I just added full cubic bezier curve support to the vector tessellation creation tool I am developing (see here). I’ll include some video at a later date (Update: here) showing this feature in action but first I want to document and make publicly available some of the code involved: namely, a bezier class in C# that supports two functions that are not the easiest things in the world to implement. Both of these functions come up as soon as you try to do anything non-trivial with beziers and the literature online can be unhelpful.

My code is here: link.

The first function is splitting a cubic bezier into two curves at a given *t* parameter such that the concatenation of the two curves is congruent to the original. This turns out to be easy; however, I could find no actual code on offer, rather just descriptions of code. There is good coverage elsewhere on the theory so I won’t go into it deeply here, but, briefly, splitting a bezier is easy because how one does so follows naturally from the sort of recursive definition of the bezier that is implied by De Casteljau’s algorithm. Anyway here is my code:

...
public Tuple<Bezier,Bezier> Split(double t) {
PointD E = Interpolate(t, CtrlPoint1, CtrlPoint2);
PointD F = Interpolate(t, CtrlPoint2, CtrlPoint3);
PointD G = Interpolate(t, CtrlPoint3, CtrlPoint4);
PointD H = Interpolate(t, E, F);
PointD J = Interpolate(t, F, G);
PointD K = Interpolate(t, H, J);
return new Tuple<Bezier, Bezier>(
new Bezier(CtrlPoint1, E, H, K),
new Bezier(K, J, G, CtrlPoint4)
);
}
static private PointD Interpolate(double t, PointD pt1, PointD pt2) {
double x = (1.0 - t) * pt1.X + t * pt2.X;
double y = (1.0 - t) * pt1.Y + t * pt2.Y;
return new PointD(x, y);
}
....

The other function is considerably harder: finding the nearest point on a bezier curve to a given point. If you scan the literature (i.e. perform the google search) you will find that what most people end up using is some code from an article in the book *Graphics Gems I, *“Solving the Nearest-Point-On-Curve Problem” by Philip J. Schneider. I, however, have a couple of problems with this code (1) I don’t really understand it, (2) it is long and I need to port it to C#, and (3) someone on the internet is claiming that it is incorrect but given (1) I can’t really evaluate the claim; see the comment that begins “There seem to be errors in the ControlPolygonFlatEnough function […]” in the C code here.

But more relevantly I just don’t get this code on a fundamental level. The nearest point on a bezier problem is difficult to *solve* mathematically but it is easy to *formulate* mathematically. I think that the Graphics Gems code obscures the straightforward aspect of how this problem can be formulated. Further, I don’t have a problem using code that I don’t understand; however, if I am going to port an opaque block of code from one language to another I’d like to keep that portion of code to a minimum and I can do this easily by structuring my code around the formulation of the problem that I understand.

The formulation I am talking about is as follows, if **B**(*t*) is a particular cubic bezier function, **B’**(*t*) is its derivative, and **P** is an arbitrary point then [**B**(*t*) – **P**] ⋅ **B’**(*t*) = 0 when *t *is the parameter of the closest point on **B** to **P**. [**B**(*t*) – **P**] is a vector that points from **P** towards some point on the curve and **B’**(*t*) is a vector that is tangent to the curve at this point, if the distance is going to be a minimum then these two vectors will be perpendicular and thus their dot product will be zero. Since **B**(*t*) is cubic, its derivative will be quadratic and thus when we take the scalar product of [**B**(*t*) – **P**] and **B’**(*t*) we will end up with a quintic equation of one variable. This leads naturally to a nearest point function that looks like the following on the top level:

public Tuple<PointD, double> FindNearestPoint(PointD pt)
{
var polyQuintic = GetNearestBezierPtQuintic(
_points[0].X, _points[0].Y,
_points[1].X, _points[1].Y,
_points[2].X, _points[2].Y,
_points[3].X, _points[3].Y,
pt.X, pt.Y
);
List<Complex> roots = FindAllRoots(polyQuintic);
// Filter out roots with nonzero imaginary parts and roots
// with real parts that are not between 0 and 1.
List<double> candidates = roots.FindAll(
root => root.Real > 0 && root.Real <= 1.0 && Math.Abs(root.Imaginary) < ROOT_EPS
).Select(
root => root.Real
).ToList();
// add t=0 and t=1 ... the edge cases.
candidates.Add(0.0);
candidates.Add(1.0);
// find the candidate that yields the closest point on the bezier to the given point.
double t = double.NaN;
PointD output = new PointD(double.NaN,double.NaN);
double minDistance = double.MaxValue;
foreach (double candidate in candidates)
{
var ptAtCandidate = GetPoint(candidate);
double distance = DistSqu(ptAtCandidate, pt);
if (distance < minDistance)
{
minDistance = distance;
t = candidate;
output = ptAtCandidate;
}
}
return new Tuple<PointD, double>(output, t);
}

which reduces the problem to implementing GetNearestBezierPtQuintic() and a numeric root finding function for polynomials, given that quintic equations cannot be solved via a closed-form formula like the quadratic equation.

GetNearestBezierPtQuintic() –which returns the coefficients of [**B**(*t*) – **P**] ⋅ **B’**(*t*) when fully expanded given **P** and the control points of **B **— turns out to be a little two much algebra to work out comfortably by hand, so I ran the open source symbolic mathematics application Sage on this Sage script yielding the following when translated into C#: (For the first time on the internet! … as far as I can tell)

static private List<Complex> GetNearestBezierPtQuintic(double x_0, double y_0, double x_1, double y_1,
double x_2, double y_2, double x_3, double y_3, double x, double y)
{
double t5 = 3 * x_0 * x_0 - 18 * x_0 * x_1 + 27 * x_1 * x_1 + 18 * x_0 * x_2 - 54 * x_1 * x_2 + 27 * x_2 * x_2 -
6 * x_0 * x_3 + 18 * x_1 * x_3 - 18 * x_2 * x_3 + 3 * x_3 * x_3 + 3 * y_0 * y_0 - 18 * y_0 * y_1 +
27 * y_1 * y_1 + 18 * y_0 * y_2 - 54 * y_1 * y_2 + 27 * y_2 * y_2 - 6 * y_0 * y_3 + 18 * y_1 * y_3 -
18 * y_2 * y_3 + 3 * y_3 * y_3;
double t4 = -15 * x_0 * x_0 + 75 * x_0 * x_1 - 90 * x_1 * x_1 - 60 * x_0 * x_2 + 135 * x_1 * x_2 -
45 * x_2 * x_2 + 15 * x_0 * x_3 - 30 * x_1 * x_3 + 15 * x_2 * x_3 - 15 * y_0 * y_0 + 75 * y_0 * y_1 -
90 * y_1 * y_1 - 60 * y_0 * y_2 + 135 * y_1 * y_2 - 45 * y_2 * y_2 + 15 * y_0 * y_3 - 30 * y_1 * y_3 +
15 * y_2 * y_3;
double t3 = 30 * x_0 * x_0 - 120 * x_0 * x_1 + 108 * x_1 * x_1 + 72 * x_0 * x_2 - 108 * x_1 * x_2 +
18 * x_2 * x_2 - 12 * x_0 * x_3 + 12 * x_1 * x_3 + 30 * y_0 * y_0 - 120 * y_0 * y_1 +
108 * y_1 * y_1 + 72 * y_0 * y_2 - 108 * y_1 * y_2 + 18 * y_2 * y_2 - 12 * y_0 * y_3 + 12 * y_1 * y_3;
double t2 = 3 * x * x_0 - 30 * x_0 * x_0 - 9 * x * x_1 + 90 * x_0 * x_1 - 54 * x_1 * x_1 + 9 * x * x_2 -
36 * x_0 * x_2 + 27 * x_1 * x_2 - 3 * x * x_3 + 3 * x_0 * x_3 + 3 * y * y_0 - 30 * y_0 * y_0 - 9 * y * y_1 +
90 * y_0 * y_1 - 54 * y_1 * y_1 + 9 * y * y_2 - 36 * y_0 * y_2 + 27 * y_1 * y_2 - 3 * y * y_3 + 3 * y_0 * y_3;
double t1 = -6 * x * x_0 + 15 * x_0 * x_0 + 12 * x * x_1 - 30 * x_0 * x_1 + 9 * x_1 * x_1 - 6 * x * x_2 +
6 * x_0 * x_2 - 6 * y * y_0 + 15 * y_0 * y_0 + 12 * y * y_1 - 30 * y_0 * y_1 + 9 * y_1 * y_1 -
6 * y * y_2 + 6 * y_0 * y_2;
double t0 = 3 * x * x_0 - 3 * x_0 * x_0 - 3 * x * x_1 + 3 * x_0 * x_1 + 3 * y * y_0 - 3 * y_0 * y_0 -
3 * y * y_1 + 3 * y_0 * y_1;
return new List<Complex> { (Complex)t0/t5, (Complex)t1/t5, (Complex)t2/t5, (Complex)t3/t5, (Complex)t4/t5, (Complex)1.0 };
}

and decided to use Laguerre’s Method to solve the equation.

I chose Laguerre’s Method because it is optimized for polynomials, is less flaky than Newton-Raphson, is relatively concise, and is easy to port from C++ given .Net’s System.Numerics’s implementation of complex numbers. I ported the implementation of Laguerre’s Method found in *Data Structures and Algorithms in C++* by Goodrich et. al. (here). This solver works well enough for my purposes; however, I think the optimal thing to do to get the optimal nearest-point-on-a-bezier algorithm would be to implement one of a class of newish algorithms that are specifically targeted at solving the quintic, e.g. “Solving the Quintic by Iteration” by Doyle and McMullen. No implementations of these algorithms seem to be publicly available, however, and it would be a research project that I don’t have time for right now to translate that paper from academese into code.