Bézier曲线拟合与SciPy

Posted 2021-04-26

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

我最初认为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曲线无法保证通过所有点;我想要一个通过两个给定端点，并且尽可能接近指定内部点的端点。

答案

这是一段用于拟合点的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

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

另一答案

这是一种用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()

另一答案

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

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

这里有很好的教程：

Cubic Splines

Parametric Cubic Splines

另一答案

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

另一答案

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

• http://schepers.cc/getting-to-the-point
• https://github.com/DmitryBaranovskiy/raphael/blob/master/dev/raphael.core.js#L1038

另一答案

@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()

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