Balancing Rotation Minimizing Frames with Additional Objectives

When moving along 3D curves, one may require local coordinate frames for visited points, such as for animating virtual cameras, controlling robotic motion, or constructing sweep surfaces. Often, consecutive coordinate frames should be similar, avoiding sharp twists. Previous work achieved this goal by using various methods to approximate rotation minimizing frames (RMFs) with respect to a curve's tangent. In this work, we use Householder transformations to construct preliminary tangent‐aligned coordinate frames and then optimize these initial frames under the constraint that they remain tangent‐aligned. This optimization minimizes the weighted sum of squared distances between selected vectors within the new frames and fixed vectors outside them (such as the axes of previous frames). By selecting different vectors for this objective function, we reproduce existing RMF approximation methods and modify them to consider additional objectives beyond rotation minimization. We also provide some example computer graphics use cases for this new frame tracking.


Introduction
In stepping along a curve, we are interested in the transition of the orthonormal frame at some point to the orthonormal frame at a neighboring point.Specifically, for a differentiable curve P(u), we are interested in using the orthonormal frame at P k = P(u k ) to construct another at P k+1 .For our purposes, the orthonormal frame at each point that we visit has the following three orthonormal vectors: a distinguished vector t, of some particular interest to us (usually the curve tangent), and two subsidiary vectors, r and s, that span the plane perpendicular to t.
Though t is fixed, unlimited orthonormal choices for r and s exist within their plane.In many applications, we would like to choose these vectors so that neighbouring frames are as similar as possible.For example, in sweep surfaces, a 2D cross section is swept along P(u) and its orientation at each u k is defined by the respective local frame.If we generate frames in a manner that can create sharp disparities between consecutive ones, as with Frenet-Serret frames, then unappealing surfaces with sharp twists result, as in Figure 1.
One way to increase neighbouring frames' similarity is minimizing the rotation between them.A rotation minimizing frame (RMF) with respect to t, also called a Bishop frame or parallel transport frame, was first proposed by Bishop [Bis75].In RMFs, every value of u has an associated orthonormal frame [t(u), r(u), s(u)] where the derivative of the vector field r(u) is parallel to t(u), and likewise for s(u).Finding r(u) and s(u) exactly requires solving a system of first-order differential equations and, thus, may not always be possi-ble or practical.In this work, we primarily reference a set of works that calculate discrete approximations of RMFs instead.We refer to these as the projection method [CW96], rotation method [Blo90], and double reflection method [WJZL08], to be consistent with the naming used in [WJZL08].
However, some of these approximations, such as the rotation method, become unstable when points are sampled too densely.Additionally, one may sometimes have additional objectives beyond the similarity of neighbouring frames; for example, if one wishes to generate a roller coaster animation for a curve, it may be desirable to not only have the cart's orientation change gradually along the curve but to also have the acceleration felt by passengers point in the direction of its floor as much as possible.These RMF approximation works do not explore how to balance rotation minimization with other constraints or objectives.Finally, approximating RMFs is not the only way to consider "similarity" between neighbouring frames; we propose to minimize some distance metric between neighbouring r and s vectors instead.By also minimizing distances between other pairs of vectors, this technique enables us to include the additional objectives beyond neighbour similarity.
A high-level visual overview of our work can be found in Figure 2. Our work starts (in Section 4) by exploring the minimization of such a distance metric.The process begins with a frame [t k , r k , s k ] at point P(u k ) = P k , having unit tangent t k at P k .We construct a following, trial frame tk+1 , rk+1 , sk+1 , at point P(u k+1 ) = P k+1 where tk+1 is aligned with the unit tangent t k+1 ; that is,  From the trial frame, we produce a preferred frame [t k+1 , r k+1 , s k+1 ] at P k+1 .The preferred frame's subsidiary vectors are chosen to be orthonormal linear combinations of the rk+1 and sk+1 that they replace.If we wish to minimize the shift from r k to r k+1 and s k to s k+1 , we can minimize, for example, More generally, we can include other objectives and minimize the weighted squared distance sum where each f R i,k+1 is a vector within the u k+1 frame (such as r k+1 ), each g i,k is a fixed vector outside of it (such as r k or something else like acceleration), and w i ∈ R is the weight.Worded another way, we solve a modified form of Wahba's problem [Wah65] where the axis of rotation is predetermined.An example f R i,k+1 and g i,k pair is shown in Figure 3.After solving this subproblem, we finish up by setting the correct sign to t k+1 = σ tk+1 so that we have the frame [t k+1 , r k+1 , s k+1 ].
The above process enables the following productions, which compose the main contributions of our work: • New derivations of existing discrete RMF approximation methods (Section 5).• A more efficient and stable calculation of the rotation method (Section 5.1).• Optimization of frames along a curve that balances RMF approximation with additional objectives (when considering Section 4 and Section 5 together).
Lastly, we produce sweep surfaces and animations that utilize such optimization, demonstrating this contribution's possibilities when applied to computer graphics (Section 6).
Figure 3: Our work involves multiple pairs of f R i,k+1 and g i,k vectors that we want to minimize the squared distance between.The f R i,k+1 vectors are part of the current frame we wish to calculate; if the frame were to rotate about the tangent, f R i,k+1 would as well.Meanwhile, the g i,k vectors are fixed in world space and would not.

RMF approximation
Frequently, smooth tracking of frames along curves utilizes RMFs.Since RMFs are defined in terms of an ordinary differential equation (ODE) [Bis75], such frame tracking may rely on an RMF approximation for cases where it cannot be assumed that the differential equation is solvable, or at least not in the time required.One may approximate the RMF by numerically solving the ODE using the Runge-Kutta method, but as Wang et al. describe [WJZL08], this is inefficient and requires a curve's second derivatives as input.This last point is a disadvantage because the RMF can still be defined for curves that are only C 1 rather than C 2 and because other approximation methods require fewer inputs (for example, just the tangents).Another approach is to approximate P(u) using simpler curves whose RMFs can be computed exactly, such as circular arcs [WJ97] or Pythagorean-Hodograph (PH) curves [JM99].
Discrete RMF approximation methods assign a frame to each of a finite set of points sampled along the curve.Multiple discrete approximation methods exist.In what we refer to as the projection method, Chung and Wang [CW96] minimize the angle between r k and r k+1 by projecting r k into the plane normal to t k+1 ; the normalization of the projection becomes r k+1 .Meanwhile, the rotation method minimizes the magnitude of the single rotation that transforms the frame at u k into the frame at u k+1 [Blo90].The axis of the resulting rotation is t k ×t k+1 , meaning the rotation method is unstable when t k ≈ t k+1 , such as if points are sampled densely or curvature is low.The double reflection method by Wang et al. [WJZL08] transforms frames at u k into the ones at u k+1 using two reflections such that, if the curve were a spherical curve, the RMF computation would be exact.The double reflection method has fourth order global approximation error, improving upon the second order error of the projection and rotation methods, and the authors provide a workaround for handling dense sampling.Other discrete approximation methods also exist, such as the discrete distance minimizing frame by Chung and Wang [CW96], a technique by Klok that projects sweep surface cross sections along the line segments of a curve's polyline approximation [Klo86], and a method by Krajnc and Vitrih utilizing quintic PH curves for interpolation [KV12].For the most part, none of these discrete approximation methods provide means to balance RMF approximation with other objectives, though a modification for closed curves to ensure the frames at the endpoints line up is provided in [WJZL08].
There can be multiple ways to derive or calculate the same approximation method.For example, Bischof et al. [BGK17] apply the constant step-size backward Euler method [HW10] to the ODE definition of RMFs and derive the projection method as the least-squares solution.Meanwhile, Chung and Wang [CW96] derive the projection method by looking for the method that produces an r k+1 whose angle with r k is minimized.Yoon et al. [YNS12] reframe the double reflection method to use quaternions instead of reflections.Sometimes, different ways of deriving the approximation methods can reveal relationships between them, such as how Chung and Wang [CW96] derive the rotation method by trying to minimize the smallest angle one can produce when comparing cos(θ)r k+1 + sin(θ)s k+1 and cos(θ)r k + sin(θ)s k for all θ ∈ R, which is similar to how they derive the projection method by only considering the angles between r vectors.Other times, there may be performance benefits; it can be shown that the double reflection method by Wang et al. is an alternative way of representing, and more efficient way of calculating, the global minimum distance intersecting frame by Chung and Wang [CW96].Likewise, Poston et al. [PFL95] demonstrate how to compute the rotation method using fewer operations than in previous works.In a similar vein, our generalized method provides alternative derivations for multiple discrete approximation methods, relating them together, and allows the rotation method to be computed with greater efficiency and stability than before.

Non-RMF frame tracking
Other works have also explored ways to track frames along curves without using or approximating RMFs, taking considerations beyond just rotation minimization about the tangent into account.Carroll et al. [CKS13] define what they call the Beta frame, whose normal is always parallel to the Frenet-Serret normal (when it exists) but may face the opposite direction to maintain smoothness.Meanwhile, Yılmaz and Turgut [YT10] introduce a frame similar to RMFs but that shares the Frenet-Serret binormal, rather than its tangent, and is rotation-minimizing with respect to the binormal instead of the tangent.Similarly, Farouki and Giannelli [FG09] use frames that are rotation-minimizing with respect to P(u)/∥P(u)∥ to define smoothly varying camera orientations.Meanwhile, Huang and Ju [HJ16] produce extrinsically smooth direction fields to produce frames for curves and surfaces that balance the goals of smoothness and shape conformity.While these works all offer alternatives or modifications to RMFs, they only take the original curve as input and, beyond that, do not offer any fine-tuned control or customization of the output frames.

Householder transformations and frame tracking
Previous works have also discussed the role Householder transformations could play in curve tracking.Lopes et al. discuss ways to apply Householder transformations to multiple problems in computer graphics and suggest that frame tracking would be another possible use, though they do not demonstrate how this could be done [LSA13].Wang et al. apply Householder transformations directly to t k to get t k+1 , and these Householder transformations are determined uniquely from the inputs [WJZL08].Meanwhile, we use Householder transformations to get trial frames and then project onto them; because we are only using the Householder transformations to create a preliminary coordinate frame at each point, we have the freedom to choose between different reference frames, and thus different Householder transformations, when creating them.Details will follow in Sections 3 and 4.

Householder transformation preliminaries
As seen in Figure 2, we need to construct orthogonal trial frames aligned with the curve tangents.We can create these frames using Householder transformations (an approach justified by works such as [LSA13]).In this section, we will quickly cover some key properties of Householder transformations that our work relies on in preparation for Section 4.1, where our technique is described in depth.

Definition and basic properties
A Householder transformation can be expressed in matrix format as: for identity matrix I and some vector h ̸ = 0.
2. H is impervious to the scaling of h; that is, h and αh yield the same H for any scalar α ̸ = 0. 3. H is orthogonal; that is H T = H −1 , since H T H = HH = I 4. Item 3. guarantees that H preserves norms and angles.
] is an orthonormal frame, then the frame [Hv 1 , Hv 2 , Hv 3 ] will also be orthonormal.5. Every Householder matrix has a determinant equal to −1. 6. Items 1. 3. and 5. imply that the Householder transformation acts on a cross product of any two vectors v and w as follows: It can also be seen that Householder transformations act as reflections across the hyperplane that passes through the origin and is orthogonal to h.The double reflection method [WJZL08] utilizes Householder transformations to perform its reflections.

Formats and computational analysis
Instead of explicitly forming and storing the matrix H, it is more efficient in storage space and computation time if we merely need to apply H to one or more vectors, w, in the following way: Here we simply form η ← 2 h T h and retain it and h as the essential part of H.
Both (1) and (3) are computationally stable and accurate under standard floating-point arithmetic.This fact has been given in numerous papers in computational linear algebra as a strong recommendation for employing Householder transformations wherever they can be used.Higham provides a comprehensive analysis of the stability and roundoff properties of Householder transformations [Hig02].Numerous computational uses of Householder transformations can be found in [GL13].

Mapping
The vector h can be constructed to map a multiple of any given vector p onto some multiple of any other vector q.In this work, we are primarily interested in the case where p and q are unit vectors, since we are working exclusively with frames that are triples of mutually orthogonal unit vectors.With unit vectors, the choice of multiples for both p and q simplifies to either +1 or −1 alone.In this case: where σ is either plus or minus 1 and represents the p and q multiples divided together onto the right.

Tracking
A visual overview of our method appears in Figure 2. Given the frame [t k , r k , s k ] at the point P k = P(u k ) of a curve, and having generated the unit tangent t k+1 at P k+1 = P(u k+1 ), we first establish an orthonormal trial frame tk+1 , rk+1 , sk+1 at P k+1 .Vector tk+1 is aligned with t k+1 : tk+1 = σ k t k+1 = ±t k+1 and the remaining two vectors, [r k+1 , sk+1 ] are basis vectors for the plane perpendicular to t k+1 , hence also to tk+1 .The details of constructing a trial frame using a Householder transformation, H k employing σ k = ±1, are discussed in subsection 4.1.
We next wish to replace [r k+1 , sk+1 ] by with the coefficients α k+1 , β k+1 , γ k+1 and δ k+1 chosen so that [r k+1 , s k+1 ] are closest in some chosen measure to [r k , s k ], they are constrained to have unit norms, and r T k+1 s k+1 = 0. From now on, we may drop the subscripts on these coefficients when they can be inferred.We discuss the optimization and finding formulas for the coefficients in subsection 4.2.
Another potentially useful condition is to require the frame at u k+1 to have the same handedness as that at u k .When r × s = +t, we call the frame right-handed, and with a negative sign on t, we call the frame left-handed.Subsection 4.3 shows how our choice of reference frame handedness and σ k , as well as α k+1 , β k+1 , γ k+1 and δ k+1 , determine the handedness of [t k+1 , r k+1 , s k+1 ], and we also show how the tracking process can be formulated a priori so that the handedness is enforced, automatically and beforehand, to that which is desired.The handedness discussion is also summarized in Figure 5.

Trial frame
To determine a trial frame where one axis is parallel to t k+1 , we first choose any Reference Frame [c k , a k , b k ] and then apply a Householder transformation to create our trial frame.Compared to other techniques for creating orthonormal frames, using Householder transformations is efficient without sacrificing numerical robustness or accuracy [LSA13].To produce the Householder transformation, we set h k = c k − σ k+1 t k+1 , where σ k+1 = ±1, and construct H k as in (3).The frame tk+1 , rk+1 is aligned with t k+1 in the sense that tk+1 = σ k t k+1 and [r k+1 , sk+1 ] are two orthonormal vectors that define/span the plane perpendicular to t k+1 .
While the Reference Frame has subscript k, this does not require that one change the Reference Frame with every tracking step.Indeed, to parallelize the computation of trial frames, one could use , where e i is the i th column of the identity matrix, that is, the i th unit vector, for all k.
Two possible uses of σ or [c k , a k , b k ] might be mentioned here.There is always a concern that the vectors c k and tk+1 which comprise the vector h k that defines the Householder matrix H k might be nearly equal, causing a problem with h k ≈ 0 and the division by h T k h k .In a typical curve tracking use case, where one would not expect t k ≈ −t k+1 , this could likely be avoided by always choosing However, if this choice cannot be used (for example, if computing the different frames in parallel, r k and s k would be unavailable), then one can instead fix the h k ≈ 0 problem by reversing the sign of σ or calling forth an alternative Reference Frame.Neither action changes the final outcome of the tracking step, since the r k+1 and s k+1 vectors of the Final Frame are always found through our formulas defining α k+1 , β k+1 , γ k+1 and δ k+1 as the optimal vectors for the k to k + 1 step.Additionally, two different reference frames with the same handedness produce trial frames of the same handedness; such frames can be generated, for example, by reversing the sign of c k and interchanging a k and b k .However, reversing the sign of σ results in the reversal of trial frame handedness that needs to be corrected subsequently.A more detailed explanation is given in subsection 4.3.
We then find the α and β that best minimize the distances between pairs of vectors of interest.Let n ≥ 1 be the number of vector pairs to consider and let R k+1 = [t k+1 , r k+1 , s k+1 ] represent a matrix that takes a vector expressed in our coordinate frame at u k+1 into world space.For 1 ] T ∈ R 3 be vectors of interest in our coordinate frame at u k+1 , making them fixed relative to said coordinate frame.Meanwhile, let g i,k ∈ R 3 be the corresponding target vectors in world space and w i ∈ R the weighting for the corresponding objective.We want to pick the α and β that determine this final frame such that the distances between R k+1 f i,k+1 and corresponding vectors g i,k ∈ R 3 are as small as possible.That is, we wish to find the R k+1 (determined by α and β) that minimizes: We may note the similarity to Wahba's problem [Wah65].The only differences are that Wahba's problem states n ≥ 2, whereas we allow n = 1, and that R k+1 can be any rotation matrix in Wahba's problem, whereas in ours the t k+1 axis is fixed.
For convenience, we can define f R i,k+1 := R k+1 f i,k+1 , which turns (6) into: To solve this minimization problem, we observe the following: The value of is constant with respect to the choice of α and β.Therefore, to minimize the error, we want to maximize: This is the dot product of two 2D vectors.Since [α, β] T is of unit length, this dot product is maximized if [α, β] T points in the same direction as the sum on the right.This means that the solution to Figure 4: A visualization of the optimization step.In (a), we see how frame choice affects distances between f R i,k+1 and g i,k vectors; we want to minimize the squared sum of these distances.In (b), we consider a specific example; on the right are values whose global positions are already known and on the left are vectors whose positions will depend on the final calculated frame.Because the distance between g i,k and f R i,k+1 along t k+1 cannot change, we can project these vectors into the plane spanned by rk+1 and sk+1 and work with corresponding 2D vectors g i,k and f R i,k+1 .For simplicity, in this example, we choose f R i,k+1 such that f R i,k+1 is of unit length and we assume that our final and trial frames have the same handedness.In (c), since jointly rotating a pair of vectors preserves the distance between them, we can rotate f R i,k+1 by R i,k+1 such that R i,k+1 f R i,k+1 is aligned with the r k+1 axis and rotate g i,k by the same R i,k+1 .The minimization problem then becomes the problem of finding r k+1 in the coordinate frame with axes rk+1 and sk+1 closest to these R i,k+1 g i,k .The solution is the unit length normalization of the weighted average of the R i,k+1 g i,k , as shown in (d).
the generalized form is the normalized weighted average of the unnormalized individual solutions for each i.
To better visualize how this method works, as in Figure 4, we can observe the following.Let us denote , and r k+1 , respectively, into the plane spanned by rk+1 and sk+1 using a 2D coordinate system with axes rk+1 and sk+1 .Because distances along t k+1 are fixed, we can minimize ∑ w i ∥g i,k − f R i,k+1 ∥ 2 instead of (7).Now, we may observe that: for some 2D rotation matrix R i,k+1 .Then we can rewrite (8) as: Using steps similar to the ones we used to find (8), one can show that maximizing (10) is the same as minimizing That is, we want to find r k+1 that, in terms of rk+1 and sk+1 , is as close as possible to our various rotated and scaled R i,k+1 f R i,k+1 g i,k vectors.This observation allows us to produce the visualization shown in Figure 4.In the figure, for the sake of simplicity, we only selected f R i,k+1 such that we have f R i,k+1

2
= 1 and we also label r k+1 as just r k+1 .This figure also ignores the ρ = −1 case and how one in general would ensure consistent frame handedness, which is the focus of subsection 4.3 and Figure 5.

Handedness
In this subsection, we discuss how to determine and control the handedness of final frames.First, we can show that the handedness of the trial frame tk+1 , rk+1 If, similar to our use of ρ and σ, we use the variable ω = ±1 to encode the reference frame's handedness, then: Next comes the impact of ρ.To save space, the below omits the  subscripts from r k+1 , s k+1 , rk+1 and sk+1 .Using properties of the cross product, we may observe that: Thus, the handedness of the pre-final frame tk+1 , r k+1 , s k+1 is the same as that of trial frame tk+1 , rk+1 , Finally comes the impact of σ.Since tk+1 = σt k+1 , it is clear that the handedness of the final frame [t k+1 , r k+1 , s k+1 ] is the same as the handedness of the pre-final frame tk+1 , r k+1 , s k+1 if σ = +1 and the opposite if σ = −1.
The above is summarized visually in Figure 5.Using our variables σ, ρ, and ω, we can also produce the following summary:

Pseudocode implementation
We have compiled the details of our frame tracking approach as pseudocode in Algorithm 1.It is important to note two things about this pseudocode.First, as mentioned before and as we will demonstrate explicitly in Section 5, for some values of i, the g i,k may be functions of r k and s k rather than constants; this enables RMF approximation.Secondly, the pseudocode covers the most general case; more efficient code is possible in specific cases, such as our rotation method calculation in Section 5.1.Supplementary material contains an operation cost breakdown.

Deriving and calculating discrete RMF approximation methods
With our generalized technique, we can derive multiple preexisting discrete RMF approximation methods.This reveals relationships between them.Additionally, the proof that we can derive the rotation method also showcases a new method of computing it that is more efficient and stable than previous ones.

Deriving and calculating the rotation method
Let us minimize (7 That is, we wish to minimize: From (8), we can see that We can show that (16) generates the same frames as the rotation method [Blo90], and we can also show how these frames can be calculated in fewer steps and with more stability than in previous works.For this proof, we assume right-handed [t, r, s] frames; the proof for left-handed frames would be nearly identical.Since both methods are coordinate invariant, we can choose to work in a coordinate system where t k and t k+1 lie in the xy-plane and The rotation method rotates frame k to frame k + 1 by the smallest angle possible.If t k = t k+1 (which, in our setup, happens if φ = 0), then the matrix representing such a rotation is just I. Otherwise, the method calculates an axis of rotation t k ×t k+1 ∥t k ×t k+1 ∥ and an angle of rotation arccos(t k • t k+1 ) and then performs this rotation on r k and s k .In our setup, the axis of rotation is clearly [0, 0, 1] T and the angle of rotation is arccos(t k • t k+1 ) = φ + 2πm for some m ∈ Z.Therefore, regardless of the value of φ used to construct t k+1 , we can represent the rotation method with the matrix Algorithm 1 Frame Tracking Input: Unit tangent vectors t k , vectors f i,k+1 and g i,k , and weights w i for k = 1, 2, . . ., m + 1 and i = 1, 2, . . ., n.Additionally, a choice of final frame handedness.Finally, g i,1 may rely on an initial frame unit vectors r 1 and s 1 , making them inputs as well.Output: ▷ Choice of reference frame is flexible (see Section 4.1). 3: 4: σ = chooseSigma(...) ▷ Either -1 or +1; also flexible (see Section 4.1). 5: ▷ Determines Householder transformation (i.e., reflection) plane. 6: ▷ Set trial frame vectors via Householder transformation.
7: for i from 1 to n do ▷ Performs summation described in Section 4.2. 11: , and f s i,k+1 . 12: end for 13: [α, β] T = sum/∥sum∥ ▷ Normalize the weighted sum variable. 14: ▷ We now calculate the final frame vectors using the trial frame ones. 15: ▷ Could also calculate s k+1 via cross-product from t k+1 and r k+1 .16: end for Then it is clear that the rotation method generates vectors For our approach, we must first choose a reference frame and σ to produce rk+1 and sk+1 .Any choice produces the same result (so long as we choose the correct ρ, based on the other choices' impact on handedness).For this discussion, we choose σ = −1 and a reference frame of [t k , r k , s k ] and we will show that these choices result in rk+1 = r k+1 and sk+1 = s k+1 , a fact one can use to save computation steps when calculating the rotation method compared to previous works.These choices do not work if t k = −t k+1 , as then we would have h = 0, but the steps below could in that case be replicated with a different reference frame so that all cases for t k+1 are covered.Our choice of reference frame means that h = t k − σt k+1 = t k + t k+1 = [1 + cos(φ), sin(φ), 0] T .Therefore, using (3) and simplifying, we have: We can see that rk+1 then matches the r k+1 given in (18).A similar set of steps shows that the sk+1 for this reference frame and σ also matches the s k+1 in (18).
While our approach and the typical way of computing the rotation method generate the same results in theory, when working with floating-point representations of the inputs, our approach can be more robust when densely sampling a curve.Behaviour under dense sampling is important to consider because, in theory, increasing the number of samples along a continuous curve improves discrete approximation of its RMFs.Specifically, let u = 0 to u = L represent the parameter range of a curve and suppose the curve is evenly sampled so that u k = (k − 1) ℓ for k = 1, 2, . . ., m + 1 and ℓ = L/m.Then prior work shows that the double reflection method's error is O(ℓ 4 ) [WJZL08] and the rotation method's error is O(ℓ 2 ) [PFL95].The idea that error would decrease as sampling increases is not surprising.For example, consider using only two samples: the start and end of the curve.This would be highly inaccurate, because many different curves have identical endpoint locations and tangents, but discrete RMF approximation would assign them all one common pair of frames unless more samples are taken.
However, because of finite floating-point precision, the rotation method lacks robustness when t k and t k+1 are too similar, such as when sampling a curve densely enough (especially in low curvature areas).If t k ≈ t k+1 , then t k × t k+1 ≈ 0, making the axis of rotation ambiguous.Meanwhile, if t k ≈ t k+1 , then our method just involves a reflection across span{h} ⊥ , where h = t k + t k+1 ≈ 2t k , which is well-defined.
To observe this, we can construct an example using a Pythagorean-hodograph (PH) curve, for which exact RMFs may be computed via the process in [Far02].We create a cubic Bézier PH curve from desired endpoints P 0 = [0, 0, 0] T and P 3 = [23, 266, 1] T and respective unit tangents [0, 1, 0] and [23/265, 264/265, 0] T , chosen so that the curve has low, but still nonzero, curvature for most of its length.We can use [JM99] to find the correct middle control points P 1 = [0, 2, 0] T and P 2 = [0, 2, 1] T .To rewrite the curve in the form required by [Far02], we follow [FGS15].Figure 6 shows how, compared to those generated by the rotation method, frames generated by the negate-and-reflect method for this PH curve are closer to the exact frames found using [Far02] 6: Results of comparing the negate-and-reflect method (ours) and rotation method for our sample curve.Error represents the max angle (in radians) between exact r k generated by [Far02] and approximate r k across all 1 ≤ k ≤ m + 1, where m is the number of curve segments in our discrete approximation.We used 32bit floats for our computations.The plot uses "Rotation*" and "Ours*" to show what happens if we take the extra step of normalizing coordinate frame vectors before calculating the next frame.With or without this extra correction, the negate-and-reflect method produces lower error as the number of segments n increases.
Another benefit of our calculation over previous ones is the number of steps required.An efficient implementation of the rotation method, which we used to generate Figure 6, can be found in [PFL95], and an analysis of this implementation in [WJZL07] shows that, when written as efficiently as possible, it uses 26 additions, 36 multiplications, and one division to compute r k+1 and s k+1 .In comparison, if one uses our negation-and-reflection method with σ = −1 and [c k , a k , b k ] = [t k , r k , s k ], then our method only requires only 13 additions, 16 multiplications, and one division (one first calculates r k+1 = rk+1 as described in (3) and then computes s k+1 via cross product; this paper's supplementary material has a detailed breakdown).This is summarized in Table 1.

Table 1:
Step counts for the efficient rotation method implementation in [PFL95] versus our negation-and-reflection method.
Note that, in general, there is nothing special about minimizing the distance between r vectors compared to minimizing the distance between s vectors.One could choose to instead project s k into the plane orthogonal to t k+1 and normalize to obtain s k+1 and then obtain r k+1 via cross product.One could even alternate between the two: for example, one could calculate ř as the normalized projection of r k onto span{t k+1 } ⊥ , calculate š as the normalized projection of s k onto span{t k+1 } ⊥ , and then choose to set r k+1 = ř or s k+1 = š based on how ∥r k − ř∥ and ∥s k − š∥ compare.
Chung and Wang [CW96] extend this last idea to derive the rotation method.Rather than just considering the difference between r and s vectors, however, they aimed to find the r k+1 and s k+1 that would minimize the smallest possible angle (directly related with distance for angles less than π radians) between any cos(θ)r k + sin(θ)s k and cos(θ)r k+1 + sin(θ)s k+1 in the consecutive frames.The smallest possible angle between such two vectors is the angle between the two planes they would lie in, which is also the angle between t k and t k+1 , and so the rotation method provides the solution.Our version of the rotation method in (16) can also be seen as extending the projection method to derive the rotation method; there, the final α and β can be interpreted as the normalized weighted average of the α and β we would receive by using the projection method on r and the projection method on s, with the weights being the lengths of the respective projections.

Approximation methods that utilize frame origins
All of the RMF approximation methods we have re-derived so far only consider the distances between vectors, but consecutive frames along a curve do not share the same origin; thus, considering the positional offset of the frames in the approximation is also worth considering.Chung and Wang, for example, consider minimizing the distance between P(u k ) + r k and P(u k+1 ) in what they call their discrete distance minimizing frame [CW96].Using Lagrange multipliers, their solution is the normalized projection of the vector r k − P(u k+1 ) + P(u k ) into the span{t k+1 } ⊥ plane.While we omit the steps to show this, as they are straightforward, one can see that our generalized setup provides the same solution when setting n = 1, w i = 1, f R 1,k+1 = r k+1 , and g 1,k = r k − P(u k+1 ) + P(u k ), which, as required, minimizes Of course, depending on the application, one may want to minimize multiple such distances.One case that gives an interesting result is minimizing the total squared distance between P(u k ) + v k and P(u k+1 ) + v k+1 for the n = 4 cases v = ±r and v = ±s.In other words, for a sweep surface with cross section vertices P(u k ) ± r k and P(u k ) ± s k (four vertices arranged in a square), this would be the sum of the squared lengths of the edges between the cross sections.By plugging the values into Equation ( 7) and simplifying, one can see that the given error is minimized if and only if the error in Equation ( 15) is minimized, meaning the same solution as in ( 16), and thus the rotation method, applies.In fact, the result holds for any polygonal cross-section with rotational symmetry of order four.
Lastly, to a more limited degree, we can relate our method to the double reflection method [WJZL08].In that method, one computes two Householder transformations, namely H 0 , with h = P(u k+1 ) − P(u k ), and H 1 , which maps H 0 t k to t k+1 .Then one sets r k+1 = H 1 H 0 r k and s k+1 = H 1 H 0 s k .Our negation-and-reflection method from Section 5.1 is identical except for replacing H 0 with a different reflection, one taking t k to tk+1 = σt k = −t k without changing r k or s k , as they lie in the hyperplane that is reflected about.Alternatively, one can say that the double reflection method is the same as the rotation method if P(u k+1 ) − P(u k ) is estimated by t k .
As for actually calculating the double reflection method, consider an orthonormal frame which we know has the same handedness as [t k , r k , s k ] because of (2).Now suppose one minimizes the following: We know from Section 5.1 that we can minimize this by finding the Householder transformation H ′ that maps −t L k to t k+1 and setting r k+1 = H ′ r L k and s k+1 = H ′ s L k .But since −t L k = H 0 t k , we can conclude that minimizing (19) yields the same result as the double reflection method.This adjustment of Section 5.1 performs the exact same computations as the original double reflection method, meaning it offers no additional computational stability or efficiency.However, this reframing of double reflection, in addition to the above insight into how it compares to the rotation method, allows augmentation with additional objectives (similar to how we augment the rotation method in Section 6).

Example applications to computer graphics
In this section, we demonstrate some sample applications of the generalized frame tracking from Section 4.2 to computer graphics.Because of the ability to balance multiple objectives, and because objectives are provided in a way that would be easy for artists to understand and visualize (vector fields along the curve), we imagine there could be many applications beyond what we showcase here.In all figures referenced in this section, the RMF approximation method that we augment with additional objectives is the rotation method, meaning that frame origins do not impact the calculation as they do with methods like double reflection; we only rely on curve tangents.Therefore, like in (15) for the rotation method, f R 1,k+1 = r k+1 , f R 2,k+1 = s k+1 , g 1,k = r k , g 2,k = s k , and w 1 = w 2 .Then we include g 3,k as our additional objective, use w 3 as its weight, and set f R 3,k+1 = s k+1 as the in-frame vector to align the additional objective with.More visuals, including animated video, can be found in this paper's supplementary materials.
When simulating a roller coaster, one generally wants to have the cart's floor aligned so that its normal force opposes the gravity and provides the acceleration of the passengers.This way, passengers are pulled down into their seats.However, in artistic situations where exact realism is not required, one might want to smooth out the sharp twists that can result from this choice of normal vector, such as the twist seen in Figure 1a.At the same time, a pure RMF might capture none of the acceleration information, as seen in Fig- ure 1c.Our version in Figure 1d offers a compromise between the two, with the contribution of each factor weighted by the user.Here, f R 3,k+1 = s k+1 and g 3,k is aligned with the direction of passengers' acceleration minus gravity.Acceleration can be calculated numerically by assigning the cart a fixed speed at the highest point of the coaster, using conservation of energy to compute the speed at each point, using the tangent vectors to compute velocities, and then using the velocities to approximate acceleration.A figure showing the effect of different values for the w 1 = w 2 influence of the RMFapproximating component can be found in this paper's supplementary material.
In Figure 7 we model a snake as a sweep surface.In this case, we set f R 3,k+1 = s k+1 and g 3,k to point directly upwards in world-space.This allows us to enforce that all parts should generally be groundlevel (i.e., that r k+1 should be parallel to the ground plane), which prevents the head from turning upside-down as it would with a pure RMF approximation.We use the same configuration to animate a drone's flight, as seen in Figure 8.One thing should be noted about closed curves.Currently, our work does not ensure that the first and last frames generated are equivalent.The solution by Wang et al. [WJZL08] for pure RMFs -to add the angle difference between the two, in an accumulating fashion, to all frames along the curve so that the last two line up -is not ideal for our proposed applications.For example, along lengthy straight sections of a roller coaster, the track would keep twisting as this difference is distributed gradually along it, whereas the carts could be pointing straight upwards to both align with acceleration and reduce the differences between frames.A workaround for now could be to apply this technique over a segment of the curve where such effect is not noticeable.

Future work
Currently, our work offers a solution to a modified version of Wahba's problem where rotation is constrained about a single axis and thus has only one degree of freedom.Meanwhile, in the original Wahba's problem, the rotation being optimized has three degrees of freedom.It would be interesting to consider the case where    there is exactly two degrees of freedom.Part of the motivation here would be for the sake of completeness, though it could still be relevant to considering different ways that RMF approximations can be derived and calculated, as Chung and Wang derive the rotation and double reflection method by minimizing the distance between a single pair of vectors where each vector has the freedom to rotate about its respective tangent, creating a problem with two degrees of freedom.In comparison, the future work we are proposing would be figuring out whether one can accommodate arbitrary pairs of vectors in these two frames.
In this work, we have focused on optimizing distances only between consecutive frames.However, one may instead be interested in minimizing the general error in (7) globally, as we imagine there may be situations in graphics, robotics, and more where it is this global error, rather than local error, that one wants to minimize.This may also offer us a better approach to handle closed curves in our generalized optimization.
Finally, another possibility for future work is to use the freedom offered by our generalized technique to track coordinate frames across entire surfaces rather than just along curves.As with the issue our current method faces for closed curves, one primary challenge here would be closed surfaces.

Conclusion
In this work, we have shown how one can track frames along curves by considering the distances between vectors in consecutive coordinate frames as well as vectors representing other objectives.We also demonstrated how to represent some existing RMF approximation techniques using this method.The benefits of this are relating the techniques together in new ways, deriving a new efficient and stable way to compute the rotation method, and allowing existing RMF approximation methods to be balanced with other objectives by adding other distances to minimize into the objective function.Finally, we have discussed some examples of how the last benefit can be applied to create 3D models and animations.

Figure 1 :
Figure 1: A comparison of our generalized system compared to other curve tracking options.For a roller coaster track example, frames based on only acceleration (a) and Frenet-Serret frames (b) have sharp twists.Meanwhile, with RMFs (c), passengers' acceleration relative to the cart has no effect on the shape.Our generalized frame tracking (d) allows us to create a "compromise" track between (a) and (c).

Figure 2 :
Figure 2: An overview of our technique.When calculating the frame at P(u k+1 ), we may assume that we know its tangent t k+1 and that we know the previous frame [t k , r k , s k ] at P(u k ).Using reference frame [c k , a k , b k ] and Householder transformation H k (justified in Section 3), we construct a trial frame tk+1 , rk+1 , sk+1 aligned with t k+1 , as discussed in Section 4.1.Then, in an optimization step described in Section 4.2, we find appropriate α, β, γ, and δ and set r k+1 = α rk+1 + β sk+1 and s k+1 = γ rk+1 + δ sk+1 .The effects of [c k , a k , b k ], H k , α, β, γ, and δ on the handedness of [t k+1 , r k+1 , s k+1 ] are discussed in Section 4.3.

Figure 5 :
Figure 5: Handedness influence of various tracking options, as discussed in Section 4.3.
final frames should be right-handed, else ω σ ▷ Logic for choice of ρ comes from Section 4.3.

Figure 7 :
Figure7: A snake sweep surface example demonstrating three frame tracking options.For each, the top row contains a side view and a birds-eye view of the render, while the bottom row contains a larger view of the twisting region highlighted in the side view.Using only RMF approximation (a), if one chooses the starting frame so that the belly lies on the ground, then the head is upside-down.Meanwhile, if we simply choose the coordinate frames so that one of the axes is parallel to the ground (b), then sharp twists can appear.Our version (c) allows one to balance smoothness with ground alignment.

Figure 8 :
Figure 8: A drone example.Using only RMF approximation (a), the drone flies upside-down for a substantial period.Meanwhile, if we simply choose the coordinate frames so that one of the axes is parallel to the ground (b), the drone's orientation rapidly flips during its climb.Our version (c) smooths out the transition.
© 2023 The Authors.Computer Graphics Forum published by Eurographics and John Wiley & Sons Ltd.