<p>For cases where the diffusion approximation is preferred over full volumetric path tracing, the method published by Solid Angle is fairly efficient:
<a href="https://www.solidangle.com/research/s2013_bssrdf_slides.pdf" rel="noreferrer">https://www.solidangle.com/research/s2013_bssrdf_slides.pdf</a></p>

<p>It is implemented in the Arnold render engine, Blender's Cycles, and pbrt, the latter being open source.</p>

<p>The file which implements it in PBRT is here:
<a href="https://github.com/mmp/pbrt-v3/blob/7095acfd8331cdefb03fe0bcae32c2fc9bd95980/src/core/bssrdf.cpp" rel="noreferrer">https://github.com/mmp/pbrt-v3/blob/7095acfd8331cdefb03fe0bcae32c2fc9bd95980/src/core/bssrdf.cpp</a></p>

CC BY-SA 3.0
2017-06-14T12:17:21.757Z
PBRT implementation seems like photon beam diffusion. Are photon beam diffusion and the what Solid Angle does the same? I ask because I thought PBD also take point samples on the geometry as a preprocessing step but Solid Angle slides don't do it.
Mustafa Işık
CC BY-SA 3.0
2017-06-21T08:36:52.163Z
There are two things to distinguish here: The function to integrate (BSSRDF) and the method of integration. The first one are functions like the dipole, multipole, photon beam diffusion, etc. The second are methods like Jensen's hierarchical octree or Solid Angle's importance sampled ray tracing.
Stefan Werner
CC BY-SA 3.0
2017-06-13T19:01:26.207Z
2017-06-13T19:01:26.207Z
Stefan Werner
<p>As mentioned in the comments, I would highly suggest starting with Full Volumetric Scattering. This is two fold:</p>

<ol>
<li>Since you are doing path tracing, adding volumetrics isn't super difficult.</li>
<li>Fully understanding how full volumetric scattering works will be a great basis for understanding the estimations. In addition, it can provide great "references" to see if your estimations are performing well / correctly.</li>
</ol>

<p>With that in mind, below is a basic introduction to how to implement Full Volumetric Scattering in a backwards path tracer.</p>

<hr>

<p>To start, let's look at the code for a backwards path tracer with only reflection, no transmission / refraction:</p>

<pre><code>void RenderPixel(uint x, uint y, UniformSampler *sampler) {
    Ray ray = m_scene-&gt;Camera-&gt;CalculateRayFromPixel(x, y, sampler);

    float3 color(0.0f);
    float3 throughput(1.0f);
    SurfaceInteraction interaction;

    // Bounce the ray around the scene
    const uint maxBounces = 15;
    for (uint bounces = 0; bounces &lt; maxBounces; ++bounces) {
        m_scene-&gt;Intersect(ray);

        // The ray missed. Return the background color
        if (ray.GeomID == INVALID_GEOMETRY_ID) {
            color += throughput * m_scene-&gt;BackgroundColor;
            break;
        }

        // Fetch the material
        Material *material = m_scene-&gt;GetMaterial(ray.GeomID);
        // The object might be emissive. If so, it will have a corresponding light
        // Otherwise, GetLight will return nullptr
        Light *light = m_scene-&gt;GetLight(ray.GeomID);

        // If this is the first bounce or if we just had a specular bounce,
        // we need to add the emmisive light
        if ((bounces == 0 || (interaction.SampledLobe &amp; BSDFLobe::Specular) != 0) &amp;&amp; light != nullptr) {
            color += throughput * light-&gt;Le();
        }

        interaction.Position = ray.Origin + ray.Direction * ray.TFar;
        interaction.Normal = normalize(m_scene-&gt;InterpolateNormal(ray.GeomID, ray.PrimID, ray.U, ray.V));
        interaction.OutputDirection = normalize(-ray.Direction);


        // Calculate the direct lighting
        color += throughput * SampleLights(sampler, interaction, material-&gt;bsdf, light);


        // Get the new ray direction
        // Choose the direction based on the bsdf        
        material-&gt;bsdf-&gt;Sample(interaction, sampler);
        float pdf = material-&gt;bsdf-&gt;Pdf(interaction);

        // Accumulate the weight
        throughput = throughput * material-&gt;bsdf-&gt;Eval(interaction) / pdf;

        // Shoot a new ray

        // Set the origin at the intersection point
        ray.Origin = interaction.Position;

        // Reset the other ray properties
        ray.Direction = interaction.InputDirection;
        ray.TNear = 0.001f;
        ray.TFar = infinity;


        // Russian Roulette
        if (bounces &gt; 3) {
            float p = std::max(throughput.x, std::max(throughput.y, throughput.z));
            if (sampler-&gt;NextFloat() &gt; p) {
                break;
            }

            throughput *= 1 / p;
        }
    }

    m_scene-&gt;Camera-&gt;FrameBufferData.SplatPixel(x, y, color);
}
</code></pre>

<h2>In English:</h2>

<ol>
<li>Shoot a ray through the scene.</li>
<li>Check if we hit anything. If not we return the skybox color and break.</li>
<li>If this is the first ray, or we just bounced off a specular surface, check if we hit a light. If so, we add the light emission to our color accumulation.</li>
<li>Sample the direct lighting.</li>
<li>Choose a new direction for the next ray. We can do this uniformly, or importance sample based on the BRDF.</li>
<li>Evaluate the BRDF and accumulate it. </li>
<li>Create a new ray based on our chosen direction and where we just came from.</li>
<li>[Optional] Use Russian Roulette to choose if we should terminate the ray.</li>
<li>Goto 1.</li>
</ol>

<p>For more details on Direct Light Sampling, see <a href="https://computergraphics.stackexchange.com/questions/5152/progressive-path-tracing-with-explicit-light-sampling/5153#5153">this answer</a>.</p>

<h2>BSDFs</h2>

<p>To add transmission, we first need to upgrade some of our BRDF materials to full BSDFs, with a BTDF function.
<a href="https://i.sstatic.net/2oBo9.png" rel="noreferrer"><img src="https://i.sstatic.net/2oBo9.png" alt="BSDF Explanation"></a></p>

<p>For example, a BSDF for an ideal specular dielectric (ie. glass) is:</p>

<pre><code>class IdealSpecularDielectric : public BSDF {
public:
    IdealSpecularDielectric(float3 albedo, float ior)
        : BSDF(BSDFLobe::Specular, albedo), 
          m_ior(ior) {
    }

private:
    float m_ior;

public:
    float3 Eval(SurfaceInteraction &amp;interaction) const override {
        return m_albedo;
    }

    void Sample(SurfaceInteraction &amp;interaction, UniformSampler *sampler) const override {
        float VdotN = dot(interaction.OutputDirection, interaction.Normal);
        float IORo = m_ior;
        if (VdotN &lt; 0.0f) {
            IORo = 1.0f;
            interaction.Normal = -interaction.Normal;
            VdotN = -VdotN;
        }

        float eta = interaction.IORi / IORo;
        float sinSquaredThetaT = SinSquaredThetaT(VdotN, eta);
        float fresnel = Fresnel(interaction.IORi, IORo, VdotN, sinSquaredThetaT);

        float rand = sampler-&gt;NextFloat();
        if (rand &lt;= fresnel) {
            // Reflect
            interaction.InputDirection = reflect(interaction.OutputDirection, interaction.Normal);
            interaction.SampledLobe = BSDFLobe::SpecularReflection;
            interaction.IORo = interaction.IORi;
        } else {
            // Refract
            interaction.InputDirection = refract(interaction.OutputDirection, interaction.Normal, VdotN, eta, sinSquaredThetaT);
            interaction.SampledLobe = BSDFLobe::SpecularTransmission;
            interaction.IORo = IORo;
        }

        if (AnyNan(interaction.InputDirection)) {
            printf("nan");
        }
    }

    float Pdf(SurfaceInteraction &amp;interaction) const override {
        return 1.0f;
    }
};
</code></pre>

<p>The interesting part of the code is <code>Sample()</code>. This is where a ray will choose whether to reflect or refract. How we do this is really up to us, as shown in <a href="https://computergraphics.stackexchange.com/questions/2482/choosing-reflection-or-refraction-in-path-tracing">this answer</a>. However, an obvious choice would be to sample based on the <a href="https://en.wikipedia.org/wiki/Fresnel_equations" rel="noreferrer">Fresnel equations</a>.</p>

<p>In the case of an ideal Specular Dielectric, there are only two potential output directions: perfect refraction or perfect reflection. So we evaluate the Fresnel and randomly choose to refract / reflect with proportion equal to the Fresnel. </p>

<h2>Media</h2>

<p>Next, we need to talk about media, and how they affect light bouncing around. As explained by Nathan Reed in <a href="https://computergraphics.stackexchange.com/a/3570/310">this answer</a>: </p>

<blockquote>
  <p>The way I like to think about volume scattering is that a photon traveling through a medium has a certain probability per unit length of interacting (getting scattered or absorbed). As long as it doesn't interact, it just goes in a straight line unimpeded and without losing energy. The greater the distance, the greater the probability that it interacts somewhere in that distance. The interaction probability per unit length is the coefficient $\sigma$ that you see in the equations. We usually have separate coefficients for scattering and absorption probabilities, so $\sigma = \sigma_s + \sigma_a$</p>
  
  <p>This probability per unit length is exactly the origin of the Beer-Lambert law. Slice a ray segment into infinitesimal intervals, treat each interval as an independent possible place to interact, then integrate along the ray; you get an exponential distribution (with rate parameter σσ) for the probability of interaction as a function of distance.</p>
  
  <p>You <em>can</em> technically choose the distance between events however you want, as long as you correctly weight the path for the probability that a photon can make it between two adjacent events without interacting with the medium. In other words, each path segment within the medium contributes a weight factor of $e^{-\sigma x}$, where $x$ is the length of the segment.</p>
  
  <p>Given this, a usually good choice for the distance is to importance-sample it from the exponential distribution. In other words, you set $x = -(\ln \xi)/\sigma$</p>
</blockquote>

<p>In summary:</p>

<ol>
<li>A ray travelling through a medium has some probability to:

<ul>
<li>Be absorbed</li>
<li>Be scattered</li>
</ul></li>
<li>Since this is Monte Carlo Integration, we are free to choose the interaction distance however we want. However, a good choice is to importance sample from an exponential distribution similar to Beer-Lambert.</li>
</ol>

<p>I chose to implement this with a scattering coefficient, and letting Beer-Lambert take care of the absorption coefficient.</p>

<pre><code>class Medium {
public:
    Medium(float3 absorptionColor, float absorptionAtDistance) 
        : m_absorptionCoefficient(-log(absorptionColor) / absorptionAtDistance) {
        // This method for calculating the absorption coefficient is borrowed from Burley's 2015 Siggraph Course Notes "Extending the Disney BRDF to a BSDF with Integrated Subsurface Scattering"
        // It's much more intutive to specify a color and a distance, then back-calculate the coefficient
    }
    virtual ~Medium() {
    }

protected:
    const float3a m_absorptionCoefficient;

public:
    virtual float SampleDistance(UniformSampler *sampler, float tFar, float *weight, float *pdf) const = 0;

    virtual float3a SampleScatterDirection(UniformSampler *sampler, float3a &amp;wo, float *pdf) const = 0;
    virtual float ScatterDirectionPdf(float3a &amp;wi, float3a &amp;wo) const = 0;

    virtual float3 Transmission(float distance) const = 0;
};
</code></pre>

<p>A Non-scattering medium is one where a photon travels straight through, but it is attenuated along the way according to Beer-Lambert:</p>

<pre><code>class NonScatteringMedium : public Medium {
public:
    NonScatteringMedium(float3 color, float atDistance)
        : Medium(color, atDistance) {
    }

public:
    float SampleDistance(UniformSampler *sampler, float tFar, float *weight, float *pdf) const override {
        *pdf = 1.0f;
        return tFar;
    }

    float3a SampleScatterDirection(UniformSampler *sampler, float3a &amp;wo, float *pdf) const override {
        return wo;
    }
    float ScatterDirectionPdf(float3a &amp;wi, float3a &amp;wo) const override {
        return 1.0f;
    }

    float3 Transmission(float distance) const override {
        return exp(-m_absorptionCoefficient * distance);
    }
};
</code></pre>

<p><a href="https://i.sstatic.net/aERCt.png" rel="noreferrer"><img src="https://i.sstatic.net/aERCt.png" alt="Non-scattering example"></a></p>

<p>An Isotropic Scattering Medium allows the ray to have scattering (reflection) events while the ray is traveling through it. The ray can be reflected in any direction in the sphere around the interaction point (hence Isotropic). </p>

<p><a href="https://i.sstatic.net/tOp55.jpg" rel="noreferrer"><img src="https://i.sstatic.net/tOp55.jpg" alt="Explanation of Isotropic scattering"></a></p>

<p>The distance between reflections is random, based on an exponential "scatter" probability.</p>

<pre><code>class IsotropicScatteringMedium : public Medium {
public:
    IsotropicScatteringMedium(float3 absorptionColor, float absorptionAtDistance, float scatteringDistance)
        : Medium(absorptionColor, absorptionAtDistance),
          m_scatteringCoefficient(1 / scatteringDistance) {
    }

private:
    float m_scatteringCoefficient;

public:
    float SampleDistance(UniformSampler *sampler, float tFar, float *weight, float *pdf) const override {
        float distance = -logf(sampler-&gt;NextFloat()) / m_scatteringCoefficient;
        // If we sample a distance farther than the next intersecting surface, clamp to the surface distance
        if (distance &gt;= tFar) {
            *pdf = 1.0f;
            return tFar;
        }

        *pdf = std::exp(-m_scatteringCoefficient * distance);
        return distance;
    }
    float3a SampleScatterDirection(UniformSampler *sampler, float3a &amp;wo, float *pdf) const override {
        *pdf = 0.25f * M_1_PI; // 1 / (4 * PI)
        return UniformSampleSphere(sampler);
    }

    float ScatterDirectionPdf(float3a &amp;wi, float3a &amp;wo) const override {
        return 0.25f * M_1_PI; // 1 / (4 * PI)
    }
    float3 Transmission(float distance) const override {
        return exp(-m_absorptionCoefficient * distance);
    }
};
</code></pre>

<p>At high scattering distances, the material behaves almost like a NonScatteringMedium, because it has a very small probability of scattering before it travels through the medium.</p>

<p><a href="https://i.sstatic.net/H7OX2.png" rel="noreferrer"><img src="https://i.sstatic.net/H7OX2.png" alt="High Scattering Distance"></a></p>

<p>At low scattering distances, the material behaves like wax or jade. </p>

<p><a href="https://i.sstatic.net/JgjlE.png" rel="noreferrer"><img src="https://i.sstatic.net/JgjlE.png" alt="Low Scattering Distance"></a></p>

<p>As the scattering distance gets lower and lower, we have a higher and higher chance to scatter while traveling through the medium. This presents a few issues:</p>

<ol>
<li>The number of bounces per ray is going to explode. </li>
<li>If Russian Roulette kills the ray, or if you go over "maxBounces", you will "lose" any color contribution from inside the medium. <a href="https://i.sstatic.net/tZ6QO.png" rel="noreferrer"><img src="https://i.sstatic.net/tZ6QO.png" alt="Over Max bounces example"></a></li>
</ol>

<p>Thus, you need to set maxBounces to some high number, and rely on Russian Roulette to terminate rays in the general case. In addition, IsotropicScatteringMedium is extremely inefficient for representing mostly opaque materials. To get better performance, we should use a non Isotropic phase function, <a href="https://www.cs.dartmouth.edu/~wjarosz/publications/dissertation/chapter4.pdf" rel="noreferrer">like Henyey-Greenstein or Schlick</a>. These bias the scattering in a certain direction. So we could set them to have high back-scattering, thereby lessening the issue of "lost" rays.</p>

<h2>Putting it all together</h2>

<p>So, with this new information, how do we modify the integrator to understand BSDFs and media? </p>

<p>First, we need to start tracking what medium we are currently in, and its Index of Refraction (IOR). </p>

<pre><code>SurfaceInteraction interaction;
interaction.IORi = 1.0f; // Vacuum 
Medium *medium = nullptr; // nullptr == vacuum
</code></pre>

<p>Then, we break the integrator up into two parts: transmission and material interaction.</p>

<pre><code>void RenderPixel(uint x, uint y, UniformSampler *sampler) const {
    Ray ray = m_scene-&gt;Camera-&gt;CalculateRayFromPixel(x, y, sampler);

    float3 color(0.0f);
    float3 throughput(1.0f);
    SurfaceInteraction interaction;
    interaction.IORi = 1.0f; // Air
    Medium *medium = nullptr;
    bool hitSurface = false;

    // Bounce the ray around the scene
    uint bounces = 0;
    const uint maxBounces = 1500;
    for (; bounces &lt; maxBounces; ++bounces) {
        m_scene-&gt;Intersect(ray);

        // The ray missed. Return the background color
        if (ray.GeomID == INVALID_GEOMETRY_ID) {
            color += throughput * m_scene-&gt;BackgroundColor;
            break;
        }

        // We hit an object
        hitSurface = true;

        // Calculate any transmission
        if (medium != nullptr) {
            float weight = 1.0f;
            float pdf = 1.0f;
            float distance = medium-&gt;SampleDistance(sampler, ray.TFar, &amp;weight, &amp;pdf);
            float3 transmission = medium-&gt;Transmission(distance);
            throughput = throughput * weight * transmission;

            if (distance &lt; ray.TFar) {
                // Create a scatter event
                hitSurface = false;

                ray.Origin = ray.Origin + ray.Direction * distance;

                // Reset the other ray properties
                float directionPdf;
                float3a wo = normalize(ray.Direction);
                ray.Direction = medium-&gt;SampleScatterDirection(sampler, wo, &amp;directionPdf);
                ray.TNear = 0.001f;
                ray.TFar = infinity;
                ray.GeomID = INVALID_GEOMETRY_ID;
                ray.PrimID = INVALID_PRIMATIVE_ID;
                ray.InstID = INVALID_INSTANCE_ID;
                ray.Mask = 0xFFFFFFFF;
                ray.Time = 0.0f;
            }
        }

        if (hitSurface) {
            // Fetch the material
            Material *material = m_scene-&gt;GetMaterial(ray.GeomID);
            // The object might be emissive. If so, it will have a corresponding light
            // Otherwise, GetLight will return nullptr
            Light *light = m_scene-&gt;GetLight(ray.GeomID);

            // If this is the first bounce or if we just had a specular bounce,
            // we need to add the emmisive light
            if ((bounces == 0 || (interaction.SampledLobe &amp; BSDFLobe::Specular) != 0) &amp;&amp; light != nullptr) {
                color += throughput * light-&gt;Le();
            }

            interaction.Position = ray.Origin + ray.Direction * ray.TFar;
            interaction.Normal = normalize(m_scene-&gt;InterpolateNormal(ray.GeomID, ray.PrimID, ray.U, ray.V));
            interaction.OutputDirection = normalize(-ray.Direction);
            interaction.IORo = 0.0f;


            // Calculate the direct lighting
            color += throughput * SampleOneLight(sampler, interaction, material-&gt;bsdf, light);


            // Get the new ray direction
            // Choose the direction based on the bsdf        
            material-&gt;bsdf-&gt;Sample(interaction, sampler);
            float pdf = material-&gt;bsdf-&gt;Pdf(interaction);

            // Accumulate the weight
            throughput = throughput * material-&gt;bsdf-&gt;Eval(interaction) / pdf;

            // Update the current IOR and medium if we refracted
            if (interaction.SampledLobe == BSDFLobe::SpecularTransmission) {
                interaction.IORi = interaction.IORo;
                medium = material-&gt;medium;
            }

            // Shoot a new ray

            // Set the origin at the intersection point
            ray.Origin = interaction.Position;

            // Reset the other ray properties
            ray.Direction = interaction.InputDirection;
            if (AnyNan(ray.Direction))
                printf("bad");
            ray.TNear = 0.001f;
            ray.TFar = infinity;
            ray.GeomID = INVALID_GEOMETRY_ID;
            ray.PrimID = INVALID_PRIMATIVE_ID;
            ray.InstID = INVALID_INSTANCE_ID;
            ray.Mask = 0xFFFFFFFF;
            ray.Time = 0.0f;
        }

        // Russian Roulette
        if (bounces &gt; 3) {
            float p = std::max(throughput.x, std::max(throughput.y, throughput.z));
            if (sampler-&gt;NextFloat() &gt; p) {
                break;
            }

            throughput *= 1 / p;
        }
    }

    if (bounces == maxBounces) {
        printf("Over max bounces");
    }

    m_scene-&gt;Camera-&gt;FrameBufferData.SplatPixel(x, y, color);
}
</code></pre>

<h2>In English:</h2>

<ol>
<li>Shoot a ray through the scene.</li>
<li>Check if we hit anything. If not we return the skybox color and break.</li>
<li>Check if we are currently in a medium (not in a vacuum).

<ol>
<li>Sample the medium for a distance.</li>
<li>Evaluate the medium's transmission and accumulate the throughput.</li>
<li>If the scatter distance is less than the distance from the ray origin to the next surface, make a scatter event.</li>
<li>Sample the medium for a scatter direction.</li>
</ol></li>
<li>Check if we hit a surface (ie, we didn't have a scatter event).

<ol>
<li>If this is the first ray, or we just bounced off a specular surface, check if we hit a light. If so, we add the light emission to our color accumulation.</li>
<li>Sample the direct lighting.</li>
<li>Choose a new direction for the next ray. We can do this uniformly, or importance sample based on the BRDF.</li>
<li>Evaluate the BRDF and accumulate it. </li>
<li>If we refracted into the material, update the IOR and medium for the next bounce.</li>
</ol></li>
<li>Create a new ray based on our chosen direction and where we just came from.</li>
<li>[Optional] Use Russian Roulette to choose if we should terminate the ray.</li>
<li>Goto 1.</li>
</ol>

CC BY-SA 3.0
2017-06-12T07:04:07.920Z
Thank you for this amazing answer! The problem is that no matter what I try for the coefficients, my outputs look like a glass. Does this mean I am doing scattering part wrong? Here is the sequence of outputs: http://imgur.com/a/gg2u1
I only changed the scattering coefficient while keeping absorption the same
Mustafa Işık
CC BY-SA 3.0
2017-06-12T13:25:28.770Z
Are you sure that you're not accidentally setting the coefficient equal to the scattering distance? It should be scatteringCoefficient = (1 / scatteringDistance). Also absorptionCoefficient = -log(absorptionColor) / absorptionAtDistance. In the end, it doesn't really matter how you come up with the coefficient. That said, a lower coefficient will mean less absorption / scattering. A higher coefficient will mean less. For artists, it's a bit easier to think of it in terms of distances, so we back-calculate the coefficients from distances.
RichieSams
CC BY-SA 3.0
2017-06-15T07:19:48.460Z
No, parameter setting part is the same. I couldn't solve the problem but if I solve it, I will comment here.
Mustafa Işık
CC BY-SA 3.0
2017-06-11T18:53:07.640Z
2017-06-12T17:56:21.303Z
2017-06-12T17:56:21.303Z
RichieSams
<p>I read "A Practical Model for Subsurface Light Transport" and "A Rapid Hierarchical Rendering Technique for Translucent Materials". If I understand correctly, the former does not consider the case when an object is lit behind and we can see the light is going through it.</p>

<p>The latter can achieve this but there is a preprocessing stage which is not suitable for my CUDA path tracer and I want it to be as natural as possible. So, do you have any paper in mind which is not obsolete or would you suggest me to solve the full radiative transfer equation?</p>

CC BY-SA 3.0
2017-06-08T17:17:30.390Z
Are you looking to estimate subsurface scattering in your path tracer, or evaluate it exactly using full-volumetric scattering?
RichieSams
CC BY-SA 3.0
2017-06-08T19:22:55.860Z
I read that simulating full-volumetric scattering is really slow. If it is not that slow I would like to implement it. Although diffusion approximations work faster, they assume lots of things or they need ugly preprocessing steps. Have you done it before? Any suggestion would be great!
Mustafa Işık
CC BY-SA 3.0
2017-06-08T19:57:59.250Z
It can be slow, depending on the implementation. I have implemented brute-force, isotropic scattering media in my hobby path tracer. Example image: http://imgur.com/a/oy9F8. (The FPS / total frame times are in the upper left corner). Using a backscattering phase function like Henyey-Greenstein might make it converge faster for more opaque materials, but I'm not sure since I haven't tried it. If nothing else, implementing full volumetric scattering is a great way to learn what the approximations are trying to approximate.
RichieSams
CC BY-SA 3.0
2017-06-08T19:59:38.570Z
I'll try to do a write up later tonight when I get home. In the mean time, here's the link to my path tracer. https://github.com/RichieSams/lantern It's not production ready at all. Rather, it's just for helping me learn the math behind everything. That said, it's open source under Apache 2.0, so feel free to copy and learn from it yourself.
RichieSams
CC BY-SA 3.0
2017-06-08T20:01:22.773Z
Lastly, here are two questions that I asked here about this topic: https://computergraphics.stackexchange.com/questions/2482/choosing-reflection-or-refraction-in-path-tracing https://computergraphics.stackexchange.com/questions/2563/full-monte-carlo-volumetric-scattering
RichieSams
CC BY-SA 3.0
2017-06-09T08:08:28.930Z
Thank you! I implemented perfect refractive and reflective surfaces already. By the way, the answers to both questions are enlightening for me but I would prefer to see something you wrote up.
Mustafa Işık
CC BY-SA 3.0
2017-06-11T13:40:01.567Z
@RichieSams I read physically based rendering and got the intuition. There is only one part that confuses me. We cannot sample the light sources directly due to refraction at the boundaries, right? If so, do we have to wait until the ray we are tracing exit the object because otherwise there would be lots of black pixels on the image?
Mustafa Işık
CC BY-SA 3.0
2017-06-11T15:24:57.320Z
Sorry, I got busy this week. I'm writing up the answer now. But yes, you are correct. Light sampling is almost impossible from within media with traditional path tracing. Photon mapping works quite well in this area, though. Vertex Connection and Merging combines the best of both worlds: http://www.smallvcm.com/
RichieSams
<raytracing><pathtracing><global-illumination><subsurface-scattering><cuda>
A recent approach for subsurface scattering
Mustafa Işık
# A recent approach for subsurface scattering

(Asked by: username_0 on 2017-06-08)

I read "A Practical Model for Subsurface Light Transport" and "A Rapid Hierarchical Rendering Technique for Translucent Materials". If I understand correctly, the former does not consider the case when an object is lit behind and we can see the light is going through it.

The latter can achieve this but there is a preprocessing stage which is not suitable for my CUDA path tracer and I want it to be as natural as possible. So, do you have any paper in mind which is not obsolete or would you suggest me to solve the full radiative transfer equation?

**Tags:** raytracing, pathtracing, global-illumination, subsurface-scattering, cuda

> **Comments:**
> **username_1**: Are you looking to estimate subsurface scattering in your path tracer, or evaluate it exactly using full-volumetric scattering?  
> **username_0**: I read that simulating full-volumetric scattering is really slow. If it is not that slow I would like to implement it. Although diffusion approximations work faster, they assume lots of things or they need ugly preprocessing steps. Have you done it before? Any suggestion would be great!  
> **username_1**: It can be slow, depending on the implementation. I have implemented brute-force, isotropic scattering media in my hobby path tracer. Example image: http://imgur.com/a/oy9F8. (The FPS / total frame times are in the upper left corner). Using a backscattering phase function like Henyey-Greenstein might make it converge faster for more opaque materials, but I'm not sure since I haven't tried it. If nothing else, implementing full volumetric scattering is a great way to learn what the approximations are trying to approximate.  
> **username_1**: I'll try to do a write up later tonight when I get home. In the mean time, here's the link to my path tracer. https://github.com/username_1/lantern It's not production ready at all. Rather, it's just for helping me learn the math behind everything. That said, it's open source under Apache 2.0, so feel free to copy and learn from it yourself.  
> **username_1**: Lastly, here are two questions that I asked here about this topic: https://computergraphics.stackexchange.com/questions/2482/choosing-reflection-or-refraction-in-path-tracing https://computergraphics.stackexchange.com/questions/2563/full-monte-carlo-volumetric-scattering  
> **username_0**: Thank you! I implemented perfect refractive and reflective surfaces already. By the way, the answers to both questions are enlightening for me but I would prefer to see something you wrote up.  
> **username_0**: @username_1 I read physically based rendering and got the intuition. There is only one part that confuses me. We cannot sample the light sources directly due to refraction at the boundaries, right? If so, do we have to wait until the ray we are tracing exit the object because otherwise there would be lots of black pixels on the image?  
> **username_1**: Sorry, I got busy this week. I'm writing up the answer now. But yes, you are correct. Light sampling is almost impossible from within media with traditional path tracing. Photon mapping works quite well in this area, though. Vertex Connection and Merging combines the best of both worlds: http://www.smallvcm.com/  

---

## 2 answers

### Answer by username_1 on 2017-06-11. Score: 13 | Accepted answer

As mentioned in the comments, I would highly suggest starting with Full Volumetric Scattering. This is two fold:

1. Since you are doing path tracing, adding volumetrics isn't super difficult.
2. Fully understanding how full volumetric scattering works will be a great basis for understanding the estimations. In addition, it can provide great "references" to see if your estimations are performing well / correctly.

With that in mind, below is a basic introduction to how to implement Full Volumetric Scattering in a backwards path tracer.

---

To start, let's look at the code for a backwards path tracer with only reflection, no transmission / refraction:

```
void RenderPixel(uint x, uint y, UniformSampler *sampler) {
    Ray ray = m_scene->Camera->CalculateRayFromPixel(x, y, sampler);

    float3 color(0.0f);
    float3 throughput(1.0f);
    SurfaceInteraction interaction;

    // Bounce the ray around the scene
    const uint maxBounces = 15;
    for (uint bounces = 0; bounces < maxBounces; ++bounces) {
        m_scene->Intersect(ray);

        // The ray missed. Return the background color
        if (ray.GeomID == INVALID_GEOMETRY_ID) {
            color += throughput * m_scene->BackgroundColor;
            break;
        }

        // Fetch the material
        Material *material = m_scene->GetMaterial(ray.GeomID);
        // The object might be emissive. If so, it will have a corresponding light
        // Otherwise, GetLight will return nullptr
        Light *light = m_scene->GetLight(ray.GeomID);

        // If this is the first bounce or if we just had a specular bounce,
        // we need to add the emmisive light
        if ((bounces == 0 || (interaction.SampledLobe & BSDFLobe::Specular) != 0) && light != nullptr) {
            color += throughput * light->Le();
        }

        interaction.Position = ray.Origin + ray.Direction * ray.TFar;
        interaction.Normal = normalize(m_scene->InterpolateNormal(ray.GeomID, ray.PrimID, ray.U, ray.V));
        interaction.OutputDirection = normalize(-ray.Direction);

        // Calculate the direct lighting
        color += throughput * SampleLights(sampler, interaction, material->bsdf, light);

        // Get the new ray direction
        // Choose the direction based on the bsdf        
        material->bsdf->Sample(interaction, sampler);
        float pdf = material->bsdf->Pdf(interaction);

        // Accumulate the weight
        throughput = throughput * material->bsdf->Eval(interaction) / pdf;

        // Shoot a new ray

        // Set the origin at the intersection point
        ray.Origin = interaction.Position;

        // Reset the other ray properties
        ray.Direction = interaction.InputDirection;
        ray.TNear = 0.001f;
        ray.TFar = infinity;

        // Russian Roulette
        if (bounces > 3) {
            float p = std::max(throughput.x, std::max(throughput.y, throughput.z));
            if (sampler->NextFloat() > p) {
                break;
            }

            throughput *= 1 / p;
        }
    }

    m_scene->Camera->FrameBufferData.SplatPixel(x, y, color);
}

```

## In English:

1. Shoot a ray through the scene.
2. Check if we hit anything. If not we return the skybox color and break.
3. If this is the first ray, or we just bounced off a specular surface, check if we hit a light. If so, we add the light emission to our color accumulation.
4. Sample the direct lighting.
5. Choose a new direction for the next ray. We can do this uniformly, or importance sample based on the BRDF.
6. Evaluate the BRDF and accumulate it.
7. Create a new ray based on our chosen direction and where we just came from.
8. [Optional] Use Russian Roulette to choose if we should terminate the ray.
9. Goto 1.

For more details on Direct Light Sampling, see [this answer](https://computergraphics.stackexchange.com/questions/5152/progressive-path-tracing-with-explicit-light-sampling/5153#5153).

## BSDFs

To add transmission, we first need to upgrade some of our BRDF materials to full BSDFs, with a BTDF function.
[![BSDF Explanation](https://i.sstatic.net/2oBo9.png)](https://i.sstatic.net/2oBo9.png)

For example, a BSDF for an ideal specular dielectric (ie. glass) is:

```
class IdealSpecularDielectric : public BSDF {
public:
    IdealSpecularDielectric(float3 albedo, float ior)
        : BSDF(BSDFLobe::Specular, albedo), 
          m_ior(ior) {
    }

private:
    float m_ior;

public:
    float3 Eval(SurfaceInteraction &interaction) const override {
        return m_albedo;
    }

    void Sample(SurfaceInteraction &interaction, UniformSampler *sampler) const override {
        float VdotN = dot(interaction.OutputDirection, interaction.Normal);
        float IORo = m_ior;
        if (VdotN < 0.0f) {
            IORo = 1.0f;
            interaction.Normal = -interaction.Normal;
            VdotN = -VdotN;
        }

        float eta = interaction.IORi / IORo;
        float sinSquaredThetaT = SinSquaredThetaT(VdotN, eta);
        float fresnel = Fresnel(interaction.IORi, IORo, VdotN, sinSquaredThetaT);

        float rand = sampler->NextFloat();
        if (rand <= fresnel) {
            // Reflect
            interaction.InputDirection = reflect(interaction.OutputDirection, interaction.Normal);
            interaction.SampledLobe = BSDFLobe::SpecularReflection;
            interaction.IORo = interaction.IORi;
        } else {
            // Refract
            interaction.InputDirection = refract(interaction.OutputDirection, interaction.Normal, VdotN, eta, sinSquaredThetaT);
            interaction.SampledLobe = BSDFLobe::SpecularTransmission;
            interaction.IORo = IORo;
        }

        if (AnyNan(interaction.InputDirection)) {
            printf("nan");
        }
    }

    float Pdf(SurfaceInteraction &interaction) const override {
        return 1.0f;
    }
};

```

The interesting part of the code is `Sample()`. This is where a ray will choose whether to reflect or refract. How we do this is really up to us, as shown in [this answer](https://computergraphics.stackexchange.com/questions/2482/choosing-reflection-or-refraction-in-path-tracing). However, an obvious choice would be to sample based on the [Fresnel equations](https://en.wikipedia.org/wiki/Fresnel_equations).

In the case of an ideal Specular Dielectric, there are only two potential output directions: perfect refraction or perfect reflection. So we evaluate the Fresnel and randomly choose to refract / reflect with proportion equal to the Fresnel.

## Media

Next, we need to talk about media, and how they affect light bouncing around. As explained by Nathan Reed in [this answer](https://computergraphics.stackexchange.com/a/3570/310):

> The way I like to think about volume scattering is that a photon traveling through a medium has a certain probability per unit length of interacting (getting scattered or absorbed). As long as it doesn't interact, it just goes in a straight line unimpeded and without losing energy. The greater the distance, the greater the probability that it interacts somewhere in that distance. The interaction probability per unit length is the coefficient $\sigma$ that you see in the equations. We usually have separate coefficients for scattering and absorption probabilities, so $\sigma = \sigma\_s + \sigma\_a$
>
> This probability per unit length is exactly the origin of the Beer-Lambert law. Slice a ray segment into infinitesimal intervals, treat each interval as an independent possible place to interact, then integrate along the ray; you get an exponential distribution (with rate parameter σσ) for the probability of interaction as a function of distance.
>
> You *can* technically choose the distance between events however you want, as long as you correctly weight the path for the probability that a photon can make it between two adjacent events without interacting with the medium. In other words, each path segment within the medium contributes a weight factor of $e^{-\sigma x}$, where $x$ is the length of the segment.
>
> Given this, a usually good choice for the distance is to importance-sample it from the exponential distribution. In other words, you set $x = -(\ln \xi)/\sigma$

In summary:

1. A ray travelling through a medium has some probability to:
   * Be absorbed
   * Be scattered
2. Since this is Monte Carlo Integration, we are free to choose the interaction distance however we want. However, a good choice is to importance sample from an exponential distribution similar to Beer-Lambert.

I chose to implement this with a scattering coefficient, and letting Beer-Lambert take care of the absorption coefficient.

```
class Medium {
public:
    Medium(float3 absorptionColor, float absorptionAtDistance) 
        : m_absorptionCoefficient(-log(absorptionColor) / absorptionAtDistance) {
        // This method for calculating the absorption coefficient is borrowed from Burley's 2015 Siggraph Course Notes "Extending the Disney BRDF to a BSDF with Integrated Subsurface Scattering"
        // It's much more intutive to specify a color and a distance, then back-calculate the coefficient
    }
    virtual ~Medium() {
    }

protected:
    const float3a m_absorptionCoefficient;

public:
    virtual float SampleDistance(UniformSampler *sampler, float tFar, float *weight, float *pdf) const = 0;

    virtual float3a SampleScatterDirection(UniformSampler *sampler, float3a &wo, float *pdf) const = 0;
    virtual float ScatterDirectionPdf(float3a &wi, float3a &wo) const = 0;

    virtual float3 Transmission(float distance) const = 0;
};

```

A Non-scattering medium is one where a photon travels straight through, but it is attenuated along the way according to Beer-Lambert:

```
class NonScatteringMedium : public Medium {
public:
    NonScatteringMedium(float3 color, float atDistance)
        : Medium(color, atDistance) {
    }

public:
    float SampleDistance(UniformSampler *sampler, float tFar, float *weight, float *pdf) const override {
        *pdf = 1.0f;
        return tFar;
    }

    float3a SampleScatterDirection(UniformSampler *sampler, float3a &wo, float *pdf) const override {
        return wo;
    }
    float ScatterDirectionPdf(float3a &wi, float3a &wo) const override {
        return 1.0f;
    }

    float3 Transmission(float distance) const override {
        return exp(-m_absorptionCoefficient * distance);
    }
};

```

[![Non-scattering example](https://i.sstatic.net/aERCt.png)](https://i.sstatic.net/aERCt.png)

An Isotropic Scattering Medium allows the ray to have scattering (reflection) events while the ray is traveling through it. The ray can be reflected in any direction in the sphere around the interaction point (hence Isotropic).

[![Explanation of Isotropic scattering](https://i.sstatic.net/tOp55.jpg)](https://i.sstatic.net/tOp55.jpg)

The distance between reflections is random, based on an exponential "scatter" probability.

```
class IsotropicScatteringMedium : public Medium {
public:
    IsotropicScatteringMedium(float3 absorptionColor, float absorptionAtDistance, float scatteringDistance)
        : Medium(absorptionColor, absorptionAtDistance),
          m_scatteringCoefficient(1 / scatteringDistance) {
    }

private:
    float m_scatteringCoefficient;

public:
    float SampleDistance(UniformSampler *sampler, float tFar, float *weight, float *pdf) const override {
        float distance = -logf(sampler->NextFloat()) / m_scatteringCoefficient;
        // If we sample a distance farther than the next intersecting surface, clamp to the surface distance
        if (distance >= tFar) {
            *pdf = 1.0f;
            return tFar;
        }

        *pdf = std::exp(-m_scatteringCoefficient * distance);
        return distance;
    }
    float3a SampleScatterDirection(UniformSampler *sampler, float3a &wo, float *pdf) const override {
        *pdf = 0.25f * M_1_PI; // 1 / (4 * PI)
        return UniformSampleSphere(sampler);
    }

    float ScatterDirectionPdf(float3a &wi, float3a &wo) const override {
        return 0.25f * M_1_PI; // 1 / (4 * PI)
    }
    float3 Transmission(float distance) const override {
        return exp(-m_absorptionCoefficient * distance);
    }
};

```

At high scattering distances, the material behaves almost like a NonScatteringMedium, because it has a very small probability of scattering before it travels through the medium.

[![High Scattering Distance](https://i.sstatic.net/H7OX2.png)](https://i.sstatic.net/H7OX2.png)

At low scattering distances, the material behaves like wax or jade.

[![Low Scattering Distance](https://i.sstatic.net/JgjlE.png)](https://i.sstatic.net/JgjlE.png)

As the scattering distance gets lower and lower, we have a higher and higher chance to scatter while traveling through the medium. This presents a few issues:

1. The number of bounces per ray is going to explode.
2. If Russian Roulette kills the ray, or if you go over "maxBounces", you will "lose" any color contribution from inside the medium. [![Over Max bounces example](https://i.sstatic.net/tZ6QO.png)](https://i.sstatic.net/tZ6QO.png)

Thus, you need to set maxBounces to some high number, and rely on Russian Roulette to terminate rays in the general case. In addition, IsotropicScatteringMedium is extremely inefficient for representing mostly opaque materials. To get better performance, we should use a non Isotropic phase function, [like Henyey-Greenstein or Schlick](https://www.cs.dartmouth.edu/~wjarosz/publications/dissertation/chapter4.pdf). These bias the scattering in a certain direction. So we could set them to have high back-scattering, thereby lessening the issue of "lost" rays.

## Putting it all together

So, with this new information, how do we modify the integrator to understand BSDFs and media?

First, we need to start tracking what medium we are currently in, and its Index of Refraction (IOR).

```
SurfaceInteraction interaction;
interaction.IORi = 1.0f; // Vacuum 
Medium *medium = nullptr; // nullptr == vacuum

```

Then, we break the integrator up into two parts: transmission and material interaction.

```
void RenderPixel(uint x, uint y, UniformSampler *sampler) const {
    Ray ray = m_scene->Camera->CalculateRayFromPixel(x, y, sampler);

    float3 color(0.0f);
    float3 throughput(1.0f);
    SurfaceInteraction interaction;
    interaction.IORi = 1.0f; // Air
    Medium *medium = nullptr;
    bool hitSurface = false;

    // Bounce the ray around the scene
    uint bounces = 0;
    const uint maxBounces = 1500;
    for (; bounces < maxBounces; ++bounces) {
        m_scene->Intersect(ray);

        // The ray missed. Return the background color
        if (ray.GeomID == INVALID_GEOMETRY_ID) {
            color += throughput * m_scene->BackgroundColor;
            break;
        }

        // We hit an object
        hitSurface = true;

        // Calculate any transmission
        if (medium != nullptr) {
            float weight = 1.0f;
            float pdf = 1.0f;
            float distance = medium->SampleDistance(sampler, ray.TFar, &weight, &pdf);
            float3 transmission = medium->Transmission(distance);
            throughput = throughput * weight * transmission;

            if (distance < ray.TFar) {
                // Create a scatter event
                hitSurface = false;

                ray.Origin = ray.Origin + ray.Direction * distance;

                // Reset the other ray properties
                float directionPdf;
                float3a wo = normalize(ray.Direction);
                ray.Direction = medium->SampleScatterDirection(sampler, wo, &directionPdf);
                ray.TNear = 0.001f;
                ray.TFar = infinity;
                ray.GeomID = INVALID_GEOMETRY_ID;
                ray.PrimID = INVALID_PRIMATIVE_ID;
                ray.InstID = INVALID_INSTANCE_ID;
                ray.Mask = 0xFFFFFFFF;
                ray.Time = 0.0f;
            }
        }

        if (hitSurface) {
            // Fetch the material
            Material *material = m_scene->GetMaterial(ray.GeomID);
            // The object might be emissive. If so, it will have a corresponding light
            // Otherwise, GetLight will return nullptr
            Light *light = m_scene->GetLight(ray.GeomID);

            // If this is the first bounce or if we just had a specular bounce,
            // we need to add the emmisive light
            if ((bounces == 0 || (interaction.SampledLobe & BSDFLobe::Specular) != 0) && light != nullptr) {
                color += throughput * light->Le();
            }

            interaction.Position = ray.Origin + ray.Direction * ray.TFar;
            interaction.Normal = normalize(m_scene->InterpolateNormal(ray.GeomID, ray.PrimID, ray.U, ray.V));
            interaction.OutputDirection = normalize(-ray.Direction);
            interaction.IORo = 0.0f;

            // Calculate the direct lighting
            color += throughput * SampleOneLight(sampler, interaction, material->bsdf, light);

            // Get the new ray direction
            // Choose the direction based on the bsdf        
            material->bsdf->Sample(interaction, sampler);
            float pdf = material->bsdf->Pdf(interaction);

            // Accumulate the weight
            throughput = throughput * material->bsdf->Eval(interaction) / pdf;

            // Update the current IOR and medium if we refracted
            if (interaction.SampledLobe == BSDFLobe::SpecularTransmission) {
                interaction.IORi = interaction.IORo;
                medium = material->medium;
            }

            // Shoot a new ray

            // Set the origin at the intersection point
            ray.Origin = interaction.Position;

            // Reset the other ray properties
            ray.Direction = interaction.InputDirection;
            if (AnyNan(ray.Direction))
                printf("bad");
            ray.TNear = 0.001f;
            ray.TFar = infinity;
            ray.GeomID = INVALID_GEOMETRY_ID;
            ray.PrimID = INVALID_PRIMATIVE_ID;
            ray.InstID = INVALID_INSTANCE_ID;
            ray.Mask = 0xFFFFFFFF;
            ray.Time = 0.0f;
        }

        // Russian Roulette
        if (bounces > 3) {
            float p = std::max(throughput.x, std::max(throughput.y, throughput.z));
            if (sampler->NextFloat() > p) {
                break;
            }

            throughput *= 1 / p;
        }
    }

    if (bounces == maxBounces) {
        printf("Over max bounces");
    }

    m_scene->Camera->FrameBufferData.SplatPixel(x, y, color);
}

```

## In English:

1. Shoot a ray through the scene.
2. Check if we hit anything. If not we return the skybox color and break.
3. Check if we are currently in a medium (not in a vacuum).
   1. Sample the medium for a distance.
   2. Evaluate the medium's transmission and accumulate the throughput.
   3. If the scatter distance is less than the distance from the ray origin to the next surface, make a scatter event.
   4. Sample the medium for a scatter direction.
4. Check if we hit a surface (ie, we didn't have a scatter event).
   1. If this is the first ray, or we just bounced off a specular surface, check if we hit a light. If so, we add the light emission to our color accumulation.
   2. Sample the direct lighting.
   3. Choose a new direction for the next ray. We can do this uniformly, or importance sample based on the BRDF.
   4. Evaluate the BRDF and accumulate it.
   5. If we refracted into the material, update the IOR and medium for the next bounce.
5. Create a new ray based on our chosen direction and where we just came from.
6. [Optional] Use Russian Roulette to choose if we should terminate the ray.
7. Goto 1.

> **Comments:**
> **username_0**: Thank you for this amazing answer! The problem is that no matter what I try for the coefficients, my outputs look like a glass. Does this mean I am doing scattering part wrong? Here is the sequence of outputs: http://imgur.com/a/gg2u1 I only changed the scattering coefficient while keeping absorption the same  
> **username_1**: Are you sure that you're not accidentally setting the coefficient equal to the scattering distance? It should be scatteringCoefficient = (1 / scatteringDistance). Also absorptionCoefficient = -log(absorptionColor) / absorptionAtDistance. In the end, it doesn't really matter how you come up with the coefficient. That said, a lower coefficient will mean less absorption / scattering. A higher coefficient will mean less. For artists, it's a bit easier to think of it in terms of distances, so we back-calculate the coefficients from distances.  
> **username_0**: No, parameter setting part is the same. I couldn't solve the problem but if I solve it, I will comment here.  

---

### Answer by username_2 on 2017-06-13. Score: 5

For cases where the diffusion approximation is preferred over full volumetric path tracing, the method published by Solid Angle is fairly efficient:
<https://www.solidangle.com/research/s2013_bssrdf_slides.pdf>

It is implemented in the Arnold render engine, Blender's Cycles, and pbrt, the latter being open source.

The file which implements it in PBRT is here:
<https://github.com/mmp/pbrt-v3/blob/7095acfd8331cdefb03fe0bcae32c2fc9bd95980/src/core/bssrdf.cpp>

> **Comments:**
> **username_0**: PBRT implementation seems like photon beam diffusion. Are photon beam diffusion and the what Solid Angle does the same? I ask because I thought PBD also take point samples on the geometry as a preprocessing step but Solid Angle slides don't do it.  
> **username_2**: There are two things to distinguish here: The function to integrate (BSSRDF) and the method of integration. The first one are functions like the dipole, multipole, photon beam diffusion, etc. The second are methods like Jensen's hierarchical octree or Solid Angle's importance sampled ray tracing.  

---

