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

我有一些地理数据(下面的图像显示了河流的路径为红点),我想用多段三次贝塞尔曲线近似.通过stackoverflow here和here的其他问题,我找到了来自“Graphics Gems”的Philip J. Schneider的算法.我成功地实现了它并且可以报告即使有数千个点它也非常快.不幸的是,速度带来了一些缺点,即装配非常不合适.请考虑以下图形：

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

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

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

我发现可以创造更好拟合的解决方案或者仅适用于单个贝塞尔曲线(并且省略了如何在多段贝塞尔曲线中找到每个贝塞尔曲线的良好起点和终点的问题)或者不允许最小曲率约束.我认为最小曲率约束是这里的棘手条件.

这是另一个例子(这是手绘而不是100％精确)：

让我们假设图1显示了两者,曲率约束(圆必须适合整个曲线)以及任何数据点与曲线的最大距离(恰好是绿色圆的半径).图2中红色路径的成功近似显示为蓝色.该近似值符合曲率条件(圆可以在整个曲线内滚动并在任何地方触摸它)以及距离条件(以绿色显示).图3显示了路径的不同近似值.虽然它符合距离条件但很明显圆圈不再适合曲率.图4显示了一条不可能用给定约束近似的路径,因为它太尖了.该示例应该说明为了正确地近似路径中的一些尖转弯,算法必须选择不属于路径的控制点.图3显示,如果选择沿路径的控制点,则不能再满足曲率约束.此示例还显示算法必须退出某些输入,因为无法使用给定的约束来近似它.

这个问题是否存在解决方案？解决方案不一定要快.如果需要一天时间来处理1000点,那就没问题了.解决方案也不必是最佳的,因为它必须导致最小二乘拟合.

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

最佳答案

我找到了满足我的criterea的解决方案.解决方案是首先找到近似于最小二乘意义的点的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参数,以减少平滑度.通过以编程方式遍历多个参数,我可以找到符合给定要求的良好参数.

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