From b010f79c0701fbcee7b5bd8c552b05acb1f91f84 Mon Sep 17 00:00:00 2001 From: Javiera Jilberto Vallejos Date: Fri, 5 Dec 2025 11:29:56 -0800 Subject: [PATCH 1/5] adding fiber generation codes --- utilities/fiber_generation/README.md | 35 + .../fiber_generation/example/ot/doste.png | 3 + .../example/ot/mesh-complete.mesh.vtu | 3 + .../example/ot/mesh-surfaces/EPI_APEX.vtp | 3 + .../example/ot/mesh-surfaces/av.vtp | 3 + .../example/ot/mesh-surfaces/endo_lv.vtp | 3 + .../example/ot/mesh-surfaces/endo_rv.vtp | 3 + .../example/ot/mesh-surfaces/epi.vtp | 3 + .../example/ot/mesh-surfaces/mv.vtp | 3 + .../example/ot/mesh-surfaces/pv.vtp | 3 + .../example/ot/mesh-surfaces/top.vtp | 3 + .../example/ot/mesh-surfaces/tv.vtp | 3 + .../example/truncated/VOLUME.vtu | 3 + .../example/truncated/bayer.png | 3 + .../example/truncated/mesh-surfaces/BASE.vtp | 3 + .../example/truncated/mesh-surfaces/EPI.vtp | 3 + .../truncated/mesh-surfaces/EPI_APEX.vtp | 3 + .../truncated/mesh-surfaces/EPI_APEX2.vtp | 3 + .../truncated/mesh-surfaces/EPI_MID.vtp | 3 + .../example/truncated/mesh-surfaces/LV.vtp | 3 + .../example/truncated/mesh-surfaces/RV.vtp | 3 + .../truncated/mesh-surfaces/exterior.vtp | 3 + utilities/fiber_generation/main_bayer.py | 64 + utilities/fiber_generation/main_doste.py | 85 ++ utilities/fiber_generation/src/FibGen.py | 1084 +++++++++++++++++ .../src/templates/solver_bayer.xml | 266 ++++ .../src/templates/solver_doste.xml | 492 ++++++++ 27 files changed, 2089 insertions(+) create mode 100644 utilities/fiber_generation/README.md create mode 100644 utilities/fiber_generation/example/ot/doste.png create mode 100644 utilities/fiber_generation/example/ot/mesh-complete.mesh.vtu create mode 100644 utilities/fiber_generation/example/ot/mesh-surfaces/EPI_APEX.vtp create mode 100644 utilities/fiber_generation/example/ot/mesh-surfaces/av.vtp create mode 100644 utilities/fiber_generation/example/ot/mesh-surfaces/endo_lv.vtp create mode 100644 utilities/fiber_generation/example/ot/mesh-surfaces/endo_rv.vtp create mode 100644 utilities/fiber_generation/example/ot/mesh-surfaces/epi.vtp create mode 100644 utilities/fiber_generation/example/ot/mesh-surfaces/mv.vtp create mode 100644 utilities/fiber_generation/example/ot/mesh-surfaces/pv.vtp create mode 100644 utilities/fiber_generation/example/ot/mesh-surfaces/top.vtp create mode 100644 utilities/fiber_generation/example/ot/mesh-surfaces/tv.vtp create mode 100644 utilities/fiber_generation/example/truncated/VOLUME.vtu create mode 100644 utilities/fiber_generation/example/truncated/bayer.png create mode 100644 utilities/fiber_generation/example/truncated/mesh-surfaces/BASE.vtp create mode 100644 utilities/fiber_generation/example/truncated/mesh-surfaces/EPI.vtp create mode 100644 utilities/fiber_generation/example/truncated/mesh-surfaces/EPI_APEX.vtp create mode 100644 utilities/fiber_generation/example/truncated/mesh-surfaces/EPI_APEX2.vtp create mode 100644 utilities/fiber_generation/example/truncated/mesh-surfaces/EPI_MID.vtp create mode 100644 utilities/fiber_generation/example/truncated/mesh-surfaces/LV.vtp create mode 100644 utilities/fiber_generation/example/truncated/mesh-surfaces/RV.vtp create mode 100644 utilities/fiber_generation/example/truncated/mesh-surfaces/exterior.vtp create mode 100644 utilities/fiber_generation/main_bayer.py create mode 100644 utilities/fiber_generation/main_doste.py create mode 100644 utilities/fiber_generation/src/FibGen.py create mode 100644 utilities/fiber_generation/src/templates/solver_bayer.xml create mode 100644 utilities/fiber_generation/src/templates/solver_doste.xml diff --git a/utilities/fiber_generation/README.md b/utilities/fiber_generation/README.md new file mode 100644 index 000000000..41379c364 --- /dev/null +++ b/utilities/fiber_generation/README.md @@ -0,0 +1,35 @@ +# SV-fibergen +Python + SVmultiphysics codes for fiber generation. Two methods are implemented: +* Bayer et al. (2012). [link](https://doi.org/10.1007/s10439-012-0593-5) +* Doste et al. (2018). [link](https://doi.org/10.1002/cnm.3185) + +### Examples +The `main_bayer.py` and `main_doste.py` are scripts to run both methods in the geometry described in the `example/truncated` and `example/ot` folders respectively. + +Results for truncated BiV (Bayer) +Results for BiV w/ outflow tracts (Doste) + +Note that the Doste methods needs a geometry with outflow tracts to be run (each valve needs to be defined as a separated surface). Bayer can be run in any biventricular geometry. + + +### Updates to the old code +* All operations are vectorized now. +* The SVmultiphysics solver now solves a Laplace equation. +* In Bayer: For the bislerp interpolation, instead of using the correction described in Bayer et al. (that returns a discontinuity), the basis are flipped to maintain a coherent fiber direction (see function `generate_fibers_BiV_Bayer_cells` in `FibGen.py`). +* In Bayer: The beta angles were not being included correctly. The second rotation was being applied respect the first vector (circumferential) when it should be respect to the second vector (longitudinal) (see function `generate_fibers_BiV_Bayer_cells` in `FibGen.py`). + +### Notes on SVmultiphysics solver + +To solve a Laplace equation directly from the transient HEAT solver in SVmultiphysics, in `` we need to set, +``` + 1 + 1 + 0. +``` +and in ``, +``` + 1.0 + 0.0 + 0.0 +``` +This will allow us to solve the Laplace equation directly in 1 timestep and 1 iteration. \ No newline at end of file diff --git a/utilities/fiber_generation/example/ot/doste.png b/utilities/fiber_generation/example/ot/doste.png new file mode 100644 index 000000000..1ea35dfb1 --- /dev/null +++ b/utilities/fiber_generation/example/ot/doste.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:4da5c28acd2e2b9745ba972ec783b7341afa2156a48244059684ca28f45652ac +size 494459 diff --git a/utilities/fiber_generation/example/ot/mesh-complete.mesh.vtu b/utilities/fiber_generation/example/ot/mesh-complete.mesh.vtu new file mode 100644 index 000000000..0afc681f0 --- /dev/null +++ b/utilities/fiber_generation/example/ot/mesh-complete.mesh.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:9347d6e7907816749cc16471130397f264c7721d029626fb1b9f7fdd6ac1e02c +size 28209696 diff --git a/utilities/fiber_generation/example/ot/mesh-surfaces/EPI_APEX.vtp b/utilities/fiber_generation/example/ot/mesh-surfaces/EPI_APEX.vtp new file mode 100644 index 000000000..0dd195548 --- /dev/null +++ b/utilities/fiber_generation/example/ot/mesh-surfaces/EPI_APEX.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1431a90e42af727af648dacde59095e9c3fd8e37d4a7dfef3fef57b633eaf139 +size 2867 diff --git a/utilities/fiber_generation/example/ot/mesh-surfaces/av.vtp b/utilities/fiber_generation/example/ot/mesh-surfaces/av.vtp new file mode 100644 index 000000000..4bbb83ab5 --- /dev/null +++ b/utilities/fiber_generation/example/ot/mesh-surfaces/av.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:ee7d4d8eb31b2639f75599c188cb2af9fa4fcc466d27e34b65dc3ec4f4f3633b +size 14489 diff --git a/utilities/fiber_generation/example/ot/mesh-surfaces/endo_lv.vtp b/utilities/fiber_generation/example/ot/mesh-surfaces/endo_lv.vtp new file mode 100644 index 000000000..5b2740a69 --- /dev/null +++ b/utilities/fiber_generation/example/ot/mesh-surfaces/endo_lv.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:fa09ea16a80da5813b133a22f42b94b7aae48cd414e9c8ce174d46a4d7eb3a19 +size 563261 diff --git a/utilities/fiber_generation/example/ot/mesh-surfaces/endo_rv.vtp b/utilities/fiber_generation/example/ot/mesh-surfaces/endo_rv.vtp new file mode 100644 index 000000000..8aef06c70 --- /dev/null +++ b/utilities/fiber_generation/example/ot/mesh-surfaces/endo_rv.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:b0c07b7ca1ad1b8c38d826d363f5572e1e6359595f5701a2b3f8f25606e30073 +size 734382 diff --git a/utilities/fiber_generation/example/ot/mesh-surfaces/epi.vtp b/utilities/fiber_generation/example/ot/mesh-surfaces/epi.vtp new file mode 100644 index 000000000..8336320c8 --- /dev/null +++ b/utilities/fiber_generation/example/ot/mesh-surfaces/epi.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:2168042f6f99283d747e5895df8bb30d443631b5be7e31b51c55e1a7be83fa29 +size 1293881 diff --git a/utilities/fiber_generation/example/ot/mesh-surfaces/mv.vtp b/utilities/fiber_generation/example/ot/mesh-surfaces/mv.vtp new file mode 100644 index 000000000..4b9b8492b --- /dev/null +++ b/utilities/fiber_generation/example/ot/mesh-surfaces/mv.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:712904ededfdd630a736abfc213758155f8dd68a7f1512c2128c2b0305eeecb6 +size 90751 diff --git a/utilities/fiber_generation/example/ot/mesh-surfaces/pv.vtp b/utilities/fiber_generation/example/ot/mesh-surfaces/pv.vtp new file mode 100644 index 000000000..ba4ef7861 --- /dev/null +++ b/utilities/fiber_generation/example/ot/mesh-surfaces/pv.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1b99698dc66260cd1c63d4612402f07e98b0a547549e42d896048bf7bdbe2afb +size 12559 diff --git a/utilities/fiber_generation/example/ot/mesh-surfaces/top.vtp b/utilities/fiber_generation/example/ot/mesh-surfaces/top.vtp new file mode 100644 index 000000000..28e22860f --- /dev/null +++ b/utilities/fiber_generation/example/ot/mesh-surfaces/top.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:edec0ae0c98ef7db7add744db0340aeba1c925613403ec4bc1e580d0c1454767 +size 87143 diff --git a/utilities/fiber_generation/example/ot/mesh-surfaces/tv.vtp b/utilities/fiber_generation/example/ot/mesh-surfaces/tv.vtp new file mode 100644 index 000000000..ebf70fd05 --- /dev/null +++ b/utilities/fiber_generation/example/ot/mesh-surfaces/tv.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:b33291fd2b529845d68d4584c2f4f5241c288dda1eb6a5e76f5bc839e21e2f62 +size 19198 diff --git a/utilities/fiber_generation/example/truncated/VOLUME.vtu b/utilities/fiber_generation/example/truncated/VOLUME.vtu new file mode 100644 index 000000000..bec495b6e --- /dev/null +++ b/utilities/fiber_generation/example/truncated/VOLUME.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:bcd011125435852298177578498ebcf2b5aac1080b824163d1facdefc48c2a27 +size 24524913 diff --git a/utilities/fiber_generation/example/truncated/bayer.png b/utilities/fiber_generation/example/truncated/bayer.png new file mode 100644 index 000000000..73a017c6f --- /dev/null +++ b/utilities/fiber_generation/example/truncated/bayer.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c01ec2a92c75b6047fa9b6602924a619d1e20b577cc26de7852a91cc06d19b6a +size 313768 diff --git a/utilities/fiber_generation/example/truncated/mesh-surfaces/BASE.vtp b/utilities/fiber_generation/example/truncated/mesh-surfaces/BASE.vtp new file mode 100644 index 000000000..6268fb43c --- /dev/null +++ b/utilities/fiber_generation/example/truncated/mesh-surfaces/BASE.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:b56eb592388c9cc2226c6372ffc8e0199952698f84ffa81330672d272503d7aa +size 191256 diff --git a/utilities/fiber_generation/example/truncated/mesh-surfaces/EPI.vtp b/utilities/fiber_generation/example/truncated/mesh-surfaces/EPI.vtp new file mode 100644 index 000000000..1aa8cca64 --- /dev/null +++ b/utilities/fiber_generation/example/truncated/mesh-surfaces/EPI.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:2c5ecc05750d086d04dfe9902e7a7e200831ed5cb11ffe6339186112d8bac0ce +size 526632 diff --git a/utilities/fiber_generation/example/truncated/mesh-surfaces/EPI_APEX.vtp b/utilities/fiber_generation/example/truncated/mesh-surfaces/EPI_APEX.vtp new file mode 100644 index 000000000..f713c97d0 --- /dev/null +++ b/utilities/fiber_generation/example/truncated/mesh-surfaces/EPI_APEX.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c6e68548080a65b62e5ea62af5d685cafa9181604f43ae45fa9680b0f6e09938 +size 2833 diff --git a/utilities/fiber_generation/example/truncated/mesh-surfaces/EPI_APEX2.vtp b/utilities/fiber_generation/example/truncated/mesh-surfaces/EPI_APEX2.vtp new file mode 100644 index 000000000..f713c97d0 --- /dev/null +++ b/utilities/fiber_generation/example/truncated/mesh-surfaces/EPI_APEX2.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c6e68548080a65b62e5ea62af5d685cafa9181604f43ae45fa9680b0f6e09938 +size 2833 diff --git a/utilities/fiber_generation/example/truncated/mesh-surfaces/EPI_MID.vtp b/utilities/fiber_generation/example/truncated/mesh-surfaces/EPI_MID.vtp new file mode 100644 index 000000000..f0b984ad4 --- /dev/null +++ b/utilities/fiber_generation/example/truncated/mesh-surfaces/EPI_MID.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:cc6a73358eff58939c32e73dfd9326b9a53cf25afa409e5c27b7a87ef3ec08ad +size 525056 diff --git a/utilities/fiber_generation/example/truncated/mesh-surfaces/LV.vtp b/utilities/fiber_generation/example/truncated/mesh-surfaces/LV.vtp new file mode 100644 index 000000000..9e9dcbb00 --- /dev/null +++ b/utilities/fiber_generation/example/truncated/mesh-surfaces/LV.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:e16a121071598c4696ea25623b52595c349bca27823c6d9cfb6e1bdbc7b77767 +size 145544 diff --git a/utilities/fiber_generation/example/truncated/mesh-surfaces/RV.vtp b/utilities/fiber_generation/example/truncated/mesh-surfaces/RV.vtp new file mode 100644 index 000000000..69f68ba06 --- /dev/null +++ b/utilities/fiber_generation/example/truncated/mesh-surfaces/RV.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c6028ab47b1af4947b9ee47cd79dfc58105d8e2bf86f1167aed877609023935c +size 128852 diff --git a/utilities/fiber_generation/example/truncated/mesh-surfaces/exterior.vtp b/utilities/fiber_generation/example/truncated/mesh-surfaces/exterior.vtp new file mode 100644 index 000000000..4f3247355 --- /dev/null +++ b/utilities/fiber_generation/example/truncated/mesh-surfaces/exterior.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:889f15170752b6310d35a52f9e1b207ab84afa1e3ce7c4be24c114403015fe3f +size 978206 diff --git a/utilities/fiber_generation/main_bayer.py b/utilities/fiber_generation/main_bayer.py new file mode 100644 index 000000000..490b0677c --- /dev/null +++ b/utilities/fiber_generation/main_bayer.py @@ -0,0 +1,64 @@ +#!/usr/bin/env python +# -*-coding:utf-8 -*- +''' +Created on 2025/11/21 20:38:14 + +@author: Javiera Jilberto Vallejos +''' + +import os +import src.FibGen as fg +from time import time + +########################################################### +############ USER INPUTS ################################ +########################################################### + +run_flag = True +svfsi_exec = "svmultiphysics " + +mesh_path = "example/truncated/VOLUME.vtu" +surfaces_dir = f"example/truncated/mesh-surfaces" +outdir = "example/truncated/output_b" + +surface_names = {'epi': 'EPI.vtp', + 'epi_apex': 'EPI_APEX.vtp', # New surface + 'base': 'BASE.vtp', + 'endo_lv': 'LV.vtp', + 'endo_rv': 'RV.vtp'} + +# Parameters for the Bayer et al. method https://doi.org/10.1007/s10439-012-0593-5. +params = { + "ALFA_END": 60.0, + "ALFA_EPI": -60.0, + "BETA_END": 20.0, + "BETA_EPI": -20.0, +} + + +########################################################### +############ FIBER GENERATION ########################### +########################################################### + +# Make sure the paths are full paths +mesh_path = os.path.abspath(mesh_path) +surfaces_dir = os.path.abspath(surfaces_dir) +outdir = os.path.abspath(outdir) + +start = time() +fg.generate_epi_apex(mesh_path, surfaces_dir, surface_names) + +# Run the Laplace solver +if run_flag: + template_file = os.path.join(os.path.dirname(os.path.abspath(__file__)), "src", "templates", "solver_bayer.xml") + laplace_results_file = fg.runLaplaceSolver(mesh_path, surfaces_dir, mesh_path, svfsi_exec, template_file, outdir, surface_names) +laplace_results_file = outdir + '/result_001.vtu' + +# Generate the fiber directions +result_mesh = fg.generate_fibers_BiV_Bayer_cells(outdir, laplace_results_file, params, return_angles=True, return_intermediate=True) + +print(f"generate fibers (Bayer method) elapsed time: {time() - start:.3f} s") + +# Optional, save the result mesh with intermediate field and angles for checking +result_mesh_path = os.path.join(outdir, "results_bayer.vtu") +result_mesh.save(result_mesh_path) diff --git a/utilities/fiber_generation/main_doste.py b/utilities/fiber_generation/main_doste.py new file mode 100644 index 000000000..6a612274f --- /dev/null +++ b/utilities/fiber_generation/main_doste.py @@ -0,0 +1,85 @@ +#!/usr/bin/env python +# -*-coding:utf-8 -*- +''' +Created on 2025/11/21 20:38:14 + +@author: Javiera Jilberto Vallejos +''' + +import os +import src.FibGen as fg +from time import time + +########################################################### +############ USER INPUTS ################################ +########################################################### + +run_flag = True +method = 'doste' +svfsi_exec = "svmultiphysics " + +mesh_path = "example/ot/mesh-complete.mesh.vtu" +surfaces_dir = f"example/ot/mesh-surfaces" +outdir = "example/ot/output_d" + +surface_names = {'epi': 'epi.vtp', + 'epi_apex': 'epi_apex.vtp', # New surface + 'av': 'av.vtp', + 'mv': 'mv.vtp', + 'tv': 'tv.vtp', + 'pv': 'pv.vtp', + 'base': 'top.vtp', # This is all the valves together, it is used to find the apex. + 'endo_lv': 'endo_lv.vtp', + 'endo_rv': 'endo_rv.vtp'} + +# Parameters from the Doste paper https://doi.org/10.1002/cnm.3185 +params = { + # A = alpha angle + 'AENDORV' : 90, + 'AEPIRV' : -25, + 'AENDOLV' : 60, + 'AEPILV' : -60, + + 'AOTENDOLV' : 90, + 'AOTENDORV' : 90, + 'AOTEPILV' : 0, + 'AOTEPIRV' : 0, + + # B = beta angle (this have an opposite sign to the Doste paper, + # but it's because the longitudinal direction is opposite) + 'BENDORV' : 20, + 'BEPIRV' : -20, + 'BENDOLV' : 20, + 'BEPILV' : -20, +} + + +########################################################### +############ FIBER GENERATION ########################### +########################################################### + +# Make sure the paths are full paths +mesh_path = os.path.abspath(mesh_path) +surfaces_dir = os.path.abspath(surfaces_dir) +outdir = os.path.abspath(outdir) + +# Generate the apex surface +start = time() + +start = time() +fg.generate_epi_apex(mesh_path, surfaces_dir, surface_names) + +# Run the Laplace solver +if run_flag: + template_file = os.path.join(os.path.dirname(os.path.abspath(__file__)), "src", "templates", "solver_doste.xml") + laplace_results_file = fg.runLaplaceSolver(mesh_path, surfaces_dir, mesh_path, svfsi_exec, template_file, outdir, surface_names) +laplace_results_file = outdir + '/result_001.vtu' + +# Generate the fiber directions +result_mesh = fg.generate_fibers_BiV_Doste_cells(outdir, laplace_results_file, params, return_angles=True, return_intermediate=False) + +print(f"generate fibers (Doste method) elapsed time: {time() - start:.3f} s") + +# Optional, save the result mesh with intermediate field and angles for checking +result_mesh_path = os.path.join(outdir, "results_doste.vtu") +result_mesh.save(result_mesh_path) diff --git a/utilities/fiber_generation/src/FibGen.py b/utilities/fiber_generation/src/FibGen.py new file mode 100644 index 000000000..3e62c19fc --- /dev/null +++ b/utilities/fiber_generation/src/FibGen.py @@ -0,0 +1,1084 @@ +#!/usr/bin/env python +# -*-coding:utf-8 -*- +''' +Created on 2025/11/21 20:43:23 + +@author: Javiera Jilberto Vallejos +''' + +import os +import re +import numpy as np +import pyvista as pv +import time +import copy + + +def normalize(x): + """ + Normalize each row of an (N, 3) array. Zero rows remain zero. + + Args: + x: array-like of shape (N, 3) + Returns: + np.ndarray of shape (N, 3) with row-wise normalized vectors. + """ + a = np.asarray(x, dtype=float) + if a.ndim != 2 or a.shape[1] != 3: + raise ValueError("normalize expects an array of shape (N, 3)") + norms = np.linalg.norm(a, axis=1, keepdims=True) + safe_norms = np.where(norms == 0.0, 1.0, norms) + out = a / safe_norms + zero_rows = (norms.squeeze() == 0.0) + if np.any(zero_rows): + out[zero_rows] = 0.0 + return out + + + +def get_normal_plane_svd(points): # Find the plane that minimizes the distance given N points + centroid = np.mean(points, axis=0) + svd = np.linalg.svd(points - centroid) + normal = svd[2][-1] + normal = normal/np.linalg.norm(normal) + return normal, centroid + + +def generate_epi_apex(mesh_path, surfaces_dir, surface_names): + ''' + Generate the epi apex and epi mid surfaces from the epi surface of the BiV. + + Parameters: + ----------- + surfaces_dir : str + Directory containing surface meshes + surface_names : list of str + List of surface mesh filenames + ''' + + # Load the epi surface + epi_name = os.path.join(surfaces_dir, surface_names['epi']) + epi_mesh = pv.read(epi_name) + epi_points = epi_mesh.points + epi_cells = epi_mesh.faces + epi_eNoN = epi_cells[0] + epi_cells = epi_cells.reshape((-1, epi_eNoN + 1)) + epi_cells = epi_cells[:, 1:] + epi_global_node_id = epi_mesh.point_data['GlobalNodeID'] + epi_global_cell_id = epi_mesh.cell_data['GlobalElementID'] + + # Load the base surface + base_name = os.path.join(surfaces_dir, surface_names['base']) + base_mesh = pv.read(base_name) + base_global_node_id = base_mesh.point_data['GlobalNodeID'] + + # Extract the boundary of the epi surface (at the top) to find the apex point + epi_base_global_node_id = np.intersect1d(epi_global_node_id, base_global_node_id) + epi_base_nodes = np.where(np.isin(epi_global_node_id, epi_base_global_node_id))[0] + + # # Get normal + base_normal, base_centroid = get_normal_plane_svd(epi_points[epi_base_nodes, :]) + + # Find the index of the apex point of the epi surface + distance = np.abs(base_normal@(epi_points - base_centroid).T) + epi_apex_point_index = np.argmax(distance) + + # Find elements containing the apex point + epi_apex_cell_index = np.where(epi_cells == epi_apex_point_index)[0] + + # Create epi_apex mesh + submesh_cells = epi_cells[epi_apex_cell_index] + submesh_xyz = np.zeros([len(np.unique(submesh_cells)), epi_points.shape[1]]) + map_mesh_submesh = np.ones(epi_points.shape[0], dtype=int)*-1 + map_submesh_mesh = np.zeros(submesh_xyz.shape[0], dtype=int) + child_elems_new = np.zeros(submesh_cells.shape, dtype=int) + + cont = 0 + for e in range(submesh_cells.shape[0]): + for i in range(submesh_cells.shape[1]): + if map_mesh_submesh[submesh_cells[e,i]] == -1: + child_elems_new[e,i] = cont + submesh_xyz[cont] = epi_points[submesh_cells[e,i]] + map_mesh_submesh[submesh_cells[e,i]] = cont + map_submesh_mesh[cont] = submesh_cells[e,i] + cont += 1 + else: + child_elems_new[e,i] = map_mesh_submesh[submesh_cells[e,i]] + + epi_apex_cells_type = np.full((child_elems_new.shape[0], 1), epi_eNoN) + epi_apex_cells = np.hstack((epi_apex_cells_type, child_elems_new)) + epi_apex_cells = np.hstack(epi_apex_cells) + + # Get global IDs + epi_apex_global_node_id = epi_global_node_id[map_submesh_mesh] + epi_apex_global_cell_id = epi_global_cell_id[epi_apex_cell_index] + + # Create and save mesh + epi_apex_mesh = pv.PolyData(submesh_xyz, epi_apex_cells) + epi_apex_mesh.point_data.set_array(epi_apex_global_node_id, 'GlobalNodeID') + epi_apex_mesh.cell_data.set_array(epi_apex_global_cell_id, 'GlobalElementID') + + epi_apex_name = os.path.join(surfaces_dir, surface_names['epi_apex']) + epi_apex_mesh.save(epi_apex_name) + + + +def runLaplaceSolver(mesh_dir, surfaces_dir, mesh_file, exec_svmultiphysics, template_file, outdir, surface_names): + xml_template_path = template_file + out_name = os.path.join(surfaces_dir, "../svFibers_BiV.xml") + + with open(xml_template_path, 'r') as svFile: + xml_content = svFile.read() + + # Update mesh file path using regex + mesh_pattern = r'()\s+[^\s<]+[^<]*()' + xml_content = re.sub(mesh_pattern, r'\1 ' + mesh_file + r' \2', xml_content) + + # Update face file paths - need to identify which face by checking context + # Read lines to determine context + lines = xml_content.split('\n') + updated_lines = [] + i = 0 + while i < len(lines): + line = lines[i] + + # Check if this line has a face name + face_match = re.search(r'name="([^"]+)"', line) + face_name = face_match.group(1) if face_match else None + + # Look ahead for Face_file_path + if face_name and i + 1 < len(lines) and "" in lines[i + 1]: + # Determine which file to use based on face name + if face_name == "epi": + new_path = os.path.join(surfaces_dir, surface_names['epi']) + elif face_name == "epi_top": + new_path = os.path.join(surfaces_dir, surface_names['base']) + elif face_name == "epi_apex": + new_path = os.path.join(surfaces_dir, surface_names['epi_apex']) + elif face_name == "endo_lv": + new_path = os.path.join(surfaces_dir, surface_names['endo_lv']) + elif face_name == "endo_rv": + new_path = os.path.join(surfaces_dir, surface_names['endo_rv']) + elif face_name == "mv": + new_path = os.path.join(surfaces_dir, surface_names['mv']) + elif face_name == "tv": + new_path = os.path.join(surfaces_dir, surface_names['tv']) + elif face_name == "av": + new_path = os.path.join(surfaces_dir, surface_names['av']) + elif face_name == "pv": + new_path = os.path.join(surfaces_dir, surface_names['pv']) + else: + new_path = None + + if new_path: + # Add current line + updated_lines.append(line) + # Replace the path in the next line + i += 1 + face_pattern = r'()\s+[^\s<]+[^<]*()' + updated_line = re.sub(face_pattern, r'\1 ' + new_path + r' \2', lines[i]) + updated_lines.append(updated_line) + i += 1 + continue + + # Add line as-is + updated_lines.append(line) + i += 1 + + xml_content = '\n'.join(updated_lines) + + # Update save results folder using regex + save_pattern = r'()\s+[^\s<]+[^<]*()' + xml_content = re.sub(save_pattern, r'\1 ' + outdir + r' \2', xml_content) + + with open(out_name, 'w') as svFileNew: + svFileNew.write(xml_content) + + print(" Running svMultiPhysics solver") + print(f" {exec_svmultiphysics + out_name}") + os.system(exec_svmultiphysics + out_name) + + return outdir + '/results_001.vtu' + + +def loadLaplaceSolnBayer(fileName): + ''' + Load a solution to a Laplace-Dirichlet problem from a .vtu file and extract + the solution and its gradients at the cells. + + ARGS: + fileName : str + Path to the .vtu file with the Laplace solution. The solution should be + defined at the nodes. The Laplace fields should be named as follows: + - Phi_BiV_EPI: Laplace field for the endocardium + - Phi_BiV_LV: Laplace field for the left ventricle + - Phi_BiV_RV: Laplace field for the right ventricle + - Phi_BiV_AB: Laplace field for the apex to base direction + ''' + + DATASTR1 = 'Phi_BiV_EPI' + DATASTR2 = 'Phi_BiV_LV' + DATASTR3 = 'Phi_BiV_RV' + DATASTR4 = 'Phi_BiV_AB' + + print(" Loading Laplace solution <--- %s" % (fileName)) + + # Read mesh with pyvista + result_mesh = pv.read(fileName) + + # Convert point-data to cell-data (keep point data passed to cells) + mesh_cells = result_mesh.point_data_to_cell_data() + + print(" Extracting solution and estimating gradients at cells") + + # Get cell centers (Nx3) and scalar cell arrays (N,) + cPhiEP = np.asarray(mesh_cells.cell_data[DATASTR1]) + cPhiLV = np.asarray(mesh_cells.cell_data[DATASTR2]) + cPhiRV = np.asarray(mesh_cells.cell_data[DATASTR3]) + cPhiAB = np.asarray(mesh_cells.cell_data[DATASTR4]) + + # Use pyvista's compute_derivative to get cell gradients + # compute_derivative will add arrays named '_grad' to the cell_data + grad_mesh = mesh_cells.compute_derivative(scalars=DATASTR1, gradient=True, preference='cell') + cGPhiEP = np.asarray(grad_mesh.cell_data['gradient']) + + grad_mesh = mesh_cells.compute_derivative(scalars=DATASTR2, gradient=True, preference='cell') + cGPhiLV = np.asarray(grad_mesh.cell_data['gradient']) + + grad_mesh = mesh_cells.compute_derivative(scalars=DATASTR3, gradient=True, preference='cell') + cGPhiRV = np.asarray(grad_mesh.cell_data['gradient']) + + grad_mesh = mesh_cells.compute_derivative(scalars=DATASTR4, gradient=True, preference='cell') + cGPhiAB = np.asarray(grad_mesh.cell_data['gradient']) + + # Use the mesh with cell-data (but without the large scalar arrays) as result_mesh + mesh_cells.cell_data[DATASTR1 + '_grad'] = cGPhiEP + mesh_cells.cell_data[DATASTR2 + '_grad'] = cGPhiLV + mesh_cells.cell_data[DATASTR3 + '_grad'] = cGPhiRV + mesh_cells.cell_data[DATASTR4 + '_grad'] = cGPhiAB + + return mesh_cells, cPhiEP, cPhiLV, cPhiRV, cPhiAB, \ + cGPhiEP, cGPhiLV, cGPhiRV, cGPhiAB + + +def bislerp(Q1, Q2, interp_func): + """ + Vectorized spherical interpolation between batches of rotation matrices. + Q1, Q2: (N, 3, 3) + interp_func: (N,) values in [0,1] + Returns Q: (N, 3, 3) + Notes: + - Uses wxyz quaternion convention internally, matching quat2rot below. + - Avoids per-element Python/Scipy objects for performance. + """ + def rotm_to_quat_batch(R): + # R: (N,3,3) -> q: (N,4) [w,x,y,z] + t = np.einsum('nii->n', R) # trace + q = np.zeros((R.shape[0], 4), dtype=float) + # Branch where trace is positive + mask_t = t > 0.0 + if np.any(mask_t): + S = np.sqrt(t[mask_t] + 1.0) * 2.0 + q[mask_t, 0] = 0.25 * S + q[mask_t, 1] = (R[mask_t, 2, 1] - R[mask_t, 1, 2]) / S + q[mask_t, 2] = (R[mask_t, 0, 2] - R[mask_t, 2, 0]) / S + q[mask_t, 3] = (R[mask_t, 1, 0] - R[mask_t, 0, 1]) / S + # For remaining, choose major diagonal + mask_f = ~mask_t + if np.any(mask_f): + Rf = R[mask_f] + m00 = Rf[:, 0, 0] + m11 = Rf[:, 1, 1] + m22 = Rf[:, 2, 2] + idx = np.argmax(np.stack([m00, m11, m22], axis=1), axis=1) + mf_idx = np.nonzero(mask_f)[0] + # Case idx==0 + m0 = idx == 0 + if np.any(m0): + S = np.sqrt(1.0 + Rf[m0, 0, 0] - Rf[m0, 1, 1] - Rf[m0, 2, 2]) * 2.0 + rows = mf_idx[m0] + q[rows, 0] = (Rf[m0, 2, 1] - Rf[m0, 1, 2]) / S + q[rows, 1] = 0.25 * S + q[rows, 2] = (Rf[m0, 0, 1] + Rf[m0, 1, 0]) / S + q[rows, 3] = (Rf[m0, 0, 2] + Rf[m0, 2, 0]) / S + # Case idx==1 + m1 = idx == 1 + if np.any(m1): + S = np.sqrt(1.0 + Rf[m1, 1, 1] - Rf[m1, 0, 0] - Rf[m1, 2, 2]) * 2.0 + rows = mf_idx[m1] + q[rows, 0] = (Rf[m1, 0, 2] - Rf[m1, 2, 0]) / S + q[rows, 1] = (Rf[m1, 0, 1] + Rf[m1, 1, 0]) / S + q[rows, 2] = 0.25 * S + q[rows, 3] = (Rf[m1, 1, 2] + Rf[m1, 2, 1]) / S + # Case idx==2 + m2 = idx == 2 + if np.any(m2): + S = np.sqrt(1.0 + Rf[m2, 2, 2] - Rf[m2, 0, 0] - Rf[m2, 1, 1]) * 2.0 + rows = mf_idx[m2] + q[rows, 0] = (Rf[m2, 1, 0] - Rf[m2, 0, 1]) / S + q[rows, 1] = (Rf[m2, 0, 2] + Rf[m2, 2, 0]) / S + q[rows, 2] = (Rf[m2, 1, 2] + Rf[m2, 2, 1]) / S + q[rows, 3] = 0.25 * S + + # Normalize for numerical safety + q /= np.linalg.norm(q, axis=1, keepdims=True) + return q + + def quat_to_rotm_batch(q): + # q: (N,4) [w,x,y,z] -> R: (N,3,3) + w = q[:, 0] + x = q[:, 1] + y = q[:, 2] + z = q[:, 3] + x2 = x * x + y2 = y * y + z2 = z * z + wx = w * x + wy = w * y + wz = w * z + xy = x * y + xz = x * z + yz = y * z + R = np.zeros((q.shape[0], 3, 3), dtype=float) + R[:, 0, 0] = 1.0 - 2.0 * y2 - 2.0 * z2 + R[:, 1, 0] = 2.0 * xy + 2.0 * wz + R[:, 2, 0] = 2.0 * xz - 2.0 * wy + R[:, 0, 1] = 2.0 * xy - 2.0 * wz + R[:, 1, 1] = 1.0 - 2.0 * x2 - 2.0 * z2 + R[:, 2, 1] = 2.0 * yz + 2.0 * wx + R[:, 0, 2] = 2.0 * xz + 2.0 * wy + R[:, 1, 2] = 2.0 * yz - 2.0 * wx + R[:, 2, 2] = 1.0 - 2.0 * x2 - 2.0 * y2 + return R + + # Prepare inputs + t = np.clip(np.asarray(interp_func, dtype=float), 0.0, 1.0) + q1 = rotm_to_quat_batch(np.asarray(Q1, dtype=float)) + q2 = rotm_to_quat_batch(np.asarray(Q2, dtype=float)) + + # Ensure shortest path on the unit 4-sphere + dot = np.sum(q1 * q2, axis=1) + neg_mask = dot < 0.0 + if np.any(neg_mask): + q2[neg_mask] = -q2[neg_mask] + dot[neg_mask] = -dot[neg_mask] + + # SLERP weights + dot_clipped = np.clip(dot, -1.0, 1.0) + theta0 = np.arccos(dot_clipped) + sin_theta0 = np.sin(theta0) + + # Threshold for linear interpolation + lin_mask = sin_theta0 < 1e-6 + q = np.empty_like(q1) + + if np.any(~lin_mask): + theta = theta0[~lin_mask] * t[~lin_mask] + s0 = np.sin(theta0[~lin_mask] - theta) / sin_theta0[~lin_mask] + s1 = np.sin(theta) / sin_theta0[~lin_mask] + q[~lin_mask] = (s0[:, None] * q1[~lin_mask]) + (s1[:, None] * q2[~lin_mask]) + + if np.any(lin_mask): + # Nearly identical orientations: perform linear interpolation and normalize later + tl = t[lin_mask][:, None] + q[lin_mask] = (1.0 - tl) * q1[lin_mask] + tl * q2[lin_mask] + + # Normalize and convert back to rotation matrices + q /= np.linalg.norm(q, axis=1, keepdims=True) + return quat_to_rotm_batch(q) + +def axis(u, v): + """ + u, v: (nelems, 3) + return Q: (nelems, 3, 3) where columns are [e0 (circ), e1 (long), e2 (trans)] per element + """ + u = np.asarray(u, dtype=float) + v = np.asarray(v, dtype=float) + ne = u.shape[0] + + # e1 = normalize rows of u + e1 = normalize(u) + + # e2 = v - proj_{e1}(v) + proj = np.sum(e1 * v, axis=1)[:, None] * e1 + e2 = v - proj + e2 = normalize(e2) + + # e0 = cross(e1, e2) normalized + e0 = np.cross(e1, e2, axisa=1, axisb=1) + e0 = normalize(e0) + + Q = np.zeros((ne, 3, 3), dtype=float) + Q[:, :, 0] = e0 + Q[:, :, 1] = e1 + Q[:, :, 2] = e2 + + return Q + +def orient(Q, alpha, beta): + """ + Given an orthogonal matrix Q (ne,3,3), rotate each Q[i] by alpha[i] about + the z-axis and then by beta[i] about the x-axis. alpha and beta are arrays + of length ne. + """ + Q = np.asarray(Q, dtype=float) + ne = Q.shape[0] + + ca = np.cos(alpha) + sa = np.sin(alpha) + cb = np.cos(beta) + sb = np.sin(beta) + + # Rotation about z (Ra) and x (Rb) for each element + Ra = np.zeros((ne, 3, 3), dtype=float) + Ra[:, 0, 0] = ca + Ra[:, 0, 1] = -sa + Ra[:, 1, 0] = sa + Ra[:, 1, 1] = ca + Ra[:, 2, 2] = 1.0 + + # Rb = np.zeros((ne, 3, 3), dtype=float) + # Rb[:, 0, 0] = 1.0 + # Rb[:, 1, 1] = cb + # Rb[:, 1, 2] = sb + # Rb[:, 2, 1] = -sb + # Rb[:, 2, 2] = cb + + # Rb = np.zeros((ne, 3, 3), dtype=float) + # Rb[:, 0, 0] = 1.0 + # Rb[:, 1, 1] = cb + # Rb[:, 1, 2] = -sb + # Rb[:, 2, 1] = sb + # Rb[:, 2, 2] = cb + + Rb = np.zeros((ne, 3, 3), dtype=float) + Rb[:, 0, 0] = cb + Rb[:, 0, 2] = sb + Rb[:, 1, 1] = 1.0 + Rb[:, 2, 0] = -sb + Rb[:, 2, 2] = cb + + # Compose rotations and apply to Q per element + RaRb = np.einsum('nij,njk->nik', Ra, Rb) + Qt = np.einsum('nij,njk->nik', Q, RaRb) + + return Qt + +def getFiberDirectionsBayer(Phi_EPI, Phi_LV, Phi_RV, + gPhi_EPI, gPhi_LV, gPhi_RV, gPhi_AB, + params, intermediate=False): + ''' + Compute the fiber directions at the center of each cell + ''' + + # Unpack parameters + ALFA_END = np.deg2rad(params["ALFA_END"]) + ALFA_EPI = np.deg2rad(params["ALFA_EPI"]) + BETA_END = np.deg2rad(params["BETA_END"]) + BETA_EPI = np.deg2rad(params["BETA_EPI"]) + + print(" Computing fiber directions at cells") + + d = Phi_RV / (Phi_LV + Phi_RV) + alfaS = ALFA_END * (1 - d) - ALFA_END * d + betaS = BETA_END * (1 - d) - BETA_END * d + alfaW = ALFA_END * (1 - Phi_EPI) + ALFA_EPI * Phi_EPI + betaW = BETA_END * (1 - Phi_EPI) + BETA_EPI * Phi_EPI + + Q_LV0 = axis(gPhi_AB, -gPhi_LV) + Q_LV = orient(Q_LV0, alfaS, betaS) + Q_RV0 = axis(gPhi_AB, gPhi_RV) # Note that gPhi_RV points the other way + Q_RV = orient(Q_RV0, alfaS, -betaS) # Therefore, we need a minus in betaS + + Q_END = bislerp(Q_LV, Q_RV, d) + Q_END[d > 0.5,:,0] = -Q_END[d > 0.5,:,0] + Q_END[d > 0.5,:,2] = -Q_END[d > 0.5,:,2] + + Q_EPI0 = axis(gPhi_AB, gPhi_EPI) + Q_EPI = orient(Q_EPI0, alfaW, betaW) + + FST = bislerp(Q_END, Q_EPI, Phi_EPI) + + F = FST[:, :, 0] + S = FST[:, :, 1] + T = FST[:, :, 2] + + if intermediate: + return F, S, T, Q_LV[:,:,0], Q_RV[:,:,0], Q_END[:,:,0], Q_EPI[:,:,0] + + return F, S, T + + +def get_alpha_beta_angles_Bayer(F, Phi_EPI, Phi_LV, Phi_RV, + gPhi_EPI, gPhi_LV, gPhi_RV, gPhi_AB, + params): + ''' + Sanity check routine + Compute alpha and beta angles at cells given fiber directions F and Laplace gradients + + ''' + + ALFA_END = np.deg2rad(params["ALFA_END"]) + ALFA_EPI = np.deg2rad(params["ALFA_EPI"]) + + d = Phi_RV / (Phi_LV + Phi_RV) + alfaS = ALFA_END * (1 - d) - ALFA_END * d + alfaW = ALFA_END * (1 - Phi_EPI) + ALFA_EPI * Phi_EPI + + # Alpha angle + Q_LV = axis(gPhi_AB, -gPhi_LV) # (N,3,3) + Q_RV = axis(gPhi_AB, gPhi_RV) # (N,3,3) + Q = np.copy(Q_LV) + Q[d > 0.5,:,0] = -Q_RV[d > 0.5,:,0] + Q[d > 0.5,:,2] = -Q_RV[d > 0.5,:,2] + C = Q[:, :, 0] # (N,3) + L = Q[:, :, 1] # (N,3) + + # Angle in radians between F and circumferential vector + # alpha is positive in the direction of the longitudinal vector + cosang = np.clip(np.sum(F * C, axis=1), -1.0, 1.0) + sinang = np.clip(np.sum(F * L, axis=1), -1.0, 1.0) + alpha_angle = np.sign(sinang) * np.arccos(np.abs(cosang)) + + # Beta angle + Q_LV = orient(axis(gPhi_AB, -gPhi_LV), + alfaS, + 0) + Q_RV = orient(axis(gPhi_AB, gPhi_RV), + alfaS, + 0) + + Q_END = bislerp(Q_LV, Q_RV, d) + Q_END[d > 0.5,:,0] = -Q_END[d > 0.5,:,0] + Q_END[d > 0.5,:,2] = -Q_END[d > 0.5,:,2] + + Q_EPI0 = axis(gPhi_AB, gPhi_EPI) + Q_EPI = orient(Q_EPI0, alfaW, 0.0) + + Q = bislerp(Q_END, Q_EPI, Phi_EPI) + Cr = Q[:, :, 0] + Tr = Q[:, :, 2] + + # Angle in radians between F and rotated circumferential vector + # beta is negative in the direction of the transmural vector + cosang = np.clip(np.sum(F * Cr, axis=1), -1.0, 1.0) + sinang = np.clip(np.sum(F * Tr, axis=1), -1.0, 1.0) + beta_angle = - np.sign(sinang) * np.arccos(np.abs(cosang)) # Note the minus sign to match definition + + return np.rad2deg(alpha_angle), np.rad2deg(beta_angle), C, Cr + +def generate_fibers_BiV_Bayer_cells(outdir, laplace_results_file, params, return_angles=False, return_intermediate=False): + ''' + Generate fiber directions on a truncated BiV ventricular geometry using the + Laplace-Dirichlet rule-based method of Bayer et al. 2012 + + ARGS: + laplace_results_file : str + Path to the .vtu mesh with Laplace fields defined at nodes + params : dict + Dictionary of parameters for fiber generation + ''' + + t1 = time.time() + print("========================================================") + + # Load Laplace solution + result_mesh, Phi_EPI, Phi_LV, Phi_RV, Phi_AB, \ + gPhi_EPI, gPhi_LV, gPhi_RV, gPhi_AB = loadLaplaceSolnBayer(laplace_results_file) + + + # Write the fiber directions to a vtu files + output_mesh = copy.deepcopy(result_mesh) + # Ensure only FIB_DIR is present + for k in list(output_mesh.cell_data.keys()): + output_mesh.cell_data.remove(k) + for k in list(output_mesh.point_data.keys()): + output_mesh.point_data.remove(k) + + # Generate fiber directions + out = getFiberDirectionsBayer(Phi_EPI, Phi_LV, Phi_RV, + gPhi_EPI, gPhi_LV, gPhi_RV, gPhi_AB, + params, intermediate=return_intermediate) + + if return_intermediate: + F, S, T, eC_LV, eC_RV, eC_END, eC_EPI = out + result_mesh.cell_data['F'] = F + result_mesh.cell_data['S'] = S + result_mesh.cell_data['T'] = T + result_mesh.cell_data['eC_LV'] = eC_LV + result_mesh.cell_data['eC_RV'] = eC_RV + result_mesh.cell_data['eC_END'] = eC_END + result_mesh.cell_data['eC_EPI'] = eC_EPI + else: + F, S, T = out + result_mesh.cell_data['F'] = F + result_mesh.cell_data['S'] = S + result_mesh.cell_data['T'] = T + + print(" Writing domains and fibers to VTK data structure") + + + fname1 = os.path.join(outdir, "fibersLong.vtu") + print(" Writing to vtu file ---> %s" % (fname1)) + output_mesh.cell_data.set_array(F, 'FIB_DIR') + output_mesh.save(fname1) + + fname1 = os.path.join(outdir, "fibersSheet.vtu") + print(" Writing to vtu file ---> %s" % (fname1)) + output_mesh.cell_data.set_array(T, 'FIB_DIR') + output_mesh.save(fname1) + + fname1 = os.path.join(outdir, "fibersNormal.vtu") + print(" Writing to vtu file ---> %s" % (fname1)) + output_mesh.cell_data.set_array(S, 'FIB_DIR') + output_mesh.save(fname1) + + t2 = time.time() + print('\n Total time: %.3fs' % (t2-t1)) + print("========================================================") + + if return_angles: + alpha_angle, beta_angle, eC, eCr = get_alpha_beta_angles_Bayer(F, Phi_EPI, Phi_LV, Phi_RV, + gPhi_EPI, gPhi_LV, gPhi_RV, gPhi_AB, + params) + result_mesh.cell_data['Alpha_Angle'] = alpha_angle + result_mesh.cell_data['Beta_Angle'] = beta_angle + result_mesh.cell_data['eC'] = eC + result_mesh.cell_data['eCr'] = eCr + + + return result_mesh + + + +def loadLaplaceSolnDoste(fileName): + ''' + Load a solution to a Laplace-Dirichlet problem from a .vtu file and extract + the solution and its gradients at the cells. + + ARGS: + fileName : str + Path to the .vtu file with the Laplace solution. The solution should be + defined at the nodes. + + Returns: + lap : dict + Dictionary of Laplace solution at cells + grad : dict + Dictionary of gradients at cells + ''' + + varnames = ['Trans_BiV', 'Long_AV', 'Long_MV', 'Long_PV', 'Long_TV', 'Weight_LV', + 'Weight_RV', 'Trans_EPI', 'Trans_LV', 'Trans_RV'] + + print(" Loading Laplace solution <--- %s" % (fileName)) + + # Read mesh with pyvista + result_mesh = pv.read(fileName) + + # Convert point-data to cell-data (keep point data passed to cells) + mesh_cells = result_mesh.point_data_to_cell_data() + + print(" Extracting solution and estimating gradients at cells") + + # Map VTU array names to internal keys expected by downstream code + name_map = { + 'Trans_BiV': 'ven_trans', + 'Long_AV': 'lv_av_long', + 'Long_MV': 'lv_mv_long', + 'Long_PV': 'rv_pv_long', + 'Long_TV': 'rv_tv_long', + 'Weight_LV': 'lv_weight', + 'Weight_RV': 'rv_weight', + 'Weight_RV_OP': 'rv_op_weight', + 'Trans_EPI': 'epi_trans', + 'Trans_LV': 'lv_trans', + 'Trans_RV': 'rv_trans', + } + + lap = {} + grad = {} + + for vname in varnames: + if vname not in mesh_cells.cell_data: + print(f" Warning: '{vname}' not found in cell_data; skipping") + continue + key = name_map[vname] + + # Cell-centered Laplace values + lap[key] = np.asarray(mesh_cells.cell_data[vname]) + + # Cell-centered gradients via PyVista + gmesh = mesh_cells.compute_derivative(scalars=vname, gradient=True, preference='cell') + grad[key] = np.asarray(gmesh.cell_data['gradient']) + + return result_mesh, lap, grad + + +def compute_basis_vectors(lap, grad): + # LV + # longitudinal + lv_glong = grad['lv_mv_long']*lap['lv_weight'][:,None] + grad['lv_av_long']*(1 - lap['lv_weight'][:,None]) + eL_lv = normalize(lv_glong) + + # transmural + lv_gtrans = grad['lv_trans'] - (eL_lv*grad['lv_trans'])*eL_lv + eT_lv = normalize(lv_gtrans) + + # circumferential + eC_lv = np.cross(eL_lv, eT_lv, axisa=1, axisb=1) + eC_lv = normalize(eC_lv) + + # Ensuring orthogonality + eT_lv = np.cross(eC_lv, eL_lv, axisa=1, axisb=1) + eT_lv = normalize(eT_lv) + + # RV + # longitudinal + rv_glong = grad['rv_tv_long']*lap['rv_weight'][:,None] + grad['rv_pv_long']*(1 - lap['rv_weight'][:,None] ) + eL_rv = normalize(rv_glong) + + # transmural + rv_gtrans = grad['rv_trans'] - (eL_rv*grad['rv_trans'])*eL_rv + eT_rv = normalize(rv_gtrans) + + # circumferential + eC_rv = np.cross(eL_rv, eT_rv, axisa=1, axisb=1) + eC_rv = normalize(eC_rv) + + # Ensuring orthogonality + eT_rv = np.cross(eC_rv, eL_rv, axisa=1, axisb=1) + eT_rv = normalize(eT_rv) + + # Write out global circumferential vector + eC = eC_rv*(1-lap['ven_trans'][:,None]) + eC_lv*lap['ven_trans'][:,None] + eC = normalize(eC) + + basis = {'eC_lv': eC_lv, + 'eT_lv': eT_lv, + 'eL_lv': eL_lv, + 'eC_rv': eC_rv, + 'eT_rv': eT_rv, + 'eL_rv': eL_rv, + 'eC': eC} + + return basis + + +def redistribute_weight(weight, up, low, strategy='centre'): + new_weight = weight.copy() + + if strategy == 'flip': + # Shift all weights + new_mean = 1 - np.mean(weight) + shift = new_mean - np.mean(weight) + new_weight = new_weight + shift + + # Cut off values outside of range 0 - 1 + new_weight[new_weight > 1] = 1 + new_weight[new_weight < 0] = 0 + + # Redistribute new tail + new_weight = (new_weight - np.min(new_weight)) / (np.max(new_weight) - np.min(new_weight)) + tmp = new_weight.copy() + + if shift > 0: + tmp[tmp >= new_mean] = np.nan + tmp = (tmp - np.nanmin(tmp)) / (new_mean - np.nanmin(tmp)) + elif shift < 0: + tmp[tmp <= new_mean] = np.nan + tmp = (tmp - new_mean) / (np.nanmax(tmp) - new_mean) + + tmp[np.isnan(tmp)] = new_weight[np.isnan(tmp)] + new_weight = tmp + + else: # cut off tails so that weights are centered + # Find upper and lower limits + upper_lim = np.quantile(weight, up) + while upper_lim == 0: + print('Upper limit is 0, increasing upper limit') + up += 0.1 + upper_lim = np.quantile(weight, up) + + lower_lim = np.quantile(weight, low) + + # Set upper and lower values to limits + new_weight[new_weight > upper_lim] = upper_lim + new_weight[new_weight < lower_lim] = lower_lim + + # Redistribute/normalize values + new_weight = (new_weight - np.min(new_weight)) / (np.max(new_weight) - np.min(new_weight)) + + return new_weight + + +def compute_alpha_beta_angles(lap, params): + # Modify weights so the effect of outflow tracts is localized + lv_weight = redistribute_weight(lap['lv_weight'], 0.7, 0.01) + rv_weight = redistribute_weight(lap['rv_weight'], 0.1, 0.001) + + # LV + alpha_lv_endo_long = params['AENDOLV'] * lv_weight + params['AOTENDOLV'] * (1 - lv_weight) # Endo + alpha_lv_epi_long = params['AEPILV'] * lv_weight + params['AOTEPILV'] * (1 - lv_weight) + + alpha_wall_lv = alpha_lv_endo_long * (1 - lap['epi_trans']) + alpha_lv_epi_long * lap['epi_trans'] + beta_wall_lv = (params['BENDOLV'] * (1 - lap['epi_trans']) + params['BEPILV'] * lap['epi_trans']) * lv_weight + + # RV + alpha_rv_endo_long = params['AENDORV'] * rv_weight + params['AOTENDORV'] * (1 - rv_weight) + alpha_rv_epi_long = params['AEPIRV'] * rv_weight + params['AOTEPIRV'] * (1 - rv_weight) + + alpha_wall_rv = alpha_rv_endo_long * (1 - lap['epi_trans']) + alpha_rv_epi_long * lap['epi_trans'] + beta_wall_rv = (params['BENDORV'] * (1 - lap['epi_trans']) + params['BEPIRV'] * lap['epi_trans']) * rv_weight + + # Septum + sep = np.abs(lap['ven_trans'] - 0.5) + sep = (sep - np.min(sep)) / (np.max(sep) - np.min(sep)) + alpha_septum = (alpha_lv_endo_long * sep * lap['lv_trans']) + (alpha_rv_endo_long * sep * lap['rv_trans']) + beta_septum = (params['BENDOLV'] * lap['lv_trans'] * lv_weight) + (params['BENDORV'] * lap['rv_trans'] * rv_weight) + + angles = {'alpha_lv_endo_long': alpha_lv_endo_long, + 'alpha_lv_epi_long': alpha_lv_epi_long, + 'alpha_wall_lv': alpha_wall_lv, + 'beta_wall_lv': beta_wall_lv, + 'alpha_rv_endo_long': alpha_rv_endo_long, + 'alpha_rv_epi_long': alpha_rv_epi_long, + 'alpha_wall_rv': alpha_wall_rv, + 'beta_wall_rv': beta_wall_rv, + 'alpha_septum': alpha_septum, + 'beta_septum': beta_septum + } + + return angles + + +def rotate_basis(eC, eL, eT, alpha, beta): + eC = normalize(eC) + eT = normalize(eT) + eL = normalize(eL) + + # Matrix of directional vectors + Q = np.stack([eC, eL, eT], axis=-1) + Q = np.transpose(Q, (2, 1, 0)) + + # Create rotation matrix - from Doste code + axis = eT + R = np.array([[np.cos(alpha) + (axis[:, 0]**2)*(1 - np.cos(alpha)), axis[:,0] * axis[:,1]*(1 - np.cos(alpha)) - axis[:,2]*np.sin(alpha), axis[:,0]*axis[:,2]*(1 - np.cos(alpha)) + axis[:,1]*np.sin(alpha)], + [axis[:,1]*axis[:,0]*(1 - np.cos(alpha)) + axis[:,2]*np.sin(alpha), np.cos(alpha) + (axis[:,1]**2)*(1 - np.cos(alpha)), axis[:,1]*axis[:, 2]*(1 - np.cos(alpha)) - axis[:, 0]*np.sin(alpha)], + [axis[:,2]*axis[:,0]*(1 - np.cos(alpha)) - axis[:,1]*np.sin(alpha), axis[:,2]*axis[:,1]*(1 - np.cos(alpha)) + axis[:, 0]*np.sin(alpha), np.cos(alpha)+(axis[:, 2]**2)*(1 - np.cos(alpha))]]) + + # Rotate the circumferential direction around the transmural direction + QX = np.zeros_like(R) + for i in range(len(eC)): + QX[:, :, i] = np.matmul(Q[:, :, i], R[:, :, i]) + + # Second rotation (beta) about QX + axis2 = QX[1, :, :].T + R2 = np.array([ + [np.cos(beta) + (axis2[:,0]**2)*(1 - np.cos(beta)), axis2[:,0]*axis2[:, 1]*(1 - np.cos(beta)) - axis2[:,2] * np.sin(beta), axis2[:,0] * axis2[:,2] * (1 - np.cos(beta)) + axis2[:,1] * np.sin(beta)], + [axis2[:,1] * axis2[:,0]*(1 - np.cos(beta)) + axis2[:,2]*np.sin(beta), np.cos(beta) + (axis2[:,1]**2)*(1 - np.cos(beta)), axis2[:,1] * axis2[:,2] * (1 - np.cos(beta)) - axis2[:,0] * np.sin(beta)], + [axis2[:,2] * axis2[:,0]*(1 - np.cos(beta)) - axis2[:,1]*np.sin(beta), axis2[:, 2] * axis2[:,1] * (1 - np.cos(beta)) + axis2[:,0] * np.sin(beta), np.cos(beta) + (axis2[:,2]**2) * (1 - np.cos(beta))] + ]) + + QX2 = np.zeros((R.shape[2], 3, 3), dtype=float) + for i in range(len(eC)): + QX2[i] = np.matmul(QX[:, :, i], R2[:, :, i]).T + + return QX2 + + +def compute_local_basis(basis, angles): + Qlv_septum = rotate_basis(basis['eC_lv'], basis['eL_lv'], basis['eT_lv'], angles['alpha_septum'], angles['beta_septum']) + Qrv_septum = rotate_basis(basis['eC_rv'], basis['eL_rv'], basis['eT_rv'], angles['alpha_septum'], angles['beta_septum']) + Qlv_epi = rotate_basis(basis['eC_lv'], basis['eL_lv'], basis['eT_lv'], angles['alpha_wall_lv'], angles['beta_wall_lv']) + Qrv_epi = rotate_basis(basis['eC_rv'], basis['eL_rv'], basis['eT_rv'], angles['alpha_wall_rv'], angles['beta_wall_rv']) + + local_basis = {'Qlv_septum': Qlv_septum, + 'Qrv_septum': Qrv_septum, + 'Qlv_epi': Qlv_epi, + 'Qrv_epi': Qrv_epi, + } + + return local_basis + + +def interpolate_local_basis(lap, local_basis): + + epi_trans = lap['epi_trans'] + + Qrv_septum = local_basis['Qrv_septum'] + Qlv_septum = local_basis['Qlv_septum'] + Qrv_epi = local_basis['Qrv_epi'] + Qlv_epi = local_basis['Qlv_epi'] + + Qepi = bislerp(Qrv_epi, Qlv_epi, lap['ven_trans']) + Qendo = bislerp(Qrv_septum, Qlv_septum, lap['ven_trans']) + Q = bislerp(Qendo, Qepi, epi_trans) + + return Q, Qepi + + +def getFiberDirectionsDoste(lap, grad, params, intermediate=False): + # Convert parameters from degrees to radians + for key in params: + params[key] = np.deg2rad(params[key]) + + print('Computing basis vectors') + basis = compute_basis_vectors(lap, grad) + + print('Computing angles') + angles = compute_alpha_beta_angles(lap, params) + + print('Computing local basis') + local_basis = compute_local_basis(basis, angles) + + print('Interpolating basis') + Q, Qepi = interpolate_local_basis(lap, local_basis) + + print('Done!') + f = Q[:, :, 0] + s = Q[:, :, 1] + n = Q[:, :, 2] + + if intermediate: + return f, s, n, basis, angles, local_basis, Qepi[:,:,0] + + return f, s, n + + +def generate_fibers_BiV_Doste_cells(outdir, laplace_results_file, params, return_angles=False, return_intermediate=False): + ''' + Generate fiber directions on a truncated BiV ventricular geometry using the + Laplace-Dirichlet rule-based method of Bayer et al. 2012 + + ARGS: + laplace_results_file : str + Path to the .vtu mesh with Laplace fields defined at nodes + params : dict + Dictionary of parameters for fiber generation + ''' + + t1 = time.time() + print("========================================================") + + # Load Laplace solution + result_mesh,lap, grad = loadLaplaceSolnDoste(laplace_results_file) + + # Write the fiber directions to a vtu files + output_mesh = copy.deepcopy(result_mesh) + # Ensure only FIB_DIR is present + for k in list(output_mesh.cell_data.keys()): + output_mesh.cell_data.remove(k) + for k in list(output_mesh.point_data.keys()): + output_mesh.point_data.remove(k) + + # Generate fiber directions + out = getFiberDirectionsDoste(lap, grad, params, intermediate=return_intermediate) + + if return_intermediate: + F, S, T, basis, angles, local_basis, Qepi = out + result_mesh.cell_data['F'] = F + result_mesh.cell_data['S'] = S + result_mesh.cell_data['T'] = T + result_mesh.cell_data['Qepi'] = Qepi + for k, v in basis.items(): + result_mesh.cell_data[k] = v + for k, v in angles.items(): + result_mesh.cell_data[k] = v + for k, v in local_basis.items(): + # Flatten local basis matrices to store + flattened = v[:, :, 0] + result_mesh.cell_data[k] = flattened + else: + F, S, T = out + result_mesh.cell_data['F'] = F + result_mesh.cell_data['S'] = S + result_mesh.cell_data['T'] = T + + print(" Writing domains and fibers to VTK data structure") + + + fname1 = os.path.join(outdir, "fibersLong.vtu") + print(" Writing to vtu file ---> %s" % (fname1)) + output_mesh.cell_data.set_array(F, 'FIB_DIR') + output_mesh.save(fname1) + + fname1 = os.path.join(outdir, "fibersSheet.vtu") + print(" Writing to vtu file ---> %s" % (fname1)) + output_mesh.cell_data.set_array(T, 'FIB_DIR') + output_mesh.save(fname1) + + fname1 = os.path.join(outdir, "fibersNormal.vtu") + print(" Writing to vtu file ---> %s" % (fname1)) + output_mesh.cell_data.set_array(S, 'FIB_DIR') + output_mesh.save(fname1) + + t2 = time.time() + print('\n Total time: %.3fs' % (t2-t1)) + print("========================================================") + + if return_angles: + alpha_angle, beta_angle, eC, eCr = get_alpha_beta_angles_Doste(F, lap, grad, params) + result_mesh.cell_data['Alpha_Angle'] = alpha_angle + result_mesh.cell_data['Beta_Angle'] = beta_angle + result_mesh.cell_data['eC'] = eC + result_mesh.cell_data['eCr'] = eCr + + + return result_mesh + + + +def get_alpha_beta_angles_Doste(F, lap, grad, params): + ''' + Sanity check routine for Doste-based fibers. + Compute alpha and beta angles at cells given fiber directions F and the + Laplace/basis fields used by the Doste method. + + Returns: + alpha_angle_deg, beta_angle_deg, eC_ref, Cr_ref + - alpha, beta in degrees + - eC_ref: reference circumferential vector (before rotations) + - Cr_ref: circumferential vector after applying only alpha rotation + ''' + + # Reconstruct base vectors used by Doste + basis = compute_basis_vectors(lap, grad) + eC_global = basis['eC'] # blended circumferential + + # Build a blended longitudinal direction for alpha sign (LV/RV mix) + ven = lap['ven_trans'][:, None] + eL_blend = normalize(basis['eL_rv'] * (1.0 - ven) + basis['eL_lv'] * ven) + + # Alpha: signed angle between F and circumferential in the tangent plane, + # sign taken along longitudinal direction (consistent with Bayer routine, + # but negative because the longitudinal direction is opposite to the Doste paper) + cos_a = np.clip(np.sum(F * eC_global, axis=1), -1.0, 1.0) + sin_a = np.clip(np.sum(F * eL_blend, axis=1), -1.0, 1.0) + alpha_angle = -np.sign(sin_a) * np.arccos(np.abs(cos_a)) + + # Build reference frame after ONLY alpha rotation (beta = 0) + angles = compute_alpha_beta_angles(lap, params) # radians + Qlv_septum_a = rotate_basis(basis['eC_lv'], basis['eL_lv'], basis['eT_lv'], + angles['alpha_septum'], 0.0) + Qrv_septum_a = rotate_basis(basis['eC_rv'], basis['eL_rv'], basis['eT_rv'], + angles['alpha_septum'], 0.0) + Qlv_epi_a = rotate_basis(basis['eC_lv'], basis['eL_lv'], basis['eT_lv'], + angles['alpha_wall_lv'], 0.0) + Qrv_epi_a = rotate_basis(basis['eC_rv'], basis['eL_rv'], basis['eT_rv'], + angles['alpha_wall_rv'], 0.0) + + Qepi_a = bislerp(Qrv_epi_a, Qlv_epi_a, lap['ven_trans']) + Qendo_a = bislerp(Qrv_septum_a, Qlv_septum_a, lap['ven_trans']) + Qa = bislerp(Qendo_a, Qepi_a, lap['epi_trans']) + + Cr = Qa[:, :, 0] # circumferential after alpha-only rotation + Tr = Qa[:, :, 2] # transmural after alpha-only rotation + + # Beta: signed angle between F and Cr, sign w.r.t. Tr (negative by convention) + cos_b = np.clip(np.sum(F * Cr, axis=1), -1.0, 1.0) + sin_b = np.clip(np.sum(F * Tr, axis=1), -1.0, 1.0) + beta_angle = np.sign(sin_b) * np.arccos(np.abs(cos_b)) + + return np.rad2deg(alpha_angle), np.rad2deg(beta_angle), eC_global, Cr + diff --git a/utilities/fiber_generation/src/templates/solver_bayer.xml b/utilities/fiber_generation/src/templates/solver_bayer.xml new file mode 100644 index 000000000..fe42dd1a7 --- /dev/null +++ b/utilities/fiber_generation/src/templates/solver_bayer.xml @@ -0,0 +1,266 @@ + + + + + 0 + 3 + 1 + 1 + 0. + STOP_SIM + + 1 + result + 1 + /Users/jjv/Research/CEPCRTProject/ct_1010/heat/fibers + 1 + + 1 + 0 + + 1 + 1 + 0 + + + + + /Users/jjv/Research/CEPCRTProject/ct_1010/VOLUME.vtu + + /Users/jjv/Research/CEPCRTProject/ct_1010/mesh-surfaces/EPI.vtp + + + /Users/jjv/Research/CEPCRTProject/ct_1010/mesh-surfaces/BASE.vtp + + + /Users/jjv/Research/CEPCRTProject/ct_1010/mesh-surfaces/EPI_APEX.vtp + + + /Users/jjv/Research/CEPCRTProject/ct_1010/mesh-surfaces/LV.vtp + + + /Users/jjv/Research/CEPCRTProject/ct_1010/mesh-surfaces/RV.vtp + + + + + 0 + 1 + 5 + 1e-6 + + 1.0 + 0.0 + 0.0 + + + 1 + + + + Phi_BiV_EPI + + + + + fsils + + 500 + 1e-6 + + + + Dir + 1.0 + 0 + + + + Dir + 0.0 + 0 + + + + Dir + 0.0 + 0 + + + + + 0 + 1 + 5 + 1e-6 + + 1.0 + 0.0 + 0.0 + + + 1 + + + + Phi_BiV_LV + + + + + fsils + + 500 + 1e-6 + + + + Dir + 1.0 + 0 + + + + Dir + 0.0 + 0 + + + + Dir + 0.0 + 0 + + + + + 0 + 1 + 5 + 1e-6 + + 1.0 + 0.0 + 0.0 + + + 1 + + + + Phi_BiV_RV + + + + + fsils + + 500 + 1e-6 + + + + Dir + 1.0 + 0 + + + + Dir + 0.0 + 0 + + + + Dir + 0.0 + 0 + + + + + 0 + 1 + 5 + 1e-6 + + 1.0 + 0.0 + 0.0 + + + 1 + + + + Phi_BiV_AB + + + + + fsils + + 500 + 1e-6 + + + + Dir + 1.0 + 0 + + + + Dir + 0.0 + 0 + + + + + 0 + 1 + 5 + 1e-6 + + 1.0 + 0.0 + 0.0 + + + 1 + + + + Phi_BiV_D + + + + + fsils + + 500 + 1e-6 + + + + Dir + 2.0 + 0 + + + + Dir + -1.0 + 0 + + + + Dir + 0.0 + 0 + + + + + \ No newline at end of file diff --git a/utilities/fiber_generation/src/templates/solver_doste.xml b/utilities/fiber_generation/src/templates/solver_doste.xml new file mode 100644 index 000000000..fbdeb8412 --- /dev/null +++ b/utilities/fiber_generation/src/templates/solver_doste.xml @@ -0,0 +1,492 @@ + + + + + 0 + 3 + 1 + 1 + 0. + STOP_SIM + + 1 + result + 1 + /Users/jjv/Research/CEPCRTProject/ct_1010/heat/fibers + 1 + + 1 + 0 + + 1 + 1 + 0 + + + + + /Users/jjv/Research/CEPCRTProject/ct_1010/VOLUME.vtu + + /Users/jjv/Research/CEPCRTProject/ct_1010/mesh-surfaces/EPI.vtp + + + /Users/jjv/Research/CEPCRTProject/ct_1010/mesh-surfaces/BASE.vtp + + + /Users/jjv/Research/CEPCRTProject/ct_1010/mesh-surfaces/BASE.vtp + + + /Users/jjv/Research/CEPCRTProject/ct_1010/mesh-surfaces/BASE.vtp + + + /Users/jjv/Research/CEPCRTProject/ct_1010/mesh-surfaces/BASE.vtp + + + /Users/jjv/Research/CEPCRTProject/ct_1010/mesh-surfaces/EPI_APEX.vtp + + + /Users/jjv/Research/CEPCRTProject/ct_1010/mesh-surfaces/LV.vtp + + + /Users/jjv/Research/CEPCRTProject/ct_1010/mesh-surfaces/RV.vtp + + + + + + + + 0 + 1 + 5 + 1e-6 + + 1.0 + 0.0 + 0.0 + + + 1 + + + + Trans_BiV + + + + + fsils + + 500 + 1e-6 + + + + Dir + 1.0 + 0 + + + + Dir + 0.0 + 0 + + + + + + + + + + 0 + 1 + 5 + 1e-6 + + 1.0 + 0.0 + 0.0 + + + 1 + + + + Long_AV + + + + + fsils + + 500 + 1e-6 + + + + Dir + 1.0 + 0 + + + + Dir + 0.0 + 0 + + + + + + + + + + 0 + 1 + 5 + 1e-6 + + 1.0 + 0.0 + 0.0 + + + 1 + + + + Long_MV + + + + + fsils + + 500 + 1e-6 + + + + Dir + 1.0 + 0 + + + + Dir + 0.0 + 0 + + + + + + + + + 0 + 1 + 5 + 1e-6 + + 1.0 + 0.0 + 0.0 + + + 1 + + + + Long_PV + + + + + fsils + + 500 + 1e-6 + + + + Dir + 1.0 + 0 + + + + Dir + 0.0 + 0 + + + + + + + + + 0 + 1 + 5 + 1e-6 + + 1.0 + 0.0 + 0.0 + + + 1 + + + + Long_TV + + + + + fsils + + 500 + 1e-6 + + + + Dir + 1.0 + 0 + + + + Dir + 0.0 + 0 + + + + + + + + + 0 + 1 + 5 + 1e-6 + + 1.0 + 0.0 + 0.0 + + + 1 + + + + Weight_LV + + + + + fsils + + 500 + 1e-6 + + + + Dir + 0.0 + 0 + + + + Dir + 1.0 + 0 + + + + + + + + 0 + 1 + 5 + 1e-6 + + 1.0 + 0.0 + 0.0 + + + 1 + + + + Weight_RV + + + + + fsils + + 500 + 1e-6 + + + + Dir + 0.0 + 0 + + + + Dir + 1.0 + 0 + + + + + + + + 0 + 1 + 5 + 1e-6 + + 1.0 + 0.0 + 0.0 + + + 1 + + + + Trans_EPI + + + + + fsils + + 500 + 1e-6 + + + + Dir + 1.0 + 0 + + + Dir + 0.0 + 0 + + + Dir + 0.0 + 0 + + + + + + + 0 + 1 + 5 + 1e-6 + + 1.0 + 0.0 + 0.0 + + + 1 + + + + Trans_LV + + + + + fsils + + 500 + 1e-6 + + + + Dir + 0.0 + 0 + + + Dir + 0.0 + 0 + + + Dir + 1.0 + 0 + + + + + + + 0 + 1 + 5 + 1e-6 + + 1.0 + 0.0 + 0.0 + + + 1 + + + + Trans_RV + + + + + fsils + + 500 + 1e-6 + + + + Dir + 0.0 + 0 + + + Dir + 0.0 + 0 + + + Dir + 1.0 + 0 + + + + \ No newline at end of file From c36a6b89212aea6b116f71462fda831710de5296 Mon Sep 17 00:00:00 2001 From: Javiera Jilberto Vallejos Date: Tue, 9 Dec 2025 11:29:03 -0800 Subject: [PATCH 2/5] increasing linear solver max iterations --- .../src/templates/solver_bayer.xml | 10 +++++----- .../src/templates/solver_doste.xml | 20 +++++++++---------- 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/utilities/fiber_generation/src/templates/solver_bayer.xml b/utilities/fiber_generation/src/templates/solver_bayer.xml index fe42dd1a7..abf7e803d 100644 --- a/utilities/fiber_generation/src/templates/solver_bayer.xml +++ b/utilities/fiber_generation/src/templates/solver_bayer.xml @@ -65,7 +65,7 @@ fsils - 500 + 2000 1e-6 @@ -110,7 +110,7 @@ fsils - 500 + 2000 1e-6 @@ -155,7 +155,7 @@ fsils - 500 + 2000 1e-6 @@ -200,7 +200,7 @@ fsils - 500 + 2000 1e-6 @@ -239,7 +239,7 @@ fsils - 500 + 2000 1e-6 diff --git a/utilities/fiber_generation/src/templates/solver_doste.xml b/utilities/fiber_generation/src/templates/solver_doste.xml index fbdeb8412..7f9dc0f4a 100644 --- a/utilities/fiber_generation/src/templates/solver_doste.xml +++ b/utilities/fiber_generation/src/templates/solver_doste.xml @@ -77,7 +77,7 @@ fsils - 500 + 2000 1e-6 @@ -121,7 +121,7 @@ fsils - 500 + 2000 1e-6 @@ -165,7 +165,7 @@ fsils - 500 + 2000 1e-6 @@ -208,7 +208,7 @@ fsils - 500 + 2000 1e-6 @@ -251,7 +251,7 @@ fsils - 500 + 2000 1e-6 @@ -294,7 +294,7 @@ fsils - 500 + 2000 1e-6 @@ -336,7 +336,7 @@ fsils - 500 + 2000 1e-6 @@ -378,7 +378,7 @@ fsils - 500 + 2000 1e-6 @@ -423,7 +423,7 @@ fsils - 500 + 2000 1e-6 @@ -468,7 +468,7 @@ fsils - 500 + 2000 1e-6 From c401a9af3c4743daaca51b943f60301d8670ed8b Mon Sep 17 00:00:00 2001 From: Javiera Jilberto Vallejos Date: Mon, 29 Dec 2025 09:40:00 -0800 Subject: [PATCH 3/5] changing paths and README.md --- utilities/fiber_generation/README.md | 2 +- .../src/templates/solver_bayer.xml | 14 ++++++------- .../src/templates/solver_doste.xml | 20 +++++++++---------- 3 files changed, 18 insertions(+), 18 deletions(-) diff --git a/utilities/fiber_generation/README.md b/utilities/fiber_generation/README.md index 41379c364..fad7dfa0b 100644 --- a/utilities/fiber_generation/README.md +++ b/utilities/fiber_generation/README.md @@ -1,5 +1,5 @@ # SV-fibergen -Python + SVmultiphysics codes for fiber generation. Two methods are implemented: +Python + svMultiPhysics codes for fiber generation. Two methods are implemented: * Bayer et al. (2012). [link](https://doi.org/10.1007/s10439-012-0593-5) * Doste et al. (2018). [link](https://doi.org/10.1002/cnm.3185) diff --git a/utilities/fiber_generation/src/templates/solver_bayer.xml b/utilities/fiber_generation/src/templates/solver_bayer.xml index abf7e803d..e98f80a08 100644 --- a/utilities/fiber_generation/src/templates/solver_bayer.xml +++ b/utilities/fiber_generation/src/templates/solver_bayer.xml @@ -12,7 +12,7 @@ 1 result 1 - /Users/jjv/Research/CEPCRTProject/ct_1010/heat/fibers + PATH_TO_OUTPUT_FLDR 1 1 @@ -25,21 +25,21 @@ - /Users/jjv/Research/CEPCRTProject/ct_1010/VOLUME.vtu + VOLUME.vtu - /Users/jjv/Research/CEPCRTProject/ct_1010/mesh-surfaces/EPI.vtp + EPI.vtp - /Users/jjv/Research/CEPCRTProject/ct_1010/mesh-surfaces/BASE.vtp + BASE.vtp - /Users/jjv/Research/CEPCRTProject/ct_1010/mesh-surfaces/EPI_APEX.vtp + EPI_APEX.vtp - /Users/jjv/Research/CEPCRTProject/ct_1010/mesh-surfaces/LV.vtp + LV.vtp - /Users/jjv/Research/CEPCRTProject/ct_1010/mesh-surfaces/RV.vtp + RV.vtp diff --git a/utilities/fiber_generation/src/templates/solver_doste.xml b/utilities/fiber_generation/src/templates/solver_doste.xml index 7f9dc0f4a..00a9c01f0 100644 --- a/utilities/fiber_generation/src/templates/solver_doste.xml +++ b/utilities/fiber_generation/src/templates/solver_doste.xml @@ -12,7 +12,7 @@ 1 result 1 - /Users/jjv/Research/CEPCRTProject/ct_1010/heat/fibers + PATH_TO_OUTPUT_FLDR 1 1 @@ -25,30 +25,30 @@ - /Users/jjv/Research/CEPCRTProject/ct_1010/VOLUME.vtu + VOLUME.vtu - /Users/jjv/Research/CEPCRTProject/ct_1010/mesh-surfaces/EPI.vtp + EPI.vtp - /Users/jjv/Research/CEPCRTProject/ct_1010/mesh-surfaces/BASE.vtp + MV.vtp - /Users/jjv/Research/CEPCRTProject/ct_1010/mesh-surfaces/BASE.vtp + AV.vtp - /Users/jjv/Research/CEPCRTProject/ct_1010/mesh-surfaces/BASE.vtp + TV.vtp - /Users/jjv/Research/CEPCRTProject/ct_1010/mesh-surfaces/BASE.vtp + PV.vtp - /Users/jjv/Research/CEPCRTProject/ct_1010/mesh-surfaces/EPI_APEX.vtp + EPI_APEX.vtp - /Users/jjv/Research/CEPCRTProject/ct_1010/mesh-surfaces/LV.vtp + LV.vtp - /Users/jjv/Research/CEPCRTProject/ct_1010/mesh-surfaces/RV.vtp + RV.vtp From f15a13629138ab87de96aaae2fed72243cec4427 Mon Sep 17 00:00:00 2001 From: Javiera Jilberto Vallejos Date: Mon, 29 Dec 2025 10:40:29 -0800 Subject: [PATCH 4/5] adding .pngs for fiber, sheet, sheetnormal directions --- utilities/fiber_generation/example/ot/doste.png | 3 --- utilities/fiber_generation/example/ot/doste_fibers.png | 3 +++ utilities/fiber_generation/example/ot/doste_sheet.png | 3 +++ utilities/fiber_generation/example/ot/doste_sheetnormal.png | 3 +++ utilities/fiber_generation/example/truncated/bayer.png | 3 --- utilities/fiber_generation/example/truncated/bayer_fiber.png | 3 +++ utilities/fiber_generation/example/truncated/bayer_sheet.png | 3 +++ .../fiber_generation/example/truncated/bayer_sheetnormal.png | 3 +++ 8 files changed, 18 insertions(+), 6 deletions(-) delete mode 100644 utilities/fiber_generation/example/ot/doste.png create mode 100644 utilities/fiber_generation/example/ot/doste_fibers.png create mode 100644 utilities/fiber_generation/example/ot/doste_sheet.png create mode 100644 utilities/fiber_generation/example/ot/doste_sheetnormal.png delete mode 100644 utilities/fiber_generation/example/truncated/bayer.png create mode 100644 utilities/fiber_generation/example/truncated/bayer_fiber.png create mode 100644 utilities/fiber_generation/example/truncated/bayer_sheet.png create mode 100644 utilities/fiber_generation/example/truncated/bayer_sheetnormal.png diff --git a/utilities/fiber_generation/example/ot/doste.png b/utilities/fiber_generation/example/ot/doste.png deleted file mode 100644 index 1ea35dfb1..000000000 --- a/utilities/fiber_generation/example/ot/doste.png +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:4da5c28acd2e2b9745ba972ec783b7341afa2156a48244059684ca28f45652ac -size 494459 diff --git a/utilities/fiber_generation/example/ot/doste_fibers.png b/utilities/fiber_generation/example/ot/doste_fibers.png new file mode 100644 index 000000000..b4cbdaa3e --- /dev/null +++ b/utilities/fiber_generation/example/ot/doste_fibers.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1b7dea227b8a8ecd5421fe295a8f7fc88f443310dd3b0a13f42d1623a1056eec +size 408585 diff --git a/utilities/fiber_generation/example/ot/doste_sheet.png b/utilities/fiber_generation/example/ot/doste_sheet.png new file mode 100644 index 000000000..158dbfb15 --- /dev/null +++ b/utilities/fiber_generation/example/ot/doste_sheet.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:d51341267ff44c808812bb9b195973c948cba68869f29d65727655647acd6cde +size 428858 diff --git a/utilities/fiber_generation/example/ot/doste_sheetnormal.png b/utilities/fiber_generation/example/ot/doste_sheetnormal.png new file mode 100644 index 000000000..bb00c2ee0 --- /dev/null +++ b/utilities/fiber_generation/example/ot/doste_sheetnormal.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:8623e33c660b6ff0463b7d17a18cf41362f72651eba1d73e484b4d09df64b53a +size 401352 diff --git a/utilities/fiber_generation/example/truncated/bayer.png b/utilities/fiber_generation/example/truncated/bayer.png deleted file mode 100644 index 73a017c6f..000000000 --- a/utilities/fiber_generation/example/truncated/bayer.png +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:c01ec2a92c75b6047fa9b6602924a619d1e20b577cc26de7852a91cc06d19b6a -size 313768 diff --git a/utilities/fiber_generation/example/truncated/bayer_fiber.png b/utilities/fiber_generation/example/truncated/bayer_fiber.png new file mode 100644 index 000000000..5953f168b --- /dev/null +++ b/utilities/fiber_generation/example/truncated/bayer_fiber.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:ffd63603ad6b0b53ea17dbfa6cdaedf866a6259b2011c77a9a8ba680e9ac73d7 +size 343872 diff --git a/utilities/fiber_generation/example/truncated/bayer_sheet.png b/utilities/fiber_generation/example/truncated/bayer_sheet.png new file mode 100644 index 000000000..6173257cb --- /dev/null +++ b/utilities/fiber_generation/example/truncated/bayer_sheet.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:490128920a15952150ba753f6c96fe401ee983a4c4fe153ad1d1f2a9022acc09 +size 368433 diff --git a/utilities/fiber_generation/example/truncated/bayer_sheetnormal.png b/utilities/fiber_generation/example/truncated/bayer_sheetnormal.png new file mode 100644 index 000000000..12c410d0e --- /dev/null +++ b/utilities/fiber_generation/example/truncated/bayer_sheetnormal.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:9ed589a695c29155c65aa2fbdb41501fd7fdffcab6a75deea9ace018a3ee0b8f +size 320327 From 6fd7d1585e914a2a4665193b25bcbbb68b221bdf Mon Sep 17 00:00:00 2001 From: Javiera Jilberto Vallejos Date: Mon, 29 Dec 2025 10:47:23 -0800 Subject: [PATCH 5/5] updating README.md --- utilities/fiber_generation/README.md | 11 +++++++++-- .../example/ot/{doste_fibers.png => doste_fiber.png} | 0 2 files changed, 9 insertions(+), 2 deletions(-) rename utilities/fiber_generation/example/ot/{doste_fibers.png => doste_fiber.png} (100%) diff --git a/utilities/fiber_generation/README.md b/utilities/fiber_generation/README.md index fad7dfa0b..0542f1284 100644 --- a/utilities/fiber_generation/README.md +++ b/utilities/fiber_generation/README.md @@ -6,8 +6,15 @@ Python + svMultiPhysics codes for fiber generation. Two methods are implemented: ### Examples The `main_bayer.py` and `main_doste.py` are scripts to run both methods in the geometry described in the `example/truncated` and `example/ot` folders respectively. -Results for truncated BiV (Bayer) -Results for BiV w/ outflow tracts (Doste) +#### Bayer results (fiber, sheet, sheet-normal) +Fiber direction for truncated BiV (Bayer) +Sheet direction for truncated BiV (Bayer) +Sheet-normal direction for truncated BiV (Bayer) + +#### Doste results (fiber, sheet, sheet-normal) +Fiber direction for BiV with OT (Doste) +Sheet direction for BiV with OT (Doste) +Sheet-normal direction for BiV with OT (Doste) Note that the Doste methods needs a geometry with outflow tracts to be run (each valve needs to be defined as a separated surface). Bayer can be run in any biventricular geometry. diff --git a/utilities/fiber_generation/example/ot/doste_fibers.png b/utilities/fiber_generation/example/ot/doste_fiber.png similarity index 100% rename from utilities/fiber_generation/example/ot/doste_fibers.png rename to utilities/fiber_generation/example/ot/doste_fiber.png