Rendering Constructive Solid Geometry With Python

What It Means for Ray Tracing

Ray tracing with CSG's is surprisingly straightforward. When a ray is interesected with a CSG object, it is really intersected with the set of surfaces that make up the object. The job of the CSG algorithm is to take those two sets of intersections and filter them out so that only valid ones remain. These hits are then returned to the renderer so it can determine the closest surface that the ray interacts with.

An example of rendering the union of two surfaces. While each sub-shape returns a full set of intersections, the filter function discards points a1 and b0 since they are inside of the new object.

Writing a function to sort valid CSG hits turns out to be nontrivial. For every ray-surface hit, the function needs to determine if the hit occurred inside the opposite surface. This is easy enough for convex surfaces that have at most two intersections per ray, but as the number of intersections grows, basic methods for filtering fall apart or slow down significantly.

Complex hit arrays for two surfaces A & B, and the expected hit arrays for CSG Operations on those surfaces

The surface_count Vector

There are two additional arrays we'll use validate hits: (1) a sorted array of all hits for both surfaces, and (2) an array that tracks how many surfaces the ray is inside of at each hit, called surface_count.

The same intersection arrays from above and the corresponding surface_count array

To create surface_count we'll do the following:

1. Create a concatenated array of both hit arrays (but do not sort it yet).

2. Take the argsort of the concatenated array.

3. Create a new array with the same dimension as the concatenated array. For each index in the new array, assign +1 if the value in the equivalent argsort index is even, and -1 if is odd.

4. Take the cumulative sum of the +/-1 array, this is the surface count array.

5. Sort the concatenated hit array.

The +/-1 array is used to indicate if a hit is entering or exiting one of the two surfaces. Since both arrays have an even number of elements, every odd value of the concatenated array must enter a surface, and every even hit must exit. NumPy's where, combined with argsort and cumsum are all we need to create surface_count. Once we have that we're ready to tackle the first of our Boolean operators.

merged_array = np.concatenate((array1, array2))
merged_argsort = np.argsort(merged_array, axis=0)
merged_array = merged_array[merged_argsort]
merged_mask = np.where(merged_argsort & 1, -1, 1)
surface_count = np.cumsum(merged_mask, axis=0)

Union

Take a look at the union operator and corresponding count values. Notice that union hits are any any index (n) where count[n-1]==0 and count[n]==1, or count[n-1]==1 and count[n]==0. This is the same as taking the exclusive or (XOR) of the array with a copy of itself shifted down by one row. From above we know that the last value of the count array has to be 0 (at \(t=\infty\) the ray must have exited all surfaces), so the zeroth row of the shifted array will always be all zeros. NumPy gives us all the function calls needed to efficiently perform xor on our count array, shown below.

surface_count = np.logical_xor(surface_count, np.roll(surface_count,1,axis=0))
csg_hits = np.where(surface_count != 0, merged_array, np.inf)

A couple things from this code might pop out: (1) why specify the axis when rolling a 1D array, and (2) why does csg_hits need the same dimension as hit_array, padded with np.inf, and not just return the valid hits. Both of these are addressed in Extending To 2D Matrices.

Intersection

The intersection operator can be handled in a similar manner. Looking at the count array, an intersection hit occurs at any index (n) where count[n] == 2 or count[n-1] == 2. This time we'll use NumPy's logical_or function to create the mask.

is_two = (surface_count == 2)
mask = np.logical_or(is_two, np.roll(is_two, 1, axis=0))
csg_hits = np.where(mask, merged_array, np.inf)

The intersection operator has an interesting "blip" at t=5. This is because the ray intersects both surfaces at that time, but enters the first while exiting the second. Using integer math this becomes a 'zero thickness' shell, but it can cause unintended results when the hits are floating points and the surfaces overlap by a small amount.

Difference

The difference operator is unique from the union and intersection operator. Both union and intersection have analogs in boolean algebra, (or the binary equivalent and and or operators). Boolean algebra, however, does not have the concept of subtraction. Not only that, but we can't use the clever tricks from the count array, since we need to know which surface is the subtracting surface, and the count array only tracks how many closed surfaces the ray is inside. Rather than creating an entirely new method just for differences, we're going to redefine what a difference means so it it behaves like a boolean operation.

Don't think of A-B as shape B cutting away from shape A, think of it as the intersection of A with the infinitely large volume of space where B does not exist, called \(\bar{B}\).

Defining the function in this way lets us reuse the the same principle as the intersection operator, but first the count array has to be redefined for an inverted shape. An inverted shape still has to follow the assumptions from above, but the ray enters the shape at \(-\infty\) and exits the shape at the first hit in the hit array. Similarly, the ray enters the shape at the last hit of the hit array, and exits at \(\infty\).

Padding the intersection array with \(+/-\infty\) is unnecessary (and NumPy arrays don't like to be resized). Instead, the following observations allow us to work with the unmodified array.

• The hit at \(-\infty\) will always be the first hit in the sorted hits array, meaning the first value in the cumulative sum will always be a 1. This is the same as adding 1 to the count array and ignoring the hit at \(-\infty\).

• The hit at \(+\infty\) will always be the last hit in the sorted hits array. If we ignore it, the cumulative sum's final value will be 1 (since the ray is still inside of the inverted surface). However, looking at the Intersection operator, we only care about finding indices where count[n]==2, so we're safe to ignore it.

• The pre-summed count array needs to be filled with +/-1 from 0:n, where n is the length of the first hit array, and -/+1 from n:-1. This will successfully "invert" the second shape, where hits that used to enter the surface now exit, and visa versa.

With those observations in hand we're ready to create the count array for the difference operator, and once again NumPy's logical_xor will simplify the task.

merged_array = np.hstack((left_array, right_array))
merged_argsort = np.argsort(merged_array,axis=0)

count_array = np.where(np.logical_xor(merged_argsort&1, merged_argsort>=left_array.shape[-1]))
count_array = np.cumsum(count_array)+1

Here the xor operator is inverting the array mask for any indices in the argsort that reference the right_array.

Extending To 2D Matrices

PyRayT can perform reasonably fast ray tracing because under the hood every ray is stored in a 2x4xn matrix that gets intersected with each surface. This allows me to bypass Python for loops in favor of heavily optimized NumPy functions. Fortunately, by specifying axes and preserving array sizes, the csg function can be readily extended to 2D matrices, where every column represents the ordered hits for an individual ray with the given surface.

The only thing we need to change is how the arrays are concatenated. If a 1D array is passed, they can be concatenated along the zero axis, but 2D arrays need to be stacked column-wise.

if array1.ndim == 1:
    # if 1D arrays were passed, concatenate
    merged_array = np.concatenate((array1, array2))
    merged_argsort = np.argsort(merged_array, axis=0)
    merged_array = merged_array[merged_argsort]

else:
    # otherwise stack them where each column represents a unique ray's hits
    merged_array = np.vstack((array1, array2))
    merged_argsort = np.argsort(merged_array, axis=0)
    merged_array = merged_array[merged_argsort, np.arange(merged_array.shape[-1])]

Verifying Test Cases

Before passing this function off as complete, we need to make sure it passes some basic unit tests. I'll be using Python's UnitTest framework to verify that the two hit arrays plotted above return the correct values for union, intersection, and difference.

import unittest
import numpy as np
import csg

class TestArrayCSGOperation(unittest.TestCase):
    def setUp(self) -> None:
        self.array1 = np.array((1, 4, 5, 10))
        self.array2 = np.array((0, 2, 3, 5, 6, 7, 8, 9, 11, 12))
        self.expected = np.full(self.array1.shape[0] + self.array2.shape[0], np.inf)

    def test_add_operation(self):
        unioned = csg.array_csg(self.array1, self.array2, csg.Operation.UNION)
        self.expected[:4] = (0, 10, 11, 12)
        self.assertTrue(np.allclose(unioned, self.expected), f"expected {self.expected} but got {unioned}")

    def test_intersection_operation(self):
        intersected = csg.array_csg(self.array1, self.array2, csg.Operation.INTERSECT)
        self.expected[:10] = (1, 2, 3, 4, 5, 5, 6, 7, 8, 9)
        self.assertTrue(np.allclose(intersected, self.expected), f"expected {self.expected} but got {intersected}")

    def test_diff_operation(self):
        diffed = csg.array_csg(self.array1, self.array2, csg.Operation.DIFFERENCE)
        self.expected[:8] = (2, 3, 5, 6, 7, 8, 9, 10)
        self.assertTrue(np.allclose(diffed, self.expected), f"expected {self.expected} but got {diffed}")


if __name__ == '__main__':
    unittest.main()

With those tests passing, we're one step closer to a fully fledged Python ray tracer!

Constructive Solid Geometry can extend the functionality of a ray tracer by building complex shapes from basic primitives, but rendering them requires an additional step to filter out which intersections are valid. Thanks to NumPy, it's easy to write this function without reverting to Python for loops, and the same function can be used to process multiple ray-surface intersections at once.

In PyRayT, all lenses are CSG intersections of two spheres (defining the focus) and a cylinder that sets the aperture. This has sped up development time and is significantly easier than writing custom functions for each optical component.