From b18aea8b953949f7a277eb8b1139b30146e6d1a3 Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Tue, 24 Oct 2023 20:58:53 +0200 Subject: [PATCH 1/2] Implement 17 segment aha markers --- src/cardiac_geometries/__init__.py | 5 + src/cardiac_geometries/_mesh.py | 25 +++++ src/cardiac_geometries/aha.py | 154 +++++++++++++++++++++++++++++ src/cardiac_geometries/cli.py | 10 ++ src/cardiac_geometries/geometry.py | 14 ++- tests/test_cli.py | 10 ++ 6 files changed, 217 insertions(+), 1 deletion(-) create mode 100644 src/cardiac_geometries/aha.py diff --git a/src/cardiac_geometries/__init__.py b/src/cardiac_geometries/__init__.py index 53a8914..10fde7b 100644 --- a/src/cardiac_geometries/__init__.py +++ b/src/cardiac_geometries/__init__.py @@ -17,6 +17,7 @@ from .fibers import _slab as slab_fibers from .fibers import _lv_ellipsoid as lv_ellipsoid_fibers from .fibers import _biv_ellipsoid as biv_ellipsoid_fibers + from .aha import lv_aha, biv_aha from . import _mesh as mesh from ._mesh import ( create_biv_ellipsoid, @@ -33,6 +34,8 @@ create_lv_ellipsoid = None # type: ignore create_slab = None # type: ignore create_biv_ellipsoid_torso = None # type: ignore + lv_aha = None # type: ignore + biv_aha = None # type: ignore if has_mshr(): from . import _mshr as mshr @@ -57,4 +60,6 @@ "create_lv_ellipsoid", "create_slab", "mesh", + "lv_aha", + "biv_aha", ] diff --git a/src/cardiac_geometries/_mesh.py b/src/cardiac_geometries/_mesh.py index a686d06..88faad5 100644 --- a/src/cardiac_geometries/_mesh.py +++ b/src/cardiac_geometries/_mesh.py @@ -427,6 +427,7 @@ def create_lv_ellipsoid( fiber_angle_endo: float = -60, fiber_angle_epi: float = +60, fiber_space: str = "P_1", + aha: bool = True, ) -> Optional[Geometry]: """Create an LV ellipsoidal geometry @@ -461,6 +462,8 @@ def create_lv_ellipsoid( Angle for the epicardium, by default +60 fiber_space : str, optional Function space for fibers of the form family_degree, by default "P_1" + aha : bool, optional + If True create 17-segment AHA regions Returns ------- @@ -501,6 +504,7 @@ def create_lv_ellipsoid( "fibers_angle_endo": fiber_angle_endo, "fibers_angle_epi": fiber_angle_epi, "fiber_space": fiber_space, + "aha": aha, "mesh_type": MeshTypes.lv_ellipsoid.value, "cardiac_geometry_version": __version__, "timestamp": datetime.datetime.now().isoformat(), @@ -533,6 +537,20 @@ def create_lv_ellipsoid( geometry = gmsh2dolfin(mesh_name, unlink=False) + if aha: + from .aha import lv_aha + + geometry = lv_aha( + geometry=geometry, + r_long_endo=r_long_endo, + r_short_endo=r_short_endo, + mu_base=mu_base_endo, + ) + from dolfin import XDMFFile + + with XDMFFile((outdir / "cfun.xdmf").as_posix()) as xdmf: + xdmf.write(geometry.marker_functions.cfun) + with open(outdir / "markers.json", "w") as f: json.dump(geometry.markers, f, default=json_serial) @@ -554,6 +572,13 @@ def create_lv_ellipsoid( ) geo = Geometry.from_folder(outdir) + if aha: + # Update schema + from .geometry import H5Path + + cfun = geo.schema["cfun"].to_dict() + cfun["fname"] = "cfun.xdmf:f" + geo.schema["cfun"] = H5Path(**cfun) if _tmpfile is not None: _tmpfile.__exit__(None, None, None) diff --git a/src/cardiac_geometries/aha.py b/src/cardiac_geometries/aha.py new file mode 100644 index 0000000..9626871 --- /dev/null +++ b/src/cardiac_geometries/aha.py @@ -0,0 +1,154 @@ +import dolfin +import numpy as np + +from .dolfin_utils import Geometry + + +def focal(r_long_endo: float, r_short_endo: float): + return np.sqrt(r_long_endo**2 - r_short_endo**2) + + +def full_arctangent(x, y): + t = np.arctan2(x, y) + if t < 0: + return t + 2 * np.pi + else: + return t + + +def cartesian_to_prolate_ellipsoidal(x, y, z, a): + b1 = np.sqrt((x + a) ** 2 + y**2 + z**2) + b2 = np.sqrt((x - a) ** 2 + y**2 + z**2) + + sigma = 1 / (2.0 * a) * (b1 + b2) + tau = 1 / (2.0 * a) * (b1 - b2) + phi = full_arctangent(z, y) + nu = np.arccosh(sigma) + mu = np.arccos(tau) + return nu, mu, phi + + +def get_level(regions, mu): + A = np.intersect1d( + np.where((regions.T[3] <= mu))[0], np.where((mu <= regions.T[0]))[0] + ) + if len(A) == 0: + return [np.shape(regions)[0] + 1] + else: + return A + + +def get_sector(regions, theta): + if not ( + np.count_nonzero(regions.T[1] <= regions.T[2]) >= 0.5 * np.shape(regions)[0] + ): + raise ValueError("Surfaces are flipped") + + sectors = [] + for i, r in enumerate(regions): + if r[1] == r[2]: + sectors.append(i) + else: + if r[1] > r[2]: + if theta > r[1] or r[2] > theta: + sectors.append(i) + + else: + if r[1] < theta < r[2]: + sectors.append(i) + + return sectors + + +class AHA_LV_17(dolfin.UserExpression): + def __init__(self, foc: float, mu_base: float) -> None: + super().__init__() + self.mu_base = abs(mu_base) + self.foc = foc + + @property + def dmu(self) -> float: + return (np.pi - self.mu_base) / 4 + + def eval_cell(self, value, x, ufc_cell): + nu, mu, phi = cartesian_to_prolate_ellipsoidal(*x, a=self.foc) + + if self.mu_base < mu <= self.mu_base + self.dmu: + # BASE + if 0 < phi <= np.pi / 3: + value[0] = 1 + elif np.pi / 3 < phi <= 2 * np.pi / 3: + value[0] = 2 + elif 2 * np.pi / 3 < phi <= np.pi: + value[0] = 3 + elif np.pi < phi <= 4 * np.pi / 3: + value[0] = 4 + elif 4 * np.pi / 3 < phi <= 5 * np.pi / 3: + value[0] = 5 + else: + value[0] = 6 + + elif self.mu_base + self.dmu < mu <= self.mu_base + 2 * self.dmu: + # MID + if 0 < phi <= np.pi / 3: + value[0] = 7 + elif np.pi / 3 < phi <= 2 * np.pi / 3: + value[0] = 8 + elif 2 * np.pi / 3 < phi <= np.pi: + value[0] = 9 + elif np.pi < phi <= 4 * np.pi / 3: + value[0] = 10 + elif 4 * np.pi / 3 < phi <= 5 * np.pi / 3: + value[0] = 11 + else: + value[0] = 12 + elif self.mu_base + 2 * self.dmu < mu <= self.mu_base + 3 * self.dmu: + # APICAL + if 0 < phi <= np.pi / 2: + value[0] = 13 + elif np.pi / 2 < phi <= np.pi: + value[0] = 14 + elif np.pi < phi <= 3 * np.pi / 2: + value[0] = 15 + else: + value[0] = 16 + else: + value[0] = 17 + + def value_shape(self): + return () + + +def lv_aha( + geometry: Geometry, r_long_endo: float, r_short_endo: float, mu_base: float +) -> Geometry: + foc = focal(r_long_endo=r_long_endo, r_short_endo=r_short_endo) + + V = dolfin.FunctionSpace(geometry.mesh, "DG", 0) + # f = dolfin.Function(V) + expr = AHA_LV_17(foc=foc, mu_base=mu_base) + f = dolfin.interpolate(expr, V) + geometry.marker_functions.cfun.array()[:] = f.vector().get_local() + i = 1 + + for level in ["BASAL", "MID", "APICAL"]: + for sector in [ + "ANTERIOR", + "ANTEROSEPTAL", + "SEPTAL", + "INFERIOR", + "POSTERIOR", + "LATERAL", + ]: + if level == "APICAL" and sector in ("ANTEROSEPTAL", "POSTERIOR"): + continue + + geometry.markers["-".join((level, sector))] = (i, 3) + i += 1 + geometry.markers["APEX"] = (17, 3) + geometry.markers.pop("MYOCARDIUM") + return geometry + + +def biv_aha(): + ... diff --git a/src/cardiac_geometries/cli.py b/src/cardiac_geometries/cli.py index 2f0682d..7eb74d5 100644 --- a/src/cardiac_geometries/cli.py +++ b/src/cardiac_geometries/cli.py @@ -129,6 +129,14 @@ def app(): help="Function space for fibers of the form family_degree", show_default=True, ) +@click.option( + "--aha", + default=True, + is_flag=True, + type=bool, + help="If True create 17-segment AHA regions", + show_default=True, +) def create_lv_ellipsoid( outdir: Path, r_short_endo: float = 7.0, @@ -144,6 +152,7 @@ def create_lv_ellipsoid( fiber_angle_endo: float = -60, fiber_angle_epi: float = +60, fiber_space: str = "P_1", + aha: bool = True, ): outdir = Path(outdir) outdir.mkdir(exist_ok=True) @@ -165,6 +174,7 @@ def create_lv_ellipsoid( fiber_angle_endo=fiber_angle_endo, fiber_angle_epi=fiber_angle_epi, fiber_space=fiber_space, + aha=aha, ) if geo is not None: geo.save(outdir / "lv_ellipsoid.h5") diff --git a/src/cardiac_geometries/geometry.py b/src/cardiac_geometries/geometry.py index 0a3d72a..2973a2c 100644 --- a/src/cardiac_geometries/geometry.py +++ b/src/cardiac_geometries/geometry.py @@ -466,12 +466,24 @@ def from_folder(cls, folder, schema: Optional[Dict[str, H5Path]] = None): folder = Path(folder) if schema is None: - schema = cls.default_schema() + # First look for a schema file + for f in folder.iterdir(): + if f.suffix == ".json": + try: + schema = load_schema(f) + except Exception: + continue + else: + break + + else: + schema = cls.default_schema() # Load mesh first data = {} for name, p in schema.items(): + print(name, p) if p.fname == "": continue if not p.is_mesh: diff --git a/tests/test_cli.py b/tests/test_cli.py index ba7312d..f60feb2 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -46,6 +46,16 @@ def test_create_lv_ellipsoid(tmp_path: Path): geo = Geometry.from_file(outfile) assert geo.f0 is not None + # Make sure all aha regions have nonzero volume + import dolfin + + dx = dolfin.dx(domain=geo.mesh, subdomain_data=geo.cfun) + for k, v in geo.markers.items(): + if v[1] != 3: + continue + vol = dolfin.assemble(dolfin.Constant(1.0) * dx(v[0])) + assert vol > 0, k + def test_create_biv_ellipsoid(tmp_path: Path): runner = CliRunner() From 81a4100307f0bb9cbd46db1b5b59f62255ad55cb Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Tue, 24 Oct 2023 21:06:59 +0200 Subject: [PATCH 2/2] Implement 17 segment aha markers --- src/cardiac_geometries/aha.py | 8 ++++++-- src/cardiac_geometries/cli.py | 2 +- src/cardiac_geometries/geometry.py | 33 +++++++++++++++++------------- 3 files changed, 26 insertions(+), 17 deletions(-) diff --git a/src/cardiac_geometries/aha.py b/src/cardiac_geometries/aha.py index 9626871..401acbf 100644 --- a/src/cardiac_geometries/aha.py +++ b/src/cardiac_geometries/aha.py @@ -30,7 +30,8 @@ def cartesian_to_prolate_ellipsoidal(x, y, z, a): def get_level(regions, mu): A = np.intersect1d( - np.where((regions.T[3] <= mu))[0], np.where((mu <= regions.T[0]))[0] + np.where((regions.T[3] <= mu))[0], + np.where((mu <= regions.T[0]))[0], ) if len(A) == 0: return [np.shape(regions)[0] + 1] @@ -120,7 +121,10 @@ def value_shape(self): def lv_aha( - geometry: Geometry, r_long_endo: float, r_short_endo: float, mu_base: float + geometry: Geometry, + r_long_endo: float, + r_short_endo: float, + mu_base: float, ) -> Geometry: foc = focal(r_long_endo=r_long_endo, r_short_endo=r_short_endo) diff --git a/src/cardiac_geometries/cli.py b/src/cardiac_geometries/cli.py index 7eb74d5..051c78d 100644 --- a/src/cardiac_geometries/cli.py +++ b/src/cardiac_geometries/cli.py @@ -130,7 +130,7 @@ def app(): show_default=True, ) @click.option( - "--aha", + "--aha/--no-aha", default=True, is_flag=True, type=bool, diff --git a/src/cardiac_geometries/geometry.py b/src/cardiac_geometries/geometry.py index 2973a2c..51e6235 100644 --- a/src/cardiac_geometries/geometry.py +++ b/src/cardiac_geometries/geometry.py @@ -112,6 +112,19 @@ def to_dict(self): return self._asdict() +def find_schema_in_folder(folder: Path) -> Optional[Dict[str, H5Path]]: + # First look for a schema file + for f in folder.iterdir(): + if f.suffix == ".json": + try: + schema = load_schema(f) + except Exception: + pass + else: + return schema + return None + + def load_schema(path: Path) -> Optional[Dict[str, H5Path]]: if not path.is_file(): return None @@ -174,7 +187,10 @@ def read( raise RuntimeError(f"Unknown file format for {fname}") -def extract_fname_group(fname: str, folder=".") -> Tuple[Path, Optional[str]]: +def extract_fname_group( + fname: str, + folder: Union[Path, str] = ".", +) -> Tuple[Path, Optional[str]]: fg = fname.split(":") if len(fg) == 1: return Path(folder) / fg[0], None @@ -466,24 +482,13 @@ def from_folder(cls, folder, schema: Optional[Dict[str, H5Path]] = None): folder = Path(folder) if schema is None: - # First look for a schema file - for f in folder.iterdir(): - if f.suffix == ".json": - try: - schema = load_schema(f) - except Exception: - continue - else: - break - - else: + if (schema := find_schema_in_folder(folder)) is None: schema = cls.default_schema() - + assert schema is not None # Load mesh first data = {} for name, p in schema.items(): - print(name, p) if p.fname == "": continue if not p.is_mesh: