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

我有一些地理数据(下面的图像显示了河流的路径为红点),我想用多段三次贝塞尔曲线近似.通过stackoverflow herehere的其他问题,我找到了来自“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救了我的一天:)

版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。

相关推荐


Python中的函数(二) 在上一篇文章中提到了Python中函数的定义和使用,在这篇文章里我们来讨论下关于函数的一些更深的话题。在学习C语言函数的时候,遇到的问题主要有形参实参的区别、参数的传递和改变、变量的作用域。同样在Python中,关于对函数的理解和使用也存在这些问题。下面来逐一讲解。一.函
Python中的字符串 可能大多数人在学习C语言的时候,最先接触的数据类型就是字符串,因为大多教程都是以"Hello world"这个程序作为入门程序,这个程序中要打印的"Hello world"就是字符串。如果你做过自然语言处理方面的研究,并且用Python
Python 面向对象编程(一) 虽然Python是解释性语言,但是它是面向对象的,能够进行对象编程。下面就来了解一下如何在Python中进行对象编程。一.如何定义一个类 在进行python面向对象编程之前,先来了解几个术语:类,类对象,实例对象,属性,函数和方法。 类是对现实世界中一些事物的封装,
Python面向对象编程(二) 在前面一篇文章中谈到了类的基本定义和使用方法,这只体现了面向对象编程的三大特点之一:封装。下面就来了解一下另外两大特征:继承和多态。 在Python中,如果需要的话,可以让一个类去继承一个类,被继承的类称为父类或者超类、也可以称作基类,继承的类称为子类。并且Pytho
Python中的函数(一) 接触过C语言的朋友对函数这个词肯定非常熟悉,无论在哪门编程语言当中,函数(当然在某些语言里称作方法,意义是相同的)都扮演着至关重要的角色。今天就来了解一下Python中的函数用法。一.函数的定义 在某些编程语言当中,函数声明和函数定义是区分开的(在这些编程语言当中函数声明
在windows下如何快速搭建web.py开发框架 用Python进行web开发的话有很多框架供选择,比如最出名的Django,tornado等,除了这些框架之外,有一个轻量级的框架使用起来也是非常方便和顺手,就是web.py。它由一名黑客所创建,但是不幸的是这位创建者于2013年自杀了。据说现在由
将Sublime Text 2搭建成一个好用的IDE 说起编辑器,可能大部分人要推荐的是Vim和Emacs,本人用过Vim,功能确实强大,但是不是很习惯,之前一直有朋友推荐SUblime Text 2这款编辑器,然后这段时间就试了一下,就深深地喜欢上这款编辑器了...
Python中的模块 有过C语言编程经验的朋友都知道在C语言中如果要引用sqrt这个函数,必须用语句"#include<math.h>"引入math.h这个头文件,否则是无法正常进行调用的。那么在Python中,如果要引用一些内置的函数,该怎么处理呢?在Python中
Python的基础语法 在对Python有了基础的认识之后,下面来了解一下Python的基础语法,看看它和C语言、java之间的基础语法差异。一.变量、表达式和语句 Python中的语句也称作命令,比如print "hello python"这就是一条语句。 表达式,顾名思义,是
Eclipse+PyDevʽjango+Mysql搭建Python web开发环境 Python的web框架有很多,目前主流的有Django、Tornado、Web.py等,最流行的要属Django了,也是被大家最看好的框架之一。下面就来讲讲如何搭建Django的开发环境。一.准备工作 需要下载的
在windows下安装配置Ulipad 今天推荐一款轻便的文本编辑器Ulipad,用来写一些小的Python脚本非常方便。 Ulipad下载地址: https://github.com/limodou/ulipad http://files.cnblogs.com/dolphin0520/u...
Python中的函数(三) 在前面两篇文章中已经探讨了函数的一些相关用法,下面一起来了解一下函数参数类型的问题。在C语言中,调用函数时必须依照函数定义时的参数个数以及类型来传递参数,否则将会发生错误,这个是严格进行规定的。然而在Python中函数参数定义和传递的方式相比而言就灵活多了。一.函数参数的
在Notepad++中搭配Python开发环境 Python在最近几年一度成为最流行的语言之一,不仅仅是因为它简洁明了,更在于它的功能之强大。它不仅能够完成一般脚本语言所能做的事情,还能很方便快捷地进行大规模的项目开发。在学习Python之前我们来看一下Python的历史由来,"Pytho
Python中的条件选择和循环语句 同C语言、Java一样,Python中也存在条件选择和循环语句,其风格和C语言、java的很类似,但是在写法和用法上还是有一些区别。今天就让我们一起来了解一下。一.条件选择语句 Python中条件选择语句的关键字为:if 、elif 、else这三个。其基本形式如
关于raw_input( )和sys.stdin.readline( )的区别 之前一直认为用raw_input( )和sys.stdin.readline( )来获取输入的效果完全相同,但是最近在写程序时有类似这样一段代码:import sysline = sys.stdin.readline()
初识Python 跟学习所有的编程语言一样,首先得了解这门语言的编程风格和最基础的语法。下面就让我们一起来了解一下Python的编程风格。1.逻辑行与物理行 在Python中有逻辑行和物理行这个概念,物理行是指在编辑器中实际看到的一行,逻辑行是指一条Python语句。在Python中提倡一个物理行只
当我们的代码是有访问网络相关的操作时,比如http请求或者访问远程数据库,经常可能会发生一些错误,有些错误可能重新去发送请求就会成功,本文分析常见可能需要重试的场景,并最后给出python代码实现。
1.经典迭代器 2.将Sentence中的__iter__改成生成器函数 改成生成器后用法不变,但更加简洁。 3.惰性实现 当列表比较大,占内存较大时,我们可以采用惰性实现,每次只读取一个元素到内存。 或者使用更简洁的生成器表达式 4.yield from itertools模块含有大量生成器函数可
本文介绍简单介绍socket的常用函数,并以python-kafka中的源码socketpair为例,来讲解python socket的运用
python实践中经常出现编码相关的异常,大多网上找资料而没有理解原理,导致一次次重复错误。本文对常用Unicode、UTF-8、GB2312编码的原理进行介绍,接着介绍了python字符类型unicode和str以及常见编解码错误UnicodeEncodeError和UnicodeDEcodeEr