Getting consistent normals along a 3D (Bezier) curve

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: • It would be better if you showed the “2-D normals” for the same curve. – amd Jul 6 '18 at 23:40
• Do you have some concrete way to characterize the “right” direction for the normal in 3-D? Which ones in your first illustration are the “correct” ones? Does the answer change if you look at it from the other side? – amd Jul 6 '18 at 23:42
• In that case you should orient and scale them so that they look the same. The curve in the first illustration bends left as you move up the page, but the second one bends to the right. – amd Jul 6 '18 at 23:53
• @Mike'Pomax'Kamermans: In my answer, I used the exact same matrices as you have, just combined them into vectors (using Maple), so they should match yours exactly. Even if they do not, don't mind the differences, as the numerical values are not the problem here. – Nominal Animal Jul 7 '18 at 2:10
• – Oppenede Jul 8 '18 at 9:11

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;

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

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

We have \begin{aligned} B(t) =& \left [ \begin{matrix} -0.54 t^3 + 1.53 t^2 - 1.14 t \\ 0.02 t^3 + 0.15 t^2 + 8.04 t \\ 0 \end{matrix} \right ] \\ B^\prime(t) =& \left [ \begin{matrix} -1.62 t^2 + 3.06 t - 1.14 \\ 0.06 t^2 + 0.30 t + 8.04 \\ 0 \end{matrix} \right ] \\ B^{\prime\prime}(t) =& \left [ \begin{matrix} -3.24 t + 3.06 \\ 0.12 t + 0.30 \\ 0 \end{matrix} \right ] \\ \; & 0 \le t \le 1\end{aligned} with $B(t)$ the point on the curve, $B^{\prime}(t)$ the direction (velocity or tangent), and $B^{\prime\prime}(t)$ is curvature (or acceleration), at $t$.

If the curve curves (changes direction, even infinitesimally) at $t$, then the direction vector rotates around vector $n(t)$ at that point: $$n(t) = B^{\prime}(t) \times B^{\prime\prime}(t)$$

If the curve is straight at $t$, then $B^{\prime\prime}(t) = 0$.

A cubic curve can be straight, but "accelerate" or "decelerate". (If you consider a line parametrized as a cubic curve, putting the two control points along the line segment between the starting and ending points, will not change the shape of the curve at all, only how the curve is formed; i.e., the function describing the ratio $\lVert B(t) - B(0)\rVert / \lVert B(1) - B(0)\rVert$.) This "acceleration" and "deceleration" occurs when $B^{\prime} \parallel B^{\prime\prime}$. Because at such points the curve does not change its direction (the direction of $B^{\prime}(t)$ does not change, only its magnitude changes), $n(t) = 0$.

Herein lies the problem.

In two dimensions, the axis of rotation (around which you rotate $B^{\prime}(t)$ to get the perpendicular vector) stays constant; it is the "implicit" third axis. (It is implicit in the way the 2D analog of a cross product is defined.)

In three dimensions, the axis vanishes at $t \approx 0.9398$. Because all $z$ components for the curve are zero, when the rotation axis exists, it is parallel to the $z$ axis. Before $t \approx 0.9398$, the axis is towards $+z$, after it is towards $-z$.

A simple answer is to use $n(t) = N$ (so that the normal vector OP is looking for is $B^{\prime}(t) \times N$), with $N$ defined as the normal to the plane the curve lies in. (A cubic curve has four control points. If these lie in the same plane, they curve lies in that plane too. $N$ is the normal to this plane.)

That yields the exact same normal vector for the 3D-extended 2D curve.

The question of what to do with curves that are not planar, and actually curve in all three dimensions, remains.

If you switch to the $n(t)$ mentioned above, the normal vectors thus calculated will behave similarly to OP's image, even when the curve is otherwise planar. You could require it is in the same halfspace as $n(0)$, but it would still leave the point-like discontinuity at $t$ where $B^{\prime}(t) \parallel B^{\prime\prime}(t)$.

The proper answer to this remaining part is that it depends on what these normals are used for.

In a practical sense, when you extend a line or curve from 2D to 3D, you should get a surface to get analogous properties. If you keep the objects dimensionality unchanged, you get an additional degree of freedom. (If you extended the 2D cubic curve into 3D cubic patch, you'd have an analogous surface normal, too.)

This additional degree of freedom means that while a line in 2D has a normal vector, a line in 3D has a normal plane (as defined by its normal, $B^{\prime}(t)$).

Similarly, while a curve in 2D has a normal vector (the tangent, $B^{\prime}(t)$, rotated 90° clockwise or counterclockwise), a curve in 3D has a normal plane (as defined by $B^{\prime}(t)$).

• It feels like there should be an approach here, somewhere, based on the fact that any non-planar curve can be described as a planar curve as long as we restrict the interval we're examining (in the same way that any curve can locally be considered equivalent to its tangent). While it is possible to catch the B''(t) transition from +z to -z, and flip the axis or rotation, that technically leaves a discontinuity at the point where B''(t)=0. I'm now trying to see if I can make things at least "work" by using the placement of the circle implied by the curve's radius of curvature. – Mike 'Pomax' Kamermans Jul 7 '18 at 17:05
• @Mike'Pomax'Kamermans: The issue is that a straight line, or a straight segment of a curve, has no normal vector; only a normal plane. If you extend from a 2D curve to a 3D curve, you need to extend the concept of a normal vector to a normal plane, to get analogous properties; that's what the extra degree of freedom does. To me, it feels like a desperate attempt of shoehorning a concept into a situation it does not fit in. Maybe you could expand in your question a bit about the purpose of these normal vectors? – Nominal Animal Jul 8 '18 at 2:02
• That is indeed the problem. Another is the nomenclature: the vector I'm looking for is not "the" normal vector as you point out, it's really more an aesthetically useful vector in the normal plane, for things like graphical design and visual artistry work - for instance, offsetting along these "normals" of a curve and getting something that yields the kind of contour designers and artists can actually make use of, rather than offsets that kink and flip wildly (as in this example). – Mike 'Pomax' Kamermans Jul 8 '18 at 5:00
• @Mike'Pomax'Kamermans: All I use are the Wikipedia articles on versors, axis-angle representation, Quaternions and spatial rotation, and Quaternion (Hamilton product). The key is that you can normalize a quaternion $(w, x, y, z)$ by dividing each component by $\sqrt{w^2 + x^2 + y^2 + z^2}$ after operations, [...] – Nominal Animal Jul 8 '18 at 17:47
• [...] and blend between two quaternions as you would scalars; just normalize the result. An unit quaternion or versor describes an orientation. You can think of it as describing an axis, and an angle around that axis. If you get your mind around those concepts, that's about it: quaternions are really much easier to use than for example Euler angles! – Nominal Animal Jul 8 '18 at 17:49

The normal vector of a parametric curve $B(t)$ can be found in the following way:

1) compute the first derivative $B'(t)$ and the 2nd derivative $B''(t)$,
2) compute the unit tangent vector $t=B'(t)/|B'(t)|$,
3) compute the bi-normal vector $$b= \frac{B'(t) \times B''(t)}{\left|B'(t) \times B''(t)\right|}$$
4) compute the normal vector as $n= b \times t$.

The 3 unit vectors $t$, $b$ and $n$ are mutually orthogonal and will form the so-called "Frenet-Serret frame".

• Those are the normals I already have, and are super unreliable for graphics work (the first image shows them based on unnormalized vectors, but their alignment is identical). I need "normals" that lie in the normal plane, but oriented to fit a regular viewer's idea of "sticking out from the curve in a relatively uniform manner", so in the stated case, 2D and 3D showing the same result. If those can't be called normals, then that's fine: I don't care about nomenclature in this case, only that the actual vectors look right. – Mike 'Pomax' Kamermans Jul 7 '18 at 21:17
• For points lying on a plane, it may be possible to check the orientation of the binormal by examining the signs of elements of the curvature vector given by $\boldsymbol{\kappa} = (B'(t) \times B''(t))/|B'(t)|^3$. It may also to useful to examine the torsion for more general situations. – Biswajit Banerjee Jul 7 '18 at 22:22
• A possibly better way is to determine whether the normal needs to be flipped is make sure that the Frenet-Serret frame is always right-handed, i.e., check whether $(t \times b) \cdot n > 0$. – Biswajit Banerjee Jul 7 '18 at 22:39
• Upon further thought, the frame will always have the same orientation. So, the only robust way to check the flip is to take two adjacent points on the curve, check $t_1 \cdot t_2 > 0$ and $b_1 \cdot b_2 > 0$ and flip when the sign is negative. – Biswajit Banerjee Jul 7 '18 at 22:54
• That's what I initially thought, too, but that's still not enough: both dot products are positive on both sides of the crossing point, only having a flipped sign where we actually cross from positive to negative second derivative. – Mike 'Pomax' Kamermans Jul 8 '18 at 5:15