ported all of the BRDF/BSDF models, some bugs remain

data/tests/test_bsdf.xml

@@ -2,26 +2,8 @@
 	 to be tested for consistency. This is done
 	 using the testcase 'test_chisquare' -->
 <scene version="0.3.0">
-	<!-- Test the coating model with the Hanrahan-Krueger model -->
+	<!-- Test the smooth plastic model -->
-	<bsdf type="coating">
+	<bsdf type="plastic"/>
-		<bsdf type="hk">
-			<rgb name="sigmaA" value="0.1 0.2 0.3"/>
-			<rgb name="sigmaS" value="1"/>
-			<float name="thickness" value="2"/>
-		</bsdf>
-	</bsdf>
-
-	<!-- Test the Hanrahan-Krueger model with an
-		 isotropic phase function -->
-	<bsdf type="hk"/>
-	
-	<!-- Test the Hanrahan-Krueger model with a 
-		 forward-scattering phase function -->
-	<bsdf type="hk">
-		<phase type="hg">
-			<float name="g" value="0.8"/>
-		</phase>
-	</bsdf>
 
 	<!-- Test the smooth diffuse model -->
 	<bsdf type="diffuse"/>
@@ -158,6 +140,18 @@
 		<bsdf type="diffuse"/>
 	</bsdf>
 
+	<!-- Test the Hanrahan-Krueger model with an
+		 isotropic phase function -->
+	<bsdf type="hk"/>
+	
+	<!-- Test the Hanrahan-Krueger model with a 
+		 forward-scattering phase function -->
+	<bsdf type="hk">
+		<phase type="hg">
+			<float name="g" value="0.8"/>
+		</phase>
+	</bsdf>
+
 	<!-- Test the smooth coating over a diffuse base material -->
 	<bsdf type="coating">
 		<float name="intIOR" value="1.5"/>
@@ -200,4 +194,13 @@
 			<spectrum name="opacity" value="0.5"/>
 		</bsdf>
 	</bsdf>
+
+	<!-- Test the coating model with the Hanrahan-Krueger model -->
+	<bsdf type="coating">
+		<bsdf type="hk">
+			<rgb name="sigmaA" value="0.1 0.2 0.3"/>
+			<rgb name="sigmaS" value="1"/>
+			<float name="thickness" value="2"/>
+		</bsdf>
+	</bsdf>
 </scene>

src/bsdfs/SConscript

@@ -2,26 +2,26 @@ Import('env', 'plugins')
 
 # Basic library of smooth and rough materials
 plugins += env.SharedLibrary('diffuse', ['diffuse.cpp'])
-#plugins += env.SharedLibrary('dielectric', ['dielectric.cpp'])
+plugins += env.SharedLibrary('dielectric', ['dielectric.cpp'])
-#plugins += env.SharedLibrary('conductor', ['conductor.cpp'])
+plugins += env.SharedLibrary('conductor', ['conductor.cpp'])
-#plugins += env.SharedLibrary('plastic', ['plastic.cpp'])
+plugins += env.SharedLibrary('plastic', ['plastic.cpp'])
-#plugins += env.SharedLibrary('roughdiffuse', ['roughdiffuse.cpp'])
+plugins += env.SharedLibrary('roughdiffuse', ['roughdiffuse.cpp'])
-#plugins += env.SharedLibrary('roughdielectric', ['roughdielectric.cpp'])
+plugins += env.SharedLibrary('roughdielectric', ['roughdielectric.cpp'])
-#plugins += env.SharedLibrary('roughconductor', ['roughconductor.cpp'])
+plugins += env.SharedLibrary('roughconductor', ['roughconductor.cpp'])
-#plugins += env.SharedLibrary('roughplastic', ['roughplastic.cpp'])
+plugins += env.SharedLibrary('roughplastic', ['roughplastic.cpp'])
 
 # Materials that act as modifiers
-#plugins += env.SharedLibrary('twosided', ['twosided.cpp'])
+plugins += env.SharedLibrary('twosided', ['twosided.cpp'])
-#plugins += env.SharedLibrary('mask', ['mask.cpp'])
+plugins += env.SharedLibrary('mask', ['mask.cpp'])
-#plugins += env.SharedLibrary('mixture', ['mixture.cpp'])
+plugins += env.SharedLibrary('mixture', ['mixture.cpp'])
-#plugins += env.SharedLibrary('coating', ['coating.cpp'])
+plugins += env.SharedLibrary('coating', ['coating.cpp'])
-#plugins += env.SharedLibrary('bump', ['bump.cpp'])
+plugins += env.SharedLibrary('bump', ['bump.cpp'])
 
 # Other materials
-#plugins += env.SharedLibrary('ward', ['ward.cpp'])
+plugins += env.SharedLibrary('ward', ['ward.cpp'])
-#plugins += env.SharedLibrary('phong', ['phong.cpp'])
+plugins += env.SharedLibrary('phong', ['phong.cpp'])
-#plugins += env.SharedLibrary('irawan', ['irawan.cpp'])
+plugins += env.SharedLibrary('irawan', ['irawan.cpp'])
-#plugins += env.SharedLibrary('difftrans', ['difftrans.cpp'])
+plugins += env.SharedLibrary('difftrans', ['difftrans.cpp'])
-#plugins += env.SharedLibrary('hk', ['hk.cpp'])
+plugins += env.SharedLibrary('hk', ['hk.cpp'])
 
 Export('plugins')

src/bsdfs/bump.cpp

@@ -191,28 +191,6 @@ public:
 		return m_nested->pdf(perturbedQuery, measure);
 	}
 
-	Spectrum sample(BSDFQueryRecord &bRec, Float &pdf, const Point2 &sample) const {
-		const Intersection& its = bRec.its;
-		Intersection perturbed;
-		perturbIntersection(its, perturbed);
-
-		BSDFQueryRecord perturbedQuery(perturbed, bRec.sampler, bRec.quantity);
-		perturbedQuery.wi = perturbed.toLocal(its.toWorld(bRec.wi));
-		perturbedQuery.typeMask = bRec.typeMask;
-		perturbedQuery.component = bRec.component;
-		Spectrum result = m_nested->sample(perturbedQuery, pdf, sample);
-
-		if (!result.isZero()) {
-			bRec.sampledComponent = perturbedQuery.sampledComponent;
-			bRec.sampledType = perturbedQuery.sampledType;
-			bRec.wo = its.toLocal(perturbed.toWorld(perturbedQuery.wo));
-			if (Frame::cosTheta(bRec.wo) * Frame::cosTheta(perturbedQuery.wo) <= 0)
-				return Spectrum(0.0f);
-		}
-
-		return result;
-	}
-
 	Spectrum sample(BSDFQueryRecord &bRec, const Point2 &sample) const {
 		const Intersection& its = bRec.its;
 		Intersection perturbed;
@@ -234,6 +212,28 @@ public:
 		return result;
 	}
 
+	Spectrum sampleXXX(BSDFQueryRecord &bRec, Float &pdf, const Point2 &sample) const {
+		const Intersection& its = bRec.its;
+		Intersection perturbed;
+		perturbIntersection(its, perturbed);
+
+		BSDFQueryRecord perturbedQuery(perturbed, bRec.sampler, bRec.quantity);
+		perturbedQuery.wi = perturbed.toLocal(its.toWorld(bRec.wi));
+		perturbedQuery.typeMask = bRec.typeMask;
+		perturbedQuery.component = bRec.component;
+		Spectrum result = m_nested->sampleXXX(perturbedQuery, pdf, sample);
+
+		if (!result.isZero()) {
+			bRec.sampledComponent = perturbedQuery.sampledComponent;
+			bRec.sampledType = perturbedQuery.sampledType;
+			bRec.wo = its.toLocal(perturbed.toWorld(perturbedQuery.wo));
+			if (Frame::cosTheta(bRec.wo) * Frame::cosTheta(perturbedQuery.wo) <= 0)
+				return Spectrum(0.0f);
+		}
+
+		return result;
+	}
+
 	Shader *createShader(Renderer *renderer) const;
 
 	std::string toString() const {

src/bsdfs/coating.cpp

@@ -307,7 +307,7 @@ public:
 		}
 	}
 
-	Spectrum sample(BSDFQueryRecord &bRec, Float &pdf, const Point2 &_sample) const {
+	Spectrum sampleXXX(BSDFQueryRecord &bRec, Float &pdf, const Point2 &_sample) const {
 		bool sampleSpecular = (bRec.typeMask & EDeltaReflection)
 			&& (bRec.component == -1 || bRec.component == (int) m_components.size()-1);
 		bool sampleNested = (bRec.typeMask & m_nested->getType() & BSDF::EAll)
@@ -341,14 +341,14 @@ public:
 			bRec.sampledType = EDeltaReflection;
 			bRec.wo = reflect(bRec.wi);
 			pdf = sampleNested ? probSpecular : 1.0f;
-			return Spectrum(R12);
+			return Spectrum(R12) / pdf;
 		} else {
 			if (R12 == 1.0f) /* Total internal reflection */
 				return Spectrum(0.0f);
 
 			Vector wiBackup = bRec.wi;
 			bRec.wi = wiPrime;
-			Spectrum result = m_nested->sample(bRec, pdf, sample);
+			Spectrum result = m_nested->sampleXXX(bRec, pdf, sample);
 			bRec.wi = wiBackup;
 			if (result.isZero()) 
 				return Spectrum(0.0f);
@@ -366,34 +366,28 @@ public:
 			if (R21 == 1.0f) /* Total internal reflection */
 				return Spectrum(0.0f);
 
-			if (sampleSpecular)
+			if (sampleSpecular) {
-				pdf *= 1 - probSpecular;
+				pdf *= 1.0f - probSpecular;
+				result /= 1.0f - probSpecular;
+			}
 
 			result *= (1 - R12) * (1 - R21);
 
 			if (BSDF::getMeasure(bRec.sampledType) == ESolidAngle) {
 				/* Solid angle compression & irradiance conversion factors */
 				Float eta = m_extIOR / m_intIOR, etaSqr = eta*eta;
-				Float temp = Frame::cosTheta(bRec.wo) / Frame::cosTheta(woPrime);
 
-				result *= etaSqr *
+				result *= Frame::cosTheta(bRec.wi) / Frame::cosTheta(wiPrime);
-					Frame::cosTheta(bRec.wi) / Frame::cosTheta(wiPrime) * temp;
+				pdf *= etaSqr * Frame::cosTheta(bRec.wo) / Frame::cosTheta(woPrime);
-				pdf *= etaSqr * temp;
 			}
 
-
 			return result;
 		}
 	}
 
 	Spectrum sample(BSDFQueryRecord &bRec, const Point2 &sample) const {
 		Float pdf;
-		Spectrum result = SmoothCoating::sample(bRec, pdf, sample);
+		return SmoothCoating::sampleXXX(bRec, pdf, sample);
-
-		if (result.isZero())
-			return Spectrum(0.0f);
-		else
-			return result / pdf;
 	}
 
 	Shader *createShader(Renderer *renderer) const;

src/bsdfs/conductor.cpp

@@ -255,7 +255,7 @@ public:
 			fresnelConductor(Frame::cosTheta(bRec.wi), m_eta, m_k);
 	}
 
-	Spectrum sample(BSDFQueryRecord &bRec, Float &pdf, const Point2 &sample) const {
+	Spectrum sampleXXX(BSDFQueryRecord &bRec, Float &pdf, const Point2 &sample) const {
 		bool sampleReflection   = (bRec.typeMask & EDeltaReflection)
 				&& (bRec.component == -1 || bRec.component == 0);
 		

src/bsdfs/dielectric.cpp

@@ -387,7 +387,7 @@ public:
 		}
 	}
 
-	Spectrum sample(BSDFQueryRecord &bRec, Float &pdf, const Point2 &sample) const {
+	Spectrum sampleXXX(BSDFQueryRecord &bRec, Float &pdf, const Point2 &sample) const {
 		bool sampleReflection   = (bRec.typeMask & EDeltaReflection)
 				&& (bRec.component == -1 || bRec.component == 0);
 		bool sampleTransmission = (bRec.typeMask & EDeltaTransmission)
@@ -433,7 +433,7 @@ public:
 				bRec.wo = reflect(bRec.wi);
 
 				pdf = Fr;
-				return m_specularReflectance->getValue(bRec.its) * Fr;
+				return m_specularReflectance->getValue(bRec.its);
 			} else {
 				bRec.sampledComponent = 1;
 				bRec.sampledType = EDeltaTransmission;
@@ -447,7 +447,7 @@ public:
 				/* When transporting radiance, account for the solid angle
 				   change at boundaries with different indices of refraction. */
 				return m_specularTransmittance->getValue(bRec.its) 
-					* (1-Fr) * (bRec.quantity == ERadiance ? (eta*eta) : (Float) 1);
+					* (bRec.quantity == ERadiance ? (eta*eta) : (Float) 1);
 			}
 		} else if (sampleReflection) {
 			bRec.sampledComponent = 0;

src/bsdfs/difftrans.cpp

@@ -99,7 +99,7 @@ public:
 		return m_transmittance->getValue(bRec.its);
 	}
 
-	Spectrum sample(BSDFQueryRecord &bRec, Float &pdf, const Point2 &sample) const {
+	Spectrum sampleXXX(BSDFQueryRecord &bRec, Float &pdf, const Point2 &sample) const {
 		if (!(bRec.typeMask & m_combinedType)) 
 			return Spectrum(0.0f);
 		bRec.wo = squareToHemispherePSA(sample);
@@ -108,8 +108,7 @@ public:
 		bRec.sampledComponent = 0;
 		bRec.sampledType = EDiffuseTransmission;
 		pdf = std::abs(Frame::cosTheta(bRec.wo)) * INV_PI;
-		return m_transmittance->getValue(bRec.its) 
+		return m_transmittance->getValue(bRec.its);
-			* (INV_PI * std::abs(Frame::cosTheta(bRec.wo)));
 	}
 
 	void addChild(const std::string &name, ConfigurableObject *child) {

src/bsdfs/hk.cpp

@@ -272,7 +272,7 @@ public:
 		return 0.0f;
 	}
 
-	inline Spectrum sample(BSDFQueryRecord &bRec, Float &_pdf, const Point2 &_sample) const {
+	inline Spectrum sampleXXX(BSDFQueryRecord &bRec, Float &_pdf, const Point2 &_sample) const {
 		AssertEx(bRec.sampler != NULL, "The BSDFQueryRecord needs to have a sampler!");
 
 		bool hasSpecularTransmission = (bRec.typeMask & EDeltaTransmission)
@@ -307,7 +307,7 @@ public:
 			bRec.wo = -bRec.wi;
 
 			_pdf = hasSingleScattering ? probSpecularTransmission : 1.0f;
-			return eval(bRec, EDiscrete);
+			return eval(bRec, EDiscrete) / _pdf;
 		} else {
 			/* The glossy transmission/scattering component should be sampled */
 			bool hasGlossyReflection = (bRec.typeMask & EGlossyReflection)
@@ -317,7 +317,7 @@ public:
 
 			/* Sample According to the phase function lobes */
 			PhaseFunctionQueryRecord pRec(MediumSamplingRecord(), bRec.wi, bRec.wo);
-			m_phase->sample(pRec, _pdf, bRec.sampler);
+			m_phase->sampleXXX(pRec, _pdf, bRec.sampler);
 
 			/* Store the sampled direction */
 			bRec.wo = pRec.wo;
@@ -337,19 +337,14 @@ public:
 			if (_pdf == 0) 
 				return Spectrum(0.0f);
 			else
-				return eval(bRec, ESolidAngle);
+				return eval(bRec, ESolidAngle) / _pdf;
 
 		}
 	}
 
 	Spectrum sample(BSDFQueryRecord &bRec, const Point2 &sample) const {
-		Float pdf = 0;
+		Float pdf;
-		Spectrum result = HanrahanKrueger::sample(bRec, pdf, sample);
+		return HanrahanKrueger::sampleXXX(bRec, pdf, sample);
-
-		if (result.isZero())
-			return Spectrum(0.0f);
-		else
-			return result / pdf;
 	}
 
 	void serialize(Stream *stream, InstanceManager *manager) const {

src/bsdfs/irawan.cpp

@@ -303,10 +303,10 @@ public:
 		bRec.wo = squareToHemispherePSA(sample);
 		bRec.sampledComponent = 0;
 		bRec.sampledType = EGlossyReflection;
-		return eval(bRec, ESolidAngle) / Frame::cosTheta(bRec.wo);
+		return eval(bRec, ESolidAngle) * M_PI / Frame::cosTheta(bRec.wo);
 	}
 
-	Spectrum sample(BSDFQueryRecord &bRec, Float &pdf, const Point2 &sample) const {
+	Spectrum sampleXXX(BSDFQueryRecord &bRec, Float &pdf, const Point2 &sample) const {
 		bool hasSpecular = (bRec.typeMask & EGlossyReflection) &&
 			(bRec.component == -1 || bRec.component == 0) && m_ksMultiplier > 0;
 		bool hasDiffuse = (bRec.typeMask & EDiffuseReflection) &&
@@ -321,7 +321,7 @@ public:
 		bRec.sampledComponent = 0;
 		bRec.sampledType = EGlossyReflection;
 		pdf = Frame::cosTheta(bRec.wo) * INV_PI;
-		return eval(bRec, ESolidAngle);
+		return eval(bRec, ESolidAngle) / pdf;
 	}
 
 	void serialize(Stream *stream, InstanceManager *manager) const {

src/bsdfs/mask.cpp

@@ -140,42 +140,6 @@ public:
 		}
 	}
 
-	Spectrum sample(BSDFQueryRecord &bRec, Float &pdf, const Point2 &_sample) const {
-		Point2 sample(_sample);
-		Spectrum result(0.0f);
-
-		Spectrum opacity = m_opacity->getValue(bRec.its);
-		Float prob = opacity.getLuminance();
-
-		bool sampleTransmission = bRec.typeMask & EDeltaTransmission
-			&& (bRec.component == -1 || bRec.component == getComponentCount()-1);
-		bool sampleNested = bRec.component == -1 || bRec.component < getComponentCount()-1;
-
-		if (sampleTransmission && sampleNested) {
-			if (sample.x <= prob) {
-				sample.x /= prob;
-				result = m_nestedBSDF->sample(bRec, pdf, sample) * opacity;
-				pdf *= prob;
-			} else {
-				bRec.wo = -bRec.wi;
-				bRec.sampledComponent = getComponentCount()-1;
-				bRec.sampledType = EDeltaTransmission;
-				pdf = 1-prob;
-				result = Spectrum(1.0f) - opacity;
-			}
-		} else if (sampleTransmission) {
-			bRec.wo = -bRec.wi;
-			bRec.sampledComponent = getComponentCount()-1;
-			bRec.sampledType = EDeltaTransmission;
-			pdf = 1;
-			result = Spectrum(1.0f) - opacity;
-		} else if (sampleNested) {
-			result = m_nestedBSDF->sample(bRec, pdf, sample) * opacity;
-		}
-
-		return result;
-	}
-
 	Spectrum sample(BSDFQueryRecord &bRec, const Point2 &_sample) const {
 		Point2 sample(_sample);
 		Spectrum opacity = m_opacity->getValue(bRec.its);
@@ -208,6 +172,42 @@ public:
 		}
 	}
 
+	Spectrum sampleXXX(BSDFQueryRecord &bRec, Float &pdf, const Point2 &_sample) const {
+		Point2 sample(_sample);
+		Spectrum result(0.0f);
+
+		Spectrum opacity = m_opacity->getValue(bRec.its);
+		Float prob = opacity.getLuminance();
+
+		bool sampleTransmission = bRec.typeMask & EDeltaTransmission
+			&& (bRec.component == -1 || bRec.component == getComponentCount()-1);
+		bool sampleNested = bRec.component == -1 || bRec.component < getComponentCount()-1;
+
+		if (sampleTransmission && sampleNested) {
+			if (sample.x <= prob) {
+				sample.x /= prob;
+				result = m_nestedBSDF->sampleXXX(bRec, pdf, sample) * opacity / prob;
+				pdf *= prob;
+			} else {
+				bRec.wo = -bRec.wi;
+				bRec.sampledComponent = getComponentCount()-1;
+				bRec.sampledType = EDeltaTransmission;
+				pdf = 1-prob;
+				result = (Spectrum(1.0f) - opacity) / pdf;
+			}
+		} else if (sampleTransmission) {
+			bRec.wo = -bRec.wi;
+			bRec.sampledComponent = getComponentCount()-1;
+			bRec.sampledType = EDeltaTransmission;
+			pdf = 1;
+			result = Spectrum(1.0f) - opacity;
+		} else if (sampleNested) {
+			result = m_nestedBSDF->sampleXXX(bRec, pdf, sample) * opacity;
+		}
+
+		return result;
+	}
+
 	void addChild(const std::string &name, ConfigurableObject *child) {
 		if (child->getClass()->derivesFrom(MTS_CLASS(Texture)) && name == "opacity") 
 			m_opacity = static_cast<Texture *>(child);

src/bsdfs/microfacet.h

@@ -240,50 +240,138 @@ public:
 	Normal sample(const Point2 &sample, Float alphaU, Float alphaV) const {
 		/* The azimuthal component is always selected 
 		   uniformly regardless of the distribution */
-		Float phiM = (2.0f * M_PI) * sample.y,
+		Float cosThetaM = 0.0f, phiM = (2.0f * M_PI) * sample.y;
-			  thetaM = 0.0f;
 
 		switch (m_type) {
-			case EBeckmann: 
+			case EBeckmann: {
-				thetaM = std::atan(std::sqrt(-alphaU*alphaU *
+					Float tanThetaMSqr = -alphaU*alphaU * std::log(1.0f - sample.x);
-						 std::log(1.0f - sample.x)));
+					cosThetaM = 1.0f / std::sqrt(1 + tanThetaMSqr);
+				}
 				break;
 
-			case EGGX: 
+			case EGGX: {
-				thetaM = std::atan(alphaU * std::sqrt(sample.x) /
+					Float tanThetaMSqr = alphaU * alphaU * sample.x / (1.0f - sample.x);
-						 std::sqrt(1.0f - sample.x));
+					cosThetaM = 1.0f / std::sqrt(1 + tanThetaMSqr);
+				}
 				break;
 
-			case EPhong:
+			case EPhong: {
-				thetaM = std::acos(std::pow(sample.x, (Float) 1 / 
+					cosThetaM = std::pow(sample.x, 1 / (alphaU + 2));
-						 (alphaU + 2)));
+				}
 				break;
 
 			case EAshikhminShirley: {
 					/* Sampling method based on code from PBRT */
-					Float phi, cosTheta;
 					if (sample.x < 0.25f) {
 						sampleFirstQuadrant(alphaU, alphaV,
-							4 * sample.x, sample.y, phi, cosTheta);
+							4 * sample.x, sample.y, phiM, cosThetaM);
 					} else if (sample.x < 0.5f) {
 						sampleFirstQuadrant(alphaU, alphaV,
-							4 * (0.5f - sample.x), sample.y, phi, cosTheta);
+							4 * (0.5f - sample.x), sample.y, phiM, cosThetaM);
-						phi = M_PI - phi;
+						phiM = M_PI - phiM;
 					} else if (sample.x < 0.75f) {
 						sampleFirstQuadrant(alphaU, alphaV,
-							4 * (sample.x - 0.5f), sample.y, phi, cosTheta);
+							4 * (sample.x - 0.5f), sample.y, phiM, cosThetaM);
-						phi += M_PI;
+						phiM += M_PI;
 					} else {
 						sampleFirstQuadrant(alphaU, alphaV,
-							4 * (1 - sample.x), sample.y, phi, cosTheta);
+							4 * (1 - sample.x), sample.y, phiM, cosThetaM);
-						phi = 2 * M_PI - phi;
+						phiM = 2 * M_PI - phiM;
 					}
-					const Float sinTheta = std::sqrt(
+				}
-						std::max((Float) 0, 1 - cosTheta*cosTheta));
+				break;
+			default: 
+				SLog(EError, "Invalid distribution function!");
+		}
+
+		const Float sinThetaM = std::sqrt(
+			std::max((Float) 0, 1 - cosThetaM*cosThetaM));
+		return Vector(
+			sinThetaM * std::cos(phiM),
+			sinThetaM * std::sin(phiM),
+			cosThetaM
+		);
+	}
+
+
+	/**
+	 * \brief Draw a sample from the microsurface normal distribution
+	 * and return the associated probability density
+	 *
+	 * \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
+	 * \param pdf     The probability density wrt. solid angles
+	 */
+	Normal sample(const Point2 &sample, Float alphaU, Float alphaV, Float &pdf) const {
+		/* The azimuthal component is always selected 
+		   uniformly regardless of the distribution */
+		Float cosThetaM = 0.0f;
+
+		switch (m_type) {
+			case EBeckmann: {
+					Float tanThetaMSqr = -alphaU*alphaU * std::log(1.0f - sample.x);
+					cosThetaM = 1.0f / std::sqrt(1 + tanThetaMSqr);
+					Float cosThetaM2 = cosThetaM * cosThetaM,
+						  cosThetaM3 = cosThetaM2 * cosThetaM;
+					pdf = (1.0f - sample.x) / (M_PI * alphaU*alphaU * cosThetaM3);
+				}
+				break;
+
+			case EGGX: {
+					Float alphaUSqr = alphaU * alphaU;
+					Float tanThetaMSqr = alphaUSqr * sample.x / (1.0f - sample.x);
+					cosThetaM = 1.0f / std::sqrt(1 + tanThetaMSqr);
+
+					Float cosThetaM2 = cosThetaM * cosThetaM,
+						  cosThetaM3 = cosThetaM2 * cosThetaM,
+						  temp = alphaUSqr + tanThetaMSqr;
+
+					pdf = INV_PI * alphaUSqr / (cosThetaM3 * temp * temp);
+				}
+				break;
+
+			case EPhong: {
+					Float exponent = 1 / (alphaU + 2);
+					cosThetaM = std::pow(sample.x, exponent);
+					pdf = (alphaU + 2) * INV_TWOPI * std::pow(sample.x, (alphaU+1) * exponent);
+				}
+				break;
+
+			case EAshikhminShirley: {
+					Float phiM;
+
+					/* Sampling method based on code from PBRT */
+					if (sample.x < 0.25f) {
+						sampleFirstQuadrant(alphaU, alphaV,
+							4 * sample.x, sample.y, phiM, cosThetaM);
+					} else if (sample.x < 0.5f) {
+						sampleFirstQuadrant(alphaU, alphaV,
+							4 * (0.5f - sample.x), sample.y, phiM, cosThetaM);
+						phiM = M_PI - phiM;
+					} else if (sample.x < 0.75f) {
+						sampleFirstQuadrant(alphaU, alphaV,
+							4 * (sample.x - 0.5f), sample.y, phiM, cosThetaM);
+						phiM += M_PI;
+					} else {
+						sampleFirstQuadrant(alphaU, alphaV,
+							4 * (1 - sample.x), sample.y, phiM, cosThetaM);
+						phiM = 2 * M_PI - phiM;
+					}
+					const Float sinThetaM = std::sqrt(
+						    std::max((Float) 0, 1 - cosThetaM*cosThetaM)),
+						  sinPhiM = std::sin(phiM),
+						  cosPhiM = std::cos(phiM);
+
+					const Float exponent = alphaU * cosPhiM*cosPhiM
+							+ alphaV * sinPhiM*sinPhiM;
+					pdf = std::sqrt((alphaU + 1) * (alphaV + 1))
+						* INV_TWOPI * std::pow(cosThetaM, exponent);
+					
 					return Vector(
-						sinTheta * std::cos(phi),
+						sinThetaM * cosPhiM,
-						sinTheta * std::sin(phi),
+						sinThetaM * sinPhiM,
-						cosTheta
+						cosThetaM
 					);
 				}
 				break;
@@ -291,33 +379,14 @@ public:
 				SLog(EError, "Invalid distribution function!");
 		}
 
-		return Normal(sphericalDirection(thetaM, phiM));
+		const Float sinThetaM = std::sqrt(
-	}
+			std::max((Float) 0, 1 - cosThetaM*cosThetaM));
-
+		Float phiM = (2.0f * M_PI) * sample.y;
-	/**
+		return Vector(
-	 * \brief Draw a sample from an isotropic microsurface normal
+			sinThetaM * std::cos(phiM),
-	 * distribution and return the magnitude of its 'z' component.
+			sinThetaM * std::sin(phiM),
-	 *
+			cosThetaM
-	 * \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;
-		}
 	}
 
 	/**

src/bsdfs/mixture.cpp

@@ -204,41 +204,6 @@ public:
 		return result;
 	}
 
-	Spectrum sample(BSDFQueryRecord &bRec, Float &pdf, const Point2 &_sample) const {
-		Point2 sample(_sample);
-		if (bRec.component == -1) {
-			/* Choose a component based on the normalized weights */
-			size_t entry = m_pdf.sampleReuse(sample.x);
-
-			Spectrum result = m_bsdfs[entry]->sample(bRec, pdf, sample);
-			if (result.isZero()) // sampling failed
-				return result;
-
-			result *= m_weights[entry];
-			pdf *= m_pdf[entry];
-
-			EMeasure measure = BSDF::getMeasure(bRec.sampledType);
-			for (size_t i=0; i<m_bsdfs.size(); ++i) {
-				if (entry == i)
-					continue;
-				pdf += m_bsdfs[i]->pdf(bRec, measure) * m_pdf[i];
-				result += m_bsdfs[i]->eval(bRec, measure) * m_weights[i];
-			}
-
-			bRec.sampledComponent += m_offsets[entry];
-			return result;
-		} else {
-			/* Pick out an individual component */
-			int requestedComponent = bRec.component;
-			int bsdfIndex = m_indices[requestedComponent].first;
-			bRec.component = m_indices[requestedComponent].second;
-			Spectrum result = m_bsdfs[bsdfIndex]->sample(bRec, pdf, sample)
-				* m_weights[bsdfIndex];
-			bRec.component = bRec.sampledComponent = requestedComponent;
-			return result;
-		}
-	}
-
 	Spectrum sample(BSDFQueryRecord &bRec, const Point2 &_sample) const {
 		Point2 sample(_sample);
 		if (bRec.component == -1) {
@@ -246,11 +211,11 @@ public:
 			size_t entry = m_pdf.sampleReuse(sample.x);
 
 			Float pdf;
-			Spectrum result = m_bsdfs[entry]->sample(bRec, pdf, sample);
+			Spectrum result = m_bsdfs[entry]->sampleXXX(bRec, pdf, sample);
 			if (result.isZero()) // sampling failed
 				return result;
 
-			result *= m_weights[entry];
+			result *= m_weights[entry] * pdf;
 			pdf *= m_pdf[entry];
 
 			EMeasure measure = BSDF::getMeasure(bRec.sampledType);
@@ -275,6 +240,41 @@ public:
 		}
 	}
 
+	Spectrum sampleXXX(BSDFQueryRecord &bRec, Float &pdf, const Point2 &_sample) const {
+		Point2 sample(_sample);
+		if (bRec.component == -1) {
+			/* Choose a component based on the normalized weights */
+			size_t entry = m_pdf.sampleReuse(sample.x);
+
+			Spectrum result = m_bsdfs[entry]->sampleXXX(bRec, pdf, sample);
+			if (result.isZero()) // sampling failed
+				return result;
+
+			result *= m_weights[entry] * pdf;
+			pdf *= m_pdf[entry];
+			
+			EMeasure measure = BSDF::getMeasure(bRec.sampledType);
+			for (size_t i=0; i<m_bsdfs.size(); ++i) {
+				if (entry == i)
+					continue;
+				pdf += m_bsdfs[i]->pdf(bRec, measure) * m_pdf[i];
+				result += m_bsdfs[i]->eval(bRec, measure) * m_weights[i];
+			}
+
+			bRec.sampledComponent += m_offsets[entry];
+			return result/pdf;
+		} else {
+			/* Pick out an individual component */
+			int requestedComponent = bRec.component;
+			int bsdfIndex = m_indices[requestedComponent].first;
+			bRec.component = m_indices[requestedComponent].second;
+			Spectrum result = m_bsdfs[bsdfIndex]->sampleXXX(bRec, pdf, sample)
+				* m_weights[bsdfIndex];
+			bRec.component = bRec.sampledComponent = requestedComponent;
+			return result;
+		}
+	}
+
 	void addChild(const std::string &name, ConfigurableObject *child) {
 		if (child->getClass()->derivesFrom(MTS_CLASS(BSDF))) {
 			BSDF *bsdf = static_cast<BSDF *>(child);

src/bsdfs/phong.cpp

@@ -175,7 +175,7 @@ public:
 			return 0.0f;
 	}
 
-	inline Spectrum sample(BSDFQueryRecord &bRec, Float &_pdf, const Point2 &_sample) const {
+	inline Spectrum sampleXXX(BSDFQueryRecord &bRec, Float &_pdf, const Point2 &_sample) const {
 		Point2 sample(_sample);
 
 		bool hasSpecular = (bRec.typeMask & EGlossyReflection)
@@ -230,17 +230,12 @@ public:
 		if (_pdf == 0) 
 			return Spectrum(0.0f);
 		else
-			return eval(bRec, ESolidAngle);
+			return eval(bRec, ESolidAngle) / _pdf;
 	}
 
 	Spectrum sample(BSDFQueryRecord &bRec, const Point2 &sample) const {
-		Float pdf = 0;
+		Float pdf;
-		Spectrum result = Phong::sample(bRec, pdf, sample);
+		return Phong::sampleXXX(bRec, pdf, sample);
-
-		if (result.isZero())
-			return Spectrum(0.0f);
-		else
-			return result / pdf;
 	}
 
 	void addChild(const std::string &name, ConfigurableObject *child) {

src/bsdfs/plastic.cpp

@@ -182,17 +182,17 @@ public:
 
 		Float probSpecular = 1.0f;
 		if (hasSpecular && hasDiffuse) {
-			probSpecular = fresnel(Frame::cosTheta(bRec.wi), m_extIOR, m_intIOR);
+			Float Fr = fresnel(Frame::cosTheta(bRec.wi), m_extIOR, m_intIOR);
-			probSpecular = (probSpecular*m_specularSamplingWeight) /
+			probSpecular = (Fr*m_specularSamplingWeight) /
-				(probSpecular*m_specularSamplingWeight + 
+				(Fr*m_specularSamplingWeight + 
-				(1-probSpecular) * (1-m_specularSamplingWeight));
+				(1-Fr) * (1-m_specularSamplingWeight));
 		}
 
 		if (measure == EDiscrete && hasSpecular) {
 			/* Check if the provided direction pair matches an ideal
 			   specular reflection; tolerate some roundoff errors */
 			if (std::abs(1 - dot(reflect(bRec.wi), bRec.wo)) < Epsilon)
-				return hasDiffuse ? probSpecular : 1.0f;
+				return probSpecular;
 		} else if (measure == ESolidAngle && hasDiffuse) {
 			return Frame::cosTheta(bRec.wo) * INV_PI *
 				(hasSpecular ? (1 - probSpecular) : 1.0f);
@@ -224,7 +224,7 @@ public:
 				bRec.wo = reflect(bRec.wi);
 
 				return m_specularReflectance->getValue(bRec.its) *
-					Fr / probSpecular;
+					(Fr / probSpecular);
 			} else {
 				bRec.sampledComponent = 1;
 				bRec.sampledType = EDiffuseReflection;
@@ -234,7 +234,7 @@ public:
 				));
 
 				return m_diffuseReflectance->getValue(bRec.its) *
-					(1-Fr) / (1-probSpecular);
+					((1-Fr) / (1-probSpecular));
 			}
 		} else if (hasSpecular) {
 			bRec.sampledComponent = 0;
@@ -254,7 +254,7 @@ public:
 		}
 	}
 
-	Spectrum sample(BSDFQueryRecord &bRec, Float &pdf, const Point2 &sample) const {
+	Spectrum sampleXXX(BSDFQueryRecord &bRec, Float &pdf, const Point2 &sample) const {
 		bool hasSpecular   = (bRec.typeMask & EDeltaReflection)
 				&& (bRec.component == -1 || bRec.component == 0);
 		bool hasDiffuse = (bRec.typeMask & EDiffuseReflection)
@@ -276,7 +276,8 @@ public:
 				bRec.wo = reflect(bRec.wi);
 
 				pdf = probSpecular;
-				return m_specularReflectance->getValue(bRec.its) * Fr;
+				return m_specularReflectance->getValue(bRec.its)
+					* Fr / probSpecular;
 			} else {
 				bRec.sampledComponent = 1;
 				bRec.sampledType = EDiffuseReflection;
@@ -288,7 +289,7 @@ public:
 				pdf = (1-probSpecular) * Frame::cosTheta(bRec.wo) * INV_PI;
 
 				return m_diffuseReflectance->getValue(bRec.its) 
-					* (INV_PI * Frame::cosTheta(bRec.wo) * (1-Fr));
+					* (1-Fr) / (1-probSpecular);
 			}
 		} else if (hasSpecular) {
 			bRec.sampledComponent = 0;
@@ -307,8 +308,7 @@ public:
 
 			pdf = Frame::cosTheta(bRec.wo) * INV_PI;
 
-			return m_diffuseReflectance->getValue(bRec.its) 
+			return m_diffuseReflectance->getValue(bRec.its) * (1-Fr);
-				* (INV_PI * Frame::cosTheta(bRec.wo) * (1-Fr));
 		}
 	}
 

src/bsdfs/roughconductor.cpp

@@ -298,7 +298,9 @@ public:
 					m_alphaV->getValue(bRec.its).average());
 
 		/* Sample M, the microsurface normal */
-		const Normal m = m_distribution.sample(sample, alphaU, alphaV);
+		Float microfacetPDF;
+		const Normal m = m_distribution.sample(sample, 
+			alphaU, alphaV, microfacetPDF);
 
 		/* Perfect specular reflection based on the microsurface normal */
 		bRec.wo = reflect(bRec.wi, m);
@@ -316,14 +318,14 @@ public:
 			* m_distribution.G(bRec.wi, bRec.wo, m, alphaU, alphaV)
 			* dot(bRec.wi, m);
 
-		Float denominator = m_distribution.pdf(m, alphaU, alphaV)
+		Float denominator = microfacetPDF
 			* Frame::cosTheta(bRec.wi);
 
 		return m_specularReflectance->getValue(bRec.its) * F
 				* (numerator / denominator);
 	}
 
-	Spectrum sample(BSDFQueryRecord &bRec, Float &_pdf, const Point2 &sample) const {
+	Spectrum sampleXXX(BSDFQueryRecord &bRec, Float &pdf, const Point2 &sample) const {
 		if (Frame::cosTheta(bRec.wi) < 0 ||
 			((bRec.component != -1 && bRec.component != 0) ||
 			!(bRec.typeMask & EGlossyReflection)))
@@ -336,7 +338,8 @@ public:
 					m_alphaV->getValue(bRec.its).average());
 
 		/* Sample M, the microsurface normal */
-		const Normal m = m_distribution.sample(sample, alphaU, alphaV);
+		const Normal m = m_distribution.sample(sample, 
+			alphaU, alphaV, pdf);
 
 		/* Perfect specular reflection based on the microsurface normal */
 		bRec.wo = reflect(bRec.wi, m);
@@ -347,15 +350,23 @@ public:
 		if (Frame::cosTheta(bRec.wo) <= 0)
 			return Spectrum(0.0f);
 
-		/* Guard against numerical imprecisions */
+		const Spectrum F = fresnelConductor(Frame::cosTheta(bRec.wi),
-		_pdf = pdf(bRec, ESolidAngle);
+				m_eta, m_k);
 
-		if (_pdf == 0) 
+		Float numerator = m_distribution.eval(m, alphaU, alphaV)
-			return Spectrum(0.0f);
+			* m_distribution.G(bRec.wi, bRec.wo, m, alphaU, alphaV)
-		else
+			* dot(bRec.wi, m);
-			return eval(bRec, ESolidAngle);
+		
+		Float denominator = pdf * Frame::cosTheta(bRec.wi);
+
+		/* Jacobian of the half-direction transform */
+		pdf /= 4.0f * dot(bRec.wo, m);
+
+		return m_specularReflectance->getValue(bRec.its) * F
+				* (numerator / denominator);
 	}
 
+
 	void addChild(const std::string &name, ConfigurableObject *child) {
 		if (child->getClass()->derivesFrom(MTS_CLASS(Texture))) {
 			if (name == "alpha")

src/bsdfs/roughdielectric.cpp

@@ -367,8 +367,8 @@ public:
 							  && (bRec.typeMask & EGlossyReflection)),
 		     hasTransmission = ((bRec.component == -1 || bRec.component == 1)
 							  && (bRec.typeMask & EGlossyTransmission)),
-		     reflect            = Frame::cosTheta(bRec.wi) 
+		     reflect         = Frame::cosTheta(bRec.wi) 
-				                * Frame::cosTheta(bRec.wo) > 0;
+				             * Frame::cosTheta(bRec.wo) > 0;
 
 		/* Determine the appropriate indices of refraction */
 		Float etaI = m_extIOR, etaT = m_intIOR;
@@ -385,7 +385,7 @@ public:
 				return 0.0f;
 
 			/* Calculate the reflection half-vector (and possibly flip it
-			   so that it lies inside the hemisphere around the normal) */
+			   so that it lies inside the hemisphere around the surface normal) */
 			H = normalize(bRec.wo+bRec.wi) 
 				* signum(Frame::cosTheta(bRec.wo));
 
@@ -460,8 +460,9 @@ public:
 #endif
 
 		/* Sample M, the microsurface normal */
+		Float microfacetPDF;
 		const Normal m = m_distribution.sample(sample,
-				sampleAlphaU, sampleAlphaV);
+				sampleAlphaU, sampleAlphaV, microfacetPDF);
 
 		Float F = fresnel(dot(bRec.wi, m), m_extIOR, m_intIOR),
 			  numerator = 1.0f;
@@ -511,13 +512,13 @@ public:
 			* m_distribution.G(bRec.wi, bRec.wo, m, alphaU, alphaV)
 			* dot(bRec.wi, m);
 
-		Float denominator = m_distribution.pdf(m, sampleAlphaU, sampleAlphaV)
+		Float denominator = microfacetPDF
 			* Frame::cosTheta(bRec.wi);
 
 		return result * std::abs(numerator / denominator);
 	}
 
-	Spectrum sample(BSDFQueryRecord &bRec, Float &_pdf, const Point2 &_sample) const {
+	Spectrum sampleXXX(BSDFQueryRecord &bRec, Float &pdf, const Point2 &_sample) const {
 		Point2 sample(_sample);
 
 		bool hasReflection = ((bRec.component == -1 || bRec.component == 0)
@@ -546,15 +547,35 @@ public:
 #endif
 
 		/* Sample M, the microsurface normal */
+		Float microfacetPDF;
 		const Normal m = m_distribution.sample(sample,
-				sampleAlphaU, sampleAlphaV);
+				sampleAlphaU, sampleAlphaV, microfacetPDF);
+
+#if 1
+		Float ref = m_distribution.pdf(m, sampleAlphaU, sampleAlphaV);
+		if (std::abs(ref-microfacetPDF)/ref > Epsilon)
+			cout << "OOPS! ref=" << ref << ", got=" << microfacetPDF << endl;
+#endif
+
+		pdf = microfacetPDF;
+
+		Float F = fresnel(dot(bRec.wi, m), m_extIOR, m_intIOR),
+			  numerator = 1.0f;
 
 		if (hasReflection && hasTransmission) {
-			Float F = fresnel(dot(bRec.wi, m), m_extIOR, m_intIOR);
+			if (bRec.sampler->next1D() > F) {
-			if (bRec.sampler->next1D() > F)
 				choseReflection = false;
+				pdf *= (1-F);
+			} else {
+				pdf *= F;
+			}
+		} else {
+			numerator = hasReflection ? F : 1-F;
 		}
 
+		Spectrum result;
+		Float dwh_dwo;
+
 		if (choseReflection) {
 			/* Perfect specular reflection based on the microsurface normal */
 			bRec.wo = reflect(bRec.wi, m);
@@ -564,6 +585,11 @@ public:
 			/* Side check */
 			if (Frame::cosTheta(bRec.wi) * Frame::cosTheta(bRec.wo) <= 0)
 				return Spectrum(0.0f);
+		
+			result = m_specularReflectance->getValue(bRec.its);
+
+			/* Jacobian of the half-direction transform */
+			dwh_dwo = 1.0f / (4.0f * dot(bRec.wo, m));
 		} else {
 			/* Determine the appropriate indices of refraction */
 			Float etaI = m_extIOR, etaT = m_intIOR;
@@ -580,17 +606,28 @@ public:
 			/* Side check */
 			if (Frame::cosTheta(bRec.wi) * Frame::cosTheta(bRec.wo) >= 0)
 				return Spectrum(0.0f);
+
+			result = m_specularTransmittance->getValue(bRec.its)
+				* ((bRec.quantity == ERadiance) ?
+					((etaI*etaI) / (etaT*etaT)) : (Float) 1);
+
+			/* Jacobian of the half-direction transform. */
+			Float sqrtDenom = etaI * dot(bRec.wi, m) + etaT * dot(bRec.wo, m);
+			dwh_dwo = (etaT*etaT * dot(bRec.wo, m)) / (sqrtDenom*sqrtDenom);
 		}
 
-		/* Guard against numerical imprecisions */
+		numerator *= m_distribution.eval(m, alphaU, alphaV)
-		_pdf = pdf(bRec, ESolidAngle);
+			* m_distribution.G(bRec.wi, bRec.wo, m, alphaU, alphaV)
+			* dot(bRec.wi, m);
 
-		if (_pdf == 0) 
+		Float denominator = microfacetPDF * Frame::cosTheta(bRec.wi);
-			return Spectrum(0.0f);
+
-		else
+		pdf *= std::abs(dwh_dwo);
-			return eval(bRec, ESolidAngle);
+
+		return result * std::abs(numerator / denominator);
 	}
 
+
 	void addChild(const std::string &name, ConfigurableObject *child) {
 		if (child->getClass()->derivesFrom(MTS_CLASS(Texture))) {
 			if (name == "alpha")

src/bsdfs/roughdiffuse.cpp

@@ -236,7 +236,7 @@ public:
 			(Frame::cosTheta(bRec.wo) * INV_PI);
 	}
 
-	Spectrum sample(BSDFQueryRecord &bRec, Float &pdf, const Point2 &sample) const {
+	Spectrum sampleXXX(BSDFQueryRecord &bRec, Float &pdf, const Point2 &sample) const {
 		if (!(bRec.typeMask & EGlossyReflection) || Frame::cosTheta(bRec.wi) <= 0)
 			return Spectrum(0.0f);
 		
@@ -244,7 +244,7 @@ public:
 		bRec.sampledComponent = 0;
 		bRec.sampledType = EGlossyReflection;
 		pdf = Frame::cosTheta(bRec.wo) * INV_PI;
-		return eval(bRec, ESolidAngle);
+		return eval(bRec, ESolidAngle) / pdf;
 	}
 
 	void addChild(const std::string &name, ConfigurableObject *child) {

src/bsdfs/roughplastic.cpp

@@ -305,7 +305,7 @@ public:
 		return result;
 	}
 
-	inline Spectrum sample(BSDFQueryRecord &bRec, Float &_pdf, const Point2 &_sample) const {
+	inline Spectrum sampleXXX(BSDFQueryRecord &bRec, Float &_pdf, const Point2 &_sample) const {
 		bool hasSpecular = (bRec.typeMask & EGlossyReflection) &&
 			(bRec.component == -1 || bRec.component == 0);
 		bool hasDiffuse = (bRec.typeMask & EDiffuseReflection) &&
@@ -359,17 +359,12 @@ public:
 		if (_pdf == 0) 
 			return Spectrum(0.0f);
 		else
-			return eval(bRec, ESolidAngle);
+			return eval(bRec, ESolidAngle) / _pdf;
 	}
 
 	Spectrum sample(BSDFQueryRecord &bRec, const Point2 &sample) const {
-		Float pdf = 0;
+		Float pdf;
-		Spectrum result = RoughPlastic::sample(bRec, pdf, sample);
+		return RoughPlastic::sampleXXX(bRec, pdf, sample);
-
-		if (result.isZero())
-			return Spectrum(0.0f);
-		else
-			return result / pdf;
 	}
 
 	void serialize(Stream *stream, InstanceManager *manager) const {

src/bsdfs/twosided.cpp

@@ -123,14 +123,14 @@ public:
 		return result;
 	}
 
-	Spectrum sample(BSDFQueryRecord &bRec, Float &pdf, const Point2 &sample) const {
+	Spectrum sampleXXX(BSDFQueryRecord &bRec, Float &pdf, const Point2 &sample) const {
 		bool flipped = false;
 		if (Frame::cosTheta(bRec.wi) < 0) {
 			bRec.wi.z *= -1;
 			flipped = true;
 		}
 
-		Spectrum result = m_nestedBRDF->sample(bRec, pdf, sample);
+		Spectrum result = m_nestedBRDF->sampleXXX(bRec, pdf, sample);
 
 		if (flipped) {
 			bRec.wi.z *= -1;

src/bsdfs/ward.cpp

@@ -262,7 +262,7 @@ public:
 			return 0.0f;
 	}
 
-	inline Spectrum sample(BSDFQueryRecord &bRec, Float &_pdf, const Point2 &_sample) const {
+	inline Spectrum sampleXXX(BSDFQueryRecord &bRec, Float &_pdf, const Point2 &_sample) const {
 		Point2 sample(_sample);
 
 		bool hasSpecular = (bRec.typeMask & EGlossyReflection)
@@ -321,17 +321,12 @@ public:
 		if (_pdf == 0) 
 			return Spectrum(0.0f);
 		else
-			return eval(bRec, ESolidAngle);
+			return eval(bRec, ESolidAngle) / _pdf;
 	}
 
 	Spectrum sample(BSDFQueryRecord &bRec, const Point2 &sample) const {
-		Float pdf=0;
+		Float pdf;
-		Spectrum result = Ward::sample(bRec, pdf, sample);
+		return Ward::sampleXXX(bRec, pdf, sample);
-
-		if (result.isZero())
-			return Spectrum(0.0f);
-		else
-			return result / pdf;
 	}
 
 	void addChild(const std::string &name, ConfigurableObject *child) {

src/tests/test_chisquare.cpp

@@ -128,6 +128,10 @@ public:
 			EMeasure measure = ESolidAngle;
 			if (!sampled.isZero())
 				measure = BSDF::getMeasure(bRec.sampledType);
+
+			if (sampled.isZero() && sampled2.isZero())
+				return boost::make_tuple(Vector(0.0f), 0.0f, measure);
+
 			Spectrum f = m_bsdf->eval(bRec, measure);
 			pdfVal = m_bsdf->pdf(bRec, measure);
 			Spectrum manual = f/pdfVal;
@@ -178,8 +182,7 @@ public:
 				disableFPExceptions();
 			#endif
 
-			return boost::make_tuple(bRec.wo, 
+			return boost::make_tuple(bRec.wo, 1.0f, measure);
-				sampled.isZero() ? 0.0f : 1.0f, measure);
 		}
  
 		Float pdf(const Vector &wo, EMeasure measure) {
@@ -393,7 +396,7 @@ public:
 
 			Log(EInfo, "Checking the model for %i incident directions and 2D sampling", wiSamples);
 			progress->reset();
-
+#if 1
 			/* Test for a number of different incident directions */
 			for (size_t j=0; j<wiSamples; ++j) {
 				Vector wi;
@@ -429,6 +432,7 @@ public:
 				progress->update(j+1);
 			}
 			Log(EInfo, "The largest encountered importance weight was = %.2f", largestWeight);
+#endif
 			largestWeight = 0;
 
 			if (bsdf->getComponentCount() > 1) {