如何绘制多段闭合三次贝塞尔曲线(参数化)

问题描述 投票:0回答:1

在我的作业中,我尝试用Python为过山车创建一条封闭的三次贝塞尔曲线。我在代码中使用三次贝塞尔公式,如下所述: https://en.wikipedia.org/wiki/B%C3%A9zier_curve#:~:text=Cubic%20B%C3%A9zier%20curve%20with%20four,that%20can%20be%20scaled%20无限期。 我有支持点(P0,P3),但没有控制点(P1,P2),并且知道控制点应该通过线性方程和高斯算法获得,但因为找不到将它们植入为代码的方法,所以操纵p0、p3中任意段做一些控制点。我知道这是错误的,但这样做至少是为了检查代码的其他部分。现在有一条曲线与我在 CAD 应用程序(rhino3d)中通过精确点植入制作的曲线不相似。我想要它是因为控制点的实现或者我的可视化有一些错误。 干杯。

import numpy as np
import matplotlib.pyplot as plt

with open("C:\\Users\poori\Desktop\support points.txt") as file:
    text = file.read()#.split("Track")[1][4:]

lines = text.split('\n')
non_empty_lines = [line for line in lines if line.strip() != '']
text = ('\n'.join(non_empty_lines))
text_list = text.replace(" ", ",").splitlines()
for i,y in enumerate(text_list):
    text_list[i]=[float(text_list[i].split(",")[0]) , float(text_list[i].split(",")[1]), float(text_list[i].split(",")[2])]

#g = [[sum(a * b for a, b in zip(aa_row, bb_col)) for bb_col in zip(*bb)] for aa_row in bb]
s = len(text_list)-1

x = np.array([])
y = np.array([])
z = np.array([])
p1x = np.array([])
p1y = np.array([])
p1z = np.array([])
p2x = np.array([])
p2y = np.array([])
p2z = np.array([])

for i in text_list:
    x = np.append(x, np.array(i[0]))
    y = np.append(y, np.array(i[1]))
    z = np.append(z, np.array(i[2]))

for i in range(len(x)):
    if i == len(x)-1:
        break
    if (i + 2) == len(x):
        p1x = np.append(p1x, np.array((x[i] + x[i + 1] + x[1]) / 6))
        p1y = np.append(p1y, np.array((y[i] + y[i + 1] + y[1]) / 6))
        p1z = np.append(p1z, np.array((z[i] + z[i + 1] + z[1]) / 6))
        p2x = np.append(p2x, np.array((x[i] + x[i + 1] + x[1]) / 3))
        p2y = np.append(p2y, np.array((y[i] + y[i + 1] + y[1]) / 3))
        p2z = np.append(p2z, np.array((z[i] + z[i + 1] + z[1]) / 3))
    else:
        p1x = np.append(p1x, np.array((x[i] + x[i + 1] + x[i + 2]) / 6))
        print(x[i])
        p1y = np.append(p1y, np.array((y[i] + y[i + 1] + y[i + 2]) / 6))
        p1z = np.append(p1z, np.array((z[i] + z[i + 1] + z[i + 2]) / 6))
        p2x = np.append(p2x, np.array((x[i] + x[i + 1] + x[i + 2]) / 3))
        p2y = np.append(p2y, np.array((y[i] + y[i + 1] + y[i + 2]) / 3))
        p2z = np.append(p2z, np.array((z[i] + z[i + 1] + z[i + 2]) / 3))


s = np.linspace(0, 1, 300, endpoint = False)
Bx = np.array([])
By = np.array([])
Bz = np.array([])

for v in range(len(x)-1):
    for i in s:
        B = (((1-i)**3) * (np.array((x[v],y[v],z[v])))) + (3* ((1-i)**2)* i * (np.array((p1x[v],p1y[v],p1z[v])))) + (3 * (1-i) * (i**2) * (np.array((p2x[v],p2y[v],p2z[v])))) + ((i**3)*(np.array((x[v+1],y[v+1],z[v+1]))))
        Bx = np.append(Bx , B[0])
        By = np.append(By, B[1])
        Bz = np.append(Bz, B[2])


CELLS = 200
nCPTS = np.size(x, 0)
n = nCPTS -1
i = 0
t = np.linspace(0,1, CELLS)
b = []

xBezier = np.zeros((1, CELLS))
yBezier = np.zeros((1, CELLS))
zBezier = np.zeros((1, CELLS))

plt.figure().add_subplot(projection='3d')
plt.plot(Bx,By,Bz)#,out[0],out[1])
plt.show()

#support points(needs to find control points by these)
#0, 0, 0
#0,-20.836,0
#0,-38.948,0
#14.277,-38.948,0
#30.064,-38.948,0
#47.36,-38.948,0
#63.942,-38.948,3
#78.327,-38.948,6
#78.327,-30.509,9
#78.327,-22.191,12
#64.542,-22.191,15
#47.662,-22.191,18
#25.768,-22.191,21
#25.768,0.98532,24
#25.768,13.24,27.261
#25.768,23.53,30
#25.768,30.216,33
#42.9,30.216,33
#57.202,30.216,33
#75.909,30.216,33
#75.909,-0.99106,33
#75.909,-2.4834,40
#75.909,-2.6028,49
#76.909,3.2814,53
#77.909,8.2948,53
#78.909,14.468,49
#78.909,14.307,40
#78.909,11.782,33
#78.909,-15.121,36
#60.112,-15.121,39
#43.57,-15.121,42
#27.405,-15.121,48
#15.597,-15.121,48
#15.597,-3.0904,48
#15.597,11.195,48
#15.597,27.638,48
#26.208,27.638,48
#34.6,27.638,48
#34.6,14.052,48
#34.6,3.8624,48
#34.6,-5.5279,45
#44.19,-5.5279,42
#57.776,-5.5279,39
#67.565,-5.5279,33
#67.565,2.9199,27
#67.565,11.511,24
#67.565,19.902,21
#67.565,27.894,18
#67.565,39.226,15
#54.779,39.226,12
#43.924,39.226,9
#30.976,39.226,6
#18.356,39.226,3
#0,39.226,0
#0,17.662,0
#0,0,0
python numpy matplotlib bezier
1个回答
0
投票

正如我之前在评论中所写,对于平滑的复合贝塞尔曲线,您需要确保每个贝塞尔曲线末尾的梯度是连续的。本质上,您需要在此处估计切线方向,并选择 P1 和 P2,以使 P1-P0 和 P3-P2 的方向与其相关切线相匹配。 (P1-P0和P3-P2的距离,这个要求并没有固定,粗略的猜测是端点之间直线距离的三分之一。)

在下面的代码中,您可以看到函数 smooth() 执行此操作,并从支撑点(又名“结”)设置所有 P0、P1、P2、P3。

我已经对您的支撑点进行了硬编码,并将它们与复合贝塞尔代码一起绘制。我稍微转动了绘图,尝试从与您的 CAD 相同的视图中查看它;似乎有一个额外的小垂直环 - 您确定支撑点完全相同吗?最后,您的 CAD 没有使用复合三次贝塞尔曲线 - 它似乎使用您所说的支撑点作为高阶贝塞尔曲线的实际控制点。

请注意,我创建了一个简单的贝塞尔曲线类,以简化代码。

import math
import numpy as np
import matplotlib.pyplot as plt

############################################################

class Bezier:
    '''class for Cubic Bezier curve'''
    def __init__( self, q0, q1, q2, q3 ):
        self.p0 = np.array( q0 )
        self.p1 = np.array( q1 )
        self.p2 = np.array( q2 )
        self.p3 = np.array( q3 )

    def pt( self, t ):
        return ( 1.0 - t ) ** 3 * self.p0 + 3.0 * t * ( 1.0 - t ) ** 2 * self.p1 + 3.0 * t ** 2 * ( 1.0 - t ) * self.p2 + t ** 3 * self.p3

    def curve( self, n ):
        crv = []
        for t in np.linspace( 0.0, 1.0, n ):
            crv.append( self.pt( t ) )
        return crv

############################################################

def draw3d( BezierData ):
    x = [];   y = [];   z = []
    for B in BezierData:
        arc = Bezier( B[0], B[1], B[2], B[3] )
        for p in arc.curve( 100 ):
            x.append( p[0] )
            y.append( p[1] )
            z.append( p[2] )
    plt.plot( x, y, z )
    
############################################################

def length( vector ): return math.sqrt( np.dot( vector, vector ) )   # find length of a vector (numpy array)

def normalise( vector ): return vector / length( vector )            # normalise a vector (numpy array)

############################################################

def smooth( xyvals ):
    '''Turns list of coordinates into list of Bezier control points by adding intermediate control points'''
    N = len( xyvals )
    result = []

    for i in range( 0, N ):                                          # use tangents to set points P1 and P2
        im, ip, ipp = i - 1, ( i + 1 ) % N, ( i + 2 ) % N            # indices of adjacent points (cyclic)
        if im < 0: im += N
        pm = np.array( xyvals[im] )                                  # earlier point
        p0 = np.array( xyvals[i  ] )                                 # start point of Bezier curve
        p3 = np.array( xyvals[ip ] )                                 # end point of Bezier curve
        p6 = np.array( xyvals[ipp] )                                 # later point
        p1 = p0 + normalise( p3 - pm ) * length( p3 - p0 ) / 3.0     # use tangent at P0 to get P1
        p2 = p3 - normalise( p6 - p0 ) * length( p3 - p0 ) / 3.0     # use tangent at P3 to get P2
        result.append( [ p0, p1, p2, p3 ] )

    return result
    
############################################################

points = [ [0, 0, 0             ], [0,-20.836,0         ], [0,-38.948,0         ], [14.277,-38.948,0    ],
           [30.064,-38.948,0    ], [47.36,-38.948,0     ], [63.942,-38.948,3    ], [78.327,-38.948,6    ],
           [78.327,-30.509,9    ], [78.327,-22.191,12   ], [64.542,-22.191,15   ], [47.662,-22.191,18   ],
           [25.768,-22.191,21   ], [25.768,0.98532,24   ], [25.768,13.24,27.261 ], [25.768,23.53,30     ],
           [25.768,30.216,33    ], [42.9,30.216,33      ], [57.202,30.216,33    ], [75.909,30.216,33    ],
           [75.909,-0.99106,33  ], [75.909,-2.4834,40   ], [75.909,-2.6028,49   ], [76.909,3.2814,53    ],
           [77.909,8.2948,53    ], [78.909,14.468,49    ], [78.909,14.307,40    ], [78.909,11.782,33    ],
           [78.909,-15.121,36   ], [60.112,-15.121,39   ], [43.57,-15.121,42    ], [27.405,-15.121,48   ],
           [15.597,-15.121,48   ], [15.597,-3.0904,48   ], [15.597,11.195,48    ], [15.597,27.638,48    ],
           [26.208,27.638,48    ], [34.6,27.638,48      ], [34.6,14.052,48      ], [34.6,3.8624,48      ],
           [34.6,-5.5279,45     ], [44.19,-5.5279,42    ], [57.776,-5.5279,39   ], [67.565,-5.5279,33   ],
           [67.565,2.9199,27    ], [67.565,11.511,24    ], [67.565,19.902,21    ], [67.565,27.894,18    ],
           [67.565,39.226,15    ], [54.779,39.226,12    ], [43.924,39.226,9     ], [30.976,39.226,6     ],
           [18.356,39.226,3     ], [0,39.226,0          ], [0,17.662,0          ], [0,0,0               ] ]

plt.figure().add_subplot( projection='3d' )

Beziers = smooth( points )                       # create list of Bezier curves from support points
draw3d( Beziers )                                # draw composite Bezier curve

x = [ pt[0] for pt in points ]
y = [ pt[1] for pt in points ]
z = [ pt[2] for pt in points ]
plt.plot( x, y, z, 'o' )                         # plot original support points

plt.show()

输出(稍微旋转以匹配 CAD):

© www.soinside.com 2019 - 2024. All rights reserved.