Skip to content

API Reference

Public entry points and main modules.

Package

Build ellipsoidal cardiac geometry: tetrahedral mesh, surfaces, fibers, Purkinje network.

plot_all(out, off_screen=False)

Plot in one image: mesh, surfaces, Purkinje, and fibers.

Source code in gdellipsoid/plot.py 
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
def plot_all(out, off_screen=False):
    """Plot in one image: mesh, surfaces, Purkinje, and fibers."""
    with gd.Container(out) as ct:
        # Required: tetra and at least one surface
        tetra = Mesh.load(ct / "tetra.vtk")
        regular_surf = Mesh.load(ct / "regular_surf.vtk")

        # Optional: endo, epi, points, fibers, purkinje
        endo = _load_optional(ct, "endo.vtk")
        epi = _load_optional(ct, "epi.vtk")
        points_mesh = _load_optional(ct, "points.vtk")
        purkinje = _load_optional(ct, "purkinje_mesh.vtk")
        # mesh_cell_fibers.vtk is PolyData (point cloud); load with PyVista
        mesh_cell_fibers_pv = None
        if Path(ct / "mesh_cell_fibers.vtk").exists():
            mesh_cell_fibers_pv = pv.read(str(ct / "mesh_cell_fibers.vtk"))

        with Plotter(off_screen=off_screen) as p:
            # Volume mesh: transparent or wireframe so interior is visible
            p.add_mesh(tetra, opacity=0.12, show_edges=False, color="lightgray")

            # Surfaces: combined regular surface + optional endo/epi
            p.add_mesh(regular_surf, opacity=0.5, show_edges=False, color="white")

            if endo is not None:
                p.add_mesh(endo, opacity=0.6, show_edges=False, color="coral")
            if epi is not None:
                p.add_mesh(epi, opacity=0.6, show_edges=False, color="lightblue")

            p.camera_position = CAM_POS
            p.render()
            if not off_screen:
                p.show()
            p.screenshot(ct / "surfaces_img.png")

        # Purkinje network
        if purkinje is not None:
            with Plotter(off_screen=off_screen) as p:
                p.add_mesh(endo)
                p.add_mesh(purkinje, color="red", line_width=3)

                p.camera_position = [
                    (0.6051986645911922, -0.46423570086746, 37.204488072156266),
                    (
                        -0.0003440849696247916,
                        -5.655034101526013e-05,
                        -5.999310324980369,
                    ),
                    (0.004839410606533647, -0.9999298485328891, -0.010811018359763551),
                ]
                p.render()
                if not off_screen:
                    p.show()
                p.screenshot(ct / "purkinje_img.png")

        # Fibers: vectors at cell centers (subsample for clarity)
        if (
            mesh_cell_fibers_pv is not None
            and "fiber" in mesh_cell_fibers_pv.point_data
        ):
            with Plotter(off_screen=off_screen) as p:
                p.add_mesh(tetra, opacity=0.12, show_edges=False, color="lightgray")
                if endo is not None:
                    p.add_mesh(endo, opacity=0.6, show_edges=False, color="coral")
                if epi is not None:
                    p.add_mesh(epi, opacity=0.6, show_edges=False, color="lightblue")

                pf = mesh_cell_fibers_pv
                n_pts = pf.n_points
                step = max(1, n_pts // 800)
                subsampled = pf.copy()
                subsampled.points = np.asarray(pf.points)[::step]
                for key in ("fiber", "sheet"):
                    if key in pf.point_data.keys():
                        subsampled.point_data[key] = np.asarray(pf.point_data[key])[
                            ::step
                        ]
                subsampled.set_active_vectors("fiber")
                arrows = subsampled.glyph(orient="fiber", factor=0.4, geom=pv.Arrow())
                p.add_mesh(arrows, color="darkred", opacity=0.8)

                p.camera_position = CAM_POS
                p.render()
                if not off_screen:
                    p.show()
                p.screenshot(ct / "fibers_img.png")

        # Points (group 0,1 = blue, 2 = red line, 3 = green)
        if points_mesh is not None and "group" in points_mesh.point_data:
            with Plotter(off_screen=off_screen) as p:
                p.add_mesh(tetra, opacity=0.12, show_edges=False, color="lightgray")
                pts = points_mesh.points
                grp = points_mesh.point_data["group"]
                for g, color in ((0, "blue"), (1, "blue"), (3, "green")):
                    mask = grp == g
                    if np.any(mask):
                        p.add_mesh(
                            pv.PolyData(pts[mask]),
                            color=color,
                            point_size=12,
                            render_points_as_spheres=True,
                        )
                red_mask = grp == 2
                if np.sum(red_mask) >= 2:
                    red_pts = pts[red_mask]
                    n = len(red_pts)
                    # PyVista lines: flat array [n_pts, i, j, n_pts, i, j, ...] per segment
                    lines = np.column_stack(
                        [np.full(n - 1, 2), np.arange(n - 1), np.arange(1, n)]
                    ).ravel()
                    line_mesh = pv.PolyData(red_pts, lines=lines)
                    p.add_mesh(line_mesh, color="red", line_width=5)

                p.camera_position = CAM_POS
                p.render()
                if not off_screen:
                    p.show()
                p.screenshot(ct / "points_img.png")

        with Plotter(off_screen=off_screen) as p:
            p.add_mesh(tetra)
            p.camera_position = CAM_POS
            p.render()
            if not off_screen:
                p.show()
            p.screenshot(ct / "mesh_img.png")

Main pipeline

Discretisation

Bases: NamedTuple

Source code in gdellipsoid/main.py 
264
265
266
class Discretisation(NamedTuple):
    nu: int
    nv: int

generate_ellipsoid(out, off_screen=False, h=2.0)

Source code in gdellipsoid/main.py 
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
def generate_ellipsoid(out, off_screen=False, h: float = 2.0):
    meshes = create_data(off_screen=off_screen, h=h)

    with gd.Container(out) as ct:
        # Convert integer arrays to float before saving to VTK for compatibility
        # SOFA's VTK reader doesn't support integer types
        meshes["regular_surf"] = _convert_int_to_float(meshes["regular_surf"])
        meshes["tetra"] = _convert_int_to_float(meshes["tetra"])
        meshes["points"] = _convert_int_to_float(meshes["points"])

        meshes["regular_surf"].save(ct / "regular_surf.vtk")
        meshes["tetra"].save(ct / "tetra.vtk")
        meshes["tetra"].save_sofa(ct / "tetra_sofa.vtk")
        meshes["points"].save(ct / "points.vtk")
        pv.PolyData(meshes["points"].points).save(ct / "points_obj.obj")
        pv.PolyData(meshes["rings"].points).save(ct / "rings.obj")
        # Force PolyData format for endo and epi (surface meshes) while keeping .vtk extension
        _mesh_to_polydata(meshes["endo"]).save(ct / "endo.vtk")
        _mesh_to_polydata(meshes["epi"]).save(ct / "epi.vtk")

        # Create clipped_endo (clipped at z=0, keep z<0)
        clipped_endo = _clip_mesh_at_z(meshes["endo"], z_value=0.0, keep_below=True)
        clipped_endo = _convert_int_to_float(clipped_endo)
        _mesh_to_polydata(clipped_endo).save(ct / "clipped_endo.vtk")

create_data(off_screen=False, h=2.0)

Source code in gdellipsoid/main.py 
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
def create_data(off_screen: bool = False, h: float = 2.0):
    n = Discretisation(31, 31)
    # Ellipsoid dimensions (kept unchanged)
    m1 = ell([7, 7, 17], n)
    m2 = ell([10, 10, 20], n)
    mf = merge(m1, m2)

    mmg = MmgOptions(hmin=1.0, hmax=2.0, hausd=0.01, hgrad=1)
    m1 = m1.remesh_surface(mmg_args=mmg)
    m2 = m2.remesh_surface(mmg_args=mmg)

    mt = _tetra_from_structured_shell(n, nr=13)
    mt = _remesh_tetra_volume(mt, h=h)
    mt = _label_tetra_points(mt, m1, m2)

    layers, labels = [], defaultdict(list)
    nu = 10
    for r in (0.1, 0.5, 0.9):
        for i, _iu in enumerate(range(nu)):
            u1 = -np.pi
            u2 = -np.arccos(5 / (17 + 3 * r))
            u = u1 + (u2 - u1) / nu * (i + 1) * 0.95
            layers.append(get_pts_r(r, u, 0))
            labels["r"].append(r)
            labels["u"].append(u)
            labels["v"].append(0)
            labels["group"].append(0)

            v = np.pi / 10
            layers.append(get_pts_r(r, u, v))
            labels["r"].append(r)
            labels["u"].append(u)
            labels["v"].append(v)
            labels["group"].append(1)

    r, v = 0.5, 0
    for u in np.linspace(-np.pi, -np.arccos(5 / (17 + 3 * r)), 30):
        layers.append(get_pts_r(r, u, v))
        labels["r"].append(r)
        labels["u"].append(u)
        labels["v"].append(v)
        labels["group"].append(2)

    u = -np.pi
    for r in (0, 1):
        layers.append(get_pts_r(r, u, v))
        labels["r"].append(r)
        labels["u"].append(u)
        labels["v"].append(v)
        labels["group"].append(3)

    mfp_pts = np.vstack(layers)
    mfp = _point_mesh(
        mfp_pts,
        point_data={
            "r": np.asarray(labels["r"], dtype=float),
            "u": np.asarray(labels["u"], dtype=float),
            "v": np.asarray(labels["v"], dtype=float),
            "group": np.asarray(labels["group"], dtype=np.int8),
        },
    )

    ring_mask = (mf.points[:, 2] >= 4.95) & (mf.points[:, 2] <= 10.0)
    mrings = _point_mesh(mf.points[ring_mask])

    # Only plot if not in offscreen mode (plotting causes segfaults in headless CI)
    if not off_screen:
        with Plotter(off_screen=False) as p:
            p.add_mesh(mt, scalars=Labels.ENDOS_EPI)

    return {
        "regular_surf": mf,
        "points": mfp,
        "tetra": mt,
        "rings": mrings,
        "endo": m1,
        "epi": m2,
    }

Fibers

create_fibers(out, angles=(80, -80))

Source code in gdellipsoid/fibers.py 
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
def create_fibers(out, angles=(80, -80)):
    with gd.Container(out) as ct:
        m = load_mesh(ct.tetra)

    n_points = m.points.shape[0]
    m.point_data["labels"] = np.zeros(n_points, dtype=np.float64)
    endo_epi = m.point_data[Labels.ENDOS_EPI]
    m.point_data["labels"][np.where(endo_epi == Labels.Endos_Epi.ENDO)] = 1
    m.point_data["labels"][np.where(endo_epi == Labels.Endos_Epi.EPI)] = 2

    with tempfile.NamedTemporaryFile(suffix=".vtk", delete=False) as f:
        try:
            m.save(f.name)
            res = get_transmural_info(f.name)
        finally:
            os.unlink(f.name)

    pts = res.points
    gradr = res.point_data["gradr"]
    gradr[np.where(np.isclose(pts[:, 2], pts[:, 2].max()))] = 0

    angles_rad = np.deg2rad(angles)
    fiber, sheet = [], []
    r_arr = res.point_data["r"]
    for i in range(res.points.shape[0]):
        r = r_arr[i]
        ang = angles_rad[0] * (1 - r) + r * angles_rad[1]
        r_ = N(gradr[i])
        t_ = -N(np.cross(r_, [0, 0, 1]))
        rot = R.from_rotvec(ang * r_)
        t_ = rot.apply(t_)
        fiber.append(t_)
        sheet.append(r_)

    res.point_data["fiber"] = np.array(fiber, dtype=float)
    res.point_data["sheet"] = np.array(sheet, dtype=float)

    _sync_point_data_to_pv(res)
    pv_centers = m.pv_mesh.cell_centers()
    pv_centers = pv_centers.sample(res.pv_mesh)

    n_cells = pv_centers.n_points
    vertices = np.arange(n_cells, dtype=np.int64).reshape(-1, 1)
    mcenters = Mesh(
        points=np.asarray(pv_centers.points),
        cells={"vertex": vertices},
        point_data={k: np.asarray(v) for k, v in pv_centers.point_data.items()},
    )
    with gd.Container(out) as ct:
        pv_centers.save(ct / "mesh_cell_fibers.vtk")
        fibs_np = pv_centers.point_data["fiber"]
        Mesh(fibs_np[:, :3], {}).save(ct / "fibs_tetra.obj")


    return res, mcenters

Purkinje network

create_purkinje(out, show_plot=True)

Source code in gdellipsoid/purkinje.py 
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
def create_purkinje(out, show_plot=True):
    with gd.Container(out) as ct:
        endo_path = ct / "endo.vtk"
        msh = _polydata_to_meshio(endo_path)
        print(msh)
        tri_cells = next(c for c in msh.cells if c.type == "triangle")
        n_points = msh.points.shape[0]
        init_idx = min(1308, n_points - 1)
        second_idx = min(3327, n_points - 1)
        if second_idx == init_idx:
            second_idx = (init_idx + 1) % n_points
        mesh = FMesh(
            verts=msh.points,
            connectivity=tri_cells.data,
            init_node=msh.points[init_idx],
        )
        param = FractalTreeParameters(
            filename=out / "purkinje_lv",
            second_node=msh.points[second_idx],
            N_it=10,
            init_length=15,
            length=10,
            l_segment=5,
            branch_angle=np.deg2rad(30),
            repulsitivity=0.3,
            generate_fascicles=True,
            fascicles_angles=(np.deg2rad(0),),
            fascicles_length=(5,),
            save=True,
            save_paraview=True,
        )
        res: FractalTreeResult = generate_fractal_tree(mesh, param)
        print(res)
        plot_purk(mesh, res, show=show_plot)
        tree_mesh = Mesh(
            points=res.nodes,
            cells={"line": np.array(res.lines, dtype=np.int64)},
        )
        tree_mesh.save(ct / "purkinje_mesh.vtk")

plot_purk(mesh, res, show=True)

Source code in gdellipsoid/purkinje.py 
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
def plot_purk(mesh, res, show=True):
    # Skip plotting when show=False to avoid segfaults in headless CI
    if not show:
        return
    surface = Mesh(
        points=np.asarray(mesh.verts),
        cells={"triangle": np.asarray(mesh.connectivity, dtype=np.int64)},
    )
    tree = Mesh(
        points=res.nodes,
        cells={"line": np.array(res.lines, dtype=np.int64)},
    )
    with Plotter(off_screen=False) as p:
        p.add_mesh(surface, opacity=0.4, show_edges=True)
        p.add_mesh(tree, color="red", line_width=2)

Plotting

Plot mesh, surfaces, Purkinje network, and fibers in a single view.

plot_all(out, off_screen=False)

Plot in one image: mesh, surfaces, Purkinje, and fibers.

Source code in gdellipsoid/plot.py 
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
def plot_all(out, off_screen=False):
    """Plot in one image: mesh, surfaces, Purkinje, and fibers."""
    with gd.Container(out) as ct:
        # Required: tetra and at least one surface
        tetra = Mesh.load(ct / "tetra.vtk")
        regular_surf = Mesh.load(ct / "regular_surf.vtk")

        # Optional: endo, epi, points, fibers, purkinje
        endo = _load_optional(ct, "endo.vtk")
        epi = _load_optional(ct, "epi.vtk")
        points_mesh = _load_optional(ct, "points.vtk")
        purkinje = _load_optional(ct, "purkinje_mesh.vtk")
        # mesh_cell_fibers.vtk is PolyData (point cloud); load with PyVista
        mesh_cell_fibers_pv = None
        if Path(ct / "mesh_cell_fibers.vtk").exists():
            mesh_cell_fibers_pv = pv.read(str(ct / "mesh_cell_fibers.vtk"))

        with Plotter(off_screen=off_screen) as p:
            # Volume mesh: transparent or wireframe so interior is visible
            p.add_mesh(tetra, opacity=0.12, show_edges=False, color="lightgray")

            # Surfaces: combined regular surface + optional endo/epi
            p.add_mesh(regular_surf, opacity=0.5, show_edges=False, color="white")

            if endo is not None:
                p.add_mesh(endo, opacity=0.6, show_edges=False, color="coral")
            if epi is not None:
                p.add_mesh(epi, opacity=0.6, show_edges=False, color="lightblue")

            p.camera_position = CAM_POS
            p.render()
            if not off_screen:
                p.show()
            p.screenshot(ct / "surfaces_img.png")

        # Purkinje network
        if purkinje is not None:
            with Plotter(off_screen=off_screen) as p:
                p.add_mesh(endo)
                p.add_mesh(purkinje, color="red", line_width=3)

                p.camera_position = [
                    (0.6051986645911922, -0.46423570086746, 37.204488072156266),
                    (
                        -0.0003440849696247916,
                        -5.655034101526013e-05,
                        -5.999310324980369,
                    ),
                    (0.004839410606533647, -0.9999298485328891, -0.010811018359763551),
                ]
                p.render()
                if not off_screen:
                    p.show()
                p.screenshot(ct / "purkinje_img.png")

        # Fibers: vectors at cell centers (subsample for clarity)
        if (
            mesh_cell_fibers_pv is not None
            and "fiber" in mesh_cell_fibers_pv.point_data
        ):
            with Plotter(off_screen=off_screen) as p:
                p.add_mesh(tetra, opacity=0.12, show_edges=False, color="lightgray")
                if endo is not None:
                    p.add_mesh(endo, opacity=0.6, show_edges=False, color="coral")
                if epi is not None:
                    p.add_mesh(epi, opacity=0.6, show_edges=False, color="lightblue")

                pf = mesh_cell_fibers_pv
                n_pts = pf.n_points
                step = max(1, n_pts // 800)
                subsampled = pf.copy()
                subsampled.points = np.asarray(pf.points)[::step]
                for key in ("fiber", "sheet"):
                    if key in pf.point_data.keys():
                        subsampled.point_data[key] = np.asarray(pf.point_data[key])[
                            ::step
                        ]
                subsampled.set_active_vectors("fiber")
                arrows = subsampled.glyph(orient="fiber", factor=0.4, geom=pv.Arrow())
                p.add_mesh(arrows, color="darkred", opacity=0.8)

                p.camera_position = CAM_POS
                p.render()
                if not off_screen:
                    p.show()
                p.screenshot(ct / "fibers_img.png")

        # Points (group 0,1 = blue, 2 = red line, 3 = green)
        if points_mesh is not None and "group" in points_mesh.point_data:
            with Plotter(off_screen=off_screen) as p:
                p.add_mesh(tetra, opacity=0.12, show_edges=False, color="lightgray")
                pts = points_mesh.points
                grp = points_mesh.point_data["group"]
                for g, color in ((0, "blue"), (1, "blue"), (3, "green")):
                    mask = grp == g
                    if np.any(mask):
                        p.add_mesh(
                            pv.PolyData(pts[mask]),
                            color=color,
                            point_size=12,
                            render_points_as_spheres=True,
                        )
                red_mask = grp == 2
                if np.sum(red_mask) >= 2:
                    red_pts = pts[red_mask]
                    n = len(red_pts)
                    # PyVista lines: flat array [n_pts, i, j, n_pts, i, j, ...] per segment
                    lines = np.column_stack(
                        [np.full(n - 1, 2), np.arange(n - 1), np.arange(1, n)]
                    ).ravel()
                    line_mesh = pv.PolyData(red_pts, lines=lines)
                    p.add_mesh(line_mesh, color="red", line_width=5)

                p.camera_position = CAM_POS
                p.render()
                if not off_screen:
                    p.show()
                p.screenshot(ct / "points_img.png")

        with Plotter(off_screen=off_screen) as p:
            p.add_mesh(tetra)
            p.camera_position = CAM_POS
            p.render()
            if not off_screen:
                p.show()
            p.screenshot(ct / "mesh_img.png")

Utilities

Labels

Source code in gdellipsoid/utils.py 
 4
 5
 6
 7
 8
 9
10
11
class Labels:
    ENDOS_EPI = "Endos_Epi"

    class Endos_Epi:
        LV = 1
        ENDO = LV
        EPI = 2
        MYO = 3

N(x)

Source code in gdellipsoid/utils.py 
14
15
16
17
def N(x):
    x = np.asarray(x, float)
    n = np.linalg.norm(x, axis=-1, keepdims=True)
    return x / np.where(n == 0, 1, n)