01 — Ellipsoid¶
Bundle Mesh → Grow → Activate → Export¶
End-to-end demo using the bundled data/ellipsoid.obj:
Load/validate the mesh
Grow a fractal tree on the ellipsoid
Wrap into a
PurkinjeTreeand run FIM activationExport VTK for ParaView and render inline
Outputs
output/examples/01_ellipsoid/ellipsoid_AT.vtuoutput/examples/01_ellipsoid/ellipsoid_meshio.vtuoutput/examples/01_ellipsoid/ellipsoid_pmj.vtu(PMJs as points)
[1]:
# Install core runtime dependencies.
%pip install -q --upgrade pip purkinje-uv plotly
[2]:
from pathlib import Path
import os
import numpy as np
from purkinje_uv import FractalTreeParameters, FractalTree, PurkinjeTree
import pyvista as pv
import plotly.graph_objects as go
import urllib.request
import logging
[3]:
SEED = int(os.getenv("EXAMPLES_SEED", "1234"))
LITE = bool(int(os.getenv("EXAMPLES_LITE", "1")))
AUTO_SEEDS = bool(int(os.getenv("AUTO_SEEDS", "0"))) # 0 = use known seeds; 1 = compute
SHOW = bool(int(os.getenv("EXAMPLES_SHOW", "1"))) # viewer cell later
OUT_DIR = Path("output") / "examples" / "01_ellipsoid"
OUT_DIR.mkdir(parents=True, exist_ok=True)
DATA_DIR = Path("data")
DATA_DIR.mkdir(parents=True, exist_ok=True)
np.random.seed(SEED)
print(f"SEED={SEED} LITE={LITE} AUTO_SEEDS={AUTO_SEEDS} SHOW={SHOW}")
print("OUT_DIR:", OUT_DIR)
logging.basicConfig(level=logging.INFO)
logger = logging.getLogger("purkinje_uv.examples.ellipsoid_demo")
SEED=1234 LITE=True AUTO_SEEDS=False SHOW=True
OUT_DIR: output/examples/01_ellipsoid
Load & Validate Mesh¶
We use the bundled data/ellipsoid.obj.
[4]:
mesh_path = DATA_DIR / "ellipsoid.obj"
mesh_path.parent.mkdir(parents=True, exist_ok=True)
# Download mesh if it's not on the DATA_DIR
if not mesh_path.exists():
url = "https://raw.githubusercontent.com/ricardogr07/purkinje-uv/main/data/ellipsoid.obj"
last_err = None
try:
print("Trying to download from:", url)
with urllib.request.urlopen(url, timeout=20) as r, open(mesh_path, "wb") as f:
f.write(r.read())
print("Downloaded:", mesh_path, "| size:", mesh_path.stat().st_size, "bytes")
except Exception as e:
last_err = e
print(" failed:", e)
# Load + quick validate (report openness; no auto-open needed)
surf = pv.read(str(mesh_path))
if not isinstance(surf, pv.PolyData):
surf = surf.extract_geometry()
surf = surf.clean().triangulate()
edges = surf.extract_feature_edges(
boundary_edges=True, feature_edges=False, manifold_edges=False, non_manifold_edges=True
)
is_open = edges.n_cells > 0
print("Ellipsoid faces:", surf.n_cells, "| points:", surf.n_points, "| is_open:", is_open)
print(f"Using mesh: {mesh_path}")
Ellipsoid faces: 1470 | points: 780 | is_open: True
Using mesh: data/ellipsoid.obj
Seed Strategy (Manual vs. Auto)¶
Manual (default): use known good IDs (
738,210) for this ellipsoid.
[5]:
init_node_id, second_node_id = 738, 210
print("Seeds → init_node_id:", init_node_id, " second_node_id:", second_node_id)
Seeds → init_node_id: 738 second_node_id: 210
Parameters¶
Balanced defaults that work well on the ellipsoid:
l_segment=0.01init_length=0.30,length=0.15branch_angle=0.15 rad,w=0.1N_it=5(LITE) or10(non-LITE)Optional fascicles for broader coverage:
[-0.4, 0.5]with lengths[20*lseg, 40*lseg].
[6]:
lseg = 0.01
N_IT_LITE, N_IT_FULL = 5, 10
n_it = N_IT_LITE if LITE else N_IT_FULL
params = FractalTreeParameters(
meshfile=str(mesh_path),
init_node_id=738,
second_node_id=210,
l_segment=lseg,
init_length=0.3,
length=0.15,
fascicles_length=[20 * lseg, 40 * lseg],
fascicles_angles=[-0.4, 0.5], # radians
)
# Save params snapshot
params.to_json_file(OUT_DIR / "params.json")
print("Saved params to:", OUT_DIR / "params.json")
print(params.to_json())
Saved params to: output/examples/01_ellipsoid/params.json
{
"meshfile": "data/ellipsoid.obj",
"init_node_id": 738,
"second_node_id": 210,
"init_length": 0.3,
"N_it": 10,
"length": 0.15,
"branch_angle": 0.15,
"w": 0.1,
"l_segment": 0.01,
"fascicles_angles": [
-0.4,
0.5
],
"fascicles_length": [
0.2,
0.4
]
}
Grow the Fractal Tree¶
FractalTree loads the mesh, performs UV mapping/scaling, and generates:
nodes_xyz(points),connectivity(edges),end_nodes(PMJs).
[7]:
ft = FractalTree(params=params)
ft.grow_tree()
print(f"Tree grown. Nodes={len(ft.nodes_xyz)} Edges={len(ft.connectivity)} PMJs={len(ft.end_nodes)}")
computing uv map
generation 0
generation 1
generation 2
generation 3
generation 4
generation 5
generation 6
generation 7
generation 8
generation 9
Tree grown. Nodes=3250 Edges=3249 PMJs=143
PurkinjeTree + FIM Activation¶
We wrap the tree, stimulate node 0 at time 0.0, and compute activation for all nodes.
[8]:
P = PurkinjeTree(
nodes=np.asarray(ft.nodes_xyz, dtype=float),
connectivity=np.asarray(ft.connectivity, dtype=int),
end_nodes=np.asarray(ft.end_nodes, dtype=int),
)
act = P.activate_fim(
x0=np.array([0], dtype=int),
x0_vals=np.array([0.0], dtype=float),
return_only_pmj=False,
)
print("Activation:", act.shape, "| min/max:", float(act.min()), float(act.max()))
print("PMJs:", len(P.pmj) if hasattr(P, "pmj") else 0)
Activation: (3250,) | min/max: 0.0 1.4677733182907104
PMJs: 143
Export Files¶
Write three files under OUT_DIR:
ellipsoid_AT.vtu— unstructured grid + activation times,ellipsoid_meshio.vtu— line mesh viameshio,ellipsoid_pmj.vtu— PMJs as points (use.vtufor broad compatibility).
[9]:
at_path = OUT_DIR / "ellipsoid_AT.vtu"
meshio_path = OUT_DIR / "ellipsoid_meshio.vtu"
pmj_path = OUT_DIR / "ellipsoid_pmj.vtu"
P.save(str(at_path))
P.save_meshio(str(meshio_path))
P.save_pmjs(str(pmj_path))
print("Wrote:")
print(" -", at_path)
print(" -", meshio_path)
print(" -", pmj_path)
Wrote:
- output/examples/01_ellipsoid/ellipsoid_AT.vtu
- output/examples/01_ellipsoid/ellipsoid_meshio.vtu
- output/examples/01_ellipsoid/ellipsoid_pmj.vtu
Inline Viewer (OBJ overlay + Tree + PMJs)¶
Renders the ellipsoid surface (slightly transparent), the tree as tubes, and PMJs as points.
[10]:
# Inline 3D viewer in Colab using Plotly
if SHOW:
# Paths from earlier steps
out_dir = Path("output/examples/01_ellipsoid")
tree_vtu = out_dir / "ellipsoid_AT.vtu"
pmj_vtu = out_dir / "ellipsoid_pmj.vtu" # optional
obj_path = Path("data/ellipsoid.obj") # optional
# ---- helpers ----
def get_point_scalars(ds: pv.DataSet):
preferred = ["Activation", "AT", "activation", "activation_time", "time_activation"]
for nm in preferred:
if nm in ds.point_data:
return ds.point_data[nm], nm
if len(ds.point_data.keys()) > 0:
k = list(ds.point_data.keys())[0]
return ds.point_data[k], k
return None, None
def polydata_from(ds: pv.DataSet) -> pv.PolyData:
if isinstance(ds, pv.PolyData):
return ds
try:
pd = ds.extract_geometry()
if isinstance(pd, pv.PolyData):
return pd
except Exception:
pass
return pv.PolyData(ds.points) if hasattr(ds, "points") else pv.PolyData()
def parse_lines(pd: pv.PolyData):
"""Return list of (u,v) index pairs from PolyData 'lines' array."""
if pd.lines.size == 0:
return []
arr = pd.lines.reshape(-1)
edges = []
i = 0
n = arr.size
while i < n:
m = int(arr[i]); i += 1
ids = arr[i:i+m]; i += m
# convert each polyline into segments (id[j], id[j+1])
for j in range(m-1):
edges.append((int(ids[j]), int(ids[j+1])))
return edges
# -----------------
# Load tree dataset (unstructured grid with line cells + activation)
tree_ds = pv.read(str(tree_vtu))
tree_pd = polydata_from(tree_ds) # PolyData with lines
tree_pts = tree_pd.points
edges = parse_lines(tree_pd)
# Scalars for coloring (point data)
scalars, scal_name = get_point_scalars(tree_ds)
# Optional PMJs
pmj_pts = None
if pmj_vtu.exists():
try:
pmj_ds = pv.read(str(pmj_vtu))
pmj_pd = polydata_from(pmj_ds)
pmj_pts = pmj_pd.points if pmj_pd.n_points else None
except Exception:
pmj_pts = None
# Optional surface
surf_pts = surf_tris = None
if obj_path.exists():
try:
surf = pv.read(str(obj_path))
if not isinstance(surf, pv.PolyData):
surf = surf.extract_geometry()
surf = surf.clean().triangulate()
if surf.n_cells:
faces = surf.faces.reshape(-1, 4)[:, 1:] # tri indices (i,j,k)
surf_pts = surf.points
surf_tris = faces
except Exception:
pass
# Build Plotly figure
fig = go.Figure()
# Surface mesh (semi-transparent)
if surf_pts is not None and surf_tris is not None:
fig.add_trace(go.Mesh3d(
x=surf_pts[:,0], y=surf_pts[:,1], z=surf_pts[:,2],
i=surf_tris[:,0], j=surf_tris[:,1], k=surf_tris[:,2],
opacity=0.18, name="Surface", color="lightgray",
lighting=dict(ambient=0.6, diffuse=0.6, specular=0.2)
))
# Tree as line segments (single color; fast)
if edges:
# Build one big line trace with NaN breaks
xs, ys, zs = [], [], []
for u, v in edges:
xs += [tree_pts[u,0], tree_pts[v,0], None]
ys += [tree_pts[u,1], tree_pts[v,1], None]
zs += [tree_pts[u,2], tree_pts[v,2], None]
fig.add_trace(go.Scatter3d(
x=xs, y=ys, z=zs, mode="lines",
line=dict(width=3),
name="Purkinje tree"
))
# Nodes colored by activation (shows the scalar field)
if scalars is not None:
fig.add_trace(go.Scatter3d(
x=tree_pts[:,0], y=tree_pts[:,1], z=tree_pts[:,2],
mode="markers",
marker=dict(
size=3,
color=scalars, colorscale="Viridis",
colorbar=dict(title=scal_name or "activation"),
opacity=0.9
),
name="Nodes (colored)"
))
# PMJs as small squares
if pmj_pts is not None:
fig.add_trace(go.Scatter3d(
x=pmj_pts[:,0], y=pmj_pts[:,1], z=pmj_pts[:,2],
mode="markers",
marker=dict(size=4, symbol="square", color="crimson", opacity=0.95),
name="PMJs"
))
fig.update_layout(
scene=dict(aspectmode="data"),
margin=dict(l=0, r=0, t=0, b=0),
showlegend=True
)
fig.show()