立方贝塞尔曲线上的最近点？

Question

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

Answer 1

经过大量的搜索,我发现了一篇论文,讨论了在贝塞尔曲线上找到最近点到给定点的方法:

Bezier曲线点投影的改进代数算法,由陈小貂,尹周,甄宇宇,苏华和Jean-Claude Paul撰写.

此外,我发现维基百科和MathWorld对Sturm序列的描述有助于理解算法的第一部分,因为论文本身在其描述中并不十分清楚.

Answer 2

我编写了一些快速而肮脏的代码,可以对任何程度的Bézier曲线进行估算.(注意:这是伪蛮力,而不是封闭形式的解决方案.)

演示: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()调用将是微不足道的bézierPoint().

调整算法

该closestPoint()功能分两个阶段:

• 首先,沿曲线计算点(t参数的均匀间隔值).记录t的哪个值与该点的距离最小.
• 然后,使用该localMinimum()函数捕获最小距离周围的区域,使用二分搜索找到产生真实最小距离的t和点.

scansin 的值closestPoint()确定第一遍中要使用的样本数.更少的扫描速度更快,但增加了错过真正最小点的机会.

?传递给localMinimum()函数的限制控制了它继续寻找最佳值的时间.的价值1e-2量化曲线进入〜100分,这样的话你可以看到返回的点closestPoint()沿线大跌眼镜.precision-每增加一个小数点1e-3,1e-4要额外6-8电话,... -costs bézierPoint().

Answer 3

根据您的公差.蛮力和接受错误.对于一些罕见的情况,该算法可能是错误的.但是,在大多数情况下,它会找到一个非常接近正确答案的点,结果会提高您设置切片的位置.它只是以规则的间隔沿着曲线尝试每个点并返回它找到的最佳点.

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,但是您还需要它来执行实际和正确的Bezier细分算法(还有其他算法,但这就是我想要的推荐),不仅要计算一个近似值,而是将其除以线段,而不是实际曲线.或者更确切地说,多边形船体肯定包含曲线.

您可以使用上面的算法来细分给定t处的曲线.因此T = 0.5将曲线切成两半(注意0.2将通过曲线将其切割20%80%).然后,您可以索引金字塔侧面的各个点以及从基座构建的金字塔的另一侧.所以例如立方体:

   9
  7 8
 4 5 6
0 1 2 3

您将算法0 1 2 3作为控制点,然后您将两个完美细分的曲线索引为0,4,7,9和9,8,6,3.请特别注意这些曲线的开始和结束在同一点上.并且最终索引9(曲线上的点)用作另一个新锚点.鉴于此,您可以完美地细分贝塞尔曲线.

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

但是,这需要我们实际上有一个强大的解决方案,并做一个肯定正确的算法,并正确找到曲线的一小部分,肯定是我们点的最接近点.它应该相对较快.

Answer 4

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

查看此处的交互式示例。

点击放大。

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

从贝塞尔曲线方程开始。

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)。如果您正在处理两个浮点变量，您可能需要多个方程（forx和y）：Bx(t) = (1-t)^3*px_0 + ...

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

您可以通过展开原始方程来写出四个系数。

可以使用勾股定理计算从点p到曲线B(t)的距离。

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

我现在不计算平方根，因为如果我们比较相对平方距离就足够了。以下所有等式均指平方距离。

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

二阶导数是：

一个desmos图让我们检查不同的功能。

D(t)有其局部最小值，其中d'(t) = 0且d''(t) >= 0。这意味着，我们必须为d'(t) = 0'找到t。

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

找到d'(t)的根。我使用 numpy 库，它采用多项式的系数。

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

删除虚根（只保留实根）并删除任何是< 0或 的根> 1。使用三次贝塞尔曲线，可能会剩下大约 0-3 个根。

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

返回最近的点。

我正在为 Python 中的 Bezier 附加课程。检查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] /sf/ask/3297063241/
        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]