使用SciPy进行Bézier曲线拟合

lnvxswe2  于 2022-11-10  发布在  其他
关注(0)|答案(8)|浏览(140)

我有一组近似于2D曲线的点,我想使用Python和numpy、scipy来找到一条近似拟合这些点的三次贝塞尔曲线,其中我指定了两个端点的精确坐标,它将返回另外两个控制点的坐标。
我最初认为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第622页)。
谢谢您的建议!
编辑:我理解一条三次贝塞尔曲线并不能保证通过所有的点;我想要一个通过两个给定的端点,并尽可能接近指定的内部点。

6vl6ewon

6vl6ewon1#

以下是使用numpy创建Bezier曲线的方法:

import numpy as np
from scipy.special 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()
zwghvu4y

zwghvu4y2#

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

u7up0aaq

u7up0aaq3#

@keynesiancross要求“对[罗兰的]代码中的变量进行注解”,而其他人则完全忽略了这个问题。(以获得完美匹配),这使得我们更难理解这个问题(至少对我来说)的解决方案。与插值的区别是更容易看到的输入,留下残差。这里既有释义的代码和非贝兹输入--和一个意想不到的结果。

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

我想,失败的原因是范数测量的是曲线上的点之间的距离,而不是一条曲线上的点到另一条曲线上最近的点之间的距离。

5lhxktic

5lhxktic4#

Resulting Plot
根据@reptilicus和@Guillaume P.的回答,以下是完整的代码:

  • 获取贝塞尔参数,即从点列表中获取控制点。
  • 从贝塞尔参数(即控制点)创建贝塞尔曲线。
  • 绘制原始点、控制点和生成的贝塞尔曲线。

获取贝塞尔参数,即从一组X、Y点或坐标中获取控制点。所需的另一个参数是近似度,生成的控制点将为(度+1)

import numpy as np
from scipy.special import comb

def get_bezier_parameters(X, Y, degree=3):
    """ Least square qbezier fit using penrose pseudoinverse.

    Parameters:

    X: array of x data.
    Y: array of y data. Y[0] is the y point for X[0].
    degree: degree of the Bézier curve. 2 for quadratic, 3 for cubic.

    Based on https://stackoverflow.com/questions/12643079/b%C3%A9zier-curve-fitting-with-scipy
    and probably on the 1998 thesis by Tim Andrew Pastva, "Bézier Curve Fitting".
    """
    if degree < 1:
        raise ValueError('degree must be 1 or greater.')

    if len(X) != len(Y):
        raise ValueError('X and Y must be of the same length.')

    if len(X) < degree + 1:
        raise ValueError(f'There must be at least {degree + 1} points to '
                         f'determine the parameters of a degree {degree} curve. '
                         f'Got only {len(X)} points.')

    def bpoly(n, t, k):
        """ Bernstein polynomial when a = 0 and b = 1. """
        return t**k * (1 - t)**(n - k) * comb(n, k)
        #return comb(n, i) * ( t**(n-i) ) * (1 - t)**i

    def bmatrix(T):
        """ Bernstein matrix for Bézier curves. """
        return np.matrix([[bpoly(degree, t, k) for k in range(degree + 1)] for t in T])

    def least_square_fit(points, M):
        M_ = np.linalg.pinv(M)
        return M_ * points

    T = np.linspace(0, 1, len(X))
    M = bmatrix(T)
    points = np.array(list(zip(X, Y)))

    final = least_square_fit(points, M).tolist()
    final[0] = [X[0], Y[0]]
    final[len(final)-1] = [X[len(X)-1], Y[len(Y)-1]]
    return final

根据给定的Bezier参数(即控制点)创建Bezier曲线。

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=50):
    """
       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

使用的样本数据(可替换为任何数据,这是GPS数据)。

points = []
xpoints = [19.21270, 19.21269, 19.21268, 19.21266, 19.21264, 19.21263, 19.21261, 19.21261, 19.21264, 19.21268,19.21274, 19.21282, 19.21290, 19.21299, 19.21307, 19.21316, 19.21324, 19.21333, 19.21342]
ypoints = [-100.14895, -100.14885, -100.14875, -100.14865, -100.14855, -100.14847, -100.14840, -100.14832, -100.14827, -100.14823, -100.14818, -100.14818, -100.14818, -100.14818, -100.14819, -100.14819, -100.14819, -100.14820, -100.14820]
for i in range(len(xpoints)):
    points.append([xpoints[i],ypoints[i]])

绘制原始点、控制点和生成的贝塞尔曲线。

import matplotlib.pyplot as plt

# Plot the original points

plt.plot(xpoints, ypoints, "ro",label='Original Points')

# Get the Bezier parameters based on a degree.

data = get_bezier_parameters(xpoints, ypoints, degree=4)
x_val = [x[0] for x in data]
y_val = [x[1] for x in data]
print(data)

# Plot the control points

plt.plot(x_val,y_val,'k--o', label='Control Points')

# Plot the resulting Bezier curve

xvals, yvals = bezier_curve(data, nTimes=1000)
plt.plot(xvals, yvals, 'b-', label='B Curve')
plt.legend()
plt.show()
m528fe3b

m528fe3b5#

简短回答:你不需要,因为贝塞尔曲线不是这样工作的。更长的答案是:看一下Catmull-Rom样条。它们非常容易形成(除起点和终点外,任意点P处的切向量都平行于直线{P-1,P+1},因此它们也很容易编程),并且总是通过定义它们的点,不像Bezier曲线,它在所有控制点建立的船体内的“某个地方”插值。

sr4lhrrt

sr4lhrrt6#

Bezier曲线不能保证通过您提供给它的每个点;控制点是任意的(从没有特定的算法来找到它们的意义上说,你只是自己选择它们),并且只在一个方向上 * 拉 * 曲线。
如果你想要一条曲线通过你提供给它的每一个点,你需要像自然三次样条曲线这样的东西,由于这些曲线的限制(你必须提供递增的x坐标,否则它会趋向无穷大),你可能会需要一个参数化的自然三次样条曲线。
这里有一些不错的教程:
Cubic Splines
Parametric Cubic Splines

7y4bm7vi

7y4bm7vi7#

我遇到了和问题中详细描述的一样的问题。我使用了罗兰Puntayer提供的代码,并能够使它工作。这里:

def get_bezier_parameters(X, Y, degree=2):
    """ Least square qbezier fit using penrose pseudoinverse.

    Parameters:

    X: array of x data.
    Y: array of y data. Y[0] is the y point for X[0].
    degree: degree of the Bézier curve. 2 for quadratic, 3 for cubic.

    Based on https://stackoverflow.com/questions/12643079/b%C3%A9zier-curve-fitting-with-scipy
    and probably on the 1998 thesis by Tim Andrew Pastva, "Bézier Curve Fitting".
    """
    if degree < 1:
        raise ValueError('degree must be 1 or greater.')

    if len(X) != len(Y):
        raise ValueError('X and Y must be of the same length.')

    if len(X) < degree + 1:
        raise ValueError(f'There must be at least {degree + 1} points to '
                         f'determine the parameters of a degree {degree} curve. '
                         f'Got only {len(X)} points.')

    def bpoly(n, t, k):
        """ Bernstein polynomial when a = 0 and b = 1. """
        return t**k * (1 - t)**(n - k) * comb(n, k)

    def bmatrix(T):
        """ Bernstein matrix for Bézier curves. """
        return np.matrix([[bpoly(degree, t, k) for k in range(degree + 1)] for t in T])

    def least_square_fit(points, M):
        M_ = np.linalg.pinv(M)
        return M_ * points

    T = np.linspace(0, 1, len(X))
    M = bmatrix(T)
    points = np.array(list(zip(X, Y)))
    return least_square_fit(points, M).tolist()

要固定曲线的端点,请忽略函数返回的第一个和最后一个参数,并使用您自己的点。

xxe27gdn

xxe27gdn8#

Mike Kamermans说的是真的,但我也想指出,据我所知,catmull-rom样条可以用三次贝塞尔曲线来定义。因此,如果你只有一个处理三次贝塞尔曲线的库,你应该仍然能够处理catmull-rom样条:

相关问题