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:

 A: 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!
A: 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)$).
A: 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". 
A: (I came across this question trying to debug my implementation of essentially the same problem. The stated answer did not work for my case. It could be that I implemented it wrong, or maybe the stated answer works for less curvy curves than I need [see below]. I did manage to fix my own initial implementation, and it is conceptually simpler than the procedure in the accepted answer, so I thought it would be a good additional answer, regardless of whether the accepted answer is sometimes inferior.)
Given a curve approximated by a sequence of points and tangents, find a sequence of normals that follows the curve and follows some intuition of non-drifting. (Actually, we need a complete "frame", which consists of a sequence of three mutually orthogonal vectors, made from the tangents, the "consistent" normals, and their cross product.)
Let $t_i$ be the sequence of tangents.
We compute the three axes $x_0$, $y_0$, and $z_0|$ of the first frame using Method 1 described in this answer. It is simply a way to find a good set of orthogonal vectors in a plane, and we use the plane at the first point of our curve with the same normal as our first tangent.
We then compute each frame iteratively as follows:
$z_i = t_i$
$y_i = t_i \times x_{i - 1}$
$x_i = y_i \times z_i$
The idea is basically we start with some normal. Then the next normal is a vector that 1) is perpendicular to the tangent, and 2) lies in the same plane as the tangent and previous normal. The cross products are just calculating such a vector.
I use this frame to calculate a "thin tube" around a curve (represented by a sequence of points and tangents).
Below on the left, you can see the result of the implementation in the accepted answer (and remember my implementation may be wrong rather than the answer itself), and on the right is the implementation described here. Notice how the seem (the black line) drifts from the top to the bottom in the left mesh, but stays on the bottom in the right mesh.

