Simple GPU Path Tracing, Part 7.1 : Volumetric materials

 

In this post, we will discuss the implementation of volumetric materials into our path tracer.

 


Here's the commit for this post.

Basic theory

A volumetric material is a material which contains "participating media", meaning it's not an empty, it contains particles. Think of it as some kind of fog or cloud that's filling up the shape.

The way it works is when a ray hits the shape, instead of getting reflected out in a direction in the hemisphere around the normal, it's getting inside of the shape ! and then all sorts of crazy things can happen, and we will get into those details later.

But what's important to understand is volumetric materials don't render the surface of the object, but rather their inside volume.


Here's a schematic explanation of what we'll be doing : 


 

We shoot a ray, and as it hits a volumetric shape, we start entering the shape and evaluating the colour of the volume at multiple positions inside the shape, and use that to modulate our weight.

Note that in the specific case that's depicted in the schema, there's no radiance accumulated within the volume (Obviously, the material is not emitting any light). The radiance will only start to accumulate when we get out of the volume on the other end of the shape.

 Code set up

First thing we need to do is add some fields into our material structs and add a new volumetric material type : 
#define MATERIAL_TYPE_VOLUMETRIC   2

And in the material and materialPoint structs :
glm::vec3 ScatteringColour = {};
float TransmissionDepth = 0.01f;

Next in our scene definition, we will get back to our old cornell box, and set the short box to be a volumetric object : 
ShortBoxMaterial.Colour = {0.8, 0.8, 0.8};
ShortBoxMaterial.MaterialType = MATERIAL_TYPE_VOLUMETRIC;
ShortBoxMaterial.ScatteringColour = {0.9, 0.2, 0.4};

we'll also add those 2 new parameters in the material gui : 
    if(Mat.MaterialType == MATERIAL_TYPE_VOLUMETRIC)
    {
        ImGui::Separator();
        Changed |= ImGui::DragFloat("TrDepth", &Mat.TransmissionDepth, 0.001f, 0, 1);
        Changed |= ImGui::ColorEdit3("Scattering", &Mat.ScatteringColour[0]);
    }

Now we're ready to change our path tracing code to take all that stuff into account. There will be quite a lot of changes in the path tracing function, but I'll explain everything.
 

Volumetric rendering done horribly wrong

The goal of this part is just to set a general architecture for volumetric rendering, and we will take lots of shortcuts. The result will not be a pretty image, but afterwards we will just have to change a few bits of code to make it better.

First thing we need to do is figure out what's the next ray direction once we hit a volumetric material.
Well, as I said, we are continuing the ray within the shape. 
 
But there's one subtlety here, which is that on the surface of the shape, the light is always going to to go in the same direction, both when it enters and when it escapes the media. Therefore it can be represented as a "delta" material, meaning its bsdf is a delta function : it only has a value of 1 at a specific direction, and no value at all elsewhere. When we use a delta material, the main difference is that we don't need to sample the lights, as we only need and want to sample the only direction that matters for the bsdf. Let's handle this case now : 
vec3 Incoming = vec3(0);
if(!IsDelta(Material))
{
    if(RandomUnilateral(Isect.RandomState) < 0.5f)
    {
        Incoming = SampleBSDFCos(Material, Normal, OutgoingDir,  
                RandomUnilateral(Isect.RandomState), Random2F(Isect.RandomState));
    }
    else
    {
        Incoming = SampleLights(Position, RandomUnilateral(Isect.RandomState),  
                    RandomUnilateral(Isect.RandomState), Random2F(Isect.RandomState));
    }
    if(Incoming == vec3(0,0,0)) break;
    Weight *= EvalBSDFCos(Material, Normal, OutgoingDir, Incoming) /
            vec3(0.5 * SampleBSDFCosPDF(Material, Normal, OutgoingDir, Incoming) + 
                 0.5f * SampleLightsPDF(Position, Incoming));
}
else
{
    Incoming = SampleDelta(Material, Normal, OutgoingDir, 
                     RandomUnilateral(Isect.RandomState));
    Weight *= EvalDelta(Material, Normal, OutgoingDir, Incoming) /
            SampleDeltaPDF(Material, Normal, OutgoingDir, Incoming);                        
}
Now the volumetric is considered a delta material when the ray hits it, because there's only 1 direction it can go to : the same direction as the previous ray.

Let's now see how the delta function evaluations look : 
FN_DECL vec3 EvalDelta(INOUT(materialPoint) Material, vec3 Normal, vec3 OutgoingDir, vec3 Incoming)
{
    if(Material.MaterialType == MATERIAL_TYPE_VOLUMETRIC)
    {
        return EvalVolumetric(Material.Colour, Normal, OutgoingDir, Incoming);
    }    
}

FN_DECL float SampleDeltaPDF(INOUT(materialPoint) Material, INOUT(vec3) Normal, 
                         INOUT(vec3) OutgoingDir, INOUT(vec3) Incoming)
{
    if(Material.MaterialType == MATERIAL_TYPE_VOLUMETRIC)
    {
        return SampleVolumetricPDF(Material.Colour, Normal, OutgoingDir, Incoming);
    }    
}

FN_DECL vec3 SampleDelta(INOUT(materialPoint) Material, INOUT(vec3) Normal,  
                    INOUT(vec3) OutgoingDir, float RNL)
{
    if(Material.MaterialType == MATERIAL_TYPE_VOLUMETRIC)
    {
        return SampleVolumetric(OutgoingDir);
    }    
}
Quite simple, as we only have one delta material so far.
 

Let's now review the 3 volumetric functions :
FN_DECL vec3 EvalVolumetric(vec3 Colour, vec3 Normal, vec3 Outgoing, vec3 Incoming)
{
    // If Incoming and outgoing are in the same direction, return 0
    // For a volume that's not the case.
    if(dot(Normal, Incoming) * dot(Normal, Outgoing) >= 0)
    {
        return vec3(0);
    }
    else
    {
        return vec3(1);
    }
}

FN_DECL float SampleVolumetricPDF(vec3 Colour, vec3 Normal, vec3 Outgoing, vec3 Incoming)
{
    if(dot(Normal, Incoming) * dot(Normal, Outgoing) >= 0)
    {
        return 0;
    }
    else
    {
        return 1;
    }
}

FN_DECL vec3 SampleVolumetric(vec3 OutgoingDir)
{
    return -OutgoingDir;
}

So as you can see, when a ray hits a volumetric material, the ray is just going to keep going inside of the material with the exact same direction, with a probability of 1.
 
Then, we will officially enter a volumetric material, which is going to lead us into another branch of the path tracing algorithm specially fitted for volumes.
 
After we have hit a shape and calculated the weight of the ray and the next direction, we check if the hit material is volumetric and set some variables accordingly :
//If the hit material is volumetric
// And the ray keeps going in the same direction (It always does for volumetric materials)
// we set the volume material to be the current material
if(IsVolumetric(Material) && dot(Normal, OutgoingDir) * dot(Normal, Incoming) < 0)
{
    VolumeMaterial = Material;
    HasVolumeMaterial = !HasVolumeMaterial;
}
 
The IsVolumetric function is simple, but might evolve later if we add other volumetric materials : 

FN_DECL bool IsVolumetric(INOUT(materialPoint) Material)
{
    return Material.MaterialType == MATERIAL_TYPE_VOLUMETRIC;
}
 
Fine. Now, we have a way of keeping track wether we're inside a volume or not, and of the material of that volume.
We will then continue the bounces loop and get to the next ray, and we want to use this variable to check if we're in a volume. If we are, we will advance in the volume in a certain direction, and evaluate the colour at this position in the volume : 
bool StayInVolume=false;
if(HasVolumeMaterial)
{
    // If we're in a volume, we sample the distance that the ray's gonna intersect.
    // The transmittance is based on the colour of the object. The higher, the thicker.
    float Distance = VolumeMaterial.TransmissionDepth;
    Weight *= VolumeMaterial.Colour;

    //If the distance is higher than the next intersection, it means that we're stepping out of the volume
    StayInVolume = Distance < Isect.Distance;

    Isect.Distance = Distance;
}

For now, we advance by TransmissionDepth, and use the scatteringColour to set the colour of the volume, but this will change later.
Because we traced a ray from inside the volume, Isect.Distance contains the distance from inside the volume to the other side. So if distance is less than the intersection distance, we stay in the volume. Otherwise, it means we're getting out : 
 
 
 
Then, depending on wether we're inside a volume or not, we will take 2 different branches.
If we're not in, we'll use the same code as before to find the next ray direction, calculate the BRDF and the PDF associated with this direction, and update the weight and radiance.
Otherwise, here's what we will do : 

else
{
    vec3 Outgoing = -Ray.Direction;
    vec3 Position = Ray.Origin + Ray.Direction * Isect.Distance;

    vec3 Incoming = vec3(0);
    if(RandomUnilateral(Isect.RandomState) < 0.5f)
    {
        // Sample a scattering direction inside the volume using the phase function
        Incoming = -Outgoing;
    }
    else
    {
        Incoming = SampleLights(Position, RandomUnilateral(Isect.RandomState), RandomUnilateral(Isect.RandomState), Random2F(Isect.RandomState));                
    }

    if(Incoming == vec3(0)) break;

    Weight *= VolumeMaterial.ScatteringColour;

    Ray.Origin = Position;
    Ray.Direction = Incoming;                    
}

We calculate a new incoming direction a bit like before, except that instead of sampling the bsdf, we just continue in the same direction as before.
then, we just modulate the weight by the volume scattering colour, and that's all.

This is going to create quite an ugly looking image, but that was just to do the global set up of volumetric rendering. Next we'll dive into more physically correct equations and we will get to better looking images.

Getting (more) physically correct

Now that we've got the basic set up for path tracing volumetric materials, we can start getting more physically correct to get better looking results.

First thing we'll change is how we calculate the distance that we're travelling inside the volume, and to do that we will be using the Beer-Lambert law that determines the attenuation of light as it passes through a medium, or what's called "transmission". What we want to know, when evaluating transmission, is how much light gets tramsmitted through a volume, meaning the amount of light that doesn't get absorbed by the particles within the volume.

To do that, we need to define a density value (or absorption coefficient) for this medium which we can calculate as follows 


As you see, we divide the log term by our transmissionDepth value which is expressed as a distance.
The unit of that term is therefore inverse meters. It gives the probability that a photon gets absorbed or scattered, given the colour and the transmission depth.
If we take the inverse of that term, we get what's called the "mean free path", which corresponds to the average distance a photon can travel through a medium before an interaction happens.

Then, we can use Beer's law to calculate how light attenuates inside the medium : 

As you see, there's an exponential dependence between the transmittance of light in a volume, and the distance * density value.
 
So let's add that into our code in EvalMaterial : 
Point.Density = vec3(0,0,0);
if(Material.MaterialType == MATERIAL_TYPE_VOLUMETRIC)
{
    Point.Density = -log(clamp(Point.Colour, 0.0001f, 1.0f)) / Point.TransmissionDepth;
}
 
 
and then in the code where we march in the volume, we'll now use this : 
 
if(HasVolumeMaterial)
{
    // If we're in a volume, we sample the distance that the ray's gonna intersect.
    float Distance = SampleTransmittance(VolumeMaterial.Density, Isect.Distance, RandomUnilateral(Isect.RandomState), RandomUnilateral(Isect.RandomState));
    Weight *= EvalTransmittance(VolumeMaterial.Density, Distance) /
            SampleTransmittancePDF(VolumeMaterial.Density, Distance, Isect.Distance);
   
   
    //If the distance is higher than the next intersection, it means that we're stepping out of the volume
    InVolume = Distance < Isect.Distance;

    Isect.Distance = Distance;
}
 
 Let's review the 3 functions we're using here, with the more in depth explanations in comments :
 
FN_DECL float SampleTransmittance(vec3 Density, float MaxDistance, float RL, float RD)
{
    // Choose a random channel to use
    int Channel = clamp(int(RL * 3), 0, 2);
   
    // Here we calculate the distance that a ray traverses inside a medium. 
    // We sample this distance using the exponential function,
    // using the mean free path 
    // (=average distance a photon can travel through a medium before an interaction, 
    // absorbtion or scattering, happens.)
    //Here, the mean free path is simply 1 / density

    // Calculate the distance we travel, 
    // using the inverse of the CDF of the exponential function * mean free path.
    float Distance = (Density[Channel] == 0) ? 1e30f : -log(1 - RD) / Density[Channel];
   
    return min(Distance, MaxDistance);
}

FN_DECL vec3 EvalTransmittance(vec3 Density, float Distance)
{
    // Beer-Lambert law : attenuation of light as it passes through a medium, 
    // as a function of the extinction coefficient and of the distance 
    // travelled inside the medium.
    return exp(-Density * Distance);
}

FN_DECL float SampleTransmittancePDF(vec3 Density, float Distance, float MaxDistance)
{
    // We use the pdf of the exponential distribution, 
    // because we're sampling distances with the exponential distribution.
    // the pdf is pdf(x, delta) = delta * exp(-delta * x)
    // x = distance 
    // delta = rate parameter, inverse of the mean free path, 
    // so delta = density.

    if(Distance < MaxDistance)
    {
        return Sum(Density * exp(-Density * Distance)) / 3;
    }
    else
    {
        return Sum(exp(-Density * MaxDistance)) / 3;
    }
}
 
Cool, now we have a better physically correct way of marching inside a volume, but we're not done yet !
We're still not using a bsdf when generating the next ray directions, so let's fix this.
Instead of setting the next ray direction as the -OutgoingDir, we'll sample from what's called a "phase" function : 
Incoming = SamplePhase(VolumeMaterial, Outgoing, RandomUnilateral(Isect.RandomState), 
                       Random2F(Isect.RandomState));
 
And we will evaluate this bsdf when modulating the weight : 

Weight *= EvalPhase(VolumeMaterial, Outgoing, Incoming) /
(
0.5f * SamplePhasePDF(VolumeMaterial, Outgoing, Incoming) +
0.5f * SampleLightsPDF(Position, Incoming)
);
 
A Phase function essentially determines the probability of scattering angles inside a volumetric medium. It's like a BRDF, except it works for the entire sphere around a point, rather than just a hemisphere. It returns how much light is scattered, given a view direction and an incoming light direction. As for BRDF, it integrates to 1 over the entire sphere, and it is reciprocal (we can interchange the 2 input vectors and still get the same result).

We will more specifically use the Henyey-Greenstein phase function which is quite popular in computer graphics.
In this model, there's an "Isotropy" parameter that determines the shape of the distribution around the sphere, a bit like the "roughness" parameter does for our BRDF. 
An isotropic phase function will scatter light uniformly in all directions around the sphere, whereas an anisotropic material will scatter light in a specific direction.

Here's the phase function equation : 
 
Where g corresponds to our Anisotropy parameter. If g is > 0, the light is scattered forward, and if g < 0, the light is scattered backwards. Note that if g = 0, the equation resolves to 1 / (4 * pi), which corresponds to a uniform distribution over an entire sphere, like an isotropic volume.
 
And there's its implementation :
 
FN_DECL vec3 EvalPhase(INOUT(materialPoint) Material, vec3 Outgoing, vec3 Incoming)
{
    if(Material.Density == vec3(0)) return vec3(0);

    float Cosine = -dot(Outgoing, Incoming);
    float Denom = pow(1 + Material.Anisotropy * Material.Anisotropy - 2 * Material.Anisotropy * Cosine, 1.5f);
    float PhaseFunction = (1 - Material.Anisotropy * Material.Anisotropy) / (4 * PI_F * Denom * sqrt(Denom));

    return Material.ScatteringColour * Material.Density * PhaseFunction;
}
 
Then, let's look at the sampling function.We don't sample a direction directly, but rather a cos(theta) value that we then use to infer a direction, using a random phi value between 0 and 2*PI.
Here's the formula :

Note that we're using the inversion technique to sample from that function, by sampling the inverse of the phase function with a random number.


FN_DECL vec3 SamplePhase(INOUT(materialPoint) Material, vec3 Outgoing, float RNL, vec2 RN)
{
    if(Material.Density == vec3(0)) return vec3(0);
    float CosTheta = 0;
    if(abs(Material.Anisotropy) < 1e-3f)
    {
        CosTheta = 1 - 2 * RN.y;
    }
    else
    {
        float Square = (1 - Material.Anisotropy * Material.Anisotropy) /
                       (1 + Material.Anisotropy - 2 * Material.Anisotropy * RN.y);
        CosTheta = (1 + Material.Anisotropy * Material.Anisotropy - Square * Square) /
                   (2 * Material.Anisotropy);
    }

    float SinTheta = sqrt(max(0.0f, 1- CosTheta * CosTheta));
    float Phi = 2 * PI_F * RN.x;
    vec3 LocalIncoming = vec3(SinTheta * cos(Phi), SinTheta * sin(Phi), CosTheta);
    return BasisFromZ(-Outgoing) * LocalIncoming;
}

And lastly, we need to implement the pdf function, which is simply equal to the phase function : 
FN_DECL float SamplePhasePDF(INOUT(materialPoint) Material, vec3 Outgoing, vec3 Incoming)
{
    if(Material.Density == vec3(0)) return 0;

    float Cosine = -dot(Outgoing, Incoming);
    float Denom = 1 + Material.Anisotropy * Material.Anisotropy - 2 * Material.Anisotropy * Cosine;

    return (1 - Material.Anisotropy * Material.Anisotropy) / (4 * PI_F * Denom * sqrt(Denom));

}
 
 and we're done ! We can now render volumetric objects in our path tracer, here's the result of our scene : 

 
 
Nice ! Next, we'll add another volumetric material. The material we've been using here doesn't have a surface to it, all the rays go through it and none of them reflects from it surface. We will implement a material that can either reflect light from its surface, and scatter light inside of it.

If you want to learn more, here are a few useful links : 
  • Volume Rendering from Scratchapixel. This is generally speaking an excellent resource for learning computer graphics in depth, it really goes in the nitty gritty details of equations and their meaning. I really recommend reading through it if you have time. 
  • The Volume Scattering a Volumetric Light Transport chapters from pbrt. As always with this book, great explanations with code.
  • Production Volume Rendering SIGGRAPH 2017 Course by Pixar. 100 pages course on the theory and implementation of volume rendering in a path tracer.

 

Next post : Simple GPU Path Tracing, Part 7.2 : Refractive material

Commentaires

Enregistrer un commentaire

Articles les plus consultés