Homework 3: Pathtracer

site: https://cal-cs184-student.github.io/hw-webpages-sp24-ashmchiu/hw3/

Overview

This homework began with generating rays and performing primitive intersections, in which we learned how to convert between camera space and world space, generated samples and intersected rays with triangles and spheres. This allowed us to progress to the next step, dealing with bounding volume hierarchies where we could greatly increase the speed of our work. By splitting primitives into BVH nodes, we no longer had to naively intersect the ray with each primitive and instead, test on broader bounding boxes–and wrap this all into one larger intersection method that was able to determine the nearest intersection for a ray within a BVH. Once our BVH was working, we could move onto direct illumination! Direct illumination, involving zero- and one-bounce illumination, allows us to show the direct emission from lights and render one-bounce illumination with uniform hemisphere sampling and importance sampling with lights. With direct illumination complete, we moved into global illumination, first writing a method to implement indirect illumination (all outgoing light past the first bounce in ray tracing). We also implemented Russian Roulette-style random termination, noting that this made previously infeasible renderings now possible by randomly terminating the recursive tracing. Finally, in an effort to dynamically determine how many times we need to sample a pixel, we implemented adaptive sampling, which, with the ability to calculate which pixels took longer to converge, would spend more samples on those difficult to render (or important) locations and sample less on those that converged faster.

We found it interesting how now, beyond just pure two-dimensional and three-dimensional work, we now were able to show how light interacts with a scene. The difference between direct and indirect lighting was particularly illuminating: prior to learning about the Monte Carlo estimations and how multi-bounce inverse paths worked, we never thought about how light could bounce back and forth between surfaces, actually lightening areas that smoothed out shadows. It’s been really rewarding to be able to render realistically lit scenes as opposed to the pure rasterizations we previously performed for two-dimensional shapes in Homework 1 and smooth meshes in Homework 2. Our images are no longer purely of one two- or three-dimensional object, but rather a more realistic interaction of lights and textures.

A few difficult debugging journeys we crossed were:

• Part 2 Task 1: we originally initialized two new vectors of Primitive *s, but the memory management here was difficult as we debugged and learned that the relative addressing wasn’t working the way it was. We eventually stumbled upon the std::stable_partition method that luckily worked on std::vectors, and this not only worked, but simplified a lot of the code we were originally had to traverse through the list of Primitive *s. The only difficulty was coming up with the syntax for the lambda function to split the groups–we hadn’t worked with lambda functions in C++ yet, so we had to learn to capture variables by reference if they weren’t going to be parameters into our lambda function.
• Part 3 Task 3: We initially thought that our issue was that our zero-bounce (get_emission) was always returning Vector3D(0, 0, 0). We spent a while debugging this, only to learn that we had accidentally overridden the Intersection parameter with another Intersection variable of the same name, isect. Lesson learned to keep track of accidentally overwriting variables :’)
• Part 4 Task 2: Namely, to add non-accumulated rendering, we were having trouble with outgoing light from previous bounces still showing up in our renders. We originally we believed we always needed one_bounce_radiance to be called to simulate the bounce, but we realized that we only needed to ensure that in our L_out calculation to not accumulate: we didn’t need to call one_bounce_radiance unless it was the last bounce we were performing–and that was the L_out we returned.

Part 1: Ray Generation and Scene Intersection

Let’s utilize Task 1 and Task 2 to discuss the ray generation and primitive intersection parts of the rendering pipeline. We will also explain primitive intersections and sphere intersections in Task 4.

Then, we’ll continue onto Task 3 to explain the triangle intersection algorithm.

Task 1: Generating Camera Rays

We first want to transform the image coordinates (x, y) to camera space by interpolating. Knowing that we’re using an axis-aligned rectangular virtual camera sensor on the Z = -1 plane, we use hFov and vFov field of view angles along the X and Y axis to transform into camera space.

Now, in 3-dimensional camera coordinations, we have the vector containing x_camera, y_camera, and -1, where the *_camera represents the respective point in camera space. From here, we transform the camera space ray to world space using the camera-to-world,c2w_pos, which is a 4x4 homogeneous coordinate system transform matrix given in lecture. We also make sure to normalize the ray’s direction after the transform. Note specifically that since we placed the camera at pos (in world space), we utilize this as column 4.

Then, to set our range for the clipping planes, we utilized nClip and fClip as provided for min_t and max_t as directed.

Task 2: Generation Pixel Samples

Then, now that we’ve generated our camera rays in world space, we now need to generate pixel samples!

We generate ns_aa random samples within the pixel, ensuring to normalize these coordinates. Then, we call camera->generate_ray, passing in these normalized (x, y) coordinates, and then estimate the radiance by calling est_radiance_global_illumination. Finally, once we’ve generated these ns_aa samples, we can average out the pixel color, and then update the sampleBuffer by calling update_pixel with that color.

Task 3: Ray-Triangle Intersection

In order to implement our ray-triangle intersection algorithm, we had to modify the Triangle::has_intersection and Triangle::intersect methods. We created a helper function, Triangle::test, to match the code structure of the ray-sphere intersection, as well as to avoid rewriting code. Our overall strategy was to use the Möller-Trumbore intersection algorithm depicted in lecture.

Our Triangle::has_intersection method allows us to test whether a given ray r intersects a triangle and if so, updates the t, u, and v passed in.

To do so, we first check whether the ray and the plane the triangle lies on is parallel by determining whether the dot product between the two gives 0. If so, we return false (since a parallel ray and plane will never intersect). Then, we calculate t, u, and v. Given the triangle vertices $p_1$, $p_2$, and $p_3$, we calculate

1. the cross product between the direction of the ray and $p_3 - p_1$
2. the cross product between the difference between the origin of the ray and $p_1$, and the difference between $p_2 - p_1$.

From here, we get

• t is the dot product between (2) and $p_3 - p_1$,
• u is the dot product between (1) and the difference between the origin of the ray and $p_1$, and
• v is the division between the dot product of (2) and the direction of the ray by the dot product of (1) and $p_2 - p_1$.

Finally, we check to ensure that

• the intersection point is within the triangle
• 0 $\leq$ min_t $\leq$ t $\leq$ max_t

and then finally update max_t to complete our intersection.

Then, if an intersection occured, we populated isect as desired, filling in the surface normal, primitive, and bsdf.

After completing Task 3, our output for ./pathtracer -r 800 600 -f CBempty.png ../dae/sky/CBempty.dae with normal shading is

Task 4: Ray-Sphere Intersection

Like with ray-triangle intersection, in order to implement our ray-triangle intersection algorithm, we had to modify the Sphere::has_intersection and Sphere::intersect methods (and inadvertently Sphere::test).

Like in our Triangle class, we use Sphere::has_intersection and Sphere::intersect as a way to call into our Sphere::test.

Namely, within the Sphere::test function, we reduce the intersection points to the roots of a quadratic equation. We check whether the discriminant is less than 0, noting that if so, this meant that the ray missed the sphere, and as such, we could immediately return false. We calculate the discriminant by determining

1. the difference between the ray’s origin and the origin of the sphere
2. the dot product between the ray’s direction with itself
3. two times the dot product of the difference between the ray’s origin and the origin of the sphere with the direction of the ray
4. and finally the dot product between the difference between the ray’s origin and the origin of the sphere and itself subtracted by the radius of the sphere, squared

From here, we calculate the discriminant by taking

(3)^2 - 4 * (2) * (4)

If the discriminant is negative, we know that the ray and the sphere do not intersect, and we can return false. If the discriminant is non-negative, we take the square root of the discriminant and select the closest t that still lies between min_t and max_t, using the quadratic formula to determine the candidate’s time of intersection.

After calculating and ensuring that the minimally valid t candidate was selected (by ensuring that only intersections of 0 $\leq$ min_t $\leq$ t $\leq$ max_t are valid), we update max_t like desired to be the chosen candidate.

Then, if an intersection occured, we populated isect as desired, filling in the surface normal, primitive, and bsdf.

After completing Task 4, our output for ./pathtracer -r 800 600 -f CBspheres.png ../dae/sky/CBspheres_lambertian.dae with normal shading is

Here are a few more rendered images with normal shading

Part 2: Bounding Volume Hierarchy

In this part of the homework, we want to implement bounding volume hierarchy (BVH) such that we can speed up the rendering of our path tracer.

We will utilize Task 1, Task 2, and Task 3 to demonstrate our BVH construction algorithm. Within Task 1, we will explain the heuristic we choose for picking the splitting point.

Task 1: Constructing the BVH

In BVHAccel::construct_bvh, we first generate the outermost bounding box. We iterate through the primitives passed in, and expand (union) that bounding box with our current one.

If the number of primitives is less than or equal to the max_leaf_size, we create a new BVHNode, initializing its start and end to be the start and end of the primitives and return.

If not, we need to recurse. In doing so, we will need to first determine the split point and axis selection to split up our bounding volume hierarchy spatially. In reading this paper, we decided to first compute the average centroid across all the primitives. Then, we take the extent of the resulting bbox, which is a vector from opposing vertices of the box (creating a diagonal across the interior of the box). We determine which axis has the longest extent component and use that as the axis to split on.

Then, using the average centroid heuristic, we take the x, y, or z component of our avg_centroid depending on the split axis, and use that as our split point. Then, we call std:stable_partition on the primitives given in, checking whether that primitive’s centroid value at the chosen axis, and checking whether that is less than the averaged centroid across the entire bounding box (for this level) at the chosen axis. If it was less than, we partitioned it in the first set and if greater than or equal two, in the second partition.

Finally, we checked that this start of the second partition was not equal to the original start or end (so one of the partitions wasn’t empty), and if so, we assigned node->l and node->r to be recursive calls to construct_bvh given our two partitions so our l attribute was the first partition and the r attribute was the second partition.

We can see running ./pathtracer -t 8 -r 800 600 ../dae/meshedit/cow.dae and then clicking the right arrow consecutively is

You can see that as we descend into right children in the bounding volume hierarchy, we are always splitting on the longest axis of the extent. For instance, for the first descent, we take the length of the cow, since that is the longest axis. Then, on the second descent, we take the height of the cow, since at that point, the height axis is longer than both the length and width axes. Finally, on the third descent, we split on the width axis of the cow, for the same reason.

Task 2: Intersecting the Bounding Box

Now that we’ve constructed our bounding volume heirarchy, we want to check whether a ray intersects a giving bounding box. Implementing BBox::intersect, we utilize the given ray and axis-aligned plane intersection and ray and axis-aligned box intersection equations given in lecture.

Namely, since we want to represent time as $t = \frac{p_x’ - o_x}{d_x}$, calculating perpendicular to x-axis, we calculate a min_t and a max_t. Then from here, we ensure that for all axes, min_t[axis] < max_t[axis], swapping if that was not the case.

Finally, we calculate the interval of intersection as the maximum of all the min_t axes and the minimum of all the max_t axes since we want to create the tighest bound. If the maximum min_t is greater than the minimum max_t, then we return false (there is no intersection since the ray missed the box), otherwise we set t_0 and t_1 accordingly and return true.

Task 3: Intersecting the BVH

Finally, now that we have the ability to intersect a ray with a bounding box, we can now test that intersection with the BVH. Namely, here, we’ve completed BVHAccel::has_intersection and BVHAccel:intersect. We implement this using this lecture slide on BVH Recursive Traversal as a guide.

For BVHAccel::has_intersection, we first check whether the ray intersects the outermost bounding box and if not, we know it won’t intersect with any primitive, so we can immediately return false. Then, we check that the intersection occurs such that 0 $\leq$ min_t $\leq$ t $\leq$ max_t, and if not, also return false. Finally, we perform a check if the current node is a leaf.

• If so, then we iterate through all the primitives of the node, incrementing the total isects, and then return true if for any of the primitives, it does intersect with the ray.
• Otherwise, we recursively check (by calling BVHAccel::has_intersection on the same ray and check whether it intersects with node->l or node->r).

BVHAccel::intersect functions similarly, except instead of the last step (when we return whether there is an intersection between the ray and node->l or node->r), we create two new Intersections, calling BVHAccel::intersect recursively on node->l and node->r, passing in these two new Intersection objects. Then, we determine which is closer of the left hit or the right hit using boolean algebra and setting the parameter Intersection i to whichever of the hits was closer. Then, like previously, we return whether there was a recursive hit, bubbling up an intersection from the leaves.

BVH Acceleration Analysis

Following the completion of Part 2, here are how multiple different .dae files render. Without BVH acceleration, these files took much longer to render (as we’ll discuss in the analysis below).

The following data was collected by calling ./pathtracer -t 8 -r 800 600 -f {filename}.png ../dae/{path to file}.dae with and without BVH acceleration.

We can see that from these data points, that the rendering speeds up from 800 to a couple thousand times from rendering without BVH acceleration to rendering with BVH acceleration. When we render based on the starter code, we know that we naively have to test intersection with the ray against all the primitives whereas when we use BVH acceleration, we test intersections with the ray against bounding boxes, which hierarchically, means that at the leaf level, we only test against a max of max_leaf_size primitives per leaf node that intersects with the ray. Furthermore, we see that this drastically increases the number of rays that we can test intersection with (since we no longer have to test for all primitives, as our bounding boxes split primitives in a tree-like structure). By enforcing that only the branches of the tree that intersect with the ray need to be traversed impacts the speed at which we can test rays, which greatly increases when we no longer have to test each ray with each primitive. Furthermore, we can see that we also decrease the total number of intersections per ray via the third table: this also makes sense because we no longer need to test a ray against all the thousands of primitives in the images, but rather a subset that will actually be intersected with by the hierarchy of the BVH. This mimics a logarithmic asymptotic, similar to data structures like binary trees. In sum, rendering with BVH acceleration extensively decreases the rendering time since we now need to perform less ray intersection tests with bounding boxes and primitives, which thus allows us to efficiently test through many more rays per second, speeding up the entire operation.

Part 3: Direct Illumination

We will separately walk through the implementations of uniform hemisphere sampling and importance sampling lights. First, though, we’ll describe how to implement the f and sample_f functions and zero-bounce illumination.

Task 1: Diffuse BSDF

In this task, we needed to implement DiffuseBSDF::f to return the evaluation of the BSDF, or reflectance in the given directions. We originally returned the reflectance of the DiffuseBSDF, but noticed this was too bright. Then, we divided this reflectance by \(2 * \pi\) since the surface area of a unit hemisphere is \(2 * \pi\), however this was too dark. We then chose to just divide reflectance by \(\pi\) and this seemed to match the output as expected (we tested this while working on Task 3). Below we’ve attached the output of running ./pathtracer -t 8 -s 16 -l 8 -m 6 -H -r 480 360 ../dae/sky/CBbunny.dae.

We see that using reflectance \(/ \pi\) gave the brightness that best matched the desired output. Dividing by $\pi$ also matches this lecture slide, which we unfortunately did not discover until after completing this task.

We also needed to implement DiffuseBSDF::sample_f, in which we sample a value based on the given pdf for wi, and then return DiffuseBSDF::f called with the passed in wo and our sampled wi.

Task 2: Zero-bounce Illumination

To implement zero-bounce illumination, we took the given Intersection object, and accessed its bsdf attribute and called get_emission() to return the light that results from no bounces of light (the raw emission).

We then updated est_radiance_global_illumination to call zero_bounce_radiance, gaining the following output as expected for running ./pathtracer -t 8 -s 16 -l 8 -m 6 -H -f CBbunny_16_8.png -r 480 360 ../dae/sky/CBbunny.dae.

Task 3: Direct Lighting with Uniform Hemisphere Sampling

Iterating through a total of num_samples samples, we first calculate the incoming light source value by sampling a w_in vector from the hemisphere via hemisphereSampler and calculating its world-space coordinates. Then, we initialized a sample Ray, setting its min_t to EPS_F.

Then, creating an Intersection, we checked whether the bounding volume hierarchy intersected the sampled Ray by calling the intersect function from Part 2 Task 3.

If there was an intersection from calling intersect, then we calculated the \(f_r\), \(L_i\), \(\cos{\theta_j}\), and \(\text{pdf}\) to compute the reflection equation to get the outgoing lighting from this intersection.

• \(f_r\): We call the function we wrote in Task 1 to calculate the BSDF.
• \(L_i\): The incoming light is given by the emission of the intersection’s BSDF.
• \(\cos{\theta_j}\): The dot product between the surface normal of the intersection and the world-space units of the ray gives the cosine angle between the two unit vectors.
• \(\text{pdf}\): Our pdf for hemisphere sampling is \(1 / (2 * \pi)\) because the surface area for a unit sphere is \(4 \pi * r^2 = 4 \pi * 1^2 = 4 \pi\), and hence, the surface area for a unit hemisphere is \(2 \pi\), so the probability of sampling any point uniformly is \(1 / (2 * \pi)\).

Given these values, for each sample Ray, when an intersection existed, we added \((f_r * L_i * \cos{\theta_j}) / \text{pdf}\) to our L_out value.

Finally, after performing num_sample samples, we normalized the outgoing light, L_out, by dividing by num_samples, returning this final L_out.

Running ./pathtracer -t 8 -s 64 -l 32 -m 6 -H -f {filename}.png -r 480 360 ../dae/sky/{filename}.dae for uniform hemisphere sampling gave these two renders.

Task 4: Direct Lighting by Importance Sampling Lights

For importance sampling, we performed a similar workflow as in Task 3’s hemisphere sampling conceptually by calculating incoming light, and generating a sample ray to intersect with, but the structure of the looping was difference since we now were iterating through all the lights.

Particularly, iterating through all of scene->lights, we first checked whether is_delta_light() was true since if we were sampling a point light, it wouldn’t matter what direction the ray was, the outgoing light would still be the same (so we would only need 1 sample, as opposed to num_samples, which we defined similarly to Task 3).

Then, using however many samples we needed for that light as an index, we would iterate that many times, first calling the light’s sample_L, which returns the \(L_i\) for this sample. Calling sample_L would populate

• w_in in world-space, so we would use the w2o matrix to also find the incoming light in local coordinates
• distToLight for the sample
• pdf that will be used to calculate outgoing light

We then utilized the resulting world-space incoming light vector and the hit_p to generate a sample Ray. We would set this Ray’s min_t to EPS_F and its max_t to distToLight - EPS_F.

Then, creating an Intersection, we checked whether the bounding volume hierarchy intersected the sampled Ray by calling the intersect function from Part 2 Task 3.

If there wasn’t an intersection from calling intersect, then we calculated the \(f_r\) and \(\text{pdf}\) (since the \(L_i\) was returned from sample_L and the \(\text{pdf}\) was populated from sample_L) to compute the reflection equation to get the outgoing lighting from this intersection. We note this difference from Task 3 since we are only calculating if there wasn’t an intersection since this means that the light source was not obstructed and thus, would actually hit this hit point. If there was an intersection, then we wouldn’t actually want to illuminate this hit point because another object intersected with it first.

• \(f_r\): We call the function we wrote in Task 1 to calculate the BSDF.
• \(\cos{\theta_j}\): The dot product between the surface normal of the intersection and the world-space units of the ray gives the cosine angle between the two unit vectors.

Given these values, for each sample Ray, when an intersection existed, we added \((f_r * L_i * \cos{\theta_j}) / \text{pdf}\) to our L_out value.

Finally, after performing num_sample samples per light in scene->lights (or 1 if the light was a point light), we normalized the outgoing light, L_out, by dividing by num_samples, adding this outgoing light through this light to our desired outgoing light L_out. After aggregating across all the lights, we would return this final L_out.

While working on this task, we ran into an interesting debugging problem where the shadows seemingly rendered lighter than the areas where light should have hit, like shown in the bunny below.

We realized this was an error as we had kept our original code from Task 3 to add to outgoing light if there was an intersection between the ray and the light–however, we soon realized we should only add to outgoing light if it were the opposite since this would mean there wouldn’t be any obstruction from the light to the desired hit_point.

Running ./pathtracer -t 8 -s 64 -l 32 -m 6 -f {filename}.png -r 480 360 ../dae/sky/{filename}.dae for importance sampling lights gave these three renders.

Now, using ../dae/sky/CBbunny.dae, we can compare the noise levels in soft shadows when rendering with 1, 4, 16, and 64 light rays, and with 1 sample per pixel by running the command ./pathtracer -t 8 -s 1 -l {num_rays} -m 6 -r 480 360 ../dae/sky/CBbunny.dae.

We can see that as the number of light rays increases, there’s less noise in our outputted render. Particularly, the shadows of the bunny becomes more smoothed out, and the edges of shadows become softer (less harsh). The reasoning for this is because when there are less light rays, each shadow point calculated is clearer–each pixelated. However, when we add more light rays, there’s more of a range (similar to when we supersampled in Homework 1). This allows for an overall smoother shadow where there are varying levels of grey that mark out the general portion of the shadow (darker) and the edges (lighter).

Uniform Hemisphere Sampling v. Lighting Sampling

Overall, when we compare uniform hemisphere sampling to lighting sampling, we see that the lighting sampling outputs are much smoother and sharper as opposed to uniform hemisphere sampling, where there was a lot more graininess and noise. Notably, if we look at the walls, we see that the light that hits the walls isn’t smoothly dispersed, but much more grainy in uniform hemisphere sampling. Since in uniform hemisphere sampling, we’re taking samples in different directions around the given point, which means that in parts of an image that don’t converge as quickly, they appear darker since they weren’t sampled as much as needed, resulting in a grainy texture. In contrast, the lighting sampling we used placed importance of samples that actually impact the final lighting: we choose to include outgoing light only when it would be unobstructed and we sample directly from the light source rather than just points within the hemisphere itself. Therefore, since not all the samples of hemisphere sampling were not in the direction of the light source, there was a lot of noice as compared to lighting sampling, where we only get light and shadow rays since we only sample from lights, effectively removing all the noise.

Part 4: Global Illumination

Task 1: Sampling with Diffuse BSDF

This task was a repeat of Task 1 from Part 3.

Task 2: Global Illumination with up to N Bounces of Light

Let’s describe our implementation of the indirect lighting function. We first update est_radiance_global_illumination to first calculate L_out as the sum of the zero_bounce_radiance and at_least_one_bounce_radiance. We also modify raytrace_pixel to mark our camera ray’s depth to be the max_ray_depth.

Within at_least_one_bounce_radiance is where the bulk of our remaining code resides. If the ray’s depth is greater than 0, then we accumulate into L_out the return value of one_bounce_radiance. If the ray’s depth is equal to 1, then this is the L_out that we actually return (if the ray’s depth is less than 1, we return Vector3D(0, 0, 0)).

Then, if the ray’s depth is greater than 1 (so we need to make a recursive call), we first take one random sample_f of the intersection’s BSDF to get the incoming light w_in and calculate its world-space value.

Then, we generate a sample Ray from hit_p and the world-space w_in, setting the Ray’s min_t to EPS_F and max_t to r.depth - 1. This reflects the fact that we previously set the depth of the ray to the number of bounces specified on the command line, and here, our recursive call is decrementing the depth for each bounce we perform.

Then, we check whether nodes in our BVH intersect the sample ray, and if so, we calculate

• \(L_i\): The recursive call to at_least_one_bounce_radiance, which will return to us the radiance from all later bounces.
• \(\cos{\theta_j}\): The dot product between the surface normal of the intersection and the world-space units of the ray gives the cosine angle between the two unit vectors.

Then, we accumulate into L_out the calculated \((f_r * L_i * \cos{\theta_j}) / \text{pdf}\), returning L_out at the end.

Using global illumination (now a combination of direct and indirect illumination), we render the following images using 1024 sampels per pixel:

Direct v. Indirect Illumination

For both ../dae/sky/CBspheres_lambertian.dae and ../dae/sky/CBbunny.dae, we’ve rendered direct illumination (0 and 1 bounce) versus indirect illumination (> 1 bounce) at 1024 samples per pixel. Note in particular that the 0 and 1 bounce illumination resembles what we did in Part 3 for zero- and one-bounce illumination while indirect illumination is what we accounted for in Part 4. Interestingly, we can see light now hitting the ceiling in the indirect illumination only renders as the light has come out of the area light, bounced off to some other surface, and then bounced to illuminate the ceiling.

These were probably the prettiest and most interesting images to render: in the indirect images, we can see that there’s somewhat of a soft glow around the sphere and bunny. It’s quite satisfiying!

Adding Non-Accumulation Rendering

In order to account for rendering non-accumulated bounces, we needed to utilize the isAccumBounces variable of the PathTracer class. Namely, we know that when isAccumBounces is false, we only want to add the outgoing light for the last bounce, and for all other bounces, return Vector3D(0, 0, 0).

To facilitate this, we updated two things, namely just instances where we add to L_out.

• Instead of always adding one_bounce_radiance to L_out, we only call this when isAccumBounces is true or the depth of the ray is one (which means we’ve hit the last bounce).
• Instead of accumulating to L_out (summing), we directly set L_out equal to \((f_r * L_i * \cos{\theta_j}) / \text{pdf}\) if isAccumBounces is false since we don’t want to accumulate the outgoing light, but instead, just return \((f_r * L_i * \cos{\theta_j}) / \text{pdf}\), noting that \(L_i\) is the output of our recursive call to at_least_one_bounce_radiance.

Here, in two columns, we shown the output of ../dae/sky/CBbunny.dae sampled 1024 times per pixel, from the 0th to 5th bounce of light, where the left column accumulates the light (so isAccumBouces is true) while the right column does not accumulate the outgoing light and instead just shows the light from that specific bounce (isAccumBounces is false).

Between the second and third bounce of light, we see that the image gets a lot darker overall (the bunny texture isn’t as clear around its legs), however, it looks like the edges of the floor still retain some light. Overall, we can understand why because after the second bounce, it makes sense that the light has already dissapated. The decrease in intensity of outgoing light experienced here could also be due to the fact that the light could have also been absorbed by the bunny or the walls. In the second bounce, light that hit the ground would’ve hit the bottom side of the bunny, explaining why its legs were most likely illuminated. By the third bounce, the only really light areas with higher density of sharpened light seems to be near the edges of the floor in the back and a little halo under the bunny.

In contrast to purely rasterized images, we know that adding bounces and global illumination definitely impacts the realism of the image since we’re now able to add complex shadows and diffuse shadows and colors more smoothly across images. Pure rasterization did not allow for the interaction of light within a scene, which meant that prior to our work here, in previous projects, we were unable to render anything out of a controlled setting–it was always just one item in empty space. Now, we have the ability to utilize lights in our staging and mimic shadows and perspective similarly how we perceive three-dimensional objects in real life. Pretty neat!

Pixel Sampling Rate Variation

Next, we render ../dae/sky/CBbunny.dae at 4 light rays with various sample-per-pixel rates.

We see that as the pixel sampling rate increases, the bunny gets smoother, clearer, and less grainier: the perks of higher sampling rates!

Task 3: Global Illumination with Russian Roulette

The intention behind Russian Roulette is to allow for unbiased random termination in our recursive model for modeling outgoing light at different bounces. In an ideal world, we would be able to continuously recurse for high depth renders, however, Russian Roulette allows you to probabilistically choose, at each step, whether to stop the recursion or not. This makes previously computationally infeasible renders that would take (down) horrendous amounts of time now feasible.

To implement Russian Roulette random termination in at_least_one_bounce_radiance, we decided on a termination probability of 0.3. Then, when we need to make a recursive call (so when r.depth > 1), then we take a coin_flip, passing in the probability of 0.7 since we want a termination probability of 0.3. Then, only if the coin_flip returns true or if the depth of the ray is equal to the maximum ray depth, then we will recurse. We have the second condition since we note that if max_ray_depth > 1, we need to at least trace one indirect ray bounce (hence, we must recurse at least once).

Here we render ../dae/sky/CBbunny.dae with Russian Roulette with 1024 samples per pixel, setting the maximum ray depth set to 0, 1, 2, 3, 4, and 100.

With Russian Roulette work, we note that similar to previous renders, the images get brighter (more illuminated) with maximum ray depth. This makes sense with our light accumulation renders above. However, we note that between a maximum ray depth of 4 and a maximum ray depth of 100, there isn’t that big of a difference between the images. We note this probably is due to the fact that with Russian Roulette random termination, the image rendered with a maximum ray depth of 100 most likely did not reach up to 100 levels of bounces, with a termination probability of 0.3 at each level of recursion.

Part 5: Adaptive Sampling

Path tracing via Monte Carlo estimation is powerful, but now, we want to solve the problem of having noise in our images: thus, we now move towards adaptive sampling, which takes a high number of samples per pixel, but concentrates samples in areas that are more difficult to render, namely, reducing the number of samples for pixels that converge quickly.

Our implementation of adaptive sampling builds off of the raytrace_pixel implementation from Part 1 Task 2. We define num_samples as ns_aa and two doubles, s1 and s2, to keep track of the sum of illuminances and the sum of squared illuminances respectively.

Iterating through num_samples times, we first check whether our index i is now a multiple of samplesPerBatch, and if so, we want to perform a convergence check.

Our convergence check composes of defining a mean mu, equal to s1 / i. We also defined the variance sigma_squared, equal to (s2 - (s1 * s1) / i) / (i - 1). Given this, we want to calculate \(I = 1.96 \cdot \frac{\sigma}{\sqrt{n}}\) to measure the pixel’s convergence. We calculate I = 1.96 * sqrt(sigma_squared / i), where i is the number of samples so far. Knowing our I, we know that if \(I \leq maxTolerance \cdot \mu\), then we can assume the pixel has converged, so we can stop tracing more rays. Thus, we set the num_samples equal to our current index i and break out of our loop.

However, if the pixel has not yet converged, then we get a sample from the grid and normalize the sample within the pixel. From there, we generate a Ray and set a sample_radiance by calling est_radiance_global_illumination, adding that to our overal pixel’s value (radiance). Then, we sample the illuminance from sample_radiance, adding that value to s1 and adding its square to s2.

Outside, once we’ve completed the iteration, we normalize the pixel’s average radiance by dividing by the number of samples we actually took, saving this value into our sampleCountBuffer, and finally call update_pixel, unmodified from our original code from Part 1 Task 2.

Running ./pathtracer -t 8 -s 2048 -a 64 0.05 -l 1 -m 5 -r 480 360 -f {filename}.png ../dae/sky/{filename}.dae shows the impact of adaptive sampling. We show on the left, the actual image rendered alongside the image on the right that shows the sample rate of every pixel. We note that red and blue colors represent high and low sampling rates respectively, where sample rates are the ratio between the actual number of samples taken (so our ending num_samples, depending on whether we converged early) and the maximum number of samples allowed (here, 2048).

For instance, if we take a look at the bunny, we see that the area around the bunny’s lower body and ear have a higher sampling rate. This makes sense with the final rendered result since that area is rather bumpy with texture.

Contributors

Ashley Chiu, Edward Park

This project was difficult to execute and coordinate with the conceptual complexity and large codebase to work with. However, we were able to divide out the work well and debug together when needed. We mostly worked asynchronously and kept each other up on what progress we’d made, however, we did work together in person for Part 2 Task 1 and to debug Part 3 Task 3. In general, this collaboration worked for us since we were able to work when we had the time, but still discuss problems that were difficult to debug. Throughout the project, we not only learned about how to represent light within images, how ray tracing is implemented, and how lights interact with surfaces, but furthermore, how to document our code well and communicate blockers. We also learned which computers had the fastest compute power and were able to render images faster, and which laptops had issues rendering larger images or would run out of memory (we had to deal with memory management to resolve this).