Impossible Bézier Calligraphy

Need for composite curves

To step back from calligraphy for a moment, consider a simpler example: drawing with a circular sponge, marker, whatever. Make it comically big in your mind.

If you draw a straight line, the composite path is very simple: the two endpoints are capped by a semicircle and are connected by a rectangle.

Now consider a curved path: you can quickly imagine that two semicircles joined by the exact path will not do the job. First of all the endpoints are wrong to connect with the endcaps, second of all the curve would look funny!

If you slow down and look at some point on the curve very closely, what points on the circle are actually doing the work of drawing? What part of the circle is extremally far from the curve at that point? The part that is tangent to the curve!

Thus we will end up offsetting each point on the curve by the radius perpendicular to the curveʼs tangent.

This thought experiment shows that we really want to find the tangent point on the pen nib that corresponds with the tangent from the pen path. If we can correlate the two for each point in time, we would get a composite path that fills out the proper area, more than the rake and stamp method.

Death

Letʼs try a simpler thing first and see why it fails: we can keep the pen path as a cubic Bézier, but restrict the pen nib to being quadratic.

Taking the tangent vector of each curve decreases degree by one: the cubic Bézier has a quadratic tangent vector, and the quadratic Bézier has a linear tangent vector. This sounds okay so far, but recall that we want the tangents to have the same slope (to be parallel). This makes us take fractions (see below for more details in the cubic case).

So the solution is a rational function (ratio of two polynomials). Bézier curves are polynomials, not rational functions, so the result will not be a Bézier curve.

Dealing with cubic curves, their tangent vector (being the derivative of their position vector) is a quadratic function. We want the two tangent vectors to be parallel, so we end up with a quadratic equation of one in terms of the other. Solving the quadratic equation introduces radicals, so it is no longer even a rational function in the cubic case.

Reckoning

So we set about approximating the exact curve by a Bézier one.

The first question to ask ourselves is: what are the tangents of the curve at the endpoints? This is a simple question, actually: since we are picking points from the curves where the tangents match, the tangent is simply what it was for the base curve. (We will work through the math below to prove why this is the case.)

If we were approximating by a quadratic curve, this would be all we need to know: the last control point would be at the intersection of the tangents from the endpoints.

But since it is cubic, we have two more points to pick, which should correspond abstractly to another parameter to control at each endpoint, in addition to the tangent there.

Youʼd think this parameter would be curvature. Youʼd think!!

2D cross product

Normally we think of it in 3D space, because the cross product of two 3D vectors is another 3D vector. But it also works in 2D space, it just produces a scalar (1D vector) instead! And it turns out to be a useful abstraction for a lot of our calculations.

$$\vec{u} \times \vec{v} = u_x v_y - u_y v_x.$$

Issues

Thereʼs a few issues we run into.

The first is that the solution doesnʼt necessarily lie on the actual Bézier segment drawn out by $0 \le q \le 1$.

Second there might be two solutions, since weʼre solving a quadratic!

The solution to both is to split things up! We need to split up the pen path so it indexes the tangents at the end of the Bézier segments of the pen nib, after first splitting the pen nib at its inflection points.

Splitting at inflection points ensures that the tangent slope is always increasing or decreasing along the segment, making there only be a single solution. Actually this requires also knowing that the Bézier segment doesnʼt rotate 180°, so we need to split it if it reaches its original tangents again.

Solving these issues means we can think of the equation above as giving us a function for $q$ in terms of $p$: $$q(p) = \frac{-b(p) \pm \sqrt{b(p)^2 - 4a(p)c(p)}}{2a(p)}.$$

This puts the functions in lock-step in terms of their tangents, giving us what we need to calculate the outside of their sweep.

Derivative of this

Weʼll need the derive of this equation soon, so letʼs calculate it while weʼre here.

My first thought was great, we have a quadratic equation, so we know the formula and can just take the derivative of it!

This was … naïve, oh so naïve. Letʼs see why.

We have our solution here:

$$q = q(p) = \frac{-b\pm\sqrt{b^2 - 4ac}}{2a}.$$

So we can take its derivative $q' = q'(p)$, using the chain rule, quotient rule, product rule … oh Iʼll spare you the gory details.

$$\begin{aligned} q' &= \frac{(-b'\pm\frac{2bb' - 4a' c - 4ac'}{\sqrt{b^2 - 4ac}})2a + 2a'(-b\pm\sqrt{b^2 - 4ac})}{4a^2}\\ &= \frac{-b'\pm\frac{2bb' - 4a' c - 4ac'}{\sqrt{b^2 - 4ac}}}{2a} + \frac{a'(-b\pm\sqrt{b^2 - 4ac})}{2a^2} \end{aligned}$$

(Recall that $a$, $b$, and $c$ are functions of $p$, so they have derivatives $a'$, $b'$, and $c'$ in terms of that variable.)

Notice any problems?

Well, first off, itʼs an ugly, messy formula! And thatʼs even with hiding the definitions of the coefficients $a$, $b$, and $c$.

The biggest problem, though, is that everything is divided by $2a$ or $4a^2$, which means it doesnʼt work when $a = 0$. That shouldnʼt be too surprising, given that the quadratic formula also fails in that case. (Itʼs the quadratic formula after all, not the quadratic-or-linear formula!)

I mean, we could solve the linear case separately:

$$bq + c = 0$$ $$q = \frac{-c}{b}$$ $$q' = \frac{-c' b + cb'}{b^2} = \frac{-c'}{b} + \frac{cb'}{b^2}$$

But that also doesnʼt work; it omits the contribution of $a'$, which does in fact influence the result of the rate of change of the quadratic formula, even when $a = 0$.I first verified this numerically, but it makes a lot of sense if you think about it.

So I took a deep breath, started texting my math advisor, and I rubber ducked myself into a much much better solution.

You see, the quadratic formula is a lie. How did we define $q$? Certainly not as a complicated quadratic formula solution. It is really defined as the implicit solution to an equation (an equation which happens to be quadratic): $$aq^2 + bq + c = 0.$$

Look, we can just take the derivative of that whole equation, even before we attempt to solve it (only takes the product rule this time!):

$$a'q^2 + 2aqq' + b'q + bq' + c' = 0.$$

And this, now this has a nicer solution: $$q' = \frac{a' q^2 + b' q + c'}{2aq + b}.$$

I think itʼs cute how the numerator is another quadratic polynomial with the derivatives of the coefficients of the original polynomial. Itʼs also convenient how we have no square roots or plusminus signs anymore – instead we write the derivative in terms of the original solution $q$.

We still have a denominator that can be zero, but this is for deeper reasons: $$2aq + b = -b \pm \sqrt{b^2 - 4ac} + b = \pm \sqrt{b^2 - 4ac} = 0.$$

Obviously this is zero exactly when $b^2 - 4ac = 0$. This quantity is called the discriminant of the quadratic, and controls much of its behavior: the basic property of how many real-valued solutions it has, as well as deeper number-theoretic properties studied in Galois theory.

Mystery coefficients

In which I attempt to write out what $a$, $b$, and $c$ actually are. Wish me luck.

We take the standard quadratic form of the tangent of $\mathbf{Q}$:

$$\mathbf{B}'_\mathbf{Q}(q) = 3(\mathbf{Q}_1 - \mathbf{Q}_0) + 6(\mathbf{Q}_2 - 2\mathbf{Q}_1 + \mathbf{Q}_0)q + 3(\mathbf{Q}_3 - 3\mathbf{Q}_2 + 3\mathbf{Q}_1 - \mathbf{Q}_0)q^2$$

(Notice how there is a constant term, a term multiplied by $q$, and a term with $q^2$.)

We want to wrangle this cross product of that with $\mathbf{B}'_\mathbf{P}(p)$ into a quadratic equation of $q$: $$\mathbf{B}'_\mathbf{Q}(q) \times \mathbf{B}'_\mathbf{P}(p) = a(p)q^2 + b(p)q + c(p)$$

So by distributivity, each coefficient we saw above is cross-producted with $\mathbf{B}'_\mathbf{P}(p)$ to obtain our mystery coefficients: $$\begin{aligned} a(p) &= 3(\mathbf{Q}_3 - 3\mathbf{Q}_2 + 3\mathbf{Q}_1 - \mathbf{Q}_0) \times \mathbf{B}'_\mathbf{P}(p),\\ b(p) &= 6(\mathbf{Q}_2 - 2\mathbf{Q}_1 + \mathbf{Q}_0) \times \mathbf{B}'_\mathbf{P}(p),\\ c(p) &= 3(\mathbf{Q}_1 - \mathbf{Q}_0) \times \mathbf{B}'_\mathbf{P}(p).\\ \end{aligned}$$

You could expand it out more into the individual components, but it would be painful, not very insightful, and waste ink.

Note that this vector cross-product cannot be cancelled out as a common term, because the $\mathbf{Q}$-vector coefficients control how the separate components of $\mathbf{B}'_\mathbf{P}(p)$ are mixed together to create the coefficients of $a$, $b$, and $c$. However, it could be divided by its norm without a problem (that is, only the direction of $\mathbf{B}'_\mathbf{P}(p)$ matters, not is magnitude – and this is by design.)

Computing compound curvature

We now know an exact formula for the idealized $\mathbf{C}_{\mathbf{P}, \mathbf{Q}}$. We can use this to get some key bits of information that will allow us to construct a good approximation to its behavior.

In particular, we want to know the slope at the endpoints and also the curvature at the endpoints. The curvature is the complicated part.

Itʼs going to get verbose very quickly, so letʼs trim down the notation a bit by leaving $p$ implicit, focusing on $t = q(p)$, and remove the extraneous parts of the Bézier notation:

$$\mathbf{C} = \mathbf{P} + \mathbf{Q}(q).$$

By construction, the slope at the endpoints is just the slope of $\mathbf{P}$ at the endpoints: $$\mathbf{C}' = \mathbf{P}' + \mathbf{Q}'(q)q' \parallel \mathbf{P}',$$ since those vectors are parallel by the construction of $q(p)$: $$\mathbf{P}' \parallel \mathbf{Q}'(q).$$

The curvature is a bit complicated, but we can work through it and just requires applying the formulas, starting with this formula for curvature of a parametric curve:

$$\mathbf{C}^\kappa = \frac{\mathbf{C}' \times \mathbf{C}''}{\|\mathbf{C}'\|^3}.$$

We already computed $\mathbf{C}'$ above, so we just need to compute $\mathbf{C}''$ and compute the cross product. $$\mathbf{C}'' = \mathbf{P}'' + \mathbf{Q}''(q)q'^2 + \mathbf{Q}'(q)q''.$$

However, I promised that we wouldnʼt need the second derivative $q''$, so letʼs see how it cancels out in the cross product. With some distributivity we can expand it:

$$\begin{aligned} \mathbf{C}' \times \mathbf{C}'' && = &\ (\mathbf{P}' + \mathbf{Q}'(q)q')\\ && \times &\ (\mathbf{P}'' + \mathbf{Q}''(q)q'^2 + \mathbf{Q}'(q)q'')\\\\ && = &\ (\mathbf{P}' + \mathbf{Q}'(q)q') \times \mathbf{P}''\\ && + &\ (\mathbf{P}' + \mathbf{Q}'(q)q') \times \mathbf{Q}''(q)q'^2\\ && + &\ \cancel{(\mathbf{P}' \times \mathbf{Q}'(q))}q''\\ && + &\ \cancel{(\mathbf{Q}'(q) \times \mathbf{Q}'(q))}q'q''\\\\ && = &\ \mathbf{C}' \times (\mathbf{P}'' + \mathbf{Q}''(q)q'^2).\\ \end{aligned}$$

Obviously $\mathbf{Q}'(q) \times \mathbf{Q}'(q) = 0$, since those are parallel, being the same vector. But $\mathbf{P}' \times \mathbf{Q}'(q) = 0$ as well, since those are parallel by construction of $q = q(p)$. That means we do not need to deal with the second derivative $q''$.

Reconstructive surgery

Now we can get down to business. How do we find the best Bézier curve to fill in for the much more complicated curve $\mathbf{C}_{\mathbf{P},\mathbf{Q}}$ that combines the two curves $\mathbf{P}$ and $\mathbf{Q}$?

Weʼll take six (6) pieces of data from $\mathbf{C}_{\mathbf{P},\mathbf{Q}}$:

1. The endpoints: $f_0$, $f_1$ (x2),
2. The tangents at the endpoints: $d_0$, $d_1$ (x2), and
3. The curvature at the endpoints: $\kappa_0$, $\kappa_1$ (x2).

This should be enough to pin down a Bézier curve, and indeed there is a way to find cubic Bézier curves that match these parameters.

We will be following the paper High accuracy geometric Hermite interpolation by Carl de Boor, Klaus Höllig, and Malcolm Sabin to answer this question. The basic sketch of the math is pretty straightforward, but the authors have done the work to come up with the right parameterizations to make it easy to compute and reason about.

The bad news is we end up with a quartic (degree 4) equation, to solve the system of two quadratic equations. So we see that there can be up to 4 solutions. But we can narrow them down a bunch, like if there are solutions with loops (self-intersections) we can rule them out, or other outlandish solutions with far-flung control points.

For example, one can take these same datapoints from a real cubic Bézier curve and reconstruct its control points from those six pieces of information. In our case, we are hoping that the curves we come across, although not technically being of that form, are very close and will still produce a similar curve to the perfect idealized solution.Note that we arenʼt doing other methods of error reduction, like minimizing area.