{ "cells": [ { "cell_type": "markdown", "id": "5f835a6a", "metadata": {}, "source": [ "# **multimoda-rs**: CCTA Module Tutorial\n", "\n", "This notebook demonstrates the full CCTA–intravascular fusion pipeline:\n", "\n", "1. Read in and label a CCTA geometry from an STL file and centerline CSV files\n", "2. Prepare centerlines, detect branches, and discretize the vessel tree\n", "3. Load an intravascular geometry and fine-align it to the CCTA point cloud\n", "4. Label the anomalous (intramural) region within the CCTA mesh\n", "5. Compute radial scaling factors for the proximal, distal, and aortic regions\n", "6. Morph the CCTA mesh to match the intravascular geometry\n", "7. Remove the intramural region and stitch the CCTA to the intravascular geometry\n", "8. Remesh and smooth the stitched geometry\n", "9. Visualise labeled regions and export section STL files\n", "10. Re-label the final stitched geometry\n", "\n", "> **Note**: Example data is loaded automatically from `examples/data/` when run from the repository root." ] }, { "cell_type": "markdown", "id": "8e4dfe8f", "metadata": {}, "source": [ "Load required packages and set the working directory. Performance is typically faster in the console than in a notebook.\n" ] }, { "cell_type": "code", "execution_count": 1, "id": "3e81c481", "metadata": { "execution": { "iopub.execute_input": "2026-05-04T15:10:57.242463Z", "iopub.status.busy": "2026-05-04T15:10:57.242255Z", "iopub.status.idle": "2026-05-04T15:11:15.856292Z", "shell.execute_reply": "2026-05-04T15:11:15.855584Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Working directory: /mnt/c/WorkingData/Documents/2_Coding/Rust/multimodars/examples/data\n" ] } ], "source": [ "import os\n", "from pathlib import Path\n", "import numpy as np\n", "import multimodars as mm\n", "\n", "cwd = Path.cwd()\n", "for candidate in [cwd, cwd.parent, cwd.parent.parent]:\n", " if (candidate / \"examples\" / \"data\").exists():\n", " os.chdir(candidate / \"examples\" / \"data\")\n", " break\n", " elif (candidate / \"data\").exists():\n", " os.chdir(candidate / \"data\")\n", " break\n", "print(f\"Working directory: {os.getcwd()}\")\n" ] }, { "cell_type": "code", "execution_count": 2, "id": "3e8428c5", "metadata": { "execution": { "iopub.execute_input": "2026-05-04T15:11:15.858537Z", "iopub.status.busy": "2026-05-04T15:11:15.858168Z", "iopub.status.idle": "2026-05-04T15:11:23.984837Z", "shell.execute_reply": "2026-05-04T15:11:23.983925Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Note: you may need to restart the kernel to use updated packages.\n" ] } ], "source": [ "# Install if needed\n", "%pip --disable-pip-version-check install trimesh plotly scipy pymeshlab | grep -v 'already satisfied'\n", "\n", "import trimesh\n", "import plotly.graph_objects as go\n", "import plotly.io as pio\n", "pio.renderers.default = \"notebook_connected\"\n", "\n", "_DARK = dict(\n", " scene=dict(aspectmode='data', bgcolor='black'),\n", " margin=dict(r=0, l=0, b=0, t=30),\n", " legend=dict(yanchor=\"top\", y=0.99, xanchor=\"left\", x=0.01),\n", " paper_bgcolor=\"black\",\n", ")\n", "\n", "\n", "def trimesh_to_mesh3d(mesh, color, name, opacity=0.6):\n", " verts, faces = mesh.vertices, mesh.faces\n", " return go.Mesh3d(\n", " x=verts[:,0], y=verts[:,1], z=verts[:,2],\n", " i=faces[:,0], j=faces[:,1], k=faces[:,2],\n", " color=color, opacity=opacity, name=name, flatshading=True,\n", " )\n", "\n", "\n", "def _to_arr(pts):\n", " return np.array(pts, dtype=np.float64) if pts else np.empty((0, 3))\n", "\n", "\n", "def plot_labeled_geometry(results):\n", " \"\"\"Visualise initial geometry labeling (aorta / RCA / LCA / removed).\"\"\"\n", " mesh = results[\"mesh\"]\n", " fig = go.Figure()\n", " fig.add_trace(go.Mesh3d(\n", " x=mesh.vertices[:,0], y=mesh.vertices[:,1], z=mesh.vertices[:,2],\n", " i=mesh.faces[:,0], j=mesh.faces[:,1], k=mesh.faces[:,2],\n", " opacity=0.1, color='lightgray', name='Mesh', hoverinfo='skip',\n", " ))\n", " for key, color, label in [\n", " (\"aorta_points\", \"yellow\", \"Aorta\"),\n", " (\"rca_points\", \"blue\", \"RCA\"),\n", " (\"lca_points\", \"green\", \"LCA\"),\n", " (\"rca_removed_points\", \"red\", \"RCA removed (intramural)\"),\n", " (\"lca_removed_points\", \"red\", \"LCA removed (intramural)\"),\n", " ]:\n", " arr = _to_arr(results.get(key, []))\n", " if len(arr):\n", " fig.add_trace(go.Scatter3d(\n", " x=arr[:,0], y=arr[:,1], z=arr[:,2],\n", " mode='markers', marker=dict(size=2, color=color), name=label,\n", " ))\n", " fig.update_layout(**_DARK)\n", " fig.show()\n", "\n", "\n", "def plot_anomalous_regions(results):\n", " \"\"\"Visualise proximal / anomalous / distal sub-regions.\"\"\"\n", " mesh = results[\"mesh\"]\n", " fig = go.Figure()\n", " fig.add_trace(go.Mesh3d(\n", " x=mesh.vertices[:,0], y=mesh.vertices[:,1], z=mesh.vertices[:,2],\n", " i=mesh.faces[:,0], j=mesh.faces[:,1], k=mesh.faces[:,2],\n", " opacity=0.1, color='lightgray', name='Mesh', hoverinfo='skip',\n", " ))\n", " for key, color, label in [\n", " (\"proximal_points\", \"cyan\", \"Proximal\"),\n", " (\"anomalous_points\", \"orange\", \"Anomalous\"),\n", " (\"distal_points\", \"magenta\", \"Distal\"),\n", " ]:\n", " arr = _to_arr(results.get(key, []))\n", " if len(arr):\n", " fig.add_trace(go.Scatter3d(\n", " x=arr[:,0], y=arr[:,1], z=arr[:,2],\n", " mode='markers', marker=dict(size=2, color=color), name=label,\n", " ))\n", " fig.update_layout(**_DARK)\n", " fig.show()\n", "\n", "\n", "def plot_results_key(results):\n", " \"\"\"Plotly replacement for mm.plot_results_key — all region labels on one figure.\"\"\"\n", " mesh = results[\"mesh\"]\n", " fig = go.Figure()\n", " fig.add_trace(go.Mesh3d(\n", " x=mesh.vertices[:,0], y=mesh.vertices[:,1], z=mesh.vertices[:,2],\n", " i=mesh.faces[:,0], j=mesh.faces[:,1], k=mesh.faces[:,2],\n", " opacity=0.15, color='lightgray', name='Mesh', hoverinfo='skip',\n", " ))\n", " for key, color, label in [\n", " (\"aorta_points\", \"yellow\", \"Aorta\"),\n", " (\"rca_points\", \"blue\", \"RCA\"),\n", " (\"lca_points\", \"green\", \"LCA\"),\n", " (\"rca_removed_points\", \"red\", \"Intramural wall\"),\n", " (\"proximal_points\", \"cyan\", \"Proximal\"),\n", " (\"distal_points\", \"magenta\", \"Distal\"),\n", " (\"anomalous_points\", \"orange\", \"Anomalous\"),\n", " ]:\n", " arr = _to_arr(results.get(key, []))\n", " if len(arr):\n", " fig.add_trace(go.Scatter3d(\n", " x=arr[:,0], y=arr[:,1], z=arr[:,2],\n", " mode='markers', marker=dict(size=2, color=color), name=label,\n", " ))\n", " fig.update_layout(**_DARK)\n", " fig.show()\n", "\n", "\n", "def plot_centerline_branches(rca_cl, lca_cl, results_dict=None):\n", " \"\"\"Visualise centerline branch assignments for RCA and LCA.\n", "\n", " Each branch gets its own toggleable legend entry. When *results_dict*\n", " is provided, labelled surface-mesh points are overlaid so you can verify\n", " that branch labels transferred correctly to the geometry.\n", "\n", " Colour coding: blues = RCA (dark = main, lighter = side branches);\n", " reds/oranges = LCA.\n", " \"\"\"\n", " _RCA = [\"#1f77b4\", \"#17becf\", \"#9467bd\", \"#2ca02c\", \"#7f7f7f\"]\n", " _LCA = [\"#d62728\", \"#ff7f0e\", \"#e377c2\", \"#bcbd22\", \"#8c564b\"]\n", "\n", " fig = go.Figure()\n", "\n", " def add_cl(cl, colors, vessel):\n", " from collections import defaultdict\n", " by_branch = defaultdict(list)\n", " for p in cl.points:\n", " by_branch[p.branch_id].append(p.contour_point)\n", " for bid in sorted(by_branch):\n", " pts = by_branch[bid]\n", " label = f\"{vessel} main\" if bid == 0 else f\"{vessel} branch {bid}\"\n", " fig.add_trace(go.Scatter3d(\n", " x=[p.x for p in pts], y=[p.y for p in pts], z=[p.z for p in pts],\n", " mode='markers', marker=dict(color=colors[bid % len(colors)], size=3),\n", " name=label,\n", " ))\n", "\n", " def add_surf(key, color, label):\n", " arr = _to_arr(results_dict.get(key, []))\n", " if len(arr):\n", " fig.add_trace(go.Scatter3d(\n", " x=arr[:,0], y=arr[:,1], z=arr[:,2],\n", " mode='markers', marker=dict(color=color, size=1.5, opacity=0.3),\n", " name=label,\n", " ))\n", "\n", " add_cl(rca_cl, _RCA, \"RCA\")\n", " add_cl(lca_cl, _LCA, \"LCA\")\n", "\n", " if results_dict is not None:\n", " add_surf(\"rca_points_main\", _RCA[0], \"RCA main pts\")\n", " i = 1\n", " while f\"rca_points_side_{i}\" in results_dict:\n", " add_surf(f\"rca_points_side_{i}\", _RCA[i % len(_RCA)], f\"RCA side {i} pts\")\n", " i += 1\n", " add_surf(\"lca_points_main\", _LCA[0], \"LCA main pts\")\n", " i = 1\n", " while f\"lca_points_side_{i}\" in results_dict:\n", " add_surf(f\"lca_points_side_{i}\", _LCA[i % len(_LCA)], f\"LCA side {i} pts\")\n", " i += 1\n", "\n", " fig.update_layout(title=\"Centerline Branch Assignment\", **_DARK)\n", " fig.show()\n", "\n", "\n", "def plot_vessel_tree(tree, pts_per_contour=24):\n", " \"\"\"Visualise a discretized vessel tree.\n", "\n", " Renders contour rings (subsampled to *pts_per_contour* points),\n", " yellow centroid dots tracing each vessel axis, and orientation\n", " reference triplets at the ostium and bifurcation sites.\n", "\n", " Colour coding: silver = aorta; steel-blue = RCA main; shades of blue =\n", " RCA side branches; coral = LCA main; shades of orange = LCA side branches.\n", " Red × = main reference; orange ▲ = CCW reference; magenta ▼ = CW reference.\n", " \"\"\"\n", " _RCA = [\"#1f77b4\", \"#4fa3e0\", \"#7ec8e3\", \"#a8d8ea\"]\n", " _LCA = [\"#d62728\", \"#e07f4f\", \"#e3a87e\", \"#eac0a8\"]\n", "\n", " fig = go.Figure()\n", "\n", " def ring_xyz(contour):\n", " pts = contour.points\n", " m = len(pts)\n", " if m == 0:\n", " return [], [], []\n", " step = max(1, m // pts_per_contour)\n", " s = pts[::step]\n", " return ([p.x for p in s] + [s[0].x],\n", " [p.y for p in s] + [s[0].y],\n", " [p.z for p in s] + [s[0].z])\n", "\n", " def add_contours(contours, color, name, opacity=0.7):\n", " rx, ry, rz, cx, cy, cz = [], [], [], [], [], []\n", " for c in contours:\n", " xs, ys, zs = ring_xyz(c)\n", " if not xs:\n", " continue\n", " rx += xs + [None]; ry += ys + [None]; rz += zs + [None]\n", " if c.centroid:\n", " cx.append(c.centroid[0]); cy.append(c.centroid[1]); cz.append(c.centroid[2])\n", " if rx:\n", " fig.add_trace(go.Scatter3d(x=rx, y=ry, z=rz, mode='lines',\n", " name=name, line=dict(color=color, width=1), opacity=opacity))\n", " if cx:\n", " fig.add_trace(go.Scatter3d(x=cx, y=cy, z=cz, mode='markers',\n", " marker=dict(color='yellow', size=2), showlegend=False,\n", " name=f\"{name} centroids\"))\n", "\n", " def add_refs(refs, label):\n", " mx, my, mz, ccx, ccy, ccz, clx, cly, clz = [], [], [], [], [], [], [], [], []\n", " for main_ref, cc_ref, clock_ref in refs:\n", " mx.append(main_ref[0]); my.append(main_ref[1]); mz.append(main_ref[2])\n", " ccx.append(cc_ref[0]); ccy.append(cc_ref[1]); ccz.append(cc_ref[2])\n", " clx.append(clock_ref[0]); cly.append(clock_ref[1]); clz.append(clock_ref[2])\n", " if mx:\n", " fig.add_trace(go.Scatter3d(x=mx, y=my, z=mz, mode='markers',\n", " name=f\"{label} main ref\", marker=dict(color='red', symbol='x', size=5)))\n", " fig.add_trace(go.Scatter3d(x=ccx, y=ccy, z=ccz, mode='markers',\n", " name=f\"{label} CCW ref\", marker=dict(color='orange', symbol='diamond', size=5)))\n", " fig.add_trace(go.Scatter3d(x=clx, y=cly, z=clz, mode='markers',\n", " name=f\"{label} CW ref\", marker=dict(color='magenta', symbol='square', size=5)))\n", "\n", " add_contours(tree.discretized_aorta, \"silver\", \"Aorta\", opacity=0.3)\n", " add_contours(tree.discretized_rca_main, \"steelblue\", \"RCA main\")\n", " for i, branch in enumerate(tree.rca_branches):\n", " add_contours(branch, _RCA[i % len(_RCA)], f\"RCA branch {i+1}\")\n", " add_contours(tree.discretized_lca_main, \"coral\", \"LCA main\")\n", " for i, branch in enumerate(tree.lca_branches):\n", " add_contours(branch, _LCA[i % len(_LCA)], f\"LCA branch {i+1}\")\n", " add_refs(tree.rca_references, \"RCA\")\n", " add_refs(tree.lca_references, \"LCA\")\n", "\n", " fig.update_layout(title=\"Discretized Vessel Tree\", **_DARK)\n", " fig.show()\n", "\n", "def plot_centerline_edges(cl, cos_threshold=0.0, title=\"Centerline Edges\"):\n", " \"\"\"Show centerline branches with sharp-angle positions as red x markers.\n", "\n", " Each branch gets a distinct colour. Large red x markers indicate positions\n", " from find_sharp_angles — use them to decide where to call\n", " split_branch or merge_branches.\n", " \"\"\"\n", " _PALETTE = [\n", " \"#1f77b4\", \"#ff7f0e\", \"#2ca02c\", \"#d62728\", \"#9467bd\",\n", " \"#8c564b\", \"#e377c2\", \"#7f7f7f\", \"#bcbd22\", \"#17becf\",\n", " ]\n", " from collections import defaultdict\n", " by_branch = defaultdict(list)\n", " for p in cl.points:\n", " by_branch[p.branch_id].append(p.contour_point)\n", "\n", " fig = go.Figure()\n", " for bid in sorted(by_branch):\n", " pts = by_branch[bid]\n", " color = _PALETTE[bid % len(_PALETTE)]\n", " label = \"main\" if bid == 0 else f\"branch {bid}\"\n", " fig.add_trace(go.Scatter3d(\n", " x=[p.x for p in pts], y=[p.y for p in pts], z=[p.z for p in pts],\n", " mode='lines+markers', name=label,\n", " line=dict(color=color, width=2),\n", " marker=dict(color=color, size=3),\n", " ))\n", " sharp_pos = cl.find_sharp_angles(bid, cos_threshold)\n", " if sharp_pos:\n", " sharp_pts = [pts[i] for i in sharp_pos if i < len(pts)]\n", " fig.add_trace(go.Scatter3d(\n", " x=[p.x for p in sharp_pts], y=[p.y for p in sharp_pts], z=[p.z for p in sharp_pts],\n", " mode='markers', name=f\"{label} sharp angles\",\n", " marker=dict(color='red', size=8, symbol='x'),\n", " ))\n", " fig.update_layout(title=title, **_DARK)\n", " fig.show()" ] }, { "cell_type": "markdown", "id": "68985b25", "metadata": {}, "source": [ "## 1. Read in and prepare CCTA geometries and centerlines\n", "\n", "The entry point for the CCTA module is `mm.label_geometry`, which reads a triangulated surface mesh\n", "(STL) together with centerline CSV files for the aorta, RCA, and LCA. It returns a labeled `results`\n", "dictionary and three `PyCenterline` objects used in all subsequent steps.\n", "\n", "The centerline CSV files must contain three columns (no header): `x`, `y`, `z` in mm.\n", "\n", "**Key parameters:**\n", "- `bounding_sphere_radius_mm`: radius of the rolling sphere for the initial vessel labeling pass.\n", " Larger values capture more distant vertices; smaller values are more conservative.\n", "- `n_points_intramural`: number of centerline points defining the end of the intramural segment.\n", "- `anomalous_rca` / `anomalous_lca`: when `True`, re-assigns incorrectly labeled intramural points.\n", " Set to `False` for normal coronary anatomy.\n", "- `control_plot`: open an interactive scene to inspect labeling. Useful when tuning the parameters above.\n", "\n", "The returned `results` dictionary contains keys `\"mesh\"`, `\"aorta_points\"`, `\"rca_points\"`,\n", "`\"lca_points\"`, `\"rca_removed_points\"`, and `\"lca_removed_points\"`.\n" ] }, { "cell_type": "code", "execution_count": 3, "id": "b33a0ed3", "metadata": { "execution": { "iopub.execute_input": "2026-05-04T15:11:23.987453Z", "iopub.status.busy": "2026-05-04T15:11:23.987223Z", "iopub.status.idle": "2026-05-04T15:11:35.123454Z", "shell.execute_reply": "2026-05-04T15:11:35.122613Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Loaded mesh: 25171 vertices, 50338 faces\n", "Loaded aorta centerline: 75 points\n", "Loaded LCA centerline: 1092 points\n", "Loaded RCA centerline: 788 points\n", "\n", "RCA points found: 4820\n", "LCA points found: 5899\n", "Applying occlusion removal for anomalous RCA...\n", "Total faces to exclude: 232\n", "Excluded 232 faces, removed 110 points (filtered from 4820 to 4710 points)\n", "RCA: relabeled 110 points in intramual course\n", "\n", "Removing LCA and RCA island points...\n", "length before: 5899\n", "length after: 5879\n", "\n", "Applying final reclassification based on adjacency map...\n", "aorta_points:14472\n", "rca_points:4710\n", "lca_points:5879\n", "rca_removed_points:110\n", "lca_removed_points:0\n" ] }, { "data": { "text/html": [ " \n", " \n", " " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "