Closest point on a cubic Bezier curve?

I've written some quick-and-dirty code that estimates this for Bézier curves of any degree. (Note: this is pseudo-brute force, not a closed-form solution.)

Demo: http://phrogz.net/svg/closest-point-on-bezier.html

/** Find the ~closest point on a Bézier curve to a point you supply.
 * out    : A vector to modify to be the point on the curve
 * curve  : Array of vectors representing control points for a Bézier curve
 * pt     : The point (vector) you want to find out to be near
 * tmps   : Array of temporary vectors (reduces memory allocations)
 * returns: The parameter t representing the location of `out`
 */
function closestPoint(out, curve, pt, tmps) {
    let mindex, scans=25; // More scans -> better chance of being correct
    const vec=vmath['w' in curve[0]?'vec4':'z' in curve[0]?'vec3':'vec2'];
    for (let min=Infinity, i=scans+1;i--;) {
        let d2 = vec.squaredDistance(pt, bézierPoint(out, curve, i/scans, tmps));
        if (d2<min) { min=d2; mindex=i }
    }
    let t0 = Math.max((mindex-1)/scans,0);
    let t1 = Math.min((mindex+1)/scans,1);
    let d2ForT = t => vec.squaredDistance(pt, bézierPoint(out,curve,t,tmps));
    return localMinimum(t0, t1, d2ForT, 1e-4);
}

/** Find a minimum point for a bounded function. May be a local minimum.
 * minX   : the smallest input value
 * maxX   : the largest input value
 * ƒ      : a function that returns a value `y` given an `x`
 * ε      : how close in `x` the bounds must be before returning
 * returns: the `x` value that produces the smallest `y`
 */
function localMinimum(minX, maxX, ƒ, ε) {
    if (ε===undefined) ε=1e-10;
    let m=minX, n=maxX, k;
    while ((n-m)>ε) {
        k = (n+m)/2;
        if (ƒ(k-ε)<ƒ(k+ε)) n=k;
        else               m=k;
    }
    return k;
}

/** Calculate a point along a Bézier segment for a given parameter.
 * out    : A vector to modify to be the point on the curve
 * curve  : Array of vectors representing control points for a Bézier curve
 * t      : Parameter [0,1] for how far along the curve the point should be
 * tmps   : Array of temporary vectors (reduces memory allocations)
 * returns: out (the vector that was modified)
 */
function bézierPoint(out, curve, t, tmps) {
    if (curve.length<2) console.error('At least 2 control points are required');
    const vec=vmath['w' in curve[0]?'vec4':'z' in curve[0]?'vec3':'vec2'];
    if (!tmps) tmps = curve.map( pt=>vec.clone(pt) );
    else tmps.forEach( (pt,i)=>{ vec.copy(pt,curve[i]) } );
    for (var degree=curve.length-1;degree--;) {
        for (var i=0;i<=degree;++i) vec.lerp(tmps[i],tmps[i],tmps[i+1],t);
    }
    return vec.copy(out,tmps[0]);
}

The code above uses the vmath library to efficiently lerp between vectors (in 2D, 3D, or 4D), but it would be trivial to replace the lerp() call in bézierPoint() with your own code.

Tuning the Algorithm

The closestPoint() function works in two phases:

• First, calculate points all along the curve (uniformly-spaced values of the t parameter). Record which value of t has the smallest distance to the point.
• Then, use the localMinimum() function to hunt the region around the smallest distance, using a binary search to find the t and point that produces the true smallest distance.

The value of scans in closestPoint() determines how many samples to use in the first pass. Fewer scans is faster, but increases the chances of missing the true minimum point.

The ε limit passed to the localMinimum() function controls how long it continues to hunt for the best value. A value of 1e-2 quantizes the curve into ~100 points, and thus you can see the points returned from closestPoint() popping along the line. Each additional decimal point of precision—1e-3, 1e-4, …—costs about 6-8 additional calls to bézierPoint().

After lots of searching I found a paper that discusses a method for finding the closest point on a Bezier curve to a given point:

Improved Algebraic Algorithm On Point Projection For Bezier Curves, by Xiao-Diao Chen, Yin Zhou, Zhenyu Shu, Hua Su, and Jean-Claude Paul.

Furthermore, I found Wikipedia and MathWorld's descriptions of Sturm sequences useful in understanding the first part of the algoritm, as the paper itself isn't very clear in its own description.