00 — Basic Tree

Programmatic Mesh → Grow → Activate → Export

Minimal, self-contained demo that:

1. Generates a small surface mesh programmatically

2. Grows a fractal tree on it

3. Wraps it into a PurkinjeTree and runs FIM activation

4. Exports VTK files for visualization.

Why start here? No external data is required. This notebook is the fastest path to confirm the pipeline works end-to-end on your machine.

Outputs

• output/examples/00_basic_tree/basic_tree_AT.vtu — unstructured grid with activation time field.

• output/examples/00_basic_tree/basic_tree_meshio.vtu — line mesh via meshio (compatible with ParaView).

• output/examples/00_basic_tree/basic_tree_pmj.vtp — PMJs as a point set.

You can keep runtimes small with EXAMPLES_LITE=1 (default in this notebook).

%pip install -q --upgrade pip
%pip install -q purkinje-uv
%pip install "pyvista[jupyter]"

Settings, Reproducibility, and Output Paths

We use:

• SEED for determinism where applicable.

• LITE toggle to keep the example small and fast by default.

• A dedicated OUT_DIR under output/examples/00_basic_tree/ to avoid clutter.

You may override these with environment variables:

• EXAMPLES_SEED=2025

• EXAMPLES_LITE=0

• EXAMPLES_SHOW=1 (optional 3D viewer at the end)

from pathlib import Path
import os
import numpy as np
import pyvista as pv
from IPython.display import HTML, display

from purkinje_uv import FractalTreeParameters, FractalTree, PurkinjeTree

SEED = int(os.getenv("EXAMPLES_SEED", "1234"))
LITE = bool(int(os.getenv("EXAMPLES_LITE", "1")))
SHOW = bool(int(os.getenv("EXAMPLES_SHOW", "1")))  # optional PyVista preview

OUT_DIR = Path("output") / "examples" / "00_basic_tree"
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}  SHOW={SHOW}")
print("OUT_DIR:", OUT_DIR)

Create a Programmatic Mesh (Sphere)

We generate a simple triangulated sphere with PyVista and save it as an OBJ under data/. This avoids any external files and ensures the example is portable.

# Generate an OPEN ellipsoid-like surface and save as OBJ
mesh_path = Path("data/basic_open_ellipsoid.obj")
mesh_path.parent.mkdir(parents=True, exist_ok=True)

# 1) Start with a triangulated sphere
sphere = pv.Sphere(radius=1.0, theta_resolution=160, phi_resolution=160).triangulate()

# 2) Scale non-uniformly → ellipsoid-like surface
ellip = sphere.copy(deep=True)
ellip.scale([1.30, 1.00, 0.75])

# 3) Clip with a plane to create an OPEN boundary (no capping)
open_surf = ellip.clip(normal=(0, 0, 1), origin=(0, 0, 0.35), invert=False)
open_surf = open_surf.clean().triangulate()

# 4) Verify it's actually open: boundary edges must be > 0
edges = open_surf.extract_feature_edges(
    boundary_edges=True, feature_edges=False, manifold_edges=False, non_manifold_edges=True
)
n_boundary = edges.n_cells
print("Boundary edges (first try):", n_boundary)

# If somehow still closed, flip the keep-side and try again
if n_boundary == 0:
    open_surf = ellip.clip(normal=(0, 0, 1), origin=(0, 0, 0.35), invert=True).clean().triangulate()
    edges = open_surf.extract_feature_edges(
        boundary_edges=True, feature_edges=False, manifold_edges=False, non_manifold_edges=True
    )
    n_boundary = edges.n_cells
    print("Boundary edges (after invert):", n_boundary)
    assert n_boundary > 0, "Surface still closed; adjust clip origin (e.g., 0.3–0.6) or invert."

# 5) Save OBJ (PyVista will use meshio under the hood if needed)
open_surf.save(mesh_path)
print(f"Saved OPEN surface: {mesh_path} | points={open_surf.n_points} faces={open_surf.n_cells}")

Fractal Tree Parameters (Small & CI-friendly)

We start with conservative values that run quickly:

• l_segment=0.01

• init_length=0.25, length=0.12 (initial and subsequent segment lengths)

• branch_angle=0.15 rad, w=0.1 (branching control)

• init_node_id=0, second_node_id=1 (work well on the generated sphere)

If you disable LITE mode, you can increase size/complexity later.

lseg = 0.01
params = FractalTreeParameters(
    meshfile=str(mesh_path),
    l_segment=lseg,
    init_length=0.30,
    length=0.16,
    branch_angle=0.22,     # wider fan-out
    w=0.08,                # small curvature tweak
    N_it=10,                # deeper branching
    fascicles_angles=[-0.6, -0.2, 0.2, 0.6],      # radians
    fascicles_length=[25*lseg, 35*lseg, 35*lseg, 25*lseg],
)

params.to_json_file(OUT_DIR / "params.json")
print("Saved params to", OUT_DIR / "params.json")
print(params.to_json())

Grow the Fractal Tree

We instantiate FractalTree with our parameters and call grow_tree(). The class loads the OBJ, performs the internal UV mapping/scaling, and generates:

• nodes_xyz (list of 3D points),

• connectivity (list of edges),

• end_nodes (PMJ indices).

ft = FractalTree(params=params)
ft.grow_tree()

n_nodes = len(ft.nodes_xyz)
n_edges = len(ft.connectivity)
n_pmj   = len(ft.end_nodes)

print(f"Tree grown.")
print(f"Nodes={n_nodes}  Edges={n_edges}  PMJs={n_pmj}")

Wrap into PurkinjeTree and Run FIM Activation

We now build a PurkinjeTree and run the Fast Iterative Method (FIM).

• We stimulate node 0 at time 0.0.

• return_only_pmj=False returns activation times for all nodes.

You can inspect PMJ activation times via act[P.pmj].

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,
)

act_min = float(np.min(act)) if act.size else float("nan")
act_max = float(np.max(act)) if act.size else float("nan")

print("Activation array:", act.shape, " min/max:", act_min, act_max)
if hasattr(P, "pmj"):
    pmj_count = len(getattr(P, "pmj"))
    print("PMJ count:", pmj_count)

Export for Visualization

We save three files into OUT_DIR:

• *_AT.vtu — Unstructured grid with activation times (native writer).

• *_meshio.vtu — Line mesh via meshio.

• *_pmj.vtp — PMJs as points.

Open them in ParaView to verify visually.

at_path    = OUT_DIR / "basic_tree_AT.vtu"
meshio_path= OUT_DIR / "basic_tree_meshio.vtu"
pmj_path   = OUT_DIR / "basic_tree_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)

Optional: Quick 3D Preview (PyVista)

If you set EXAMPLES_SHOW=1, the cell below displays:

• Tree lines,

• PMJ points.

This is optional and may require a local 3D backend.

# --- Colab-safe PyVista visualization: OBJ surface + tree + PMJs ---

# Resolve paths from notebook globals or sensible defaults
out_dir   = Path(globals().get("OUT_DIR", "output/examples/00_basic_tree"))
at_path   = Path(globals().get("at_path", out_dir / "basic_tree_AT.vtu"))
pmj_path  = Path(globals().get("pmj_path", out_dir / "basic_tree_pmj.vtu"))

# Try to reuse the exact mesh you used to grow the tree
mesh_path = globals().get("mesh_path", None)
if mesh_path is None:
    # Fallback chain (adjust to your generation cell name)
    for candidate in [
        "data/basic_open_ellipsoid.obj",
        "data/basic_open_surface.obj",
        "data/basic_sphere.obj",
    ]:
        if Path(candidate).exists():
            mesh_path = candidate
            break
if mesh_path is None:
    raise SystemExit("Could not locate the OBJ used for growth. Set `mesh_path` before running.")

mesh_path = Path(mesh_path)

print("Reading tree VTU:", at_path)
tree_ds = pv.read(str(at_path))

pmj_ds = None
if pmj_path.exists():
    try:
        pmj_ds = pv.read(str(pmj_path))
        print("Reading PMJs:", pmj_path)
    except Exception as e:
        print("PMJ read failed; overlay will be skipped:", e)

print("Reading surface OBJ:", mesh_path)
surface = pv.read(str(mesh_path))  # expect PolyData
if not isinstance(surface, pv.PolyData):
    surface = surface.extract_geometry()
surface = surface.clean()
try:
    surface = surface.compute_normals(auto_orient_normals=True, splitting=False)
except Exception:
    pass

# Optional thinning of heavy surfaces (0.0 keeps full detail)
DECIMATE = float(os.getenv("EXAMPLES_DECIMATE", "0.0"))
if DECIMATE > 0.0 and surface.n_cells > 0:
    try:
        surface = surface.decimate_proportion(DECIMATE)
        print(f"Decimated surface to ~{int((1-DECIMATE)*100)}% faces.")
    except Exception:
        pass

def pick_scalar(ds: pv.DataSet):
    preferred = ["Activation", "AT", "activation", "activation_time", "time_activation"]
    names = set(list(ds.point_data.keys()) + list(ds.cell_data.keys()))
    for nm in preferred:
        if nm in names:
            return nm
    if len(ds.point_data.keys()) > 0:
        return list(ds.point_data.keys())[0]
    if len(ds.cell_data.keys()) > 0:
        return list(ds.cell_data.keys())[0]
    return None

def ensure_polydata(ds: pv.DataSet) -> pv.PolyData:
    try:
        pd = ds.extract_geometry()
        if isinstance(pd, pv.PolyData) and (pd.n_cells > 0 or pd.n_points > 0):
            return pd
    except Exception:
        pass
    try:
        edges = ds.extract_all_edges()
        if isinstance(edges, pv.PolyData) and edges.n_cells > 0:
            return edges
    except Exception:
        pass
    return pv.PolyData(ds.points) if hasattr(ds, "points") else pv.PolyData()

def to_tubes(pd: pv.PolyData, radius=0.003):
    try:
        if pd.n_cells > 0:
            return pd.tube(radius=radius)
    except Exception:
        pass
    return pd

# Build tree geometry as tubes for better visibility
tree_geom = ensure_polydata(tree_ds)
tree_tubes = to_tubes(tree_geom, radius=0.003)
tree_scalar = pick_scalar(tree_ds)

# Extract PMJ points if present
pmj_geom = None
if pmj_ds is not None:
    pmj_geom = ensure_polydata(pmj_ds) if getattr(pmj_ds, "n_points", 0) > 0 else None

# Also show the boundary edge loop of the OBJ (helps confirm it's OPEN)
try:
    boundary_edges = surface.extract_feature_edges(
        boundary_edges=True, feature_edges=False, manifold_edges=False, non_manifold_edges=True
    )
    has_boundary = boundary_edges.n_cells > 0
except Exception:
    has_boundary = False
    boundary_edges = None

# Try inline HTML backend first
try:
    pv.global_theme.jupyter_backend = "html"
    print("Using PyVista backend:", pv.global_theme.jupyter_backend)

    p = pv.Plotter(notebook=True)
    # Surface first (slight opacity so tubes are visible)
    p.add_mesh(surface, opacity=0.2, smooth_shading=True)
    if has_boundary:
        p.add_mesh(boundary_edges, line_width=2)

    # Tree tubes with scalar if present
    p.add_mesh(
        tree_tubes,
        scalars=tree_scalar if tree_scalar else None,
        show_scalar_bar=bool(tree_scalar),
        scalar_bar_args={"title": tree_scalar or ""},
    )

    # PMJs overlay as points
    if pmj_geom is not None and pmj_geom.n_points > 0:
        p.add_mesh(pmj_geom, style="points", point_size=12)

    p.show(title="OBJ Surface + Purkinje Tree + PMJs")
except Exception as e_html_backend:
    print("Inline HTML backend render failed:", e_html_backend)

    # Fallback 1: standalone HTML export and embed
    try:
        html_path = out_dir / "basic_tree_with_surface.html"

        try:
            pv.start_xvfb()
        except Exception:
            pass

        p = pv.Plotter(off_screen=True)
        p.add_mesh(surface, opacity=0.2, smooth_shading=True)
        if has_boundary:
            p.add_mesh(boundary_edges, line_width=2)
        p.add_mesh(
            tree_tubes,
            scalars=tree_scalar if tree_scalar else None,
            show_scalar_bar=bool(tree_scalar),
            scalar_bar_args={"title": tree_scalar or ""},
        )
        if pmj_geom is not None and pmj_geom.n_points > 0:
            p.add_mesh(pmj_geom, style="points", point_size=12)

        p.export_html(str(html_path))
        display(HTML(html_path.read_text(encoding="utf-8")))
        print("Embedded exported HTML:", html_path)
    except Exception as e_export:
        print("HTML export/embed failed:", e_export)

        # Fallback 2: screenshot
        try:
            img_path = out_dir / "basic_tree_with_surface.png"
            p = pv.Plotter(off_screen=True)
            p.add_mesh(surface, opacity=0.2, smooth_shading=True)
            if has_boundary:
                p.add_mesh(boundary_edges, line_width=2)
            p.add_mesh(
                tree_tubes,
                scalars=tree_scalar if tree_scalar else None,
                show_scalar_bar=bool(tree_scalar),
                scalar_bar_args={"title": tree_scalar or ""},
            )
            if pmj_geom is not None and pmj_geom.n_points > 0:
                p.add_mesh(pmj_geom, style="points", point_size=12)
            p.show(screenshot=str(img_path))
            from IPython.display import Image
            display(Image(filename=str(img_path)))
            print("Screenshot fallback saved:", img_path)
        except Exception as e_ss:
            print("Screenshot fallback failed:", e_ss)

Verification Checklist

• This notebook produced three files in output/examples/00_basic_tree/.

• Activation array has the same length as the number of tree nodes.

• PMJ file exists and contains points.

Reproducibility tips

• Keep SEED fixed to stabilize stochastic parts.

• Use EXAMPLES_LITE=1 for CI or quick runs; set EXAMPLES_LITE=0 to explore larger trees.

Troubleshooting

• ``vtk`` wheel installation errors (Windows/Linux): ensure you’re using a supported Python version and recent pip (python -m pip install --upgrade pip).

• OBJ writer not available: the notebook falls back to meshio to write OBJ or VTU. Ensure meshio is installed (the first cell does this).

• Blank preview window: set EXAMPLES_SHOW=1 to enable the viewer and ensure your environment supports interactive windows (e.g., not running headless).