In this article I introduce an algorithm to approximate free hand drawn curves with cubic Bezier curves in real time.

The motivation

If you follow me you probably know that I'm developping a little web app to edit files from the KanjiVG database: KanjiVGEditor. Currently it's only possible to edit the existing strokes of a kanji, or to add new ones by copy pasting strokes from another kanji. A desirable new functionality would be to allow the user to draw the strokes freely. In KanjiVG, strokes are encoded with cubic Bezier curves. Then I needed a way to convert free hand drawn curves into cubic Bezier curves.

As far as I know all vector graphic capable softwares (Adobe Illustrator, Inkscape, Libre Office Draw, Gimp, Photoshop, ...) have tools to do that, it's generally called the "pen" tool. The details vary with each software but it always works as a variation on the same theme: "point the anchor, drag the handle". They also often provide another tool, generally called the "pencil" tool, which allows to convert free hand drawn curves to some kind of parametric curves (not necessarily Bezier curve). The "pencil" tool lets the user draw the whole curve, and once it's done it calculates the approximating curve from the registered positions.

I took some time to search for other methods used in vector graphic softwares outside those I already know. In Affinity Designer the approximating curve is calculated with a delay relative to the drawing. I don't understand well what's the algorithm used, but in term of user experience it gave me the impression it "pulls" the curve toward the pointer rather than drawing it. Strange... In AutoCAD one doesn't draw the curve itself but gives a direction from the last point and an approximated curve is calculated according to that direction and previous points. One can indeed find some other methods, but pen and pencil really seem to be overwhelmingly dominant over them.

The algorithm for the "pen" tool is straightforward: it edits directly the coordinates of the controls with the pointer coordinates and extends the curve with a new control at each "down" event. The algorithm for the "pencil" tool is more complex: it memorises all the registered positions from the pointer between a "down" and an "up" event, and then solves an optimisation problem to find the controls best fitting these positions.

I did quite a lot of vector graphic drawing when I worked for an archeological survey company, both on the user side (drawing relics) and the software side (programming tools to automate part of these drawings). I found the "pen" and "pencil" tools both have inconvenients.

The "pen" tool is unnatural, you're not drawing the curve you want, you're editing the controls of a Bezier to make it match the curve you want. I happen to know the "mechanic" of a Bezier curve, but those who don't and use that tool for the first time needs to get the hang of where to place controls. The "pencil" tool is more natural, but you don't see the result curve before you've completely drawn it (note that Libre Office Draw does a better job on that point, it converts piecewise so you don't wait completely up to the end to see the result). Also, the optimisation process introduces a small delay at the end of the drawing of each curve.

These are really not terrible inconvenients, if it were it wouldn't be the methods used by all major softwares. But when you spend full days using them even a small improvement can make a big difference in your workflow. Hence the motivation to try to find a better method (beside that pathologic need to always try to reinvent the wheel myself!).

The algorithm

The algorithm I was aiming for is, in loose words, "the input of a pencil with the output of a pen". More precisely, the user draws naturally the desired curve, and the result curve approximating the drawn curve is available immediately (updated after each registered position from the pointer). To ensure an immediate response I wanted a closed form solution (no optimisation algorithm), and of course the simpler to implement the better. I went through several iterations before reaching one that please me, I'm introducing here only the final one.

There are two subproblems: 1) how to calculate a Bezier curve approximating a set of points, 2) how to decide when/where to split a curve into several ones.

The first subproblem is well known and easy to solve. Given \(A\),\(B\),\(C\),\(D\) the controls of the curve, and \(t\) the parameter of the curve, the equation of the curve is \(\mathcal{B}(t)=(1-t)^3A+3(1-t)^2tB+3(1-t)t^2C+t^3D\). Given \(n\) positions \(P_i\) and their matching parameters \(t_i\), one gets a system of linear equations \(P_i=(1-t_i)^3A+3(1-t_i)^2t_iB+3(1-t_i)t_i^2C+t_i^3D\). Forcing \(A=P_0\) and \(D=P_{n-1}\) (ie the curve must pass through the first and last points), and solving per coordinate \((B_x,C_x)\) and \((B_y,C_y)\) because they are independant, it becomes two systems of linear equations with two unknows each \(MX_.=Y_.\), $$ M_i=[3(1-t_i)^2t_i,3(1-t_i)t^2_i] $$ $$ X_.=[B_.,C_.]^t $$ $$ Y_{.i}=[P_{.i}-((1-t_i)^3P_{.0}+t_i^3P_{.(n-1)}))] $$ of which an approximate solution is obtained using the pseudo-inverse \(X_.=(M^tM)^{-1}M^tY_.\).

Here the reader may think "ok, \(P\) are the registered positions, but what about \(t\)?". Well, they aren't. Just keep reading for now. The reader may also think "and what about continuity?". Well, in that case you can define \(B=k\vec{u}\) where \(\vec{u}\) is a vector colinear to the previous curve's last anchor-handle, and solve one linear system with three unknowns \((k, C_x, C_y)\). But that's not how I'll do it. Just keep reading, again.

For the second subproblem we need a measure of fidelity, how much the approximating curve matches the set of points. That measure associated with a user defined threshold gives us a condition to start a new curve: when the measure passes over the threshold the set of points can't be approximated "well enough" with one single curve then use the previous curve to approximate the previous set of points and start a new Bezier curve with a new set of points. Such a measure is easy to define, it's the sum of squared distance between the curve and the registered points \(\sum_i||\mathcal{B}(t_i)-P_i||^2\).

So far, we need to keep in memory the set of registered points (at least those for the current Bezier curve) and the approximating Bezier curves. I didn't like the idea of keeping the registered points. In my specification the approximating curve is updated each time a new registered position is added. So, the approximating curve in itself encodes all the previously registered positions, up to the approximation error. Then, what about throwing away the previously registered positions and use the approximating curve to evaluate them when we need them?

In the first subproblem \(Y_{.i}\) becomes \([\mathcal{B}_.(t_i)-((1-t_i)^3\mathcal{B}_.(0)+t_i^3\mathcal{B}_.(1)))]\), where \(t_i\) can be freely choosen. Yes we loose time computing back \(P\) but it's actually faster most of the time because we don't need to solve a huge linear system with all the registered positions. For a cubic Bezier curve 5 \(P\) are enough (I call them "supports"), while there will be almost always many more registered positions. And of course we had the new registered position as a 6th \(P\) to calculate the new approximating curve.

In the second subproblem, we could keep as described above with \(P\) reevaluated, but there is actually a better way. If \(\mathcal{B}(t)\) is the previous approximating curve and \(\mathcal{B}'(t')\) is the new approximating curve, the fidelity measurement becomes \(\int_0^1||\mathcal{B}'(t')-\mathcal{B}(t)||^2dt\), and I choose to i \(t'=\alpha t\) where \(\alpha=length(\mathcal{B})/length(\mathcal{B'})\). The length of a Bezier curve isn't calculable in a closed form, but rememember we are creating these curve step by step, almost pixel by pixel. Approximating the length of the curve from the cumulated distance between successive registered points gives us a far sufficient approximation. Then, developping the integral above we get a neat formula for the fidelity measurement: $$ E=(1-\alpha^3)A+3\alpha^3B'-3B-3\alpha^3C'+3C+\alpha^3D'-D $$ $$ F=3(\alpha^2-1)A-6\alpha^2B'+6B+3\alpha^2C'-3C $$ $$ G=3(1-\alpha)A+3\alpha B'-3B $$ $$ H=E_x^2+E_y^2 $$ $$ I=2(E_xF_x+E_yF_y) $$ $$ J=2(E_xG_x+E_yG_y)+F_x^2+F_y^2 $$ $$ K=2(F_xG_x+F_yG_y) $$ $$ L=G_x^2+G_y^2 $$ $$ \int_0^1||\mathcal{B}'(\alpha t)-\mathcal{B}(t)||^2dt=H/7+I/6+J/5+K/4+L/3 $$

Note that here we need to take care of cumulating what I call the "drift" (the result of the integral above) because we are not calculating the fidelity relative to the original registered position but relative to the previous approximation.

Now, how to manage new segments? If it's the very first one there is not much to do, we don't know anything yet about the curve except for one point. Then set all the first segment' controls to the first registered position. If it's not the first segment we can do a bit better. We know it must start at the end of the previous segment, and it must end at the new registered position. We can't presume much about the second handle, so I choose to set it at the new registered position. In the other hand, I like to impose by default some continuity for the first handle, then I choose to set it at a quarter of the distance of the new segment and colinear with previous segment \(CD\).

Last point to clarify: what to do about the continuity in that context? I've tried several ways and what gave me the best (to my personal opinion) results is to calculate the new approximating curve without worrying about continuity, and then to correct the first anchor \(A'\) (same as \(D\)) to keep C1 continuity while respecting the ratio of distance \(CD\) and \(A'B'\). First it keeps everything simple, and second there is a nice way to think about it. As the approximating curves are created step by step one segment can't include information about registered positions of the following segment. When we correct \(D\) (the last anchor of the previous segment) we reintroduced information a posteriori in the previous segment, hence a more fidel representation of the hand drawn curve.

This way of doing also tames a problem I had with other solutions where I instead forced \(B'\) according to \(CD\), which I called "the ripples problem". As I said, the approximating curve can't account for information from the next segment. Hence \(A'B'\) forced from \(CD\) may "overshoot" the registered positions in the new segment. When this happens, \(C'\) starts to overshoot in the opposite direction to compensate, which of course makes the problem even worse when we move to the following segment. It results a "ripple" effect very unpleasing. An example is shown below. In the final version of my algorithm ripples still happen when you're moving the pointer frenetically, but I don't see it under a reasonable usage, so I consider it problem solved.

Then, the algorithm becomes:

on down event at position pos: start a new segment F() reduced to a single point at pos on move event at position pos: calculate l the distance from the end of current segment to pos calculate alpha according to l and current segment's approximated length calculate 5 equally spaced supports P from F() calculate the supports' parameter from alpha add pos to the supports with its parameter set to 1 calculate the new segment F'()'s B',C' (A' is same as A and D' is pos) if we are in continuous mode and there is a previous segment: replace the previous segment's D and current and new segments' A and A' to keep C1 continuity calculate the drift between F() and F'() if the cumulated drift is over a threshold: start a new segment with A set to the current segment's D, B set to be colinear and 1/4 of A'D' length, C and D set to D' else update the controls, cumulated drift, and approximated length of the current segment

Finally we get everything we need and the pseudocode is as follow.

STRUCT FBH: lengthLastSegment: real cumulDrift: real thresholdDrift: real, default set to 1.0 flagContinuity: bool nbSegment: int ctrls: array of pair of real, the Bezier curve's controls (nbSegment*3+1) coordinates FUNCTION EventDown(fbh: FBH, pos: pair of real): fbh.lengthLastSegment = 0.0 fbh.cumulDrift = 0.0 fbh.thresholdDrift = ... user defined fbh.flagContinuity = ... user defined fbh.nbSegment = 1 fbh.ctrls = [] FOR(i=0; i<4; i+=1): fbh.ctrls[i] = [pos[0], pos[1]] FUNCTION Eval(fbh: FBH, iSegment: int, t: real): pos = [0, 0] s = 1.0 - t FOR(i=0; i<2; i+=1): pos[i] = s * s * s * fbh.ctrls[iSegment][i] + 3.0 * s * s * t * fbh.ctrls[iSegment+1][i] + 3.0 * s * t * t * fbh.ctrls[iSegment+2][i] + t * t * t * fbh.ctrls[iSegment+3][i] return pos FUNCTION CalculateBestFitting( fbh: FBH, supports: array of pair of real, ts: array of real, ctrls: array of pair of real): FOR(iCoord=0; iCoord<2; iCoord+=1): M = [[0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0]] N = [0, 0, 0, 0, 0, 0] FOR(i=0; i<6; i+=1): M[i][0] = 3.0 * (1.0 - ts[i]) * (1.0 - ts[i]) * ts[i] M[i][1] = 3.0 * (1.0 - ts[i]) * ts[i] * ts[i] N[i] = supports[i][iCoord] - ( ctrls[0][iCoord] * (1.0 - ts[i]) * (1.0 - ts[i]) * (1.0 - ts[i]) + ctrls[3][iCoord] * ts[i] * ts[i] * ts[i]) A = [[0, 0], [0, 0]] FOR(i=0; i<2; i+=1): FOR(j=0; j<2; j+=1): FOR(k=0; k<6; k+=1): A[i][j] += M[k][j] * M[k][i] det = A[0][0] * A[1][1] - A[0][1] * A[1][0] IF(abs(det) == 0): ctrls[1] = [ctrls[0][0], ctrls[0][1]] ctrls[2] = [ctrls[3][0], ctrls[3][1]] else: B = [[0, 0], [0, 0]] B[0][0] = A[1][1] / det B[0][1] = -A[0][1] / det B[1][0] = -A[1][0] / det B[1][1] = A[0][0] / det Mp = [[0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0]] FOR(i=0; i<6; i+=1): Mp[0][i] = B[0][0] * M[i][0] + B[0][1] * M[i][1] Mp[1][i] = B[1][0] * M[i][0] + B[1][1] * M[i][1] FOR(i=0; i<6; i+=1): ctrls[1][iCoord] += Mp[0][i] * N[i] ctrls[2][iCoord] += Mp[1][i] * N[i] FUNCTION CalculateDrift( fbh: FBH, iSegment: int, ctrls: array of pair of real, alpha: real): a = fbh.ctrls[iSegment] b = fbh.ctrls[iSegment+1] c = fbh.ctrls[iSegment+2] d = fbh.ctrls[iSegment+3] bp = ctrls[1] cp = ctrls[2] dp = ctrls[3] e = [0, 0] f = [0, 0] g = [0, 0] FOR(iCoord=0; iCoord<2; iCoord+=1): e[iCoord] = (-alpha * alpha * alpha + 1.0) * a[iCoord] + 3.0 * (alpha * alpha * alpha * bp[iCoord] - b[iCoord] - alpha * alpha * alpha * cp[iCoord] + c[iCoord]) + alpha * alpha * alpha * dp[iCoord] - d[iCoord] f[iCoord] = 3.0 * ((alpha * alpha - 1.0) * a[iCoord] - 2.0 * alpha * alpha * bp[iCoord] + 2.0 * b[iCoord] + alpha * alpha * cp[iCoord] - c[iCoord]) g[iCoord] = 3.0 * ((-alpha + 1.0) * a[iCoord] + alpha * bp[iCoord] - b[iCoord]) h = e[0] * e[0] + e[1] * e[1] i = 2.0 * (e[0] * f[0] + e[1] * f[1]) j = 2.0 * (e[0] * g[0] + e[1] * g[1]) + f[0] * f[0] + f[1] * f[1] k = 2.0 * (f[0] * g[0] + f[1] * g[1]) l = g[0] * g[0] + g[1] * g[1] drift = h / 7.0 + i / 6.0 + j / 5.0 + k / 4.0 + l / 3.0 return drift FUNCTION AddSegment(fbh: FBH, pos: pair of real, distNewPos: real): fbh.nbSegment += 1 fbh.lengthLastSegment = distNewPos fbh.cumulDrift = 0.0 iNewSegment = 3 * (fbh.nbSegment - 1) iPrevSegment = 3 * (fbh.nbSegment - 2) fbh.ctrls[iNewSegment+1] = [0.0, 0.0] fbh.ctrls[iNewSegment+2] = [pos[0], pos[1]] fbh.ctrls[iNewSegment+3] = [pos[0], pos[1]] u = [0, 0] FOR(i=0; i<2; i+=1): u[i] = fbh.ctrls[iPrevSegment+3][i] - fbh.ctrls[iPrevSegment+2][i] l = sqrt(u[0] * u[0] + u[1] * u[1]) IF(l != 0): FOR(i=0; i<2; i+=1): u[i] /= l FOR(i=0; i<2; i+=1): fbh.ctrls[iNewSegment+1][i] = fbh.ctrls[iNewSegment][i] + u[i] * distNewPos / 4.0 FUNCTION EventMove(fbh: FBH, pos: pair of real): iSegment = 3 * (fbh.nbSegment - 1) distNewPos = sqrt( (fbh.ctrls[iSegment+3][0] - pos[0]) * (fbh.ctrls[iSegment+3][0] - pos[0]) + (fbh.ctrls[iSegment+3][1] - pos[1]) * (fbh.ctrls[iSegment+3][1] - pos[1])) lengthNewSegment = fbh.lengthLastSegment + distNewPos ratioLength = fbh.lengthLastSegment / lengthNewSegment supports = [] ts = [] FOR(iPos=0; iPos<5; iPos+=1): t = iPos / 4.0 supports[iPos] = Eval(fbh, iSegment, t) ts[iPos] = t * ratioLength supports[5] = [pos[0], pos[1]] ts[5] = 1.0 ctrls = [ [fbh.ctrls[iSegment][0],fbh.ctrls[iSegment][1]], [0.0,0.0], [0.0,0.0], [pos[0],pos[1]]] CalculateBestFitting(fbh, supports, ts, ctrls) IF(fbh.flagContinuity && fbh.nbSegment > 1): iPrevSegment = iSegment - 3 dA = sqrt( (fbh.ctrls[iPrevSegment+3][0] - fbh.ctrls[iPrevSegment+2][0]) * (fbh.ctrls[iPrevSegment+3][0] - fbh.ctrls[iPrevSegment+2][0]) + (fbh.ctrls[iPrevSegment+3][1] - fbh.ctrls[iPrevSegment+2][1]) * (fbh.ctrls[iPrevSegment+3][1] - fbh.ctrls[iPrevSegment+2][1])) dB = sqrt( (ctrls[1][0] - ctrls[0][0]) * (ctrls[1][0] - ctrls[0][0]) + (ctrls[1][1] - ctrls[0][1]) * (ctrls[1][1] - ctrls[0][1])) t = dA / (dA + dB) FOR(i=0; i<2; i+=1): fbh.ctrls[iPrevSegment+3][i] = fbh.ctrls[iPrevSegment+2][i] * (1.0 - t) + ctrls[1][i] * t ctrls[0][i] = fbh.ctrls[iPrevSegment+3][i] drift = CalculateDrift(fbh, iSegment, ctrls, ratioLength) IF(fbh.cumulDrift + drift < fbh.thresholdDrift): FOR(i=0; i<4; i+=1): fbh.ctrls[iSegment+i] = [ctrls[i][0], ctrls[i][1]] fbh.lengthLastSegment = lengthNewSegment fbh.cumulDrift += drift else: AddSegment(fbh, pos, distNewPos)

The results

To test my algorithm I've implemented it in C and made a small demo which allows me to trace automatically a spiral (to have repeatable tests) or draw freely. The image below shows the result of the automatic spiral for 3 threshold values and C1 continuity turned on/off. Click to enlarge, grey dots are registered pointer positions, blue lines are the approximating curves, red lines are the anchor-handle segments.

As expected, when the threshold increases the number of controls decreases but the fidelity to the drawn curve decreases. If C1 continuity is turned on, it is indeed respected, but the fidelity decreases slightly (obviously it's harder to find a fitting curve with more constraints). Controls are well spaced, and it happens that a threshold value of 1 seems to be a good value. Very neat and not expected at all!

The image below shows the result for free hand drawing, also with varying threshold and continuity.

And a video recording:

I personnally find it pleasant to use and I'm happy with the result. However it's a subjective judgement and one would have to try it to know if one prefers it over the usual algorithms. Then I've also implemented it in Javascript and integrated below so everyone can try. The range controls the threshold for the fidelity (from 0.1 to 2.0), the checkbox turns on/off C1. It's a minimalist demo, only one curve is allowed, but you can draw as many time you want (the previous curve is deleted). Have fun!

threshold: C1:

The code is also available here, in C and Javascript, each with a ready to use demo.