首页 文章

Bézier曲线拟合与SciPy

提问于
浏览 次
10

我有一组近似于2D曲线的点 . 我想使用带有numpy和scipy的Python来找到一个大致适合点的三次Bézier路径，其中我指定了两个 endpoints 的精确坐标，并返回其他两个控制点的坐标 .

我最初认为 scipy.interpolate.splprep() 可能会做我想要的，但它似乎迫使曲线通过每个数据点（因为我想你想要插值） . 我会假设自己走错了路 .

我的问题类似于这个：How can I fit a Bézier curve to a set of data?，除了他们说他们不想使用numpy . 我的偏好是找到我需要已经在scipy或numpy中实现的东西 . 否则，我计划使用numpy实现从该问题的答案之一链接的算法：An algorithm for automatically fitting digitized curves（pdf.page 622） .

谢谢你的任何建议！

编辑：据我所知，立方Bézier曲线无法保证通过所有点;我想要一个通过两个给定 endpoints ，并且尽可能接近指定内部点的 endpoints .

python scipy

6 回答

  • 2

    这是一段用于拟合点的python代码：

    '''least square qbezier fit using penrose pseudoinverse
    >>> V=array
    >>> E,  W,  N,  S =  V((1,0)), V((-1,0)), V((0,1)), V((0,-1))
    >>> cw = 100
    >>> ch = 300
    >>> cpb = V((0, 0))
    >>> cpe = V((cw, 0))
    >>> xys=[cpb,cpb+ch*N+E*cw/8,cpe+ch*N+E*cw/8, cpe]            
    >>> 
    >>> ts = V(range(11), dtype='float')/10
    >>> M = bezierM (ts)
    >>> points = M*xys #produces the points on the bezier curve at t in ts
    >>> 
    >>> control_points=lsqfit(points, M)
    >>> linalg.norm(control_points-xys)<10e-5
    True
    >>> control_points.tolist()[1]
    [12.500000000000037, 300.00000000000017]

'''
from numpy import array, linalg, matrix
from scipy.misc import comb as nOk
Mtk = lambda n, t, k: t**(k)*(1-t)**(n-k)*nOk(n,k)
bezierM = lambda ts: matrix([[Mtk(3,t,k) for k in range(4)] for t in ts])
def lsqfit(points,M):
    M_ = linalg.pinv(M)
    return M_ * points

    通常在贝塞尔曲线上查看Animated bezier和bezierinfo

    回复于
  • 14

    这是一种用numpy做Bezier曲线的方法：

    import numpy as np
from scipy.misc import comb

def bernstein_poly(i, n, t):
    """
     The Bernstein polynomial of n, i as a function of t
    """

    return comb(n, i) * ( t**(n-i) ) * (1 - t)**i


def bezier_curve(points, nTimes=1000):
    """
       Given a set of control points, return the
       bezier curve defined by the control points.

       points should be a list of lists, or list of tuples
       such as [ [1,1], 
                 [2,3], 
                 [4,5], ..[Xn, Yn] ]
        nTimes is the number of time steps, defaults to 1000

        See http://processingjs.nihongoresources.com/bezierinfo/
    """

    nPoints = len(points)
    xPoints = np.array([p[0] for p in points])
    yPoints = np.array([p[1] for p in points])

    t = np.linspace(0.0, 1.0, nTimes)

    polynomial_array = np.array([ bernstein_poly(i, nPoints-1, t) for i in range(0, nPoints)   ])

    xvals = np.dot(xPoints, polynomial_array)
    yvals = np.dot(yPoints, polynomial_array)

    return xvals, yvals


if __name__ == "__main__":
    from matplotlib import pyplot as plt

    nPoints = 4
    points = np.random.rand(nPoints,2)*200
    xpoints = [p[0] for p in points]
    ypoints = [p[1] for p in points]

    xvals, yvals = bezier_curve(points, nTimes=1000)
    plt.plot(xvals, yvals)
    plt.plot(xpoints, ypoints, "ro")
    for nr in range(len(points)):
        plt.text(points[nr][0], points[nr][1], nr)

    plt.show()
    回复于
  • 0

    贝塞尔曲线不能保证通过您提供的每个点;控制点是任意的（在某种意义上说，没有特定的算法可以找到它们，你只需自己选择它们）并且只在一个方向上拉动曲线 .

    如果你想要一条能够通过你提供的每一个点的曲线，你需要像天然三次样条这样的东西，并且由于它们的局限性（你必须为它们提供增加的x坐标，或者它倾向于无穷大） ，你可能想要一个参数化的自然三次样条 .

    这里有很好的教程：

    Cubic Splines

    Parametric Cubic Splines

    回复于
  • 8

    简短的回答：你没有，因为这不是Bezier曲线的工作方式 . 更长的答案：看看Catmull-Rom splines . 它们很容易形成（任何点P处的切向量，禁止开始和结束，与线{P-1，P 1}平行，所以它们也很容易编程）并且总是通过定义它们的点，与贝塞尔曲线不同，贝塞尔曲线在所有控制点设置的凸包内插入“某处” .

    回复于
  • 1

    Mike Kamermans所说的是真的，但我也想指出，据我所知，catmull-rom样条可以用立方贝塞尔数来定义 . 因此，如果您只有一个可以使用Cubic的库，那么您仍然可以使用catmull-rom样条：

    回复于
  • 0

    @keynesiancross要求“[罗兰的代码中关于变量是什么的评论”，而其他人完全错过了陈述的问题 . Roland以Bézier曲线作为输入开始（以获得完美匹配），这使得更难理解问题和（至少对我而言）解决方案 . 对于留下残差的输入，插值的差异更容易看出 . 这是释义代码和非贝塞尔输入 - 以及意想不到的结果 .

    import matplotlib.pyplot as plt
import numpy as np
from scipy.special import comb as n_over_k
Mtk = lambda n, t, k: t**k * (1-t)**(n-k) * n_over_k(n,k)
BézierCoeff = lambda ts: [[Mtk(3,t,k) for k in range(4)] for t in ts]

fcn = np.log
tPlot = np.linspace(0. ,1. , 81)
xPlot = np.linspace(0.1,2.5, 81)
tData = tPlot[0:81:10]
xData = xPlot[0:81:10]
data = np.column_stack((xData, fcn(xData))) # shapes (9,2)

Pseudoinverse = np.linalg.pinv(BézierCoeff(tData)) # (9,4) -> (4,9)
control_points = Pseudoinverse.dot(data)     # (4,9)*(9,2) -> (4,2)
Bézier = np.array(BézierCoeff(tPlot)).dot(control_points)
residuum = fcn(Bézier[:,0]) - Bézier[:,1]

fig, ax = plt.subplots()
ax.plot(xPlot, fcn(xPlot),   'r-')
ax.plot(xData, data[:,1],    'ro', label='input')
ax.plot(Bézier[:,0],
        Bézier[:,1],         'k-', label='fit')
ax.plot(xPlot, 10.*residuum, 'b-', label='10*residuum')
ax.plot(control_points[:,0],
        control_points[:,1], 'ko:', fillstyle='none')
ax.legend()
fig.show()

    这适用于 fcn = np.cos 但不适用于 log . 我有点期望拟合会使用控制点的t分量作为额外的自由度，就像拖动控制点一样：

    manual_points = np.array([[0.1,np.log(.1)],[.27,-.6],[.82,.23],[2.5,np.log(2.5)]])
Bézier = np.array(BézierCoeff(tPlot)).dot(manual_points)
residuum = fcn(Bézier[:,0]) - Bézier[:,1]

fig, ax = plt.subplots()
ax.plot(xPlot, fcn(xPlot),   'r-')
ax.plot(xData, data[:,1],    'ro', label='input')
ax.plot(Bézier[:,0],
        Bézier[:,1],         'k-', label='fit')
ax.plot(xPlot, 10.*residuum, 'b-', label='10*residuum')
ax.plot(manual_points[:,0],
        manual_points[:,1],  'ko:', fillstyle='none')
ax.legend()
fig.show()

    我猜测，失败的原因是标准测量曲线上点之间的距离，而不是一条曲线上的点与另一条曲线上最近点之间的距离 .

    回复于

相关问题