In our last post, we looked at the SSE optimized BVH traversal of cranberray. BVH construction and traversal is a fairly complicated portion of cranberray; cranberray’s sampling strategy however is quite a bit simpler.
Sampling in cranberray is simple. We first split our image into a grid representing every pixel of our image. (Note that a rectangular pixel is a simplifying assumption ) Then we take our pixel box and sample it X number of times using X different sample points and average the results. Our sampling strategy determines how we spread out those X different sample points within our pixel rectangle.
The simplest strategy is simply to select a point randomly within our pixel box for every sample. However, this has the undesirable property of potentially introducing clumping causing us to miss some desirable information. 
Instead, cranberray makes use of a simple technique called N-Rooks sampling.
N-Rooks is what’s known as a low discrepancy sequence. It has the desirable property of avoiding the clumping of samples that random sampling has and avoiding the aliasing issues of regular sampling. We won’t be explaining low discrepancy sequences here, but I refer the interested reader to  and  for excellent treatments on the subject.
N-Rooks sampling is very simple.
You first generate jittered samples along the diagonal of an X by X box.
And then you shuffle the columns of said box by starting from the leftmost column and randomly swapping with the columns to the right of our current column. (Including our current column). This is known as the Fisher-Yates shuffle.  It has the desirable property of uniformly selecting from any of the X! permutations of our columns.
And that’s it! We now sample our scene using these sample points and average their results to get our final output.
In our next post, we’ll be looking at BRDF multiple importance sampling using the balance heuristic as well as the GGX-Smith microfacet specular BRDF. Until then, happy coding!
While writing this, I was considering the possibility of selecting the value of our pixel using a different strategy than the mean filter used with Monte Carlo integration. I would have to do some further digging to see if this has been presented anywhere. Something along the lines of a Gaussian filter or potentially something else entirely.
I would also like to implement other low discrepancy sequences in cranberray if ever time permits. N-Rooks was very simple, it would be interesting to be able to select between different sampling strategies.
In our last post, we looked at our SIMD-friendly math library. In this post, we’ll be taking a look at how cranberray builds its BVH. For a primer to BVHs, I recommend going to PBRT’s treatment on the topic here.
Canberray’s BVH is built recursively as a binary tree. It starts with the top level node containing all nodes and recursively splits it into 2 children at each level. Originally, the BVH was split simply by selecting the left and right halves of the children nodes. The next approach was to split according to the median of a particular axis (typically the largest one), and now it makes use of the surface area heuristic.
The Surface Area Heuristic or SAH is a well known strategy to build an efficient BVH that tries to minimize the number of intersection tests applied on a tree.
PBRT has an excellent treatment on the topic , and I’ve linked one of the original papers referencing the idea . I will try to instead provide an intuitive picture of why and how SAH works.
The idea behind SAH is to discover an efficient partition of a collection that will try to minimize the cost of traversing the tree for any ray.
The SAH function is presented as:
This cost model is based on probabilities. is the cost of traversing the parent node to determine if any of the two children nodes intersect the ray. In cranberray, this would be the cost of testing the ray against both children AABBs and determining if it intersects either. is the cost of testing the shape itself against the ray. In cranberray, this would be the cost of testing the triangles against the ray. Finally, and are the probabilities that the bounding volumes containing the children nodes intersect the ray.
Now, well make a slight modification and reformulate and as and where C is our parent node. This states that represents the probability of the event A happening given the event C has already happened. In this context, this means that represents the probability that we intersect our volume A given that we’ve already intersected our parent volume C.
With this change, our function is now:
Now, let’s imagine that we’re a lonesome ray traveling through our BVH. We arrive at a node and we test both of it’s children to see if we intersect it. This takes us 10 units of time (). At this point we would intersect A with a probability of and we would intersect B with a probability of . In this example, we only intersect with A, and so we test our ray against it’s contained shapes. Let’s say that A has 10 shapes, as a result, we have to pay the cost of testing the intersection with all 10 shapes. This is . Finally, let’s say that the cost of is 15 for every child shape.
This gives us the final actual cost of:
What our original cost function is doing, is calculating the expected cost of doing this operation for any possible ray.
As a result, if we were to build with these values:
Then our expected cost for this node would be:
Which is exactly the cost of our example above because if our A node has a probability of 1, we would always experience the scenario above.
Of course, if we were to take some more likely value such as:
Which would be the expected value of our node. Notice that our probabilities can add up to more than or less than 1 because it is reasonable that we might be able to intersect with both volumes or none at all.
Now we need to determine how we want to define . This probability is the probability of a random ray intersecting with our volume given that we’ve intersected with our parent volume.
is given as the surface area of our volume A divided by the surface area of our parent volume C. As a result, we can define as the probability of intersecting A given that we’ve intersected C.
Derivation In 2D
In 2 dimensions, this makes sense, as we can image a rectangle, embedded in another rectangle. The probability of a random point being placed in the child rectangle given that it is contained in the parent rectangle is equal to the area of the embedded rectangle divided by the area of the parent rectangle.
Now, instead we’ll imagine uniformally random ray instead of a random line. To select a random ray, we can select a random direction and then select a random ray along the plane defined by that direction.
With this in mind, we want to calculate the probability that our ray will intersect the children volume given that it has intersected the parent volume.
To do this, we want to calculate the average projected area of our shape.
Imagine that we’ve selected our random direction for our ray, we can see that our ray will have intersected our shape if it’s projection onto the plane intersects our random ray.
Now if we add a second volume embedded in the first volume, we can see that the probability of intersecting the first volume would be given by the projected area of the child volume divided by the projected area of the parent volume.
This example only works for a single direction. However, since we’re looking for a probability given any ray direction, we need to calculate the average of this projection across the sphere of all directions.
In 2 dimensions, the projected area of a circle is always it’s diameter. We will use this fact to test our work.
To calculate the average projected area of an arbitrary convex shape, we first want to split that shape into a series of small discrete patches.
This will be our integration domain.
Once we’ve split our shape into discrete patches, we want to calculate the average projected area of a single patch.
Now we’ll integrate across the circle of directions containing the patch and divide it by the area of our integration. This will effectively give us the average projected area for a single patch of our shape. This works because we will have calculated the projection for every single direction.
Integrating with an absolute value is tricky. Instead due to the symmetry of the problem, we can look at only a quarter of our circle and multiply by 4 and remove our absolute term.
Since our expression reduces to
Finally, we want to divide our value by 2, this is because we’ve now counted our projected area twice.
Using this formulation, we can calculate the average projected area of our shape by adding all of our average projected areas together.
And there you have it, the average projected area of a 2 dimensional convex object is the area of said object.
Notice that if we apply this to the surface area of a circle with the formula we get which is our diameter which matches up with our expectation.
Derivation In 3D
We can take the same steps in 3D.
Calculate the average projection over our sphere of directions:
Integrating across the positive hemisphere and multiplying by 2
Finally, dividing by 2 for our double projection
And plugging into our surface area calculation
Putting It Together
Finally we can see that our average projected area in 3D is it’s surface area.
To calculate our probability, we simply want to divide our parent’s projected area (C) divided by our child’s projected area (A)
Where A and C is the surface area of our volumes. And voila, that’s how we got that original formula.
In Depth View
Now that we have all the pieces, we can take a look at the construction of our BVH.
Cranberray builds it’s BVH from the top down. Starting from a containing bounding volume and splitting it in 2 recursively.
Cranberray keeps a queue of bounding volumes to process as a ring buffer. This makes management of front popping and back pushing very simple but makes resizing the queue trickier. As a result, Cranberray simply allocates a large buffer. We could likely allocate the maximum possible number of elements in the ring buffer instead (something along the lines of where n is the next largest power of 2)
Cranberray then selects the axis with the widest breadth of centroids. The code for this looks like so:
The final special format of our BVH is that all the leaf nodes in our tree are stored at the end of our array. This allows us to test if a node is a leaf node by simply comparing the index in our array against the size of our array minus the number of leaves contained in the tree. This allows us to use data that we’ve already loaded into memory instead of requiring us to load extra data to use in our branch introducing a data dependency.
I’ve recently been diving into the world of importance sampling and I decided to share my derivation for importance sampling the cosine lobe.
When we’re shading a point in path tracing, we typically shoot a ray from our surface in a uniformity random direction contained on a hemisphere centered about our normal. This has the downside of introducing quite a bit of variance into our renders.
Imagine that we have a very bright light that only occupies a very small projected area from our shaded point.
We would be very likely to miss this light with most of our samples and our render could turn out much darker than we would expect.
This is where importance sampling comes in.
If you imagine that your illumination is a polar function
If we were to sample a random variable with a distribution that matches this function, we would be much more likely to hit the important points of our function. (Hence importance sampling)
I won’t dive deeply into this topic, as there are a variety of excellent resources detailing this topic. 
The essence of it however, is that you want to find a Probability Density Function (PDF) that matches the shape of your illumination function. Once you’ve define this PDF, you can sample it using the Cumulative Density Function (CDF).
Since our cosine lobe illumination will look like this:
We will use this as the basis to derive our distribution since we’re most likely to get the most light from directions arriving parallel to our normal.
Thankfully, our cosine lobe has an analytical formula that we can use as our PDF.
Our PDF must integrate to 1, we integrate the PDF across our hemisphere
Plug in (1)
Plug in (2)
Since our PDF has to integrate to 1,
Plug in (3),
Finally, plug in (4) into our PDF,
Now that we have our PDF, we can calculate our PDF in terms of and .
Now we integrate with respect to to get
And then to get ,
Now we want to calculate the CDF of each function,
Now we want to invert our CDF to sample it using our random variable ,
Now we have our CDFs and our PDFs, we can finally calculate our direction.
In pseudo code you can simply do:
With these directions, you can now sample your scene:
Plug in (5)
That’s it! These formulas will sample the hemisphere where it is receiving more light defined by the cosine lobe. The results are pretty awesome.
Recently, I’ve been wrestling with attempting to define the mathematics behind the illumination of a triangular area light. Although I have not succeeded yet, I figured I would chronicle my recent attempts. Here goes!