how to draw closed cubic Bezier curve with multiple segments (parametrized)

Question

In my assignment I tried to create a closed cubic Bezier curve for a roller coaster in python. I use cubic Bezier Formula on my code as mentioned here: https://en.wikipedia.org/wiki/B%C3%A9zier_curve#:~:text=Cubic%20B%C3%A9zier%20curve%20with%20four,that%20can%20be%20scaled%20indefinitely. I have support points (P0, P3) but not control points (P1, P2) and know that control points should get with linear equations and Gaussian Algorithm, but because can't find a way to implant them as a code, so manipulate the p0, p3 in any segments to make some control points. I know that it's wrong but did that to at least check the other parts of the code. now have a curve that is not similar to my curve that made by exact points implant in a cad application (rhino3d). I want to its because of control points implementation or my visualization has some mistakes. cheers.

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

Answer 1

As I wrote earlier in the comments, for a smooth composite Bezier curve you need to ensure that the gradient at the end of each Bezier is continuous.

If - as in your post - you intend to use the given "support points" as the end points of each cubic Bezier (P0 and P3) then you need to make an estimate of the tangent direction at those support points here and choose P1 and P2 to match the directions of P1-P0 and P3-P2 to their relevant tangents. (The distances, i.e. lengths of P1-P0 and P3-P2, are not fixed by this requirement. A rough guess is to take them as one third of the straight line distance between end points.) Note: there is now what I think is a better approach at the end of this post.

In the code below you can see function smooth() doing that, and setting all of P0, P1, P2, P3 from the support points (aka "knots").

I have hard-coded your support points and also plotted them together with the composite Bezier code. I have turned the plot a little to try and see it from the same view as your CAD; there appears to be one small extra vertical loop - are you sure the support points are exactly the same? Finally, your CAD is not using a composite cubic Bezier curve - it appears to be using what you are calling support points as the actual control points of a higher-order Bezier curve.

Note that I have created a simple Bezier class, to simplify the code.

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

Output (slightly rotated to match CAD):

EDIT:

A better approach (and one that looks more like your rhino3d CAD) is to treat the given points as the P1 and P2 points of the Bezier curves. P3 points are then chosen by averaging the P2 point of one segment with the P1 point of the next (and similarly for the P0 points). This ensures common tangency at the joined cubic Beziers. Note that, in this circumstance the curve doesn't necessarily go through any of the given points - they are all Bezier control points.

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

This time the roller-coaster is a little smoother: