multimoda-rs: CCTA Module Tutorial
This notebook demonstrates the full CCTA–intravascular fusion pipeline:
Read in and label a CCTA geometry from an STL file and centerline CSV files
Prepare centerlines, detect branches, and discretize the vessel tree
Load an intravascular geometry and fine-align it to the CCTA point cloud
Label the anomalous (intramural) region within the CCTA mesh
Compute radial scaling factors for the proximal, distal, and aortic regions
Morph the CCTA mesh to match the intravascular geometry
Remove the intramural region and stitch the CCTA to the intravascular geometry
Remesh and smooth the stitched geometry
Visualise labeled regions and export section STL files
Re-label the final stitched geometry
Note: Example data is loaded automatically from
examples/data/when run from the repository root.
Load required packages and set the working directory. Performance is typically faster in the console than in a notebook.
import os
from pathlib import Path
import numpy as np
import multimodars as mm
cwd = Path.cwd()
for candidate in [cwd, cwd.parent, cwd.parent.parent]:
if (candidate / "examples" / "data").exists():
os.chdir(candidate / "examples" / "data")
break
elif (candidate / "data").exists():
os.chdir(candidate / "data")
break
print(f"Working directory: {os.getcwd()}")
Working directory: /mnt/c/WorkingData/Documents/2_Coding/Rust/multimodars/examples/data
# Install if needed
%pip --disable-pip-version-check install trimesh plotly scipy pymeshlab | grep -v 'already satisfied'
import trimesh
import plotly.graph_objects as go
import plotly.io as pio
pio.renderers.default = "notebook_connected"
_DARK = dict(
scene=dict(aspectmode='data', bgcolor='black'),
margin=dict(r=0, l=0, b=0, t=30),
legend=dict(yanchor="top", y=0.99, xanchor="left", x=0.01),
paper_bgcolor="black",
)
def trimesh_to_mesh3d(mesh, color, name, opacity=0.6):
verts, faces = mesh.vertices, mesh.faces
return go.Mesh3d(
x=verts[:,0], y=verts[:,1], z=verts[:,2],
i=faces[:,0], j=faces[:,1], k=faces[:,2],
color=color, opacity=opacity, name=name, flatshading=True,
)
def _to_arr(pts):
return np.array(pts, dtype=np.float64) if pts else np.empty((0, 3))
def plot_labeled_geometry(results):
"""Visualise initial geometry labeling (aorta / RCA / LCA / removed)."""
mesh = results["mesh"]
fig = go.Figure()
fig.add_trace(go.Mesh3d(
x=mesh.vertices[:,0], y=mesh.vertices[:,1], z=mesh.vertices[:,2],
i=mesh.faces[:,0], j=mesh.faces[:,1], k=mesh.faces[:,2],
opacity=0.1, color='lightgray', name='Mesh', hoverinfo='skip',
))
for key, color, label in [
("aorta_points", "yellow", "Aorta"),
("rca_points", "blue", "RCA"),
("lca_points", "green", "LCA"),
("rca_removed_points", "red", "RCA removed (intramural)"),
("lca_removed_points", "red", "LCA removed (intramural)"),
]:
arr = _to_arr(results.get(key, []))
if len(arr):
fig.add_trace(go.Scatter3d(
x=arr[:,0], y=arr[:,1], z=arr[:,2],
mode='markers', marker=dict(size=2, color=color), name=label,
))
fig.update_layout(**_DARK)
fig.show()
def plot_anomalous_regions(results):
"""Visualise proximal / anomalous / distal sub-regions."""
mesh = results["mesh"]
fig = go.Figure()
fig.add_trace(go.Mesh3d(
x=mesh.vertices[:,0], y=mesh.vertices[:,1], z=mesh.vertices[:,2],
i=mesh.faces[:,0], j=mesh.faces[:,1], k=mesh.faces[:,2],
opacity=0.1, color='lightgray', name='Mesh', hoverinfo='skip',
))
for key, color, label in [
("proximal_points", "cyan", "Proximal"),
("anomalous_points", "orange", "Anomalous"),
("distal_points", "magenta", "Distal"),
]:
arr = _to_arr(results.get(key, []))
if len(arr):
fig.add_trace(go.Scatter3d(
x=arr[:,0], y=arr[:,1], z=arr[:,2],
mode='markers', marker=dict(size=2, color=color), name=label,
))
fig.update_layout(**_DARK)
fig.show()
def plot_results_key(results):
"""Plotly replacement for mm.plot_results_key — all region labels on one figure."""
mesh = results["mesh"]
fig = go.Figure()
fig.add_trace(go.Mesh3d(
x=mesh.vertices[:,0], y=mesh.vertices[:,1], z=mesh.vertices[:,2],
i=mesh.faces[:,0], j=mesh.faces[:,1], k=mesh.faces[:,2],
opacity=0.15, color='lightgray', name='Mesh', hoverinfo='skip',
))
for key, color, label in [
("aorta_points", "yellow", "Aorta"),
("rca_points", "blue", "RCA"),
("lca_points", "green", "LCA"),
("rca_removed_points", "red", "Intramural wall"),
("proximal_points", "cyan", "Proximal"),
("distal_points", "magenta", "Distal"),
("anomalous_points", "orange", "Anomalous"),
]:
arr = _to_arr(results.get(key, []))
if len(arr):
fig.add_trace(go.Scatter3d(
x=arr[:,0], y=arr[:,1], z=arr[:,2],
mode='markers', marker=dict(size=2, color=color), name=label,
))
fig.update_layout(**_DARK)
fig.show()
def plot_centerline_branches(rca_cl, lca_cl, results_dict=None):
"""Visualise centerline branch assignments for RCA and LCA.
Each branch gets its own toggleable legend entry. When *results_dict*
is provided, labelled surface-mesh points are overlaid so you can verify
that branch labels transferred correctly to the geometry.
Colour coding: blues = RCA (dark = main, lighter = side branches);
reds/oranges = LCA.
"""
_RCA = ["#1f77b4", "#17becf", "#9467bd", "#2ca02c", "#7f7f7f"]
_LCA = ["#d62728", "#ff7f0e", "#e377c2", "#bcbd22", "#8c564b"]
fig = go.Figure()
def add_cl(cl, colors, vessel):
from collections import defaultdict
by_branch = defaultdict(list)
for p in cl.points:
by_branch[p.branch_id].append(p.contour_point)
for bid in sorted(by_branch):
pts = by_branch[bid]
label = f"{vessel} main" if bid == 0 else f"{vessel} branch {bid}"
fig.add_trace(go.Scatter3d(
x=[p.x for p in pts], y=[p.y for p in pts], z=[p.z for p in pts],
mode='markers', marker=dict(color=colors[bid % len(colors)], size=3),
name=label,
))
def add_surf(key, color, label):
arr = _to_arr(results_dict.get(key, []))
if len(arr):
fig.add_trace(go.Scatter3d(
x=arr[:,0], y=arr[:,1], z=arr[:,2],
mode='markers', marker=dict(color=color, size=1.5, opacity=0.3),
name=label,
))
add_cl(rca_cl, _RCA, "RCA")
add_cl(lca_cl, _LCA, "LCA")
if results_dict is not None:
add_surf("rca_points_main", _RCA[0], "RCA main pts")
i = 1
while f"rca_points_side_{i}" in results_dict:
add_surf(f"rca_points_side_{i}", _RCA[i % len(_RCA)], f"RCA side {i} pts")
i += 1
add_surf("lca_points_main", _LCA[0], "LCA main pts")
i = 1
while f"lca_points_side_{i}" in results_dict:
add_surf(f"lca_points_side_{i}", _LCA[i % len(_LCA)], f"LCA side {i} pts")
i += 1
fig.update_layout(title="Centerline Branch Assignment", **_DARK)
fig.show()
def plot_vessel_tree(tree, pts_per_contour=24):
"""Visualise a discretized vessel tree.
Renders contour rings (subsampled to *pts_per_contour* points),
yellow centroid dots tracing each vessel axis, and orientation
reference triplets at the ostium and bifurcation sites.
Colour coding: silver = aorta; steel-blue = RCA main; shades of blue =
RCA side branches; coral = LCA main; shades of orange = LCA side branches.
Red × = main reference; orange ▲ = CCW reference; magenta ▼ = CW reference.
"""
_RCA = ["#1f77b4", "#4fa3e0", "#7ec8e3", "#a8d8ea"]
_LCA = ["#d62728", "#e07f4f", "#e3a87e", "#eac0a8"]
fig = go.Figure()
def ring_xyz(contour):
pts = contour.points
m = len(pts)
if m == 0:
return [], [], []
step = max(1, m // pts_per_contour)
s = pts[::step]
return ([p.x for p in s] + [s[0].x],
[p.y for p in s] + [s[0].y],
[p.z for p in s] + [s[0].z])
def add_contours(contours, color, name, opacity=0.7):
rx, ry, rz, cx, cy, cz = [], [], [], [], [], []
for c in contours:
xs, ys, zs = ring_xyz(c)
if not xs:
continue
rx += xs + [None]; ry += ys + [None]; rz += zs + [None]
if c.centroid:
cx.append(c.centroid[0]); cy.append(c.centroid[1]); cz.append(c.centroid[2])
if rx:
fig.add_trace(go.Scatter3d(x=rx, y=ry, z=rz, mode='lines',
name=name, line=dict(color=color, width=1), opacity=opacity))
if cx:
fig.add_trace(go.Scatter3d(x=cx, y=cy, z=cz, mode='markers',
marker=dict(color='yellow', size=2), showlegend=False,
name=f"{name} centroids"))
def add_refs(refs, label):
mx, my, mz, ccx, ccy, ccz, clx, cly, clz = [], [], [], [], [], [], [], [], []
for main_ref, cc_ref, clock_ref in refs:
mx.append(main_ref[0]); my.append(main_ref[1]); mz.append(main_ref[2])
ccx.append(cc_ref[0]); ccy.append(cc_ref[1]); ccz.append(cc_ref[2])
clx.append(clock_ref[0]); cly.append(clock_ref[1]); clz.append(clock_ref[2])
if mx:
fig.add_trace(go.Scatter3d(x=mx, y=my, z=mz, mode='markers',
name=f"{label} main ref", marker=dict(color='red', symbol='x', size=5)))
fig.add_trace(go.Scatter3d(x=ccx, y=ccy, z=ccz, mode='markers',
name=f"{label} CCW ref", marker=dict(color='orange', symbol='diamond', size=5)))
fig.add_trace(go.Scatter3d(x=clx, y=cly, z=clz, mode='markers',
name=f"{label} CW ref", marker=dict(color='magenta', symbol='square', size=5)))
add_contours(tree.discretized_aorta, "silver", "Aorta", opacity=0.3)
add_contours(tree.discretized_rca_main, "steelblue", "RCA main")
for i, branch in enumerate(tree.rca_branches):
add_contours(branch, _RCA[i % len(_RCA)], f"RCA branch {i+1}")
add_contours(tree.discretized_lca_main, "coral", "LCA main")
for i, branch in enumerate(tree.lca_branches):
add_contours(branch, _LCA[i % len(_LCA)], f"LCA branch {i+1}")
add_refs(tree.rca_references, "RCA")
add_refs(tree.lca_references, "LCA")
fig.update_layout(title="Discretized Vessel Tree", **_DARK)
fig.show()
def plot_centerline_edges(cl, cos_threshold=0.0, title="Centerline Edges"):
"""Show centerline branches with sharp-angle positions as red x markers.
Each branch gets a distinct colour. Large red x markers indicate positions
from find_sharp_angles — use them to decide where to call
split_branch or merge_branches.
"""
_PALETTE = [
"#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd",
"#8c564b", "#e377c2", "#7f7f7f", "#bcbd22", "#17becf",
]
from collections import defaultdict
by_branch = defaultdict(list)
for p in cl.points:
by_branch[p.branch_id].append(p.contour_point)
fig = go.Figure()
for bid in sorted(by_branch):
pts = by_branch[bid]
color = _PALETTE[bid % len(_PALETTE)]
label = "main" if bid == 0 else f"branch {bid}"
fig.add_trace(go.Scatter3d(
x=[p.x for p in pts], y=[p.y for p in pts], z=[p.z for p in pts],
mode='lines+markers', name=label,
line=dict(color=color, width=2),
marker=dict(color=color, size=3),
))
sharp_pos = cl.find_sharp_angles(bid, cos_threshold)
if sharp_pos:
sharp_pts = [pts[i] for i in sharp_pos if i < len(pts)]
fig.add_trace(go.Scatter3d(
x=[p.x for p in sharp_pts], y=[p.y for p in sharp_pts], z=[p.z for p in sharp_pts],
mode='markers', name=f"{label} sharp angles",
marker=dict(color='red', size=8, symbol='x'),
))
fig.update_layout(title=title, **_DARK)
fig.show()
Note: you may need to restart the kernel to use updated packages.
1. Read in and prepare CCTA geometries and centerlines
The entry point for the CCTA module is mm.label_geometry, which reads a triangulated surface mesh
(STL) together with centerline CSV files for the aorta, RCA, and LCA. It returns a labeled results
dictionary and three PyCenterline objects used in all subsequent steps.
The centerline CSV files must contain three columns (no header): x, y, z in mm.
Key parameters:
bounding_sphere_radius_mm: radius of the rolling sphere for the initial vessel labeling pass. Larger values capture more distant vertices; smaller values are more conservative.n_points_intramural: number of centerline points defining the end of the intramural segment.anomalous_rca/anomalous_lca: whenTrue, re-assigns incorrectly labeled intramural points. Set toFalsefor normal coronary anatomy.control_plot: open an interactive scene to inspect labeling. Useful when tuning the parameters above.
The returned results dictionary contains keys "mesh", "aorta_points", "rca_points",
"lca_points", "rca_removed_points", and "lca_removed_points".
rca_cl_raw = np.genfromtxt("./centerline_rca_short.csv", delimiter=',')
lca_cl_raw = np.genfromtxt("./centerline_lca.csv", delimiter=',')
aorta_cl_raw = np.genfromtxt("./centerline_aorta.csv", delimiter=',')
rca_cl = mm.numpy_to_centerline(rca_cl_raw)
lca_cl = mm.numpy_to_centerline(lca_cl_raw)
aorta_cl = mm.numpy_to_centerline(aorta_cl_raw)
results, (rca_cl, lca_cl, ao_cl) = mm.label_geometry(
path_ccta_geometry="./NARCO_119_noside.stl",
path_centerline_aorta="./centerline_aorta.csv",
path_centerline_rca="./centerline_rca_short.csv",
path_centerline_lca="./centerline_lca.csv",
bounding_sphere_radius_mm=3.0,
n_points_intramural=100,
anomalous_rca=True,
anomalous_lca=False,
control_plot=False,
)
plot_labeled_geometry(results)
Loaded mesh: 25171 vertices, 50338 faces
Loaded aorta centerline: 75 points
Loaded LCA centerline: 1092 points
Loaded RCA centerline: 788 points
RCA points found: 4820
LCA points found: 5899
Applying occlusion removal for anomalous RCA...
Total faces to exclude: 232
Excluded 232 faces, removed 110 points (filtered from 4820 to 4710 points)
RCA: relabeled 110 points in intramual course
Removing LCA and RCA island points...
length before: 5899
length after: 5879
Applying final reclassification based on adjacency map...
aorta_points:14472
rca_points:4710
lca_points:5879
rca_removed_points:110
lca_removed_points:0
2. Discretize CCTA and map reference points
prepare_centerlines detects branches on both coronary centerlines, validates topology, and labels
the surface-mesh points by branch — all in one call. The returned results dictionary gains keys
rca_points_main, rca_points_side_1, …, lca_points_main, lca_points_side_1, ….
After verifying the branch assignment in the interactive plot, discretize_vessel_tree slices each
vessel into evenly-spaced cross-sectional contours and computes orientation reference triplets at
the ostium and every side-branch bifurcation.
Key parameters:
branch_sigma: Gaussian smoothing radius (mm) for branch detection — increase if the algorithm over-segments noisy centerlines.step_size: arc-length spacing between consecutive cross-sections (mm).n_points: number of evenly-spaced points per output contour ring.b_spline: replace each discretized contour with a closed periodic B-spline fit — useful for smoothing jagged contours before computing reference points. Control strength withbspline_smoothing(0 = exact; ≈n_points= gentle; ≈5 × n_points= strong).
rca_cl, lca_cl, results = mm.prepare_centerlines(
rca_cl, lca_cl, results,
branch_sigma=2.0,
)
# Inspect each centerline for sharp angles — red × marks a position that may
# need a split_branch / merge_branches fix. Tune cos_threshold toward negative
# values (e.g. -0.3) to suppress false positives on gently curved segments.
plot_centerline_edges(lca_cl, cos_threshold=0.0, title="LCA edges")
# Manual corrections (uncomment and adapt as needed):
list_edges = lca_cl.find_sharp_angles(branch_id=0, cos_threshold=0.0)
lca_cl = lca_cl.split_branch(0, list_edges[4])
lca_cl = lca_cl.merge_branches(0, 4)
lca_cl = lca_cl.check_centerline()
results = mm.label_branches(lca_cl, results, results_key="lca_points")
# Verify final branch assignment after any corrections.
plot_centerline_branches(rca_cl, lca_cl, results)
Branch labeling for 'rca_points' (branch_ids=[0]):
rca_points_main: 3749
rca_points_side: 961
rca_points_side_1: 515
rca_points_side_2: 423
rca_points_side_3: 23
Branch labeling for 'lca_points' (branch_ids=[0]):
lca_points_main: 5228
lca_points_side: 651
lca_points_side_1: 248
lca_points_side_2: 277
lca_points_side_3: 114
lca_points_side_4: 12
Branch labeling for 'lca_points' (branch_ids=[0]):
lca_points_main: 2623
lca_points_side: 3256
lca_points_side_1: 496
lca_points_side_2: 277
lca_points_side_3: 114
lca_points_side_4: 2617
tree = mm.discretize_vessel_tree(
ao_cl, rca_cl, lca_cl,
results,
step_size=1.0,
n_points=100,
b_spline=True, # set True + tune bspline_smoothing to smooth noisy contours
bspline_smoothing=5.0,
)
plot_vessel_tree(tree)
3. Load and align intravascular geometry
Load the intravascular segmentation with mm.from_file_singlepair (see the intravascular tutorial
for the full range of loading options and parameter tuning). Then align it to the CCTA centerline
and point cloud with mm.align_combined, which first performs a coarse three-point alignment using
anatomical landmarks and then refines the rotation by minimising Hausdorff distances.
Landmark points (in mm, CCTA coordinate system):
Aortic reference — centre of vessel on the aortic side
Superior reference — proximal end of the intramural segment
Inferior reference — distal end of the intramural segment
Key parameters for align_combined:
angle_range_deg: angular search window (±degrees) for Hausdorff refinement.write/watertight/output_dir: export OBJ meshes whenwrite=True.
rest, (dia_logs, sys_logs) = mm.from_file_singlepair(
input_path="ivus_rest",
labels=["aligned_dia", "aligned_sys"],
output_path="output/rest",
)
ref_points = tree.rca_references[0]
rca_cl_main = rca_cl.get_branch(0) # alignment needs single-branch CL
aligned, resampled_cl = mm.align_combined(
rca_cl_main,
rest,
ref_points[0], # aortic reference point
ref_points[1], # superior reference point
ref_points[2], # inferior reference point
results['rca_points'], # CCTA point cloud for Hausdorff refinement
angle_range_deg=30.0,
write=True,
watertight=False,
output_dir="test",
align_wall_anomalous=True,
)
✅ Successfully built geometry from path
-----------------------------------------
✅ Lumen
❌ Eem
❌ Calcification
❌ Sidebranch
✅ Catheter
-----------------------------------------
Label: aligned_dia
Diastole phase: Yes
✅ Successfully built geometry from path
-----------------------------------------
✅ Lumen
❌ Eem
❌ Calcification
❌ Sidebranch
✅ Catheter
-----------------------------------------
Label: aligned_sys
Diastole phase: No
+--------------------------------------------------------------------+
| ✅ Finished aligning 'aligned_sys' (anomalous: true) |
+---------+------------+---------------+-------+-------+-------------+
| Contour | Matched To | Rotation (°) | Tx | Ty | Centroid |
+---------+------------+---------------+-------+-------+-------------+
| 1 | 0 | -2.50 | -0.19 | -0.01 | (3.63,5.12) |
| 2 | 1 | 12.00 | -0.54 | 0.15 | (3.63,5.12) |
| 3 | 2 | -10.00 | -0.44 | 0.53 | (3.63,5.12) |
| 4 | 3 | -4.00 | -0.72 | 0.94 | (3.63,5.12) |
| 5 | 4 | 1.00 | -0.75 | 0.26 | (3.63,5.12) |
| 6 | 5 | 5.00 | -0.79 | 0.65 | (3.63,5.12) |
| 7 | 6 | -6.50 | -1.21 | 0.69 | (3.63,5.12) |
| 8 | 7 | 15.00 | -1.08 | 0.81 | (3.63,5.12) |
| 9 | 8 | -35.00 | -1.03 | 0.24 | (3.63,5.12) |
| 10 | 9 | 25.50 | -1.46 | 0.28 | (3.63,5.12) |
| 11 | 10 | 40.00 | -2.40 | 0.73 | (3.63,5.12) |
| 12 | 11 | -21.00 | -1.95 | 0.68 | (3.63,5.12) |
| 13 | 12 | 5.00 | -2.06 | 0.66 | (3.63,5.12) |
| 14 | 13 | -7.00 | -2.14 | 0.52 | (3.63,5.12) |
| 15 | 14 | 23.00 | -2.27 | 0.46 | (3.63,5.12) |
| 16 | 15 | 7.50 | -2.28 | 0.53 | (3.63,5.12) |
+---------+------------+---------------+-------+-------+-------------+
+--------------------------------------------------------------------+
| ✅ Finished aligning 'aligned_dia' (anomalous: true) |
+---------+------------+---------------+-------+-------+-------------+
| Contour | Matched To | Rotation (°) | Tx | Ty | Centroid |
+---------+------------+---------------+-------+-------+-------------+
| 1 | 0 | 7.50 | 0.12 | -0.36 | (3.72,5.25) |
| 2 | 1 | 1.00 | 0.04 | -0.16 | (3.72,5.25) |
| 3 | 2 | 6.00 | -0.28 | 0.03 | (3.72,5.25) |
| 4 | 3 | -13.00 | -0.16 | 0.45 | (3.72,5.25) |
| 5 | 4 | -24.00 | -0.32 | 0.89 | (3.72,5.25) |
| 6 | 5 | -5.00 | -0.48 | 1.10 | (3.72,5.25) |
| 7 | 6 | -29.50 | -0.94 | 0.82 | (3.72,5.25) |
| 8 | 7 | -18.00 | -1.29 | 0.87 | (3.72,5.25) |
| 9 | 8 | 1.50 | -1.36 | 0.86 | (3.72,5.25) |
| 10 | 9 | -1.00 | -1.68 | 1.27 | (3.72,5.25) |
| 11 | 10 | 47.00 | -1.54 | 1.87 | (3.72,5.25) |
| 12 | 11 | 3.50 | -1.45 | 1.95 | (3.72,5.25) |
| 13 | 12 | 1.00 | -1.51 | 2.08 | (3.72,5.25) |
| 14 | 13 | 30.50 | -0.96 | 2.12 | (3.72,5.25) |
| 15 | 14 | -7.00 | -1.11 | 1.93 | (3.72,5.25) |
| 16 | 15 | 19.50 | -0.68 | 2.47 | (3.72,5.25) |
| 17 | 16 | 3.50 | -0.50 | 2.35 | (3.72,5.25) |
| 18 | 17 | -1.00 | -0.58 | 2.24 | (3.72,5.25) |
| 19 | 18 | -11.00 | -0.95 | 2.37 | (3.72,5.25) |
+---------+------------+---------------+-------+-------+-------------+
✅ Aligned geometry 'aligned_sys' to 'aligned_dia'
-----------------------------------------
Applied initial translation: (0.09, 0.13, -3.94) mm
Found best rotation of 3.50° with parameters:
range: 90.00°
step size: 0.5°
Applied final translation: ( 0, 0.00, 0.00) mm
-----------------------------------------
Saving files for 'aligned_dia - aligned_sys' to 'output/rest'
LUMEN .obj files: 2/2 written successfully
CATHETER .obj files: 2/2 written successfully
WALL .obj files: 2/2 written successfully
Step 1: Finding initial rotation via three-point method
---------------------Centerline alignment: Finding optimal rotation---------------------
✅ Best angle found: 219.00°
Step 2: Refining with Hausdorff distance
---------------------Refining alignment with Hausdorff---------------------
Initial rotation: 0.00°, Initial CL index: 22
Refined rotation: 5.00°, Refined CL index: 21, Hausdorff: 5.12
---------------------Applying final transformation---------------------
Total rotation (initial + delta): 224.00°
Moving ostium by 1 centerline points
Saving files for 'None' to 'test'
LUMEN .obj files: 2/2 written successfully
CATHETER .obj files: 2/2 written successfully
WALL .obj files: 2/2 written successfully
sidebranch file not found, skipping: "ivus_rest/branch_diastolic_contours.csv"
eem file not found, skipping: "ivus_rest/eem_diastolic_contours.csv"
process_directory: unknown mapping name 'catheter', skipping
calcification file not found, skipping: "ivus_rest/calcium_diastolic_contours.csv"
process_directory: unknown mapping name 'catheter', skipping
eem file not found, skipping: "ivus_rest/eem_systolic_contours.csv"
sidebranch file not found, skipping: "ivus_rest/branch_systolic_contours.csv"
calcification file not found, skipping: "ivus_rest/calcium_systolic_contours.csv"
resample_centerline_by_contours: centroid_count=14, centroid_mean_spacing=Some(1.2955758552631584), centerline_length=185.4959091541807, spacing=1.295576
resample_centerline_by_contours: produced 144 points
4. Label the anomalous region
mm.label_anomalous_region subdivides the RCA points into three sub-regions — proximal, anomalous
(intramural), and distal — based on spatial overlap between the aligned intravascular frames and the
CCTA mesh. The results dictionary is extended with "proximal_points", "anomalous_points",
and "distal_points".
Set debug_plot=True to open an interactive scene that shows the sub-region boundaries — useful
when a boundary appears misplaced.
results = mm.label_anomalous_region(
centerline=rca_cl,
frames=aligned.geom_a.frames,
results=results,
results_key='rca_points',
debug_plot=False,
)
plot_anomalous_regions(results)
Applying anomalous labeling based on aligned intravascular frames...
proximal_points: 36
distal_points: 4273
anomalous_points: 401
5. Compute scaling factors
Before morphing, the optimal radial scaling factor for each region is computed independently. Each function searches for the scale that minimises distance between the CCTA mesh and the corresponding portion of the aligned intravascular geometry. All return values are signed floats in mm (+ expands, − contracts).
Function |
Region |
|---|---|
|
proximal and distal segments |
|
aortic vertices, minimising distance to outer intravascular wall |
|
aortic wall vertices, targeting the intramural→free-segment transition |
Note:
find_aortic_wall_scalingraisesValueErrorif no frame with elliptic ratio < 1.3 is found (nearly circular geometry). Omit this step in that case.
prox_scaling, distal_scaling = mm.find_distal_and_proximal_scaling(
frames=aligned.geom_a.frames,
centerline=rca_cl,
results=results,
)
aortic_scaling = mm.find_aorta_scaling(
frames=aligned.geom_a.frames,
cl_aorta=ao_cl,
results=results,
)
aortic_wall_scaling = mm.find_aortic_wall_scaling(
frames=aligned.geom_a.frames,
cl_aorta=ao_cl,
results=results,
)
print(f"Proximal scaling: {prox_scaling:.3f} mm")
print(f"Distal scaling: {distal_scaling:.3f} mm")
print(f"Aortic scaling: {aortic_scaling:.3f} mm")
print(f"Aortic wall scaling: {aortic_wall_scaling:.3f} mm")
Finding best proximal/distal radial scaling factors...
Proximal scaling: -0.4 mm
Distal scaling: 0.4 mm
Finding best aortic radial scaling factor...
Aortic scaling: -0.4 mm
Finding best aortic wall radial scaling factor...
elliptic ratio <1.3 for frame index 7
Aortic wall scaling: 1.16 mm
Proximal scaling: -0.400 mm
Distal scaling: 0.400 mm
Aortic scaling: -0.400 mm
Aortic wall scaling: 1.158 mm
6. Morph CCTA to intravascular geometry
mm.scale_region_centerline_morphing moves mesh vertices radially along the local centerline normal
so that the region diameter changes by diameter_adjustment_mm. After each call,
mm.sync_results_to_mesh updates all coordinate lists in results to reflect the new vertex
positions. Always sync before the next morphing step.
The three regions are morphed in order: distal → aortic → proximal.
# 1. Scale the distal segment along the RCA centerline
scaled_distal = mm.scale_region_centerline_morphing(
mesh=results['mesh'],
region_points=results['distal_points'],
centerline=rca_cl,
diameter_adjustment_mm=distal_scaling,
)
results = mm.sync_results_to_mesh(results, results['mesh'], scaled_distal)
# 2. Scale the aortic region (aorta + intramural wall) along the aortic centerline
scaled_distal_aortic = mm.scale_region_centerline_morphing(
mesh=results['mesh'],
region_points=results['aorta_points'] + results['rca_removed_points'],
centerline=aorta_cl,
diameter_adjustment_mm=aortic_scaling,
)
results = mm.sync_results_to_mesh(results, results['mesh'], scaled_distal_aortic)
# 3. Scale the proximal segment along the RCA centerline
scaled_proximal = mm.scale_region_centerline_morphing(
mesh=results['mesh'],
region_points=results['proximal_points'],
centerline=rca_cl,
diameter_adjustment_mm=prox_scaling,
)
results = mm.sync_results_to_mesh(results, results['mesh'], scaled_proximal)
# Overlay scaled CCTA with the aligned intravascular geometry
anomaly_mesh = trimesh.load("test/lumen_000_None.obj")
anomaly_wall_mesh = trimesh.load("test/wall_000_None.obj")
fig = go.Figure(data=[
trimesh_to_mesh3d(results['mesh'], "darkred", "CCTA scaled", opacity=0.5),
trimesh_to_mesh3d(anomaly_mesh, "royalblue", "Intravascular lumen", opacity=0.7),
trimesh_to_mesh3d(anomaly_wall_mesh, "lightblue", "Intravascular wall", opacity=0.3),
])
fig.update_layout(scene=dict(aspectmode="data"), margin=dict(l=0, r=0, t=30, b=0))
fig.show()
Scaling 4273 vertices around Centerline(len=786, spacing=0.57 mm)
Diameter adjustment: 0.4 mm
Scaling 14582 vertices around Centerline(len=75, spacing=0.75 mm)
Diameter adjustment: -0.4 mm
Scaling 36 vertices around Centerline(len=786, spacing=0.57 mm)
Diameter adjustment: -0.4 mm
7. Remove intramural region and stitch geometries
mm.remove_labeled_points_from_mesh deletes the anomalous and proximal vertices, opening a boundary
ring at the proximal end of the intravascular segment. The results dictionary gains a
"boundary_points" key.
mm.stitch_ccta_to_intravascular then triangulates a patch connecting this boundary ring to the
proximal contour of the aligned intravascular geometry.
prox_start_mode:
"nearest_iv"(default) — rotate to the boundary vertex closest to intravascular point 0."highest_z"— rotate to the boundary vertex with the largest z-coordinate; prefer for straight intramural segments where the pullback axis is nearly aligned with the image z-axis.
Anomalous ostium clamping (clamp_overshoot)
For anomalous coronary arteries the intramural course runs roughly perpendicular to the aortic wall. This means the plane of the intravascular (IV) ring at the ostium is nearly at 90° to the plane of the aortic boundary ring. Without correction, some aortic boundary vertices end up on the wrong side of the ostium plane (i.e. inside the coronary lumen), creating a jagged seam.
When the angle between the two planes exceeds 45°, the following three-step correction is applied automatically:
Clamp — any boundary vertex that has crossed to the intravascular side of the ostium plane is projected back onto that plane.
Minimum gap (
clamp_overshoot, default 0.5 mm) — every boundary vertex (including freshly clamped ones) is pushed at leastclamp_overshootmm away from the ostium plane on the aortic side. This prevents the stitching patch from meeting the aortic wall at a razor-sharp angle at the ostium edge.Layer propagation — the two rings of aortic mesh vertices adjacent to the boundary ring are shifted radially outward from the coronary axis (within the ostium plane), by 0.1 mm and 0.2 mm respectively. This prevents those second-layer vertices from sitting closer to the coronary axis than the boundary ring itself, which would otherwise create a visible ridge just outside the seam.
clamp_overshoot can be tuned per case. Set it to 0 to disable steps 2 and 3 (pure clamping
only).
updated_results = mm.remove_labeled_points_from_mesh(
results,
["anomalous_points", "proximal_points"],
)
stitched = mm.stitch_ccta_to_intravascular(
aligned.geom_a,
updated_results['mesh'],
updated_results,
prox_start_mode="highest_z",
clamp_overshoot=0.5,
)
stitched['mesh'].export("prefixed_mesh.stl")
print("Raw stitched mesh exported → prefixed_mesh.stl")
Applying removal of '['anomalous_points', 'proximal_points']'
Removed 437
Created 56 boundary points
Stitching: 100/100 triangles created (n_boundary=37, n_iv=100, step=2, remainder=26)
Stitching: 100/100 triangles created (n_boundary=19, n_iv=100, step=5, remainder=5)
Raw stitched mesh exported → prefixed_mesh.stl
# Visualise the stitching seam:
# boundary ring (red→blue by index) + IV frame-0 lumen (red→blue by index)
boundary_pts = np.array(stitched['prox_boundary_points'], dtype=np.float64)
n_bnd = len(boundary_pts)
iv_viz = aligned.geom_a.downsample(100).sort_frame_points()
iv_pts = iv_viz.frames[0].lumen.points
iv_xyz = np.array([[p.x, p.y, p.z] for p in iv_pts])
n_iv = len(iv_xyz)
fig = go.Figure()
fig.add_trace(trimesh_to_mesh3d(stitched['mesh'], "lightgray", "Stitched mesh", opacity=0.2))
fig.add_trace(go.Scatter3d(
x=boundary_pts[:,0], y=boundary_pts[:,1], z=boundary_pts[:,2],
mode='markers',
marker=dict(size=4, color=list(range(n_bnd)), colorscale='RdBu', showscale=False),
name='CCTA boundary ring',
))
fig.add_trace(go.Scatter3d(
x=iv_xyz[:,0], y=iv_xyz[:,1], z=iv_xyz[:,2],
mode='markers',
marker=dict(size=5, color=list(range(n_iv)), colorscale='RdBu', showscale=False),
name='IV frame-0 lumen',
))
fig.update_layout(scene=dict(aspectmode='data'), margin=dict(l=0, r=0, t=30, b=0))
fig.show()
8. Remesh and smooth
mm.fix_and_remesh_stitched_mesh applies a three-step repair pipeline:
non-manifold repair → hole filling → isotropic remesh (via pymeshlab).
Taubin smoothing then reduces surface noise while preserving overall shape.
Requires:
pip install 'multimodars[meshlab]'
Key parameters:
target_edge_length_mm: desired edge length after remeshing (Noneuses the 25th-percentile edge length of the input mesh).verbose: print per-step vertex/face counts and watertightness.
remeshed = stitched.copy()
remeshed['mesh'] = mm.fix_and_remesh_stitched_mesh(
stitched['mesh'],
target_edge_length_mm=0.5,
verbose=True,
)
print(f"Watertight? {remeshed['mesh'].is_watertight}")
trimesh.smoothing.filter_taubin(remeshed['mesh'], lamb=0.6)
remeshed['mesh'].export("fixed_mesh.stl")
print("Remeshed and smoothed geometry exported → fixed_mesh.stl")
[input ] verts= 26,134 faces= 52,264 watertight=False
non-manifold edges/vertices repaired
holes closed
[after hole fill ] verts= 26,134 faces= 52,264 watertight=True
target edge=0.5000 mm (0.2586% of bbox diag=193.34 mm)
[after remesh ] verts= 58,004 faces=116,004 watertight=True
Watertight? True
Remeshed and smoothed geometry exported → fixed_mesh.stl
9. Visualise and export sections
Inspect all labeled regions of the stitched result, then export individual anatomical sections as
STL files with mm.export_section_stl.
Colour |
Region |
|---|---|
Yellow |
Aorta |
Blue |
RCA |
Green |
LCA |
Red |
Intramural wall (removed) |
Cyan |
Proximal |
Magenta |
Distal |
Orange |
Anomalous |
mesh = stitched['mesh']
verts, faces = mesh.vertices, mesh.faces
intensity = mesh.vertex_normals[:, 2]
fig = go.Figure(data=[go.Mesh3d(
x=verts[:,0], y=verts[:,1], z=verts[:,2],
i=faces[:,0], j=faces[:,1], k=faces[:,2],
intensity=intensity,
colorscale="RdBu",
opacity=1.0,
flatshading=False,
name="Stitched geometry",
)])
fig.update_layout(scene=dict(aspectmode="data"), margin=dict(l=0, r=0, t=30, b=0))
fig.show()
mm.export_section_stl(stitched, "all")
mm.export_section_stl(stitched, "aorta")
mm.export_section_stl(stitched, "rca")
mm.export_section_stl(stitched, "lca")
print("Exported: all.stl aorta.stl rca.stl lca.stl")
Exported: all.stl aorta.stl rca.stl lca.stl
10. Re-label the final geometry
After remeshing and smoothing, vertex coordinates have changed and the stored point lists are no
longer valid. Re-running mm.label_geometry on the exported fixed mesh produces a fresh, consistent
labeling for downstream biomechanical simulation or further analysis.
results_final, (rca_cl_f, lca_cl_f, ao_cl_f) = mm.label_geometry(
path_ccta_geometry="fixed_mesh.stl",
path_centerline_aorta="../data/centerline_aorta.csv",
path_centerline_rca="../data/centerline_rca_short.csv",
path_centerline_lca="../data/centerline_lca.csv",
bounding_sphere_radius_mm=3.0,
n_points_intramural=100,
anomalous_rca=True,
anomalous_lca=False,
control_plot=False,
)
plot_labeled_geometry(results_final)
mm.export_section_stl(results_final, "all")
mm.export_section_stl(results_final, "aorta")
mm.export_section_stl(results_final, "lca")
mm.export_section_stl(results_final, "rca")
print("Final labeled sections exported.")
Loaded mesh: 58004 vertices, 116004 faces
Loaded aorta centerline: 75 points
Loaded LCA centerline: 1092 points
Loaded RCA centerline: 788 points
RCA points found: 13674
LCA points found: 14130
Applying occlusion removal for anomalous RCA...
Total faces to exclude: 322
Excluded 322 faces, removed 191 points (filtered from 13674 to 13483 points)
RCA: relabeled 191 points in intramual course
Removing LCA and RCA island points...
length before: 14130
length after: 14087
Applying final reclassification based on adjacency map...
aorta_points:30243
rca_points:13483
lca_points:14087
rca_removed_points:191
lca_removed_points:0
Final labeled sections exported.