14 KiB
Canonical form (for cubic curves)
While quadratic curves are relatively simple curves to analyze, the same cannot be said of the cubic curve. As a curvature is controlled by more than one control point, it exhibits all kinds of features like loops, cusps, odd colinear features, and as many as two inflection points because the curvature can change direction up to three times. Now, knowing what kind of curve we're dealing with means that some algorithms can be run more efficiently than if we have to implement them as generic solvers, so is there a way to determine the curve type without lots of work?
As it so happens, the answer is yes, and the solution we're going to look at was presented by Maureen C. Stone from Xerox PARC and Tony D. deRose from the University of Washington in their joint paper "A Geometric Characterization of Parametric Cubic curves". It was published in 1989, and defines curves as having a "canonical" form (i.e. a form that all curves can be reduced to) from which we can immediately tell what features a curve will have. So how does it work?
The first observation that makes things work is that if we have a cubic curve with four points, we can apply a linear transformation to these points such that three of the points end up on (0,0), (0,1) and (1,1), with the last point then being "somewhere". After applying that transformation, the location of that last point can then tell us what kind of curve we're dealing with. Specifically, we see the following breakdown:
This is a fairly funky image, so let's see how it breaks down. We see the three fixed points at (0,0), (0,1) and (1,1), and then the fourth point is somewhere. Depending on where it is, our curve will have certain features. Namely, if the fourth point is...
- anywhere on and in the red zone, the curve will either be self-intersecting (yielding a loop), or it will have a sharp discontinuity (yielding a cusp). Anywhere inside the red zone, this will be a loop. We won't know where that loop is (in terms of t values), but we are guaranteed that there is one.
- on the left (red) edge, the curve will have a cusp. We again don't know where, just that it has one. This edge is described by the function:
[ y = \frac{-x^2 + 2x + 3}{4}, { x \leq 1 } ]
- on the lower right (pink) edge, the curve will have a loop at t=1, so we know the end coordinate of the curve also lies on the curve. This edge is described by the function:
[ y = \frac{\sqrt{3(4x - x^2)} - x}{2}, { 0 \leq x \leq 1 } ]
- on the top (blue) edge, the curve will have a loop at t=0, so we know the start coordinate of the curve also lies on the curve. This edge is described by the function:
[ y = \frac{-x^2 + 3x}{3}, { x \leq 0 } ]
- inside the green zone, the curve will have a single inflection, switching concave/convex once.
- between the red and green zones, the curve has two inflections, meaning its curvature switches between concave/convex form twice.
- anywhere on the right of the red zone, the curve will have no inflections. It'll just be a well-behaved arch.
Of course, this map is fairly small, but the regions extend to infinity, with well defined boundaries.
Wait, where do those lines come from?
Without repeating the paper mentioned at the top of this section, the loop-boundaries come from rewriting the curve into canonical form, and then solving the formulae for which constraints must hold for which possible curve properties. In the paper these functions yield formulae for where you will find cusp points, or loops where we know t=0 or t=1, but those functions are derived for the full cubic expression, meaning they apply to t=-∞ to t=∞... For Bézier curves we only care about the "clipped interval" t=0 to t=1, so some of the properties that apply when you look at the curve over an infinite interval simply don't apply to the Bézier curve interval.
The right bound for the loop region, indicating where the curve switches from "having inflections" to "having a loop", for the general cubic curve, is actually mirrored over x=1, but for Bézier curves this right half doesn't apply, so we don't need to pay attention to it. Similarly, the boundaries for t=0 and t=1 loops are also nice clean curves but get "cut off" when we only look at what the general curve does over the interval t=0 to t=1.
For the full details, head over to the paper and read through sections 3 and 4. If you still remember your high school precalculus, you can probably follow along with this paper, although you might have to read it a few times before all the bits "click".
So now the question becomes: how do we manipulate our curve so that it fits this canonical form, with three fixed points, and one "free" point? Enter linear algebra. Don't worry, I'll be doing all the math for you, as well as show you what the effect is on our curves, but basically we're going to be using linear algebra, rather than calculus, because "it's way easier". Sometimes a calculus approach is very hard to work with, when the equivalent geometrical solution is super obvious.
The approach is going to start with a curve that doesn't have all-colinear points (so we need to make sure the points don't all fall on a straight line), and then applying four graphics operations that you will probably have heard of: translation (moving all points by some fixed x- and y-distance), scaling (multiplying all points by some x and y scale factor), and shearing (an operation that turns rectangles into parallelograms).
Step 1: we translate any curve by -p1.x and -p1.y, so that the curve starts at (0,0). We're going to make use of an interesting trick here, by pretending our 2D coordinates are 3D, with the z coordinate simply always being 1. This is an old trick in graphics to overcome the limitations of 2D transformations: without it, we can only turn (x,y) coordinates into new coordinates of the form (ax + by, cx + dy), which means we can't do translation, since that requires we end up with some kind of (x + a, y + b). If we add a bogus z coordinate that is always 1, then we can suddenly add arbitrary values. For example:
[ \left [ \begin{array}{ccc} 01 & 0 & a \ 0 & 1 & b \ 0 & 0 & 1 \end{array} \right ] \cdot \left [ \begin{matrix} x \ y \ z=1 \end{matrix} \right ]
\left [ \begin{matrix} 1 \cdot x + 0 \cdot y + a \cdot z \ 0 \cdot x + 1 \cdot y + b \cdot z \ 0 \cdot x + 0 \cdot y + 1 \cdot z \end{matrix} \right ]
\left [ \begin{matrix} x + a \cdot 1 \ y + b \cdot 1 \ 1 \cdot z \end{matrix} \right ]
\left [ \begin{matrix} x + a \ y + b \ z=1 \end{matrix} \right ] ]
Sweet! z stays 1, so we can effectively ignore it entirely, but we added some plain values to our x and y coordinates. So, if we want to subtract p1.x and p1.y, we use:
[ T_1 = \left [ \begin{array}{ccc} 01 & 0 & -{P_1}_x \ 0 & 1 & -{P_1}_y \ 0 & 0 & 1 \end{array} \right ] \cdot \left [ \begin{matrix} x \ y \ 1 \end{matrix} \right ]
\left [ \begin{matrix} 1 \cdot x + 0 \cdot y - {P_1}_x \cdot 1 \ 0 \cdot x + 1 \cdot y - {P_1}_y \cdot 1 \ 0 \cdot x + 0 \cdot y + 1 \cdot 1 \end{matrix} \right ]
\left [ \begin{matrix} x - {P_1}_x \ y - {P_1}_y \ 1 \end{matrix} \right ] ]
Running all our coordinates through this transformation gives a new set of coordinates, let's call those U, where the first coordinate lies on (0,0), and the rest is still somewhat free. Our next job is to make sure point 2 ends up lying on the x=0 line, so what we want is a transformation matrix that, when we run it, subtracts x from whatever x we currently have. This is called shearing, and the typical x-shear matrix and its transformation looks like this:
[ \left [ \begin{matrix} 1 & S & 0 \ 0 & 1 & 0 \ 0 & 0 & 1 \end{matrix} \right ] \cdot \left [ \begin{matrix} x \ y \ 1 \end{matrix} \right ]
\left [ \begin{matrix} x + S \cdot y \ y \ 1 \end{matrix} \right ] ]
So we want some shearing value that, when multiplied by y, yields -x, so our x coordinate becomes zero. That value is simply -x/y, because -x/y * y = -x. Done:
[ T_2 = \left [ \begin{matrix} 1 & -\frac{ {U_2}_x }{ {U_2}_y } & 0 \ 0 & 1 & 0 \ 0 & 0 & 1 \end{matrix} \right ] ]
Now, running this on all our points generates a new set of coordinates, let's call those V, which now have point 1 on (0,0) and point 2 on (0, some-value), and we wanted it at (0,1), so we need to do some scaling to make sure it ends up at (0,1). Additionally, we want point 3 to end up on (1,1), so we can also scale x to make sure its x-coordinate will be 1 after we run the transform. That means we'll be x-scaling by 1/point3x, and y-scaling by point2y. This is really easy:
[ T_3 = \left [ \begin{matrix} \frac{1}{ {V_3}_x } & 0 & 0 \ 0 & \frac{1}{ {V_2}_y } & 0 \ 0 & 0 & 1 \end{matrix} \right ] ]
Then, finally, this generates a new set of coordinates, let's call those W, of which point 1 lies on (0,0), point 2 lies on (0,1), and point three lies on (1, ...) so all that's left is to make sure point 3 ends up at (1,1) - but we can't scale! Point 2 is already in the right place, and y-scaling would move it out of (0,1) again, so our only option is to y-shear point three, just like how we x-sheared point 2 earlier. In this case, we do the same trick, but with y/x rather than x/y because we're not x-shearing but y-shearing. Additionally, we don't actually want to end up at zero (which is what we did before) so we need to shear towards an offset, in this case 1:
[ T_4 = \left [ \begin{matrix} 1 & 0 & 0 \ \frac{1 - {W_3}_y}{ {W_3}_x } & 1 & 0 \ 0 & 0 & 1 \end{matrix} \right ] ]
And this generates our final set of four coordinates. Of these, we already know that points 1 through 3 are (0,0), (0,1) and (1,1), and only the last coordinate is "free". In fact, given any four starting coordinates, the resulting "transformation mapped" coordinate will be:
[ mapped_4 = \left ( \begin{matrix} x = \left ( \frac { -x_1 + x_4 - \frac{(-x_1+x_2)(-y_1+y_4)}{-y_1+y_2} } { -x_1+x_3-\frac{(-x_1+x_2)(-y_1+y_3)}{-y_1+y_2} } \right ) \ y = \left ( \frac{(-y_1+y_4)}{-y_1+y_2} + \frac { \left ( 1 - \frac{-y_1+y_3}{-y_1+y_2} \right ) \left ( -x_1 + x_4 - \frac{(-x_1+x_2)(-y_1+y_4)}{-y_1+y_2} \right ) } { -x_1+x_3-\frac{(-x_1+x_2)(-y_1+y_3)}{-y_1+y_2} } \right ) \end{matrix} \right ) ]
That looks very complex, but notice that every coordinate value is being offset by the initial translation, and a lot of terms in there repeat: it's pretty easy to calculate this fast, since there's so much we can cache and reuse while we compute this mapped coordinate!
First, let's just do that translation step as a "preprocessing" operation so we don't have to subtract the values all the time. What does that leave?
[ ... = \left ( \begin{matrix} x = \left ( x_4 - \frac{x_2 \cdot y_4}{y_2} \middle/ x_3-\frac{x_2 \cdot y_3}{y_2} \right ) \ y = \frac{y_4}{y_2} + \left ( 1 - \frac{y_3}{y_2} \right ) \cdot \left ( x_4 - \frac{x_2 \cdot y_4}{y_2} \middle/ x_3-\frac{x_2 \cdot y_3}{y_2} \right ) \end{matrix} \right ) ]
Suddenly things look a lot simpler: the mapped x is fairly straight forward to compute, and we see that the mapped y actually contains the mapped x in its entirety, so we'll have that part already available when we need to evaluate it. In fact, let's pull out all those common factors to see just how simple this is:
[ ... = \left ( \begin{matrix} x = (x_4 - x_2 \cdot f_{42}) / ( x_3- x_2 \cdot f_{32} ) \ y = f_{42} + \left ( 1 - f_{32} \right ) \cdot x \end{matrix} \right ), f_{32} = \frac{y_3}{y_2}, f_{42} = \frac{y_4}{y_2} ]
That's kind of super-simple to write out in code, I think you'll agree. Coding math tends to be easier than the formulae initially make it look!
How do you track all that?
Doing maths can be a pain, so whenever possible, I like to make computers do the work for me. Especially for things like this, I simply use Mathematica. Tracking all this math by hand is insane, and we invented computers, literally, to do this for us. I have no reason to use pen and paper when I can write out what I want to do in a program, and have the program do the math for me. And real math, too, with symbols, not with numbers. In fact, here's the Mathematica notebook if you want to see how this works for yourself.
Now, I know, you're thinking "but Mathematica is super expensive!" and that's true, it's $295 for home use, but it's also free when you buy a $35 raspberry pi. Obviously, I bought a raspberry pi, and I encourage you to do the same. With that, as long as you know what you want to do, Mathematica can just do it for you. And we don't have to be geniuses to work out what the maths looks like. That's what we have computers for.
So, let's write up a sketch that'll show us the canonical form for any curve drawn in blue, overlaid on our canonical map, so that we can immediately tell which features our curve must have, based on where the fourth coordinate is located on the map: