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/d/00_coding/multimoda-rs/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.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. Ignored whenvtp_data=True(see section below).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
Alternative: VTP centerlines from 3D-Slicer / VMTK
If centerlines were exported as .vtp files (ASCII VTK PolyData) by 3D-Slicer or VMTK, the
entire preparation workflow simplifies significantly. read_centerline_vtp parses the file and
identifies branches automatically; cleanup_vtp_data removes the run-alongside-main-branch
prefix that VTP files attach to every side branch and optionally smooths the result.
prepare_centerlines(vtp_data=True) then skips calculate_branches (branch indices are already
set) and only runs check_centerline to normalise ordering — no manual corrections needed.
# Load and clean VTP centerlines
rca_cl = mm.read_centerline_vtp("./rca_cl.vtp").cleanup_vtp_data(smooth=True)
lca_cl = mm.read_centerline_vtp("./lca_cl.vtp").cleanup_vtp_data(smooth=True)
ao_cl = mm.read_centerline_vtp("./ao_cl.vtp").cleanup_vtp_data(smooth=True)
# Re-run label_geometry with PyCenterline objects instead of file paths
results, (rca_cl, lca_cl, ao_cl) = mm.label_geometry(
path_ccta_geometry="./NARCO_119.stl",
path_centerline_aorta=ao_cl,
path_centerline_rca=rca_cl,
path_centerline_lca=lca_cl,
bounding_sphere_radius_mm=3.0,
n_points_intramural=200,
anomalous_rca=True,
anomalous_lca=False,
control_plot=False,
)
# Branch indices already set from VTP — skip calculate_branches
rca_cl, lca_cl, results = mm.prepare_centerlines(
rca_cl, lca_cl, results, vtp_data=True
)
plot_centerline_branches(rca_cl, lca_cl, results)
Loaded mesh: 25171 vertices, 50338 faces
Using provided aorta centerline: 150 points
Using provided LCA centerline: 2074 points
Using provided RCA centerline: 1098 points
RCA points found: 4806
LCA points found: 5901
Applying occlusion removal for anomalous RCA...
Total faces to exclude: 243
Excluded 243 faces, removed 111 points (filtered from 4806 to 4695 points)
RCA: relabeled 111 points in intramual course
Removing LCA and RCA island points...
length before: 5901
length after: 5882
Applying final reclassification based on adjacency map...
aorta_points:14483
rca_points:4695
lca_points:5882
rca_removed_points:111
lca_removed_points:0
Branch labeling for 'rca_points' (branch_ids=[0]):
rca_points_main: 3734
rca_points_side: 961
rca_points_side_1: 23
rca_points_side_2: 515
rca_points_side_3: 423
Branch labeling for 'lca_points' (branch_ids=[0]):
lca_points_main: 2627
lca_points_side: 3255
lca_points_side_1: 2616
lca_points_side_2: 2275
lca_points_side_3: 278
lca_points_side_4: 115
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: 222.00°
Step 2: Refining with Hausdorff distance
---------------------Refining alignment with Hausdorff---------------------
Initial rotation: 0.00°, Initial CL index: 31
Refined rotation: -5.00°, Refined CL index: 32, Hausdorff: 4.99
---------------------Applying final transformation---------------------
Total rotation (initial + delta): 217.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
calcification file not found, skipping: "ivus_rest/calcium_diastolic_contours.csv"
process_directory: unknown mapping name 'catheter', skipping
sidebranch file not found, skipping: "ivus_rest/branch_diastolic_contours.csv"
eem file not found, skipping: "ivus_rest/eem_diastolic_contours.csv"
eem file not found, skipping: "ivus_rest/eem_systolic_contours.csv"
sidebranch file not found, skipping: "ivus_rest/branch_systolic_contours.csv"
process_directory: unknown mapping name 'catheter', skipping
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=199.93001601905405, spacing=1.295576
resample_centerline_by_contours: produced 155 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: 38
distal_points: 4243
anomalous_points: 414
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.5 mm
Distal scaling: 0.4 mm
Finding best aortic radial scaling factor...
Aortic scaling: -0.7 mm
Finding best aortic wall radial scaling factor...
elliptic ratio <1.3 for frame index 7
Aortic wall scaling: 1.35 mm
Proximal scaling: -0.500 mm
Distal scaling: 0.400 mm
Aortic scaling: -0.700 mm
Aortic wall scaling: 1.353 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 4243 vertices around Centerline(len=1098, spacing=0.26 mm, branches=4)
Diameter adjustment: 0.4 mm
Scaling 14594 vertices around Centerline(len=75, spacing=0.75 mm, branches=1)
Diameter adjustment: -0.7 mm
Scaling 38 vertices around Centerline(len=1098, spacing=0.26 mm, branches=4)
Diameter adjustment: -0.5 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 452
Created 55 boundary points
Stitching: 100/100 triangles created (n_boundary=37, n_iv=100, step=2, remainder=26)
Stitching: 100/100 triangles created (n_boundary=18, n_iv=100, step=5, remainder=10)
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,119 faces= 52,234 watertight=False
non-manifold edges/vertices repaired
holes closed
[after hole fill ] verts= 26,119 faces= 52,236 watertight=False
target edge=0.5000 mm (0.2588% of bbox diag=193.22 mm)
[after remesh ] verts= 57,247 faces=114,492 watertight=False
[after post-remesh fix ] verts= 57,248 faces=114,492 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: 57248 vertices, 114492 faces
Loaded aorta centerline: 75 points
Loaded LCA centerline: 1092 points
Loaded RCA centerline: 788 points
RCA points found: 13612
LCA points found: 14147
Applying occlusion removal for anomalous RCA...
Total faces to exclude: 191
Excluded 191 faces, removed 135 points (filtered from 13612 to 13477 points)
RCA: relabeled 135 points in intramual course
Removing LCA and RCA island points...
length before: 14147
length after: 14111
Applying final reclassification based on adjacency map...
aorta_points:29525
rca_points:13477
lca_points:14111
rca_removed_points:135
lca_removed_points:0
Final labeled sections exported.