Closest point on a cubic Bezier curve?

Question

How can I find the point B(t) along a cubic Bezier curve that is closest to an arbitrary point P in the plane?

Ary

Here is a good reference with an implementation in basic: tinaja.com/glib/bezdist.pdf – Harborage

After looking at the linked PDF, I think I'm looking for something more descriptive -- more like an academic paper. As it is, I'm not sure I understand what the algorithm being described really does. – Ary

Nothing about distance, but this is a fun page to read if you are just interested in bezier curves: redpicture.com/bezier/bezier-01.html – Harborage

This pomax.github.io/bezierjs is also a relevant implementation of the closest point in javascript. I ended up researching this and the last answer to this question for a similar problem that I have. – Tradesfolk

Answer 1

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.

Ary

This looks like a great resource. Have you translated either of the 2 algorithms in there into code that can be shared? – Illa

Anyone gotten around to implementing the described algorithm? Or know of an implementation available somewhere? I read the paper, but find some parts of it very terse, such as when the the use of Sturm sequences and continued fractions is described. I'm actually unsure whether the authors used Sturm sequences in the end or not. They also don't describe how they determined the tau parameter when bisecting (Algorithm 1) (the cut-off interval length at which they start using Newton's method). – Redwood

[This is a bit late, just for clarification] The method described in the paper uses strum sequences (this is an "old" n' well studied method for root finding -probably why the paper didn't elaborate) to find the roots of the polynomial. However, you can use any method that you would like. It is only a degree-5 polynomial. I have tried strum sequences, VCA (Vincent-Collins-Akritas), and Isolator Polynomials. Finally I opted to use isolator polynomials for its simplicity and performance. citeseerx.ist.psu.edu/viewdoc/… – Syrinx

Both VCA and Sturm sequences require the polynomial to be square-free (simple roots only). Can it be shown that the projection polynomial (which is being solved to find the closest point) (p-B(t)) dot B'(t) = 0 is square-free? – Effable

Somebody has translated it into code. – Bryonbryony

@Syrinx I've tried to implement the isolator polynomials for performance reasons, but am running into this problem (math.stackexchange.com/questions/4003082/…). Did you experience the same? – Cloudscape

Answer 2

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

Densimeter

Answer 3

Seeing as the other methods on this page seem to be approximation, this answer will provide a simple numerical solution. It is a python implementation depending on the numpy library to supply Bezier class. In my tests, this approach performed about three times better than my brute-force implementation (using samples and subdivision).

Look at the interactive example here.

Click to enlarge.

I used basic algebra to solve this minimal problem.

---

Start with the bezier curve equation.

B(t) = (1 - t)^3 * p0 + 3 * (1 - t)^2 * t * p1 + 3 * (1 - t) * t^2 * p2 + t^3 * p3

Since I'm using numpy, my points are represented as numpy vectors (matrices). This means that p0 is a one-dimensional, e.g. (1, 4.2). If you are handling two floating point variables you probably need mutliple equations (for x and y): Bx(t) = (1-t)^3*px_0 + ...

Convert it to a standard form with four coefficients.

You can write the four coefficients by expanding the original equation.

The distance from a point p to the curve B(t) can be calculated using the pythagorean theorem.

Here a and b are the two dimensions of our points x and y. This means that the squared distance D(t) is:

I'm not calculating a square root just now, because it is enough if we compare relative squared distances. All following equation will refer to the squared distance.

This function D(t) describes the distance between the graph and the points. We are interested in the minima in the range of t in [0, 1]. To find them, we have to derive the function twice. The first derivative of the distance function is a 5 order polynomial:

The second derivative is:

A desmos graph let's us examine the different functions.

D(t) has its local minima where d'(t) = 0 and d''(t) >= 0. This means, that we have to find the t for d'(t) = 0'.

^{black: the bezier curve, green: d(t), purple: d'(t), red:d''(t)}

Find the roots of d'(t). I use the numpy library, which takes the coefficients of a polynomial.

dcoeffs = np.stack([da, db, dc, dd, de, df])
roots = np.roots(dcoeffs)

Remove the imaginary roots (keep only the real roots) and remove any roots which are < 0 or > 1. With a cubic bezier, there will probably be about 0-3 roots left.

Next, check the distances of each |B(t) - pt| for each t in roots. Also check the distances for B(0) and B(1) since start and end of the Bezier curve could be the closest points (although they aren't local minima of the distance function).

Return the closest point.

I am attaching the class for the Bezier in python. Check the github link for a usage example.

import numpy as np

# Bezier Class representing a CUBIC bezier defined by four
# control points.
# 
# at(t):            gets a point on the curve at t
# distance2(pt)      returns the closest distance^2 of
#                   pt and the curve
# closest(pt)       returns the point on the curve
#                   which is closest to pt
# maxes(pt)         plots the curve using matplotlib
class Bezier(object):
    exp3 = np.array([[3, 3], [2, 2], [1, 1], [0, 0]], dtype=np.float32)
    exp3_1 = np.array([[[3, 3], [2, 2], [1, 1], [0, 0]]], dtype=np.float32)
    exp4 = np.array([[4], [3], [2], [1], [0]], dtype=np.float32)
    boundaries = np.array([0, 1], dtype=np.float32)

    # Initialize the curve by assigning the control points.
    # Then create the coefficients.
    def __init__(self, points):
        assert isinstance(points, np.ndarray)
        assert points.dtype == np.float32
        self.points = points
        self.create_coefficients()
    
    # Create the coefficients of the bezier equation, bringing
    # the bezier in the form:
    # f(t) = a * t^3 + b * t^2 + c * t^1 + d
    #
    # The coefficients have the same dimensions as the control
    # points.
    def create_coefficients(self):
        points = self.points
        a = - points[0] + 3*points[1] - 3*points[2] + points[3]
        b = 3*points[0] - 6*points[1] + 3*points[2]
        c = -3*points[0] + 3*points[1]
        d = points[0]
        self.coeffs = np.stack([a, b, c, d]).reshape(-1, 4, 2)

    # Return a point on the curve at the parameter t.
    def at(self, t):
        if type(t) != np.ndarray:
            t = np.array(t)
        pts = self.coeffs * np.power(t, self.exp3_1)
        return np.sum(pts, axis = 1)

    # Return the closest DISTANCE (squared) between the point pt
    # and the curve.
    def distance2(self, pt):
        points, distances, index = self.measure_distance(pt)
        return distances[index]

    # Return the closest POINT between the point pt
    # and the curve.
    def closest(self, pt):
        points, distances, index = self.measure_distance(pt)
        return points[index]

    # Measure the distance^2 and closest point on the curve of 
    # the point pt and the curve. This is done in a few steps:
    # 1     Define the distance^2 depending on the pt. I am 
    #       using the squared distance because it is sufficient
    #       for comparing distances and doesn't have the over-
    #       head of an additional root operation.
    #       D(t) = (f(t) - pt)^2
    # 2     Get the roots of D'(t). These are the extremes of 
    #       D(t) and contain the closest points on the unclipped
    #       curve. Only keep the minima by checking if
    #       D''(roots) > 0 and discard imaginary roots.
    # 3     Calculate the distances of the pt to the minima as
    #       well as the start and end of the curve and return
    #       the index of the shortest distance.
    #
    # This desmos graph is a helpful visualization.
    # https://www.desmos.com/calculator/ktglugn1ya
    def measure_distance(self, pt):
        coeffs = self.coeffs

        # These are the coefficients of the derivatives d/dx and d/(d/dx).
        da = 6*np.sum(coeffs[0][0]*coeffs[0][0])
        db = 10*np.sum(coeffs[0][0]*coeffs[0][1])
        dc = 4*(np.sum(coeffs[0][1]*coeffs[0][1]) + 2*np.sum(coeffs[0][0]*coeffs[0][2]))
        dd = 6*(np.sum(coeffs[0][0]*(coeffs[0][3]-pt)) + np.sum(coeffs[0][1]*coeffs[0][2]))
        de = 2*(np.sum(coeffs[0][2]*coeffs[0][2])) + 4*np.sum(coeffs[0][1]*(coeffs[0][3]-pt))
        df = 2*np.sum(coeffs[0][2]*(coeffs[0][3]-pt))

        dda = 5*da
        ddb = 4*db
        ddc = 3*dc
        ddd = 2*dd
        dde = de
        dcoeffs = np.stack([da, db, dc, dd, de, df])
        ddcoeffs = np.stack([dda, ddb, ddc, ddd, dde]).reshape(-1, 1)
        
        # Calculate the real extremes, by getting the roots of the first
        # derivativ of the distance function.
        extrema = np_real_roots(dcoeffs)
        # Remove the roots which are out of bounds of the clipped range [0, 1].
        # [future reference] https://mcmap.net/q/339983/-deleting-every-3rd-element-of-a-tensor-in-tensorflow
        dd_clip = (np.sum(ddcoeffs * np.power(extrema, self.exp4)) >= 0) & (extrema > 0) & (extrema < 1)
        minima = extrema[dd_clip]

        # Add the start and end position as possible positions.
        potentials = np.concatenate((minima, self.boundaries))

        # Calculate the points at the possible parameters t and 
        # get the index of the closest
        points = self.at(potentials.reshape(-1, 1, 1))
        distances = np.sum(np.square(points - pt), axis = 1)
        index = np.argmin(distances)

        return points, distances, index


    # Point the curve to a matplotlib figure.
    # maxes         ... the axes of a matplotlib figure
    def plot(self, maxes):
        import matplotlib.path as mpath
        import matplotlib.patches as mpatches
        Path = mpath.Path
        pp1 = mpatches.PathPatch(
            Path(self.points, [Path.MOVETO, Path.CURVE4, Path.CURVE4, Path.CURVE4]),
            fc="none")#, transform=ax.transData)
        pp1.set_alpha(1)
        pp1.set_color('#00cc00')
        pp1.set_fill(False)
        pp2 = mpatches.PathPatch(
            Path(self.points, [Path.MOVETO, Path.LINETO , Path.LINETO , Path.LINETO]),
            fc="none")#, transform=ax.transData)
        pp2.set_alpha(0.2)
        pp2.set_color('#666666')
        pp2.set_fill(False)

        maxes.scatter(*zip(*self.points), s=4, c=((0, 0.8, 1, 1), (0, 1, 0.5, 0.8), (0, 1, 0.5, 0.8),
                                                  (0, 0.8, 1, 1)))
        maxes.add_patch(pp2)
        maxes.add_patch(pp1)

# Wrapper around np.roots, but only returning real
# roots and ignoring imaginary results.
def np_real_roots(coefficients, EPSILON=1e-6):
    r = np.roots(coefficients)
    return r.real[abs(r.imag) < EPSILON]

Trotskyism

The Bezier curve equation isn't right (the one after Start with the bezier curve equation.). It should be B(t) = (1 - t)^3 * p0 + 3 * (1 - t)^2 * t * p1 + 3 * (1 - t) * t^2 * p2 + t^3 * p3 – Suffumigate

It was fixed. =) – Trotskyism

Answer 4

Depending on your tolerances. Brute force and being accepting of error. This algorithm could be wrong for some rare cases. But, in the majority of them it will find a point very close to the right answer and the results will improve the higher you set the slices. It just tries each point along the curve at regular intervals and returns the best one it found.

public double getClosestPointToCubicBezier(double fx, double fy, int slices, double x0, double y0, double x1, double y1, double x2, double y2, double x3, double y3)  {
    double tick = 1d / (double) slices;
    double x;
    double y;
    double t;
    double best = 0;
    double bestDistance = Double.POSITIVE_INFINITY;
    double currentDistance;
    for (int i = 0; i <= slices; i++) {
        t = i * tick;
        //B(t) = (1-t)**3 p0 + 3(1 - t)**2 t P1 + 3(1-t)t**2 P2 + t**3 P3
        x = (1 - t) * (1 - t) * (1 - t) * x0 + 3 * (1 - t) * (1 - t) * t * x1 + 3 * (1 - t) * t * t * x2 + t * t * t * x3;
        y = (1 - t) * (1 - t) * (1 - t) * y0 + 3 * (1 - t) * (1 - t) * t * y1 + 3 * (1 - t) * t * t * y2 + t * t * t * y3;

        currentDistance = Point.distanceSq(x,y,fx,fy);
        if (currentDistance < bestDistance) {
            bestDistance = currentDistance;
            best = t;
        }
    }
    return best;
}

You can get a lot better and faster by simply finding the nearest point and recursing around that point.

public double getClosestPointToCubicBezier(double fx, double fy, int slices, int iterations, double x0, double y0, double x1, double y1, double x2, double y2, double x3, double y3) {
    return getClosestPointToCubicBezier(iterations, fx, fy, 0, 1d, slices, x0, y0, x1, y1, x2, y2, x3, y3);
}

private double getClosestPointToCubicBezier(int iterations, double fx, double fy, double start, double end, int slices, double x0, double y0, double x1, double y1, double x2, double y2, double x3, double y3) {
    if (iterations <= 0) return (start + end) / 2;
    double tick = (end - start) / (double) slices;
    double x, y, dx, dy;
    double best = 0;
    double bestDistance = Double.POSITIVE_INFINITY;
    double currentDistance;
    double t = start;
    while (t <= end) {
        //B(t) = (1-t)**3 p0 + 3(1 - t)**2 t P1 + 3(1-t)t**2 P2 + t**3 P3
        x = (1 - t) * (1 - t) * (1 - t) * x0 + 3 * (1 - t) * (1 - t) * t * x1 + 3 * (1 - t) * t * t * x2 + t * t * t * x3;
        y = (1 - t) * (1 - t) * (1 - t) * y0 + 3 * (1 - t) * (1 - t) * t * y1 + 3 * (1 - t) * t * t * y2 + t * t * t * y3;


        dx = x - fx;
        dy = y - fy;
        dx *= dx;
        dy *= dy;
        currentDistance = dx + dy;
        if (currentDistance < bestDistance) {
            bestDistance = currentDistance;
            best = t;
        }
        t += tick;
    }
    return getClosestPointToCubicBezier(iterations - 1, fx, fy, Math.max(best - tick, 0d), Math.min(best + tick, 1d), slices, x0, y0, x1, y1, x2, y2, x3, y3);
}

In both cases you can do the quad just as easily:

x = (1 - t) * (1 - t) * x0 + 2 * (1 - t) * t * x1 + t * t * x2; //quad.
y = (1 - t) * (1 - t) * y0 + 2 * (1 - t) * t * y1 + t * t * y2; //quad.

By switching out the equation there.

While the accepted answer is right, and you really can figure out the roots and compare that stuff. If you really just need to find the nearest point on the curve, this will do it.

---

In regard to Ben in the comments. You cannot short hand the formula in the many hundreds of control point range, like I did for cubic and quad forms. Because the amount demanded by each new addition of a bezier curve means that you build a Pythagorean pyramids for them, and we're basically dealing with even more and more massive strings of numbers. For quad you go 1, 2, 1, for cubic you go 1, 3, 3, 1. You end up building bigger and bigger pyramids, and end up breaking it down with Casteljau's algorithm, (I wrote this for solid speed):

/**
 * Performs deCasteljau's algorithm for a bezier curve defined by the given control points.
 *
 * A cubic for example requires four points. So it should get at least an array of 8 values
 *
 * @param controlpoints (x,y) coord list of the Bezier curve.
 * @param returnArray Array to store the solved points. (can be null)
 * @param t Amount through the curve we are looking at.
 * @return returnArray
 */
public static float[] deCasteljau(float[] controlpoints, float[] returnArray, float t) {
    int m = controlpoints.length;
    int sizeRequired = (m/2) * ((m/2) + 1);
    if (returnArray == null) returnArray = new float[sizeRequired];
    if (sizeRequired > returnArray.length) returnArray = Arrays.copyOf(controlpoints, sizeRequired); //insure capacity
    else System.arraycopy(controlpoints,0,returnArray,0,controlpoints.length);
    int index = m; //start after the control points.
    int skip = m-2; //skip if first compare is the last control point.
    for (int i = 0, s = returnArray.length - 2; i < s; i+=2) {
        if (i == skip) {
            m = m - 2;
            skip += m;
            continue;
        }
        returnArray[index++] = (t * (returnArray[i + 2] - returnArray[i])) + returnArray[i];
        returnArray[index++] = (t * (returnArray[i + 3] - returnArray[i + 1])) + returnArray[i + 1];
    }
    return returnArray;
}

You basically need to use the algorithm directly, not just for the calculation of the x,y which occur on the curve itself, but you also need it to perform actual and proper Bezier subdivision algorithm (there are others but that is what I'd recommend), to calculate not just an approximation as I give by dividing it into line segments, but of the actual curves. Or rather the polygon hull that is certain to contain the curve.

You do this by using the above algorithm to subdivide the curves at the given t. So T=0.5 to cut the curves in half (note 0.2 would cut it 20% 80% through the curve). Then you index the various points at the side of the pyramid and the other side of the pyramid as built from the base. So for example in cubic:

   9
  7 8
 4 5 6
0 1 2 3

You would feed the algorithm 0 1 2 3 as control points, then you would index the two perfectly subdivided curves at 0, 4, 7, 9 and 9, 8, 6, 3. Take special note to see that these curves start and end at the same point. and the final index 9 which is the point on the curve is used as the other new anchor point. Given this you can perfectly subdivide a bezier curve.

Then to find the closest point you'd want to keep subdividing the curve into different parts noting that it is the case that the entire curve of a bezier curve is contained within the hull of the control points. Which is to say if we turn points 0, 1, 2, 3 into a closed path connecting 0,3 that curve must fall completely within that polygon hull. So what we do is define our given point P, then we continue to subdivide curves until such time as we know that the farthest point of one curve is closer than the closest point of another curve. We simply compare this point P to all the control and anchor points of the curves. And discard any curve from our active list whose closest point (whether anchor or control) is further away than the farthest point of another curve. Then we subdivide all the active curves and do this again. Eventually we will have very subdivided curves discarding about half each step (meaning it should be O(n log n)) until our error is basically negligible. At this point we call our active curves the closest point to that point (there could be more than one), and note that the error in that highly subdivided bit of curve is basically equal to a point. Or simply decide the issue by saying whichever of the two anchor point is closest is the closest point to our point P. And we know the error to a very specific degree.

This, though, requires that we actually have a robust solution and do a certainly correct algorithm and correctly find the tiny fraction of curve that will certainly be the closest point to our point. And it should be relatively fast still.

Bogor

How does this work? What do fx and fy signify? Why are there exactly 4 other x/y coordinates and no more or less? How can I apply this to an arbitrary Beziér curve with up to hundreds of points? – Ketcham

Basically because of the Bezout theorem. There's a few somewhat problematic bits in higher orders. Just as in the fundamental theorem algebra there are as many roots to a Bezier curve (potentially changing the given direction, as there are properties. The question asked for a cubic bezier curve directly, and this is a solid approximation but if you want higher order curves, you're going to want something that is mathematically certain to produce the correct answer at least within certain bounds. fx, fy, find x, findy. – Bogor

@BenLeggiero, updated the answer to explain how you would do it fully and properly. Also, you could simply apply the decasteljau's and find a couple good anchor points on the curve and bounce around in a binary search fashion there too. But, I'm not sure how good of an answer that might be. While I am certain it will work fairly well at pretty low orders. – Bogor

Beautiful, thank you so much for the extra work and information :) – Ketcham

Answer 5

There is also DOM SVG specific implementations of the closest point algorithms from Mike Bostock:

https://bl.ocks.org/mbostock/8027637

https://bl.ocks.org/mbostock/8027835

Pamphlet

Answer 6

A solution to this problem would be to get all the possible points on the bezier curve and compare each distance. The number of points can be controlled by the detail variable.

Here is a implementation made in Unity (C#):

public Vector2 FindNearestPointOnBezier(Bezier bezier, Vector2 point)
{
    float detail = 100;
    List<Vector2> points = new List<Vector2>();

    for (float t = 0; t < 1f; t += 1f / detail)
    {
        // this function can be exchanged for any bezier curve
        points.Add(Functions.CalculateBezier(bezier.a, bezier.b, bezier.c, bezier.d, t));
    }

    Vector2 closest = Vector2.zero;
    float minDist = Mathf.Infinity;

    foreach (Vector2 p in points)
    {
        // use sqrMagnitude as it is faster
        float dist = (p - point).sqrMagnitude;

        if (dist < minDist)
        {
            minDist = dist;
            closest = p;
        }
    }

    return closest;
}

Note that the Bezier class just holds 4 points.

Probably not the best way as it can become very slow depending on the detail.

Acicula