I am trying to implement a research paper where part of my task is to project tooth models onto a curve surface to which I am currently stuck on.

The way it should work is by projecting each vertice on the tooth model to the surface. Please check the image below:

The full paper is provided here (https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=6070973).

As you can see, the teeth are presented on the image with the projection surface Pr. Imagine that this is placed on a 3D space, the curve should be a surface to which all the teeth will be projected on.

So far, this is my code:

import trimesh
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import bisect
from mpl_toolkits.mplot3d import Axes3D
import cv2

# Load the STL file
mesh = trimesh.load('teeth.stl')

# Polynomial coefficients for Pr and Pc
coeffs_Pr = [-6.12e-6, 0, -1.46e-2, 0, 0]  # Updated to reflect a curve
coeffs_Pc = [-1.08e-4, 0, -2.23e-2, 0, -7.07]

# Polynomial surface function
def polynomial_surface(coeffs, x):
    # print(sum(c * x**i for i, c in enumerate(coeffs)))
    return sum(c * x**i for i, c in enumerate(coeffs))

# Angle curve phi based on the graph provided (piecewise linear approximation)
def angle_curve_phi(distance):
    if distance <= 20:
        return 90 + (distance / 20) * 10
    elif distance <= 40:
        return 100 + ((distance - 20) / 20) * 10
    elif distance <= 60:
        return 110 - ((distance - 40) / 20) * 10
    elif distance <= 80:
        return 100 - ((distance - 60) / 20) * 20
    else:
        return 80 - ((distance - 80) / 20) * 30


# Root finding function
def find_root(f, a, b, max_iter=100):
    for _ in range(max_iter):
        try:
            return bisect(f, a, b)
        except ValueError:
            a -= 10
            b += 10
    raise ValueError("Root not found within the allowed number of iterations")

# Project vertex to surface function
def project_vertex_to_surface(vertex, coeffs_Pr, coeffs_Pc):
    xc = vertex[0]
    zc = vertex[2]
    
    try:
        x_c_prime = find_root(lambda x: polynomial_surface(coeffs_Pc, x) - xc, -10, 10)
    except ValueError:
        return None
    
    distance_from_midline = np.linalg.norm(vertex[:2])
    angle_phi = np.deg2rad(angle_curve_phi(distance_from_midline))
    ray_direction = np.array([np.cos(angle_phi), 0, np.sin(angle_phi)])
    
    try:
        x_r = find_root(lambda x: polynomial_surface(coeffs_Pr, x) - (vertex + x * ray_direction)[2], -10000, 10000)
    except ValueError:
        return None
    
    projected_vertex = vertex + x_r * ray_direction
    return projected_vertex

# Align mesh to the correct orientation (if necessary)
rotation_matrix = trimesh.transformations.rotation_matrix(
    angle=np.radians(-90), direction=[1, 0, 0], point=mesh.centroid
)
mesh.apply_transform(rotation_matrix)

rotation_matrix = trimesh.transformations.rotation_matrix(
    angle=np.radians(-90), direction=[1, 0, 0], point=mesh.centroid
)
mesh.apply_transform(rotation_matrix)

# Calculate the centroid of the mesh
centroid = mesh.centroid

# Translate the mesh to center it along the X-axis
translation_vector = [-centroid[0], 0, 0]
translation_matrix = trimesh.transformations.translation_matrix(translation_vector)
mesh.apply_transform(translation_matrix)

# Project all vertices
projected_vertices = []
for vertex in mesh.vertices:
    projected = project_vertex_to_surface(vertex, coeffs_Pr, coeffs_Pc)
    if projected is not None:
        projected_vertices.append(projected)

# Convert to numpy array for plotting
projected_vertices = np.array(projected_vertices)

print(projected_vertices)

# Define the range for the projection surfaces
x_values = np.linspace(-100, 100, 400)
z_values_Pr = polynomial_surface(coeffs_Pr, x_values)
z_values_Pc = polynomial_surface(coeffs_Pc, x_values)
y_values = np.zeros_like(x_values)

# Create a plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# Plot the mesh vertices
ax.scatter(mesh.vertices[:, 0], mesh.vertices[:, 1], mesh.vertices[:, 2], color='blue', s=1, label='Teeth Mesh')

# Plot the projection surfaces on the XY plane
ax.plot(x_values, y_values, z_values_Pr, color='green', label='Projection Surface Pr')
ax.plot(x_values, y_values, z_values_Pc, color='orange', label='Central Surface Pc')

# Plot the projected vertices
if len(projected_vertices) > 0:
    ax.scatter(projected_vertices[:, 0], projected_vertices[:, 1], projected_vertices[:, 2], color='red', s=1, label='Projected Vertices')

# Set labels and legend
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.legend()

# Show the plot
plt.show()

The projected_vertices should be an array to which it outputs an empty array indicating that no vertices are projected on the surface.

Any help will be appreciated!

For the teeth, you may use this model (https://cults3d.com/en/3d-model/art/full-anatomy-upper-and-lower-teeth-1-husseinhuzam3) for testing.