当前位置：文江博客话题详情

三次贝塞尔曲线上的最近点？

发布于 2024-08-30 11:11:19 字数 40 浏览 8 评论 0原文

如何沿着三次贝塞尔曲线找到最接近平面上任意点 P 的点 B(t)？

原文

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

泪眸﹌ 2024-09-06 11:11:19

我编写了一些快速而简单的代码来估计任意阶贝塞尔曲线的这一点。 （注意：这是伪暴力，而不是封闭式解决方案。）

演示：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]);
}

上面的代码使用 vmath 库可以在向量（2D、3D 或 4D）之间有效地进行 lerp，但是替换 lerp() 就很简单了使用您自己的代码调用 bézierPoint()。

调整算法

closestPoint() 函数分两个阶段工作：

首先，计算沿曲线的所有点（t 参数的均匀间隔值）。记录哪个 t 值到该点的距离最小。
然后，使用 localMinimum() 函数寻找最小距离周围的区域，使用二分搜索来查找产生真正最小距离的 t 和点。

closestPoint() 中的 scans 值决定在第一遍中使用多少个样本。扫描次数越少速度越快，但会增加错过真正最小点的机会。

传递给 localMinimum() 函数的 ε 限制控制它继续寻找最佳值的时间。值 1e-2 将曲线量化为约 100 个点，因此您可以看到从 closestPoint() 返回的点沿着线条弹出。每个额外的小数点精度 - 1e-3、1e-4、... - 大约需要额外调用 6-8 次 bézierPoint()。

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

回复收藏 0 原文

北渚 2024-09-06 11:11:19

鉴于本页上的其他方法似乎都是近似值，因此该答案将提供一个简单的数值解。它是一个 Python 实现，依赖于 numpy 库来提供 Bezier 类。在我的测试中，这种方法的性能比我的暴力实现（使用样本和细分）好大约三倍。

请查看此处的交互式示例。

点击放大。

我使用基本代数来解决这个最小的问题。

从贝塞尔曲线方程开始。

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

由于我使用 numpy，因此我的点表示为 numpy 向量（矩阵）。这意味着p0 是一维的，例如(1, 4.2)。如果您正在处理两个浮点变量，您可能需要多个方程（对于 x 和 y）：Bx(t) = (1-t)^3* px_0 + ...

将其转换为具有四个系数的标准形式。

将原方程展开即可写出这四个系数。

从点p到曲线B(t)的距离可以使用以下公式计算毕达哥拉斯定理。

这里 a 和 b 是两个维度我们的点x和y。这意味着平方距离 D(t) 为：

我现在不计算平方根，因为如果我们比较相对平方距离。以下所有方程均指距离的平方。

该函数D(t)描述了图形和点之间的距离。我们感兴趣的是 [0, 1] 中 t 范围内的最小值。为了找到它们，我们必须对函数进行两次推导。距离函数的一阶导数是 5 阶多项式：

二阶导数是：

通过

desmos 图，我们可以检查不同的函数。

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

^{黑色：贝塞尔曲线，绿色：d(t)，紫色：d'(t)， >red:d''(t)}

求 d'(t) 的根。我使用 numpy 库，它获取多项式的系数。

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

删除虚根（仅保留实根）并删除 < 的所有根。 0 或 > 1..对于三次贝塞尔曲线，可能会剩下大约 0-3 个根。

接下来，检查根中每个 t 的每个 |B(t) - pt| 的距离。还要检查 B(0) 和 B(1) 的距离，因为贝塞尔曲线的起点和终点可能是最近的点（尽管它们不是距离函数）。

返回最近点。

我正在用 python 附加贝塞尔曲线的类。有关使用示例，请查看 github 链接。

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://stackoverflow.com/questions/47100903/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]

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://stackoverflow.com/questions/47100903/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]

回复收藏 0 原文

度的依靠╰つ 2024-09-06 11:11:19

经过大量搜索后，我找到了一篇论文，讨论了一种在贝塞尔曲线上查找与给定点最接近的点的方法：

改进的代数算法
贝塞尔曲线的投影，作者
陈小雕、周寅、舒振宇、
苏华、让-克洛德·保罗。

此外，我还找到了维基百科和MathWorld 对 Sturm 序列的描述有助于理解算法的第一部分，因为论文本身的描述不是很清楚。

回复收藏 0 原文

究竟谁懂我的在乎 2024-09-06 11:11:19

取决于你的容忍度。蛮力和接受错误。对于某些罕见的情况，该算法可能是错误的。但是，在大多数情况下，它会找到一个非常接近正确答案的点，并且将切片设置得越高，结果就会越好。它只是定期尝试沿曲线的每个点，然后返回找到的最佳点。

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;
}

只需找到最近的点并围绕该点递归，您就可以获得更好更快的结果。

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);
}

在这两种情况下，您都可以轻松地完成四边形：

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.

通过切换那里的方程。

虽然公认的答案是正确的，但你确实可以找出根源并进行比较。如果你真的只需要找到曲线上最近的点，这就可以了。

关于评论中的本。您无法在数百个控制点范围内简写公式，就像我对立方体和四边形形式所做的那样。因为每次新添加贝塞尔曲线所需的数量意味着您需要为它们建造毕达哥拉斯金字塔，而我们基本上正在处理越来越多的数字串。对于四边形，你选择 1, 2, 1，对于立方体，你选择 1, 3, 3, 1。你最终会建造越来越大的金字塔，并最终用 Casteljau 的算法将其分解，（我写这个是为了稳定的速度）：

/**
 * 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;
}

您基本上需要直接使用该算法，不仅仅是计算曲线本身上出现的 x,y，而且您还需要它来执行实际且正确的贝塞尔细分算法（还有其他算法，但这就是我想要的）推荐），不仅计算我通过将其划分为线段给出的近似值，而且还计算实际曲线。或者更确切地说，多边形外壳肯定包含曲线。

您可以通过使用上述算法在给定 t 处细分曲线来实现此目的。因此 T=0.5 将曲线切成两半（注意 0.2 会将曲线切成 20% 80%）。然后，您可以索引金字塔一侧和金字塔另一侧（从底部构建）的各个点。例如，在三次方中：

您将向算法提供 0 1 2 3 作为控制点，然后您将在 0, 4, 7, 9 和 9, 8, 6, 3 处索引两条完美细分的曲线。请特别注意查看这些曲线在同一点开始和结束。最终索引9（即曲线上的点）用作另一个新的锚点。鉴于此，您可以完美地细分贝塞尔曲线。

然后，为了找到最近的点，您需要继续将曲线细分为不同的部分，请注意，贝塞尔曲线的整个曲线都包含在控制点的外壳内。也就是说，如果我们将点 0、1、2、3 转变为连接 0,3 的闭合路径，则曲线必须完全落在该多边形外壳内。所以我们要做的就是定义给定的点 P，然后继续细分曲线，直到我们知道一条曲线的最远点比另一条曲线的最近点更近。我们简单地将点 P 与曲线的所有控制点和锚点进行比较。并丢弃活动列表中最近点（无论是锚点还是控制点）比另一条曲线的最远点更远的任何曲线。然后我们细分所有活动曲线并再次执行此操作。最终我们将得到非常细分的曲线，每一步丢弃大约一半（这意味着它应该是 O(n log n)），直到我们的错误基本上可以忽略不计。此时，我们将活动曲线称为最接近该点的点（可能有多个），并注意该高度细分的曲线位中的误差基本上等于一个点。或者简单地通过说两个锚点中最接近的那个是最接近我们的点 P 的点来确定问题。并且我们在非常具体的程度上知道误差。

然而，这要求我们实际上有一个强大的解决方案，并执行一个绝对正确的算法，并正确找到曲线的一小部分，该部分肯定是最接近我们点的点。而且它应该还是相对较快的。

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:

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.

回复收藏 0 原文

杯别 2024-09-06 11:11:19

Mike Bostock 还提供了 DOM SVG 特定最近点算法实现：

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

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

回复收藏 0 原文

温馨耳语 2024-09-06 11:11:19

该问题的解决方案是获取贝塞尔曲线上所有可能的点并比较每个距离。点数可以通过细节变量控制。

下面是在 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;
}

请注意，Bezier 类仅包含 4 个点。

可能不是最好的方法，因为根据细节它可能会变得非常慢。

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.

回复收藏 0 原文

~没有更多了~