From ef48c84915dca73e9f96299a63ccb0dac8f8ea62 Mon Sep 17 00:00:00 2001 From: johannes hanika Date: Mon, 3 Nov 2014 17:46:29 +0100 Subject: [PATCH] hslt: improve numerical robustness --- src/bsdfs/microfacet.h | 5 +++-- src/bsdfs/roughconductor.cpp | 8 +++---- src/integrators/mlt/mlt_proc.cpp | 3 +++ src/libbidir/mut_bidir.cpp | 2 +- src/libbidir/mut_caustic.cpp | 6 ++++- src/libbidir/mut_lens.cpp | 6 ++++- src/libbidir/mut_manifold.cpp | 15 ++++++++++--- src/libbidir/vertex.cpp | 38 +++++++++++++++----------------- 8 files changed, 51 insertions(+), 32 deletions(-) diff --git a/src/bsdfs/microfacet.h b/src/bsdfs/microfacet.h index 9bc6af36..6ace9bf5 100644 --- a/src/bsdfs/microfacet.h +++ b/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; diff --git a/src/bsdfs/roughconductor.cpp b/src/bsdfs/roughconductor.cpp index 2146c60b..296506e3 100644 --- a/src/bsdfs/roughconductor.cpp +++ b/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.wo) < 0 || + Frame::cosTheta(bRec.wi) <= 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.wo) < 0 || + Frame::cosTheta(bRec.wi) <= 0 || + Frame::cosTheta(bRec.wo) <= 0 || ((bRec.component != -1 && bRec.component != 0) || !(bRec.typeMask & EGlossyReflection))) return 0.0f; diff --git a/src/integrators/mlt/mlt_proc.cpp b/src/integrators/mlt/mlt_proc.cpp index 6b42f3a0..402120c4 100644 --- a/src/integrators/mlt/mlt_proc.cpp +++ b/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(); diff --git a/src/libbidir/mut_bidir.cpp b/src/libbidir/mut_bidir.cpp index a531f2bc..7cf67be0 100644 --- a/src/libbidir/mut_bidir.cpp +++ b/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; } diff --git a/src/libbidir/mut_caustic.cpp b/src/libbidir/mut_caustic.cpp index 749f7b8f..bf351c23 100644 --- a/src/libbidir/mut_caustic.cpp +++ b/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 &) { diff --git a/src/libbidir/mut_lens.cpp b/src/libbidir/mut_lens.cpp index 047d0ef2..e2f938e5 100644 --- a/src/libbidir/mut_lens.cpp +++ b/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 &) { diff --git a/src/libbidir/mut_manifold.cpp b/src/libbidir/mut_manifold.cpp index 8c42c908..0ce9de7b 100644 --- a/src/libbidir/mut_manifold.cpp +++ b/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; } diff --git a/src/libbidir/vertex.cpp b/src/libbidir/vertex.cpp index 03dc2342..1cc440cd 100644 --- a/src/libbidir/vertex.cpp +++ b/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], - weightBkw = pdf[1-mode] == 0 ? 0.0f : 1.0f / pdf[1-mode]; + Float weightFwd = pdf[mode] <= RCPOVERFLOW ? 0 : 1 / pdf[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->weight[EImportance].isZero()) + if (!vs->update(scene, pred, vt, EImportance)) return false; - vt->update(scene, succ, vs, ERadiance); - if (vt->weight[ERadiance].isZero()) + if (!vt->update(scene, succ, vs, ERadiance)) 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->weight[EImportance].isZero() || vs->pdf[EImportance] == 0.f) + if (!vs->update(scene, pred, vt, EImportance, vsMeasure)) return false; - vt->update(scene, succ, vs, ERadiance, vtMeasure); - if (vt->weight[ERadiance].isZero() || vt->pdf[ERadiance] == 0.f) + if (!vt->update(scene, succ, vs, ERadiance, vtMeasure)) return false; return edge->connect(scene, predEdge, vs, vt, succEdge);