VolPathIntegrator: track beta and 1 / path MIS weights.

(Rather than T_maj and path PDFs.)

This is "just" a matter of dividing the old T_maj by the PDF for sampling
the path and then dividing the old path PDFs by that same factor, so the
same mathematically in the end.

However, it a) better aligns the VolPathIntegrator with the other
integrators, which track beta, and b) fixes the issue of T_maj and the PDFs
coming to have very large or very small magnitudes and being succeptible to
over/underflow.  Now, all of the Rescale() stuff is no longer necessary for
them...

src/pbrt/cpu/integrators.cpp

@@ -954,7 +954,7 @@ SampledSpectrum VolPathIntegrator::Li(RayDifferential ray, SampledWavelengths &l
                                       Sampler sampler, ScratchBuffer &scratchBuffer,
                                       VisibleSurface *visibleSurf) const {
     // Declare state variables for volumetric path sampling
-    SampledSpectrum L(0.f), T_hat(1.f), uniPathPDF(1.f), lightPathPDF(1.f);
+    SampledSpectrum L(0.f), beta(1.f), inv_w_u(1.f), inv_w_l(1.f);
     bool specularBounce = false, anyNonSpecularBounces = false;
     int depth = 0;
     Float etaScale = 1;
@@ -964,8 +964,8 @@ SampledSpectrum VolPathIntegrator::Li(RayDifferential ray, SampledWavelengths &l
     while (true) {
         // Sample segment of volumetric scattering path
         PBRT_DBG("%s\n",
-                 StringPrintf("Path tracer depth %d, current L = %s, T_hat = %s\n", depth,
-                              L, T_hat)
+                 StringPrintf("Path tracer depth %d, current L = %s, beta = %s\n", depth,
+                              L, beta)
                      .c_str());
         pstd::optional<ShapeIntersection> si = Intersect(ray);
         if (ray.medium) {
@@ -982,25 +982,22 @@ SampledSpectrum VolPathIntegrator::Li(RayDifferential ray, SampledWavelengths &l
                 [&](Point3f p, MediumProperties mp, SampledSpectrum sigma_maj,
                     SampledSpectrum T_maj) {
                     // Handle medium scattering event for ray
-                    if (!T_hat) {
+                    if (!beta) {
                         terminated = true;
                         return false;
                     }
                     ++volumeInteractions;
-                    // Rescale path throughput and PDFs if necessary
-                    Rescale(&T_hat, &uniPathPDF, &lightPathPDF);
 
                     // Add emission from medium scattering event
                     if (depth < maxDepth && mp.Le) {
-                        // Compute $\hat{P}$ at new path vertex
-                        SampledSpectrum P_hat = T_hat * T_maj * mp.sigma_a * mp.Le;
+                        Float pdf = sigma_maj[0] * T_maj[0];
+                        SampledSpectrum betap = beta * T_maj * mp.sigma_a / pdf;
 
-                        // Compute PDF for absorption at path vertex
-                        SampledSpectrum emitPDF = uniPathPDF * sigma_maj * T_maj;
+                        SampledSpectrum inv_w_e = inv_w_u * sigma_maj * T_maj / pdf;
 
                         // Update _L_ for medium emission
-                        if (emitPDF)
-                            L += P_hat / emitPDF.Average();
+                        if (inv_w_e)
+                            L += betap * mp.Le / inv_w_e.Average();
                     }
 
                     // Compute medium event probabilities for interaction
@@ -1025,14 +1022,15 @@ SampledSpectrum VolPathIntegrator::Li(RayDifferential ray, SampledWavelengths &l
                             return false;
                         }
 
-                        T_hat *= T_maj * mp.sigma_s;
-                        uniPathPDF *= T_maj * mp.sigma_s;
-                        if (T_hat && uniPathPDF) {
+                        Float pdf = T_maj[0] * mp.sigma_s[0];
+                        beta *= T_maj * mp.sigma_s / pdf;
+                        inv_w_u *= T_maj * mp.sigma_s / pdf;
+                        if (beta && inv_w_u) {
                             // Sample direct lighting at volume scattering event
                             MediumInteraction intr(p, -ray.d, ray.time, ray.medium,
                                                    mp.phase);
-                            L += SampleLd(intr, nullptr, lambda, sampler, T_hat,
-                                          uniPathPDF);
+                            L += SampleLd(intr, nullptr, lambda, sampler, beta,
+                                          inv_w_u);
 
                             // Sample new direction at real scattering event
                             Point2f u = sampler.Get2D();
@@ -1042,9 +1040,8 @@ SampledSpectrum VolPathIntegrator::Li(RayDifferential ray, SampledWavelengths &l
                                 terminated = true;
                             else {
                                 // Update ray path state for indirect volume scattering
-                                T_hat *= ps->p;
-                                lightPathPDF = uniPathPDF;
-                                uniPathPDF *= ps->pdf;
+                                beta *= ps->p / ps->pdf;
+                                inv_w_l = inv_w_u / ps->pdf;
                                 prevIntrContext = LightSampleContext(intr);
                                 scattered = true;
                                 ray.o = p;
@@ -1059,21 +1056,21 @@ SampledSpectrum VolPathIntegrator::Li(RayDifferential ray, SampledWavelengths &l
                         // Handle null scattering along ray path
                         SampledSpectrum sigma_n =
                             ClampZero(sigma_maj - mp.sigma_a - mp.sigma_s);
-                        T_hat *= T_maj * sigma_n;
-                        uniPathPDF *= T_maj * sigma_n;
-                        lightPathPDF *= T_maj * sigma_maj;
-                        Rescale(&T_hat, &uniPathPDF, &lightPathPDF);
-                        return T_hat && uniPathPDF;
+                        Float pdf = T_maj[0] * sigma_n[0];
+                        beta *= T_maj * sigma_n / pdf;
+                        inv_w_u *= T_maj * sigma_n / pdf;
+                        inv_w_l *= T_maj * sigma_maj / pdf;
+                        return beta && inv_w_u;
                     }
                 });
             // Handle terminated, scattered, and unscattered rays after medium sampling
-            if (terminated || !T_hat || !uniPathPDF)
+            if (terminated || !beta || !inv_w_u)
                 return L;
             if (scattered)
                 continue;
-            T_hat *= T_maj;
-            uniPathPDF *= T_maj;
-            lightPathPDF *= T_maj;
+            beta *= T_maj / T_maj[0];
+            inv_w_u *= T_maj / T_maj[0];
+            inv_w_l *= T_maj / T_maj[0];
         }
         // Handle surviving unscattered rays
         // Add emitted light at volume path vertex or from the environment
@@ -1082,13 +1079,13 @@ SampledSpectrum VolPathIntegrator::Li(RayDifferential ray, SampledWavelengths &l
             for (const auto &light : infiniteLights) {
                 if (SampledSpectrum Le = light.Le(ray, lambda); Le) {
                     if (depth == 0 || specularBounce)
-                        L += T_hat * Le / uniPathPDF.Average();
+                        L += beta * Le / inv_w_u.Average();
                     else {
                         // Add infinite light contribution using both PDFs with MIS
                         Float lightPDF = lightSampler.PMF(prevIntrContext, light) *
                                          light.PDF_Li(prevIntrContext, ray.d, true);
-                        lightPathPDF *= lightPDF;
-                        L += T_hat * Le / (uniPathPDF + lightPathPDF).Average();
+                        inv_w_l *= lightPDF;
+                        L += beta * Le / (inv_w_u + inv_w_l).Average();
                     }
                 }
             }
@@ -1099,14 +1096,14 @@ SampledSpectrum VolPathIntegrator::Li(RayDifferential ray, SampledWavelengths &l
         if (SampledSpectrum Le = isect.Le(-ray.d, lambda); Le) {
             // Add contribution of emission from intersected surface
             if (depth == 0 || specularBounce)
-                L += T_hat * Le / uniPathPDF.Average();
+                L += beta * Le / inv_w_u.Average();
             else {
                 // Add surface light contribution using both PDFs with MIS
                 Light areaLight(isect.areaLight);
                 Float lightPDF = lightSampler.PMF(prevIntrContext, areaLight) *
                                  areaLight.PDF_Li(prevIntrContext, ray.d, true);
-                lightPathPDF *= lightPDF;
-                L += T_hat * Le / (uniPathPDF + lightPathPDF).Average();
+                inv_w_l *= lightPDF;
+                L += beta * Le / (inv_w_u + inv_w_l).Average();
             }
         }
 
@@ -1154,7 +1151,7 @@ SampledSpectrum VolPathIntegrator::Li(RayDifferential ray, SampledWavelengths &l
 
         // Sample illumination from lights to find attenuated path contribution
         if (IsNonSpecular(bsdf.Flags())) {
-            L += SampleLd(isect, &bsdf, lambda, sampler, T_hat, uniPathPDF);
+            L += SampleLd(isect, &bsdf, lambda, sampler, beta, inv_w_u);
             DCHECK(IsInf(L.y(lambda)) == false);
         }
         prevIntrContext = LightSampleContext(isect);
@@ -1166,20 +1163,16 @@ SampledSpectrum VolPathIntegrator::Li(RayDifferential ray, SampledWavelengths &l
         if (!bs)
             break;
         // Update _T_hat_ and PDFs for BSDF scattering
-        T_hat *= bs->f * AbsDot(bs->wi, isect.shading.n);
-        lightPathPDF = uniPathPDF;
-        if (bs->pdfIsProportional) {
-            Float pdf = bsdf.PDF(wo, bs->wi);
-            T_hat *= pdf / bs->pdf;
-            uniPathPDF *= pdf;
-        } else
-            uniPathPDF *= bs->pdf;
-        Rescale(&T_hat, &uniPathPDF, &lightPathPDF);
+        beta *= bs->f * AbsDot(bs->wi, isect.shading.n) / bs->pdf;
+        if (bs->pdfIsProportional)
+            inv_w_l = inv_w_u / bsdf.PDF(wo, bs->wi);
+        else
+            inv_w_l = inv_w_u / bs->pdf;
 
-        PBRT_DBG("%s\n", StringPrintf("Sampled BSDF, f = %s, pdf = %f -> T_hat = %s",
-                                      bs->f, bs->pdf, T_hat)
+        PBRT_DBG("%s\n", StringPrintf("Sampled BSDF, f = %s, pdf = %f -> beta = %s",
+                                      bs->f, bs->pdf, beta)
                              .c_str());
-        DCHECK(IsInf(T_hat.y(lambda)) == false);
+        DCHECK(IsInf(beta.y(lambda)) == false);
         // Update volumetric integrator path state after surface scattering
         specularBounce = bs->IsSpecular();
         anyNonSpecularBounces |= !bs->IsSpecular();
@@ -1225,8 +1218,10 @@ SampledSpectrum VolPathIntegrator::Li(RayDifferential ray, SampledWavelengths &l
                 break;
 
             // Update path state for subsurface scattering
-            T_hat *= bssrdfSample.Sp;
-            uniPathPDF *= interactionSampler.SampleProbability() * bssrdfSample.pdf;
+            Float pdf = interactionSampler.SampleProbability() * bssrdfSample.pdf[0];
+            beta *= bssrdfSample.Sp / pdf;
+            inv_w_u *= bssrdfSample.pdf / bssrdfSample.pdf[0];
+
             SurfaceInteraction pi = ssi;
             pi.wo = bssrdfSample.wo;
             prevIntrContext = LightSampleContext(pi);
@@ -1240,26 +1235,25 @@ SampledSpectrum VolPathIntegrator::Li(RayDifferential ray, SampledWavelengths &l
                 ++totalBSDFs;
 
             // Account for attenuated direct illumination subsurface scattering
-            L += SampleLd(pi, &Sw, lambda, sampler, T_hat, uniPathPDF);
+            L += SampleLd(pi, &Sw, lambda, sampler, beta, inv_w_u);
 
             // Sample ray for indirect subsurface scattering
             Float u = sampler.Get1D();
             pstd::optional<BSDFSample> bs = Sw.Sample_f(pi.wo, u, sampler.Get2D());
             if (!bs)
                 break;
-            T_hat *= bs->f * AbsDot(bs->wi, pi.shading.n);
-            lightPathPDF = uniPathPDF;
-            uniPathPDF *= bs->pdf;
+            beta *= bs->f * AbsDot(bs->wi, pi.shading.n) / bs->pdf;
+            inv_w_l = inv_w_u / bs->pdf;
             // Don't increment depth this time...
-            DCHECK(!IsInf(T_hat.y(lambda)));
+            DCHECK(!IsInf(beta.y(lambda)));
             specularBounce = bs->IsSpecular();
             ray = RayDifferential(pi.SpawnRay(bs->wi));
         }
 
         // Possibly terminate volumetric path with Russian roulette
-        if (!T_hat)
+        if (!beta)
             break;
-        SampledSpectrum rrBeta = T_hat * etaScale / uniPathPDF.Average();
+        SampledSpectrum rrBeta = beta * etaScale / inv_w_u.Average();
         Float uRR = sampler.Get1D();
         PBRT_DBG("%s\n",
                  StringPrintf("etaScale %f -> rrBeta %s", etaScale, rrBeta).c_str());
@@ -1267,8 +1261,7 @@ SampledSpectrum VolPathIntegrator::Li(RayDifferential ray, SampledWavelengths &l
             Float q = std::max<Float>(0, 1 - rrBeta.MaxComponentValue());
             if (uRR < q)
                 break;
-            uniPathPDF *= 1 - q;
-            lightPathPDF *= 1 - q;
+            beta /= 1 - q;
         }
     }
     return L;
@@ -1276,8 +1269,8 @@ SampledSpectrum VolPathIntegrator::Li(RayDifferential ray, SampledWavelengths &l
 
 SampledSpectrum VolPathIntegrator::SampleLd(const Interaction &intr, const BSDF *bsdf,
                                             SampledWavelengths &lambda, Sampler sampler,
-                                            SampledSpectrum T_hat,
-                                            SampledSpectrum pathPDF) const {
+                                            SampledSpectrum beta,
+                                            SampledSpectrum inv_w_p) const {
     // Estimate light-sampled direct illumination at _intr_
     // Initialize _LightSampleContext_ for volumetric light sampling
     LightSampleContext ctx;
@@ -1329,7 +1322,7 @@ SampledSpectrum VolPathIntegrator::SampleLd(const Interaction &intr, const BSDF
 
     // Declare path state variables for ray to light source
     Ray lightRay = intr.SpawnRayTo(ls->pLight);
-    SampledSpectrum T_ray(1.f), lightPathPDF(1.f), uniPathPDF(1.f);
+    SampledSpectrum T_ray(1.f), inv_w_l(1.f), inv_w_u(1.f);
     RNG rng(Hash(lightRay.o), Hash(lightRay.d));
 
     while (lightRay.d != Vector3f(0, 0, 0)) {
@@ -1351,31 +1344,29 @@ SampledSpectrum VolPathIntegrator::SampleLd(const Interaction &intr, const BSDF
                     // Update _T_ray_ and PDFs using ratio-tracking estimator
                     SampledSpectrum sigma_n =
                         ClampZero(sigma_maj - mp.sigma_a - mp.sigma_s);
-                    T_ray *= T_maj * sigma_n;
-                    lightPathPDF *= T_maj * sigma_maj;
-                    uniPathPDF *= T_maj * sigma_n;
+                    Float pdf = T_maj[0] * sigma_maj[0];
+                    T_ray *= T_maj * sigma_n / pdf;
+                    inv_w_l *= T_maj * sigma_maj / pdf;
+                    inv_w_u *= T_maj * sigma_n / pdf;
 
                     // Possibly terminate transmittance computation using Russian roulette
-                    SampledSpectrum Tr = T_ray / (lightPathPDF + uniPathPDF).Average();
+                    SampledSpectrum Tr = T_ray / (inv_w_l + inv_w_u).Average();
                     if (Tr.MaxComponentValue() < 0.05f) {
                         Float q = 0.75f;
                         if (rng.Uniform<Float>() < q)
                             T_ray = SampledSpectrum(0.);
-                        else {
-                            lightPathPDF *= 1 - q;
-                            uniPathPDF *= 1 - q;
-                        }
+                        else
+                            T_ray /= 1 - q;
                     }
 
                     if (!T_ray)
                         return false;
-                    Rescale(&T_ray, &lightPathPDF, &uniPathPDF);
                     return true;
                 });
             // Update transmittance estimate for final segment
-            T_ray *= T_maj;
-            lightPathPDF *= T_maj;
-            uniPathPDF *= T_maj;
+            T_ray *= T_maj / T_maj[0];
+            inv_w_l *= T_maj / T_maj[0];
+            inv_w_u *= T_maj / T_maj[0];
         }
 
         // Generate next ray segment or return final transmittance
@@ -1386,12 +1377,12 @@ SampledSpectrum VolPathIntegrator::SampleLd(const Interaction &intr, const BSDF
         lightRay = si->intr.SpawnRayTo(ls->pLight);
     }
     // Return path contribution function estimate for direct lighting
-    lightPathPDF *= pathPDF * lightPDF;
-    uniPathPDF *= pathPDF * scatterPDF;
+    inv_w_l *= inv_w_p * lightPDF;
+    inv_w_u *= inv_w_p * scatterPDF;
     if (IsDeltaLight(light.Type()))
-        return T_hat * f_hat * T_ray * ls->L / lightPathPDF.Average();
+        return beta * f_hat * T_ray * ls->L / inv_w_l.Average();
     else
-        return T_hat * f_hat * T_ray * ls->L / (lightPathPDF + uniPathPDF).Average();
+        return beta * f_hat * T_ray * ls->L / (inv_w_l + inv_w_u).Average();
 }
 
 std::string VolPathIntegrator::ToString() const {

src/pbrt/cpu/integrators.h

@@ -284,27 +284,7 @@ class VolPathIntegrator : public RayIntegrator {
     // VolPathIntegrator Private Methods
     SampledSpectrum SampleLd(const Interaction &intr, const BSDF *bsdf,
                              SampledWavelengths &lambda, Sampler sampler,
-                             SampledSpectrum T_hat, SampledSpectrum pathPDF) const;
-
-    static void Rescale(SampledSpectrum *T_hat, SampledSpectrum *uniPathPDF,
-                        SampledSpectrum *lightPathPDF) {
-        if (T_hat->MaxComponentValue() > 0x1p24f ||
-            lightPathPDF->MaxComponentValue() > 0x1p24f ||
-            uniPathPDF->MaxComponentValue() > 0x1p24f) {
-            // Downscale _T_hat_, _lightPathPDF_, and _uniPathPDF_
-            *T_hat *= 1.f / 0x1p24f;
-            *lightPathPDF *= 1.f / 0x1p24f;
-            *uniPathPDF *= 1.f / 0x1p24f;
-        }
-        // Upscale _T_hat_, _lightPathPDF_, and _uniPathPDF_ if necessary
-        if (T_hat->MaxComponentValue() < 0x1p-24f ||
-            lightPathPDF->MaxComponentValue() < 0x1p-24f ||
-            uniPathPDF->MaxComponentValue() < 0x1p-24f) {
-            *T_hat *= 0x1p24f;
-            *lightPathPDF *= 0x1p24f;
-            *uniPathPDF *= 0x1p24f;
-        }
-    }
+                             SampledSpectrum beta, SampledSpectrum inv_w_u) const;
 
     // VolPathIntegrator Private Members
     int maxDepth;