I am trying to create a mesh (.stl file) that follows a bezier curve. I have succesfully created the side face but I have trouble creating the top and bottom faces. The problem is that the bezier curve is quite complex (the shape can be concave or convex) and the faces i have managed to create don’t follow the curve exactly (I will add a picture). Is there a way to make the top and bottom faces follow the curve exactly and also fill the shape (i don’t want it to be hollow).

This is the code I use to create the mesh.

import numpy as np
import trimesh
import shape
from scipy.spatial import Delaunay

def create_mesh(x, y, z_bottom, z_top):
    """ Create a mesh from the bezier curve points """
    # Create vertices for the bottom and top surfaces
    vertices = np.vstack((np.column_stack((x, y, z_bottom)), np.column_stack((x, y, z_top))))
    num_points = len(x)

    # Create faces
    faces = []
    
    for i in range(num_points - 1):
        # Side face
        faces.append([i, i + 1, num_points + i + 1])
        faces.append([i, num_points + i + 1, num_points + i])

    # Closing the side face
    faces.append([num_points - 1, 0, num_points])
    faces.append([num_points - 1, num_points, num_points + (num_points - 1)])

    # Top face
    top_vertices = np.column_stack((x, y, z_top))
    top = Delaunay(top_vertices[:, :2])
    top_faces = top.simplices + num_points

    # Bottom face
    bottom_vertices = np.column_stack((x, y, z_bottom))
    bottom = Delaunay(bottom_vertices[:, :2])
    bottom_faces = bottom.simplices
    
    # Combine faces
    faces = np.array(faces + top_faces.tolist() + bottom_faces.tolist())
    
    # Create the mesh
    mesh = trimesh.Trimesh(vertices=vertices, faces=faces)

    return mesh

And this is the shape.py code i use to get the curve.

import numpy as np
from scipy.special import binom
from scipy.interpolate import UnivariateSpline
import matplotlib.pyplot as plt

bernstein = lambda n, k, t: binom(n,k)* t**k * (1.-t)**(n-k)

def bezier(points, num=200):
    N = len(points)
    t = np.linspace(0, 1, num=num)
    curve = np.zeros((num, 2))
    for i in range(N):
        curve += np.outer(bernstein(N - 1, i, t), points[i])
    return curve

class Segment():
    def __init__(self, p1, p2, angle1, angle2, **kw):
        self.p1 = p1; self.p2 = p2
        self.angle1 = angle1; self.angle2 = angle2
        self.numpoints = kw.get("numpoints", 100)
        r = kw.get("r", 0.3)
        d = np.sqrt(np.sum((self.p2-self.p1)**2))
        self.r = r*d
        self.p = np.zeros((4,2))
        self.p[0,:] = self.p1[:]
        self.p[3,:] = self.p2[:]
        self.calc_intermediate_points(self.r)

    def calc_intermediate_points(self,r):
        self.p[1,:] = self.p1 + np.array([self.r*np.cos(self.angle1),
                                    self.r*np.sin(self.angle1)])
        self.p[2,:] = self.p2 + np.array([self.r*np.cos(self.angle2+np.pi),
                                    self.r*np.sin(self.angle2+np.pi)])
        self.curve = bezier(self.p,self.numpoints)


def get_curve(points, **kw):
    segments = []
    for i in range(len(points)-1):
        seg = Segment(points[i,:2], points[i+1,:2], points[i,2],points[i+1,2],**kw)
        segments.append(seg)
    curve = np.concatenate([s.curve[:-1] for s in segments])
    # curve = curve[::-1]
    return segments, curve

def ccw_sort(p):
    d = p-np.mean(p,axis=0)
    s = np.arctan2(d[:,0], d[:,1])
    return p[np.argsort(s),:]

def get_bezier_curve(rad=0.2, edgy=0, height=1.0):
    """ Given an array of points *a*, create a 3D curve through those points.
    *rad* controls the distance of control points.
    *edgy* controls how "edgy" the curve is.
    *height* is the fixed height to extend in the z direction."""
    
    a = get_random_points(n=7, scale=1)
    p = np.arctan(edgy) / np.pi + 0.5
    a = ccw_sort(a)
    a = np.append(a, np.atleast_2d(a[0, :]), axis=0)
    d = np.diff(a, axis=0)
    ang = np.arctan2(d[:, 1], d[:, 0])
    f = lambda ang: (ang >= 0) * ang + (ang < 0) * (ang + 2 * np.pi)
    ang = f(ang)
    ang1 = ang
    ang2 = np.roll(ang, 1)
    ang = p * ang1 + (1 - p) * ang2 + (np.abs(ang2 - ang1) > np.pi) * np.pi
    ang = np.append(ang, [ang[0]])
    a = np.append(a, np.atleast_2d(ang).T, axis=1)

    # Generate 2D curve
    s, c = get_curve(a, r=rad, method="var")
    x, y = c.T

    # Create 3D vertices by extending in the z direction
    z = np.zeros_like(x)
    z_top = np.full_like(x, height)

    # Combine into a 3D array
    return x, y, z, z_top, a


def get_random_points(n=5, scale=0.8, mindst=None, rec=0):
    """ create n random points in the unit square, which are *mindst*
    apart, then scale them."""
    mindst = mindst or .7/n
    a = np.random.rand(n,2)
    d = np.sqrt(np.sum(np.diff(ccw_sort(a), axis=0), axis=1)**2)
    if np.all(d >= mindst) or rec>=200:
        return a*scale
    else:
        return get_random_points(n=n, scale=scale, mindst=mindst, rec=rec+1)


def get_curvature(x, y=None, error=0.1):
    if y is None:
        x, y = x.real, x.imag

    t = np.arange(x.shape[0])
    std = error * np.ones_like(x)

    fx = UnivariateSpline(t, x, k=4, w=1 / np.sqrt(std))
    fy = UnivariateSpline(t, y, k=4, w=1 / np.sqrt(std))

    xˈ = fx.derivative(1)(t)
    xˈˈ = fx.derivative(2)(t)
    yˈ = fy.derivative(1)(t)
    yˈˈ = fy.derivative(2)(t)
    curvature = (xˈ* yˈˈ - yˈ* xˈˈ) / np.power(xˈ** 2 + yˈ** 2, 3 / 2)
    return curvature

This is an example of what i get (just adding the top and bottom faces and not filling the object). As you can see the top and bottom faces don’t match exactly.