How to handle the situation where the result generated by buildCylinderBetweenTwoHoles collides between two holes?

[alpinebuster]

Bonjour 👋

The situation where the result generated by buildCylinderBetweenTwoHoles collides between two holes, as shown in the image below, which represents the cross-sectional view of the maxilla and mandible teeth, where there is an intersection during generation.

[Fedr]

Hello,

buildCylinderBetweenTwoHoles behavior depends a lot on the metric passed via StitchHolesParams.metric parameter. Have you tried different options?

If nothing helps, could you please attach your mesh for us?

[alpinebuster]

I tried several FillHoleMetric options but didn't get the expected results (FYI: the code I am using is generateOrthodonticBiteImpl and generateOrthodonticBiteWithFillHoleMetricImpl). I feel I can handle this using lower-level interfaces, but it's a bit complicated and time-consuming. According to the intended use case of buildCylinderBetweenTwoHoles, it should handle this situation well 🤔.

Here is the Mesh files used: bite-binary.zip.

[Grantim]

Hello!
Could you please demonstrate what result do you want to achieve on these meshes?

[alpinebuster]

First, invert the occlusal-facing surfaces of the maxillary and mandibular teeth, and then create a water-tight mesh (an orthodontic occlusal retainer) by connecting the peripheral margins of the cavities that join the upper and lower jaws. In the example image, some portions of the generated connections between the upper and lower margins intersect and penetrate the tooth structures themselves.

The desired final result is to form a water-tight mesh, but only the yellow areas in the image below.

[Grantim]

I see. It seems that in this case it will be hard to get nice result without providing additional guidelines. Please note that both fillHole and buildCylinderBetweenTwoHoles only take boundary into account, and does not deepend on inner surface of the mesh. It is also true for positionVertsSmoothly functions.
I think in this case it will be easier to add guildelines in the middle of the mesh for better result (or even several layers) like this:

from meshlib import mrmeshpy as mm
import math

meshA = mm.loadMesh(__file__ + "/../mandible-binary.stl")
meshB = mm.loadMesh(__file__ + "/../maxilla-binary.stl")

# combine meshes
mesh = mm.Mesh()
mesh.addMesh(meshA)
mesh.addMesh(meshB)

# flip normals
mesh.topology.flipOrientation()

avgEdgeLength = mesh.averageEdgeLength() # store it for future subdivision

# prepare guidelines (circular for simplicity)
holes = mesh.topology.findHoleRepresentiveEdges()
assert holes.size() == 2
box = mesh.computeBoundingBox()
center = box.center()
radius = box.diagonal()*0.5
# Z axis is known to be occlusal?
numGuides = 2 # play with range to get better result
for i in range(numGuides): 
    zRatio = float(i+1)/float(numGuides+1)
    zLevel = box.min.z * (1.0-zRatio) + box.max.z * zRatio
    points = mm.Contour3f()
    points.resize(200+1)
    for j in range(points.size()-1):
        angle = math.pi*2*j/(points.size()-1)
        points[j].x = center.x + radius * math.cos(angle)
        points[j].y = center.y + radius * math.sin(angle)
        points[j].z = zLevel
    points[points.size()-1] = points[0] # close points
    loopId = mesh.addSeparateEdgeLoop(points)
    holes.insert(holes.size() - 1, loopId.sym())
    holes.insert(holes.size() - 1, loopId)

# stitch holes
sParams = mm.StitchHolesParams()
sParams.metric = mm.getCircumscribedMetric(mesh)
oldFaces = mesh.topology.getValidFaces()
for i in range(int(holes.size()/2)):
    mm.buildCylinderBetweenTwoHoles(mesh,holes[i*2],holes[i*2+1],sParams)

# post process
newFaces = mesh.topology.getValidFaces() - oldFaces
sSettings = mm.SubdivideSettings()
sSettings.maxEdgeLen = avgEdgeLength * 5
sSettings.maxEdgeSplits = 1000000 # big number not to stop early 
sSettings.region = newFaces
mm.subdivideMesh(mesh,sSettings)

vertRegion = mm.getIncidentVerts(mesh.topology,newFaces)
mm.shrink(mesh.topology,vertRegion,3)

mm.positionVertsSmoothly(mesh,vertRegion,mm.EdgeWeights.Cotan,mm.VertexMass.NeiArea)

mm.saveMesh(mesh,__file__+"/../result.stl")

it will prodice result like this:

So, better guidelines will produce betteer result, but getting them might be tricky

[alpinebuster]

Got it, thks~

[alpinebuster]

Hello, I have a feedback. I have achieved the desired result based on the guidelines method~

from meshlib import mrmeshpy as mm
import math
import os
import pdb
import numpy as np


def extract_hole_boundary_points(mesh, hole_edges):
    """
    Extract the 3D coordinates of vertices around a hole boundary
    given a list of edges (hole_edges). Returns the vertex positions
    corresponding to the origin of each edge.
    """
    boundary_points = mm.Contour3f()

    for edge in hole_edges:
        vertex_id = mesh.topology.org(edge)
        vertex_pos = mesh.points[vertex_id]
        boundary_points.append(vertex_pos)

    return boundary_points

def compute_hole_centroid_and_normal(boundary_points):
    """
    Calculate the centroid (center point) and approximate normal direction
    of a hole boundary. This helps better understanding the hole's orientation.
    """
    n = boundary_points.size()
    if n < 3:
        return mm.Vector3f(0, 0, 0), mm.Vector3f(0, 0, 1)
    
    # Calculate centroid by averaging all boundary points
    # 
    # centroid = mm.Vector3f(0, 0, 0)
    # for i in range(n):
    #     centroid.x += boundary_points[i].x
    #     centroid.y += boundary_points[i].y
    #     centroid.z += boundary_points[i].z
    
    # centroid.x /= n
    # centroid.y /= n
    # centroid.z /= n
    # 
    pts = np.array([[p.x, p.y, p.z] for p in boundary_points])
    centroid_np = pts.mean(axis=0)
    centroid = mm.Vector3f(*centroid_np)

    
    centroid_np = pts.mean(axis=0)
    centroid = mm.Vector3f(*centroid_np)



    ###
    # PCA for the points
    # pdb.set_trace()
    pts_centered = pts - centroid_np
    cov = np.dot(pts_centered.T, pts_centered) / n
    eigvals, eigvecs = np.linalg.eigh(cov)
    normal_np = eigvecs[:, 0]  # the normal direction corresponding to the smallest eigenvalue

    normal = mm.Vector3f(*normal_np).normalized()
    ###

    
    return centroid, normal


def scale_boundary_to_bbox(original_points, original_centroid, target_centroid, bbox, custom_scale=1.):
    """Scale the original_points uniformly so they at least touch one face of the BoundingBox.
        margin: optional distance to scale slightly beyond the BoundingBox.
    """
    # Compute the vector from the original centroid to the target centroid
    translation_vector = mm.Vector3f(
        target_centroid.x - original_centroid.x,
        target_centroid.y - original_centroid.y,
        target_centroid.z - original_centroid.z
    )
    # Translate all points so centroid is exactly target_centroid
    translated_points = mm.Contour3f()
    for i in range(original_points.size()):
        translated_points.append(mm.Vector3f(
            original_points[i].x + translation_vector.x,
            original_points[i].y + translation_vector.y,
            # NOTE: this makes things complicated
            # original_points[i].z + translation_vector.z
            target_centroid.z
        ))


    # Find required scale_factor
    min_scale = float("inf")
    for p in translated_points:
        # Direction from new centroid
        dir_vec = mm.Vector3f(
            p.x - target_centroid.x,
            p.y - target_centroid.y,
            p.z - target_centroid.z
        )

        if dir_vec.x == dir_vec.y == dir_vec.z == 0:
            continue

        possible_scales = []
        for d, c, bmin, bmax in zip(
            (dir_vec.x, dir_vec.y, dir_vec.z),
            (target_centroid.x, target_centroid.y, target_centroid.z),
            (bbox.min.x, bbox.min.y, bbox.min.z),
            (bbox.max.x, bbox.max.y, bbox.max.z),
        ):
            if abs(d) < 1e-6:
                continue

            if d > 0:
                scale = (bmax - c) / d
            else:
                scale = (bmin - c) / d

            if scale > 0:
                possible_scales.append(scale)

        if possible_scales:
            min_scale = min(min_scale, min(possible_scales))

    scale_factor = min_scale * custom_scale
    # pdb.set_trace()


    # Scale up translated_points
    scaled_points = mm.Contour3f()
    scaled_points.resize(translated_points.size() + 1)
    for i in range(translated_points.size()):
        offset = mm.Vector3f(
            translated_points[i].x - target_centroid.x, 
            translated_points[i].y - target_centroid.y, 
            translated_points[i].z - target_centroid.z,
        )
        scaled_points[i] = mm.Vector3f(
            target_centroid.x + offset.x * scale_factor,
            target_centroid.y + offset.y * scale_factor,
            target_centroid.z
        )

    scaled_points[translated_points.size()] = scaled_points[0]  # close points

    return scaled_points



# Main
base_dir = os.path.dirname(__file__)
meshA_path = os.path.abspath(os.path.join(base_dir, "..", "_assets", "mandible-binary.stl"))
meshB_path = os.path.abspath(os.path.join(base_dir, "..", "_assets", "maxilla-binary.stl"))

# Load and combine the meshes
meshA = mm.loadMesh(meshA_path)
meshB = mm.loadMesh(meshB_path)

mesh = mm.Mesh()
mesh.addMesh(meshA)
mesh.addMesh(meshB)
# Flip normals to ensure consistent orientation
mesh.topology.flipOrientation()

# Store average edge length for later use in subdivision
avgEdgeLength = mesh.averageEdgeLength()


###
# Find the hole boundaries - these represent the openings we need to connect
holes = mesh.topology.findHoleRepresentiveEdges()
assert holes.size() == 2, f"Expected 2 holes, found {holes.size()}"

print(f"Found {holes.size()} holes to connect")
###


###
# Extract boundary information from both holes
holes_with_all_edges = mm.findRightBoundary(mesh.topology)
hole1_boundary = extract_hole_boundary_points(mesh, holes_with_all_edges[0])
hole2_boundary = extract_hole_boundary_points(mesh, holes_with_all_edges[1])

print(f"Hole 1 has {hole1_boundary.size()} boundary points")
print(f"Hole 2 has {hole2_boundary.size()} boundary points")

# Calculate centroids and normals for both holes
centroid1, normal1 = compute_hole_centroid_and_normal(hole1_boundary)
centroid2, normal2 = compute_hole_centroid_and_normal(hole2_boundary)

print(f"Hole 1 centroid: ({centroid1.x:.2f}, {centroid1.y:.2f}, {centroid1.z:.2f}), normal: {normal1}")
print(f"Hole 2 centroid: ({centroid2.x:.2f}, {centroid2.y:.2f}, {centroid2.z:.2f}), normal: {normal2}")
###


###
# Determine number of guide curves - more curves give smoother transitions
# but require more conservative scaling to avoid intersections
numGuides = 5
scale = 1.7
# Create intermediate guide curves by interpolating between the two holes
guide_holes = mm.std_vector_Id_EdgeTag()
bbox = mesh.computeBoundingBox()
for i in range(numGuides):
    # Calculate interpolation ratio (0 = closer to hole1, 1 = closer to hole2)
    zRatio = float(i + 1) / float(numGuides + 1)
    print(f"Creating guide curve {i+1}/{numGuides} at interpolation ratio {zRatio:.2f}")
    
    # Interpolate position between the two centroids
    intermediate_centroid = mm.Vector3f(
        centroid1.x * (1.0 - zRatio) + centroid2.x * zRatio,
        centroid1.y * (1.0 - zRatio) + centroid2.y * zRatio,
        centroid1.z * (1.0 - zRatio) + centroid2.z * zRatio
    )
    
    # Choose which boundary to use as template (alternate or use the larger one)
    if zRatio < 1/2:
        # Use hole1 boundary as template
        template_boundary = hole1_boundary
        template_centroid = centroid1
    else:
        # Use hole2 boundary as template
        template_boundary = hole2_boundary
        template_centroid = centroid2
    
    # Generate dynamically scaled and positioned guide curve
    guide_points = scale_boundary_to_bbox(
        template_boundary,
        template_centroid,
        intermediate_centroid,
        bbox,
        scale,
    )
    
    # Add this guide curve to the mesh as a separate edge loop
    loopId = mesh.addSeparateEdgeLoop(guide_points)
    if zRatio < 1/2:
        holes.insert(holes.size() - 1, loopId.sym())
        holes.insert(holes.size() - 1, loopId)
    else:
        holes.insert(holes.size() - 1, loopId)
        holes.insert(holes.size() - 1, loopId.sym())

print(f"Total curves to stitch (original two + all guidelines): {holes.size()}")
###


###
# Stitch all the curves together to create the connecting surface
sParams = mm.StitchHolesParams()
sParams.metric = mm.getCircumscribedMetric(mesh)
# Remember which faces existed before stitching
oldFaces = mesh.topology.getValidFaces()

# Stitch adjacent pairs of curves
for i in range(int(holes.size() / 2)):
    mm.buildCylinderBetweenTwoHoles(mesh,holes[i*2],holes[i*2+1],sParams)
###


# Post-processing: improve the quality of the newly created surface
print("Post-processing the new surface...")
# Identify the newly created faces
newFaces = mesh.topology.getValidFaces() - oldFaces

###
# Subdivide the new faces to improve surface quality
sSettings = mm.SubdivideSettings()
sSettings.maxEdgeLen = avgEdgeLength * 3  # Smaller multiplier for finer detail
sSettings.maxEdgeSplits = 1000000
sSettings.region = newFaces
mm.subdivideMesh(mesh, sSettings)
###


###
# Get vertices in the new region for smoothing
vertRegion = mm.getIncidentVerts(mesh.topology, newFaces)
# Shrink the region slightly to avoid affecting the original boundaries
mm.shrink(mesh.topology, vertRegion, 2)
# Smooth the vertices for a more natural surface
mm.positionVertsSmoothly(mesh, vertRegion, mm.EdgeWeights.Cotan, mm.VertexMass.NeiArea)
###


# Save the result
result_path = os.path.abspath(os.path.join(base_dir, "MeshFillHoleWithGuidelines.stl"))
mm.saveMesh(mesh, result_path)
print(f"Successfully created mesh connection and saved to: {result_path}")

[alpinebuster]

I have an ancillary question (may be off topic). I’ve exported part of MeshLib’s API to WebAssembly using Emscripten’s embind and open‑sourced the work as meshsdk (thanks — MeshLib performs excellently under WASM and the API is as convenient to use as the Python bindings). However, all the WASM binding code was handwritten and I want to automate the binding generation. I noticed MeshLib’s Python and C bindings were produced with mrbind. Would it be feasible to leverage mrbind to parse MeshLib and then emit embind bindings? Which aspects, steps, keywords, or design considerations should I focus on to implement automated WASM bindings? Any guidance would be greatly appreciated 🆘.

[adalisk-emikhaylov]

It should be possible, and we plan to do it eventually, but currently we plan to automate C# bindings first, and only then WASM.

It takes several man-months for each new language, so probably not feasible to do yourself. I haven't looked into embind too closely yet, not sure if we're going to use it or something else.