feat: new aha strain

doc/source/changelog/1212.added.md

@@ -0,0 +1 @@
+New aha strain

src/ansys/health/heart/post/auto_process.py

@@ -227,6 +227,13 @@ def mech_post(directory: str, model: HeartModel) -> None:
     out_dir = os.path.join(directory, "post", "lrc_strain")
     os.makedirs(out_dir, exist_ok=True)
     aha_strain = AhaStrainCalculator(model, d3plot_file=os.path.join(directory, "d3plot"))
-    aha_strain.compute_aha_strain(out_dir, write_vtk=True, t_to_keep=last_cycle_duration)
+    # aha_strain.compute_aha_strain(out_dir, write_vtk=True, t_to_keep=last_cycle_duration)
+
+    l_strain, c_strain = aha_strain.compute_longitudinal_radial_strain(
+        time_array=exporter.save_time, vtk_dir=out_dir
+    )
+
+    np.savetxt(os.path.join(out_dir, "longitudinal_strain.csv"), l_strain)
+    np.savetxt(os.path.join(out_dir, "circumferential_strain.csv"), c_strain)
 
     return

src/ansys/health/heart/post/strain_calculator.py

@@ -24,12 +24,19 @@
 
 import pathlib
 
+from deprecated import deprecated
 import numpy as np
 import pyvista as pv
 
+from ansys.health.heart import LOG as LOGGER
 from ansys.health.heart.models import BiVentricle, FourChamber, FullHeart, HeartModel, LeftVentricle
 from ansys.health.heart.post.dpf_utils import D3plotReader
-from ansys.health.heart.utils.landmark_utils import compute_aha17, compute_element_cs
+from ansys.health.heart.utils.landmark_utils import (
+    compute_aha17,
+    compute_aha17_lines,
+    compute_aha17_points,
+    compute_element_cs,
+)
 from ansys.health.heart.utils.vtk_utils import find_corresponding_points, generate_thickness_lines
 
 
@@ -54,6 +61,140 @@ def __init__(self, model: HeartModel, d3plot_file):
 
         self.d3plot = D3plotReader(d3plot_file)
 
+    def compute_longitudinal_radial_strain(
+        self, time_array: np.ndarray | list = None, vtk_dir: pathlib.Path | str = None
+    ) -> tuple[np.ndarray, np.ndarray]:
+        """
+        Compute longitudinal and radial strain.
+
+        Parameters
+        ----------
+        time_array: np.ndarray | list, optional
+            Array of time points to compute strain at. If None, use all time points.
+        vtk_dir: pathlib.Path | str, optional
+            Directory to save VTK files. If None, do not save VTK files.
+
+        Returns
+        -------
+        tuple[np.ndarray, np.ndarray]
+            Longitudinal and radial strain arrays of 17 segments.
+        """
+        if time_array is None:
+            time_array = self.d3plot.time
+
+        surface_endo: pv.PolyData = self.model.left_ventricle.endocardium.copy()
+
+        # Find AHA points on the endocardium surface.
+        points = compute_aha17_points(
+            self.model, self.model.short_axis, self.model.l2cv_axis, surface=surface_endo
+        )
+        ids = [surface_endo.find_closest_point(p) for p in points]
+
+        hlength_array = []
+        vlength_array = []
+        for it, t in enumerate(time_array):
+            coordinates = (
+                self.d3plot.model.results.coordinates.on_time_scoping(float(t)).eval()[0].data
+            )
+            surface_endo.points = coordinates[surface_endo["_global-point-ids"]]
+            try:
+                # points = compute_aha17_points(
+                #     self.model, self.model.short_axis, self.model.l2cv_axis, surface=surface_endo
+                # )
+                points = [surface_endo.points[id] for id in ids]
+                hlines, vlines = compute_aha17_lines(surface_endo, points)
+
+            except ValueError as e:
+                LOGGER.error(f"Error in computing AHA strain at time {t}: {e}")
+                exit()
+
+            hlength_array.append([hl.length for hl in hlines])
+            vlength_array.append([vl.length for vl in vlines])
+
+            if vtk_dir is not None:
+                vtk_dir = pathlib.Path(vtk_dir)
+                surface_endo.save(vtk_dir / f"endo_{it}.vtp")
+                pv.merge(hlines).save(vtk_dir / f"hlines_{it}.vtp")
+                pv.merge(vlines).save(vtk_dir / f"vlines_{it}.vtp")
+
+        hlength_array = np.array(hlength_array)
+        vlength_array = np.array(vlength_array)
+
+        # length variation of horizontal and vertical lines
+        # compared to the first image
+        h_strain = (hlength_array - hlength_array[0]) / hlength_array[0]
+        v_strain = (vlength_array - vlength_array[0]) / vlength_array[0]
+
+        """
+        save strain values with the following format:
+
+        P1---H1---P2---H2---P3----H3--P4---H4---P5---H5---P6---H6---P1
+        |          |         |         |         |         |         |
+        |          |         |         |         |         |         |
+        V1   S4   V2    S5  V3    S6  V4    S1  V5   S2   V6   S3   V1
+        |          |         |         |         |         |         |
+        |          |         |         |         |         |         |
+        P7---H7---P8---H8---P9----H9--P10--H10--P11--H11--P12--H12--P7
+        |          |         |         |         |         |         |
+        |          |         |         |         |         |         |
+        V7  S10   V8   S11  V9   S12  V10  S7   V11  S8   V12  S9   V7
+        |          |         |         |         |         |         |
+        |          |         |         |         |         |         |
+        P13--H13--P14--H14--P15--H15--P16--H16--P17--H17--P18--H18--P13
+            P19---H19---P20---H20---P21---H21---P22---H22---P19
+            |           |           |           |           |
+            |           |           |           |           |
+            V13  S15   V14   S16   V15   S13   V16   S14   V13
+            |           |           |           |           |
+            |           |           |           |           |
+            P23---H23---P24---H24---P25---H25---P26---H26---P23
+        """
+        # longitudinal strain
+        # seg 17 is always 0
+        aha_l_strain = np.zeros((len(time_array), 17))
+        aha_l_strain[:, 0] = 0.5 * (v_strain[:, 3] + v_strain[:, 4])
+        aha_l_strain[:, 1] = 0.5 * (v_strain[:, 4] + v_strain[:, 5])
+        aha_l_strain[:, 2] = 0.5 * (v_strain[:, 5] + v_strain[:, 0])
+        aha_l_strain[:, 3] = 0.5 * (v_strain[:, 0] + v_strain[:, 1])
+        aha_l_strain[:, 4] = 0.5 * (v_strain[:, 1] + v_strain[:, 2])
+        aha_l_strain[:, 5] = 0.5 * (v_strain[:, 2] + v_strain[:, 3])
+
+        aha_l_strain[:, 6] = 0.5 * (v_strain[:, 9] + v_strain[:, 10])
+        aha_l_strain[:, 7] = 0.5 * (v_strain[:, 10] + v_strain[:, 11])
+        aha_l_strain[:, 8] = 0.5 * (v_strain[:, 11] + v_strain[:, 6])
+        aha_l_strain[:, 9] = 0.5 * (v_strain[:, 6] + v_strain[:, 7])
+        aha_l_strain[:, 10] = 0.5 * (v_strain[:, 7] + v_strain[:, 8])
+        aha_l_strain[:, 11] = 0.5 * (v_strain[:, 8] + v_strain[:, 9])
+
+        aha_l_strain[:, 12] = 0.5 * (v_strain[:, 14] + v_strain[:, 15])
+        aha_l_strain[:, 13] = 0.5 * (v_strain[:, 15] + v_strain[:, 12])
+        aha_l_strain[:, 14] = 0.5 * (v_strain[:, 12] + v_strain[:, 13])
+        aha_l_strain[:, 15] = 0.5 * (v_strain[:, 13] + v_strain[:, 14])
+
+        # circumferential strain
+        # seg 17 is always 0
+        aha_c_strain = np.zeros((len(time_array), 17))
+        aha_c_strain[:, 0] = 0.5 * (h_strain[:, 3] + h_strain[:, 9])
+        aha_c_strain[:, 1] = 0.5 * (h_strain[:, 4] + h_strain[:, 10])
+        aha_c_strain[:, 2] = 0.5 * (h_strain[:, 5] + h_strain[:, 11])
+        aha_c_strain[:, 3] = 0.5 * (h_strain[:, 0] + h_strain[:, 6])
+        aha_c_strain[:, 4] = 0.5 * (h_strain[:, 1] + h_strain[:, 7])
+        aha_c_strain[:, 5] = 0.5 * (h_strain[:, 2] + h_strain[:, 8])
+
+        aha_c_strain[:, 6] = 0.5 * (h_strain[:, 9] + h_strain[:, 15])
+        aha_c_strain[:, 7] = 0.5 * (h_strain[:, 10] + h_strain[:, 16])
+        aha_c_strain[:, 8] = 0.5 * (h_strain[:, 11] + h_strain[:, 17])
+        aha_c_strain[:, 9] = 0.5 * (h_strain[:, 6] + h_strain[:, 12])
+        aha_c_strain[:, 10] = 0.5 * (h_strain[:, 7] + h_strain[:, 13])
+        aha_c_strain[:, 11] = 0.5 * (h_strain[:, 8] + h_strain[:, 14])
+
+        aha_c_strain[:, 12] = 0.5 * (h_strain[:, 20] + h_strain[:, 24])
+        aha_c_strain[:, 13] = 0.5 * (h_strain[:, 21] + h_strain[:, 25])
+        aha_c_strain[:, 14] = 0.5 * (h_strain[:, 18] + h_strain[:, 22])
+        aha_c_strain[:, 15] = 0.5 * (h_strain[:, 19] + h_strain[:, 23])
+
+        return aha_l_strain, aha_c_strain
+
     def _compute_thickness_lines(self, time_array: np.ndarray | list = None) -> list[pv.PolyData]:
         """Compute ventricular myocardium thickness.
 
@@ -144,6 +285,10 @@ def _compute_thickness(
             res.append(thickness_lines)
         return res
 
+    @deprecated(
+        reason="""This method will be deprecated in the future. Use
+                "`compute_longitudinal_radial_strain` instead."""
+    )
     def compute_aha_strain(
         self, out_dir: str = None, write_vtk: bool = False, t_to_keep: float = 10e10
     ) -> np.ndarray:
@@ -195,6 +340,10 @@ def compute_aha_strain(
 
         return strain
 
+    @deprecated(
+        reason="""This method will be deprecated in the future. Use
+                "`compute_longitudinal_radial_strain` instead."""
+    )
     def compute_aha_strain_at(self, frame: int = 0, out_dir: pathlib.Path = None) -> np.ndarray:
         """
         Export AHA strain and/or save a VTK file for a given frame.
@@ -232,6 +381,10 @@ def compute_aha_strain_at(self, frame: int = 0, out_dir: pathlib.Path = None) ->
 
         return aha_lrc
 
+    @deprecated(
+        reason="""This method will be deprecated in the future. Use
+                "`compute_longitudinal_radial_strain` instead."""
+    )
     def _compute_myocardial_strain(
         self, at_frame, reference=None
     ) -> tuple[list[float], list[float], list[float]]: