various cosmetic changes involving dielectric materials, still debugging roughglass..

include/mitsuba/core/util.h

@@ -423,30 +423,30 @@ extern MTS_EXPORT_CORE Point2 squareToTriangle(const Point2 &sample);
  * Calculates the unpolarized fresnel reflection coefficient for a 
  * dielectric material
  *
- * \param cosTheta1
+ * \param cosThetaI
  * 		Cosine of the angle between the normal and the incident ray
- * \param cosTheta2
+ * \param cosThetaT
  * 		Cosine of the angle between the normal and the transmitted ray
- * \param etaExt
+ * \param etaI
- * 		Refraction coefficient outside of the material
+ * 		Refraction coefficient at the incident direction
- * \param etaInt
+ * \param etaT
- * 		Refraction coefficient inside the material
+ * 		Refraction coefficient at the transmitted direction
  */
-extern MTS_EXPORT_CORE Float fresnelDielectric(Float cosTheta1, 
+extern MTS_EXPORT_CORE Float fresnelDielectric(Float cosThetaI, 
-		Float cosTheta2, Float etaExt, Float etaInt);
+		Float cosThetaT, Float etaI, Float etaT);
 
 /**
  * Calculates the unpolarized fresnel reflection coefficient for a 
  * dielectric material. Handles incidence from either sides.
  *
- * \param cosTheta1
+ * \param cosThetaI
  * 		Cosine of the angle between the normal and the incident ray
  * \param etaExt
  * 		Refraction coefficient outside of the material
  * \param etaInt
  * 		Refraction coefficient inside the material
  */
-extern MTS_EXPORT_CORE Float fresnel(Float cosTheta1, Float etaExt,
+extern MTS_EXPORT_CORE Float fresnel(Float cosThetaI, Float etaExt,
 		Float etaInt);
 
 /**

src/bsdfs/dielectric.cpp

@@ -97,37 +97,36 @@ public:
 
 	Float refract(Float intIOR, Float extIOR, 
 			const Vector &wi, Vector &wo, ETransportQuantity quantity) const {
-		Float cosTheta1 = Frame::cosTheta(wi);
+		Float cosThetaI = Frame::cosTheta(wi),
-		bool entering = cosTheta1 > 0.0f;
+			  etaI = extIOR, etaT = intIOR;
+		bool entering = cosThetaI > 0.0f;
 
-		/* Swap the indices of refraction if the interaction starts
-		   at the inside of the object */
 		if (!entering)
-			std::swap(intIOR, extIOR);
+			std::swap(etaT, etaI);
 
-		Float eta = extIOR/intIOR;
+		Float eta = etaI / etaT;
 
 		/* Using Snell's law, calculate the squared sine of the
 		   angle between the normal and the transmitted ray */
-		Float sinTheta2Sqr = eta*eta * Frame::sinTheta2(wi);
+		Float sinThetaTSqr = eta*eta * Frame::sinTheta2(wi);
 
-		if (sinTheta2Sqr > 1.0f) /* Total internal reflection! */
+		if (sinThetaTSqr > 1.0f) /* Total internal reflection! */
 			return 0.0f;
 
 		/* Compute the cosine, but guard against numerical imprecision */
-		Float cosTheta2 = std::sqrt(std::max((Float) 0.0f, 1.0f - sinTheta2Sqr));
+		Float cosThetaT = std::sqrt(1.0f - sinThetaTSqr);
 		if (entering)
-			cosTheta2 = -cosTheta2;
+			cosThetaT = -cosThetaT;
 
-		/* With cos(N, transmittedRay) on tap, calculating the 
+		/* With cos(N, transmittedRay) avilable, calculating the 
 		   transmission direction is straightforward */
-		wo = Vector(-eta*wi.x, -eta*wi.y, cosTheta2);
+		wo = Vector(-eta*wi.x, -eta*wi.y, cosThetaT);
 
 		/* Finally compute transmission coefficient. When transporting
-		   radiance, account for the change at boundaries with different 
+		   radiance, account for the solid angle change at boundaries 
-		   indices of refraction. */
+		   with different indices of refraction. */
 		if (quantity == ERadiance)
-			return (extIOR*extIOR)/(intIOR*intIOR);
+			return (etaI*etaI) / (etaT*etaT);
 		else
 			return 1.0f;
 	}

src/bsdfs/roughglass.cpp

@@ -155,8 +155,8 @@ public:
 		}
 
 		/* Prevent numerical issues in other stages of the model */
-		if (result < 1e-10)
+		//if (result < 1e-40)
-			result = 0;
+		//	result = 0;
 
 		return result;
 	}
@@ -244,7 +244,9 @@ public:
 				Log(EError, "Invalid distribution function!");
 		}
 
-		return Normal(sphericalDirection(thetaM, phiM));
+		Normal result(sphericalDirection(thetaM, phiM));
+		Assert(result.z >= 0);
+		return result;
 	}
 
 
@@ -256,43 +258,33 @@ public:
 		return (value < 0) ? -1.0f : 1.0f;
 	}
 
-	Float refract(const Vector &wi, Vector &wo, ETransportQuantity quantity) const {
+	bool refract(const Vector &wi, Vector &wo) const {
-		Float cosTheta1 = Frame::cosTheta(wi);
+		Float cosThetaI = Frame::cosTheta(wi),
-		Float intIOR = m_intIOR, extIOR = m_extIOR;
+			  etaI = m_extIOR, etaT = m_intIOR;
-		bool entering = cosTheta1 > 0.0f;
+		bool entering = cosThetaI > 0.0f;
 
-		/* Swap the indices of refraction if the interaction starts
-		   at the inside of the object */
 		if (!entering)
-			std::swap(intIOR, extIOR);
+			std::swap(etaT, etaI);
 
-		Float eta = extIOR/intIOR;
+		Float eta = etaI / etaT;
 
 		/* Using Snell's law, calculate the squared sine of the
 		   angle between the normal and the transmitted ray */
-		Float sinTheta2Sqr = eta*eta * Frame::sinTheta2(wi);
+		Float sinThetaTSqr = eta*eta * Frame::sinTheta2(wi);
 
-		if (sinTheta2Sqr > 1.0f) /* Total internal reflection! */
+		if (sinThetaTSqr > 1.0f) /* Total internal reflection! */
-			return 0.0f;
+			return false;
 
-		/* Use the sin^2+cos^2=1 identity - max() guards against
+		/* Compute the cosine, but guard against numerical imprecision */
-		   numerical imprecision*/
+		Float cosThetaT = std::sqrt(1.0f - sinThetaTSqr);
-		Float cosTheta2 = std::sqrt(std::max((Float) 0.0f, 1.0f - sinTheta2Sqr));
 		if (entering)
-			cosTheta2 = -cosTheta2;
+			cosThetaT = -cosThetaT;
 
-		/* Having cos(N, transmittedRay), calculating the actual
+		/* With cos(N, transmittedRay) avilable, calculating the 
-		   direction becomes easy. */
+		   transmission direction is straightforward */
-		wo = Vector(-eta*wi.x, -eta*wi.y, cosTheta2);
+		wo = Vector(-eta*wi.x, -eta*wi.y, cosThetaT);
 
-		/* Finally compute transmission coefficient. When transporting
+		return true;
-		   radiance, account for the change at boundaries with different 
-		   indices of refraction. */
-		if (quantity == ERadiance)
-			return (extIOR*extIOR)/(intIOR*intIOR) * 
-				(1.0f - fresnel(Frame::cosTheta(wi), extIOR, intIOR));
-		else
-			return 1.0f - fresnel(Frame::cosTheta(wi), extIOR, intIOR);
 	}
 
 	inline Spectrum fReflection(const BSDFQueryRecord &bRec) const {
@@ -396,19 +388,19 @@ public:
 	}
 
 	inline Float pdfTransmission(const BSDFQueryRecord &bRec, Float alpha) const {
-		Float etaI = m_extIOR, etaO = m_intIOR;
+		Float etaI = m_extIOR, etaT = m_intIOR;
 
 		if (Frame::cosTheta(bRec.wi) * Frame::cosTheta(bRec.wo) >= 0)
 			return 0.0f;
 
 		if (Frame::cosTheta(bRec.wi) < 0)
-			std::swap(etaI, etaO);
+			std::swap(etaI, etaT);
 
-		Vector Ht = -normalize(bRec.wi*etaI+bRec.wo*etaO);
+		Vector Ht = -normalize(bRec.wi*etaI + bRec.wo*etaT);
 
 		/* Jacobian of the half-direction transform. */
-		Float sqrtDenom = etaI * dot(bRec.wi, Ht) + etaO * dot(bRec.wo, Ht);
+		Float sqrtDenom = etaI * dot(bRec.wi, Ht) + etaT * dot(bRec.wo, Ht);
-		Float dwht_dwo = (etaO*etaO * absDot(bRec.wo, Ht)) / (sqrtDenom*sqrtDenom);
+		Float dwht_dwo = (etaT*etaT * absDot(bRec.wo, Ht)) / (sqrtDenom*sqrtDenom);
 
 		return evalD(Ht, alpha) * std::abs(Frame::cosTheta(Ht)) * dwht_dwo;
 	}
@@ -422,11 +414,12 @@ public:
 		/* Suggestion by Bruce Walter: sample using a slightly different 
 		   value of alpha. This in practice limits the weights to 
 		   values <= 4. See also \ref sample() */
-		Float alpha = m_alpha * (1.2f - 0.2f * std::sqrt(
+//		Float alpha = m_alpha * (1.2f - 0.2f * std::sqrt(
-				std::abs(Frame::cosTheta(bRec.wi))));
+//				std::abs(Frame::cosTheta(bRec.wi))));
+		Float alpha = m_alpha;
 
 		if (hasReflection && hasTransmission) {
-			/* PDF for importance sampling according to approximate F
+			/* PDF for importance sampling according to approximate 
 			   Fresnel coefficients (approximate, because we don't know 
 			   the microfacet normal at this point) */
 			Float fr = fresnel(Frame::cosTheta(bRec.wi), m_extIOR, m_intIOR);
@@ -466,14 +459,31 @@ public:
 	}
 
 	inline Spectrum sampleTransmission(BSDFQueryRecord &bRec, Float alpha, const Point2 &sample) const {
+		Vector m = sampleD(sample, alpha);
+		cout << endl;
+#if 1
+		Float etaI = m_extIOR, etaT = m_intIOR;
+
+		if (Frame::cosTheta(bRec.wi) < 0)
+			std::swap(etaI, etaT);
+
+		Float eta = etaI / etaT, c = dot(bRec.wi, m);
+
+		Vector wo2 = m * (eta*c - signum(bRec.wi.z) * std::sqrt(1 + eta * (c*c-1)))
+			- bRec.wi * eta;
+#endif
+		cout << "wo2: " << wo2.toString() << endl;
+#if 1
 		/* Sample M, the microsurface normal */
-		Frame mFrame(sampleD(sample, alpha));
+		Frame mFrame(m);
 
 		/* Perfect specular reflection along the microsurface normal */
-		if (refract(mFrame.toLocal(bRec.wi), bRec.wo, bRec.quantity) == 0)
+		if (!refract(mFrame.toLocal(bRec.wi), bRec.wo))
 			return Spectrum(0.0f);
 
 		bRec.wo = mFrame.toWorld(bRec.wo);
+#endif
+		cout << "bRec.wo: " << bRec.wo.toString() << endl;
 
 		bRec.sampledComponent = 1;
 		bRec.sampledType = EGlossyTransmission;
@@ -502,11 +512,13 @@ public:
 		   value of alpha. This in practice limits the weights to 
 		   values <= 4. The change is of course also accounted for 
 		   in \ref pdf(), hence no error is introduced. */
-		Float alpha = m_alpha * (1.2f - 0.2f * std::sqrt(
+		//Float alpha = m_alpha * (1.2f - 0.2f * std::sqrt(
-				std::abs(Frame::cosTheta(bRec.wi))));
+		//		std::abs(Frame::cosTheta(bRec.wi))));
+		
+		Float alpha = m_alpha;
 
 		if (hasReflection && hasTransmission) {
-			/* PDF for importance sampling according to approximate F
+			/* PDF for importance sampling according to approximate 
 			   Fresnel coefficients (approximate, because we don't know 
 			   the microfacet normal at this point) */
 			Float fr = fresnel(Frame::cosTheta(bRec.wi), m_extIOR, m_intIOR);
@@ -515,7 +527,7 @@ public:
 				sample.x /= fr;
 				return sampleReflection(bRec, alpha, sample);
 			} else {
-				sample.x = (sample.x - fr) / (1-fr);
+				sample.x = (sample.x - fr) / (1 - fr);
 				return sampleTransmission(bRec, alpha, sample);
 			}
 		} else if (hasReflection) {

src/libcore/chisquare.cpp

@@ -91,7 +91,7 @@ void ChiSquare::fill(
 	Float min[2], max[2];
 	size_t idx = 0;
 
-	NDIntegrator integrator(1, 2, 100000, 0, 1e-6f);
+	NDIntegrator integrator(1, 2, 1000000, 1e-8f, 1e-8f);
 	Float maxError = 0, integral = 0;
 	for (int i=0; i<m_thetaBins; ++i) {
 		min[0] = i * factor.x;

src/libcore/util.cpp

@@ -678,11 +678,11 @@ Float lanczosSinc(Float t, Float tau) {
 	PBRT and the paper "Derivation of Refraction Formulas" 
 	by Paul S. Heckbert. */
 Float fresnelDielectric(Float cosTheta1, Float cosTheta2, 
-						Float etaExt, Float etaInt) {
+						Float etaI, Float etaT) {
-	Float Rs = (etaExt * cosTheta1 - etaInt * cosTheta2)
+	Float Rs = (etaI * cosTheta1 - etaT * cosTheta2)
-			/ (etaExt * cosTheta1 + etaInt * cosTheta2);
+			/ (etaI * cosTheta1 + etaT * cosTheta2);
-	Float Rp = (etaInt * cosTheta1 - etaExt * cosTheta2)
+	Float Rp = (etaT * cosTheta1 - etaI * cosTheta2)
-			/ (etaInt * cosTheta1 + etaExt * cosTheta2);
+			/ (etaT * cosTheta1 + etaI * cosTheta2);
 
 	return (Rs * Rs + Rp * Rp) / 2.0f;
 }
@@ -701,28 +701,27 @@ Spectrum fresnelConductor(Float cosTheta, const Spectrum &eta, const Spectrum &k
 	return (rParl2 + rPerp2) / 2.0f;
 }
 
-Float fresnel(Float cosTheta1, Float etaExt, Float etaInt) {
+Float fresnel(Float cosThetaI, Float etaExt, Float etaInt) {
+	Float etaI = etaExt, etaT = etaInt;
+
 	/* Swap the indices of refraction if the interaction starts
 		at the inside of the object */
-	if (cosTheta1 < 0.0f)
+	if (cosThetaI < 0.0f)
-		std::swap(etaInt, etaExt);
+		std::swap(etaI, etaT);
 
 	/* Using Snell's law, calculate the sine of the angle
 		between the transmitted ray and the surface normal */
-	Float sinTheta2 = etaExt/etaInt * 
+	Float sinThetaT = etaI / etaT * 
-		std::sqrt(std::max((Float) 0.0f, 1.0f - cosTheta1*cosTheta1));
+		std::sqrt(std::max((Float) 0.0f, 1.0f - cosThetaI*cosThetaI));
 
-	if (sinTheta2 > 1.0f)
+	if (sinThetaT > 1.0f)
 		return 1.0f;  /* Total internal reflection! */
 
-	/* Use the sin^2+cos^2=1 identity - max() guards against
+	Float cosThetaT = std::sqrt(1.0f - sinThetaT*sinThetaT);
-		numerical imprecision*/
-	Float cosTheta2 = std::sqrt(std::max((Float) 0.0f, 
-		1.0f - sinTheta2*sinTheta2));
 
 	/* Finally compute the reflection coefficient */
-	return fresnelDielectric(std::abs(cosTheta1), cosTheta2, 
+	return fresnelDielectric(std::abs(cosThetaI),
-		etaInt, etaExt);
+		cosThetaT, etaI, etaT);
 }
 
 Float radicalInverse(int b, size_t i) {

src/tests/test_chisquare.cpp

@@ -38,7 +38,8 @@ public:
 	class BSDFAdapter {
 	public:
 		BSDFAdapter(const BSDF *bsdf, Random *random, const Vector &wi, int component = -1)
-				: m_bsdf(bsdf), m_random(random), m_wi(wi), m_component(component) { }
+			: m_bsdf(bsdf), m_random(random), m_wi(wi), m_component(component),
+			  m_largestWeight(0) { }
 
 		std::pair<Vector, Float> generateSample() {
 			Point2 sample(m_random->nextFloat(), m_random->nextFloat());
@@ -69,8 +70,10 @@ public:
 
 			for (int i=0; i<SPECTRUM_SAMPLES; ++i) {
 				Float a = sampled[i], b=sampled2[i];
-				Float min = std::min(std::abs(a), std::abs(b));
+				SAssert(a >= 0 && b >= 0);
-				Float err = std::abs(a-b);
+				Float min = std::min(a, b);
+				Float err = std::abs(a - b);
+				m_largestWeight = std::max(m_largestWeight, a);
 
 				if (min < Epsilon && err > Epsilon) // absolute error threshold
 					mismatch = true;
@@ -95,11 +98,14 @@ public:
 				return 0.0f;
 			return m_bsdf->pdf(bRec);
 		}
+
+		inline Float getLargestWeight() const { return m_largestWeight; }
 	private:
 		ref<const BSDF> m_bsdf;
 		ref<Random> m_random;
 		Vector m_wi;
 		int m_component;
+		Float m_largestWeight;
 	};
 	
 	void test01_BSDF() {
@@ -116,8 +122,10 @@ public:
 				continue;
 
 			const BSDF *bsdf = static_cast<const BSDF *>(objects[i]);
+			Float largestWeight = 0;
 
 			Log(EInfo, "Processing BSDF model %s", bsdf->toString().c_str());
+#if 0
 			Log(EInfo, "Checking the combined model for %i incident directions", wiSamples);
 
 			/* Test for a number of different incident directions */
@@ -135,8 +143,8 @@ public:
 
 				// Initialize the tables used by the chi-square test
 				chiSqr->fill(
-					boost::bind(&BSDFAdapter::generateSample, adapter),
+					boost::bind(&BSDFAdapter::generateSample, &adapter),
-					boost::bind(&BSDFAdapter::pdf, adapter, _1)
+					boost::bind(&BSDFAdapter::pdf, &adapter, _1)
 				);
 
 				// (the following assumes that the distribution has 1 parameter, e.g. exponent value)
@@ -150,8 +158,10 @@ public:
 				} else {
 					succeed();
 				}
+				largestWeight = std::max(largestWeight, adapter.getLargestWeight());
 				++testCount;
 			}
+#endif
 
 			if (bsdf->getComponentCount() > 1) {
 				for (int comp=0; comp<bsdf->getComponentCount(); ++comp) {
@@ -167,13 +177,14 @@ public:
 							wi = squareToHemispherePSA(Point2(random->nextFloat(), random->nextFloat()));
 
 						BSDFAdapter adapter(bsdf, random, wi, comp);
+
 						ref<ChiSquare> chiSqr = new ChiSquare(thetaBins, 2*thetaBins, wiSamples);
 						chiSqr->setLogLevel(EDebug);
 
 						// Initialize the tables used by the chi-square test
 						chiSqr->fill(
-							boost::bind(&BSDFAdapter::generateSample, adapter),
+							boost::bind(&BSDFAdapter::generateSample, &adapter),
-							boost::bind(&BSDFAdapter::pdf, adapter, _1)
+							boost::bind(&BSDFAdapter::pdf, &adapter, _1)
 						);
 
 						// (the following assumes that the distribution has 1 parameter, e.g. exponent value)
@@ -187,10 +198,12 @@ public:
 						} else {
 							succeed();
 						}
+						largestWeight = std::max(largestWeight, adapter.getLargestWeight());
 						++testCount;
 					}
 				}
 			}
+			Log(EInfo, "Done with this model. The largest encountered sampling weight was = %f", largestWeight);
 		}
 		Log(EInfo, "%i/%i BSDF checks succeeded", testCount-failureCount, testCount);
 	}