Closest point on Bezier spline

Given a two-dimensional cubic bezier spline defined by 4 control-points as described here, is there a way to solve analytically the parameter along the curve (0.0 to 1.0 parameter domain) which is closest to an arbitrary point in space?

B(t) = (1-t)3 P0 + 3(1-t)2 tP1 + 3(1-t) t2P2 + t3P3, t E [0,1]

where P0, P1, P2 and P3 are the 4 control-points of the curve.

I can solve it pretty reliably and quickly with a divide-and-conquer algorithm, but it makes me feel dirty...

-- David Rutten

• You may have more luck with this kind of question in the places mentioned in the FAQ: mathoverflow.net/faq#whatnot – Loop Space Dec 15 '09 at 13:03
• A reliable and fast numerical algorithm is your best bet. Analytic solutions, even if they exist (which I doubt in this case), are not always the best choice. Even finding the roots of a cubic polynomial with Cardano's formula is messy and numerically unstable. – lhf Dec 15 '09 at 15:23
• Thanks lhf. Suppose I'll stick with my current solution then. – David Rutten Dec 15 '09 at 15:32
• If you want to see an implementation of calculating the distance from one point to a Bezier Curve(the closest point) , you can check out the "Runtime Curve Editor" assetstore.unity3d.com/en/#!/content/11835 , that's an Unity package(you perhaps need to install Unity) ,all the code is available,is C#,but math is the same, the package is doing much more than just calculating that distance(projection) , the price of package is 45$. – Varaughe Feb 28 '15 at 16:06 • Which algorithm does it use, for those of us who do not have 45$ to spare? – Federico Poloni Feb 28 '15 at 16:50

If you have a Bezier curve $(x(t),y(t))$, the closest point to the origin (say) is given by the minimum of $f(t) = x(t)^2 + y(t)^2$. By calculus, this minimum is either at the endpoints or when the derivative vanishes, $f'(t) = 0$. This latter condition is evidently a quintic polynomial. Now, there is no exact formula in radicals for solving the quintic. However, there is a really nifty new iterative algorithm based on the symmetry group of the icosahedron due to Doyle and McMullen. They make the point that you use a dynamical iteration anyway to find radicals via Newton's method; if you think of a quintic equation as a generalized radical, then it has an iteration that it just as robust numerically as finding radicals with Newton's method.
There is also a more ordinary approach to finding real roots of a quintic polynomial. (Like Cardano's formula, the Doyle-McMullen solution requires complex numbers and finds the complex roots equally easily.) Namely, you can use a cutoff procedure to switch from divide-and-conquer to Newton's method. For example, if your quintic $q(x)$ on a unit interval $[0,1]$ is $40-100x+x^5$, then it is clearly close enough to linear that Newton's method will work; you don't need divide-and-conquer. So if you have cut down the solution space to any interval, you can change the interval to $[0,1]$ (or maybe better $[-1,1]$), and then in the new variable decide whether the norms of the coefficients guarantee that Newton's method will converge. This method should only make you feel "a little dirty", because for general high-degree polynomials it's a competitive numerical algorithm. (Higher than quintic, maybe; Doyle-McMullen is really pretty good.)