fancy roughplastic sampling strategy, just before cleanup

data/tests/test_bsdf.xml

@@ -6,7 +6,7 @@
 		 Beckmann microfacet distribution -->
 	<bsdf type="roughplastic">
 		<string name="distribution" value="beckmann"/>
-		<float name="alpha" value=".3"/>
+		<float name="alpha" value=".7"/>
 	</bsdf>
 
 	<!-- Test the smooth diffuse model -->

include/mitsuba/render/mipmap.h

@@ -56,7 +56,8 @@ public:
 	/// Construct a mip map from a HDR bitmap
 	static ref<MIPMap> fromBitmap(Bitmap *bitmap, 
 		EFilterType filterType = EEWA, EWrapMode wrapMode = ERepeat,
-		Float maxAnisotropy = 8.0f);
+		Float maxAnisotropy = 8.0f, 
+		Spectrum::EConversionIntent intent = Spectrum::EReflectance);
 
 	/// Do a mip-map lookup at the appropriate level
 	Spectrum getValue(Float u, Float v,

src/bsdfs/microfacet.h

@@ -469,7 +469,11 @@ public:
 
 	/**
 	 * \brief Compute a spline representation that gives the probability
-	 * of sampling a transmission event e.g. in the plugin 'roughdielectric'.
+	 * of choosing a reflection event when importance sampling wrt. the
+	 * Fresnel coefficient between a sampled microsurface normal and the
+	 * incident direction.
+	 *
+	 * This function is currently used by the plugin 'roughplastic'.
 	 *
 	 * Like \ref computeRoughTransmittance, the spline is parameterized by the
 	 * cosine of the angle between the indident direction and the (macro-) 
@@ -477,7 +481,8 @@ public:
 	 * 
 	 * \remark This function only works for isotropic microfacet distributions
 	 */
-	CubicSpline *computeTransmissionProbability(Float extIOR, Float intIOR, Float alpha, size_t resolution) const {
+	CubicSpline *computeTransmissionProbability(Float extIOR, Float intIOR, 
+			Float alpha, Float specularSamplingWeight, size_t resolution) const {
 		if (isAnisotropic())
 			SLog(EError, "MicrofacetDistribution::computeTransmissionProbability(): only "
 				"supports isotropic distributions!");
@@ -496,7 +501,7 @@ public:
 
 			integrator.integrateVectorized(
 				boost::bind(&MicrofacetDistribution::integrand2, this,
-					wi, extIOR, intIOR, alpha, _1, _2, _3),
+					wi, extIOR, intIOR, alpha, specularSamplingWeight, _1, _2, _3),
 				min, max, &integral, &error, &nEvals
 			);
 
@@ -545,17 +550,14 @@ protected:
 
 	/// Integrand helper function called by \ref computeTransmissionProbability
 	void integrand2(const Vector &wi, Float extIOR, Float intIOR, Float alpha,
-			size_t nPts, const Float *in, Float *out) const {
+			Float specularSamplingWeight, 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;
+			Float probSpecular = fresnel(dot(wi, m), extIOR, intIOR);
-			if (Frame::cosTheta(wo) <= 0) {
+			probSpecular = (probSpecular*specularSamplingWeight) /
-				out[i] = 0;
+				(probSpecular*specularSamplingWeight + 
-				continue;
+				(1-probSpecular) * (1-specularSamplingWeight));
-			}
+			out[i] = 1-probSpecular;
-
-			/* Calculate the specular reflection component */
-			out[i] = 1 - fresnel(dot(wi, m), extIOR, intIOR);
 		}
 	}
 

src/bsdfs/roughplastic.cpp

@@ -24,6 +24,8 @@
 
 MTS_NAMESPACE_BEGIN
 
+#define SPLINE_PRECOMP_NODES 200
+
 /*!\plugin{roughplastic}{Rough plastic material}
  * \order{8}
  * \parameters{
@@ -44,7 +46,7 @@ MTS_NAMESPACE_BEGIN
  *              \vspace{-4mm}
  *       \end{enumerate}
  *     }
- *     \parameter{alpha}{\Float\Or\Texture}{
+ *     \parameter{alpha}{\Float}{
  *         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. 
@@ -91,7 +93,7 @@ public:
 			Log(EError, "The 'roughplastic' plugin does not support "
 				"anisotropic microfacet distributions!");
 
-		m_alpha = new ConstantFloatTexture(
+		m_alpha = m_distribution.transformRoughness(
 			props.getFloat("alpha", 0.1f));
 
 		m_usesRayDifferentials = false;
@@ -102,15 +104,15 @@ public:
 		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_diffuseProb = static_cast<CubicSpline *>(manager->getInstance(stream));
+		m_alpha = stream->readFloat();
 		m_intIOR = stream->readFloat();
 		m_extIOR = stream->readFloat();
 
 		m_usesRayDifferentials = 
-			m_alpha->usesRayDifferentials() ||
 			m_specularReflectance->usesRayDifferentials() ||
 			m_diffuseReflectance->usesRayDifferentials();
 
@@ -119,8 +121,8 @@ public:
 
 	void configure() {
 		m_components.clear();
-		m_components.push_back(EGlossyReflection | EFrontSide);
+		m_components.push_back(EGlossyReflection | ECanUseSampler | EFrontSide);
-		m_components.push_back(EDiffuseReflection | EFrontSide);
+		m_components.push_back(EDiffuseReflection | ECanUseSampler | EFrontSide);
 
 		/* Verify the input parameters and fix them if necessary */
 		m_specularReflectance = ensureEnergyConservation(
@@ -128,10 +130,22 @@ public:
 		m_diffuseReflectance = ensureEnergyConservation(
 			m_diffuseReflectance, "diffuseReflectance", 1.0f);
 
-		if (m_roughTransmittance == NULL) {
+		/* Compute weights that further steer samples towards
-			Float alpha = m_distribution.transformRoughness(m_alpha->getValue(Intersection()).average());
+		   the specular or diffuse components */
-			m_roughTransmittance = m_distribution.computeRoughTransmittance(m_extIOR, m_intIOR, alpha, 200);
+		Float dAvg = m_diffuseReflectance->getAverage().getLuminance(),
-		}
+			  sAvg = m_specularReflectance->getAverage().getLuminance();
+		m_specularSamplingWeight = sAvg / (dAvg + sAvg);
+
+		/* Precompute the rough transmittance through the interface */
+		m_roughTransmittance = m_distribution.computeRoughTransmittance(
+				m_extIOR, m_intIOR, m_alpha, SPLINE_PRECOMP_NODES);
+
+		/* Precompute a spline that specifies the probability of
+		   sampling the diffuse component for different angles
+		   of incidence. */
+		m_diffuseProb = m_distribution.computeTransmissionProbability(
+				m_extIOR, m_intIOR, m_alpha, m_specularSamplingWeight,
+				SPLINE_PRECOMP_NODES);
 
 		BSDF::configure();
 	}
@@ -157,21 +171,17 @@ public:
 	
 		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);
+			const Float D = m_distribution.eval(H, m_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);
+			const Float G = m_distribution.G(bRec.wi, bRec.wo, H, m_alpha);
 
 			/* Calculate the specular reflection component */
 			Float value = F * D * G / 
@@ -201,32 +211,48 @@ public:
 			return 0.0f;
 
 		/* Calculate the reflection half-vector */
-		Vector H = normalize(bRec.wo+bRec.wi);
+		const Vector H = normalize(bRec.wo+bRec.wi);
 
-		Float roughTransmittance = 0.0f;
+		Float probSpecular, probDiffuse;
-		if (sampleDiffuse && sampleSpecular)
+		if (sampleSpecular && sampleDiffuse) {
-			roughTransmittance = m_roughTransmittance->eval(
+			if (bRec.sampler && false) {
-				Frame::cosTheta(bRec.wi));
+				/* Fancy sampling strategy */
+				probSpecular = fresnel(dot(bRec.wi, H), m_extIOR, m_intIOR);
+				
+				/* Reallocate samples */
+				probSpecular = (probSpecular*m_specularSamplingWeight) /
+					(probSpecular*m_specularSamplingWeight + 
+					(1-probSpecular) * (1-m_specularSamplingWeight));
+	
+				probDiffuse = m_diffuseProb->eval(Frame::cosTheta(bRec.wi));
+			} else {
+				/* Basic sampling strategy that only needs 2 random numbers */
+				probSpecular = 1 - m_roughTransmittance->eval(Frame::cosTheta(bRec.wi));
+
+				/* Reallocate samples */
+				probSpecular = (probSpecular*m_specularSamplingWeight) /
+					(probSpecular*m_specularSamplingWeight + 
+					(1-probSpecular) * (1-m_specularSamplingWeight));
+
+				probDiffuse = 1 - probSpecular;
+			}
+		} else {
+			probDiffuse = probSpecular = 1.0f;
+		}
 
 		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));
+			const Float dwh_dwo = 1.0f / (4.0f * dot(bRec.wo, H));
 
 			/* Evaluate the microsurface normal distribution */
-			Float prob = m_distribution.pdf(H, alpha);
+			const Float prob = m_distribution.pdf(H, m_alpha);
 
-			result += prob * dwh_dwo *
+			result = prob * dwh_dwo * probSpecular;
-				(sampleDiffuse ? (1-roughTransmittance) : 1.0f);
 		}
 
 		if (sampleDiffuse) 
-			result += Frame::cosTheta(bRec.wo) * INV_PI *
+			result += Frame::cosTheta(bRec.wo) * INV_PI * probDiffuse;
-				(sampleSpecular ? roughTransmittance : 1.0f);
 
 		return result;
 	}
@@ -240,30 +266,57 @@ public:
 		if (Frame::cosTheta(bRec.wi) <= 0 || (!sampleSpecular && !sampleDiffuse))
 			return Spectrum(0.0f);
 
-		bool choseReflection = sampleSpecular;
 
+		bool choseSpecular = sampleSpecular;
+		Normal m;
 		Point2 sample(_sample);
+
 		if (sampleSpecular && sampleDiffuse) {
-			Float roughTransmittance = m_roughTransmittance->eval(
+			if (bRec.sampler && false) {
-				Frame::cosTheta(bRec.wi));
+				/**
-			if (sample.x < roughTransmittance) {
+				 * We have access to a sampler -- use a good sampling 
-				sample.x /= roughTransmittance;
+				 * technique, which is somewhat wasteful in terms of
-				choseReflection = false;
+				 * random numbers
+				 */
+				m = m_distribution.sample(sample, m_alpha);
+
+				Float probSpecular = fresnel(dot(bRec.wi, m), m_extIOR, m_intIOR);
+
+				/* Reallocate samples */
+				probSpecular = (probSpecular*m_specularSamplingWeight) /
+					(probSpecular*m_specularSamplingWeight + 
+					(1-probSpecular) * (1-m_specularSamplingWeight));
+
+				if (bRec.sampler->next1D() > probSpecular) {
+					choseSpecular = false;
+					sample = bRec.sampler->next2D();
+				}
 			} else {
-				sample.x = (sample.x - roughTransmittance)
+				/**
-					/ (1 - roughTransmittance);
+				 * Basic strategy -- use a clamped Fresnel coefficient 
+				 * wrt. the macro-surface normal to choose between 
+				 * diffuse and specular component. 
+				 */
+				Float probSpecular = 1 - m_roughTransmittance->eval(Frame::cosTheta(bRec.wi));
+
+				/* Reallocate samples */
+				probSpecular = (probSpecular*m_specularSamplingWeight) /
+					(probSpecular*m_specularSamplingWeight + 
+					(1-probSpecular) * (1-m_specularSamplingWeight));
+				
+				if (sample.x < probSpecular) {
+					sample.x /= probSpecular;
+					m = m_distribution.sample(sample, m_alpha);
+				} else {
+					sample.x = (sample.x - probSpecular) / (1 - probSpecular);
+					choseSpecular = false;
+				}
 			}
+		} else if (choseSpecular) {
+			m = m_distribution.sample(sample, m_alpha);
 		}
 
-		if (choseReflection) {
+		if (choseSpecular) {
-			/* 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;
@@ -288,10 +341,7 @@ public:
 	}
 
 	void addChild(const std::string &name, ConfigurableObject *child) {
-		if (child->getClass()->derivesFrom(MTS_CLASS(Texture)) && name == "alpha") {
+		if (child->getClass()->derivesFrom(MTS_CLASS(Texture)) && name == "specularReflectance") {
-			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") {
@@ -316,10 +366,11 @@ public:
 		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());
+		manager->serialize(stream, m_diffuseProb.get());
+		stream->writeFloat(m_alpha);
 		stream->writeFloat(m_intIOR);
 		stream->writeFloat(m_extIOR);
 	}
@@ -329,9 +380,11 @@ public:
 		oss << "RoughPlastic[" << endl
 			<< "  name = \"" << getName() << "\"," << endl
 			<< "  distribution = " << m_distribution.toString() << "," << endl
-			<< "  alpha = " << indent(m_alpha->toString()) << "," << endl
+			<< "  alpha = " << m_alpha << "," << endl
 			<< "  specularReflectance = " << indent(m_specularReflectance->toString()) << "," << endl
 			<< "  diffuseReflectance = " << indent(m_diffuseReflectance->toString()) << "," << endl
+			<< "  specularSamplingWeight = " << m_specularSamplingWeight << "," << endl
+			<< "  diffuseSamplingWeight = " << (1-m_specularSamplingWeight) << "," << endl
 			<< "  intIOR = " << m_intIOR << "," << endl
 			<< "  extIOR = " << m_extIOR << endl
 			<< "]";
@@ -344,10 +397,11 @@ public:
 private:
 	MicrofacetDistribution m_distribution;
 	ref<CubicSpline> m_roughTransmittance;
+	ref<CubicSpline> m_diffuseProb;
 	ref<Texture> m_diffuseReflectance;
 	ref<Texture> m_specularReflectance;
-	ref<Texture> m_alpha;
+	Float m_alpha, m_intIOR, m_extIOR;
-	Float m_intIOR, m_extIOR;
+	Float m_specularSamplingWeight;
 };
 
 /* Fake plastic shader -- it is really hopeless to visualize

src/librender/mipmap.cpp

@@ -154,7 +154,8 @@ Spectrum MIPMap::getMaximum() const {
 }
 	
 ref<MIPMap> MIPMap::fromBitmap(Bitmap *bitmap, EFilterType filterType,
-		EWrapMode wrapMode, Float maxAnisotropy) {
+		EWrapMode wrapMode, Float maxAnisotropy,
+		Spectrum::EConversionIntent intent) {
 	int width = bitmap->getWidth();
 	int height = bitmap->getHeight();
 	float *data = bitmap->getFloatData();
@@ -165,9 +166,9 @@ ref<MIPMap> MIPMap::fromBitmap(Bitmap *bitmap, EFilterType filterType,
 			float r = data[(y*width+x)*4+0];
 			float g = data[(y*width+x)*4+1];
 			float b = data[(y*width+x)*4+2];
-			s.fromLinearRGB(r, g, b);
-			s.clampNegative();
 			/* Convert to a spectral representation */
+			s.fromLinearRGB(r, g, b, intent);
+			s.clampNegative();
 			pixels[y*width+x] = s;
 		}
 	}

src/luminaires/envmap.cpp

@@ -44,7 +44,8 @@ public:
 		ref<Stream> is = new FileStream(m_path, FileStream::EReadOnly);
 		ref<Bitmap> bitmap = new Bitmap(Bitmap::EEXR, is);
 
-		m_mipmap = MIPMap::fromBitmap(bitmap);
+		m_mipmap = MIPMap::fromBitmap(bitmap, MIPMap::ETrilinear,
+				MIPMap::ERepeat, 0.0f, Spectrum::EIlluminant);
 		m_average = m_mipmap->triangle(m_mipmap->getLevels()-1, 0, 0) * m_intensityScale;
 		m_type = EOnSurface;
 	}
@@ -60,7 +61,8 @@ public:
 		stream->copyTo(mStream, size);
 		mStream->setPos(0);
 		ref<Bitmap> bitmap = new Bitmap(Bitmap::EEXR, mStream);
-		m_mipmap = MIPMap::fromBitmap(bitmap);
+		m_mipmap = MIPMap::fromBitmap(bitmap, MIPMap::ETrilinear,
+				MIPMap::ERepeat, 0.0f, Spectrum::EIlluminant);
 		m_average = m_mipmap->triangle(m_mipmap->getLevels()-1, 0, 0) * m_intensityScale;
 		m_surfaceArea = 4 * m_bsphere.radius * m_bsphere.radius * M_PI;
 		m_invSurfaceArea = 1/m_surfaceArea;