bugfix in volpath regarding index-matched participating media. improved the handling of the strictNormals parameter in all path tracers

include/mitsuba/render/scene.h

@@ -167,12 +167,17 @@ public:
 	 *    extent information, as well as a time (which applies when
 	 *    the shapes are animated)
 	 *
+	 * \param medium
+	 *    Initial medium containing the ray \c ray
+	 *
 	 * \param its
 	 *    A detailed intersection record, which will be filled by the
 	 *    intersection query
 	 *
-	 * \param medium
+	 * \param indexMatchedMediumTransition
-	 *    Initial medium containing the ray \c ray.
+	 *    This parameter is used to return whether or not this routine
+	 *    passed through an index-matched medium-to-medium transition 
+	 *    while computing the attenuation.
 	 *
 	 * \param transmittance
 	 *    Transmittance from the ray origin to the returned intersection
@@ -180,8 +185,9 @@ public:
 	 *
 	 * \return \c true if an intersection was found
 	 */
-	bool attenuatedRayIntersect(const Ray &ray, const Medium *medium,
+	bool attenuatedRayIntersect(const Ray &ray, const Medium * medium,
-		Intersection &its, Spectrum &transmittance, Sampler *sampler = NULL) const;
+		Intersection &its, bool &indexMatchedMediumTransition,
+		Spectrum &transmittance, Sampler *sampler = NULL) const;
 
 	/**
 	 * \brief Intersect a ray against all primitives stored in the scene

src/integrators/path/path.cpp

@@ -84,7 +84,7 @@ public:
 			/*                          Luminaire sampling                          */
 			/* ==================================================================== */
 
-			/* Prevent light leaks due to the use of shading normals -- [Veach, p. 158] */
+			/* Prevent light leaks due to the use of shading normals */
 			Float wiDotGeoN = -dot(its.geoFrame.n, ray.d),
 				  wiDotShN  = Frame::cosTheta(its.wi);
 			if (wiDotGeoN * wiDotShN < 0 && m_strictNormals) 
@@ -95,12 +95,17 @@ public:
 			if (rRec.type & RadianceQueryRecord::EDirectSurfaceRadiance && 
 				scene->sampleLuminaire(its.p, ray.time, lRec, rRec.nextSample2D())) {
 				/* Allocate a record for querying the BSDF */
-				const BSDFQueryRecord bRec(its, its.toLocal(-lRec.d));
+				const Vector wo = -lRec.d;
-
+				const BSDFQueryRecord bRec(its, its.toLocal(wo));
+	
 				/* Evaluate BSDF * cos(theta) */
 				const Spectrum bsdfVal = bsdf->fCos(bRec);
 
-				if (!bsdfVal.isZero()) {
+				Float woDotGeoN = dot(its.geoFrame.n, wo);
+
+				/* Prevent light leaks due to the use of shading normals */
+				if (!bsdfVal.isZero() && (!m_strictNormals
+						|| woDotGeoN * Frame::cosTheta(bRec.wo) > 0)) {
 					/* Calculate prob. of having sampled that direction
 					   using BSDF sampling */
 					Float bsdfPdf = (lRec.luminaire->isIntersectable() 
@@ -125,7 +130,7 @@ public:
 				break;
 			bsdfVal /= bsdfPdf;
 	
-			/* Prevent light leaks due to the use of shading normals -- [Veach, p. 158] */
+			/* Prevent light leaks due to the use of shading normals */
 			const Vector wo = its.toWorld(bRec.wo);
 			Float woDotGeoN = dot(its.geoFrame.n, wo);
 			if (woDotGeoN * Frame::cosTheta(bRec.wo) <= 0 && m_strictNormals)
@@ -178,8 +183,9 @@ public:
 				break;
 			rRec.type = RadianceQueryRecord::ERadianceNoEmission;
 
-			/* Russian roulette - Possibly stop the recursion. Don't use for transmission
+			/* Russian roulette - Possibly stop the recursion. Don't do this when
-			   due to IOR weighting factors, which throw the heuristic off */
+			   dealing with a transmission component, since solid angle compression
+			   factors cause problems with the heuristic below */
 			if (rRec.depth >= m_rrDepth && !(bRec.sampledType & BSDF::ETransmission)) {
 				/* Assuming that BSDF importance sampling is perfect,
 				   'bsdfVal.max()' should equal the maximum albedo

src/integrators/path/volpath.cpp

@@ -75,10 +75,11 @@ public:
 				/* ==================================================================== */
 				/*                          Luminaire sampling                          */
 				/* ==================================================================== */
-
+				
 				/* Estimate the single scattering component if this is requested */
 				if (rRec.type & RadianceQueryRecord::EDirectMediumRadiance && 
-					scene->sampleAttenuatedLuminaire(mRec.p, ray.time, rRec.medium, lRec, rRec.nextSample2D(), rRec.sampler)) {
+					scene->sampleAttenuatedLuminaire(mRec.p, ray.time, rRec.medium, 
+						lRec, rRec.nextSample2D(), rRec.sampler)) {
 					/* Evaluate the phase function */
 					Float phaseVal = phase->f(PhaseFunctionQueryRecord(mRec, -ray.d, -lRec.d));
 
@@ -94,7 +95,7 @@ public:
 						Li += pathThroughput * lRec.value * phaseVal * weight;
 					}
 				}
-
+				
 				/* ==================================================================== */
 				/*                         Phase function sampling                      */
 				/* ==================================================================== */
@@ -110,10 +111,11 @@ public:
 				ray = Ray(mRec.p, pRec.wo, ray.time);
 				ray.mint = 0;
 
-				bool hitLuminaire = false;
+				bool hitLuminaire = false, indexMatchedMediumTransition = false;
 				Spectrum transmittance;
 
-				if (scene->attenuatedRayIntersect(ray, rRec.medium, its, transmittance, rRec.sampler)) {
+				if (scene->attenuatedRayIntersect(ray, rRec.medium, its, 
+						indexMatchedMediumTransition, transmittance, rRec.sampler)) {
 					/* Intersected something - check if it was a luminaire */
 					if (its.isLuminaire()) {
 						lRec = LuminaireSamplingRecord(its, -ray.d);
@@ -130,7 +132,7 @@ public:
 						hitLuminaire = true;
 					}
 				}
-
+				
 				/* If a luminaire was hit, estimate the local illumination and
 				   weight using the power heuristic */
 				if (hitLuminaire && (rRec.type & RadianceQueryRecord::EDirectMediumRadiance)) {
@@ -140,6 +142,14 @@ public:
 					Li += pathThroughput * lRec.value * phaseVal * weight * transmittance;
 				}
 
+				if (indexMatchedMediumTransition) {
+					/* The previous ray intersection code passed through an index-matched
+					   medium transition while looking for a luminaire. For the recursion,
+					   we need to rewind and account for this transition -- therefore,
+					   another ray intersection call is neccessary */
+					scene->rayIntersect(ray, its);
+				}
+
 				/* ==================================================================== */
 				/*                         Multiple scattering                          */
 				/* ==================================================================== */
@@ -174,7 +184,7 @@ public:
 						Li += pathThroughput * scene->LeBackground(ray);
 					break;
 				}
-
+				
 				/* Possibly include emitted radiance if requested */
 				if (its.isLuminaire() && (rRec.type & RadianceQueryRecord::EEmittedRadiance))
 					Li += pathThroughput * its.Le(-ray.d);
@@ -185,19 +195,19 @@ public:
 
 				if (rRec.depth == m_maxDepth && m_maxDepth > 0)
 					break;
-
+				
 				const BSDF *bsdf = its.getBSDF(ray);
 				if (!bsdf) {
 					/* Pass right through the surface (there is no BSDF) */
-					if (its.isMediumTransition())
+					if (its.isMediumTransition()) 
 						rRec.medium = its.getTargetMedium(ray.d);
 					ray.setOrigin(its.p);
 					ray.mint = Epsilon;
 					scene->rayIntersect(ray, its);
 					continue;
 				}
-			
+
-				/* Prevent light leaks due to the use of shading normals -- [Veach, p. 158] */
+				/* Prevent light leaks due to the use of shading normals */
 				Float wiDotGeoN = -dot(its.geoFrame.n, ray.d),
 					  wiDotShN  = Frame::cosTheta(its.wi);
 				if (wiDotGeoN * wiDotShN < 0 && m_strictNormals) 
@@ -209,14 +219,21 @@ public:
 
 				/* Estimate the direct illumination if this is requested */
 				if (rRec.type & RadianceQueryRecord::EDirectSurfaceRadiance && 
-					scene->sampleAttenuatedLuminaire(its, rRec.medium, lRec, rRec.nextSample2D(), rRec.sampler)) {
+					scene->sampleAttenuatedLuminaire(its, rRec.medium, lRec, 
-					/* Allocate a record for querying the BSDF */
+						rRec.nextSample2D(), rRec.sampler)) {
-					const BSDFQueryRecord bRec(its, its.toLocal(-lRec.d));
+					const Vector wo = -lRec.d;
 
+					/* Allocate a record for querying the BSDF */
+					const BSDFQueryRecord bRec(its, its.toLocal(wo));
+	
 					/* Evaluate BSDF * cos(theta) */
 					const Spectrum bsdfVal = bsdf->fCos(bRec);
 
-					if (!bsdfVal.isZero()) {
+					Float woDotGeoN = dot(its.geoFrame.n, wo);
+
+					/* Prevent light leaks due to the use of shading normals */
+					if (!bsdfVal.isZero() && (!m_strictNormals
+							|| woDotGeoN * Frame::cosTheta(bRec.wo) > 0)) {
 						/* Calculate prob. of having sampled that direction
 						   using BSDF sampling */
 						Float bsdfPdf = (lRec.luminaire->isIntersectable() 
@@ -242,7 +259,7 @@ public:
 	
 				bsdfVal /= bsdfPdf;
 
-				/* Prevent light leaks due to the use of shading normals -- [Veach, p. 158] */
+				/* Prevent light leaks due to the use of shading normals */
 				const Vector wo = its.toWorld(bRec.wo);
 				Float woDotGeoN = dot(its.geoFrame.n, wo);
 				if (woDotGeoN * Frame::cosTheta(bRec.wo) <= 0 && m_strictNormals)
@@ -254,10 +271,11 @@ public:
 				if (its.isMediumTransition())
 					rRec.medium = its.getTargetMedium(ray.d);
 
-				bool hitLuminaire = false;
+				bool hitLuminaire = false, indexMatchedMediumTransition = false;
 				Spectrum transmittance;
 
-				if (scene->attenuatedRayIntersect(ray, rRec.medium, its, transmittance, rRec.sampler)) {
+				if (scene->attenuatedRayIntersect(ray, rRec.medium, its,
+						indexMatchedMediumTransition, transmittance, rRec.sampler)) {
 					/* Intersected something - check if it was a luminaire */
 					if (its.isLuminaire()) {
 						lRec = LuminaireSamplingRecord(its, -ray.d);
@@ -285,6 +303,14 @@ public:
 					Li += pathThroughput * lRec.value * bsdfVal * weight * transmittance;
 				}
 
+				if (indexMatchedMediumTransition) {
+					/* The previous ray intersection code passed through an index-matched
+					   medium transition while looking for a luminaire. For the recursion,
+					   we need to rewind and account for this transition -- therefore,
+					   another ray intersection call is neccessary */
+					scene->rayIntersect(ray, its);
+				}
+
 				/* ==================================================================== */
 				/*                         Indirect illumination                        */
 				/* ==================================================================== */
@@ -294,8 +320,9 @@ public:
 					break;
 				rRec.type = RadianceQueryRecord::ERadianceNoEmission;
 
-				/* Russian roulette - Possibly stop the recursion. Don't use for transmission
+				/* Russian roulette - Possibly stop the recursion. Don't do this when
-				   due to IOR weighting factors, which throw the heuristic off */
+				   dealing with a transmission component, since solid angle compression
+				   factors cause problems with the heuristic below */
 				if (rRec.depth >= m_rrDepth && !(bRec.sampledType & BSDF::ETransmission)) {
 					/* Assuming that BSDF importance sampling is perfect,
 					   'bsdfVal.max()' should equal the maximum albedo
@@ -317,8 +344,7 @@ public:
 	}
 
 	inline Float miWeight(Float pdfA, Float pdfB) const {
-		pdfA *= pdfA;
+		pdfA *= pdfA; pdfB *= pdfB;
-		pdfB *= pdfB;
 		return pdfA / (pdfA + pdfB);
 	}
 

src/integrators/path/volpath_simple.cpp

@@ -78,7 +78,6 @@ public:
 				if (rRec.type & RadianceQueryRecord::EDirectMediumRadiance && 
 					scene->sampleAttenuatedLuminaire(mRec.p, ray.time,
 						rRec.medium, lRec, rRec.nextSample2D(), rRec.sampler)) {
-
 					Li += pathThroughput * lRec.value * phase->f(
 							PhaseFunctionQueryRecord(mRec, -ray.d, -lRec.d));
 				}
@@ -153,7 +152,7 @@ public:
 					continue;
 				}
 				
-				/* Prevent light leaks due to the use of shading normals -- [Veach, p. 158] */
+				/* Prevent light leaks due to the use of shading normals */
 				Float wiDotGeoN = -dot(its.geoFrame.n, ray.d),
 					  wiDotShN  = Frame::cosTheta(its.wi);
 				if (wiDotGeoN * wiDotShN < 0 && m_strictNormals) 
@@ -168,8 +167,14 @@ public:
 					scene->sampleAttenuatedLuminaire(its, rRec.medium, lRec, 
 						rRec.nextSample2D(), rRec.sampler)) {
 					/* Allocate a record for querying the BSDF */
-					const BSDFQueryRecord bRec(its, its.toLocal(-lRec.d));
+					const Vector wo = -lRec.d;
-					Li += pathThroughput * lRec.value * bsdf->fCos(bRec);
+					const BSDFQueryRecord bRec(its, its.toLocal(wo));
+					
+					Float woDotGeoN = dot(its.geoFrame.n, wo);
+					/* Prevent light leaks due to the use of shading normals */
+					if (!m_strictNormals ||
+						woDotGeoN * Frame::cosTheta(bRec.wo) > 0)
+						Li += pathThroughput * lRec.value * bsdf->fCos(bRec);
 				}
 
 				/* ==================================================================== */
@@ -182,7 +187,7 @@ public:
 				if (bsdfVal.isZero()) 
 					break;
 	
-				/* Prevent light leaks due to the use of shading normals -- [Veach, p. 158] */
+				/* Prevent light leaks due to the use of shading normals */
 				const Vector wo = its.toWorld(bRec.wo);
 				Float woDotGeoN = dot(its.geoFrame.n, wo);
 				if (woDotGeoN * Frame::cosTheta(bRec.wo) <= 0 && m_strictNormals)
@@ -231,8 +236,9 @@ public:
 					}
 				}
 
-				/* Russian roulette - Possibly stop the recursion. Don't use for transmission
+				/* Russian roulette - Possibly stop the recursion. Don't do this when
-				   due to IOR weighting factors, which throw the heuristic off */
+				   dealing with a transmission component, since solid angle compression
+				   factors cause problems with the heuristic below */
 				if (rRec.depth >= m_rrDepth && !(bRec.sampledType & BSDF::ETransmission)) {
 					/* Assuming that BSDF importance sampling is perfect,
 					   'bsdfVal.max()' should equal the maximum albedo

src/librender/scene.cpp

@@ -431,7 +431,8 @@ Spectrum Scene::getTransmittance(const Point &p1, const Point &p2,
 }
 
 bool Scene::attenuatedRayIntersect(const Ray &_ray, const Medium *medium,
-		Intersection &its, Spectrum &transmittance, Sampler *sampler) const {
+		Intersection &its, bool &indexMatchedMediumTransition,
+		Spectrum &transmittance, Sampler *sampler) const {
 	Ray ray(_ray);
 	transmittance = Spectrum(1.0f);
 	int iterations = 0;
@@ -446,10 +447,12 @@ bool Scene::attenuatedRayIntersect(const Ray &_ray, const Medium *medium,
 			return false;
 		else if (its.shape->isOccluder())
 			return true;
-		else if (its.shape->isMediumTransition()) 
+		else if (its.shape->isMediumTransition()) {
 			medium = dot(its.geoFrame.n, ray.d) > 0 ?
 				  its.shape->getExteriorMedium()
 				: its.shape->getInteriorMedium();
+			indexMatchedMediumTransition = true;
+		}
 
 		ray.o = ray(its.t);
 		ray.mint = Epsilon;

src/librender/trimesh.cpp

@@ -23,6 +23,7 @@
 #include <mitsuba/core/timer.h>
 #include <mitsuba/core/properties.h>
 #include <mitsuba/render/subsurface.h>
+#include <mitsuba/render/medium.h>
 #include <mitsuba/render/bsdf.h>
 #include <mitsuba/render/luminaire.h>
 
@@ -803,8 +804,12 @@ std::string TriMesh::toString() const {
 		<< "  surfaceArea = " << m_surfaceArea << "," << endl
 		<< "  aabb = " << m_aabb.toString() << "," << endl
 		<< "  bsdf = " << indent(m_bsdf.toString()) << "," << endl
-		<< "  subsurface = " << indent(m_subsurface.toString()) << "," << endl
+		<< "  subsurface = " << indent(m_subsurface.toString()) << "," << endl;
-		<< "  luminaire = " << indent(m_luminaire.toString()) << endl
+	if (isMediumTransition()) {
+		oss << "  interiorMedium = " << indent(m_interiorMedium.toString()) << "," << endl
+			<< "  exteriorMedium = " << indent(m_exteriorMedium.toString()) << "," << endl;
+	}
+	oss << "  luminaire = " << indent(m_luminaire.toString()) << endl
 		<< "]";
 	return oss.str();
 }