Past exam of the mathematics course of the University of Cambridge / 2008 / iii / Paper 77 (Computer Aided Geometric Design)

Let the triangles have vertices $A_0,A_1,A_2$ and $B_0,B_1,B_2$. Work initially with nonzero-area triangles. A triangle-triangle intersection algorithm begins with the unnormalized normal vector $N_A=(A_1-A_0)\times (A_2-A_0)$ and the three signed plane evaluations $d_j=N_A\cdot (B_j-A_0)$. If all $d_j$ are strictly positive, or all strictly negative, the second triangle lies strictly on one side of the first one's plane, so the triangles are disjoint. Otherwise compute $N_B$ and $e_i=N_B\cdot (A_i-B_0)$ and apply the symmetric rejection test. Strict inequalities are important: contact on a vertex or edge counts as intersection.

If the planes are not parallel, put $L=N_A\times N_B$. Each triangle meets the other one's plane in an empty set, a point or a line segment. Its endpoints are obtained from zero-valued vertices and the edges whose endpoint plane evaluations have opposite signs. For example, an edge $B_j{B}_k$ contributesAll these points lie on the common line of the two planes. Choose a coordinate $r$ for which $L_r\ne 0$ and form the coordinate intervals $I_A=[a_{-},a_{+}]$, $I_B=[b_{-},b_{+}]$ of the two slices. This coordinate is one-to-one along the common line, so the exact final test isOnly that coordinate of each slice endpoint needs to be evaluated; no full three-dimensional line parameterization is necessary.

Parallel planes which have survived the strict plane tests are coplanar. Drop the coordinate corresponding to the largest component of $N_A$ and test the two projected planar triangles. For each edge of either projected triangle, project both vertex sets onto the perpendicular direction. A strictly disjoint pair of projection intervals certifies separation. If none of the six edge directions separates them, the planar convex sets overlap: a separating line for two disjoint convex polygons can be chosen parallel to an edge of one of them. This includes containment and boundary contact.

A zero normal vector signifies a degenerate triangle; reduce its convex hull to a point or its extreme line segment and use point containment or segment intersection instead. For floating-point input, reliable sign predicates and scale-aware error bounds are preferable to treating a tiny signed distance as an arbitrary exact zero. Plane rejection followed by interval overlap, with a planar test for coplanar input, determines all nondegenerate cases.

Count scalar additions, subtractions, multiplications and divisions for the specific triangle-triangle intersection algorithm above; comparisons and index selection are separate. Computing one plane normal vector uses two vector differences, costing six subtractions, and one cross product, costing six multiplications and three subtractions. The three vertex-to-plane evaluations use nine further coordinate subtractions and three dot products, each with three multiplications and two additions. The first rejection test therefore costsFor randomly located triangles of characteristic size $a$ in a box of size $H\gg a$, a second triangle straddles the first triangle's plane only when its position is in a slab of relative thickness $O(a/H)$. Consequently the first plane test rejects almost every pair: the expected arithmetic count for this implementation is $39+O(a/H)$, with the constant in the remainder determined by the subsequent branches. This is an average count, rather than the assertion that every input needs 39 operations.

For a generic surviving noncoplanar pair, the symmetric plane test costs another 39 operations, the cross product $N_A\times N_B$ costs nine, and each of the four slice endpoints needs five scalar operations in the selected coordinate: one subtraction and division for the interpolation fraction, then one subtraction, multiplication and addition for the coordinate. Thus this branch usesA fully executed coplanar test can use six planar edge directions. Each needs two subtractions to form its perpendicular and six two-dimensional dot products, costing $2+6(3)=20$ operations. Including both plane tests and the parallelism cross product gives an upper count $78+9+120=207$ before any early exit. Contact and degeneracy can shorten or change these counts. Precomputed normal vectors or a preliminary axis-aligned bounding box test also change the answer: the latter uses only comparisons after its coordinate extrema are available. The relevant small-triangle average is 39 operations for the specified plane-first implementation; a generic complete test uses 107.

Store a bounding volume hierarchy for each object, using axis-aligned bounding boxes or tighter enclosing volumes. A leaf contains one triangle or a small group; every internal node encloses all its descendants. Start with the two roots. If their bounding volumes are disjoint, discard that pair. Otherwise descend into children of one or both nodes, usually splitting the larger node, until overlapping leaf pairs are reached. Apply the exact triangle-triangle intersection algorithm only at these leaves, and stop at the first hit if only a yes/no answer is required. The pruning is sound because disjoint enclosing volumes imply disjoint enclosed triangles.

A balanced bounding volume hierarchy typically takes $O(n\log n)$ work to construct and $O(n)$ storage. It can be reused for repeated queries and refitted when vertex positions move; a rigid transformation can also be handled in a common coordinate frame. Spatial grids or octrees provide another broad-phase strategy: place triangles in every cell they meet and test only pairs from common cells, suppressing duplicates. A box test is especially valuable because the small randomly distributed triangles almost always fail it.

For $n,m$ of order $10^5$, the naive method has about $10^{10}$ candidate pairs. Hierarchical pruning replaces this by the number of overlapping node pairs and actual candidate leaf pairs. No universal subquadratic worst-case bound follows: interpenetrating or very poorly separated meshes can still generate $nm$ candidates. If “intersection of objects” means overlap of solid volumes rather than intersection of their boundaries, add a point-in-solid containment test when no surface intersections are found; one closed object can lie wholly inside the other. Use spatial enclosure to prune pairs before exact triangle tests.

At $X=P(u,v)$, orient the unit normal vector $n=(P_u\times P_v)/\Vert P_u\times P_v\Vert$ towards the eye. Define the unit directions from the reflecting point to the source and eye byThe incident propagation direction is $-a$. The law of specular reflection is $b=-a+2(a\cdot n)n$, so $a+b$ is normal to the surface. Conversely, if the tangential components of $a+b$ vanish and $a\cdot n,b\cdot n>0$, its unit-length summands have equal positive normal components and obey the reflection law. Therefore solve the two scalar equationsThe stated orientation condition gives $b\cdot n>0$. It removes the back-facing ambiguity, but not multiple images or caustics. Exclude $Q(t)=P(u,v)$ and $E=P(u,v)$, where these directions are undefined.

At a root where the $2\times 3$ Jacobian matrix $DF$ has rank two, the implicit function theorem makes the solution set a local curve in parameter space. Compute a unit tangent $w$ in its nullspace; equivalently take the normalized cross product of the two rows of $DF$. A continuation method predicts $x_p=x+hw$, where $x=(u,v,t)$, then corrects with the Newton method applied to the three equationsThe last equation fixes a transverse section and allows tracing through a turning point of $t$, where solving only for $(u,v)$ at successive fixed source parameters would fail. The corresponding space point is $R=P(u,v)$, with derivative $R^{\prime }=P_u{u}^{\prime }+P_v{v}^{\prime }$ along the traced parameter-space curve.

Find initial roots by subdividing the parameter cube, using interval arithmetic bounds on $F$ to discard impossible boxes and validated local solves in the remaining boxes. Trace in both directions, stop at the boundary or a recognized closed loop, and retain a record of traced boxes to avoid duplicates. Search the remaining boxes as well: starting from just one seed would miss disconnected closed reflection curves. Where $DF$ loses rank, reduce the step and subdivide locally to find the outgoing branches instead of continuing a single guessed tangent. If physical visibility is required, discard roots for which either open segment $XE$ or $XQ(t)$ is occluded; test these segments against the surface or any supplied scene geometry.

Control the geometric approximation rather than just the parameter step. For any regular traced arc $R(s)$ with a bound $\Vert R^{\prime \prime }(s)\Vert \le M$, the chord error bound gives an arc-to-chord Hausdorff distance at most $M(\Delta s{)}^2/8$. Budget, for example, half of $ϵ$ for corrected-vertex error and half for chord error; accept a segment only whenBound the vertex error through the inverse transverse Jacobian matrix and a bound on $DP$, or use validated root enclosures. If analytic second-derivative bounds are unavailable, recursively subdivide and enclose each arc in a geometric tube of radius at most the remaining tolerance around its chord. A midpoint test alone is a useful heuristic but cannot certify the requested error for an arbitrary oscillatory curve.

The returned vertices $P(u_k,v_k)$ and their ordered connecting segments give the reflection curve to the prescribed tolerance on every resolved branch. Finite certified termination requires regular evaluable patches, derivative or enclosure bounds, and finitely resolvable branches; eye-facing normals alone do not guarantee these properties. Singular components or an identically vanishing reflection constraint must be handled as such rather than silently reported as an ordinary curve.

Write a degree-$p$ B-spline as $C(t)={\sum }_i{P}_i{N}_{i,p}(t)$ and locate its spline knot span. There are several equivalent exact evaluation methods.

Evaluate the at most $p+1$ nonzero basis functions by the Cox-de Boor recurrence, then form their weighted sum of control points. Alternatively apply De Boor's algorithm directly to the $p+1$ active control points. In a span $[t_k,t_{k+1})$, initialize $d_j^{(0)}=P_{k-p+j}$ and repeatedly interpolate:for $r=1,\dots ,p$ and $j=p,p-1,\dots ,r$. The last point $d_p^{(p)}$ is $C(t)$; zero knot intervals are handled by the standard limiting convention or an appropriate neighboring nonempty span.

A third method precomputes the local polynomial coefficients on each span and evaluates them by Horner's method. For a uniform cubic B-spline, a fixed four-by-four basis matrix gives these coefficients. A fourth method uses knot insertion to express the relevant span as a Bézier curve, followed by De Casteljau's algorithm. Both interpolation algorithms also split the curve at the requested parameter. Repeated subdivision supplies arbitrarily accurate polygonal evaluation, although a finite subdivision mesh is generally an approximation rather than the exact requested point. Basis summation, de Boor interpolation, local polynomial evaluation and Bézier extraction are equivalent exact representations.

The Cox-de Boor recurrence and De Boor's algorithm handle nonuniform knots and repeated knots without converting the whole curve. With a tabulated recurrence they require $O(p^2)$ arithmetic after locating the span and use only local control points; a naive recursive basis evaluation wastes work by recomputing the same terms. De Boor's algorithm evaluates the vector directly using affine combinations. Inside the span these are convex combinations, which makes the geometric construction transparent and avoids the cancellation that can occur in a power-basis representation. Basis evaluation is especially useful when the same knots and parameter are shared by many coordinate functions or many curves, since the weights can be reused.

Precomputed power coefficients followed by Horner's method need only $O(p)$ arithmetic per coordinate and are attractive for many evaluations within an unchanged span, particularly for fixed-degree uniform B-splines. Coefficient construction has an initial cost, control-point edits require updates to affected spans, and poor scaling of high-degree power coefficients can damage accuracy. Rescaling each span to $[0,1]$ helps. Bézier curve extraction has a similar initial conversion cost; De Casteljau's algorithm then gives stable local evaluation, subdivision, a convex hull bound and geometric flatness tests useful for rendering and intersection.

Subdivision is convenient when an entire polygonal approximation is wanted: refine only where needed and reuse the hierarchy. It is less economical for a single prescribed parameter and requires a geometric error bound before its approximate point can replace exact evaluation. Span lookup itself costs $O(\log n)$ for sorted arbitrary knots, constant work for a uniform interior span, or amortized constant work for an ordered sequence of sample parameters. Use de Boor for robust general knots, shared basis weights for batches, and precomputed local polynomials or Bézier pieces for repeated interrogation.

Use unit knot spacing and the centered Cardinal cubic B-spline convention $C(t)={\sum }_i{P}_i{M}_3(t-i)$, where $M_3$ is supported on $[-2,2]$. On the span $[i,i+1]$, with $s=t-i$, the active control points are $P_{i-1},P_i,P_{i+1},P_{i+2}$, and their B-spline weights areFor $t=11/4$, take $i=2$, $s=3/4$. The weights and their derivatives areThe PDF has $P_3=(0,0,0)$ and $P_4=(0,1,0)$, so the point and tangent vector areThus, under this centered convention,The original PDF specifies neither the knot origin nor its indexing relative to $P_i$. If instead the convention associates span $[i,i+1]$ with $P_i,P_{i+1},P_{i+2},P_{i+3}$, the identical local weights multiply $P_2,P_3,P_4,P_5$ at the requested parameter, givingThese are parameter translations of the same uniform construction, rather than contradictory evaluations with the same knots. A specified knot convention is necessary to select the numerical pair uniquely; changing the knot spacing would also rescale the derivative.

A parallel surface at signed distance $d$ from a design surface describes a constant-thickness skin, shell, coating or clearance envelope. A particularly important application is machining with a spherical cutter of radius $r$: when the cutter touches $P$ without penetration, its center lies at $P+rn$, where $n$ is the outward unit normal vector. Thus a tool-center path can be computed on the surface offset rather than on the desired contact surface. Positive and negative surface offsets also model material addition and removal, and pairs of offsets describe the two faces of a thin manufactured wall. Offsets convert contact, thickness and clearance specifications into geometric surfaces which can be interrogated.

The construction is initially local. Where the offset develops self-intersections or reaches a curvature radius, not every point of the raw parallel surface is a usable physical boundary. Tool paths and solid envelopes may require trimming those pieces and checking global collisions.

For a regular parametric surface $P(u,v)$, form $N=P_u\times P_v$, $n=N/\Vert N\Vert$ and define the offset evaluatorIt is already a parametric surface with the same parameter domain, even if its coordinates no longer belong to the original polynomial or rational patch family. Point interrogation requires evaluations of $P,P_u,P_v$ and one normalization. For tangent or normal vector interrogation, differentiate the evaluator. For example,and similarly for $v$. The offset normal vector is the normalized cross product $(P_d{)}_u\times (P_d{)}_v$ wherever it is nonzero. These derivatives also support intersection, projection and curvature calculations by the usual methods for parametric surfaces. For a variable normal offset, add $d_{u}n$ and $d_{v}n$ to the corresponding tangent derivatives.

With shape operator $W=-Dn$, the constant-distance offset differential is $D{P}_d=(I-dW)DP$. Along a principal direction its factor is $1-d{\kappa }_i$, where ${\kappa }_i$ is a principal curvature. Consequently regular interrogation fails at $d{\kappa }_i=1$; on regular pieces the offset principal curvatures are ${\kappa }_i/(1-d{\kappa }_i)$ with the continued normal orientation. Global self-intersections can occur even when both local factors are nonzero. Evaluate and differentiate the base surface and its normalized normal; an exact conversion back to the original patch type is unnecessary.

Treat the subdivision surface as an evaluable limit together with its refinement hierarchy. At a requested location $\xi$, evaluate the limit position $S(\xi )$ and two independent limit tangent vectors, then setOn regular regions, local spline patches supply position and derivative evaluations. Near extraordinary vertices, use the actual subdivision rule and its limit-position and tangent evaluation machinery, or a convergent local refinement with controlled position and normal vector errors. If a tangent pair degenerates, the regular offset evaluator is not defined there; a crease requires separate one-sided offsets and an explicitly chosen join.

For polygonal interrogation, subdivide the base mesh adaptively. At every new sample, evaluate the base limit and its limit normal vector, and displace that sample by $dn$. Connect these samples with the inherited connectivity, refining until flatness and normal-variation bounds meet the desired tolerance. If the position error is $\delta$ and unit-normal error is $\eta$, the offset-position error is at most $\delta +|d|\eta$, in addition to the interpolation error between samples. A variable normal offset can be sampled at the same locations, with its approximation error included as well. The result retains the base subdivision hierarchy for locating points and refining patches.

Displacing the original control vertices along estimated vertex normals and applying the original subdivision mask just once does not, in general, give the exact surface offset. If $A$ is one linear refinement step, then $A(P+dn(P))=AP+dAn(P)$, while displacing the refined base gives $AP+dn(AP)$; normalization and the normal calculation are nonlinear, so $An(P)\ne n(AP)$ in general. Thus the exact offset is an interrogation wrapper around the base limit, and a conventional subdivision control mesh for it is a fitted approximation requiring refinement or error control. Refine and interrogate the base limit, then offset its evaluated points; do not assume offsetting commutes with subdivision.

A practical family consists of constant-distance parallel surfaces and smooth variable normal offsetsConstants cover wall thickness, cutter radius and clearance. A scalar B-spline or piecewise polynomial displacement field permits thickness variation, local bumps and smooth blending; compact support makes local edits economical. Use the base parameterization or compatible data on its subdivision hierarchy so that evaluating the displacement and its derivatives requires only a small local stencil. Choose displacement size and gradients to avoid unintended folds, and trim global self-intersections where the physical application requires an envelope. A constant displacement smaller than the local curvature radii prevents the corresponding local focal singularities, although global clearance still needs a separate check.

When remaining within a polynomial or rational representation is more useful than prescribing an exact distance, a related practical form is $P+\lambda (u,v)(P_u\times P_v)$ with a low-degree scalar field $\lambda$. Its displacement is normal but has magnitude $|\lambda |\Vert P_u\times P_v\Vert$, so it is not a constant-distance parallel surface. For polynomial $P$ and $\lambda$, it stays polynomial; for rational input it stays rational. More general low-degree vector displacement fields similarly preserve efficient evaluation, while relinquishing the normal-offset condition. These forms make explicit the tradeoff between an exact metric offset, which involves normalization, and an easily represented deformation. Simple scalar normal-displacement fields are useful and directly evaluable; algebraic displacement fields are useful when representation closure is required.

The supplied subdivision mask describes the regular, valence-six triangular grid. Its coefficients split into four parity classes after binary refinement: one class for old vertices and three classes for the three edge directions. For an old vertex $V$ with cyclic neighbors $V_1,\dots ,V_6$, the even stencil isFor an edge with endpoints $A,B$ and opposite vertices $C,D$ in its two incident triangles, the new odd vertex isRotating this edge stencil supplies the other two odd parity classes. All stencil weights sum to one, giving affine invariance; their first moments place an affine regular grid at the old positions and edge midpoints. Apply the same weights separately to $x,y,z$. Boundary edges and extraordinary vertices need additional rules, which are not specified by this uniform mask. The subsequent arguments concern the regular interior grid and compatible extruded boundary treatment.

Choose one grid-edge vector $a$ as the extrusion direction and an independent edge vector $b$. Label grid points $r_{i,j}=ia+jb$, so extruded height data have $z_{i,j}=g_j$, independent of $i$. The Loop subdivision stencils are affine-exact in the $x,y$ coordinates, so the refined horizontal grid has spacing one-half and the same row structure.

At an old vertex, two neighbors lie in its own row, two in the row above and two in the row below. Its refined height is thereforeA new vertex on an edge parallel to $a$ has both endpoints in row $j$ and its opposite vertices in rows $j-1,j+1$, giving exactly the same value $g_{2j}^{\prime }$. A new vertex on either of the other edge directions has endpoints in rows $j,j+1$; its two opposite vertices are also one in each of these rows. Its height isEvery refined vertex in a given row consequently has the same height. Induction proves this at every refinement level. The piecewise planar interpolating mesh is also constant along $a$ between rows, since adjacent rows have constant heights and form parallel straight lines. Taking its convergent limit preserves this invariance. Thuswhere the points lie in the domain. The limit subdivision surface is an extrusion of the univariate limit curve. This is a statement about functional data on the regular grid, not about an arbitrary spatial mesh with extraordinary vertices or independently imposed boundary stencils.

The transverse row sequence in the extrusion invariance of Loop subdivision obeysIn the convention $g_k^{\prime }={\sum }_j{a}_{k-2j}{g}_j$, the centered univariate subdivision mask isThe even coefficients sum to one and the odd coefficients sum to one; the whole binary mask sums to two. Its generating expression is $a(z)=z^{-2}(1+z{)}^4/8$, identifying the Cardinal cubic B-spline refinement rule. Its limit is $G(y)={\sum }_j{g}_j{M}_3(y-j)$ for the centered unit-grid cubic basis $M_3$. The curve across rows is this cubic B-spline; extrusion generates the whole limit surface. The along-extrusion direction itself has constant height, so its profile is a straight line rather than an additional nonconstant cubic curve.

Interpret polynomial reproduction as a guarantee for every polynomial of degree at most $d$, with the planar grid coordinates fixed by the affine-exact refinement. Constants are unchanged because every stencil sums to one. Linear height functions are unchanged because each stencil's weighted position is its refined grid location. Therefore the Loop subdivision limit reproduces all affine polynomials exactly.

Quadratic data already disprove the next degree. Use extruded data $g_j=j^2$ on a unit-spaced transverse row coordinate $y$. At spacing $h$, both even and odd rules from the induced subdivision mask give the value of $y^2$ at the refined location plus $h^2/4$. For the even stencil this follows fromand for the odd stencil from the average of the values at $y-h/2$ and $y+h/2$. Constants are reproduced, so the accumulated vertical displacement over all refinement levels isThe limit is $y^2+1/3$, not the original quadratic graph. HenceThe word “every” matters: some special higher-degree bivariate polynomials have cancelling second moments and can be reproduced. The conclusion concerns the entire polynomial space, which already fails at degree two.

Let the three grid-edge vectors be $v_1=a$, $v_2=b$, $v_3=b-a$, and define the constant-coefficient differential operator $D={\sum }_{r=1}^3(v_r\cdot \nabla {)}^2$. For a stencil at grid spacing $h$, its displacement vectors $\delta$ about the new vertex have vanishing first and third moments by central symmetry. All four stencil types have the same second momentFor the old-vertex stencil, each pair $\pm h{v}_r$ contributes $h^2{v}_r{v}_r^T/8$. For an edge parallel to $a$, the endpoint displacements are $\pm ha/2$ with weights $3/8$, and the opposite displacements are $\pm h(b-a/2)$ with weights $1/8$. Their moment is $h^2[a{a}^T/4+b{b}^T/4-(a{b}^T+b{a}^T)/8]$, precisely the same expression; the other edge orientations follow by symmetry.

If $f$ is a multivariate polynomial of total degree at most three, its Taylor expansion is exact through the third-order term. Thus every refined vertex samples the single polynomialThe polynomial $Df$ has degree at most one, so subsequent refinement reproduces it exactly. Starting at unit spacing and successively halving $h$, the limit therefore isThe correction lowers degree by two, leaving the leading homogeneous term unchanged. Hence every sampled polynomial of degree at most three generates a polynomial of that same degree, even though quadratics and cubics need not be reproduced identically. For the extruded cases this gives $y^2\mapsto y^2+1/3$ and $y^3\mapsto y^3+y$.

Degree four is not guaranteed. Extruded data $g_j=j^4$ reduce to the Cardinal cubic B-spline $G(y)={\sum }_j{j}^4{M}_3(y-j)$, which is piecewise cubic and hence cannot be a single polynomial of degree four. It is not a lower-degree polynomial either: at integer rows, $G(k)=[(k-1{)}^4+4k^4+(k+1{)}^4]/6=k^4+2k^2+1/3$, which cannot agree with a cubic at all integers. This counterexample excludes every guaranteed degree $d\ge 4$. ConsequentlyThis distinguishes polynomial generation from exact reproduction and again refers to all polynomials through the stated degree.

Use the limit subdivision surface and limit subdivision curve, rather than treating their initial control meshes as the reflecting geometry. Build adaptive refinement hierarchies for both. Each surface cell supplies a limit-position evaluator, local tangent vectors and a normal vector; each curve interval supplies a limit-position evaluator and tangent. Regular regions may be evaluated as local splines, while extraordinary neighborhoods require the evaluation machinery of the specified subdivision rules. A piecewise smooth surface is handled patch by patch, with explicit boundary or crease joins.

A useful initial polygonal construction refines the surface to small triangles and the source curve to small line segments. For a candidate mirror triangle in a plane $n\cdot (X-X_0)=0$, reflect the eye across that plane:For a source edge $Y(t)=(1-t)Y_0+t{Y}_1$, a physical reflection point is the plane intersection of the line $E^{\ast }Y(t)$:Keep only $0\le t\le 1$, $0<\lambda <1$, points inside the mirror triangle, and eye/source points on the appropriate reflecting side. Unfolding the reflected path across the mirror proves the specular reflection condition. As $Y$ moves along a straight edge, $X$ lies on the intersection of the mirror plane with the plane through $E^{\ast },Y_0,Y_1$, so the retained locus is a straight segment, a point or the empty set. Plane and triangle clipping computes its endpoints. If these defining points are collinear, or a denominator vanishes, use the corresponding limiting geometric case rather than a generic plane formula. A bounding volume hierarchy accelerates candidate search where geometric bounds exclude pairs; refined facet segments provide seeds and approximate connectivity.

To obtain the reflection curve of the actual limits, correct these seeds. In a local surface chart write $S(u,v)$, and on a source interval write $C(t)$. SetThe exact limit reflection locus is $F=0$ together with the front-side conditions. For a two-sided mirror, orient $n$ locally so that $b\cdot n>0$ and require $a\cdot n>0$; for a one-sided mirror use its prescribed orientation and discard the wrong side. At grazing incidence this orientation condition ceases to be strict and the case needs separate handling. Unlike the eye-facing assumption in the parametric question, no global orientation condition is supplied here.

Apply a continuation method: at a rank-two root of $DF$, predict in a nullspace tangent direction and correct by the Newton method on a transverse section. Map corrected parameters to $S(u,v)$, crossing patch and curve-interval boundaries through their hierarchy adjacency. Trace all branches, recognizing closed components and endpoints. Supplement facet seeds by recursive searches of the surface-cell/source-interval product, with interval arithmetic or other certified enclosures of $F$; otherwise a small loop absent from the initial polygonal approximation could be missed. Split singular neighborhoods where the Jacobian matrix loses rank, and test the open source-to-mirror and mirror-to-eye segments for occlusion when selecting actually visible components.

Allocate a geometric error budget to limit evaluation, root correction and polygonal chords. For example, require combined position-evaluation and corrected-root error at most $ϵ/2$ and, on every regular arc, require the chord error bound $M(\Delta s{)}^2/8\le ϵ/2$. Use subdivision position and normal bounds for evaluations: a tiny position error alone is insufficient, since a wrong normal vector changes the reflected ray. Away from singularities an inverse transverse Jacobian matrix bound converts reflection residual into parameter error and a bound on $DS$ converts this to space error. Refine until these bounds hold, or use validated enclosures of entire arcs and their distances to the chords.

Near a caustic the inverse bound can become large, so the hierarchy must refine more strongly or isolate the singular branch geometry. Mesh flatness alone does not certify a reflection tolerance. With regular convergent evaluators and finite resolvable components, the resulting ordered vertices and chords satisfyThe problem supplies no particular subdivision rules or smoothness guarantees, so this construction specifies what their evaluators and bounds must provide; an arbitrary nonsmooth or degenerate limit need not have an ordinary finite reflection curve. Refine for seeds, enforce reflection on the actual limit geometry, then refine the traced arcs to the required spatial tolerance.