[Math] Getting consistent normals along a 3D (Bezier) curve

Question

I'm trying to get consistent normals along a 3D Bezier curve $B(t)$, where for any point I compute the normal as:

$$
\begin{align}
\vec{a} &= B'(t) \\
\vec{b} &= B''(t) \\
\vec{c} &= \vec{a} + \vec{b} \\
\vec{r} &= \vec{c} × \vec{a} \\
\vec{n} &= \vec{r} × \vec{a} \\
\end{align}
$$

So, get the derivative at a point for time value $t$, and implicitly get the plane of curvature at the point by computing the cross product of the derivative vector at the point, and the "next" derivative vector we get from moving the derivative by the amount dictated by the second derivative. The cross product yields the axis of rotation, so to then form the normal at the point for time value $t$ I take the cross product of the axis of rotation, and the original derivative vector, since these three vectors are by definition perpendicular.

The problem is that normals computed this way are not consistent: they will "flip" around inflections, and I'm not sure what the right way is to go about making sure that does not happen.

As visual illustration, consider the following 3D cubic Bezier curve:

$$
B(t) =
\left[\begin{matrix}1&t&t^2&t^3\end{matrix}\right]
\left[\begin{matrix}1&0&0&0\\-3&3&0&0\\3&-6&3&0\\-1&3&-3&1\end{matrix}\right]
\left[\begin{matrix}
0 & 0 & 0\\
-0.38 & 2.68 & 0\\
-0.25 & 5.41 & 0\\
-0.15 & 8.21 & 0
\end{matrix}\right]
$$

Now, this happens to be a 3D curve that lies entirely on the x/y plane, but it illustrates the problem rather well. The above procedure yields the following normals:

However, this is rather different from the 2D normals we get when taking advantage of the 2D plane, where a normal can be constructed by simply rotating the (normalised) derivative vector a quarter turn clockwise, setting $(x,y)$ as $(-y,x)$:

I'd like to get something similar to the 2D case for the 3D case, but I don't know how to ensure that the cross products are unaffected by "which direction" the second derivative moves the derivative across its plane of curvature

(Effectively, how do I ensure that, when considering the triplet {normal,derivative,axis of rotation} that these always map to the local {x,y,z} axes, rather than sometimes mapping to {x,y,z} and somethings mapping to {y,x,z} axes)

Edit

While more "algorithmic" than I'd like, the only workable solution I've found so far is to compute the normals for two points $B(t)$ and $B(t+\varepsilon )$, then computing the angular difference in the plane for those two normals,

$$
\theta = \textit{acos} \left ( \frac{n_1 \cdot n_2 }{||n_1|| ||n_2||} \right )
$$

and then check whether that value is close to $\pi$ or not. Even in fast-changing curves, the angle between two "reasonable" normals is a relatively small value, so if the angle suddenly flips to "nearly $\pi$" then as of that time value the "desired normals" are negative actual normal.

While that works, it feels kind of hacky.

Without algorithmic flipping:

With algorithmic flipping:

Note this does not affect cuves with "reasonable twisting", e.g. when we set the $z$ values to $\{0,200,-200,600\}$ for the first, second, third and fourth control point respectively:

Answer 1

Thanks to @Oppenede commenting on the question post, it turns out that what I was looking for in this case is called the "Rotation minimizing frame" of a point, also known as the "parallel transport frame", or "bishop frame".

This is an algorithmic procedure, where you compute the orthogonal vector triplet $\{\textit{tangent}, \textit{rotation axis}, \textit{normal}\}$ for the point at time value 0, and then compute subsequent frames based on "the previous frame", using an ever so slightly modified version of the procedure explained in section 4 or "Computation of Rotation Minimizing Frames" (Wenping Wang, Bert Jüttler, Dayue Zheng, and Yang Liu, 2008):

ArrayList<VectorFrame> getRMF(int steps) {
  ArrayList<VectorFrame> frames = new ArrayList<VectorFrame>();
  double c1, c2, step = 1.0/steps, t0, t1;
  PointVector v1, v2, riL, tiL, riN, siN;
  VectorFrame x0, x1;

  // Start off with the standard tangent/axis/normal frame
  // associated with the curve just prior the Bezier interval.
  t0 = -step;
  frames.add(getFrenetFrame(t0));

  // start constructing RM frames
  for (; t0 < 1.0; t0 += step) {
    // start with the previous, known frame
    x0 = frames.get(frames.size() - 1);

    // get the next frame: we're going to throw away its axis and normal
    t1 = t0 + step;
    x1 = getFrenetFrame(t1);

    // First we reflect x0's tangent and axis onto x1, through
    // the plane of reflection at the point midway x0--x1
    v1 = x1.o.minus(x0.o);
    c1 = v1.dot(v1);
    riL = x0.r.minus(v1.scale( 2/c1 * v1.dot(x0.r) ));
    tiL = x0.t.minus(v1.scale( 2/c1 * v1.dot(x0.t) ));

    // Then we reflection a second time, over a plane at x1
    // so that the frame tangent is aligned with the curve tangent:
    v2 = x1.t.minus(tiL);
    c2 = v2.dot(v2);
    riN = riL.minus(v2.scale( 2/c2 * v2.dot(riL) ));
    siN = x1.t.cross(riN);
    x1.n = siN;
    x1.r = riN;

    // we record that frame, and move on
    frames.add(x1);
  }

  // and before we return, we throw away the very first frame,
  // because it lies outside the Bezier interval.
  frames.remove(0);

  return frames;
}

(This uses Java syntax but should be easy enough to parse for porting to any other language)

The result of this procedure leads to rather aesthetically pleasing normals. For the original planar curve, we see the following, with the RMF normals in green and the original normals in blue:

And for the non-planar curve, again with RMF normals in green and original normals in blue:

So this works really well. The downside of course is that this means normals can no longer be computed "on demand", as each frame relies on the previous frame, necessitating a full RMF computation pass and then interpolating for missing normals. But, based on the literature available, there does not appear to be a way to get nice, consistent looking normals without an iterative approach like this.

So Rotation Minimizing Frames it is!