## 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 – Andrew Stacey Dec 15 2009 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 2009 at 15:23 Thanks lhf. Suppose I'll stick with my current solution then. – David Rutten Dec 15 2009 at 15:32

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.

Contrary to what lhf said, Cardano's formula for the cubic polynomial is perfectly stable numerically. You just need arithmetic with complex numbers even if, indeed exactly when, all three roots are real.

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.)

See also this related MO question on the multivariate situation, which you would encounter for bicubic patches in 3D. The multivariate situation is pretty much the same: You have a choice between polynomial algebra and divide-and-conquer plus Newton's method. The higher the dimension, the more justification there is for the latter over the former.

-
Perhaps I was too harsh. I meant this: linus.socs.uts.edu.au/~don/pubs/solving.html – lhf Dec 15 2009 at 19:39
That is an interesting analysis that could well be useful for high-performance graphics software. However, the author mentions in passing that you can clean up any rounding problems with a single Newton's method step at the end. Good enough for government work. – Greg Kuperberg Dec 15 2009 at 19:51