A recent approach for subsurface scattering

Question

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?

Answer 1

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.

BSDFs

To add transmission, we first need to upgrade some of our BRDF materials to full BSDFs, with a BTDF function.

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. However, an obvious choice would be to sample based on the 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:

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);
    }
};

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).

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.

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

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.

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. 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.

Answer 2

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