6.7 Curves

While triangles or bilinear patches can be used to represent thin shapes for modeling fine geometry like hair, fur, or fields of grass, it is worthwhile to have a specialized Shape in order to more efficiently render these sorts of objects, since many individual instances of them are often present. The Curve shape, introduced in this section, represents thin geometry modeled with cubic Bézier curves, which are defined by four control points, normal p 0 , normal p 1 , normal p 2 , and normal p 3 . The Bézier spline passes through the first and last control points. Points along it are given by the polynomial

normal p Subscript Baseline left-parenthesis u right-parenthesis equals left-parenthesis 1 minus u right-parenthesis cubed normal p 0 plus 3 left-parenthesis 1 minus u right-parenthesis squared u normal p 1 plus 3 left-parenthesis 1 minus u right-parenthesis u squared normal p 2 plus u cubed normal p 3 period

(See Figure 6.29.) Curves specified using another basis (e.g., Hermite splines or b-splines) must therefore be converted to the Bézier basis to be used with this Shape.

Figure 6.29: A cubic Bézier curve is defined by four control points, normal p Subscript Baseline Subscript i . The curve normal p Subscript Baseline left-parenthesis u right-parenthesis , defined in Equation (6.16), passes through the first and last control points at u equals 0 and u equals 1 , respectively.

The Curve shape is defined by a 1D Bézier curve along with a width that is linearly interpolated from starting and ending widths along its extent. Together, these define a flat 2D surface (Figure 6.30). It is possible to directly intersect rays with this representation without tessellating it, which in turn makes it possible to efficiently render smooth curves without using too much storage.

Figure 6.30: Basic Geometry of the Curve Shape. A 1D Bézier curve is offset by half of the specified width in both the directions orthogonal to the curve at each point along it. The resulting area represents the curve’s surface.

Figure 6.31 shows a bunny model with fur modeled with over one million Curves.

Figure 6.31: Furry Bunny. Bunny model with over one million Curve shapes used to model fur. Here, we have used unrealistically long curves to better show off the Curve’s capabilities, giving an unrealistically shaggy bunny. (Underlying bunny mesh courtesy of the Stanford Computer Graphics Laboratory.)

<<Curve Definition>>= 
class Curve { public: <<Curve Public Methods>> 
static pstd::vector<Shape> Create(const Transform *renderFromObject, const Transform *objectFromRender, bool reverseOrientation, const ParameterDictionary &parameters, const FileLoc *loc, Allocator alloc); PBRT_CPU_GPU Bounds3f Bounds() const; pstd::optional<ShapeIntersection> Intersect(const Ray &ray, Float tMax) const; bool IntersectP(const Ray &ray, Float tMax) const; PBRT_CPU_GPU Float Area() const; PBRT_CPU_GPU pstd::optional<ShapeSample> Sample(Point2f u) const; PBRT_CPU_GPU Float PDF(const Interaction &) const; PBRT_CPU_GPU pstd::optional<ShapeSample> Sample(const ShapeSampleContext &ctx, Point2f u) const; PBRT_CPU_GPU Float PDF(const ShapeSampleContext &ctx, Vector3f wi) const; std::string ToString() const; Curve(const CurveCommon *common, Float uMin, Float uMax) : common(common), uMin(uMin), uMax(uMax) {} DirectionCone NormalBounds() const { return DirectionCone::EntireSphere(); }
private: <<Curve Private Methods>> 
bool IntersectRay(const Ray &r, Float tMax, pstd::optional<ShapeIntersection> *si) const; bool RecursiveIntersect(const Ray &r, Float tMax, pstd::span<const Point3f> cp, const Transform &ObjectFromRay, Float u0, Float u1, int depth, pstd::optional<ShapeIntersection> *si) const;
<<Curve Private Members>> 
const CurveCommon *common; Float uMin, uMax;
};

There are three types of curves that the Curve shape can represent, shown in Figure 6.32.

  • Flat: Curves with this representation are always oriented to face the ray being intersected with them; they are useful for modeling fine swept cylindrical shapes like hair or fur.
  • Cylinder: For curves that span a few pixels on the screen (like spaghetti seen from not too far away), the Curve shape can compute a shading normal that makes the curve appear to actually be a cylinder.
  • Ribbon: This variant is useful for modeling shapes that do not actually have a cylindrical cross section (such as a blade of grass).

Figure 6.32: The Three Types of Curves That the Curve Shape Can Represent. On the top is a flat curve that is always oriented to be perpendicular to a ray approaching it. The middle is a variant of this curve where the shading normal is set so that the curve appears to be cylindrical. On the bottom is a ribbon, which has a fixed orientation at its starting and ending points; intermediate orientations are smoothly interpolated between them.

The CurveType enumerator records which of them a given Curve instance models.

The flat and cylinder curve variants are intended to be used as convenient approximations of deformed cylinders. It should be noted that intersections found with respect to them do not correspond to a physically realizable 3D shape, which can potentially lead to minor inconsistencies when taking a scene with true cylinders as a reference.

<<CurveType Definition>>= 
enum class CurveType { Flat, Cylinder, Ribbon };

Given a curve specified in a pbrt scene description file, it can be worthwhile to split it into a few segments, each covering part of the u parametric range of the curve. (One reason for doing so is that axis-aligned bounding boxes do not tightly bound wiggly curves, but subdividing Bézier curves makes them less wiggly—the variation diminishing property of polynomials.) Therefore, the Curve constructor takes a parametric range of u values, left-bracket u Subscript normal m normal i normal n Baseline comma u Subscript normal m normal a normal x Baseline right-bracket , as well as a pointer to a CurveCommon structure, which stores the control points and other information about the curve that is shared across curve segments. In this way, the memory footprint for individual curve segments is reduced, which makes it easier to keep many of them in memory.

<<Curve Public Methods>>= 
Curve(const CurveCommon *common, Float uMin, Float uMax) : common(common), uMin(uMin), uMax(uMax) {}

<<Curve Private Members>>= 
const CurveCommon *common; Float uMin, uMax;

The CurveCommon constructor initializes member variables with values passed into it for the control points, the curve width, etc. The control points provided to it should be in the curve’s object space.

For Ribbon curves, CurveCommon stores a surface normal to orient the curve at each endpoint. The constructor precomputes the angle between the two normal vectors and one over the sine of this angle; these values will be useful when computing the orientation of the curve at arbitrary points along its extent.

<<CurveCommon Definition>>= 
struct CurveCommon { <<CurveCommon Public Methods>> 
CurveCommon(pstd::span<const Point3f> c, Float w0, Float w1, CurveType type, pstd::span<const Normal3f> norm, const Transform *renderFromObject, const Transform *objectFromRender, bool reverseOrientation); std::string ToString() const;
<<CurveCommon Public Members>> 
CurveType type; Point3f cpObj[4]; Float width[2]; Normal3f n[2]; Float normalAngle, invSinNormalAngle; const Transform *renderFromObject, *objectFromRender; bool reverseOrientation, transformSwapsHandedness;
};

<<CurveCommon Public Members>>= 
CurveType type; Point3f cpObj[4]; Float width[2]; Normal3f n[2]; Float normalAngle, invSinNormalAngle; const Transform *renderFromObject, *objectFromRender; bool reverseOrientation, transformSwapsHandedness;

6.7.1 Bounding Curves

The object-space bound of a curve can be found by first bounding the spline along the center of the curve and then expanding that bound by half the maximum width the curve takes on over its extent. The Bounds() method then transforms that bound to rendering space before returning it.

<<Curve Method Definitions>>= 
Bounds3f Curve::Bounds() const { pstd::span<const Point3f> cpSpan(common->cpObj); Bounds3f objBounds = BoundCubicBezier(cpSpan, uMin, uMax); <<Expand objBounds by maximum curve width over u range>> 
Float width[2] = {Lerp(uMin, common->width[0], common->width[1]), Lerp(uMax, common->width[0], common->width[1])}; objBounds = Expand(objBounds, std::max(width[0], width[1]) * 0.5f);
return (*common->renderFromObject)(objBounds); }

<<Expand objBounds by maximum curve width over u range>>= 
Float width[2] = {Lerp(uMin, common->width[0], common->width[1]), Lerp(uMax, common->width[0], common->width[1])}; objBounds = Expand(objBounds, std::max(width[0], width[1]) * 0.5f);

The Curve shape cannot be used as an area light, as it does not provide implementations of the required sampling methods. It does provide a NormalBounds() method that returns a conservative bound.

<<Curve Public Methods>>+= 
DirectionCone NormalBounds() const { return DirectionCone::EntireSphere(); }

6.7.2 Intersection Tests

Both of the intersection methods required by the Shape interface are implemented via another Curve method, IntersectRay(). Rather than returning an optional ShapeIntersection, it takes a pointer to one.

<<Curve Method Definitions>>+=  
pstd::optional<ShapeIntersection> Curve::Intersect(const Ray &ray, Float tMax) const { pstd::optional<ShapeIntersection> si; IntersectRay(ray, tMax, &si); return si; }

IntersectP() passes nullptr to IntersectRay(), which indicates that it can return immediately if an intersection is found.

<<Curve Method Definitions>>+=  
bool Curve::IntersectP(const Ray &ray, Float tMax) const { return IntersectRay(ray, tMax, nullptr); }

The Curve intersection algorithm is based on discarding curve segments as soon as it can be determined that the ray definitely does not intersect them and otherwise recursively splitting the curve in half to create two smaller segments that are then tested. Eventually, the curve is linearly approximated for an efficient intersection test. That process starts after some initial preparation and early culling tests in IntersectRay().

<<Curve Method Definitions>>+=  
bool Curve::IntersectRay(const Ray &r, Float tMax, pstd::optional<ShapeIntersection> *si) const { <<Transform Ray to curve’s object space>> 
Ray ray = (*common->objectFromRender)(r);
<<Get object-space control points for curve segment, cpObj>> 
pstd::array<Point3f, 4> cpObj = CubicBezierControlPoints(pstd::span<const Point3f>(common->cpObj), uMin, uMax);
<<Project curve control points to plane perpendicular to ray>> 
Vector3f dx = Cross(ray.d, cpObj[3] - cpObj[0]); if (LengthSquared(dx) == 0) { Vector3f dy; CoordinateSystem(ray.d, &dx, &dy); } Transform rayFromObject = LookAt(ray.o, ray.o + ray.d, dx); pstd::array<Point3f, 4> cp = { rayFromObject(cpObj[0]), rayFromObject(cpObj[1]), rayFromObject(cpObj[2]), rayFromObject(cpObj[3]) };
<<Test ray against bound of projected control points>> 
Float maxWidth = std::max(Lerp(uMin, common->width[0], common->width[1]), Lerp(uMax, common->width[0], common->width[1])); Bounds3f curveBounds = Union(Bounds3f(cp[0], cp[1]), Bounds3f(cp[2], cp[3])); curveBounds = Expand(curveBounds, 0.5f * maxWidth); Bounds3f rayBounds(Point3f(0, 0, 0), Point3f(0, 0, Length(ray.d) * tMax)); if (!Overlaps(rayBounds, curveBounds)) return false;
<<Compute refinement depth for curve, maxDepth>> 
Float L0 = 0; for (int i = 0; i < 2; ++i) L0 = std::max( L0, std::max(std::max(std::abs(cp[i].x - 2 * cp[i + 1].x + cp[i + 2].x), std::abs(cp[i].y - 2 * cp[i + 1].y + cp[i + 2].y)), std::abs(cp[i].z - 2 * cp[i + 1].z + cp[i + 2].z))); int maxDepth = 0; if (L0 > 0) { Float eps = std::max(common->width[0], common->width[1]) * .05f; // width / 20 // Compute log base 4 by dividing log2 in half. int r0 = Log2Int(1.41421356237f * 6.f * L0 / (8.f * eps)) / 2; maxDepth = Clamp(r0, 0, 10); }
<<Recursively test for ray–curve intersection>> 
pstd::span<const Point3f> cpSpan(cp); return RecursiveIntersect(ray, tMax, cpSpan, Inverse(rayFromObject), uMin, uMax, maxDepth, si);
}

<<Transform Ray to curve’s object space>>= 

The CurveCommon class stores the control points for the full curve, but a Curve instance generally needs the four control points that represent the Bézier curve for its u extent. The CubicBezierControlPoints() utility function performs this computation.

<<Get object-space control points for curve segment, cpObj>>= 
pstd::array<Point3f, 4> cpObj = CubicBezierControlPoints(pstd::span<const Point3f>(common->cpObj), uMin, uMax);

Like the ray–triangle intersection algorithm from Section 6.5.3, the ray–curve intersection test is based on transforming the curve to a coordinate system with the ray’s origin at the origin of the coordinate system and the ray’s direction aligned to be along the plus z axis. Performing this transformation at the start greatly reduces the number of operations that must be performed for intersection tests.

For the Curve shape, we will need an explicit representation of the transformation, so the LookAt() function is used to generate it here. The origin is the ray’s origin and the “look at” point is a point offset from the origin along the ray’s direction. The “up” direction is set to be perpendicular to both the ray’s direction and the vector from the first to the last control point. Doing so helps orient the curve to be roughly parallel to the x axis in the ray coordinate system, which in turn leads to tighter bounds in y (see Figure 6.33). This improvement in the fit of the bounds often makes it possible to terminate the recursive intersection tests earlier than would be possible otherwise.

Figure 6.33: 2D Bounding Boxes of a Bézier Curve. (a) Bounding box computed using the curve’s control points as given. (b) The effect of rotating the curve so that the vector from its first to last control point is aligned with the x axis before computing bounds. The resulting bounding box is a much tighter fit.

If the ray and the vector between the first and last control points are parallel, dx will be degenerate. In that case we find an arbitrary “up” vector direction so that intersection tests can proceed in this unusual case.

<<Project curve control points to plane perpendicular to ray>>= 
Vector3f dx = Cross(ray.d, cpObj[3] - cpObj[0]); if (LengthSquared(dx) == 0) { Vector3f dy; CoordinateSystem(ray.d, &dx, &dy); } Transform rayFromObject = LookAt(ray.o, ray.o + ray.d, dx); pstd::array<Point3f, 4> cp = { rayFromObject(cpObj[0]), rayFromObject(cpObj[1]), rayFromObject(cpObj[2]), rayFromObject(cpObj[3]) };

Along the lines of the implementation in Curve::Bounds(), a conservative bounding box for a curve segment can be found by taking the bounds of the curve’s control points and expanding by half of the maximum width of the curve over the u range being considered.

Figure 6.34: Ray–Curve Bounds Test. In the ray coordinate system, the ray’s origin is at left-parenthesis 0 comma 0 comma 0 right-parenthesis and its direction is aligned with the plus z axis. Therefore, if the 2D point left-parenthesis x comma y right-parenthesis equals left-parenthesis 0 comma 0 right-parenthesis is outside the x y bounding box of the curve segment, then it is impossible that the ray intersects the curve.

Because the ray’s origin is at left-parenthesis 0 comma 0 comma 0 right-parenthesis and its direction is aligned with the plus z axis in the intersection space, its bounding box only includes the origin in x and y (Figure 6.34); its z extent is given by the z range that its parametric extent covers. Before proceeding with the recursive intersection testing algorithm, the ray’s bounding box is tested for intersection with the curve’s bounding box. The method can return immediately if they do not intersect.

<<Test ray against bound of projected control points>>= 
Float maxWidth = std::max(Lerp(uMin, common->width[0], common->width[1]), Lerp(uMax, common->width[0], common->width[1])); Bounds3f curveBounds = Union(Bounds3f(cp[0], cp[1]), Bounds3f(cp[2], cp[3])); curveBounds = Expand(curveBounds, 0.5f * maxWidth); Bounds3f rayBounds(Point3f(0, 0, 0), Point3f(0, 0, Length(ray.d) * tMax)); if (!Overlaps(rayBounds, curveBounds)) return false;

The maximum number of times to subdivide the curve is computed so that the maximum distance from the eventual linearized curve at the finest refinement level is bounded to be less than a small fixed distance. We will not go into the details of this computation, which is implemented in the fragment <<Compute refinement depth for curve, maxDepth>>. With the culling tests passed and that value in hand, the recursive intersection tests begin.

<<Recursively test for ray–curve intersection>>= 
pstd::span<const Point3f> cpSpan(cp); return RecursiveIntersect(ray, tMax, cpSpan, Inverse(rayFromObject), uMin, uMax, maxDepth, si);

The RecursiveIntersect() method then tests whether the given ray intersects the given curve segment over the given parametric range left-bracket monospace u Baseline monospace 0 monospace comma monospace u Baseline monospace 1 monospace right-bracket . It assumes that the ray has already been tested against the curve’s bounding box and found to intersect it.

<<Curve Method Definitions>>+= 
bool Curve::RecursiveIntersect( const Ray &ray, Float tMax, pstd::span<const Point3f> cp, const Transform &objectFromRay, Float u0, Float u1, int depth, pstd::optional<ShapeIntersection> *si) const { Float rayLength = Length(ray.d); if (depth > 0) { <<Split curve segment into subsegments and test for intersection>> 
pstd::array<Point3f, 7> cpSplit = SubdivideCubicBezier(cp); Float u[3] = {u0, (u0 + u1) / 2, u1}; for (int seg = 0; seg < 2; ++seg) { <<Check ray against curve segment’s bounding box>> 
Float maxWidth = std::max(Lerp(u[seg], common->width[0], common->width[1]), Lerp(u[seg + 1], common->width[0], common->width[1])); pstd::span<const Point3f> cps = pstd::MakeConstSpan(&cpSplit[3 * seg], 4); Bounds3f curveBounds = Union(Bounds3f(cps[0], cps[1]), Bounds3f(cps[2], cps[3])); curveBounds = Expand(curveBounds, 0.5f * maxWidth); Bounds3f rayBounds(Point3f(0, 0, 0), Point3f(0, 0, Length(ray.d) * tMax)); if (!Overlaps(rayBounds, curveBounds)) continue;
<<Recursively test ray-segment intersection>> 
bool hit = RecursiveIntersect(ray, tMax, cps, objectFromRay, u[seg], u[seg + 1], depth - 1, si); if (hit && !si) return true;
} return si ? si->has_value() : false;
} else { <<Intersect ray with curve segment>> 
<<Test ray against segment endpoint boundaries>> 
<<Test sample point against tangent perpendicular at curve start>> 
Float edge = (cp[1].y - cp[0].y) * -cp[0].y + cp[0].x * (cp[0].x - cp[1].x); if (edge < 0) return false;
<<Test sample point against tangent perpendicular at curve end>> 
edge = (cp[2].y - cp[3].y) * -cp[3].y + cp[3].x * (cp[3].x - cp[2].x); if (edge < 0) return false;
<<Find line w that gives minimum distance to sample point>> 
Vector2f segmentDir = Point2f(cp[3].x, cp[3].y) - Point2f(cp[0].x, cp[0].y); Float denom = LengthSquared(segmentDir); if (denom == 0) return false; Float w = Dot(-Vector2f(cp[0].x, cp[0].y), segmentDir) / denom;
<<Compute u coordinate of curve intersection point and hitWidth>> 
Float u = Clamp(Lerp(w, u0, u1), u0, u1); Float hitWidth = Lerp(u, common->width[0], common->width[1]); Normal3f nHit; if (common->type == CurveType::Ribbon) { <<Scale hitWidth based on ribbon orientation>> 
if (common->normalAngle == 0) nHit = common->n[0]; else { Float sin0 = std::sin((1 - u) * common->normalAngle) * common->invSinNormalAngle; Float sin1 = std::sin(u * common->normalAngle) * common->invSinNormalAngle; nHit = sin0 * common->n[0] + sin1 * common->n[1]; } hitWidth *= AbsDot(nHit, ray.d) / rayLength;
}
<<Test intersection point against curve width>> 
Vector3f dpcdw; Point3f pc = EvaluateCubicBezier(pstd::span<const Point3f>(cp), Clamp(w, 0, 1), &dpcdw); Float ptCurveDist2 = Sqr(pc.x) + Sqr(pc.y); if (ptCurveDist2 > Sqr(hitWidth) * 0.25f) return false; if (pc.z < 0 || pc.z > rayLength * tMax) return false;
if (si) { <<Initialize ShapeIntersection for curve intersection>> 
<<Compute tHit for curve intersection>> 
Float tHit = pc.z / rayLength; if (si->has_value() && tHit > si->value().tHit) return false;
<<Initialize SurfaceInteraction intr for curve intersection>> 
<<Compute v coordinate of curve intersection point>> 
Float ptCurveDist = std::sqrt(ptCurveDist2); Float edgeFunc = dpcdw.x * -pc.y + pc.x * dpcdw.y; Float v = (edgeFunc > 0) ? 0.5f + ptCurveDist / hitWidth : 0.5f - ptCurveDist / hitWidth;
<<Compute partial-differential normal p slash partial-differential u and partial-differential normal p slash partial-differential v for curve intersection>> 
Vector3f dpdu, dpdv; EvaluateCubicBezier(pstd::MakeConstSpan(common->cpObj), u, &dpdu); if (common->type == CurveType::Ribbon) dpdv = Normalize(Cross(nHit, dpdu)) * hitWidth; else { <<Compute curve partial-differential normal p slash partial-differential v for flat and cylinder curves>> 
Vector3f dpduPlane = objectFromRay.ApplyInverse(dpdu); Vector3f dpdvPlane = Normalize(Vector3f(-dpduPlane.y, dpduPlane.x, 0)) * hitWidth; if (common->type == CurveType::Cylinder) { <<Rotate dpdvPlane to give cylindrical appearance>> 
Float theta = Lerp(v, -90, 90); Transform rot = Rotate(-theta, dpduPlane); dpdvPlane = rot(dpdvPlane);
} dpdv = objectFromRay(dpdvPlane);
}
<<Compute error bounds for curve intersection>> 
Vector3f pError(hitWidth, hitWidth, hitWidth);
bool flipNormal = common->reverseOrientation ^ common->transformSwapsHandedness; Point3fi pi(ray(tHit), pError); SurfaceInteraction intr(pi, {u, v}, -ray.d, dpdu, dpdv, Normal3f(), Normal3f(), ray.time, flipNormal); intr = (*common->renderFromObject)(intr);
*si = ShapeIntersection{intr, tHit};
} return true;
} }

If the maximum depth has not been reached, a call to SubdivideCubicBezier() gives the control points for the two Bézier curves that result in splitting the Bézier curve given by cp in half. The last control point of the first curve is the same as the first control point of the second, so 7 values are returned for the total of 8 control points. The u array is then initialized so that it holds the parametric range of the two curves before each one is processed in turn.

<<Split curve segment into subsegments and test for intersection>>= 
pstd::array<Point3f, 7> cpSplit = SubdivideCubicBezier(cp); Float u[3] = {u0, (u0 + u1) / 2, u1}; for (int seg = 0; seg < 2; ++seg) { <<Check ray against curve segment’s bounding box>> 
Float maxWidth = std::max(Lerp(u[seg], common->width[0], common->width[1]), Lerp(u[seg + 1], common->width[0], common->width[1])); pstd::span<const Point3f> cps = pstd::MakeConstSpan(&cpSplit[3 * seg], 4); Bounds3f curveBounds = Union(Bounds3f(cps[0], cps[1]), Bounds3f(cps[2], cps[3])); curveBounds = Expand(curveBounds, 0.5f * maxWidth); Bounds3f rayBounds(Point3f(0, 0, 0), Point3f(0, 0, Length(ray.d) * tMax)); if (!Overlaps(rayBounds, curveBounds)) continue;
<<Recursively test ray-segment intersection>> 
bool hit = RecursiveIntersect(ray, tMax, cps, objectFromRay, u[seg], u[seg + 1], depth - 1, si); if (hit && !si) return true;
} return si ? si->has_value() : false;

The bounding box test in the <<Check ray against curve segment’s bounding box>> fragment is essentially the same as the one in <<Test ray against bound of projected control points>> except that it takes u values from the u array when computing the curve’s maximum width over the u range and it uses control points from cpSplit. Therefore, it is not included here.

If the ray does intersect the bounding box, the corresponding segment is given to a recursive call of RecursiveIntersect(). If an intersection is found and the ray is a shadow ray, si will be nullptr and an intersection can immediately be reported. For non-shadow rays, even if an intersection has been found, it may not be the closest intersection, so the other segment still must be considered.

<<Recursively test ray-segment intersection>>= 
bool hit = RecursiveIntersect(ray, tMax, cps, objectFromRay, u[seg], u[seg + 1], depth - 1, si); if (hit && !si) return true;

The intersection test is made more efficient by using a linear approximation of the curve; the variation diminishing property allows us to make this approximation without introducing too much error.

<<Intersect ray with curve segment>>= 
<<Test ray against segment endpoint boundaries>> 
<<Test sample point against tangent perpendicular at curve start>> 
Float edge = (cp[1].y - cp[0].y) * -cp[0].y + cp[0].x * (cp[0].x - cp[1].x); if (edge < 0) return false;
<<Test sample point against tangent perpendicular at curve end>> 
edge = (cp[2].y - cp[3].y) * -cp[3].y + cp[3].x * (cp[3].x - cp[2].x); if (edge < 0) return false;
<<Find line w that gives minimum distance to sample point>> 
Vector2f segmentDir = Point2f(cp[3].x, cp[3].y) - Point2f(cp[0].x, cp[0].y); Float denom = LengthSquared(segmentDir); if (denom == 0) return false; Float w = Dot(-Vector2f(cp[0].x, cp[0].y), segmentDir) / denom;
<<Compute u coordinate of curve intersection point and hitWidth>> 
Float u = Clamp(Lerp(w, u0, u1), u0, u1); Float hitWidth = Lerp(u, common->width[0], common->width[1]); Normal3f nHit; if (common->type == CurveType::Ribbon) { <<Scale hitWidth based on ribbon orientation>> 
if (common->normalAngle == 0) nHit = common->n[0]; else { Float sin0 = std::sin((1 - u) * common->normalAngle) * common->invSinNormalAngle; Float sin1 = std::sin(u * common->normalAngle) * common->invSinNormalAngle; nHit = sin0 * common->n[0] + sin1 * common->n[1]; } hitWidth *= AbsDot(nHit, ray.d) / rayLength;
}
<<Test intersection point against curve width>> 
Vector3f dpcdw; Point3f pc = EvaluateCubicBezier(pstd::span<const Point3f>(cp), Clamp(w, 0, 1), &dpcdw); Float ptCurveDist2 = Sqr(pc.x) + Sqr(pc.y); if (ptCurveDist2 > Sqr(hitWidth) * 0.25f) return false; if (pc.z < 0 || pc.z > rayLength * tMax) return false;
if (si) { <<Initialize ShapeIntersection for curve intersection>> 
<<Compute tHit for curve intersection>> 
Float tHit = pc.z / rayLength; if (si->has_value() && tHit > si->value().tHit) return false;
<<Initialize SurfaceInteraction intr for curve intersection>> 
<<Compute v coordinate of curve intersection point>> 
Float ptCurveDist = std::sqrt(ptCurveDist2); Float edgeFunc = dpcdw.x * -pc.y + pc.x * dpcdw.y; Float v = (edgeFunc > 0) ? 0.5f + ptCurveDist / hitWidth : 0.5f - ptCurveDist / hitWidth;
<<Compute partial-differential normal p slash partial-differential u and partial-differential normal p slash partial-differential v for curve intersection>> 
Vector3f dpdu, dpdv; EvaluateCubicBezier(pstd::MakeConstSpan(common->cpObj), u, &dpdu); if (common->type == CurveType::Ribbon) dpdv = Normalize(Cross(nHit, dpdu)) * hitWidth; else { <<Compute curve partial-differential normal p slash partial-differential v for flat and cylinder curves>> 
Vector3f dpduPlane = objectFromRay.ApplyInverse(dpdu); Vector3f dpdvPlane = Normalize(Vector3f(-dpduPlane.y, dpduPlane.x, 0)) * hitWidth; if (common->type == CurveType::Cylinder) { <<Rotate dpdvPlane to give cylindrical appearance>> 
Float theta = Lerp(v, -90, 90); Transform rot = Rotate(-theta, dpduPlane); dpdvPlane = rot(dpdvPlane);
} dpdv = objectFromRay(dpdvPlane);
}
<<Compute error bounds for curve intersection>> 
Vector3f pError(hitWidth, hitWidth, hitWidth);
bool flipNormal = common->reverseOrientation ^ common->transformSwapsHandedness; Point3fi pi(ray(tHit), pError); SurfaceInteraction intr(pi, {u, v}, -ray.d, dpdu, dpdv, Normal3f(), Normal3f(), ray.time, flipNormal); intr = (*common->renderFromObject)(intr);
*si = ShapeIntersection{intr, tHit};
} return true;

Figure 6.35: Curve Segment Boundaries. The intersection test for a segment of a larger curve computes edge functions for the lines that are perpendicular to the segment endpoints (dashed lines). If a potential intersection point (solid dot) is on the other side of the edge than the segment, it is rejected; another curve segment (if present on that side) should account for this intersection instead.

It is important that the intersection test only accepts intersections that are on the Curve’s surface for the u segment currently under consideration. Therefore, the first step of the intersection test is to compute edge functions for lines perpendicular to the curve starting point and ending point and to classify the potential intersection point against them (Figure 6.35).

<<Test ray against segment endpoint boundaries>>= 
<<Test sample point against tangent perpendicular at curve start>> 
Float edge = (cp[1].y - cp[0].y) * -cp[0].y + cp[0].x * (cp[0].x - cp[1].x); if (edge < 0) return false;
<<Test sample point against tangent perpendicular at curve end>> 
edge = (cp[2].y - cp[3].y) * -cp[3].y + cp[3].x * (cp[3].x - cp[2].x); if (edge < 0) return false;

Projecting the curve control points into the ray coordinate system makes this test more efficient for two reasons. First, because the ray’s direction is oriented with the plus z axis, the problem is reduced to a 2D test in x and y . Second, because the ray origin is at the origin of the coordinate system, the point we need to classify is left-parenthesis 0 comma 0 right-parenthesis , which simplifies evaluating the edge function, just like the ray–triangle intersection test.

Edge functions were introduced for ray–triangle intersection tests in Equation (6.5); see also Figure 6.14. To define the edge function here, we need any two points on the line perpendicular to the curve going through the starting point. The first control point, normal p 0 , is a fine choice for the first. For the second, we will compute the vector perpendicular to the curve’s tangent and add that offset to the control point.

Differentiation of Equation (6.16) shows that the tangent to the curve at the first control point normal p 0 is 3 left-parenthesis normal p 1 minus normal p 0 right-parenthesis . The scaling factor does not matter here, so we will use bold t equals normal p 1 minus normal p 0 here. Computing the vector perpendicular to the tangent is easy in 2D: it is just necessary to swap the x and y coordinates and negate one of them. (To see why this works, consider the dot product left-parenthesis x comma y right-parenthesis dot left-parenthesis y comma negative x right-parenthesis equals x y plus minus y x equals 0 . Because the cosine of the angle between the two vectors is zero, they must be perpendicular.) Thus, the second point on the edge is

normal p 0 plus left-parenthesis normal p 1 Subscript y Baseline minus normal p 0 Subscript y Baseline comma minus left-parenthesis normal p 1 Subscript x Baseline minus normal p 0 Subscript x Baseline right-parenthesis right-parenthesis equals normal p 0 plus left-parenthesis normal p 1 Subscript y Baseline minus normal p 0 Subscript y Baseline comma normal p 0 Subscript x Baseline minus normal p 1 Subscript x Baseline right-parenthesis period

Substituting these two points into the definition of the edge function, Equation (6.5), and simplifying gives

e left-parenthesis normal p Subscript Baseline right-parenthesis equals left-parenthesis normal p 1 Subscript y Baseline minus normal p 0 Subscript y Baseline right-parenthesis left-parenthesis normal p Subscript Baseline Subscript y Baseline minus normal p 0 Subscript y Baseline right-parenthesis minus left-parenthesis normal p Subscript Baseline Subscript x Baseline minus normal p 0 Subscript x Baseline right-parenthesis left-parenthesis normal p 0 Subscript x Baseline minus normal p 1 Subscript x Baseline right-parenthesis period

Finally, substituting normal p Subscript Baseline equals left-parenthesis 0 comma 0 right-parenthesis gives the final expression to test:

e left-parenthesis left-parenthesis 0 comma 0 right-parenthesis right-parenthesis equals left-parenthesis normal p 1 Subscript y Baseline minus normal p 0 Subscript y Baseline right-parenthesis left-parenthesis minus normal p 0 Subscript y Baseline right-parenthesis plus normal p 0 Subscript x Baseline left-parenthesis normal p 0 Subscript x Baseline minus normal p 1 Subscript x Baseline right-parenthesis period

<<Test sample point against tangent perpendicular at curve start>>= 
Float edge = (cp[1].y - cp[0].y) * -cp[0].y + cp[0].x * (cp[0].x - cp[1].x); if (edge < 0) return false;

The <<Test sample point against tangent perpendicular at curve end>> fragment, not included here, does the corresponding test at the end of the curve.

The next part of the test is to determine the u value along the curve segment where the point left-parenthesis 0 comma 0 right-parenthesis is closest to the curve. This will be the intersection point, if it is no farther than the curve’s width away from the center at that point. Determining this distance for a cubic Bézier curve requires a significant amount of computation, so instead the implementation here approximates the curve with a linear segment to compute this u value.

Figure 6.36: Approximation of a Cubic Bézier Curve with a Linear Segment. For this part of the ray–curve intersection test, we approximate the Bézier with a linear segment (dashed line) passing through its starting and ending points. (In practice, after being subdivided, the curve will be already nearly linear, so the error is less than this figure suggests.)

We linearly approximate the Bézier curve with a line segment from its starting point normal p 0 to its endpoint normal p 3 that is parameterized by w . In this case, the position is normal p 0 at w equals 0 and normal p 3 at w equals 1 (Figure 6.36). Our task is to compute the value of w along the line corresponding to the point on the line normal p prime that is closest to the point normal p Subscript . The key insight to apply is that at normal p prime , the vector from the corresponding point on the line to normal p Subscript will be perpendicular to the line (Figure 6.37(a)).

Figure 6.37: (a) Given an infinite line and a point normal p , the vector from the point to the closest point on the line, normal p prime , is then perpendicular to the line. (b) Because this vector is perpendicular, we can compute the distance from the first point of the line to the point of closest approach, normal p prime , as d equals double-vertical-bar normal p Subscript Baseline minus normal p 0 double-vertical-bar cosine theta .

Equation (3.1) gives us a relationship between the dot product of two vectors, their lengths, and the cosine of the angle between them. In particular, it shows us how to compute the cosine of the angle between the vector from normal p 0 to normal p Subscript and the vector from normal p 0 to normal p 3 :

cosine theta equals StartFraction left-parenthesis normal p Subscript Baseline minus normal p 0 right-parenthesis dot left-parenthesis normal p 3 minus normal p 0 right-parenthesis Over double-vertical-bar normal p Subscript Baseline minus normal p 0 double-vertical-bar double-vertical-bar normal p 3 minus normal p 0 double-vertical-bar EndFraction period

Because the vector from normal p prime to normal p Subscript is perpendicular to the line (Figure 6.37(b)), we can compute the distance along the line from normal p 0 to normal p prime as

d equals double-vertical-bar normal p Subscript Baseline minus normal p 0 double-vertical-bar cosine theta equals StartFraction left-parenthesis normal p Subscript Baseline minus normal p 0 right-parenthesis dot left-parenthesis normal p 3 minus normal p 0 right-parenthesis Over double-vertical-bar normal p 3 minus normal p 0 double-vertical-bar EndFraction period

Finally, the parametric offset w along the line is the ratio of d to the line’s length,

w equals StartFraction d Over double-vertical-bar normal p 3 minus normal p 0 double-vertical-bar EndFraction equals StartFraction left-parenthesis normal p Subscript Baseline minus normal p 0 right-parenthesis dot left-parenthesis normal p 3 minus normal p 0 right-parenthesis Over double-vertical-bar normal p 3 minus normal p 0 double-vertical-bar squared EndFraction period

The computation of the value of w is in turn slightly simplified from the fact that normal p Subscript Baseline equals left-parenthesis 0 comma 0 right-parenthesis in the intersection coordinate system.

<<Find line w that gives minimum distance to sample point>>= 
Vector2f segmentDir = Point2f(cp[3].x, cp[3].y) - Point2f(cp[0].x, cp[0].y); Float denom = LengthSquared(segmentDir); if (denom == 0) return false; Float w = Dot(-Vector2f(cp[0].x, cp[0].y), segmentDir) / denom;

The parametric u coordinate of the (presumed) closest point on the Bézier curve to the candidate intersection point is computed by linearly interpolating along the u range of the segment. Given this u value, the width of the curve at that point can be computed.

<<Compute u coordinate of curve intersection point and hitWidth>>= 
Float u = Clamp(Lerp(w, u0, u1), u0, u1); Float hitWidth = Lerp(u, common->width[0], common->width[1]); Normal3f nHit; if (common->type == CurveType::Ribbon) { <<Scale hitWidth based on ribbon orientation>> 
if (common->normalAngle == 0) nHit = common->n[0]; else { Float sin0 = std::sin((1 - u) * common->normalAngle) * common->invSinNormalAngle; Float sin1 = std::sin(u * common->normalAngle) * common->invSinNormalAngle; nHit = sin0 * common->n[0] + sin1 * common->n[1]; } hitWidth *= AbsDot(nHit, ray.d) / rayLength;
}

For Ribbon curves, the curve is not always oriented to face the ray. Rather, its orientation is interpolated between two surface normals given at each endpoint. Here, spherical linear interpolation is used to interpolate the normal at u . The curve’s width is then scaled by the cosine of the angle between the normalized ray direction and the ribbon’s orientation so that it corresponds to the visible width of the curve from the given direction.

<<Scale hitWidth based on ribbon orientation>>= 
if (common->normalAngle == 0) nHit = common->n[0]; else { Float sin0 = std::sin((1 - u) * common->normalAngle) * common->invSinNormalAngle; Float sin1 = std::sin(u * common->normalAngle) * common->invSinNormalAngle; nHit = sin0 * common->n[0] + sin1 * common->n[1]; } hitWidth *= AbsDot(nHit, ray.d) / rayLength;

To finally classify the potential intersection as a hit or miss, the Bézier curve must still be evaluated at u . (Because the control points cp represent the curve segment currently under consideration, it is important to use w rather than u in the function call, however, since w is in the range left-bracket 0 comma 1 right-bracket .) The derivative of the curve at this point will be useful shortly, so it is recorded now.

We would like to test whether the distance from normal p Subscript to this point on the curve pc is less than half the curve’s width. Because normal p Subscript Baseline equals left-parenthesis 0 comma 0 right-parenthesis , we can equivalently test whether the distance from pc to the origin is less than half the width or whether the squared distance is less than one quarter the width squared. If this test passes, the last thing to check is if the intersection point is in the ray’s parametric t range.

<<Test intersection point against curve width>>= 
Vector3f dpcdw; Point3f pc = EvaluateCubicBezier(pstd::span<const Point3f>(cp), Clamp(w, 0, 1), &dpcdw); Float ptCurveDist2 = Sqr(pc.x) + Sqr(pc.y); if (ptCurveDist2 > Sqr(hitWidth) * 0.25f) return false; if (pc.z < 0 || pc.z > rayLength * tMax) return false;

For non-shadow rays, the ShapeIntersection for the intersection can finally be initialized. Doing so requires computing the ray t value for the intersection as well as its SurfaceInteraction.

<<Initialize ShapeIntersection for curve intersection>>= 
<<Compute tHit for curve intersection>> 
Float tHit = pc.z / rayLength; if (si->has_value() && tHit > si->value().tHit) return false;
<<Initialize SurfaceInteraction intr for curve intersection>> 
<<Compute v coordinate of curve intersection point>> 
Float ptCurveDist = std::sqrt(ptCurveDist2); Float edgeFunc = dpcdw.x * -pc.y + pc.x * dpcdw.y; Float v = (edgeFunc > 0) ? 0.5f + ptCurveDist / hitWidth : 0.5f - ptCurveDist / hitWidth;
<<Compute partial-differential normal p slash partial-differential u and partial-differential normal p slash partial-differential v for curve intersection>> 
Vector3f dpdu, dpdv; EvaluateCubicBezier(pstd::MakeConstSpan(common->cpObj), u, &dpdu); if (common->type == CurveType::Ribbon) dpdv = Normalize(Cross(nHit, dpdu)) * hitWidth; else { <<Compute curve partial-differential normal p slash partial-differential v for flat and cylinder curves>> 
Vector3f dpduPlane = objectFromRay.ApplyInverse(dpdu); Vector3f dpdvPlane = Normalize(Vector3f(-dpduPlane.y, dpduPlane.x, 0)) * hitWidth; if (common->type == CurveType::Cylinder) { <<Rotate dpdvPlane to give cylindrical appearance>> 
Float theta = Lerp(v, -90, 90); Transform rot = Rotate(-theta, dpduPlane); dpdvPlane = rot(dpdvPlane);
} dpdv = objectFromRay(dpdvPlane);
}
<<Compute error bounds for curve intersection>> 
Vector3f pError(hitWidth, hitWidth, hitWidth);
bool flipNormal = common->reverseOrientation ^ common->transformSwapsHandedness; Point3fi pi(ray(tHit), pError); SurfaceInteraction intr(pi, {u, v}, -ray.d, dpdu, dpdv, Normal3f(), Normal3f(), ray.time, flipNormal); intr = (*common->renderFromObject)(intr);
*si = ShapeIntersection{intr, tHit};

After the tHit value has been computed, it is compared against the tHit of a previously found ray–curve intersection, if there is one. This check ensures that the closest intersection is returned.

<<Compute tHit for curve intersection>>= 
Float tHit = pc.z / rayLength; if (si->has_value() && tHit > si->value().tHit) return false;

A variety of additional quantities need to be computed in order to be able to initialize the intersection’s SurfaceInteraction.

<<Initialize SurfaceInteraction intr for curve intersection>>= 
<<Compute v coordinate of curve intersection point>> 
Float ptCurveDist = std::sqrt(ptCurveDist2); Float edgeFunc = dpcdw.x * -pc.y + pc.x * dpcdw.y; Float v = (edgeFunc > 0) ? 0.5f + ptCurveDist / hitWidth : 0.5f - ptCurveDist / hitWidth;
<<Compute partial-differential normal p slash partial-differential u and partial-differential normal p slash partial-differential v for curve intersection>> 
Vector3f dpdu, dpdv; EvaluateCubicBezier(pstd::MakeConstSpan(common->cpObj), u, &dpdu); if (common->type == CurveType::Ribbon) dpdv = Normalize(Cross(nHit, dpdu)) * hitWidth; else { <<Compute curve partial-differential normal p slash partial-differential v for flat and cylinder curves>> 
Vector3f dpduPlane = objectFromRay.ApplyInverse(dpdu); Vector3f dpdvPlane = Normalize(Vector3f(-dpduPlane.y, dpduPlane.x, 0)) * hitWidth; if (common->type == CurveType::Cylinder) { <<Rotate dpdvPlane to give cylindrical appearance>> 
Float theta = Lerp(v, -90, 90); Transform rot = Rotate(-theta, dpduPlane); dpdvPlane = rot(dpdvPlane);
} dpdv = objectFromRay(dpdvPlane);
}
<<Compute error bounds for curve intersection>> 
Vector3f pError(hitWidth, hitWidth, hitWidth);
bool flipNormal = common->reverseOrientation ^ common->transformSwapsHandedness; Point3fi pi(ray(tHit), pError); SurfaceInteraction intr(pi, {u, v}, -ray.d, dpdu, dpdv, Normal3f(), Normal3f(), ray.time, flipNormal); intr = (*common->renderFromObject)(intr);

We have gotten this far without computing the v coordinate of the intersection point, which is now needed. The curve’s v coordinate ranges from 0 to 1, taking on the value 0.5 at the center of the curve; here, we classify the intersection point, left-parenthesis 0 comma 0 right-parenthesis , with respect to an edge function going through the point on the curve pc and a point along its derivative to determine which side of the center the intersection point is on and in turn how to compute v .

<<Compute v coordinate of curve intersection point>>= 
Float ptCurveDist = std::sqrt(ptCurveDist2); Float edgeFunc = dpcdw.x * -pc.y + pc.x * dpcdw.y; Float v = (edgeFunc > 0) ? 0.5f + ptCurveDist / hitWidth : 0.5f - ptCurveDist / hitWidth;

The partial derivative partial-differential normal p slash partial-differential u comes directly from the derivative of the underlying Bézier curve. The second partial derivative, partial-differential normal p slash partial-differential v , is computed in different ways based on the type of the curve. For ribbons, we have partial-differential normal p slash partial-differential u and the surface normal, and so partial-differential normal p slash partial-differential v must be the vector such that partial-differential normal p slash partial-differential u times partial-differential normal p slash partial-differential v equals bold n and has length equal to the curve’s width.

<<Compute partial-differential normal p slash partial-differential u and partial-differential normal p slash partial-differential v for curve intersection>>= 
Vector3f dpdu, dpdv; EvaluateCubicBezier(pstd::MakeConstSpan(common->cpObj), u, &dpdu); if (common->type == CurveType::Ribbon) dpdv = Normalize(Cross(nHit, dpdu)) * hitWidth; else { <<Compute curve partial-differential normal p slash partial-differential v for flat and cylinder curves>> 
Vector3f dpduPlane = objectFromRay.ApplyInverse(dpdu); Vector3f dpdvPlane = Normalize(Vector3f(-dpduPlane.y, dpduPlane.x, 0)) * hitWidth; if (common->type == CurveType::Cylinder) { <<Rotate dpdvPlane to give cylindrical appearance>> 
Float theta = Lerp(v, -90, 90); Transform rot = Rotate(-theta, dpduPlane); dpdvPlane = rot(dpdvPlane);
} dpdv = objectFromRay(dpdvPlane);
}

For flat and cylinder curves, we transform partial-differential normal p slash partial-differential u to the intersection coordinate system. For flat curves, we know that partial-differential normal p slash partial-differential v lies in the x y plane, is perpendicular to partial-differential normal p slash partial-differential u , and has length equal to hitWidth. We can find the 2D perpendicular vector using the same approach as was used earlier for the perpendicular curve segment boundary edges.

<<Compute curve partial-differential normal p slash partial-differential v for flat and cylinder curves>>= 
Vector3f dpduPlane = objectFromRay.ApplyInverse(dpdu); Vector3f dpdvPlane = Normalize(Vector3f(-dpduPlane.y, dpduPlane.x, 0)) * hitWidth; if (common->type == CurveType::Cylinder) { <<Rotate dpdvPlane to give cylindrical appearance>> 
Float theta = Lerp(v, -90, 90); Transform rot = Rotate(-theta, dpduPlane); dpdvPlane = rot(dpdvPlane);
} dpdv = objectFromRay(dpdvPlane);

The partial-differential normal p slash partial-differential v vector for cylinder curves is rotated around the dpduPlane axis so that its appearance resembles a cylindrical cross-section.

<<Rotate dpdvPlane to give cylindrical appearance>>= 
Float theta = Lerp(v, -90, 90); Transform rot = Rotate(-theta, dpduPlane); dpdvPlane = rot(dpdvPlane);