Update from book source. No functional changes.

src/pbrt/bssrdf.cpp

@@ -50,7 +50,7 @@ Float BeamDiffusionMS(Float sigma_s, Float sigma_a, Float g, Float eta, Float r)
 
     for (int i = 0; i < nSamples; ++i) {
         // Sample real point source depth $\depthreal$
-        Float zr = -std::log(1 - (i + 0.5f) / nSamples) / sigmap_t;
+        Float zr = SampleExponential((i + 0.5f) / nSamples, sigmap_t);
 
         // Evaluate dipole integrand $E_{\roman{d}}$ at $\depthreal$ and add to _Ed_
         Float zv = -zr + 2 * ze;
@@ -83,7 +83,7 @@ Float BeamDiffusionSS(Float sigma_s, Float sigma_a, Float g, Float eta, Float r)
     const int nSamples = 100;
     for (int i = 0; i < nSamples; ++i) {
         // Evaluate single-scattering integrand and add to _Ess_
-        Float ti = tCrit - std::log(1 - (i + 0.5f) / nSamples) / sigma_t;
+        Float ti = tCrit + SampleExponential((i + 0.5f) / nSamples, sigma_t);
         // Determine length $d$ of connecting segment and $\cos\theta_\roman{o}$
         Float d = std::sqrt(Sqr(r) + Sqr(ti));
         Float cosTheta_o = ti / d;

src/pbrt/bssrdf.h

@@ -233,7 +233,7 @@ class TabulatedBSSRDF {
 
     PBRT_CPU_GPU
     SampledSpectrum PDF_Sp(Point3f pi, Normal3f ni) const {
-        // Express $\pti-\pto$ and $\bold{n}_i$ with respect to local coordinates at
+        // Express $\pti-\pto$ and $\N{}_\roman{i}$ with respect to local coordinates at
         // $\pto$
         Vector3f d = pi - po;
         Frame f = Frame::FromZ(ns);

src/pbrt/cpu/integrators.cpp

@@ -847,12 +847,12 @@ SampledSpectrum SimpleVolPathIntegrator::Li(RayDifferential ray,
         pstd::optional<ShapeIntersection> si = Intersect(ray);
         bool scattered = false, terminated = false;
         if (ray.medium) {
-            // Initialize RNG for delta tracking
+            // Initialize _RNG_ for sampling the majorant transmittance
             uint64_t hash0 = Hash(sampler.Get1D());
             uint64_t hash1 = Hash(sampler.Get1D());
             RNG rng(hash0, hash1);
 
-            // Sample medium scattering using delta tracking
+            // Sample medium using delta tracking
             Float tMax = si ? si->tHit : Infinity;
             Float u = sampler.Get1D();
             Float uMode = sampler.Get1D();
@@ -867,21 +867,20 @@ SampledSpectrum SimpleVolPathIntegrator::Li(RayDifferential ray,
                             // Randomly sample medium scattering event for delta tracking
                             int mode = SampleDiscrete({pAbsorb, pScatter, pNull}, uMode);
                             if (mode == 0) {
-                                // Handle absorption event for delta tracking
+                                // Handle absorption event for medium sample
                                 L += beta * mp.Le;
                                 terminated = true;
                                 return false;
 
                             } else if (mode == 1) {
-                                // Handle regular scattering event for delta tracking
+                                // Handle regular scattering event for medium sample
                                 // Stop path sampling if maximum depth has been reached
                                 if (depth++ >= maxDepth) {
                                     terminated = true;
                                     return false;
                                 }
 
-                                // Sample phase function for delta tracking scattering
-                                // event
+                                // Sample phase function for medium scattering event
                                 Point2f u{rng.Uniform<Float>(), rng.Uniform<Float>()};
                                 pstd::optional<PhaseFunctionSample> ps =
                                     mp.phase.Sample_p(-ray.d, u);
@@ -890,7 +889,7 @@ SampledSpectrum SimpleVolPathIntegrator::Li(RayDifferential ray,
                                     return false;
                                 }
 
-                                // Update state for recursive $L_\roman{i}$ evaluation
+                                // Update state for recursive evaluation of $L_\roman{i}$
                                 beta *= ps->p / ps->pdf;
                                 ray.o = p;
                                 ray.d = ps->wi;
@@ -898,7 +897,7 @@ SampledSpectrum SimpleVolPathIntegrator::Li(RayDifferential ray,
                                 return false;
 
                             } else {
-                                // Handle null scattering event for delta tracking
+                                // Handle null scattering event for medium sample
                                 uMode = rng.Uniform<Float>();
                                 return true;
                             }
@@ -909,7 +908,7 @@ SampledSpectrum SimpleVolPathIntegrator::Li(RayDifferential ray,
             return L;
         if (scattered)
             continue;
-        // Add emission to unscattered ray
+        // Add emission to surviving ray
         if (si)
             L += beta * si->intr.Le(-ray.d, lambda);
         else {
@@ -918,7 +917,7 @@ SampledSpectrum SimpleVolPathIntegrator::Li(RayDifferential ray,
             return L;
         }
 
-        // Handle surface intersection along delta tracking path
+        // Handle surface intersection along ray path
         BSDF bsdf = si->intr.GetBSDF(ray, lambda, camera, buf, sampler);
         if (!bsdf)
             si->intr.SkipIntersection(&ray, si->tHit);
@@ -973,7 +972,7 @@ SampledSpectrum VolPathIntegrator::Li(RayDifferential ray, SampledWavelengths &l
             // Sample the participating medium
             bool scattered = false, terminated = false;
             Float tMax = si ? si->tHit : Infinity;
-            // Initialize RNG for delta tracking
+            // Initialize _RNG_ for sampling the majorant transmittance
             uint64_t hash0 = Hash(sampler.Get1D());
             uint64_t hash1 = Hash(sampler.Get1D());
             RNG rng(hash0, hash1);
@@ -989,14 +988,14 @@ SampledSpectrum VolPathIntegrator::Li(RayDifferential ray, SampledWavelengths &l
                     }
                     ++volumeInteractions;
                     // Rescale path throughput and PDFs if necessary
-                    Rescale(T_hat, uniPathPDF, lightPathPDF);
+                    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;
 
-                        // Compute PDF for absorption event at path vertex
+                        // Compute PDF for absorption at path vertex
                         SampledSpectrum emitPDF = uniPathPDF * sigma_maj * T_maj;
 
                         // Update _L_ for medium emission
@@ -1035,21 +1034,20 @@ SampledSpectrum VolPathIntegrator::Li(RayDifferential ray, SampledWavelengths &l
                         Point2f u = sampler.Get2D();
                         pstd::optional<PhaseFunctionSample> ps =
                             intr.phase.Sample_p(-ray.d, u);
-                        if (!ps || ps->pdf == 0) {
+                        if (!ps || ps->pdf == 0)
                             terminated = true;
-                            return false;
+                        else {
+                            // Update ray path state for indirect volume scattering
+                            T_hat *= ps->p;
+                            lightPathPDF = uniPathPDF;
+                            uniPathPDF *= ps->pdf;
+                            prevIntrContext = LightSampleContext(intr);
+                            scattered = true;
+                            ray.o = p;
+                            ray.d = ps->wi;
+                            specularBounce = false;
+                            anyNonSpecularBounces = true;
                         }
-                        // Update ray path state for indirect volume scattering
-                        T_hat *= ps->p;
-                        lightPathPDF = uniPathPDF;
-                        uniPathPDF *= ps->pdf;
-                        prevIntrContext = LightSampleContext(intr);
-                        scattered = true;
-                        ray.o = p;
-                        ray.d = ps->wi;
-                        specularBounce = false;
-                        anyNonSpecularBounces = true;
-
                         return false;
 
                     } else {
@@ -1059,7 +1057,8 @@ SampledSpectrum VolPathIntegrator::Li(RayDifferential ray, SampledWavelengths &l
                         T_hat *= T_maj * sigma_n;
                         uniPathPDF *= T_maj * sigma_n;
                         lightPathPDF *= T_maj * sigma_maj;
-                        Rescale(T_hat, uniPathPDF, lightPathPDF);
+                        Rescale(&T_hat, &uniPathPDF, &lightPathPDF);
+
                         return true;
                     }
                 });
@@ -1171,7 +1170,7 @@ SampledSpectrum VolPathIntegrator::Li(RayDifferential ray, SampledWavelengths &l
             uniPathPDF *= pdf;
         } else
             uniPathPDF *= bs->pdf;
-        Rescale(T_hat, uniPathPDF, lightPathPDF);
+        Rescale(&T_hat, &uniPathPDF, &lightPathPDF);
 
         PBRT_DBG("%s\n", StringPrintf("Sampled BSDF, f = %s, pdf = %f -> T_hat = %s",
                                       bs->f, bs->pdf, T_hat)
@@ -1214,18 +1213,19 @@ SampledSpectrum VolPathIntegrator::Li(RayDifferential ray, SampledWavelengths &l
             if (!interactionSampler.HasSample())
                 break;
 
-            // Convert probe intersection to _BSSRDFSample_ and update path state
+            // Convert probe intersection to _BSSRDFSample_
             SubsurfaceInteraction ssi = interactionSampler.GetSample();
             BSSRDFSample bssrdfSample =
                 bssrdf.ProbeIntersectionToSample(ssi, scratchBuffer);
             if (!bssrdfSample.Sp || !bssrdfSample.pdf)
                 break;
+
+            // Update path state for subsurface scattering
             T_hat *= bssrdfSample.Sp;
             uniPathPDF *= interactionSampler.SampleProbability() * bssrdfSample.pdf;
             SurfaceInteraction pi = ssi;
             pi.wo = bssrdfSample.wo;
             prevIntrContext = LightSampleContext(pi);
-
             // Possibly regularize subsurface BSDF
             BSDF &Sw = bssrdfSample.Sw;
             anyNonSpecularBounces = true;
@@ -1235,7 +1235,7 @@ SampledSpectrum VolPathIntegrator::Li(RayDifferential ray, SampledWavelengths &l
             } else
                 ++totalBSDFs;
 
-            // Account for attenuated direct subsurface scattering
+            // Account for attenuated direct illumination subsurface scattering
             L += SampleLd(pi, &Sw, lambda, sampler, T_hat, uniPathPDF);
 
             // Sample ray for indirect subsurface scattering
@@ -1365,10 +1365,10 @@ SampledSpectrum VolPathIntegrator::SampleLd(const Interaction &intr, const BSDF
 
                     if (!T_ray)
                         return false;
-                    Rescale(T_ray, lightPathPDF, uniPathPDF);
+                    Rescale(&T_ray, &lightPathPDF, &uniPathPDF);
                     return true;
                 });
-            // Update transmittance estimate for final unsampled segment
+            // Update transmittance estimate for final segment
             T_ray *= T_maj;
             lightPathPDF *= T_maj;
             uniPathPDF *= T_maj;

src/pbrt/cpu/integrators.h

@@ -286,23 +286,23 @@ class VolPathIntegrator : public RayIntegrator {
                              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) {
+    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;
+            *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;
+        if (T_hat->MaxComponentValue() < 0x1p-24f ||
+            lightPathPDF->MaxComponentValue() < 0x1p-24f ||
+            uniPathPDF->MaxComponentValue() < 0x1p-24f) {
+            *T_hat *= 0x1p24f;
+            *lightPathPDF *= 0x1p24f;
+            *uniPathPDF *= 0x1p24f;
         }
     }
 

src/pbrt/media.h

@@ -737,9 +737,19 @@ inline MediumProperties Medium::SamplePoint(Point3f p,
 }
 
 // Medium Sampling Function Definitions
+template <typename F>
+PBRT_CPU_GPU SampledSpectrum SampleT_maj(Ray ray, Float tMax, Float u, RNG &rng,
+                                         const SampledWavelengths &lambda, F callback) {
+    auto sample = [&](auto medium) {
+        using M = typename std::remove_reference_t<decltype(*medium)>;
+        return SampleT_maj<M>(ray, tMax, u, rng, lambda, callback);
+    };
+    return ray.medium.Dispatch(sample);
+}
+
 template <typename ConcreteMedium, typename F>
 PBRT_CPU_GPU SampledSpectrum SampleT_maj(Ray ray, Float tMax, Float u, RNG &rng,
-                                         const SampledWavelengths &lambda, F func) {
+                                         const SampledWavelengths &lambda, F callback) {
     // Normalize ray direction and update _tMax_ accordingly
     tMax *= Length(ray.d);
     ray.d = Normalize(ray.d);
@@ -749,7 +759,7 @@ PBRT_CPU_GPU SampledSpectrum SampleT_maj(Ray ray, Float tMax, Float u, RNG &rng,
     typename ConcreteMedium::MajorantIterator iter;
     medium->SampleRay(ray, tMax, lambda, &iter);
 
-    // Generate ray majorant samples until termination event
+    // Generate ray majorant samples until termination
     SampledSpectrum T_maj(1.f);
     bool done = false;
     while (!done) {
@@ -768,16 +778,17 @@ PBRT_CPU_GPU SampledSpectrum SampleT_maj(Ray ray, Float tMax, Float u, RNG &rng,
             continue;
         }
 
+        // Generate samples along current majorant segment
         Float tMin = seg->tMin;
         while (true) {
             // Try to generate sample along current majorant segment
             Float t = tMin + SampleExponential(u, seg->sigma_maj[0]);
             u = rng.Uniform<Float>();
             if (t < seg->tMax) {
-                // Call callback function for sample along majorant segment
+                // Call callback function for sample within segment
                 T_maj *= FastExp(-(t - tMin) * seg->sigma_maj);
                 MediumProperties mp = medium->SamplePoint(ray(t), lambda);
-                if (!func(ray(t), mp, seg->sigma_maj, T_maj)) {
+                if (!callback(ray(t), mp, seg->sigma_maj, T_maj)) {
                     // Returning out of doubly-nested while loop is not as good perf. wise
                     // on the GPU vs using "done" here.
                     done = true;
@@ -801,16 +812,6 @@ PBRT_CPU_GPU SampledSpectrum SampleT_maj(Ray ray, Float tMax, Float u, RNG &rng,
     return SampledSpectrum(1.f);
 }
 
-template <typename F>
-PBRT_CPU_GPU SampledSpectrum SampleT_maj(Ray ray, Float tMax, Float u, RNG &rng,
-                                         const SampledWavelengths &lambda, F func) {
-    auto sample = [&](auto medium) {
-        using Medium = typename std::remove_reference_t<decltype(*medium)>;
-        return SampleT_maj<Medium>(ray, tMax, u, rng, lambda, func);
-    };
-    return ray.medium.Dispatch(sample);
-}
-
 inline RayMajorantIterator Medium::SampleRay(Ray ray, Float tMax,
                                              const SampledWavelengths &lambda,
                                              ScratchBuffer &buf) const {

src/pbrt/wavefront/media.cpp

@@ -8,27 +8,21 @@
 
 namespace pbrt {
 
-// It's not unususal for these values to have very large or very small
-// magnitudes after multiple (null) scattering events, even though in the
-// end ratios like T_hat/uniPathPDF are generally around 1.  To avoid overflow,
-// we rescale all three of them by the same factor when they become large.
 PBRT_CPU_GPU
-static inline void rescale(SampledSpectrum &T_hat, SampledSpectrum &lightPathPDF,
-                           SampledSpectrum &uniPathPDF) {
-    // Note that no precision is lost in the rescaling since we're always
-    // multiplying by an exact power of 2.
-    if (T_hat.MaxComponentValue() > 0x1p24f ||
-        lightPathPDF.MaxComponentValue() > 0x1p24f ||
-        uniPathPDF.MaxComponentValue() > 0x1p24f) {
-        T_hat *= 1.f / 0x1p24f;
-        lightPathPDF *= 1.f / 0x1p24f;
-        uniPathPDF *= 1.f / 0x1p24f;
-    } else if (T_hat.MaxComponentValue() < 0x1p-24f ||
-               lightPathPDF.MaxComponentValue() < 0x1p-24f ||
-               uniPathPDF.MaxComponentValue() < 0x1p-24f) {
-        T_hat *= 0x1p24f;
-        lightPathPDF *= 0x1p24f;
-        uniPathPDF *= 0x1p24f;
+static inline void rescale(SampledSpectrum *T_hat, SampledSpectrum *lightPathPDF,
+                           SampledSpectrum *uniPathPDF) {
+    if (T_hat->MaxComponentValue() > 0x1p24f ||
+        lightPathPDF->MaxComponentValue() > 0x1p24f ||
+        uniPathPDF->MaxComponentValue() > 0x1p24f) {
+        *T_hat *= 1.f / 0x1p24f;
+        *lightPathPDF *= 1.f / 0x1p24f;
+        *uniPathPDF *= 1.f / 0x1p24f;
+    } else if (T_hat->MaxComponentValue() < 0x1p-24f ||
+               lightPathPDF->MaxComponentValue() < 0x1p-24f ||
+               uniPathPDF->MaxComponentValue() < 0x1p-24f) {
+        *T_hat *= 0x1p24f;
+        *lightPathPDF *= 0x1p24f;
+        *uniPathPDF *= 0x1p24f;
     }
 }
 
@@ -86,7 +80,7 @@ void WavefrontPathIntegrator::SampleMediumInteraction(int wavefrontDepth) {
                 ray, tMax, uDist, rng, lambda,
                 [&](Point3f p, MediumProperties mp, SampledSpectrum sigma_maj,
                     SampledSpectrum T_maj) {
-                    rescale(T_hat, uniPathPDF, lightPathPDF);
+                    rescale(&T_hat, &uniPathPDF, &lightPathPDF);
 
                     PBRT_DBG("Medium event T_maj %f %f %f %f sigma_a %f %f %f %f sigma_s "
                              "%f %f "

src/pbrt/wavefront/surfscatter.cpp

@@ -22,20 +22,20 @@
 namespace pbrt {
 
 PBRT_CPU_GPU
-static inline void rescale(SampledSpectrum &T_hat, SampledSpectrum &lightPathPDF,
-                           SampledSpectrum &uniPathPDF) {
-    if (T_hat.MaxComponentValue() > 0x1p24f ||
-        lightPathPDF.MaxComponentValue() > 0x1p24f ||
-        uniPathPDF.MaxComponentValue() > 0x1p24f) {
-        T_hat *= 1.f / 0x1p24f;
-        lightPathPDF *= 1.f / 0x1p24f;
-        uniPathPDF *= 1.f / 0x1p24f;
-    } else if (T_hat.MaxComponentValue() < 0x1p-24f ||
-               lightPathPDF.MaxComponentValue() < 0x1p-24f ||
-               uniPathPDF.MaxComponentValue() < 0x1p-24f) {
-        T_hat *= 0x1p24f;
-        lightPathPDF *= 0x1p24f;
-        uniPathPDF *= 0x1p24f;
+static inline void rescale(SampledSpectrum *T_hat, SampledSpectrum *lightPathPDF,
+                           SampledSpectrum *uniPathPDF) {
+    if (T_hat->MaxComponentValue() > 0x1p24f ||
+        lightPathPDF->MaxComponentValue() > 0x1p24f ||
+        uniPathPDF->MaxComponentValue() > 0x1p24f) {
+        *T_hat *= 1.f / 0x1p24f;
+        *lightPathPDF *= 1.f / 0x1p24f;
+        *uniPathPDF *= 1.f / 0x1p24f;
+    } else if (T_hat->MaxComponentValue() < 0x1p-24f ||
+               lightPathPDF->MaxComponentValue() < 0x1p-24f ||
+               uniPathPDF->MaxComponentValue() < 0x1p-24f) {
+        *T_hat *= 0x1p24f;
+        *lightPathPDF *= 0x1p24f;
+        *uniPathPDF *= 0x1p24f;
     }
 }
 
@@ -199,7 +199,7 @@ void WavefrontPathIntegrator::EvaluateMaterialAndBSDF(MaterialEvalQueue *evalQue
                 } else
                     uniPathPDF *= bsdfSample->pdf;
 
-                rescale(T_hat, uniPathPDF, lightPathPDF);
+                rescale(&T_hat, &uniPathPDF, &lightPathPDF);
                 // Update _etaScale_ accounting for BSDF scattering
                 Float etaScale = w.etaScale;
                 if (bsdfSample->IsTransmission())