Merge pull request #18 from ComputationalPhysiology/aha-regions

Aha regions

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",
 ]

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)

src/cardiac_geometries/aha.py

@@ -0,0 +1,158 @@
+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():
+    ...

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/--no-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")

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,8 +482,9 @@ def from_folder(cls, folder, schema: Optional[Dict[str, H5Path]] = None):
         folder = Path(folder)
 
         if schema is None:
-            schema = cls.default_schema()
-
+            if (schema := find_schema_in_folder(folder)) is None:
+                schema = cls.default_schema()
+        assert schema is not None
         # Load mesh first
         data = {}
 

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