Camera rays, intersections, and the first BVH

This article continues through the first commit. It follows the conversion from a pixel sample to a world-space ray, then concentrates on the choices made by HdRestir’s first CPU geometry path. Standard intersection algorithms are linked rather than derived again.

Pixel samples and camera rays

For a pixel sample \((p_x,p_y)\) in an image of width \(W\) and height \(H\), the renderer first maps image coordinates to normalized device coordinates:

The inverse projection maps the point from normalized device coordinates into camera space. The inverse view transform then maps both camera origin and the near-plane point into world space. The first-commit implementation is:

CameraRay GenerateCameraRay(
    const GfMatrix4d& inverseViewMatrix,
    const GfMatrix4d& inverseProjMatrix,
    float px, float py, int width, int height)
{
    const float ndcX = (2.0f * px / width)  - 1.0f;
    const float ndcY = (2.0f * py / height) - 1.0f;

    const GfVec3f origin{inverseViewMatrix.Transform(GfVec3d(0.0, 0.0, 0.0))};
    const GfVec3f nearCam{inverseProjMatrix.Transform(GfVec3d(ndcX, ndcY, -1.0))};
    const GfVec3f nearWorld{inverseViewMatrix.Transform(GfVec3d(nearCam))};

    return {origin, (nearWorld - origin).GetNormalized()};
}

This is a literal function from camera_ray.cpp. The code assumes that GfMatrix4d::Transform performs the homogeneous operation appropriate for the supplied point. The ray is

with \(o\) equal to origin and \(d\) equal to the normalized difference. Jittering \(p_x\) and \(p_y\) inside the pixel changes the sampled primary ray without changing this mapping.

The overload in the same file adds lens distortion and depth of field. For a circular aperture it samples a disk with

The square root is required for uniform area density: the area inside radius \(r\) is proportional to \(r^2\), so setting \(r^2=U_1\) makes the cumulative probability proportional to area. The sampled lens point is aimed at the configured focal-plane point.

Ray–triangle intersection

The initial BVH uses the algorithm described by Möller and Trumbore in Fast, Minimum Storage Ray–Triangle Intersection. That paper derives the determinant and barycentric-coordinate test. Here the relevant question is how the commit specializes it:

GfVec3f edge1 = v1 - v0;
GfVec3f edge2 = v2 - v0;
GfVec3f pvec = GfCross(rayDir, edge2);
float det = GfDot(edge1, pvec);
if (std::abs(det) < 1e-8) return false;
float invDet = 1.0f / det;
GfVec3f tvec = rayOrigin - v0;
float u = GfDot(tvec, pvec) * invDet;
if (u < 0.0f || u > 1.0f) return false;
GfVec3f qvec = GfCross(tvec, edge1);
float v = GfDot(rayDir, qvec) * invDet;
if (v < 0.0f || u + v > 1.0f) return false;
t = GfDot(edge2, qvec) * invDet;
outU = u;
outV = v;
return (t > 1e-4);

The two thresholds are policies, not parts of the general algorithm. 1e-8 classifies a determinant as numerically zero, while the 1e-4 lower bound rejects a hit close to the ray origin. Both are expressed relative to the floating-point values and scene units used by this implementation. The latter can be too large for a small scene or too small for a large one; it is not a scale-independent solution to self-intersection. Wächter and Binder give a more robust error-bound construction in A Fast and Robust Method for Avoiding Self-Intersection.

AABB slab test

The node test is the standard ray–axis-aligned-box slab test. The PBRT discussion of bounding volume hierarchies provides the broader derivation and explains how the predicate is used during traversal. This commit implements the interval update as follows:

float tmin = -1e30f;
float tmax = 1e30f;

for (int i = 0; i < 3; ++i) {
    if (std::abs(rayDir[i]) < 1e-8) {
        if (rayOrigin[i] < min[i] || rayOrigin[i] > max[i]) return false;
    } else {
        float invDir = 1.0f / rayDir[i];
        float t1 = (min[i] - rayOrigin[i]) * invDir;
        float t2 = (max[i] - rayOrigin[i]) * invDir;
        if (t1 > t2) std::swap(t1, t2);
        tmin = std::max(tmin, t1);
        tmax = std::min(tmax, t2);
        if (tmin > tmax) return false;
    }
}

The implementation-specific branch treats abs(rayDir[i]) < 1e-8 as parallel. It rejects the box when the origin lies outside that slab; otherwise it avoids the unstable reciprocal. As with the triangle test, the fixed threshold is a numerical policy that should be reviewed against the renderer’s coordinate range.

BVH construction in this revision

A BVH and common construction strategies are covered by the same PBRT reference. The useful information here is the exact builder introduced by this commit. Each node owns a contiguous triangle range; four or fewer triangles form a leaf. A larger range splits triangle centroids at the midpoint of the node’s largest axis:

GfVec3f size = node.bounds.GetSize();
int axis = 0;
if (size[1] > size[0]) axis = 1;
if (size[2] > size[axis]) axis = 2;

float splitPos = node.bounds.GetMin()[axis] + size[axis] * 0.5f;

int i = start;
int j = end - 1;
while (i <= j) {
    if (_triangles[i].centroid[axis] < splitPos) {
        i++;
    } else {
        std::swap(_triangles[i], _triangles[j]);
        j--;
    }
}

int leftCount = i - start;
if (leftCount == 0 || leftCount == count) {
    i = start + count / 2;
}

The fallback guarantees non-empty child ranges, but it does not reorder the triangles by centroid when the midpoint partition fails. It therefore guarantees progress and balanced counts, not spatially good child bounds. Traversal can still approach linear work when both child boxes overlap many rays. The commonly quoted \(O(\log N)\) traversal cost is not a worst-case property of this builder.

Reconstructing the hit

At a leaf, the closest accepted triangle updates distance, geometric normal, UV coordinates, smooth normal, tangent frame, vertex color, and material index. Barycentric interpolation uses

The code applies that same expression to UVs, normals, and colors. Normals are renormalized after interpolation; UVs and colors are not. The material index selects the shading model in the next article.