Implement non-exponential transport

src/core/grids/Grid.hpp

@@ -20,9 +20,8 @@ class Grid : public JsonSerializable
     virtual Box3f bounds() const;
 
     virtual float density(Vec3f p) const = 0;
-    virtual Vec3f transmittance(PathSampleGenerator &sampler, Vec3f p, Vec3f w, float t0, float t1, Vec3f sigmaT) const = 0;
-    virtual Vec2f inverseOpticalDepth(PathSampleGenerator &sampler, Vec3f p, Vec3f w, float t0, float t1,
-            float sigmaT, float xi) const = 0;
+    virtual float opticalDepth(PathSampleGenerator &sampler, Vec3f p, Vec3f w, float t0, float t1) const = 0;
+    virtual Vec2f inverseOpticalDepth(PathSampleGenerator &sampler, Vec3f p, Vec3f w, float t0, float t1, float xi) const = 0;
 };
 
 }

src/core/grids/VdbGrid.cpp

@@ -237,7 +237,7 @@ float VdbGrid::density(Vec3f p) const
     return gridAt(_grid->tree(), p);
 }
 
-Vec3f VdbGrid::transmittance(PathSampleGenerator &sampler, Vec3f p, Vec3f w, float t0, float t1, Vec3f sigmaT) const
+float VdbGrid::opticalDepth(PathSampleGenerator &sampler, Vec3f p, Vec3f w, float t0, float t1) const
 {
     auto accessor = _grid->getConstAccessor();
 
@@ -249,7 +249,7 @@ Vec3f VdbGrid::transmittance(PathSampleGenerator &sampler, Vec3f p, Vec3f w, flo
             integral += accessor.getValue(voxel)*(tb - ta);
             return false;
         });
-        return std::exp(-integral*sigmaT);
+        return integral;
     } else if (_integrationMethod == IntegrationMethod::ExactLinear) {
         VdbRaymarcher<openvdb::FloatGrid::TreeType, 3> dda;
 
@@ -261,40 +261,37 @@ Vec3f VdbGrid::transmittance(PathSampleGenerator &sampler, Vec3f p, Vec3f w, flo
             fa = fb;
             return false;
         });
-        return std::exp(-integral*sigmaT);
+        return integral;
     } else if (_integrationMethod == IntegrationMethod::ResidualRatio) {
         VdbRaymarcher<Vec2fGrid::TreeType, 3> dda;
 
         float scale = _supergridSubsample;
         float invScale = 1.0f/scale;
-        sigmaT *= scale;
-
-        float sigmaTc = sigmaT.max();
 
         auto superAccessor =  _superGrid->getConstAccessor();
 
         UniformSampler &generator = sampler.uniformGenerator();
 
         float controlIntegral = 0.0f;
-        Vec3f Tr(1.0f);
+        float Tr = 1.0f;
         dda.march(DdaRay(p*invScale + 0.5f, w), t0*invScale, t1*invScale, superAccessor, [&](openvdb::Coord voxel, float ta, float tb) {
             openvdb::Vec2s v = superAccessor.getValue(voxel);
             float muC = v.x();
             float muR = v.y();
-            muR *= sigmaTc;
+            muR *= scale;
 
             controlIntegral += muC*(tb - ta);
 
             while (true) {
                 ta -= BitManip::normalizedLog(generator.nextI())/muR;
                 if (ta >= tb)
                     break;
-                Tr *= 1.0f - sigmaT*((gridAt(accessor, p + w*ta*scale) - muC)/muR);
+                Tr *= 1.0f - scale*((gridAt(accessor, p + w*ta*scale) - muC)/muR);
             }
 
             return false;
         });
-        return std::exp(-controlIntegral*sigmaT)*Tr;
+        return controlIntegral - std::log(Tr);
     } else {
         float ta = t0;
         float fa = gridAt(accessor, p + w*t0);
@@ -308,12 +305,11 @@ Vec3f VdbGrid::transmittance(PathSampleGenerator &sampler, Vec3f p, Vec3f w, flo
             fa = fb;
             dT = _stepSize;
         } while (ta < t1);
-        return std::exp(-integral*sigmaT);
+        return integral;
     }
 }
 
-Vec2f VdbGrid::inverseOpticalDepth(PathSampleGenerator &sampler, Vec3f p, Vec3f w, float t0, float t1,
-        float sigmaT, float xi) const
+Vec2f VdbGrid::inverseOpticalDepth(PathSampleGenerator &sampler, Vec3f p, Vec3f w, float t0, float t1, float tau) const
 {
     auto accessor = _grid->getConstAccessor();
 
@@ -324,9 +320,9 @@ Vec2f VdbGrid::inverseOpticalDepth(PathSampleGenerator &sampler, Vec3f p, Vec3f
         Vec2f result(t1, 0.0f);
         bool exited = !dda.march(DdaRay(p + 0.5f, w), t0, t1, accessor, [&](openvdb::Coord voxel, float ta, float tb) {
             float v = accessor.getValue(voxel);
-            float delta = v*sigmaT*(tb - ta);
-            if (integral + delta >= xi) {
-                result = Vec2f(ta + (tb - ta)*(xi - integral)/delta, v);
+            float delta = v*(tb - ta);
+            if (integral + delta >= tau) {
+                result = Vec2f(ta + (tb - ta)*(tau - integral)/delta, v);
                 return true;
             }
             integral += delta;
@@ -341,11 +337,11 @@ Vec2f VdbGrid::inverseOpticalDepth(PathSampleGenerator &sampler, Vec3f p, Vec3f
         Vec2f result(t1, 0.0f);
         bool exited = !dda.march(DdaRay(p + 0.5f, w), t0, t1, accessor, [&](openvdb::Coord /*voxel*/, float ta, float tb) {
             float fb = gridAt(accessor, p + tb*w);
-            float delta = (fb + fa)*0.5f*sigmaT*(tb - ta);
-            if (integral + delta >= xi) {
-                float a = (fb - fa)*sigmaT;
-                float b = fa*sigmaT;
-                float c = (integral - xi)/(tb - ta);
+            float delta = (fb + fa)*0.5f*(tb - ta);
+            if (integral + delta >= tau) {
+                float a = (fb - fa);
+                float b = fa;
+                float c = (integral - tau)/(tb - ta);
                 float mantissa = max(b*b - 2.0f*a*c, 0.0f);
                 float x1 = (-b + std::sqrt(mantissa))/a;
                 result = Vec2f(ta + (tb - ta)*x1, fa + (fb - fa)*x1);
@@ -364,11 +360,11 @@ Vec2f VdbGrid::inverseOpticalDepth(PathSampleGenerator &sampler, Vec3f p, Vec3f
         do {
             float tb = min(ta + dT, t1);
             float fb = gridAt(accessor, p + w*tb);
-            float delta = (fa + fb)*sigmaT*0.5f*(tb - ta);
-            if (integral + delta >= xi) {
-                float a = (fb - fa)*sigmaT;
-                float b = fa*sigmaT;
-                float c = (integral - xi)/(tb - ta);
+            float delta = (fa + fb)*0.5f*(tb - ta);
+            if (integral + delta >= tau) {
+                float a = (fb - fa);
+                float b = fa;
+                float c = (integral - tau)/(tb - ta);
                 float mantissa = max(b*b - 2.0f*a*c, 0.0f);
                 float x1 = (-b + std::sqrt(mantissa))/a;
                 return Vec2f(ta + (tb - ta)*x1, fa + (fb - fa)*x1);

src/core/grids/VdbGrid.hpp

@@ -68,9 +68,8 @@ class VdbGrid : public Grid
     virtual Box3f bounds() const override;
 
     float density(Vec3f p) const override;
-    Vec3f transmittance(PathSampleGenerator &sampler, Vec3f p, Vec3f w, float t0, float t1, Vec3f sigmaT) const override;
-    Vec2f inverseOpticalDepth(PathSampleGenerator &sampler, Vec3f p, Vec3f w, float t0, float t1,
-            float sigmaT, float xi) const override;
+    float opticalDepth(PathSampleGenerator &sampler, Vec3f p, Vec3f w, float t0, float t1) const override;
+    Vec2f inverseOpticalDepth(PathSampleGenerator &sampler, Vec3f p, Vec3f w, float t0, float t1, float tau) const override;
 };
 
 }

src/core/integrators/TraceBase.cpp

@@ -108,7 +108,7 @@ inline Vec3f TraceBase::generalizedShadowRayImpl(PathSampleGenerator &sampler,
                 pdfForward *= forward;
                 pdfBackward *= backward;
             } else {
-                throughput *= medium->transmittance(sampler, ray);
+                throughput *= medium->transmittance(sampler, ray, startsOnSurface, endsOnSurface);
             }
         }
         if (info.primitive == nullptr || info.primitive == endCap)
@@ -125,10 +125,11 @@ inline Vec3f TraceBase::generalizedShadowRayImpl(PathSampleGenerator &sampler,
 }
 
 Vec3f TraceBase::generalizedShadowRay(PathSampleGenerator &sampler, Ray &ray, const Medium *medium,
-            const Primitive *endCap, int bounce) const
+            const Primitive *endCap, bool startsOnSurface, bool endsOnSurface, int bounce) const
 {
     float dummyA, dummyB;
-    return generalizedShadowRayImpl<false>(sampler, ray, medium, endCap, bounce, false, false, dummyA, dummyB);
+    return generalizedShadowRayImpl<false>(sampler, ray, medium, endCap, bounce,
+            startsOnSurface, endsOnSurface, dummyA, dummyB);
 }
 
 Vec3f TraceBase::generalizedShadowRayAndPdfs(PathSampleGenerator &sampler, Ray &ray, const Medium *medium,
@@ -147,6 +148,7 @@ Vec3f TraceBase::attenuatedEmission(PathSampleGenerator &sampler,
                          IntersectionTemporary &data,
                          IntersectionInfo &info,
                          int bounce,
+                         bool startsOnSurface,
                          Ray &ray,
                          Vec3f *transmittance)
 {
@@ -162,7 +164,7 @@ Vec3f TraceBase::attenuatedEmission(PathSampleGenerator &sampler,
     info.w = ray.dir();
     light.intersectionInfo(data, info);
 
-    Vec3f shadow = generalizedShadowRay(sampler, ray, medium, &light, bounce);
+    Vec3f shadow = generalizedShadowRay(sampler, ray, medium, &light, startsOnSurface, true, bounce);
     if (transmittance)
         *transmittance = shadow;
     if (shadow == 0.0f)
@@ -192,7 +194,7 @@ bool TraceBase::volumeLensSample(const Camera &camera,
     ray.setPrimaryRay(false);
     ray.setFarT(lensSample.dist);
 
-    Vec3f transmittance = generalizedShadowRay(sampler, ray, medium, nullptr, bounce);
+    Vec3f transmittance = generalizedShadowRay(sampler, ray, medium, nullptr, false, true, bounce);
     if (transmittance == 0.0f)
         return false;
 
@@ -231,7 +233,7 @@ bool TraceBase::surfaceLensSample(const Camera &camera,
     ray.setPrimaryRay(false);
     ray.setFarT(sample.dist);
 
-    Vec3f transmittance = generalizedShadowRay(*event.sampler, ray, medium, nullptr, bounce);
+    Vec3f transmittance = generalizedShadowRay(*event.sampler, ray, medium, nullptr, true, true, bounce);
     if (transmittance == 0.0f)
         return false;
 
@@ -270,7 +272,7 @@ Vec3f TraceBase::lightSample(const Primitive &light,
 
     IntersectionTemporary data;
     IntersectionInfo info;
-    Vec3f e = attenuatedEmission(*event.sampler, light, medium, sample.dist, data, info, bounce, ray, transmittance);
+    Vec3f e = attenuatedEmission(*event.sampler, light, medium, sample.dist, data, info, bounce, true, ray, transmittance);
     if (e == 0.0f)
         return Vec3f(0.0f);
 
@@ -306,7 +308,7 @@ Vec3f TraceBase::bsdfSample(const Primitive &light,
 
     IntersectionTemporary data;
     IntersectionInfo info;
-    Vec3f e = attenuatedEmission(*event.sampler, light, medium, -1.0f, data, info, bounce, ray, nullptr);
+    Vec3f e = attenuatedEmission(*event.sampler, light, medium, -1.0f, data, info, bounce, true, ray, nullptr);
 
     if (e == Vec3f(0.0f))
         return Vec3f(0.0f);
@@ -338,7 +340,7 @@ Vec3f TraceBase::volumeLightSample(PathSampleGenerator &sampler,
 
     IntersectionTemporary data;
     IntersectionInfo info;
-    Vec3f e = attenuatedEmission(sampler, light, medium, lightSample.dist, data, info, bounce, ray, nullptr);
+    Vec3f e = attenuatedEmission(sampler, light, medium, lightSample.dist, data, info, bounce, false, ray, nullptr);
     if (e == 0.0f)
         return Vec3f(0.0f);
 
@@ -366,7 +368,7 @@ Vec3f TraceBase::volumePhaseSample(const Primitive &light,
 
     IntersectionTemporary data;
     IntersectionInfo info;
-    Vec3f e = attenuatedEmission(sampler, light, medium, -1.0f, data, info, bounce, ray, nullptr);
+    Vec3f e = attenuatedEmission(sampler, light, medium, -1.0f, data, info, bounce, false, ray, nullptr);
 
     if (e == Vec3f(0.0f))
         return Vec3f(0.0f);

src/core/integrators/TraceBase.hpp

@@ -64,6 +64,7 @@ class TraceBase
                              IntersectionTemporary &data,
                              IntersectionInfo &info,
                              int bounce,
+                             bool startsOnSurface,
                              Ray &ray,
                              Vec3f *transmittance);
 
@@ -148,6 +149,8 @@ class TraceBase
                                Ray &ray,
                                const Medium *medium,
                                const Primitive *endCap,
+                               bool startsOnSurface,
+                               bool endsOnSurface,
                                int bounce) const;
     Vec3f generalizedShadowRayAndPdfs(PathSampleGenerator &sampler,
                                Ray &ray,

src/core/integrators/bidirectional_path_tracer/LightPath.cpp

@@ -96,19 +96,25 @@ void LightPath::toAreaMeasure()
 float LightPath::misWeight(const LightPath &camera, const LightPath &emitter,
             const PathEdge &edge, int s, int t, float *ratios)
 {
-    float *pdfForward  = reinterpret_cast<float *>(alloca((s + t)*sizeof(float)));
-    float *pdfBackward = reinterpret_cast<float *>(alloca((s + t)*sizeof(float)));
-    bool  *connectable = reinterpret_cast<bool  *>(alloca((s + t)*sizeof(bool)));
+    float *pdfForward           = reinterpret_cast<float *>(alloca((s + t)*sizeof(float)));
+    float *pdfBackward          = reinterpret_cast<float *>(alloca((s + t)*sizeof(float)));
+    bool  *connectable          = reinterpret_cast<bool  *>(alloca((s + t)*sizeof(bool)));
+    const PathVertex **vertices = reinterpret_cast<const PathVertex **>(alloca((s + t)*sizeof(const PathVertex *)));
+
+    if (!camera[t - 1].segmentConnectable(emitter[s - 1]))
+        return 0.0f;
 
     for (int i = 0; i < s; ++i) {
         pdfForward [i] = emitter[i].pdfForward();
         pdfBackward[i] = emitter[i].pdfBackward();
         connectable[i] = !emitter[i].isDirac();
+        vertices   [i] = &emitter[i];
     }
     for (int i = 0; i < t; ++i) {
         pdfForward [s + t - (i + 1)] = camera[i].pdfBackward();
         pdfBackward[s + t - (i + 1)] = camera[i].pdfForward();
         connectable[s + t - (i + 1)] = !camera[i].isDirac();
+        vertices   [s + t - (i + 1)] = &camera[i];
     }
     connectable[s - 1] = connectable[s] = true;
 
@@ -137,7 +143,7 @@ float LightPath::misWeight(const LightPath &camera, const LightPath &emitter,
         ratios[s] = 1.0f;
     for (int i = s + 1; i < s + t; ++i) {
         pi *= pdfForward[i - 1]/pdfBackward[i - 1];
-        if (connectable[i - 1] && connectable[i]) {
+        if (connectable[i - 1] && connectable[i] && vertices[i - 1]->segmentConnectable(*vertices[i])) {
             weight += pi;
             if (ratios)
                 ratios[i] = pi;
@@ -149,7 +155,7 @@ float LightPath::misWeight(const LightPath &camera, const LightPath &emitter,
     pi = 1.0f;
     for (int i = s - 1; i >= 1; --i) {
         pi *= pdfBackward[i]/pdfForward[i];
-        if (connectable[i - 1] && connectable[i]) {
+        if (connectable[i - 1] && connectable[i] && vertices[i - 1]->segmentConnectable(*vertices[i])) {
             weight += pi;
             if (ratios)
                 ratios[i] = pi;
@@ -166,7 +172,7 @@ float LightPath::misWeight(const LightPath &camera, const LightPath &emitter,
     } else {
         if (ratios)
             ratios[0] = 0.0f;
-     }
+    }
 
     return 1.0f/weight;
 }

src/core/integrators/bidirectional_path_tracer/PathVertex.cpp

@@ -160,7 +160,7 @@ bool PathVertex::sampleNextVertex(const TraceableScene &scene, TraceBase &tracer
         hitSurface = mediumRecord.mediumSample.exited;
         edgePdfForward = mediumRecord.mediumSample.pdf;
         Ray reverseRay(mediumRecord.mediumSample.p, -state.ray.dir(), 0.0f, mediumRecord.mediumSample.t);
-        edgePdfBackward = state.medium->pdf(state.sampler, reverseRay, onSurface());
+        edgePdfBackward = state.medium->pdf(state.sampler, reverseRay, hitSurface, onSurface());
         weight *= mediumRecord.mediumSample.weight;
         if (hitSurface && !didHit)
             return false;
@@ -327,6 +327,14 @@ void PathVertex::evalPdfs(const PathVertex *prev, const PathEdge *prevEdge, cons
     }}
 }
 
+bool PathVertex::segmentConnectable(const PathVertex &next) const
+{
+    if (onSurface() || next.onSurface())
+        return true;
+
+    return !_medium->isDirac();
+}
+
 void PathVertex::pointerFixup()
 {
     // Yuck. It's best not to ask.

src/core/integrators/bidirectional_path_tracer/PathVertex.hpp

@@ -132,6 +132,7 @@ class PathVertex
     Vec3f eval(const Vec3f &d, bool adjoint) const;
     void evalPdfs(const PathVertex *prev, const PathEdge *prevEdge, const PathVertex &next,
             const PathEdge &nextEdge, float *forward, float *backward) const;
+    bool segmentConnectable(const PathVertex &next) const;
 
     void pointerFixup();
 

src/core/integrators/light_tracer/LightTracer.cpp

@@ -29,7 +29,7 @@ void LightTracer::traceSample(PathSampleGenerator &sampler)
         Ray shadowRay(point.p, splat.d);
         shadowRay.setFarT(splat.dist);
 
-        Vec3f transmission = generalizedShadowRay(sampler, shadowRay, medium, nullptr, 0);
+        Vec3f transmission = generalizedShadowRay(sampler, shadowRay, medium, nullptr, true, true, 0);
         if (transmission != 0.0f) {
             Vec3f value = throughput*transmission*splat.weight
                     *light->evalDirectionalEmission(point, DirectionSample(splat.d));

src/core/integrators/path_tracer/PathTracer.cpp

@@ -50,6 +50,7 @@ Vec3f PathTracer::traceSample(Vec2u pixel, PathSampleGenerator &sampler)
     while ((didHit || medium) && bounce < _settings.maxBounces) {
         bool hitSurface = true;
         if (medium) {
+            mediumSample.continuedWeight = throughput;
             if (!medium->sampleDistance(sampler, ray, state, mediumSample))
                 return emission;
             throughput *= mediumSample.weight;