python - 使用多段三次贝塞尔曲线和距离以及曲率约束逼近数据

Question

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



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

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

多段贝塞尔曲线与数据点的距离不得超过一定距离(这是由 Schneider 的算法提供的)

多段贝塞尔曲线绝不能产生过于尖锐的曲率。检查此标准的一种方法是沿多段贝塞尔曲线滚动具有最小曲率半径的圆，并检查它是否沿其路径接触曲线的所有部分。尽管似乎有更好的方法涉及 cross product of the first and second derivative

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

这是另一个例子(这是手绘的，不是 100% 精确的):



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

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

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

Answer 1

我找到了满足我的标准的解决方案。解决方案是首先找到一个在最小二乘意义上逼近点的 B 样条，然后将该样条转换为多段贝塞尔曲线。 B 样条确实具有与贝塞尔曲线相比的优点，它们不会通过控制点，并且提供了一种指定近似曲线所需“平滑度”的方法。生成此类样条所需的功能在 scipy 提供了 python 绑定(bind)的 FITPACK 库中实现。假设我将数据读入列表 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以编程方式参数我可以找到一个符合给定要求的好参数。

但问题是如何将该结果转换为贝塞尔曲线。答案在 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 拯救了我的一天:)