partially working implementation of the rough diffuse model, added a class for representing cubic splines

data/tests/test_bsdf.xml

@@ -2,6 +2,13 @@
 	 to be tested for consistency. This is done
 	 using the testcase 'test_chisquare' -->
 <scene>
+	<!-- Test the rough plastic model with the 
+		 Beckmann microfacet distribution -->
+	<bsdf type="roughplastic">
+		<string name="distribution" value="beckmann"/>
+		<float name="alpha" value=".3"/>
+	</bsdf>
+
 	<!-- Test the smooth diffuse model -->
 	<bsdf type="diffuse"/>
 

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;

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 <http://www.gnu.org/licenses/>.
+*/
+
+#if !defined(__SPLINE_H)
+#define __SPLINE_H
+
+#include <mitsuba/core/serialization.h>
+
+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<Float> &x, std::vector<Float> &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<Float> m_x, m_y;
+	std::vector<Float> m_deriv; /// 2nd derivatives
+};
+
+MTS_NAMESPACE_END
+
+#endif /* __SPLINE_H */

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
+# Materials that act as modifiers
-plugins += env.SharedLibrary('difftrans', ['difftrans.cpp'])
+#plugins += env.SharedLibrary('coating', ['coating.cpp'])
-
-# Plugins that act as modifiers
-plugins += env.SharedLibrary('mixture', ['mixture.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'])

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}}

src/bsdfs/microfacet.h

@@ -19,8 +19,11 @@
 #if !defined(__MICROFACET_H)
 #define __MICROFACET_H
 
-#include <mitsuba/mitsuba.h>
+#include <mitsuba/core/quad.h>
+#include <mitsuba/core/timer.h>
+#include <mitsuba/core/spline.h>
 #include <boost/algorithm/string.hpp>
+#include <boost/bind.hpp>
 
 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.
@@ -87,8 +95,18 @@ public:
 	 * \brief Implements the microfacet distribution function D
 	 *
 	 * \param m The microsurface normal
-	 * \param alphaU  Surface roughness in the tangent directoin
+	 * \param alpha  The surface roughness
-	 * \param alphaV  Surface roughness in the bitangent direction
+	 */
+	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  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 alpha  The surface roughness 
-	 * \param alphaV  Surface roughness in the bitangent direction
+	 */
+	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 
@@ -253,6 +294,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
 	 * of the supported microfacet distributions
@@ -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> timer = new Timer();
+
+		Float stepSize = (1.0f-2*Epsilon)/(resolution-1);
+		for (size_t i=0; i<resolution; ++i) {
+			Float z = stepSize * i + Epsilon;
+			Vector wi(std::sqrt(std::max((Float) 0, 1-z*z)), 0, z);
+			Float min[2] = {0, 0}, max[2] = {1, 1},
+				  integral = 0, error = 0;
+
+			integrator.integrateVectorized(
+				boost::bind(&MicrofacetDistribution::integrand, this,
+					wi, extIOR, intIOR, alpha, _1, _2, _3),
+				min, max, &integral, &error, &nEvals
+			);
+
+			spline->append(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;
 };
 

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,
+		const Normal m = m_distribution.sample(sample, alphaU, alphaV);
-				sampleAlphaU, sampleAlphaV);
 
 		/* 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)
+		Float numerator = m_distribution.G(bRec.wi, bRec.wo, m, alphaU, alphaV)
-			* 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,
+		const Normal m = m_distribution.sample(sample, alphaU, alphaV);
-				sampleAlphaU, sampleAlphaV);
 
 		/* 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(wo) / dot(wo, m))," << endl
-			<< "        abs(2 * nDotM * cosTheta(wi) * tmp)));" << endl
+			<< "        abs(2 * nDotM * cosTheta(wi) / dot(wi, m))));" << endl
 			<< "}" << endl
 			<< endl
 			<< "vec3 " << evalName << "_schlick(vec3 wi) {" << endl

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
+				/* We have access to extra random numbers and are
-				   the exact Fresnel term with respect to the
+				   thus able to sample the exact Fresnel term 
-				   microfacet normal is sampled */
+				   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

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 

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 <http://www.gnu.org/licenses/>.
+*/
+
+#include <mitsuba/render/bsdf.h>
+#include <mitsuba/render/sampler.h>
+#include <mitsuba/hw/basicshader.h>
+#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<Texture *>(manager->getInstance(stream));
+		m_specularReflectance = static_cast<Texture *>(manager->getInstance(stream));
+		m_diffuseReflectance = static_cast<Texture *>(manager->getInstance(stream));
+		m_roughTransmittance = static_cast<CubicSpline *>(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<Texture *>(child);
+			m_usesRayDifferentials |= m_alpha->usesRayDifferentials();
+		} else if (child->getClass()->derivesFrom(MTS_CLASS(Texture)) && name == "specularReflectance") {
+			m_specularReflectance = static_cast<Texture *>(child);
+			m_usesRayDifferentials |= m_specularReflectance->usesRayDifferentials();
+		} else if (child->getClass()->derivesFrom(MTS_CLASS(Texture)) && name == "diffuseReflectance") {
+			m_diffuseReflectance = static_cast<Texture *>(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<CubicSpline> m_roughTransmittance;
+	ref<Texture> m_diffuseReflectance;
+	ref<Texture> m_specularReflectance;
+	ref<Texture> 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<std::string> &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

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

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;

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 {
+		const Float *max, Float *result, Float *error, size_t *_evals) const {
-	return mitsuba::integrate((unsigned int) m_fdim, f, (unsigned int) m_dim,
+	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

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 <http://www.gnu.org/licenses/>.
+*/
+
+#include <mitsuba/core/spline.h>
+
+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<size-1; i++){
+		Float sig = (m_x[i] - m_x[i-1]) / (m_x[i+1] - m_x[i-1]);
+		Float invP = 1.0f / (sig * m_deriv[i-1] + 2);
+		m_deriv[i]=(sig-1) * invP;
+		temp[i]=(m_y[i+1]-m_y[i])  / (m_x[i+1]-m_x[i]) 
+		       - (m_y[i]-m_y[i-1]) / (m_x[i]-m_x[i-1]);
+		temp[i]=(6*temp[i] / (m_x[i+1]-m_x[i-1]) - sig*temp[i-1]) * invP;
+	}
+	m_deriv[size-1] = 0.0f;
+
+	/* Backsubstitute */
+	for(ptrdiff_t k=size-2; k>=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_x.size(); ++i) {
+		oss << "    " << m_x[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

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);