diff --git a/data/tests/test_bsdf.xml b/data/tests/test_bsdf.xml index 97428d7f..bc41638e 100644 --- a/data/tests/test_bsdf.xml +++ b/data/tests/test_bsdf.xml @@ -2,6 +2,13 @@ to be tested for consistency. This is done using the testcase 'test_chisquare' --> + + + + + + diff --git a/include/mitsuba/core/quad.h b/include/mitsuba/core/quad.h index d204b4cd..40e7edc2 100644 --- a/include/mitsuba/core/quad.h +++ b/include/mitsuba/core/quad.h @@ -181,7 +181,7 @@ public: * results of the evaluation into the \c out array using \c fDim entries. */ EResult integrate(const Integrand &f, const Float *min, const Float *max, - Float *result, Float *error, size_t &evals) const; + Float *result, Float *error, size_t *evals = NULL) const; /** * \brief Integrate the function \c f over the rectangular domain @@ -211,7 +211,7 @@ public: * several hundred. */ EResult integrateVectorized(const VectorizedIntegrand &f, const Float *min, - const Float *max, Float *result, Float *error, size_t &evals) const; + const Float *max, Float *result, Float *error, size_t *evals = NULL) const; protected: size_t m_fdim, m_dim, m_maxEvals; Float m_absError, m_relError; diff --git a/include/mitsuba/core/spline.h b/include/mitsuba/core/spline.h new file mode 100644 index 00000000..62fcb8ff --- /dev/null +++ b/include/mitsuba/core/spline.h @@ -0,0 +1,89 @@ +/* + This file is part of Mitsuba, a physically based rendering system. + + Copyright (c) 2007-2011 by Wenzel Jakob and others. + + Mitsuba is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License Version 3 + as published by the Free Software Foundation. + + Mitsuba is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +*/ + +#if !defined(__SPLINE_H) +#define __SPLINE_H + +#include + +MTS_NAMESPACE_BEGIN + +/** + * \brief Simple natural cubic spline interpolation + */ +class MTS_EXPORT_CORE CubicSpline : public SerializableObject { +public: + /** + * \brief Initialize a cubic spline with the given set of + * points (must be in ascending order, with no duplicates) + */ + inline CubicSpline(std::vector &x, std::vector &y) + : m_x(x), m_y(y) { } + + /** + * \brief Initialize an empty cubic spline and reserve memory + * for the specified number of points + */ + inline CubicSpline(size_t nPoints) { + m_x.reserve(nPoints); + m_y.reserve(nPoints); + m_deriv.reserve(nPoints); + } + + /** + * \brief Unserialize a cubic spline from a + * binary data stream + */ + CubicSpline(Stream *stream, InstanceManager *manager); + + /** + * \brief Append a point at the end -- must be called with + * increasing values of \a x + */ + inline void append(Float x, Float y) { + m_x.push_back(x); + m_y.push_back(y); + } + + /// Return the number of control points + inline size_t getSize() const { return m_x.size(); } + + /// Clear the internal representation + void clear(); + + /// Compute the derivatives -- must be called prior to \ref eval + void build(); + + /// Evaluate the spline at \c x + Float eval(Float x) const; + + /// Serialize this spline to a binary data stream + void serialize(Stream *stream, InstanceManager *manager) const; + + /// Return a human-readable representation + std::string toString() const; + + MTS_DECLARE_CLASS() +private: + std::vector m_x, m_y; + std::vector m_deriv; /// 2nd derivatives +}; + +MTS_NAMESPACE_END + +#endif /* __SPLINE_H */ diff --git a/src/bsdfs/SConscript b/src/bsdfs/SConscript index 8693b6d1..73fd0a86 100644 --- a/src/bsdfs/SConscript +++ b/src/bsdfs/SConscript @@ -1,25 +1,24 @@ Import('env', 'plugins') -# Basic materials (smooth & rough versions of each) +# Basic library of smooth and rough materials plugins += env.SharedLibrary('diffuse', ['diffuse.cpp']) plugins += env.SharedLibrary('dielectric', ['dielectric.cpp']) plugins += env.SharedLibrary('conductor', ['conductor.cpp']) plugins += env.SharedLibrary('plastic', ['plastic.cpp']) - plugins += env.SharedLibrary('roughdiffuse', ['roughdiffuse.cpp']) plugins += env.SharedLibrary('roughdielectric', ['roughdielectric.cpp']) plugins += env.SharedLibrary('roughconductor', ['roughconductor.cpp']) -#plugins += env.SharedLibrary('roughplastic', ['roughplastic.cpp']) +plugins += env.SharedLibrary('roughplastic', ['roughplastic.cpp']) -# Other -plugins += env.SharedLibrary('difftrans', ['difftrans.cpp']) - -# Plugins that act as modifiers -plugins += env.SharedLibrary('mixture', ['mixture.cpp']) +# Materials that act as modifiers +#plugins += env.SharedLibrary('coating', ['coating.cpp']) #plugins += env.SharedLibrary('twosided', ['twosided.cpp']) #plugins += env.SharedLibrary('mask', ['mask.cpp']) -#plugins += env.SharedLibrary('coating', ['coating.cpp']) +plugins += env.SharedLibrary('mixture', ['mixture.cpp']) + +# Other materials +plugins += env.SharedLibrary('difftrans', ['difftrans.cpp']) #plugins += env.SharedLibrary('ward', ['ward.cpp']) #plugins += env.SharedLibrary('phong', ['phong.cpp']) #plugins += env.SharedLibrary('irawan', ['irawan.cpp']) diff --git a/src/bsdfs/coating.cpp b/src/bsdfs/coating.cpp index b77da784..b82b669b 100644 --- a/src/bsdfs/coating.cpp +++ b/src/bsdfs/coating.cpp @@ -21,7 +21,7 @@ MTS_NAMESPACE_BEGIN -/*! \plugin{coating}{Smooth coating} +/*! \plugin{coating}{Smooth dieletric coating} * * \parameters{ * \parameter{intIOR}{\Float}{Interior index of refraction \default{1.5046}} diff --git a/src/bsdfs/microfacet.h b/src/bsdfs/microfacet.h index 661adda6..50421704 100644 --- a/src/bsdfs/microfacet.h +++ b/src/bsdfs/microfacet.h @@ -19,8 +19,11 @@ #if !defined(__MICROFACET_H) #define __MICROFACET_H -#include +#include +#include +#include #include +#include MTS_NAMESPACE_BEGIN @@ -51,7 +54,7 @@ public: * \brief Create a microfacet distribution of the specified name * (ggx/phong/beckmann/as) */ - MicrofacetDistribution(const std::string &name) { + MicrofacetDistribution(const std::string &name) : m_type(EBeckmann) { std::string distr = boost::to_lower_copy(name); if (distr == "beckmann") @@ -70,6 +73,11 @@ public: /// Return the distribution type inline EType getType() const { return m_type; } + /// Is this an anisotropic microfacet distribution? + bool isAnisotropic() const { + return m_type == EAshikhminShirley; + } + /** * \brief Convert the roughness values so that they behave similarly to the * Beckmann distribution. @@ -82,13 +90,23 @@ public: value = 2 / (value * value) - 2; return std::max(value, (Float) 1e-4f); } + + /** + * \brief Implements the microfacet distribution function D + * + * \param m The microsurface normal + * \param alpha The surface roughness + */ + inline Float eval(const Vector &m, Float alpha) const { + return eval(m, alpha, alpha); + } /** * \brief Implements the microfacet distribution function D * * \param m The microsurface normal - * \param alphaU Surface roughness in the tangent directoin - * \param alphaV Surface roughness in the bitangent direction + * \param alphaU The surface roughness in the tangent direction + * \param alphaV The surface roughness in the bitangent direction */ Float eval(const Vector &m, Float alphaU, Float alphaV) const { if (Frame::cosTheta(m) <= 0) @@ -149,7 +167,20 @@ public: /** * \brief Returns the density function associated with - * the \ref{sample} function. + * the \ref sample() function. + * \param m The microsurface normal + * \param alpha The surface roughness + */ + inline Float pdf(const Vector &m, Float alpha) const { + return pdf(m, alpha, alpha); + } + + /** + * \brief Returns the density function associated with + * the \ref sample() function. + * \param m The microsurface normal + * \param alphaU The surface roughness in the tangent direction + * \param alphaV The surface roughness in the bitangent direction */ Float pdf(const Vector &m, Float alphaU, Float alphaV) const { /* Usually, this is just D(m) * cos(theta_M) */ @@ -193,8 +224,18 @@ public: * \brief Draw a sample from the microsurface normal distribution * * \param sample A uniformly distributed 2D sample - * \param alphaU Surface roughness in the tangent directoin - * \param alphaV Surface roughness in the bitangent direction + * \param alpha The surface roughness + */ + inline Normal sample(const Point2 &sample, Float alpha) const { + return MicrofacetDistribution::sample(sample, alpha, alpha); + } + + /** + * \brief Draw a sample from the microsurface normal distribution + * + * \param sample A uniformly distributed 2D sample + * \param alphaU The surface roughness in the tangent direction + * \param alphaV The surface roughness in the bitangent direction */ Normal sample(const Point2 &sample, Float alphaU, Float alphaV) const { /* The azimuthal component is always selected @@ -252,6 +293,32 @@ public: return Normal(sphericalDirection(thetaM, phiM)); } + + /** + * \brief Draw a sample from an isotropic microsurface normal + * distribution and return the magnitude of its 'z' component. + * + * \param sample A uniformly distributed number on [0,1] + * \param alphaU The surface roughness + */ + Float sampleIsotropic(Float sample, Float alpha) const { + switch (m_type) { + case EBeckmann: + return 1.0f / std::sqrt(1 + + std::abs(-alpha*alpha * std::log(1.0f - sample))); + + case EGGX: + return 1.0f / std::sqrt(1 + + alpha * alpha * sample / (1.0f - sample)); + + case EPhong: + return std::pow(sample, (Float) 1 / (alpha + 2)); + + default: + SLog(EError, "Invalid distribution function!"); + return 0.0f; + } + } /** * \brief Smith's shadow-masking function G1 for each @@ -316,6 +383,20 @@ public: * \param m The microsurface normal * \param alpha The surface roughness */ + inline Float G(const Vector &wi, const Vector &wo, const Vector &m, Float alpha) const { + return G(wi, wo, m, alpha, alpha); + } + + /** + * \brief Shadow-masking function for each of the supported + * microfacet distributions + * + * \param wi The incident direction + * \param wo The exitant direction + * \param m The microsurface normal + * \param alphaU The surface roughness in the tangent direction + * \param alphaV The surface roughness in the bitangent direction + */ Float G(const Vector &wi, const Vector &wo, const Vector &m, Float alphaU, Float alphaV) const { if (m_type != EAshikhminShirley) { return smithG1(wi, m, alphaU) @@ -330,14 +411,60 @@ public: const Float nDotM = Frame::cosTheta(m), nDotWo = Frame::cosTheta(wo), nDotWi = Frame::cosTheta(wi), - woDotM = dot(wo, m); + woDotM = dot(wo, m), + wiDotM = dot(wi, m); return std::min((Float) 1, std::min(std::abs(2 * nDotM * nDotWo / woDotM), - std::abs(2 * nDotM * nDotWi / woDotM))); + std::abs(2 * nDotM * nDotWi / wiDotM))); } } + /** + * \brief Precompute the aggregate Fresnel transmittance through + * a rough interface + * + * This function efficiently computes the integral of + * \int_{S^2} D(m) * (1 - Fr(m<->wi)) dm + * for incident directions 'wi' with a range of different inclinations + * (where Fr denotes the fresnel reflectance). It returns a cubic spline + * interpolation parameterized by the cosine of the angle between 'wi' + * and the (macro-) surface normal. + */ + CubicSpline *getRoughTransmittance(Float extIOR, Float intIOR, Float alpha, size_t resolution) const { + if (isAnisotropic()) + SLog(EError, "MicrofacetDistribution::getRoughTransmission(): only " + "supports isotropic distributions!"); + + NDIntegrator integrator(1, 2, 5000, 0, 1e-5f); + CubicSpline *spline = new CubicSpline(resolution); + size_t nEvals, nEvalsTotal = 0; + ref timer = new Timer(); + + Float stepSize = (1.0f-2*Epsilon)/(resolution-1); + for (size_t i=0; iappend(z, 1-integral); + + nEvalsTotal += nEvals; + } + SLog(EInfo, "Created a " SIZE_T_FMT "-node cubic spline approximation to the rough Frensel " + "transmittance (integration took %i ms and " SIZE_T_FMT " function evaluations)", + resolution, timer->getMilliseconds(), nEvalsTotal); + spline->build(); + return spline; + } + std::string toString() const { switch (m_type) { case EBeckmann: return "beckmann"; break; @@ -349,7 +476,26 @@ public: return ""; } } -private: +protected: + /// Integrand helper function called by \ref getRoughTransmission + void integrand(const Vector &wi, Float extIOR, Float intIOR, Float alpha, + size_t nPts, const Float *in, Float *out) const { + for (int i=0; i<(int) nPts; ++i) { + Normal m = sample(Point2(in[2*i], in[2*i+1]), alpha); + Vector wo = 2 * dot(wi, m) * Vector(m) - wi; + if (Frame::cosTheta(wi) <= 0 || Frame::cosTheta(wo) <= 0) { + out[i] = 0; + continue; + } + + /* Calculate the specular reflection component */ + out[i] = std::abs(fresnel(dot(wi, m), extIOR, intIOR) + * G(wi, wo, m, alpha) * dot(wi, m) / + (Frame::cosTheta(wi) * Frame::cosTheta(m))); + } + } + +protected: EType m_type; }; diff --git a/src/bsdfs/roughconductor.cpp b/src/bsdfs/roughconductor.cpp index 1ef66372..10ca7f39 100644 --- a/src/bsdfs/roughconductor.cpp +++ b/src/bsdfs/roughconductor.cpp @@ -24,15 +24,6 @@ MTS_NAMESPACE_BEGIN -/* Suggestion by Bruce Walter: sample the model using a slightly - wider density function. This in practice limits the importance - weights to values <= 4. - - Turned off by default, since it seems to increase the variance - of the reflection component. -*/ -#define ENLARGE_LOBE_TRICK 0 - /*!\plugin{roughconductor}{Rough conductor material} * \order{6} * \parameters{ @@ -45,7 +36,6 @@ MTS_NAMESPACE_BEGIN * \item \code{ggx}: New distribution proposed by * Walter et al. \cite{Walter07Microfacet}, which is meant to better handle * the long tails observed in measurements of ground surfaces. - * Renderings with this distribution may converge slowly. * \item \code{phong}: Classical $\cos^p\theta$ distribution. * Due to the underlying microfacet theory, * the use of this distribution here leads to more realistic @@ -289,12 +279,6 @@ public: alphaV = m_distribution.transformRoughness( m_alphaV->getValue(bRec.its).average()); -#if ENLARGE_LOBE_TRICK == 1 - Float factor = (1.2f - 0.2f * std::sqrt( - std::abs(Frame::cosTheta(bRec.wi)))); - alphaU *= factor; alphaV *= factor; -#endif - return m_distribution.pdf(H, alphaU, alphaV) / (4 * absDot(bRec.wo, H)); } @@ -311,19 +295,8 @@ public: alphaV = m_distribution.transformRoughness( m_alphaV->getValue(bRec.its).average()); -#if ENLARGE_LOBE_TRICK == 1 - Float factor = (1.2f - 0.2f * std::sqrt( - std::abs(Frame::cosTheta(bRec.wi)))); - Float sampleAlphaU = alphaU * factor, - sampleAlphaV = alphaV * factor; -#else - Float sampleAlphaU = alphaU, - sampleAlphaV = alphaV; -#endif - /* Sample M, the microsurface normal */ - const Normal m = m_distribution.sample(sample, - sampleAlphaU, sampleAlphaV); + const Normal m = m_distribution.sample(sample, alphaU, alphaV); /* Perfect specular reflection based on the microsurface normal */ bRec.wo = reflect(bRec.wi, m); @@ -337,11 +310,10 @@ public: const Spectrum F = fresnelConductor(Frame::cosTheta(bRec.wi), m_eta, m_k); - Float numerator = m_distribution.eval(m, alphaU, alphaV) - * m_distribution.G(bRec.wi, bRec.wo, m, alphaU, alphaV) + Float numerator = m_distribution.G(bRec.wi, bRec.wo, m, alphaU, alphaV) * dot(bRec.wi, m); - Float denominator = m_distribution.pdf(m, sampleAlphaU, sampleAlphaV) + Float denominator = Frame::cosTheta(m) * Frame::cosTheta(bRec.wi); return m_specularReflectance->getValue(bRec.its) * F @@ -360,19 +332,8 @@ public: alphaV = m_distribution.transformRoughness( m_alphaV->getValue(bRec.its).average()); -#if ENLARGE_LOBE_TRICK == 1 - Float factor = (1.2f - 0.2f * std::sqrt( - std::abs(Frame::cosTheta(bRec.wi)))); - Float sampleAlphaU = alphaU * factor, - sampleAlphaV = alphaV * factor; -#else - Float sampleAlphaU = alphaU, - sampleAlphaV = alphaV; -#endif - /* Sample M, the microsurface normal */ - const Normal m = m_distribution.sample(sample, - sampleAlphaU, sampleAlphaV); + const Normal m = m_distribution.sample(sample, alphaU, alphaV); /* Perfect specular reflection based on the microsurface normal */ bRec.wo = reflect(bRec.wi, m); @@ -512,10 +473,10 @@ public: << " if ((dot(wi, m) * cosTheta(wi)) <= 0 || " << endl << " (dot(wo, m) * cosTheta(wo)) <= 0)" << endl << " return 0.0;" << endl - << " float nDotM = cosTheta(m), tmp = 1.0 / dot(wo, m);" << endl + << " float nDotM = cosTheta(m);" << endl << " return min(1.0, min(" << endl - << " abs(2 * nDotM * cosTheta(wo) * tmp)," << endl - << " abs(2 * nDotM * cosTheta(wi) * tmp)));" << endl + << " abs(2 * nDotM * cosTheta(wo) / dot(wo, m))," << endl + << " abs(2 * nDotM * cosTheta(wi) / dot(wi, m))));" << endl << "}" << endl << endl << "vec3 " << evalName << "_schlick(vec3 wi) {" << endl diff --git a/src/bsdfs/roughdielectric.cpp b/src/bsdfs/roughdielectric.cpp index 8a0082ac..6ffac812 100644 --- a/src/bsdfs/roughdielectric.cpp +++ b/src/bsdfs/roughdielectric.cpp @@ -426,9 +426,9 @@ public: following seems confusing */ Float F; if (bRec.sampler) { - /* We have access to extra random numbers, hence - the exact Fresnel term with respect to the - microfacet normal is sampled */ + /* We have access to extra random numbers and are + thus able to sample the exact Fresnel term + with respect to the microfacet normal */ F = fresnel(dot(bRec.wi, H), m_extIOR, m_intIOR); } else { /* Only a simple 2D sample is given, and hence the diff --git a/src/bsdfs/roughdiffuse.cpp b/src/bsdfs/roughdiffuse.cpp index 5465c298..3340f909 100644 --- a/src/bsdfs/roughdiffuse.cpp +++ b/src/bsdfs/roughdiffuse.cpp @@ -51,7 +51,7 @@ MTS_NAMESPACE_BEGIN * reflectance at grazing angles, as well as an overall reduced contrast.}\vspace{3mm} * } * - * This reflectance model describes the interaction of light with a rough + * This reflectance model describes the interaction of light with a \emph{rough} * diffuse material, such as plaster, sand, clay, or concrete. * The underlying theory was developed by Oren and Nayar * \cite{Oren1994Generalization}, who model the microscopic surface structure as diff --git a/src/bsdfs/roughplastic.cpp b/src/bsdfs/roughplastic.cpp new file mode 100644 index 00000000..ba52d36f --- /dev/null +++ b/src/bsdfs/roughplastic.cpp @@ -0,0 +1,384 @@ +/* + This file is part of Mitsuba, a physically based rendering system. + + Copyright (c) 2007-2011 by Wenzel Jakob and others. + + Mitsuba is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License Version 3 + as published by the Free Software Foundation. + + Mitsuba is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +*/ + +#include +#include +#include +#include "microfacet.h" +#include "ior.h" + +MTS_NAMESPACE_BEGIN + +/*!\plugin{roughplastic}{Rough plastic material} + * \order{8} + * \parameters{ + * \parameter{distribution}{\String}{ + * Specifies the type of microfacet normal distribution + * used to model the surface roughness. + * \begin{enumerate}[(i)] + * \item \code{beckmann}: Physically-based distribution derived from + * Gaussian random surfaces. This is the default. + * \item \code{ggx}: New distribution proposed by + * Walter et al. \cite{Walter07Microfacet}, which is meant to better handle + * the long tails observed in measurements of ground surfaces. + * Renderings with this distribution may converge slowly. + * \item \code{phong}: Classical $\cos^p\theta$ distribution. + * Due to the underlying microfacet theory, + * the use of this distribution here leads to more realistic + * behavior than the separately available \pluginref{phong} plugin. + * \vspace{-4mm} + * \end{enumerate} + * } + * \parameter{alpha}{\Float\Or\Texture}{ + * Specifies the roughness of the unresolved surface microgeometry. + * When the Beckmann distribution is used, this parameter is equal to the + * \emph{root mean square} (RMS) slope of the microfacets. + * \default{0.1}. + * } + * \parameter{intIOR}{\Float\Or\String}{Interior index of refraction specified + * numerically or using a known material name. \default{\texttt{bk7} / 1.5046}} + * \parameter{extIOR}{\Float\Or\String}{Exterior index of refraction specified + * numerically or using a known material name. \default{\texttt{air} / 1.000277}} + * \parameter{specular\showbreak Reflectance}{\Spectrum\Or\Texture}{Optional + * factor used to modulate the specular reflectance component\default{1.0}} + * \parameter{diffuse\showbreak Reflectance}{\Spectrum\Or\Texture}{Optional + * factor used to modulate the diffuse reflectance component\default{0.5}} + * } + * \renderings{ + * \rendering{Beckmann, $\alpha=0.1$}{bsdf_roughplastic_ggx} + * \rendering{GGX, $\alpha=0.3$}{bsdf_roughplastic_ggx} + * } + * + */ +class RoughPlastic : public BSDF { +public: + RoughPlastic(const Properties &props) : BSDF(props) { + m_specularReflectance = new ConstantSpectrumTexture( + props.getSpectrum("specularReflectance", Spectrum(1.0f))); + m_diffuseReflectance = new ConstantSpectrumTexture( + props.getSpectrum("diffuseReflectance", Spectrum(0.5f))); + + /* Specifies the internal index of refraction at the interface */ + m_intIOR = lookupIOR(props, "intIOR", "bk7"); + + /* Specifies the external index of refraction at the interface */ + m_extIOR = lookupIOR(props, "extIOR", "air"); + + if (m_intIOR < 0 || m_extIOR < 0 || m_intIOR == m_extIOR) + Log(EError, "The interior and exterior indices of " + "refraction must be positive and differ!"); + + m_distribution = MicrofacetDistribution( + props.getString("distribution", "beckmann") + ); + + if (m_distribution.isAnisotropic()) + Log(EError, "The 'roughplastic' plugin does not support " + "anisotropic microfacet distributions!"); + + m_alpha = new ConstantFloatTexture( + props.getFloat("alpha", 0.1f)); + + m_usesRayDifferentials = false; + } + + RoughPlastic(Stream *stream, InstanceManager *manager) + : BSDF(stream, manager) { + m_distribution = MicrofacetDistribution( + (MicrofacetDistribution::EType) stream->readUInt() + ); + m_alpha = static_cast(manager->getInstance(stream)); + m_specularReflectance = static_cast(manager->getInstance(stream)); + m_diffuseReflectance = static_cast(manager->getInstance(stream)); + m_roughTransmittance = static_cast(manager->getInstance(stream)); + m_intIOR = stream->readFloat(); + m_extIOR = stream->readFloat(); + + m_usesRayDifferentials = + m_alpha->usesRayDifferentials() || + m_specularReflectance->usesRayDifferentials() || + m_diffuseReflectance->usesRayDifferentials(); + + configure(); + } + + void configure() { + m_components.clear(); + m_components.push_back(EGlossyReflection | EFrontSide); + m_components.push_back(EDiffuseReflection | EFrontSide); + + /* Verify the input parameters and fix them if necessary */ + m_specularReflectance = ensureEnergyConservation( + m_specularReflectance, "specularReflectance", 1.0f); + m_diffuseReflectance = ensureEnergyConservation( + m_diffuseReflectance, "diffuseReflectance", 1.0f); + + if (m_roughTransmittance == NULL) { + Float alpha = m_distribution.transformRoughness(m_alpha->getValue(Intersection()).average()); + m_roughTransmittance = m_distribution.getRoughTransmittance(m_extIOR, m_intIOR, alpha, 200); + } + + BSDF::configure(); + } + + virtual ~RoughPlastic() { } + + /// Helper function: reflect \c wi with respect to a given surface normal + inline Vector reflect(const Vector &wi, const Normal &m) const { + return 2 * dot(wi, m) * Vector(m) - wi; + } + + Spectrum eval(const BSDFQueryRecord &bRec, EMeasure measure) const { + bool sampleSpecular = (bRec.typeMask & EGlossyReflection) && + (bRec.component == -1 || bRec.component == 0); + bool sampleDiffuse = (bRec.typeMask & EDiffuseReflection) && + (bRec.component == -1 || bRec.component == 1); + + if (measure != ESolidAngle || + Frame::cosTheta(bRec.wi) <= 0 || + Frame::cosTheta(bRec.wo) <= 0 || + (!sampleSpecular && !sampleDiffuse)) + return Spectrum(0.0f); + + Spectrum result(0.0f); + if (sampleSpecular) { + /* Evaluate the roughness */ + Float alpha = m_distribution.transformRoughness( + m_alpha->getValue(bRec.its).average()); + + /* Calculate the reflection half-vector */ + const Vector H = normalize(bRec.wo+bRec.wi); + + /* Evaluate the microsurface normal distribution */ + const Float D = m_distribution.eval(H, alpha); + + /* Fresnel term */ + const Float F = fresnel(dot(bRec.wi, H), m_extIOR, m_intIOR); + + /* Smith's shadow-masking function */ + const Float G = m_distribution.G(bRec.wi, bRec.wo, H, alpha); + + /* Calculate the specular reflection component */ + Float value = F * D * G / + (4.0f * Frame::cosTheta(bRec.wi)); + + result += m_specularReflectance->getValue(bRec.its) * value; + } + + if (sampleDiffuse) + result += m_diffuseReflectance->getValue(bRec.its) * (INV_PI + * m_roughTransmittance->eval(Frame::cosTheta(bRec.wi)) + * Frame::cosTheta(bRec.wo)); + + return result; + } + + Float pdf(const BSDFQueryRecord &bRec, EMeasure measure) const { + bool sampleSpecular = (bRec.typeMask & EGlossyReflection) && + (bRec.component == -1 || bRec.component == 0); + bool sampleDiffuse = (bRec.typeMask & EDiffuseReflection) && + (bRec.component == -1 || bRec.component == 1); + + if (measure != ESolidAngle || + Frame::cosTheta(bRec.wi) <= 0 || + Frame::cosTheta(bRec.wo) <= 0 || + (!sampleSpecular && !sampleDiffuse)) + return 0.0f; + + /* Calculate the reflection half-vector */ + Vector H = normalize(bRec.wo+bRec.wi); + + Float roughTransmittance = 0.0f; + if (sampleDiffuse && sampleSpecular) + roughTransmittance = m_roughTransmittance->eval( + Frame::cosTheta(bRec.wi)); + + Float result = 0.0f; + if (sampleSpecular) { + /* Evaluate the roughness */ + Float alpha = m_distribution.transformRoughness( + m_alpha->getValue(bRec.its).average()); + + /* Jacobian of the half-direction transform */ + Float dwh_dwo = 1.0f / (4.0f * dot(bRec.wo, H)); + + /* Evaluate the microsurface normal distribution */ + Float prob = m_distribution.pdf(H, alpha); + + result += prob * dwh_dwo * + (sampleDiffuse ? (1-roughTransmittance) : 1.0f); + } + + if (sampleDiffuse) + result += Frame::cosTheta(bRec.wo) * INV_PI * + (sampleSpecular ? roughTransmittance : 1.0f); + + return result; + } + + inline Spectrum sample(BSDFQueryRecord &bRec, Float &_pdf, const Point2 &_sample) const { + bool sampleSpecular = (bRec.typeMask & EGlossyReflection) && + (bRec.component == -1 || bRec.component == 0); + bool sampleDiffuse = (bRec.typeMask & EDiffuseReflection) && + (bRec.component == -1 || bRec.component == 1); + + if (Frame::cosTheta(bRec.wi) <= 0 || (!sampleSpecular && !sampleDiffuse)) + return Spectrum(0.0f); + + bool choseReflection = sampleSpecular; + + Point2 sample(_sample); + if (sampleSpecular && sampleDiffuse) { + Float roughTransmittance = m_roughTransmittance->eval( + Frame::cosTheta(bRec.wi)); + if (sample.x < roughTransmittance) { + sample.x /= roughTransmittance; + choseReflection = false; + } else { + sample.x = (sample.x - roughTransmittance) + / (1 - roughTransmittance); + } + } + + if (choseReflection) { + /* Evaluate the roughness */ + Float alpha = m_distribution.transformRoughness( + m_alpha->getValue(bRec.its).average()); + + /* Sample M, the microsurface normal */ + const Normal m = m_distribution.sample(sample, alpha); + + + /* Perfect specular reflection based on the microsurface normal */ + bRec.wo = reflect(bRec.wi, m); + bRec.sampledComponent = 0; + bRec.sampledType = EGlossyReflection; + + /* Side check */ + if (Frame::cosTheta(bRec.wo) <= 0) + return Spectrum(0.0f); + } else { + bRec.sampledComponent = 1; + bRec.sampledType = EDiffuseReflection; + bRec.wo = squareToHemispherePSA(sample); + } + + /* Guard against numerical imprecisions */ + _pdf = pdf(bRec, ESolidAngle); + + if (_pdf == 0) + return Spectrum(0.0f); + else + return eval(bRec, ESolidAngle); + } + + void addChild(const std::string &name, ConfigurableObject *child) { + if (child->getClass()->derivesFrom(MTS_CLASS(Texture)) && name == "alpha") { + m_alpha = static_cast(child); + m_usesRayDifferentials |= m_alpha->usesRayDifferentials(); + } else if (child->getClass()->derivesFrom(MTS_CLASS(Texture)) && name == "specularReflectance") { + m_specularReflectance = static_cast(child); + m_usesRayDifferentials |= m_specularReflectance->usesRayDifferentials(); + } else if (child->getClass()->derivesFrom(MTS_CLASS(Texture)) && name == "diffuseReflectance") { + m_diffuseReflectance = static_cast(child); + m_usesRayDifferentials |= m_diffuseReflectance->usesRayDifferentials(); + } else { + BSDF::addChild(name, child); + } + } + + Spectrum sample(BSDFQueryRecord &bRec, const Point2 &sample) const { + Float pdf; + Spectrum result = RoughPlastic::sample(bRec, pdf, sample); + + if (result.isZero()) + return Spectrum(0.0f); + else + return result / pdf; + } + + void serialize(Stream *stream, InstanceManager *manager) const { + BSDF::serialize(stream, manager); + + stream->writeUInt((uint32_t) m_distribution.getType()); + manager->serialize(stream, m_alpha.get()); + manager->serialize(stream, m_specularReflectance.get()); + manager->serialize(stream, m_diffuseReflectance.get()); + manager->serialize(stream, m_roughTransmittance.get()); + stream->writeFloat(m_intIOR); + stream->writeFloat(m_extIOR); + } + + std::string toString() const { + std::ostringstream oss; + oss << "RoughPlastic[" << endl + << " name = \"" << getName() << "\"," << endl + << " distribution = " << m_distribution.toString() << "," << endl + << " alpha = " << indent(m_alpha->toString()) << "," << endl + << " specularReflectance = " << indent(m_specularReflectance->toString()) << "," << endl + << " diffuseReflectance = " << indent(m_diffuseReflectance->toString()) << "," << endl + << " intIOR = " << m_intIOR << "," << endl + << " extIOR = " << m_extIOR << endl + << "]"; + return oss.str(); + } + + Shader *createShader(Renderer *renderer) const; + + MTS_DECLARE_CLASS() +private: + MicrofacetDistribution m_distribution; + ref m_roughTransmittance; + ref m_diffuseReflectance; + ref m_specularReflectance; + ref m_alpha; + Float m_intIOR, m_extIOR; +}; + +/* Fake plastic shader -- it is really hopeless to visualize + this material in the VPL renderer, so let's try to do at least + something that suggests the presence of a translucent boundary */ +class RoughPlasticShader : public Shader { +public: + RoughPlasticShader(Renderer *renderer) : + Shader(renderer, EBSDFShader) { + m_flags = ETransparent; + } + + void generateCode(std::ostringstream &oss, + const std::string &evalName, + const std::vector &depNames) const { + oss << "vec3 " << evalName << "(vec2 uv, vec3 wi, vec3 wo) {" << endl + << " return vec3(0.08);" << endl + << "}" << endl + << endl + << "vec3 " << evalName << "_diffuse(vec2 uv, vec3 wi, vec3 wo) {" << endl + << " return " << evalName << "(uv, wi, wo);" << endl + << "}" << endl; + } + MTS_DECLARE_CLASS() +}; + +Shader *RoughPlastic::createShader(Renderer *renderer) const { + return new RoughPlasticShader(renderer); +} + +MTS_IMPLEMENT_CLASS(RoughPlasticShader, false, Shader) +MTS_IMPLEMENT_CLASS_S(RoughPlastic, false, BSDF) +MTS_EXPORT_PLUGIN(RoughPlastic, "Rough plastic BSDF"); +MTS_NAMESPACE_END diff --git a/src/libcore/SConscript b/src/libcore/SConscript index e6471440..e88b988d 100644 --- a/src/libcore/SConscript +++ b/src/libcore/SConscript @@ -33,7 +33,7 @@ libcore_objects = [ 'serialization.cpp', 'sstream.cpp', 'cstream.cpp', 'mstream.cpp', 'sched.cpp', 'sched_remote.cpp', 'sshstream.cpp', 'wavelet.cpp', 'zstream.cpp', 'shvector.cpp', 'fresolver.cpp', 'quad.cpp', 'mmap.cpp', - 'chisquare.cpp' + 'chisquare.cpp', 'spline.cpp' ] # Add some platform-specific components diff --git a/src/libcore/chisquare.cpp b/src/libcore/chisquare.cpp index 870f37f5..b077055a 100644 --- a/src/libcore/chisquare.cpp +++ b/src/libcore/chisquare.cpp @@ -147,11 +147,10 @@ m_table[thetaBin * m_phiBins + phiBin] += boost::get<1>(sample); min[1] = j * factor.y; max[1] = (j+1) * factor.y; Float result, error; - size_t evals; integrator.integrateVectorized( boost::bind(&ChiSquare::integrand, pdfFn, _1, _2, _3), - min, max, &result, &error, evals + min, max, &result, &error ); integral += result; diff --git a/src/libcore/quad.cpp b/src/libcore/quad.cpp index af85b900..46cae28e 100644 --- a/src/libcore/quad.cpp +++ b/src/libcore/quad.cpp @@ -1148,17 +1148,25 @@ NDIntegrator::NDIntegrator(size_t fDim, size_t dim, m_relError(relError) { } NDIntegrator::EResult NDIntegrator::integrate(const Integrand &f, const Float *min, - const Float *max, Float *result, Float *error, size_t &evals) const { + const Float *max, Float *result, Float *error, size_t *_evals) const { VectorizationAdapter adapter(f, m_fdim, m_dim); - return mitsuba::integrate((unsigned int) m_fdim, boost::bind( + size_t evals = 0; + EResult retval = mitsuba::integrate((unsigned int) m_fdim, boost::bind( &VectorizationAdapter::f, &adapter, _1, _2, _3), (unsigned int) m_dim, min, max, m_maxEvals, m_absError, m_relError, result, error, evals, false); + if (_evals) + *_evals = evals; + return retval; } NDIntegrator::EResult NDIntegrator::integrateVectorized(const VectorizedIntegrand &f, const Float *min, - const Float *max, Float *result, Float *error, size_t &evals) const { - return mitsuba::integrate((unsigned int) m_fdim, f, (unsigned int) m_dim, + const Float *max, Float *result, Float *error, size_t *_evals) const { + size_t evals = 0; + EResult retval = mitsuba::integrate((unsigned int) m_fdim, f, (unsigned int) m_dim, min, max, m_maxEvals, m_absError, m_relError, result, error, evals, true); + if (_evals) + *_evals = evals; + return retval; } MTS_NAMESPACE_END diff --git a/src/libcore/spline.cpp b/src/libcore/spline.cpp new file mode 100644 index 00000000..308ac09a --- /dev/null +++ b/src/libcore/spline.cpp @@ -0,0 +1,114 @@ +/* + This file is part of Mitsuba, a physically based rendering system. + + Copyright (c) 2007-2011 by Wenzel Jakob and others. + + Mitsuba is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License Version 3 + as published by the Free Software Foundation. + + Mitsuba is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +*/ + +#include + +MTS_NAMESPACE_BEGIN + +CubicSpline::CubicSpline(Stream *stream, InstanceManager *manager) { + size_t size = stream->readSize(); + m_x.resize(size); m_y.resize(size); m_deriv.resize(size); + stream->readFloatArray(&m_x[0], size); + stream->readFloatArray(&m_y[0], size); + stream->readFloatArray(&m_deriv[0], size); +} + +void CubicSpline::build() { + if (m_x.size() != m_y.size()) + Log(EError, "build(): encountered a different number " + " of X and Y components!"); + size_t size = m_x.size(); + if (size < 3) + Log(EError, "build(): need at least three points!"); + + Float *temp = new Float[size]; + m_deriv.resize(size); + + /* Solve a tridiagonal system (based on Numerical Recipes) */ + m_deriv[0] = temp[0] = 0.0f; + for(size_t i=1; i=0; k--) + m_deriv[k] = m_deriv[k] * m_deriv[k+1] + temp[k]; + + delete[] temp; +} + +void CubicSpline::clear() { + m_x.clear(); + m_y.clear(); + m_deriv.clear(); +} + +Float CubicSpline::eval(Float x) const { + size_t lo = 0, hi = m_x.size()-1; + + /* Binary search */ + while (hi-lo > 1){ + size_t k = (hi+lo) >> 1; + if (m_x[k] > x) + hi = k; + else + lo = k; + } + + const Float h = m_x[hi]-m_x[lo], + invH = 1.0f / h, + a = (m_x[hi]-x) * invH, + b = (x- m_x[lo]) * invH; + + return a*m_y[lo] + b*m_y[hi] + (1.0f/6.0f) * + ((a*a*a-a)*m_deriv[lo] + (b*b*b-b)*m_deriv[hi]) * (h*h); +} + +void CubicSpline::serialize(Stream *stream, InstanceManager *manager) const { + Assert(m_x.size() == m_y.size() && m_x.size() == m_deriv.size()); + stream->writeSize(m_x.size()); + stream->writeFloatArray(&m_x[0], m_x.size()); + stream->writeFloatArray(&m_y[0], m_x.size()); + stream->writeFloatArray(&m_deriv[0], m_x.size()); +} + +std::string CubicSpline::toString() const { + std::ostringstream oss; + Assert(m_x.size() == m_y.size()); + oss << "CubicSpline[" << endl + << " nodeCount = " << m_x.size() << "," << endl + << " nodes = {" << endl; + for (size_t i=0; i " << m_y[i]; + if (i+1 < m_x.size()) + oss << ","; + oss << endl; + } + oss << " }" << endl + << "]"; + return oss.str(); +} + +MTS_IMPLEMENT_CLASS(CubicSpline, false, SerializableObject) +MTS_NAMESPACE_END diff --git a/src/tests/test_quad.cpp b/src/tests/test_quad.cpp index e820582e..72d3cb6f 100644 --- a/src/tests/test_quad.cpp +++ b/src/tests/test_quad.cpp @@ -72,7 +72,7 @@ public: Float min = 0, max = 10, result, err; size_t evals; assertTrue(quad.integrate(boost::bind( - &TestQuadrature::testF2, this, _1, _2), &min, &max, &result, &err, evals) == NDIntegrator::ESuccess); + &TestQuadrature::testF2, this, _1, _2), &min, &max, &result, &err, &evals) == NDIntegrator::ESuccess); Float ref = 2 * std::pow(std::sin(5.0f), 2.0f); Log(EInfo, "test02_nD_01(): used " SIZE_T_FMT " function evaluations, " "error=%f", evals, err); @@ -84,7 +84,7 @@ public: size_t evals; Float min[3] = { -1, -1, -1 } , max[3] = { 1, 1, 1 }, result[2], err[2]; assertTrue(quad.integrateVectorized(boost::bind( - &TestQuadrature::testF3, this, _1, _2, _3), min, max, result, err, evals) == NDIntegrator::ESuccess); + &TestQuadrature::testF3, this, _1, _2, _3), min, max, result, err, &evals) == NDIntegrator::ESuccess); Log(EInfo, "test02_nD_02(): used " SIZE_T_FMT " function evaluations, " "error=[%f, %f]", evals, err[0], err[1]); assertEqualsEpsilon(result[0], 1, 1e-5f);