如何沿着三次贝塞尔曲线找到最接近平面中任意点 P 的点 B(t)？
问问题
34934 次
6 回答
24
我已经编写了一些快速而肮脏的代码来估计任何程度的贝塞尔曲线。（注意：这是伪蛮力，不是封闭形式的解决方案。）
演示： 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和点。
scansin的值closestPoint()决定了在第一遍中使用多少样本。更少的扫描速度更快，但会增加错过真正最小点的机会。
ε传递给函数的限制localMinimum()控制它继续寻找最佳值的时间。值1e-2将曲线量化为约 100 个点，因此您可以看到从closestPoint()沿线弹出的点返回。每个额外的小数点精度 - <code>1e-3, 1e-4, ... - 花费大约 6-8 次额外调用bézierPoint().
于 2017-07-09T06:39:31.393 回答
21
经过大量搜索后，我找到了一篇论文，其中讨论了一种在贝塞尔曲线上找到离给定点最近的点的方法：
改进的 Bezier 曲线点投影代数算法，作者 Xiao-Dio Chen、Yin Zhou、Zhenyu Shu、Hua Su 和 Jean-Claude Paul。
此外，我发现Wikipedia和MathWorld对 Sturm 序列的描述有助于理解算法的第一部分，因为论文本身的描述不是很清楚。
于 2010-05-01T17:09:00.443 回答
11
取决于你的容忍度。蛮力和接受错误。对于一些罕见的情况，此算法可能是错误的。但是，在大多数情况下，它会找到一个非常接近正确答案的点，并且结果会随着您设置的切片越高而提高。它只是定期尝试沿曲线的每个点并返回它找到的最佳点。
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%）。然后，您索引金字塔一侧和金字塔另一侧的各个点，这些点是从底部构建的。因此，例如在立方中：
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,3 的闭合路径，则该曲线必须完全落在那个多边形外壳内。所以我们要做的是定义给定的点 P，然后我们继续细分曲线，直到我们知道一条曲线的最远点比另一条曲线的最近点更近。我们简单地将这个点 P 与曲线的所有控制点和锚点进行比较。并从我们的活动列表中丢弃任何最近点（无论是锚点还是控制点）比另一条曲线的最远点更远的曲线。然后我们细分所有活动曲线并再次执行此操作。最终，我们将得到非常细分的曲线，每一步丢弃大约一半（这意味着它应该是 O(n log n)），直到我们的错误基本上可以忽略不计。在这一点上，我们称我们的活动曲线为离该点最近的点（可能不止一个），并注意那条高度细分的曲线中的误差基本上等于一个点。或者简单地通过说两个锚点中最接近我们的点 P 的点来决定问题。并且我们在非常具体的程度上知道错误。
然而，这要求我们实际上有一个稳健的解决方案，并执行一个肯定正确的算法，并正确找到肯定是最接近我们点的曲线的一小部分。而且应该还是比较快的。
于 2015-12-29T23:17:10.333 回答
11
鉴于此页面上的其他方法似乎是近似值，此答案将提供一个简单的数值解决方案。它是一个 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)。如果您正在处理两个浮点变量，您可能需要多个方程（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] 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]
于 2019-08-01T18:55:13.550 回答
1
Mike Bostock还提供了最近点算法的DOM SVG 特定实现：
于 2020-08-08T19:32:11.043 回答
0
解决此问题的方法是获取贝塞尔曲线上所有可能的点并比较每个距离。点的数量可以通过细节变量来控制。
这是在 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;
}
请注意，贝塞尔类只有 4 分。
可能不是最好的方法，因为它可能会变得非常缓慢，具体取决于细节。
于 2021-10-08T12:33:12.560 回答