hslt: improve numerical robustness

src/bsdfs/microfacet.h

@@ -550,11 +550,12 @@ protected:
 
 	/// Compute the interpolated roughness for the Phong model
 	inline Float interpolatePhongExponent(const Vector &v) const {
-		Float invSinTheta2 = 1 / Frame::sinTheta2(v);
+		const Float sinTheta2 = Frame::sinTheta2(v);
 
-		if (isIsotropic() || invSinTheta2 <= 0)
+		if (isIsotropic() || sinTheta2 <= RCPOVERFLOW)
 			return m_exponentU;
 
+		Float invSinTheta2 = 1 / sinTheta2;
 		Float cosPhi2 = v.x * v.x * invSinTheta2;
 		Float sinPhi2 = v.y * v.y * invSinTheta2;
 

src/bsdfs/roughconductor.cpp

@@ -257,8 +257,8 @@ public:
 	Spectrum eval(const BSDFSamplingRecord &bRec, EMeasure measure) const {
 		/* Stop if this component was not requested */
 		if (measure != ESolidAngle ||
-			Frame::cosTheta(bRec.wi) < 0 ||
+			Frame::cosTheta(bRec.wi) <= 0 ||
-			Frame::cosTheta(bRec.wo) < 0 ||
+			Frame::cosTheta(bRec.wo) <= 0 ||
 			((bRec.component != -1 && bRec.component != 0) ||
 			!(bRec.typeMask & EGlossyReflection)))
 			return Spectrum(0.0f);
@@ -295,8 +295,8 @@ public:
 
 	Float pdf(const BSDFSamplingRecord &bRec, EMeasure measure) const {
 		if (measure != ESolidAngle ||
-			Frame::cosTheta(bRec.wi) < 0 ||
+			Frame::cosTheta(bRec.wi) <= 0 ||
-			Frame::cosTheta(bRec.wo) < 0 ||
+			Frame::cosTheta(bRec.wo) <= 0 ||
 			((bRec.component != -1 && bRec.component != 0) ||
 			!(bRec.typeMask & EGlossyReflection)))
 			return 0.0f;

src/integrators/mlt/mlt_proc.cpp

@@ -206,7 +206,10 @@ public:
 
 				Float a;
 				if (!m_config.importanceMap) {
+					if(Qxy > RCPOVERFLOW)
 					a = std::min((Float) 1, Qyx / Qxy);
+					else
+						a = 0.f;
 				} else {
 					const Float *luminanceValues = m_config.importanceMap->getFloatData();
 					const Point2 &curPos = current->getSamplePosition();

src/libbidir/mut_bidir.cpp

@@ -260,7 +260,7 @@ Float BidirectionalMutator::Q(const Path &source, const Path &proposal,
 
 		Float luminance = weight.getLuminance();
 
-		if (luminance <= 0 || !std::isfinite(luminance)) {
+		if (luminance <= RCPOVERFLOW || !std::isfinite(luminance)) {
 			Log(EWarn, "Internal error: luminance = %f!", luminance);
 			continue;
 		}

src/libbidir/mut_caustic.cpp

@@ -212,7 +212,11 @@ Float CausticPerturbation::Q(const Path &source, const Path &proposal,
 					source.edge(i), edge);
 	}
 
-	return solidAngleDensity / weight.getLuminance();
+	const Float lumWeight = weight.getLuminance();
+	if(lumWeight <= RCPOVERFLOW)
+		return 0.f;
+
+	return solidAngleDensity / lumWeight;
 }
 
 void CausticPerturbation::accept(const MutationRecord &) {

src/libbidir/mut_lens.cpp

@@ -196,7 +196,11 @@ Float LensPerturbation::Q(const Path &source, const Path &proposal,
 					source.edge(i-1), edge);
 	}
 
-	return 1.0f / weight.getLuminance();
+	const Float lumWeight = weight.getLuminance();
+	if(lumWeight <= RCPOVERFLOW)
+		return 0.f;
+
+	return 1.0f / lumWeight;
 }
 
 void LensPerturbation::accept(const MutationRecord &) {

src/libbidir/mut_manifold.cpp

@@ -663,10 +663,19 @@ Float ManifoldPerturbation::Q(const Path &source, const Path &proposal,
 		/* Convert to area density at x_b */
 		prob *= m_manifold->G(proposal, a, b);
 
-		if (prob == 0)
+		if (prob <= RCPOVERFLOW)
 			return 0.0f;
+
+#if defined(MTS_DEBUG_FP)
+		disableFPExceptions();
+#endif
+
 		weight /= prob;
 
+#if defined(MTS_DEBUG_FP)
+		enableFPExceptions();
+#endif
+
 		/* Catch very low probabilities which round to +inf in the above division operation */
 		if (!std::isfinite(weight.average()))
 			return 0.0f;
@@ -677,7 +686,7 @@ Float ManifoldPerturbation::Q(const Path &source, const Path &proposal,
 			source.vertex(a)->pdf[mode] * m_probFactor * m_probFactor);
 
 		Float pdf = source.vertex(a+step)->perturbPositionPdf(proposal.vertex(a+step), stddev);
-		if (pdf == 0)
+		if (pdf <= RCPOVERFLOW)
 			return 0.0f;
 
 		weight /= pdf;
@@ -743,7 +752,7 @@ Float ManifoldPerturbation::Q(const Path &source, const Path &proposal,
 
 	Float lum = weight.getLuminance();
 
-	if (lum <= 0 || !std::isfinite(lum)) {
+	if (lum <= RCPOVERFLOW || !std::isfinite(lum)) {
 		Log(EWarn, "Internal error in manifold perturbation: luminance = %f!", lum);
 		return 0.f;
 	}

src/libbidir/vertex.cpp

@@ -192,7 +192,7 @@ bool PathVertex::sampleNext(const Scene *scene, Sampler *sampler,
 				/* Compute the reverse quantities */
 				bRec.reverse();
 				pdf[1-mode] = bsdf->pdf(bRec, (EMeasure) measure);
-				if (pdf[1-mode] == 0) {
+				if (pdf[1-mode] <= RCPOVERFLOW) {
 					/* This can happen rarely due to roundoff errors -- be strict */
 					return false;
 				}
@@ -508,7 +508,7 @@ bool PathVertex::perturbDirection(const Scene *scene, const PathVertex *pred,
 				Spectrum value = emitter->evalDirection(dRec, pRec);
 				Float prob = emitter->pdfDirection(dRec, pRec);
 
-				if (value.isZero() || prob == 0)
+				if (value.isZero() || prob <= RCPOVERFLOW)
 					return false;
 
 				weight[EImportance] = value/prob;
@@ -531,7 +531,7 @@ bool PathVertex::perturbDirection(const Scene *scene, const PathVertex *pred,
 				Spectrum value = sensor->evalDirection(dRec, pRec);
 				Float prob = sensor->pdfDirection(dRec, pRec);
 
-				if (value.isZero() || prob == 0)
+				if (value.isZero() || prob <= RCPOVERFLOW)
 					return false;
 
 				weight[EImportance] = value * (
@@ -556,7 +556,7 @@ bool PathVertex::perturbDirection(const Scene *scene, const PathVertex *pred,
 				Spectrum value = bsdf->eval(bRec);
 				Float prob = bsdf->pdf(bRec);
 
-				if (value.isZero() || prob == 0)
+				if (value.isZero() || prob <= RCPOVERFLOW)
 					return false;
 
 				weight[mode] = value/prob;
@@ -592,7 +592,7 @@ bool PathVertex::perturbDirection(const Scene *scene, const PathVertex *pred,
 				/* Compute the reverse quantities */
 				bRec.reverse();
 				pdf[1-mode] = bsdf->pdf(bRec, ESolidAngle);
-				if (pdf[1-mode] == 0) {
+				if (pdf[1-mode] <= RCPOVERFLOW) {
 					/* This can happen rarely due to roundoff errors -- be strict */
 					return false;
 				}
@@ -628,7 +628,7 @@ bool PathVertex::perturbDirection(const Scene *scene, const PathVertex *pred,
 				Float value = phase->eval(pRec);
 				Float prob = phase->pdf(pRec);
 
-				if (value == 0 || prob == 0)
+				if (value == 0 || prob <= RCPOVERFLOW)
 					return false;
 
 				pdf[mode] = prob;
@@ -710,7 +710,7 @@ bool PathVertex::propagatePerturbation(const Scene *scene, const PathVertex *pre
 
 	bRec.typeMask = BSDF::EAll;
 	Float prob = bsdf->pdf(bRec, EDiscrete);
-	if (prob == 0) {
+	if (prob <= RCPOVERFLOW) {
 		SLog(EWarn, "Unable to recreate specular vertex in perturbation (bsdf=%s)",
 			bsdf->toString().c_str());
 		return false;
@@ -721,7 +721,7 @@ bool PathVertex::propagatePerturbation(const Scene *scene, const PathVertex *pre
 	measure = EDiscrete;
 	componentType = componentType_;
 
-	if (weight[mode].isZero() || prob == 0)
+	if (weight[mode].isZero() || prob <= RCPOVERFLOW)
 		return false;
 
 	/* Account for medium changes if applicable */
@@ -747,7 +747,7 @@ bool PathVertex::propagatePerturbation(const Scene *scene, const PathVertex *pre
 	/* Compute the reverse quantities */
 	bRec.reverse();
 	pdf[1-mode] = bsdf->pdf(bRec, EDiscrete);
-	if (pdf[1-mode] == 0) {
+	if (pdf[1-mode] <= RCPOVERFLOW) {
 		/* This can happen rarely due to roundoff errors -- be strict */
 		return false;
 	}
@@ -1128,6 +1128,7 @@ bool PathVertex::cast(const Scene *scene, EVertexType desired) {
 		pRec.object = emitter;
 		pRec.pdf = 0.0f;
 		getPositionSamplingRecord() = pRec;
+		measure = pRec.measure;
 		degenerate = emitter->getType() & Emitter::EDeltaDirection;
 
 		return true;
@@ -1149,6 +1150,7 @@ bool PathVertex::cast(const Scene *scene, EVertexType desired) {
 		Vector2i size = sensor->getFilm()->getSize();
 		pRec.uv.x *= size.x; pRec.uv.y *= size.y;
 		getPositionSamplingRecord() = pRec;
+		measure = pRec.measure;
 		degenerate = sensor->getType() & Sensor::EDeltaDirection;
 
 		return true;
@@ -1167,11 +1169,11 @@ bool PathVertex::update(const Scene *scene, const PathVertex *pred,
 	weight[mode]    = eval(scene, pred, succ, mode, measure);
 	weight[1-mode]  = eval(scene, succ, pred, (ETransportMode) (1-mode), measure);
 
-	if (weight[mode].isZero() || pdf[mode] == 0)
+	if (weight[mode].isZero() || pdf[mode] <= RCPOVERFLOW)
 		return false;
 
-	Float weightFwd = pdf[mode]   == 0 ? 0.0f : 1.0f / pdf[mode],
+	Float weightFwd = pdf[mode]   <= RCPOVERFLOW ? 0 : 1 / pdf[mode],
-	      weightBkw = pdf[1-mode] == 0 ? 0.0f : 1.0f / pdf[1-mode];
+		  weightBkw = pdf[1-mode] <= RCPOVERFLOW ? 0 : 1 / pdf[1-mode];
 
 	this->measure = measure;
 
@@ -1331,12 +1333,10 @@ bool PathVertex::connect(const Scene *scene,
 	if (vs->isDegenerate() || vt->isDegenerate())
 		return false;
 
-	vs->update(scene, pred, vt, EImportance);
+	if (!vs->update(scene, pred, vt, EImportance))
-	if (vs->weight[EImportance].isZero())
 		return false;
 
-	vt->update(scene, succ, vs, ERadiance);
+	if (!vt->update(scene, succ, vs, ERadiance))
-	if (vt->weight[ERadiance].isZero())
 		return false;
 
 	return edge->connect(scene, predEdge, vs, vt, succEdge);
@@ -1357,12 +1357,10 @@ bool PathVertex::connect(const Scene *scene,
 			return false;
 	}
 
-	vs->update(scene, pred, vt, EImportance, vsMeasure);
+	if (!vs->update(scene, pred, vt, EImportance, vsMeasure))
-	if (vs->weight[EImportance].isZero() || vs->pdf[EImportance] == 0.f)
 		return false;
 
-	vt->update(scene, succ, vs, ERadiance, vtMeasure);
+	if (!vt->update(scene, succ, vs, ERadiance, vtMeasure))
-	if (vt->weight[ERadiance].isZero() || vt->pdf[ERadiance] == 0.f)
 		return false;
 
 	return edge->connect(scene, predEdge, vs, vt, succEdge);