Simple GPU Path Tracing, Part. 3.0 : Path Tracing Basics

 

Now that we can ray trace our scene, we're ready to start implementing a path tracing algorithm!

In this post we will be implementing a very naive approach and the result will still be far from our goal of photorealistic imagery, but we will still make some progress.

The code for this post will be on this branch of the github repo.
 

 Some theory

Let's first quickly review some overly simplified theory around how we will be calculating lighting in the scene with path tracing.

In a scene with a light, we can think of light as millions of rays that are emitted from the light source and that scatter in the scene : 

 Then, some of those rays will escape from the surface they hit, and they will get reflected into other directions : 


 

Note that each ray has a colour here. The ray that's coming from the light is white, which is the colour of the light. Then, when it hits the red box, all the rays that will get reflected from that surface will become red too, because the ray's colour gets multiplied by the colour of the shape it hits, so (1,1,1) * (1, 0, 0) becomes (1,0,0).

That's a big (huge) simplification, but it helps understanding that the rays of light change colour as they bounce across the scene. 

Some of those reflected red rays will eventually hit a camera, and form a red colour on the camera film. 

Now, if we were to generate lots of ray from the light source, and bounce them around the scene until one reaches the camera, we would have to wait for a very long time for the camera film to be filled with colours, because only a small portions of all the light rays eventually hit the camera.

Instead, in path tracing, we do the same process backwards : we generate lots of rays from the camera and we check if they intersect the scene. When an intersection is found, we want to know if it's being illuminated by a light, and what colour that shape is at this intersection. We can then do the same calculation as described before, where we multiply the colour coming from the light ray with the colour of the shape.

 So you can think of the ray coming from the camera as a black ray, meaning it's not transporting any light on it. If it hits a red shape, (that would be intersection #1), it would still be black, as the red shape doesn't emit any light as far as we know. But if we shoot another ray from that red shape, and that ray hits the light (intersection #2), then we know that intersection #1 on the red shape can be illuminated by that light, and we can apply the colour multiplication, and the camera ray will become red.

So we can return that red colour to the camera, and set the pixel to red on the film.

 

Again, that's very simplified. In theory, when we reach intersection #1, we want to shoot rays in the entire hemisphere around the surface, to know how much light comes from all directions to that point.

The only problem is there doesn't exist a way of calculating that mathematically.

I'm not going to go into too much details, there are lots of resources that explain the low level maths much better than I could possibly do, but basically that process of calculating how much light comes in to a specific point boils down to the calculation of an integral over the hemisphere surrounding the point.

This integral is called the "Rendering Equation", and doesn't have a closed form, meaning we cannot calculate it. Here's the formula of this equation : 

Let's break down each term of that equation : 

  •  L(x, wo) : How much of the light coming in the position X in the scene reflects to the direction Wo (Direction from this point to the camera), i.e what's the value of the pixel whose ray hits the scene at position x.

You can see that as I said, it's an integral over a hemisphere. Let's analyse the 4 terms of this integral : 

  • f(x, wo, wi) : This is called the BRDF (Bidirectionnal reflectance distribution function). It's a function that takes 2 input directions, wo and wi. What it calculates is, at the given point x, how much of the light coming from direction wi is reflected to direction wo. We will see more on that later when we implement BRDF functions, but essentially it defines how a material reacts to lighting. A mirror BRDF would only reflect light to a very specific direction (the "perfect" reflection direction). A diffuse or lambertian brdf would reflect incoming light to all directions in the hemisphere. You can learn more about brdf's on this course by Jakub Boksansky
  • L(x, wi) : This is the rendering equation itself ! And it's the reason why there's no closed form to this integral, it's because it's a kind of recursive integral. It accounts for how much light is coming to this point x from the direction wi, which is the direction that we're integrating with.
  • cos(theta): That's to account for a physics rule that states that the the more light arrives perpendicular to a surface, the higher the light intensity will be. If the light hits a surface at a grazing angle, it will less intense. Here's an illustration of that : 

If the light hits perpendicularly to the surface, you see that all the power of the light is concentrated in a small area part of the surface : 

 

contrarily, if the light hits at an angle, the power of the light becomes spread on a larger surface, meaning each point on that surface receives smaller amounts of light : 

 

  • d(w) : This represents the solid angle over which we're integrating.

But again, we can't solve this integral analytically :  What we can do though is using something called Monte-Carlo integration to approximate the result of this integral, by taking random numbers, also called "samples", evaluating the functions at those samples, and summing all the samples to obtain an approximation of the result.

This website has an in-depth explanation of how and why this works.

In practice, what happens is that when we shade the point X, and we want to calculate how much light comes to this point from the surrounding scene (L(X, wi), we will shoot many rays in all directions around this point, and evaluate how much light is coming. This is how Monte-Carlo integration works : To evaluate an integral, we take many many samples of that function and kind of compute an average of all those values.

In this post, we will use a much simpler version of the rendering equation. We will not use the cos(theta) term, and the BRDF will always be equal to the colour of the shape.

Let's build up an example to illustrate all of that : 

As before, we will be generating one ray for each pixel in the image.

we will test if the ray intersect the scene, keeping track of 2 colour values : 

- Radiance, which starts at (0,0,0)

- Attenuation, which starts at (1,1,1)

If the ray intersects an object, we will modulate the attenuation with the colour of that object :

So in this case, the attenuation becomes (0.8, 0.8, 0.8) (the colour of the floor) and the radiance still is (0,0,0)

We then trace a ray from this intersection point to another random location, and do the same :

 

If it hits an object, we modulate the attenuation with the colour of the this object : 

The attenuation becomes ((0.8, 0.8, 0.8) * (0.9, 0.1, 0.1) = (0.72, 0.08, 0.08) and the radiance still is (0,0,0) 

We then shoot another random ray from the new position, which is going to hit the sky : 


If it hits the sky, we set the global radiance variable to emission of the sky times the current attenuation, and we step out of the ray recursion. We set the value of the pixel to this radiance value. 

Here, let's say that the sky's emission is (1,1,1), which means it's emitting a white colour from all directions.

therefore, Radiance = (0.72, 0.08, 0.08) * (1,1,1)

the pixel value becomes (0.72, 0.08, 0.08).

 

But here, the floor is supposed to be white, why is the pixel red then ? 

Well, that's just one random sample. Another random sample could be : 

Ray hits the floor, Attenuation = (0.8, 0.8, 0.8)

Then, the secondary ray hits the sky directly, which means Radiance = (0.8, 0.8, 0.8) * (1,1,1) = ( (0.8, 0.8, 0.8)

When we accumulate a lot of samples, some will be red because they hit the red wall next to the floor, and some will be white because they hit the sky. Because we blend the result of each sample together, this will give the floor at this position a red tinted white colour, which makes sense because the wall kinda reflects its own colour on the floor.

 

Here's an outline of the algorithm that we will implement : 

For each pixel 

    Ray = GenerateRay(pixel)

    Radiance = 0

    Attenuation = 1

    for(bounce in maxBounces):

        intersection = IntersectScene(ray)

        if(!intersection.hit):

            Radiance += EnvLight()

            break out of the loop;

         else:

            NextDir = RandVectorInHemisphere()

            Attenuation *=HitShapeColour

            Ray = {HitPosition, NextDir}

    Image[Pixel] = Radiance


We will repeat this process many times, and the image will slowly converge to a good looking render.

Note that this is not a physically correct computation, but it's a good start and will give us some hints on how to improve the system later on.

Implementation

So, let's go through the code and explain each step : 

So at the start, we just calculate UV coordinates based on the thread ID. We will use that for primary ray generation.

MAIN()
{
    INIT()
   
    ivec2 ImageSize = IMAGE_SIZE(RenderImage);
    int Width = ImageSize.x;
    int Height = ImageSize.y;

    uvec2 GlobalID = GLOBAL_ID();
    vec2 UV = vec2(GlobalID) / vec2(ImageSize);
    UV.y = 1 - UV.y;
       
    if (GlobalID.x < Width && GlobalID.y < Height) {

We then enter a loop that's gonna run 100 times. For each sample, we generate a ray with the same function as before. We then initialize the Radiance and the Weight (Attenuation).

        vec4 PrevCol = imageLoad(RenderImage, ivec2(GlobalID));
        vec4 NewCol;
        for(int Sample=0; Sample < 100; Sample++)
        {
            ray Ray = GetRay(UV);
           

            vec3 Radiance = vec3(0,0,0);
            vec3 Weight = vec3(1,1,1);
       

 Then, we'll enter the tracing part : 

We need define a maximum amount of bounces, here we set that to 3.

We then test the intersection of the ray with the scene. If it doesn't hit something, it means that it hit the sky (Or the environment). 

Here we set the sky to have an emission value of 1, which is a white sky.

So when the ray hits the sky, we add the sky emission to the radiance of this pixel, multiplied by the current attenuation (or weight).

When we hit the sky, we don't wanna bounce again, so we can break out of the loop.

            for(int Bounce=0; Bounce < 3; Bounce++)
            {
                sceneIntersection Isect;
                Isect.Distance = 1e30f;
                Isect.RandomState = CreateRNG(uint(uint(GLOBAL_ID().x) * uint(1973)  
                    + uint(GLOBAL_ID().y) * uint(9277)  +  uint(Bounce + Sample) *  
                    uint(117191)) | uint(1), 371213);

                IntersectTLAS(Ray, Isect);
                if(Isect.Distance == 1e30f)
                {
                     // Environment radiance 
                    Radiance += Weight * vec3(1);
                    break;
                }

 

If we did hit a scene object, we need to calculate some values to keep going.

We first find the transform of the instance that we hit. We then get the normal of the point that we hit and transform it with the instance transform.
We then find the colour of the shape that we hit.

Finally, we calculate the position of the shape that we hit.

Note that all those calculations are based on the barycentric coordinates that allow us to interpolate any value accross the surface of a triangle.

uint Element = Isect.PrimitiveIndex;
triangle Tri = TriangleBuffer[Element];
triangleExtraData ExtraData = TriangleExBuffer[Element];    
Isect.InstanceTransform = TLASInstancesBuffer[Isect.InstanceIndex].Transform;
               
mat4 NormalTransform = TLASInstancesBuffer[Isect.InstanceIndex].NormalTransform;
vec3 Normal = ExtraData.Normal1 * Isect.U +  
              ExtraData.Normal2 * Isect.V +
              ExtraData.Normal0 * (1 - Isect.U - Isect.V);
Isect.Normal = TransformDirection(NormalTransform, Normal);
               
vec4 Colour = ExtraData.Colour1 * Isect.U + 
              ExtraData.Colour2 * Isect.V + 
              ExtraData.Colour0 * (1 - Isect.U - Isect.V);                

vec3 Position = Tri.v1 * Isect.U + Tri.v2 * Isect.V + Tri.v0 * (1 - Isect.U - Isect.V);

Next, we can modulate the weight by the colour of the shape :

Weight *= vec3(Colour);

and finally, we can start building the next ray. The ray origin will be the point we hit on the shape. 

Now we need to calculate the next ray direction. 

We need to consider the light incoming from all directions, but that's quite a lot of directions!

So for now we'll just generate a random direction in the upper hemisphere, with higher chance for directions closer to the normal. That is called a cosine distribution, because the probability of sampling a direction is proportional to the cosine between that direction and the normal.

We will heavily use that function later, so may as well define it now.

vec3 OutgoingDir = -Ray.Direction;
vec3 Incoming = SampleHemisphereCosine(Normal, Random2F(Isect.RandomState));
               
Ray.Origin = Position;
Ray.Direction = Incoming;

 

Once we've gone through all the bounces, we step out of the bounces loop, and we can calculate the new colour for the current pixel : 

float SampleWeight = 1.0f / (float(Sample) + 1);

NewCol = mix(PrevCol, vec4(Radiance.x, Radiance.y, Radiance.z, 1.0f), SampleWeight);
PrevCol = NewCol;

 We're blending between the previous value of the pixel, and the newly calculated value. mix() actually does a linear interpolation between 2 values.

The weight of the interpolation is calculated with SampleWeight : For example in the first sampe, SampleWeight will be 1, so we're taking the entire value of new colour.

In the second sample, SampleWeight is 1 / 2 = 0.5. We take half of the previous colour, and add half of the new colour. and so on...

We also set PreviousColour to the new colour.

after we've finished the sample loop, we store the NewColour value in the image, and finally exit the kernel.

imageStore(RenderImage, ivec2(GlobalID), NewCol);

 

And let's look at the result : 


So... That's a noisy and not so great image, but it's a start !

We can guess the shadows behind the cubes. we can see the green and red walls reflecting on the sides of the boxes.


Progressive rendering

At the moment, we're running 100 samples in a single kernel. The resulting image is always exactly the same at each kernel pass, because we're using deterministic random number generation. It's also very slow because it takes quite a long time to process 100 samples.

What we want is the image to converge over time, and we want to run only a few samples in a kernel pass so that the application can become interactive. 
We will then run the kernel over and over, and we will seed the random number generator with the current sample number, so that each time we have a new set of random numbers.

In order to achieve that, we will pass to the kernel some information : 
- How many samples does the kernel must run in one pass
- The value of the first sample it must process. 

To send this information to the kernel, we will use a new buffer that will store informations that will be needed by the path tracing algorithm:
 

struct tracingParameters
{
    int CurrentSample;
    int TotalSamples;
    int Batch;
    int Pad0;    
};

inline tracingParameters GetTracingParameters()
{
    tracingParameters Params;
    Params.CurrentSample = 0;
    Params.Batch = 1;
    Params.TotalSamples = 256;
    return Params;
}
 
We call this function in application::Init() : 
    Params =  GetTracingParameters();  
 
Inside of InitGpuObjects(), we then create a buffer containing this structure : 
#if API==API_GL
    PathTracingShader = std::make_shared<shaderGL>("resources/shaders/PathTrace.glsl");
    RenderTexture = std::make_shared<textureGL>(Window->Width, Window->Height, 4);
    TracingParamsBuffer = std::make_shared<uniformBufferGL>(sizeof(tracingParameters), 
                                                             &Params);
#elif API==API_CU
    RenderTexture = std::make_shared<textureGL>(Window->Width, Window->Height, 4);
    RenderBuffer = std::make_shared<bufferCu>(Window->Width * Window->Height * 4 
                                                                     * sizeof(float));
    RenderTextureMapping = CreateMapping(RenderTexture);    
    TracingParamsBuffer = std::make_shared<bufferCu>(sizeof(tracingParameters), &Params);
#endif
 
note that for OpenGL, we use a new class called uniformBufferGL. that's because we will be using uniform buffers instead of shader storage buffers, that are more suited for passing small structs to the gpu. I've added it to the bufferGL file. Again, it's a very simple wrapper.

I've now refactored all the kernel execution code into a Trace() function that we will run every frame.
Inside the function, we run the kernel and increment the CurrentSample value of the struct, and upload again to the gpu.
We also need to pass the tracingParamsBuffer to the kernels : 
void application::Trace()
{
    if(Params.CurrentSample < Params.TotalSamples)
    {
        TracingParamsBuffer->updateData(&Params, sizeof(tracingParameters));
       
    #if API==API_GL
        PathTracingShader->Use();
        PathTracingShader->SetTexture(0, RenderTexture->TextureID, GL_READ_WRITE);
        PathTracingShader->SetSSBO(BVH->TrianglesBuffer, 1);
        PathTracingShader->SetSSBO(BVH->TrianglesExBuffer, 2);
        PathTracingShader->SetSSBO(BVH->BVHBuffer, 3);
        PathTracingShader->SetSSBO(BVH->IndicesBuffer, 4);
        PathTracingShader->SetSSBO(BVH->IndexDataBuffer, 5);
        PathTracingShader->SetSSBO(BVH->TLASInstancesBuffer, 6);
        PathTracingShader->SetSSBO(BVH->TLASNodeBuffer, 7);        
        PathTracingShader->SetSSBO(Scene->CamerasBuffer, 8);
        PathTracingShader->SetUBO(TracingParamsBuffer, 9);
        PathTracingShader->Dispatch(Window->Width / 16 + 1, Window->Height / 16 +1, 1);
#elif API==API_CU
        dim3 blockSize(16, 16);
        dim3 gridSize((Window->Width / blockSize.x)+1, (Window->Height / blockSize.y) + 1);
        TraceKernel<<<gridSize, blockSize>>>((glm::vec4*)RenderBuffer->Data, Window->Width, 
        Window->Height,(triangle*)BVH->TrianglesBuffer->Data, 
        (triangleExtraData*) BVH->TrianglesExBuffer->Data, (bvhNode*) BVH->BVHBuffer->Data, 
        (u32*) BVH->IndicesBuffer->Data, (indexData*) BVH->IndexDataBuffer->Data, 
        (bvhInstance*)BVH->TLASInstancesBuffer->Data, (tlasNode*) BVH->TLASNodeBuffer->Data,
(camera*)Scene->CamerasBuffer->Data, (tracingParameters*)TracingParamsBuffer->Data);
         
        cudaMemcpyToArray(RenderTextureMapping->CudaTextureArray, 0, 0,  
                RenderBuffer->Data, Window->Width * Window->Height * sizeof(glm::vec4),
                cudaMemcpyDeviceToDevice);
#endif
        Params.CurrentSample+= Params.Batch;
    }    
}
 
 Now in the PathTracingCode, we can use those informations : 
the samples loop now becomes : 
for(int Sample=0; Sample < Parameters.Batch; Sample++)

in the random state creation, we now calculate the global sample number : 
Isect.RandomState = CreateRNG(uint(uint(GLOBAL_ID().x) * uint(1973)  
                    + uint(GLOBAL_ID().y) * uint(9277)  +  
                     uint(Bounce + GET_ATTR(Parameters,CurrentSample) + Sample)
 * uint(117191)) | uint(1), 371213);

And when we calculate the weight of the current frame, we use that as well : 
float SampleWeight = 1.0f / (float(GET_ATTR(Parameters,CurrentSample) + Sample) + 1);

And that's pretty much it, now we have a nice progressive rendering of the scene.

Links

 

Commentaires

Enregistrer un commentaire

Articles les plus consultés