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

Question

我有一些地理数据(下面的图像显示了河流的路径为红点),我想用多段三次贝塞尔曲线近似.通过对计算器等问题,在这里和这里我发现由Philip J.施耐德从"图形宝石"的算法.我成功地实现了它并且可以报告即使有数千个点它也非常快.不幸的是,速度带来了一些缺点,即装配非常不合适.请考虑以下图形:

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

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

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

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

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

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

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

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

Answer 1

我找到了满足我的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参数,以减少平滑度.通过以s编程方式查看多个参数,我可以找到符合给定要求的良好参数.

但问题是如何将该结果转换为贝塞尔曲线.Zachary Pincus 在这封电子邮件中的答案.我将在这里复制他的解决方案,以完整回答我的问题:

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;

如果您需要反向转换,请参阅:

• 贝塞尔曲线与曲线内的控制点