在我的作业中,我尝试用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
正如我之前在评论中所写,对于平滑的复合贝塞尔曲线,您需要确保每个贝塞尔曲线末尾的梯度是连续的。本质上,您需要在此处估计切线方向,并选择 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()