Getting consistent normals from a 3D cubic bezier

Question

I'm writing a BezierPath class that contains a list of BezierPoints. Each BezierPoint has a position, an inTangent and an outTangent:

BezierPath contains functions for getting linear positions and tangents from the path. My next step is to provide functionality for getting the normals from the path.

I'm aware that any given line in 3D is going to have an unlimited number of lines that are perpendicular, so there isn't going to be a set answer.

My aim is for the user to be able to specify normals (or a roll angle?) at each BezierPoint which I will interpolate between to get normals along the path. My problem is that I don't know how to choose a starting tangent (what should the default tangent be?).

My first attempt at getting the starting tangents is using the Unity3D method Quaternion.LookRotation:

Quaternion lookAt = Quaternion.LookRotation(tangent);
Vector3 normal = lookAt * Vector3.up;
Handles.DrawLine(position, position + normal * 10.0f);

This results in the following (green lines are tangents, blue are normals):

For the most part this seems like a good base result, but it looks like there are sudden changes in certain orientations:

So my question is: Is there a good way for getting consistent default normals for lines in 3D?

Thanks, Ves

Answer 1

Getting the normal for a point on a Bezier curve is actually pretty straight forward, as normals are simply perpendicular to a function's tangent (oriented in the plane of the direction of travel for the curve), and the tangent function of a Bezier curve is actually just another Bezier curve, 1 order lower. Let's find the normal for a cubic Bezier curve. The regular function, with (a,b,c,d) being the curve coordinates in a single dimension:

function computeBezier (t,a,b,c,d) {
  return a * (1-t)³ + 3 * b * (1-t)² * t + 3 * c * (1-t) * t² + d * t³
}

Note that Bezier curves are symmetrical, the only difference between t vs 1-t is which end of the curve represents "the start". Using a * (1-t)³ means the curve starts at a. Using a * ³t would make it start at d instead.

So let's define a quick curve with the following coordinates:

a = (-100,100,0)
b = (200,-100,100)
c = (0,100,-500)
d = (-100,-100,100)

In order to get the normal for this function, we first need the derivative:

function computeBezierDerivative (t,a,b,c,d) {
  a = 3*(b−a)
  b = 3*(c-b)
  c = 3*(d-c)
  return a * (1-t)² + 2 * b * (1-t) * t + 3 * c * t²
}

Done. Computing the derivative is stupidly simple (a fantastic property of Bezier curves).

Now, in order to get the normal, we need to take the normalised tangent vector at some value t, and rotate it by a quarter turn. We can turn it in quite a few directions, so a further restriction is that we want to turn it only in the plane that is defined by the tangent vector, and the tangent vector "right next to it", an infinitesimally small interval apart.

The tangent vector for any Bezier curve is formed simply by taking however-many-dimensions you have, and evaluating them separately, so for a 3D curve:

             | computeBezierDerivative(t, x values) |    |x'|
Tangent(t) = | computeBezierDerivative(t, y values) | => |y'|
             | computeBezierDerivative(t, z values) |    |z'|

Again, quite simple to compute. To normalise this vector (or in fact any vector), we simply perform a vector division by its length:

                   |x'|
NormalTangent(t) = |y'| divided by sqrt(x'² + y'² + z'²)
                   |z'|

So let's draw those in green:

The only trick is now to find the plane in which to rotate the tangent vector, to turn the tangent into a normal. We know we can use another t value arbitrarily close to the one we want, and turn that into a second tangent vector damn near on the same point, for finding the plane with arbitrary correctness, so we can do that:

Given an original point f(t1)=p we take a point f(t2)=q with t2=t1+e, where e is some small value like 0.001 -- this point q has a derivative q' = pointDerivative(t2), and in order to make things easier for us, we move that tangent vector a tiny bit by p-q so that the two vectors both "start" at p. Pretty simple.

However, this is equivalent to computing the first and second derivative at p, and then forming the second vector by adding those two together, as the second derivative gives us the change of the tangent at a point, so adding the second derivative vector to the first derivative vector gets us two vectors in the plane at p without having to find an adjacent point. This can be useful in curves where there are discontinuities in the derivative, i.e. curves with cusps.

We now have two vectors, departing at the same coordinate: our real tangent, and the "next" point's tangent, which is so close it might as well be the same point. Thankfully, due to how Bezier curves work, this second tangent is never the same, but slightly different, and "slightly different" is all we need: If we have two normalised vectors, starting at the same point but pointing in different directions, we can find the axis over which we need to rotate one to get the other simply by taking the cross product between them, and thus we can find the plane that they both go through.

Order matters: we compute c = tangent₂ × tangent₁, because if we compute c = tangent₁ × tangent₂ we'll be computing the rotation axis and resulting normals in the "wrong" direction. Correcting that is literally just a "take vector, multiply by -1" at the end, but why correct after the fact when we can get it right, here. Let's see those axes of rotation in blue:

Now we have everything we need: in order to turn our normalised tangent vectors into normal vectors, all we have to do is rotate them about the axes we just found by a quarter turn. If we turn them one way, we get normals, if we turn them the other, we get backfacing normals.

For arbitrary rotation about an axis in 3D, that job is perhaps laborious, but not difficult, and the quarter turns are generally special in that they greatly simplify the maths: to rotate a point over our rotation axis c, the rotation matrix turns out to be:

    |     c₁²     c₁*c₂ - c₃  c₁*c₃ + c₂ |
R = | c₁*c₂ + c₃      c₂²     c₂*c₃ - c₁ |
    | c₁*c₃ - c₂  c₂*c₃ + c₁      c₃²    |

Where the 1, 2 and 3 subscripts are really just the x, y, and z components of our vector. So that's still easy, and all that's left is to matrix-rotate our normalised tangent:

n = R * Tangent "T"

Which is:

    | T₁ * R₁₁ + T₂ * R₁₂ + T₃ * R₁₃ |    |nx|
n = | T₁ * R₂₁ + T₂ * R₂₂ + T₃ * R₂₃ | => |ny|
    | T₁ * R₃₁ + T₂ * R₃₂ + T₃ * R₃₃ |    |nz|

And we have the normal vector(s) we need. Perfect!

Except we can do better: since we're not working with arbitrary angles but with right angles, there's a significant shortcut we can use. In the same way that the vector c was perpendicular to both tangents, our normal n is perpendicular to both c and the regular tangent, so we can use the cross product a second time to find the normal:

                    |nx|
n = c × tangent₁ => |ny|
                    |nz|

This will give us exactly the same vector, with less work.

And if we want internal normals, it's the same vector, just multiply by -1:

Pretty easy once you know the tricks! And finally, because code is always useful this gist is the Processing program I used to make sure I was telling the truth.

What if the normals behave really weird?

For example, what if we're using a 3D curve but it's planar (e.g. all z coordinates at 0)? Things suddenly do horrible things. For instance, let's look at a curve with coordinates (0,0,0), (-38,260,0), (-25,541,0), and (-15,821,0):

Similarly, particularly curvy curves may yield rather twisting normals. Looking at a curve with coordinates (0,0,0), (-38,260,200), (-25,541,-200), and (-15,821,600):

In this case, we want normals that rotate and twist as little as possible, which can be found using a Rotation Minimising Frame algorithm, such as explained in section 4 or "Computation of Rotation Minimizing Frames" (Wenping Wang, Bert Jüttler, Dayue Zheng, and Yang Liu, 2008).

Implementing their 9 line algorithm takes a little more work in a normal programming language, such as Java/Processing:

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

Still, this works really well. With the note that the Frenet frame is the "standard" tangent/axis/normal frame:

VectorFrame getFrenetFrame(double t) {
  PointVector origin = get(t);
  PointVector tangent = derivative.get(t).normalise();
  PointVector normal = getNormal(t).normalise();
  return new VectorFrame(origin, tangent, normal);
}

For our planar curve, we now see perfectly behaved normals:

And in the non-planar curve, there is minimal rotation:

And finally, these normals can be uniformly reoriented by rotating all vectors around their associated tangent vectors.