正如我在前面的评论中所写的，对于一个平滑的复合贝塞尔曲线，你需要确保每个贝塞尔曲线末端的梯度是连续的。
如果--就像在你的帖子里一样--你打算使用给定的“支撑点”作为每个三次贝塞尔曲线的端点（P0和P3），那么您需要估计这些支撑点处的切线方向，并选择P1和P2，以将P1-P0和P3-P2的方向与它们的相关切线相匹配。（距离，即P1-P0和P3-P2的长度，不固定的这个要求.一个粗略的猜测是把它们作为端点之间的直线距离的三分之一.）**注意：**现在有什么我认为是一个更好的方法在这篇文章的结尾.
在下面的代码中，你可以看到函数smooth（）在做这件事，并从支撑点（又名“结”）设置所有的P0，P1，P2，P3。
我已经硬编码了您的支撑点，并将它们与复合Bezier代码一起绘制。我已经将图稍微旋转了一下，试图从与您的CAD相同的视图查看它;似乎有一个小的额外垂直环-您确定支撑点完全相同吗？最后，您的CAD并没有使用复合三次Bezier曲线-它似乎使用您所谓的支持点作为高阶Bezier曲线的实际控制点。
请注意，我创建了一个简单的Bezier类，以简化代码。

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）：x1c 0d1x

编辑：

更好的方法（一个看起来更像你的rhino 3d CAD）是把给定的点作为Bezier曲线的P1和P2点，然后通过将一个段的P2点与下一个段的P1点平均来选择P3点。（对于P0点也是类似的）。这确保了在连接的立方Bezier处的公共相切。注意，在这种情况下，曲线不一定要通过任何给定的点-他们都是贝塞尔控制点。

def smooth( xyzvals ):
    '''Turns list of coordinates (pairwise) into list of cubic Bezier control points'''
    N = len( xyzvals )
    result = []

    for i in range( 0, N, 2 ):                                       # set cubic Beziers s as P1 and P2
        p1 = np.array( xyzvals[i] )                                  # use original points as P1 and P2
        p2 = np.array( xyzvals[i+1] )

        im, ip = ( N + i - 1 ) % N, ( i + 2 ) % N                    # indices of adjacent points (cyclic)
        pm = np.array( xyzvals[im] )                                 
        p4 = np.array( xyzvals[ip] )
        p0 = 0.5 * ( pm + p1 )                                       # assign end points to ensure tangency
        p3 = 0.5 * ( p2 + p4 )

        result.append( [ p0, p1, p2, p3 ] )

    return result

这一次过山车稍微平稳一些：