6.7 曲线 - Physically Based Rendering From Theory To Implementation(PBRT)

Physically Based Rendering From Theory To Implementation

翻译：https://pbr-book.org/4ed/Shapes/Curves

总目录： Physically Based Rendering From Theory To Implementation - 基于物理的渲染从理论到实践第四版原书翻译

6.7 Curves 曲线 ※

虽然三角形或双直线曲面分块（bilinear patches）可以用来表示细条形状，以建模精细的几何形状，如头发、毛皮或草地，但为了更有效地渲染这些对象，有必要使用专门的形状，因为它们通常存在许多单独的实例。本节介绍的曲线形状用三次Bézier曲线建模的薄几何形状表示，这些曲线由四个控制点$p_0$，$p_1$，$p_2$，和$p_3$定义。Bézier样条通过第一个和最后一个控制点。沿着它的点由多项式给出

$$p(u)=(1-u{)}^3{p}_0+3(1-u{)}^{2}u{p}_1+3(1-u)u^2{p}_2+u^3{p}_3$$

(见图6.29)因此，为了在Shape中使用其他基(例如，Hermite样条或b样条)指定的曲线，必须转换为 Bézier 基。

图6.29:一个三次 Bézier 曲线由四个控制点 $p_i$ 定义。曲线 $p(u)$ 如式(6.16)所示，分别经过第一个控制点 $u=0$ 和最后一个控制点 $u=1$ 。

曲线形状(Curve shape)由一条1D Bézier曲线定义，其宽度从起始宽度和结束宽度沿其范围线性插值。它们共同定义了一个平面2D表面(参见图6 - 30)。可以直接将光线与这种表示形式相交，而无需对其进行细分，这反过来使高效地渲染平滑曲线而无需使用太多存储空间成为可能。

图6.30:曲线形状的基本几何形状。在与曲线正交的两个方向上，沿着曲线的每一点，一维 Bézier 曲线被偏移指定宽度的一半。得到的面积表示曲线的表面。

图6.31显示了一个兔子模型，它的皮毛有超过一百万条曲线。

图6.31:毛茸茸的兔子。兔子模型与超过一百万个曲线形状用于建模皮毛。在这里，我们使用了不切实际的长曲线来更好地展示曲线的功能，给出了一个不切实际的毛茸茸的兔子。(底层的兔子网格由斯坦福计算机图形实验室提供。)

class Curve {
public:
    <<Curve Public Methods>> 
private:
    <<Curve Private Methods>> 
    <<Curve Private Members>> 
};

曲线形状可以表示三种曲线类型，如图6.32所示。

• 平坦（Flat）：这种表示形式的曲线总是面向与之相交的光线;他们对精细斜切圆柱形（如头发或毛皮）来建模是十分有用的。
• 圆柱体（Cylinder）：对于在屏幕上跨越几个像素的曲线(比如从不远的地方看到的意大利面)，曲线形状可以计算一个阴影法线，使曲线看起来实际上是一个圆柱体。
• 丝带（Ribbon）:这个变体对于实际上没有圆柱形横截面的形状(如草叶)建模很有用。

图6.32:曲线形状可以表示的三种曲线类型。顶部是一条平坦的曲线，它的方向总是垂直于接近它的射线。中间是这条曲线的变体，其中设置了阴影法线，使曲线看起来是圆柱形的。底部为带状，其起点和终点方向固定;中间方向在它们之间平滑地插值。

CurveType枚举器记录了给定的曲线实例模型。

平面曲线和圆柱曲线的变体是用来方便的近似变形圆柱。值得注意的是，与它们的交叉点并不是物理上可模拟的3D形状，这可能会与在使用真实圆柱体作为参考的场景中导致轻微的不一致。

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

给定在pbrt场景描述文件中指定的曲线，可以将其分割为几个段，每个段覆盖曲线的部分u参数范围。(这样做的一个原因是，轴对齐包围盒不会紧密地束缚抖动曲线，但细分Bézier曲线使它们不那么抖动——多项式的变化递减特性【variation diminishing property】。)因此，Curve构造函数接受一个参数范围的 $u$ 值，$[u_{min},u_{max}]$以及一个指向 CurveCommon 结构的指针，后者存储控制点和有关曲线的其他信息，这些信息在曲线段之间共享。通过这种方式，单个曲线段的内存占用减少了，这使得更容易将更多曲线保存在内存中。

class Curve {
public:
    Curve(const CurveCommon *common, Float uMin, Float uMax)
    	: common(common), uMin(uMin), uMax(uMax) {}
private:
	const CurveCommon *common;
	Float uMin, uMax;
};

CurveCommon构造函数初始化成员变量，并传递给它控制点、曲线宽度等值。提供给它的控制点应该在曲线的对象空间中。

对于带状【Ribbon】曲线，CurveCommon存储一个曲面法线，用于在每个端点定位曲线。构造函数预先计算了两个法向量之间的夹角和这个角的正弦值；这些值在计算沿其范围任意点的曲线方向时非常有用。

<<CurveCommon Definition>>= 
struct CurveCommon {
       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;

       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 包围曲线

要确定曲线在对象空间中的包围盒，可以先将样条曲线沿曲线的中心包围，然后将该包围盒扩展为曲线在其范围内最大宽度的一半。然后，Bounds()方法将该绑定转换为渲染空间，然后返回它。

class Curve {
public:
    Curve(const CurveCommon *common, Float uMin, Float uMax)
    	: common(common), uMin(uMin), uMax(uMax) {}

	Bounds3f Curve::Bounds() const {
	    pstd::span<const Point3f> cpSpan(common->cpObj);
	    Bounds3f objBounds = BoundCubicBezier(cpSpan, uMin, uMax);
	    // <<Expand objBounds by maximum curve width over  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);
	}

	DirectionCone NormalBounds() const { return DirectionCone::EntireSphere();} 
	
private:
	const CurveCommon *common;
	Float uMin, uMax;
};

曲线形状不能用作面积光，因为它不提供所需采样方法的实现。它确实提供了一个NormalBounds()方法来返回一个保守的边界。

6.7.2 Intersection Tests 交叉测试

Shape接口需要的两个求交方法都是通过曲线的另一个方法IntersectRay()实现的。它接收一个指向ShapeIntersection的指针，而不是返回一个可选的ShapeIntersection。

<<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()将nullptr传递给IntersectRay()，这表明如果找到交集，它可以立即返回。

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

曲线相交算法是基于丢弃曲线段，只要它可以确定射线肯定不会与它们相交，否则递归地将曲线一分为二，以创建两个更小的段，然后进行测试。最后，曲线被线性逼近以进行有效的相交测试。这个过程需要在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>> 
    <<Project curve control points to plane perpendicular to ray>> 
    <<Test ray against bound of projected control points>> 
    <<Compute refinement depth for curve, maxDepth>> 
    <<Recursively test for ray–curve intersection>> 
}

CurveCommon 类存储整个曲线的控制点，但是Curve实例通常需要表示Bézier曲线u范围上的四个控制点。实用函数CubicBezierControlPoints()完成了这项计算。

<<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>> 
    <<Test ray against bound of projected control points>> 
    <<Compute refinement depth for curve, maxDepth>> 
    <<Recursively test for ray–curve intersection>> 
}

与6.5.3节中的射线-三角形求交算法类似，射线-曲线求交的方法是将曲线转换为一个坐标系，射线的原点在该坐标系的原点，射线的方向沿+z坐标轴对齐。在开始时执行这种转换会大大减少了交叉测试必须执行的操作数量。

对于曲线形状，我们需要显式地表示转换，因此这里使用LookAt()函数来生成它。原点是射线的原点，“看”点是沿着射线方向与原点偏移的点。“向上”方向被设置为垂直于两个方向，光线的方向和从第一个控制点到最后一个控制点的方向。这样做有助于使曲线的方向大致平行于射线坐标系中的x坐标轴，从而使y方向边界更紧密(参见图6 - 33)。这种边界匹配方面的改进通常可以使递归交集测试提前结束。

图6.33：Bézier 曲线的二维包围盒。(a)根据给定的曲线控制点计算包围盒。(b)旋转曲线的效果，使矢量从第一个控制点到最后一个控制点在计算边界之前与轴对齐。得到的边界框更加紧密。

如果第一个控制点和最后一个控制点之间的向量和光照射线平行，则dx是退化的。在这种情况下，我们找到一个任意的“向上”矢量方向，以便在这种不寻常的情况下进行相交测试。

<<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>> 
    <<Compute refinement depth for curve, maxDepth>> 
    <<Recursively test for ray–curve intersection>> 
}

沿着在Curve::Bounds()中实现的路线，可以通过取曲线控制点的边界并在考虑的范围$u$内扩展曲线最大宽度的一半来找到曲线段的保守包围盒。

图6.34：射线-曲线包围测试。在射线坐标系中，射线的原点在（0,0,0），方向与+z 轴对齐。因此，如果二维点$(x,y)=(0,0)$在曲线段的xy包围盒之外，则射线不可能与曲线相交。

由于光线的原点在$(0,0,0)$，方向与交点空间中的 $+z$ 轴对齐，因此其包围盒只包含在x和y轴的原点(图6.34)；其z范围由其参数z范围所涵盖的范围给出。在进行递归相交测试算法之前，先对射线的包围盒与曲线的包围盒进行相交测试。如果它们不相交，该方法可以立即返回。

<<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>> 
    <<Recursively test for ray–curve intersection>> 
}

计算细分曲线的最大次数，以便在最精细的层次上，最终线性化曲线的最大距离被限定为小于一个小的固定距离。我们不会深入这个计算的细节，它是在片段<<Compute refine depth for curve, maxDepth>>中实现的。通过筛选测试并获得该值后，开始递归交集测试。

<<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>> 
    // <<Recursively test for ray–curve intersection>> 
	pstd::span<const Point3f> cpSpan(cp);
	return RecursiveIntersect(ray, tMax, cpSpan, Inverse(rayFromObject),
	                          uMin, uMax, maxDepth, si);
}

然后，RecursiveIntersect()方法测试给定的光线是否与给定的参数范围$[u0,u1]$内的曲线段相交。它假设射线已经针对曲线的包围盒进行了测试，并发现它与曲线的包围盒相交。

<<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>> 
    } 
    else 
    {
        <<Intersect ray with curve segment>> 
    }
}

如果还没有达到最大深度，则调用SubdivideCubicBezier()将给出Bézier曲线的控制点，从而将cp给出的Bézier曲线一分为二。第一条曲线的最后一个控制点与第二条曲线的第一个控制点相同，因此总共8个控制点返回了7个值。然后初始化u数组，以便在依次处理两条曲线之前保存两条曲线的参数范围。

<<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>> 
		    <<Recursively test ray-segment intersection>> 
		}
		return si ? si->has_value() : false;
    } 
    else 
    {
        <<Intersect ray with curve segment>> 
    }
}

在<<Check ray against curve segment’s bounding box>>片段的包围盒测试在本质上与<<Test ray against bound of projected control points>>一样，除了它在计算曲线在范围u内的最大宽度时从u数组中获取u值，并且它使用来自cpSplit的控制点。因此，这里不包括它。

如果射线确实与包围盒相交，则相应的段将被递归调用RecursiveIntersect()。如果找到了交点，并且光线是阴影光线，则si将为nullptr，并且可以立即报告交点。对于非阴影光线，即使找到了一个交点，它也可能不是最近的交点，因此仍然必须考虑另一段。

<<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>> 
		    
		    // <<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>> 
    }
}

通过使用曲线的线性近似，使相交检验更有效；变化递减特性使我们可以在不引入太多误差的情况下进行近似。

<<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>> 
		    
		    // <<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>> 
		<<Find line  that gives minimum distance to sample point>> 
		<<Compute  coordinate of curve intersection point and hitWidth>> 
		<<Test intersection point against curve width>> 
		if (si) {
		    <<Initialize ShapeIntersection for curve intersection>> 
		}
		return true;
    }
}

图6.35：曲线分段边界。较大曲线的段的相交测试计算垂直于段端点的线(虚线)的边缘函数。如果一个潜在的交点(实心点)在边缘的另一边而不是线段，则拒绝；另一个曲线段(如果在那一边)应该代替这个交集。

重要的是，交点测试只接受当前正在考虑的u段的Curve表面上的交点。因此，交点测试的第一步是计算垂直于曲线起点和终点的直线的边缘函数，并根据它们对潜在交点进行分类(图6.35)。

<<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) 
    {
        ...
    } 
    else 
    {
        // <<Intersect ray with curve segment>> 
        {
			// <<Test ray against segment endpoint boundaries>> 
			{
				<<Test sample point against tangent perpendicular at curve start>> 
				<<Test sample point against tangent perpendicular at curve end>> 
			}
			<<Find line  that gives minimum distance to sample point>> 
			<<Compute  coordinate of curve intersection point and hitWidth>> 
			<<Test intersection point against curve width>> 
			if (si) {
			    <<Initialize ShapeIntersection for curve intersection>> 
			}
		}
		return true;
    }
}

将曲线控制点投影到射线坐标系中可以使这个测试更高效，原因有二。首先，因为射线的方向是与+z轴相关的，所以问题被简化为在x和y中的2D测试。其次，因为射线原点在坐标系的原点，所以我们需要分类的点是(0,0)，这简化了对边缘函数的评估，就像射线-三角形相交测试一样。

在式(6.5)中引入边缘函数进行射线-三角形相交测试;参见图6.14。为了定义这个边缘函数，我们需要直线上任意两点垂直于经过起点的曲线。第一个控制点，$p_0$，对于第一个控制点来说是一个很好的选择。对于第二个，我们将计算垂直于曲线切线的向量，并将该偏移量添加到控制点。

对式(6.16)求导可知，曲线在第一个控制点$p_0$处的切线为$3(p_1-p_0)$。缩放因子在这里并不重要，所以我们将在这里使用$t=p_1-p_0$。在2D中计算垂直于切线的向量很简单：只需要交换x和y坐标，并使其中一个向量的值为负即可。(为了了解为什么这样做，考虑点积$(x,y)\cdot (y,-x)=xy+-yx=0$。因为这两个向量夹角的余弦为0，所以它们一定是垂直的。)因此，边缘上的第二个点是

$$\begin{gathered} p_0+(p_{1_y}-p_{0_y},-(p_{1_x}-p_{0_x})) \\ =p_0+(p_{1_y}-p_{0_y},p_{0_x}-p_{1_x})) \end{gathered}$$

将这两点代入边缘函数的定义式(6.5)，化简得

$$e(p)=(p_{1_y}-p_{0_y})(p_y-p_{0_y})-(p_x-p_{0_x})(p_{0_x}-p_{1_x})$$

最后，代入 $p=(0,0)$ 得到要测试的最终表达式:

$$e((0,0))=(p_{1_y}-p_{0_y})(-p_{0_y})+p_{0_x}(p_{0_x}-p_{1_x})$$

<<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) 
    {
        ...
    } 
    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>> 
			}
			<<Find line  that gives minimum distance to sample point>> 
			<<Compute  coordinate of curve intersection point and hitWidth>> 
			<<Test intersection point against curve width>> 
			if (si) {
			    <<Initialize ShapeIntersection for curve intersection>> 
			}
		}
		return true;
    }
}

<<Test sample point against tangent perpendicular at curve end>>片段同理，在曲线的末端进行相应的测试。

<<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) 
    {
        ...
    } 
    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  that gives minimum distance to sample point>> 
			<<Compute  coordinate of curve intersection point and hitWidth>> 
			<<Test intersection point against curve width>> 
			if (si) {
			    <<Initialize ShapeIntersection for curve intersection>> 
			}
		}
		return true;
    }
}

测试的下一部分是确定曲线段上最接近点(0,0)的曲线上的u值。如果该点离中心点的距离不超过曲线的宽度，这个点就是交点。为一个三次Bézier曲线确定这个距离需要大量的计算，因此这里的实现用一个线段近似曲线来计算这个u值。

图6 - 36：用线段逼近三次Bézier曲线。对于光线曲线相交测试的这一部分，我们用通过其起点和终点的线段(虚线)近似Bézier。(实际上，经过细分后，曲线已经接近线性，因此误差小于图中所示的误差。)

我们用从起点 $p_0$ 到终点 $p_3$ 的线段线性逼近Bézier曲线，参数为$w$。在这个例子中，位置为 $w=0$ 的 $p_0$ 和 $w=1$ 的 $p_3$ (如图6 - 36所示)。我们的任务是计算与 $p$ 点最接近的直线上的点 $p^{\prime }$ 对应的直线上的值 $w$。这里的关键是$p^{\prime }$，即从直线上对应的点到点 $p$ 的向量垂直于直线(如图6 - 37(a)所示)。

图6.37:(a)给定一条无限直线和一个点$p$，从该点到直线上最近点$p^{\prime }$的向量垂直于直线。(b)因为这个向量是垂直的，所以我们可以计算直线的第一个点到最接近点 $p^{\prime }$ 的距离为$d=\Vert p-p_0\Vert cos\theta$。

公式(3.1)给出了两个向量的点积、它们的长度和它们之间夹角的余弦值之间的关系。特别是，它向我们展示了如何计算$p_0$到$p$向量和$p_0$到$p_3$的向量之间夹角的余弦值:

$$cos\theta =\frac{(p-p_0)\cdot (p_3-p_0)}{\Vert p-p_0\Vert \Vert p_3-p_0\Vert }$$

因为从$p^{\prime }$到$p$的向量垂直于直线(图6.37(b))，我们可以计算沿直线从$p_0$到$p^{\prime }$的距离为

$$\begin{gathered} d=\Vert p-p_0\Vert cos\theta =\Vert p-p_0\Vert \frac{(p-p_0)\cdot (p_3-p_0)}{\Vert p-p_0\Vert \Vert p_3-p_0\Vert } \\ =\frac{(p-p_0)\cdot (p_3-p_0)}{\Vert p_3-p_0\Vert } \end{gathered}$$

最后，沿直线的参数偏移量$w$是$d$与直线长度的比值

$$w=\frac{d}{\Vert p_3-p_0\Vert }=\frac{(p-p_0)\cdot (p_3-p_0)}{\Vert p_3-p_0{\Vert }^2}$$

$w$ 值的计算反过来又稍微简化了相交坐标系中 $p=(0,0)$ 的事实。

<<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) 
    {
        ...
    } 
    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  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  coordinate of curve intersection point and hitWidth>> 
			<<Test intersection point against curve width>> 
			if (si) {
			    <<Initialize ShapeIntersection for curve intersection>> 
			}
		}
		return true;
    }
}

Bézier曲线上(假定的)与候选交点最近的点的参数u坐标是通过沿着线段的范围进行线性插值计算的。有了这个值，就可以计算出曲线在该点处的宽度。

<<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) 
    {
        ...
    } 
    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  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  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>> 
			}
			
			<<Test intersection point against curve width>> 
			if (si) {
			    <<Initialize ShapeIntersection for curve intersection>> 
			}
		}
		return true;
    }
}

对于Ribbon曲线，曲线并不总是面向光线。相反，它的方向是在每个端点给出的两个表面法线之间插值。本文采用球面线性插值对法向量进行插值。然后，曲线的宽度被归一化的射线方向和带状方向之间的夹角的余弦缩放，以便它对应于从给定方向可见的曲线宽度。

<<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) 
    {
        ...
    } 
    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  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  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>> 
			if (si) {
			    <<Initialize ShapeIntersection for curve intersection>> 
			}
		}
		return true;
    }
}

为了最终将潜在的交点分类为命中或未命中，Bézier曲线仍然必须在u内评估。(由于控制点cp表示当前考虑的曲线段，因此重要的是在函数调用中使用$w$而不是在函数调用$u$中，因为是在范围 $(0,1)$ 内。)在这一点上，曲线的导数很快就会有用，所以现在记录下来。

我们想在曲线pc上测试$p$到这一点的距离是否小于曲线宽度的一半。因为$p=(0,0)$，我们可以等效地测试pc到原点的距离是否小于宽度的一半，或者距离的平方是否小于宽度的平方的1 / 4。如果这个测试通过了，最后要检查交点是否在光线的参数 $t$ 范围内。

<<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) 
    {
        ...
    } 
    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  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  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>> 
			}
		}
		return true;
    }
}

对于非阴影光线，最终交点的ShapeIntersection可以被初始化。这样做需要计算交点的光线 $t$ 值以及它的SurfaceInteraction。

<<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) 
    {
        ...
    } 
    else 
    {
        ...

		if (si) {
		    // <<Initialize ShapeIntersection for curve intersection>> 
		    <<Compute tHit for curve intersection>> 
			<<Initialize SurfaceInteraction intr for curve intersection>> 
			*si = ShapeIntersection{intr, tHit};
		}
		
		return true;
    }
}

计算出tHit值后，将其与先前发现的光线曲线相交点(如果有)的tHit进行比较。这个检查确保返回最近的交集。

<<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) 
    {
        ...
    } 
    else 
    {
        ...

		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>> 
			*si = ShapeIntersection{intr, tHit};
		}
		
		return true;
    }
}

为了能够初始化交集的SurfaceInteraction，需要计算各种额外的量。

<<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) {...} 
    else {
        ...

		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>> 
			<<Compute $\dpdu$ and $\dpdv$ for curve intersection>> 
			<<Compute error bounds for curve intersection>> 
			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;
    }
}

到目前为止，我们没有计算交点的 $v$ 坐标，现在需要计算交点的 $v$ 坐标。曲线的 $v$ 坐标范围为 0 ~ 1，取曲线中心的值 0.5；在这里，我们根据经过曲线PC上的点和沿其导数的点的边缘函数对交点进行分类，以确定交点在中心的哪一边，以及如何计算 $v$ 。

<<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) {...} 
    else {
        ...

		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 $\dpdu$ and $\dpdv$ for curve intersection>> 
			<<Compute error bounds for curve intersection>> 
			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;
    }
}

偏导数$\frac{\partial p}{\partial u}$直接来自底层 Bézier 曲线的导数。第二个偏导数$\frac{\partial p}{\partial v}$，根据曲线的类型有不同的计算方法。
For ribbons, we have θp/θu and the surface normal, and so θp/θv must be the vector such that θp/θu×θp/θυ = n and has length equal to the curve’s width.
对于带状，我们有$\frac{\partial p}{\partial u}$和表面法线，所以$\frac{\partial p}{\partial v}$必须是矢量，使得$\frac{\partial p}{\partial u}\times \frac{\partial p}{\partial v}=n$，并且长度等于曲线的宽度。

<<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) {...} 
    else {
        ...

		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 $\dpdu$ and $\dpdv$ 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  for flat and cylinder curves>> 
			}
			<<Compute error bounds for curve intersection>> 
			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;
    }
}

对于平面曲线和圆柱曲线，我们转换$\frac{\partial p}{\partial u}$为交点坐标系。对于平坦曲线，我们知道$\frac{\partial p}{\partial v}$在 $xy$ 平面上，垂直于$\frac{\partial p}{\partial u}$，长度等于hitWidth。我们可以使用与之前用于垂直曲线段边界边缘相同的方法来找到二维垂直向量。

<<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) {...} 
    else {
        ...

		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 dp/du and dp/dv 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 dp/dv 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>> 
				}
				dpdv = objectFromRay(dpdvPlane); 
			}
			<<Compute error bounds for curve intersection>> 
			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;
    }
}

圆柱曲线的$\partial p/\partial v$矢量绕dpduPlane轴旋转，使其外观类似于圆柱形截面。

<<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) {...} 
    else {
        ...

		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 dp/du and dp/dv 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 dp/dv 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>> 
			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;
    }
}

最终整体代码：

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
            // FIXME: this tHit isn't quite right for ribbons...
            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 $\dpdu$ and $\dpdv$ for curve intersection
            Vector3f dpdu, dpdv;
            EvaluateCubicBezier(pstd::MakeConstSpan(common->cpObj), u, &dpdu);
            CHECK_NE(Vector3f(0, 0, 0), dpdu);
            if (common->type == CurveType::Ribbon)
                dpdv = Normalize(Cross(nHit, dpdu)) * hitWidth;
            else {
                // Compute curve $\dpdv$ 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};
        }
#ifndef PBRT_IS_GPU_CODE
        ++nCurveHits;
#endif
        return true;
    }
}