PathTracer Part 1

In this project, I explored the core routines of a physically-based renderer through implementation of a pathtracing algorithm. I delved into concepts such as ray-scene intersection, acceleration structures, and physically-based lighting, culminating in the generation of complex, realistic pictures. Over the course of the project, I learned a lot about the ideas behind modelling and approximating real world phenomena, and how introducing a bit of complexity can drastically improve both performance and the end result.

Part 1: Ray Generation and Intersection

Ray generation begins with the eponymous generation of a sample of random rays through each pixel. We can trace the path of each these rays to determine the radiance associated with each one, which we can then average to estimate each pixel's overall irradiance.

Ray Generation Implementation Details

• Note that that all rays generated originate from the camera; therefore, given a sample point on a pixel, we must first convert it to a point on the sensor plane whose bottom left and top right corners can be calculated by

  Vector3D(-tan(radians(hFov)*.5), -tan(radians(vFov)*.5),-1)
  Vector3D(tan(radians(hFov)*.5), tan(radians(vFov)*.5),-1)

  where hFov and vFov are the field of view angles. We can then map this result from the camera space to the world space to get the direction of the ray. The origin of the ray is simply the position of the camera.

With the myriad of rays that are generated, we need an efficient way to determine intersection so we can actually calculate features such as visibility, shadows, and lighting for objects in the scene. Since each object's surface can be represented by one or a mesh of primitives, we can distill the problem of general intersection into the question of whether a ray intersects with a primitive in the scene. In this particular project, we took a look at two primitives in particular, triangles and spheres. To calculate triangle intersection, I used the Möller Trumbore algorithm. Given a triangle and a ray, the idea is as follows:

1. Using barycentric coordinates, we can parametrize the intersection point P using the vertices of a given triangle as

   P = wA + uB + vC

   Note that since w + u + v = 1, this can also be rewritten as

   P = (1 - u - v)A + uB + vC 
   = A + u(B−A) + v(C−A)

2. The intersection point P can also be written using the ray's equation as

   P = O + tD

3. Setting the two equations for the intersection equal to each other, we get

O + tD = A + u(B − A) + v(C − A)
 O - A = -tD + u(B − A) + v(C − A)

4. We can think about equation above as a change of basis to the tuv-space. The O - A term on the left can be thought as a recentering of the triangle to the new space's origin, and the terms on the left perform the transformation of the intersection point from xy-space to tuv-space.
5. Finally, note that since we're dealing with three dimensional vectors, the equation is actually a system of three linear equations, which can be solved using Cramer's rule. The resulting solution can be seen below.

Upon determining a valid intersection exists, we update the max time of the ray to reflect this new information, given that the intersection does not occur after the current max time of the ray.

Screenshots

Part 2: Bounding Volume Hierarchy

While the motivation for calculating intersections in Part 1 was well-founded, the methodology could've used a bit of improvement. Determining intersections by iterating over every single primitive in the scene just isn't sustainable for any moderately complex image. To remedy the sitaution, enter bounding volume hierarchies (BVH). At a very high level, BVH is an object partition akin to a binary tree. The recursive algorithm to construct the BVH is as follows:

1. Given a list of a primitives, start by computing their bounding box.
2. Using the computed bounding box, initialize a new BVH node.
3. If the size of the list of primitives is less than a predetermined limit, the node is a leaf node. Associate the list of primitives with the node and return the node
4. Else, we need to split the primitives into a left and right child node. We split the primitives based on the axis associated with largest dimension of the bounding box's extent.
5. The actual split point along the axis is determined by the corresponding dimension of the average centroid of each of the primitives' bounding boxes. We use the average centroid since it allows every primitive to make a contribution to help determine where the "center" is.
6. Finally, split the primitives into two lists, left and right depending on whether the value of their bounding box's relevant centroid coordinate is less than or greater than the average value computed in the previous step.
7. Finally, assign the left and right children of this BVH node to the results from recursing on the left and right lists generated.
8. Return the ensuing node.

To calculate intersections with this new data structure, we do the following:

1. To start, we first check if the ray intersects with the BVH node's bounding box at all. If not, we can easily short circuit and return false. If we image the box as the intersection of 3 slabs (one for each axis), we can calculate the intersection of a ray with said box by combining the results from running intersection tests on the three slabs. Since each slab has two faces, we'll get two values for the time t of intersection, a smaller one for the face closer to the ray, and a larger one for the face farther to the ray. After taking the maximum of the smaller values and the minimum of the larger values and making sure the max smaller value is less than or equal to the max larger value, we can assert intersection if the resulting time period overlaps with the valid time period for the ray. Intersection with face of each slab was computed as follows:
2. If we've determined that the ray intersects with the BVH node's bounding box, we have two options. If the node is a leaf, we iterate over the included primitive list and determine if any intersections occur, returning of the closest one. if no intersections occur, we return false.
3. If the node is not a leaf, we recurse on the left and right children of the node, and return the nearest of the two intersections if applicable.

Screenshots

Rendering Speed Experiments

As apparent in the sizable differences in rendering speed (22s vs 0.05s for Cow; 335s vs 0.05s vs Beast), BVH is a significant improvement over the naive linear alternative. By organizing the primitives into a binary tree like structure, we can drastically decrease our search space each layer we progress downwards (log vs linear running time). The short-circuiting ensures no time is wasted individually checking primitives when we don't have to. Notice that this all comes at the cost of increasing the pre-processing time by around 10-50x in order to build the BVH; however, relative to the decrease in rendering time, the increase in preprocessing time is trivial. Choosing the right data structure can mean the difference between waiting an afternoon or a few minutes. On a side note, there appears to be a baseline of 0.05s of rendering time for both Cow and Beast, despite beast having over 10x the primitives, which is probably just the overhead of the program.

Part 3: Direct Illumination

We preface this section by nothing that both types of direct lighting are examples of Monte Carlo Estimation, or inverse probability weighting.

Uniform Hemisphere Sampling

The core idea behind uniform hemisphere sampling is that we can estimate irradiance at a point by sampling a preset number of random solid angles over the enclosing hemisphere. Given a point-of-interest p, we perform the following steps:

1. For each of the desired samples, we begin by generating a random direction. Since each of the directions is equally likely (hence uniform hemisphere sampling), the pdf for any direction is equal to 1/2π (half of surface area of full sphere, 4π).
2. We then trace the resulting ray to check for any intersections.
3. If no intersection occurs, we move onto the next sample.
4. If an intersection does occur, we take the emission from the intersection and convert it to irradiance by projecting it onto the desired surface using the cosine of the direction vector and normal vector, dividing by the pdf, and multiplying by BSDF.
5. Finally, we average over all of the samples to get our direct lighting estimate.

Importance Sampling by Sampling over Lights

The core idea behind importance sampling is that we can use knowledge of what the distribution might look like (priors) to improve our estimate. In this particular case, we know that the only nonzero contributions to irradiance must come from light sources. Why not start from the light sources instead of randomly picking directions to check for said lights? Give a point-of-interest p, we perform the following steps:

1. For each light, we take a preset number of samples (1 sample if its a delta light, since all samples from the light will be the same).
2. Upon sampling the light passing p as a parameter, we get four quantities: the incoming radiance, the direction from p to the light source, the distance to the light, and the pdf evaluated at the given direction.
3. After appling the necessary world space to object space conversions and short-circuiting if the object lies behind the light source, we can then check for any objects in between the p and the light source using a shadow ray originating from p and heading towards the light. If an intersection exists (besides that of with the light), we move onto the next sample. Otherwise, we move to the next step.
4. If no intersection exists, we take the incoming radiance returned by the sampling of the light and convert it to irradiance by projecting it onto the desired surface using the cosine of the direction vector and normal vector, dividing by the pdf returned by the sampling of the light, and multiplying by the BSDF.
5. Finally, we average within the samples for each light, and add the results together (a final average should not be taken!).

Screenshots

General renders using global illumination

Variable light rays, 1 sample per pixel

Comparison

Due to the nature of uniform hemisphere sampling selecting samples with no prior knowledge, we get sizable amount of noise from samples that don't quite find the light source, resulting in a recognizable, but speckled scene. Importance sampling helps us mitigate a sizable amount of this effect by only focusing on rays originating from light sources making use of our priors, but as seen by the images of the dragon with respect to the number of light rays, more samples are almost always better. This latter point also applies to uniform hemisphere sampling, as seen below.

Part 4: Global Illumination

Global illumination takes what we learned in part 3 and applies it recursively to estimate higher bounces. To solve the problem of infinite bounces, we use Russian roulette as an unbiased random terminator. In addition, we incorporate zero bounce radiance in this part. My implementation of total radiance arriving at a given point p is as follows:

• For zero bounce, we just check if p emits light on its own.
• For one bounce, refer to the implementation of part 3.
• For at least one bounce or indirect lighting if max ray depth is greater than 1, the following steps are taken:
  1. We start by calculating the one bounce radiance at the point.
  2. We then take a sample from the surface BSDF, which takes the outgoing radiance direction and gives us the BSDF value, the incoming radiance direction, and the pdf associated with the incoming direction.
  3. We then use Russian roulette with a probability of 0.7 of not terminating to determine whether or not to terminate the current ray. However, we do guarantee one indirect bounce assuming there exists an intersection (the pdf for this specific bounce and direct bounce before it are equal to 1).
  4. If we do not terminate the current ray, we create a ray in the direction of the incoming radiance direction converted to the world space, check if intersects with anything, and recursively call the at least one bounce function to trace this ray.
  5. Note that each ray is associated with a depth counter, which is initialized at the max depth and counts down to 1. A ray does not continue after it has reached its max depth.
  6. If applicable, we take the incoming radiance returned by the recursive call and project it onto the desired surface using the cosine of the direction vector and normal vector, divide by both the pdf returned by the sampling of the surface BSDF and one minus the Russian roulette termination probability, and scale by the BSDF. We then add this outgoing radiance estimate to the initial one bounce radiance.
  7. Finally, we return the summed amount (or just the one bounce radiance amount if no recursive step was taken).

Screenshots

General renders using global illumination

Indirect vs direct lighting

Variable Max Ray Depth

Changes in sample rates

Part 5: Adaptive Sampling

Adaptive sampling applies the ideas behind confidence intervals to give us an estimate of pixel convergence. We can implement it for a given pixel as follows:

1. To start, we introduce the idea of convergence factor of sorts, defined below:
   
   
   
   We can determine if a pixel has converged if:
   
   
   
   where maxTolerance is set to 0.05.


2. To make it easier to work with, we can convert a Spectrum class to a illuminance by calling Spectrum::illum(). We can use a streaming algorithm as show below to keep track of just the sum and squared sum of the illuminance.


3. In order to actually calculate the mean and standard deviation, we use the following formulas:
   
   
4. Since checking convergence after every sample can be costly, we opt to check for convergence after a batch of samples (32 to be specific) have been made.
5. Finally, after convergence is achieved with 95% confidence, we return the mean of all of the gathered samples as the estimated irradiance for the pixel.

Screenshots

Slides and equations from the cs184 website.