具有多段三次贝塞尔曲线和距离以及曲率约束的近似数据

Question

我有一些地理数据（下图将河流的路径显示为红点），我想使用多段三次贝塞尔曲线对其进行近似。通过关于 stackoverflow here 和 here 的其他问题，我从“Graphics Gems”中找到了 Philip J. Schneider 的算法。我成功地实现了它，并且可以报告即使有数千个点它也非常快。不幸的是，速度带来了一些缺点，即拟合非常草率。考虑下图：

红点是我的原始数据，蓝线是施耐德算法创建的多段贝塞尔曲线。如您所见，算法的输入是一个容差，至少与绿线所示的一样高。然而，该算法创建了一个有太多急转弯的贝塞尔曲线。您也可以在图像中看到这些不必要的急转弯。很容易想象所显示数据的贝塞尔曲线具有较少的急转弯，同时仍保持最大容差条件（只需将贝塞尔曲线稍微推向洋红色箭头的方向）。问题似乎是该算法从我的原始数据中选择数据点作为各个贝塞尔曲线的端点（洋红色箭头点表示一些嫌疑人）。由于贝塞尔曲线的端点被这样限制，很明显该算法有时会产生相当尖锐的曲率。

我正在寻找的是一种算法，它可以用具有两个约束的多段贝塞尔曲线来近似我的数据：

• 多段贝塞尔曲线与数据点的距离不得超过一定距离（这是由 Schneider 算法提供的）
• 多段贝塞尔曲线绝不能产生过于尖锐的曲率。检查此标准的一种方法是沿多段贝塞尔曲线滚动具有最小曲率半径的圆，并检查它是否沿其路径接触曲线的所有部分。虽然似乎有一个更好的方法涉及cross product of the first and second derivative

遗憾的是，我发现创建更好拟合的解决方案要么仅适用于单个贝塞尔曲线（并且省略了如何为多段贝塞尔曲线中的每条贝塞尔曲线找到良好起点和终点的问题）或不允许最小值曲率约束。我觉得最小曲率约束是这里的棘手条件。

这里是另一个例子（这是手绘的，不是 100% 精确的）：

假设图一显示了曲率约束（圆必须沿着整条曲线拟合）以及任何数据点与曲线的最大距离（恰好是绿色圆的半径） .图二中红色路径的成功近似显示为蓝色。该近似值尊重曲率条件（圆可以在整个曲线内滚动并在任何地方接触它）以及距离条件（以绿色显示）。图三显示了路径的不同近似值。虽然它尊重距离条件，但很明显圆不再适合曲率。图四显示了一条无法用给定约束近似的路径，因为它太尖了。这个例子应该说明为了正确近似路径中的一些尖角转弯，算法必须选择不属于路径的控制点。图 3 表明，如果选择了沿路径的控制点，曲率约束将不再满足。此示例还表明，算法必须在某些输入上退出，因为不可能用给定的约束来逼近它。

这个问题有解决方案吗？解决方案不必很快。如果处理 1000 个点需要一天时间，那很好。解决方案也不必是最优的，因为它必须导致最小二乘拟合。

最后我会用 C 和 Python 实现这个，但我也可以阅读大多数其他语言。

Answer 1

我找到了满足我的标准的解决方案。解决方案是首先找到一个在最小二乘意义上逼近点的 B 样条，然后将该样条转换为多段贝塞尔曲线。 B样条曲线确实具有与贝塞尔曲线相反的优点，它们不会通过控制点，并且提供了一种指定近似曲线所需“平滑度”的方法。生成这种样条所需的功能在 FITPACK 库中实现，scipy 提供了一个 python 绑定。假设我将数据读入列表x 和y，那么我可以这样做：

import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate
tck,u = interpolate.splprep([x,y],s=3)
unew = np.arange(0,1.01,0.01)
out = interpolate.splev(unew,tck)
plt.figure()
plt.plot(x,y,out[0],out[1])
plt.show()

结果如下所示：

如果我想让曲线更平滑，那么我可以将s参数增加到splprep。如果我想要更接近数据的近似值，我可以减小s 参数以降低平滑度。通过以编程方式检查多个s 参数，我可以找到一个符合给定要求的好参数。

问题是如何将该结果转换为贝塞尔曲线。 Zachary Pincus 在this email 中的答案。我将在这里复制他的解决方案，以完整回答我的问题：

def b_spline_to_bezier_series(tck, per = False):
  """Convert a parametric b-spline into a sequence of Bezier curves of the same degree.

  Inputs:
    tck : (t,c,k) tuple of b-spline knots, coefficients, and degree returned by splprep.
    per : if tck was created as a periodic spline, per *must* be true, else per *must* be false.

  Output:
    A list of Bezier curves of degree k that is equivalent to the input spline. 
    Each Bezier curve is an array of shape (k+1,d) where d is the dimension of the
    space; thus the curve includes the starting point, the k-1 internal control 
    points, and the endpoint, where each point is of d dimensions.
  """
  from fitpack import insert
  from numpy import asarray, unique, split, sum
  t,c,k = tck
  t = asarray(t)
  try:
    c[0][0]
  except:
    # I can't figure out a simple way to convert nonparametric splines to 
    # parametric splines. Oh well.
    raise TypeError("Only parametric b-splines are supported.")
  new_tck = tck
  if per:
    # ignore the leading and trailing k knots that exist to enforce periodicity 
    knots_to_consider = unique(t[k:-k])
  else:
    # the first and last k+1 knots are identical in the non-periodic case, so
    # no need to consider them when increasing the knot multiplicities below
    knots_to_consider = unique(t[k+1:-k-1])
  # For each unique knot, bring it's multiplicity up to the next multiple of k+1
  # This removes all continuity constraints between each of the original knots, 
  # creating a set of independent Bezier curves.
  desired_multiplicity = k+1
  for x in knots_to_consider:
    current_multiplicity = sum(t == x)
    remainder = current_multiplicity%desired_multiplicity
    if remainder != 0:
      # add enough knots to bring the current multiplicity up to the desired multiplicity
      number_to_insert = desired_multiplicity - remainder
      new_tck = insert(x, new_tck, number_to_insert, per)
  tt,cc,kk = new_tck
  # strip off the last k+1 knots, as they are redundant after knot insertion
  bezier_points = numpy.transpose(cc)[:-desired_multiplicity]
  if per:
    # again, ignore the leading and trailing k knots
    bezier_points = bezier_points[k:-k]
  # group the points into the desired bezier curves
  return split(bezier_points, len(bezier_points) / desired_multiplicity, axis = 0)

所以 B-Splines、FITPACK、numpy 和 scipy 拯救了我的一天 :)

Answer 2

1. 多边形化数据

   找到点的顺序，这样您就可以找到彼此最近的点并尝试将它们“按线”连接。避免循环回到原点

2. 沿路径计算推导

   这是你达到局部最小值或最大值的“线”方向的变化，那里是你的控制点......这样做可以减少你的输入数据（只留下控制点）。

3. 曲线

   现在使用这些点作为控制点。我强烈建议对 x 和 y 分别使用插值多项式，例如：

   x=a0+a1*t+a2*t*t+a3*t*t*t
   y=b0+b1*t+b2*t*t+b3*t*t*t
   

   a0..a3 的计算方式如下：

   d1=0.5*(p2.x-p0.x);
   d2=0.5*(p3.x-p1.x);
   a0=p1.x;
   a1=d1;
   a2=(3.0*(p2.x-p1.x))-(2.0*d1)-d2;
   a3=d1+d2+(2.0*(-p2.x+p1.x));
   

   • b0 .. b3 的计算方式相同，但当然使用 y 坐标
   • p0..p3 是三次插值曲线的控制点
   • t =<0.0,1.0> 是从p1 到p2 的曲线参数

   这样可以确保位置和一阶导数是连续的 (c1)，您也可以使用 BEZIER，但它不会像这样匹配。

[edit1] 边缘太锋利是个大问题

要解决它，您可以在获取控制点之前从数据集中删除点。我现在可以想到两种方法...选择更适合您的方法

1. 从数据集中删除一阶导数过高的点

   dx/dl 或 dy/dl 其中x,y 是坐标，l 是曲线长度（沿其路径）。从曲线推导中精确计算曲率半径很棘手

2. 从数据集中移除导致曲率半径过小的点

   计算相邻线段（黑线）中点的交点。像图像上的垂直轴（红线）它的距离和连接点（蓝线）是你的曲率半径。当曲率半径小于你的极限时，删除那个点......

   现在，如果您真的只需要 BEZIER 三次方，那么您可以像这样将我的插值三次方转换为 BEZIER 三次方：

//  ---------------------------------------------------------------------------
//  x=cx[0]+(t*cx[1])+(tt*cx[2])+(ttt*cx[3]); // cubic x=f(t), t = <0,1>
//  ---------------------------------------------------------------------------
//  cubic matrix                           bz4 = it4
//  ---------------------------------------------------------------------------
//  cx[0]=                            (    x0) =                    (    X1)
//  cx[1]=                   (3.0*x1)-(3.0*x0) =           (0.5*X2)         -(0.5*X0)
//  cx[2]=          (3.0*x2)-(6.0*x1)+(3.0*x0) = -(0.5*X3)+(2.0*X2)-(2.5*X1)+(    X0)
//  cx[3]= (    x3)-(3.0*x2)+(3.0*x1)-(    x0) =  (0.5*X3)-(1.5*X2)+(1.5*X1)-(0.5*X0)
//  ---------------------------------------------------------------------------
    const double m=1.0/6.0;
    double x0,y0,x1,y1,x2,y2,x3,y3;
    x0 = X1;           y0 = Y1;
    x1 = X1-(X0-X2)*m; y1 = Y1-(Y0-Y2)*m;
    x2 = X2+(X1-X3)*m; y2 = Y2+(Y1-Y3)*m;
    x3 = X2;           y3 = Y2;

如果您需要反向转换，请参阅：

• Bezier curve with control points within the curve

Answer 3

问题是很久以前发布的，但这里有一个基于 splprep 的简单解决方案，找到 s 的最小值，允许满足最小曲率半径标准。

路线是输入点的集合，第一维是点的数量。

import numpy as np
from scipy.interpolate import splprep, splev

#The minimum curvature radius we want to enforce
minCurvatureConstraint = 2000
#Relative tolerance on the radius
relTol = 1.e-6

#Initial values for bisection search, should bound the solution
s_0 = 0
minCurvature_0 = 0
s_1 = 100000000 #Should be high enough to produce curvature radius larger than constraint
s_1 *= 2
minCurvature_1 = np.float('inf')

while np.abs(minCurvature_0 - minCurvature_1)>minCurvatureConstraint*relTol:
    s = 0.5 * (s_0 + s_1)
    tck, u  = splprep(np.transpose(route), s=s)
    smoothed_route = splev(u, tck)
    
    #Compute radius of curvature
    derivative1 = splev(u, tck, der=1)
    derivative2 = splev(u, tck, der=2) 
    xprim = derivative1[0]
    xprimprim = derivative2[0]
    yprim = derivative1[1]
    yprimprim = derivative2[1]
    curvature = 1.0 / np.abs((xprim*yprimprim - yprim* xprimprim) / np.power(xprim*xprim + yprim*yprim, 3 / 2))
    minCurvature = np.min(curvature)
    print("s is %g => Minimum curvature radius is %g"%(s,np.min(curvature)))

    #Perform bisection
    if minCurvature > minCurvatureConstraint:
        s_1 = s
        minCurvature_1 = minCurvature
    else:
        s_0 = s
        minCurvature_0 = minCurvature

它可能需要一些改进（例如迭代）才能找到合适的 s_1，但可以。