From 3613dbcf14fec070da24f095690b35b37acce7a6 Mon Sep 17 00:00:00 2001
From: wenfengye <wenfeng.ye@ansys.com>
Date: Wed, 13 Aug 2025 13:19:47 +0200
Subject: [PATCH 01/14] add method to extract aha lines

---
 .../health/heart/utils/landmark_utils.py      | 200 ++++++++++++++++++
 tests/heart/utils/test_landmarks.py           |  21 ++
 2 files changed, 221 insertions(+)

diff --git a/src/ansys/health/heart/utils/landmark_utils.py b/src/ansys/health/heart/utils/landmark_utils.py
index 31ad2844c..9229b1066 100644
--- a/src/ansys/health/heart/utils/landmark_utils.py
+++ b/src/ansys/health/heart/utils/landmark_utils.py
@@ -26,6 +26,7 @@
 
 from deprecated import deprecated
 import numpy as np
+import pyvista as pv
 from scipy.spatial.transform import Rotation
 
 from ansys.health.heart.models import HeartModel
@@ -79,6 +80,187 @@ def compute_anatomy_axis(
     return (l4cv_axis, l2cv_axis, short_axis)
 
 
+def compute_aha17_landmarks(
+    model: HeartModel,
+    short_axis: dict,
+    long_axis: dict,
+):
+    """_summary_.
+
+    P1--H1----P2---H2---P3----H3--P4---H4---P5---H5---P6---H6---P1
+    |          |         |         |         |         |         |
+    |          |         |         |         |         |         |
+    V1        V2        V3        V4        V5        V6        V1
+    |          |         |         |         |         |         |
+    |          |         |         |         |         |         |
+    P7---H7---P8---H8---P9----H9--P10--H10--P11--H11--P12--H12--P7
+    |          |         |         |         |         |         |
+    |          |         |         |         |         |         |
+    V7        V8        V9        V10       V11       V12       V7
+    |          |         |         |         |         |         |
+    |          |         |         |         |         |         |
+    P13--H13--P14--H14--P15--H15--P16--H16--P17--H17--P18--H18--P13
+          P19---H19---P20---H20---P21---H21---P22---H22---P19
+           |           |           |           |           |
+           |           |           |           |           |
+           V13        V14         V15         V16         V17
+           |           |           |           |           |
+           |           |           |           |           |
+          P23---H23---P24---H24---P25---H25---P26---H26---P23
+
+    Parameters
+    ----------
+    model : HeartModel
+        _description_.
+    short_axis : dict
+        _description_.
+    l4cv_axis : dict
+        _description_.
+
+    Returns
+    -------
+    np.ndarray
+        _description_.
+    """
+    p_basal, p_mid, p_apical, apex_endo, apex_epi = _calculate_longitudinal_points(
+        model, short_axis
+    )
+    axe_60, axe_120, axe_180, axe_45, axe_135 = _calculate_rotation_axis(short_axis, long_axis)
+
+    # Note: epicardium is not supported because it's incomplete for more than left ventricle model
+    surf: pv.PolyData = model.left_ventricle.endocardium
+
+    coords = []
+
+    for center in [p_basal, p_mid, p_apical]:
+        for normal in [-axe_60, -axe_120, axe_180, axe_60, axe_120, -axe_180]:
+            point, _ = surf.ray_trace(center, center + 1e6 * normal, first_point=True)
+            coords.append(point)
+    for center in [p_apical, apex_endo]:
+        for normal in [-axe_45, -axe_135, axe_45, axe_135]:
+            point, _ = surf.ray_trace(center, center + 1e6 * normal, first_point=True)
+            coords.append(point)
+
+    coords = np.array(coords)
+
+    hline = [
+        _project_line_segment(surf, p_basal, coords[0], coords[1]),
+        _project_line_segment(surf, p_basal, coords[1], coords[2]),
+        _project_line_segment(surf, p_basal, coords[2], coords[3]),
+        _project_line_segment(surf, p_basal, coords[3], coords[4]),
+        _project_line_segment(surf, p_basal, coords[4], coords[5]),
+        _project_line_segment(surf, p_basal, coords[5], coords[0]),
+        _project_line_segment(surf, p_mid, coords[6], coords[7]),
+        _project_line_segment(surf, p_mid, coords[7], coords[8]),
+        _project_line_segment(surf, p_mid, coords[8], coords[9]),
+        _project_line_segment(surf, p_mid, coords[9], coords[10]),
+        _project_line_segment(surf, p_mid, coords[10], coords[11]),
+        _project_line_segment(surf, p_mid, coords[11], coords[6]),
+        _project_line_segment(surf, p_apical, coords[12], coords[13]),
+        _project_line_segment(surf, p_apical, coords[13], coords[14]),
+        _project_line_segment(surf, p_apical, coords[14], coords[15]),
+        _project_line_segment(surf, p_apical, coords[15], coords[16]),
+        _project_line_segment(surf, p_apical, coords[16], coords[17]),
+        _project_line_segment(surf, p_apical, coords[17], coords[12]),
+        _project_line_segment(surf, p_apical, coords[18], coords[19]),
+        _project_line_segment(surf, p_apical, coords[19], coords[20]),
+        _project_line_segment(surf, p_apical, coords[20], coords[21]),
+        _project_line_segment(surf, p_apical, coords[21], coords[18]),
+        _project_line_segment(surf, apex_endo, coords[22], coords[23]),
+        _project_line_segment(surf, apex_endo, coords[23], coords[24]),
+        _project_line_segment(surf, apex_endo, coords[24], coords[25]),
+        _project_line_segment(surf, apex_endo, coords[25], coords[22]),
+    ]
+
+    hlines = pv.MultiBlock(hline)
+    hlines.save("hlines.vtm")
+
+    vline = [
+        _project_line_segment(surf, p_basal, coords[0], coords[6]),
+        _project_line_segment(surf, p_basal, coords[1], coords[7]),
+        _project_line_segment(surf, p_basal, coords[2], coords[8]),
+        _project_line_segment(surf, p_basal, coords[3], coords[9]),
+        _project_line_segment(surf, p_basal, coords[4], coords[10]),
+        _project_line_segment(surf, p_basal, coords[5], coords[11]),
+        _project_line_segment(surf, p_mid, coords[6], coords[12]),
+        _project_line_segment(surf, p_mid, coords[7], coords[13]),
+        _project_line_segment(surf, p_mid, coords[8], coords[14]),
+        _project_line_segment(surf, p_mid, coords[9], coords[15]),
+        _project_line_segment(surf, p_mid, coords[10], coords[16]),
+        _project_line_segment(surf, p_mid, coords[11], coords[17]),
+        _project_line_segment(surf, p_apical, coords[18], coords[22]),
+        _project_line_segment(surf, p_apical, coords[19], coords[23]),
+        _project_line_segment(surf, p_apical, coords[20], coords[24]),
+        _project_line_segment(surf, p_apical, coords[21], coords[25]),
+    ]
+
+    vlines = pv.MultiBlock()
+    for h in vline:
+        vlines.append(h)
+    vlines.save("vlines.vtm")
+    return surf, coords
+
+
+def _project_line_segment(
+    surf: pv.PolyData, center: np.ndarray, p1: np.ndarray, p2: np.ndarray
+) -> pv.PolyData:
+    """Project a line segment onto a surface.
+
+    Parameters
+    ----------
+    surf : pv.PolyData
+        The surface to project onto.
+    center : np.ndarray
+        The center point of the projection.
+    p1 : np.ndarray
+        The first point of the line segment.
+    p2 : np.ndarray
+        The second point of the line segment.
+
+    Returns
+    -------
+    pv.PolyData
+        The projected line segment.
+    """
+    segment_points = np.linspace(p1, p2, num=100)
+
+    # Project each point
+    projected_points = []
+    for pt in segment_points:
+        start = center
+        end = center + 1e9 * (pt - center)
+        intersection = surf.ray_trace(start, end, first_point=True)[0]
+        if intersection.size > 0:
+            projected_points.append(intersection)
+
+    # Convert to array
+    projected_points = np.array(projected_points)
+
+    # Create curve
+    projected_line = pv.Spline(projected_points, n_points=len(projected_points))
+
+    return projected_line
+
+
+def _calculate_longitudinal_points(model, short_axis):
+    """Define landmarks along the short axis."""
+    for apex in model.left_ventricle.apex_points:
+        if "endocardium" in apex.name:
+            apex_ed = apex.xyz
+        elif "epicardium" in apex.name:
+            apex_ep = apex.xyz
+
+    p_basal = short_axis["center"]
+    p_mid = 1 / 3 * (apex_ep - p_basal) + p_basal
+    p_apical = 2 / 3 * (apex_ep - p_basal) + p_basal
+
+    # to have a flat segment 17, project endocardical apex point on short axis
+    x = apex_ed - apex_ep
+    y = p_basal - apex_ep
+    apex_ed = y * np.dot(x, y) / np.dot(y, y) + apex_ep
+    return p_basal, p_mid, p_apical, apex_ed, apex_ep
+
+
 def compute_aha17(
     model: HeartModel,
     short_axis: dict,
@@ -237,6 +419,24 @@ def compute_aha17(
     return aha_ids
 
 
+def _calculate_rotation_axis(short_axis, l4cv_axis):
+    short_normal = short_axis["normal"]
+
+    # define reference cut plane
+    # default: rotate 60 from long axis
+    long_axis = l4cv_axis["normal"]
+    axe_60 = Rotation.from_rotvec(np.radians(60) * short_normal).apply(  # noqa:E501
+        long_axis
+    )
+
+    axe_120 = Rotation.from_rotvec(np.radians(60) * short_normal).apply(axe_60)
+    axe_180 = -Rotation.from_rotvec(np.radians(60) * short_normal).apply(axe_120)
+    axe_45 = Rotation.from_rotvec(np.radians(-15) * short_normal).apply(axe_60)
+    axe_135 = Rotation.from_rotvec(np.radians(90) * short_normal).apply(axe_45)
+
+    return axe_60, axe_120, axe_180, axe_45, axe_135
+
+
 @deprecated(reason="Using gradient from UVC to get better results.")
 def compute_element_cs(
     model: HeartModel, short_axis: dict, aha_element: np.ndarray
diff --git a/tests/heart/utils/test_landmarks.py b/tests/heart/utils/test_landmarks.py
index d4f5fadc0..fc8688cf6 100644
--- a/tests/heart/utils/test_landmarks.py
+++ b/tests/heart/utils/test_landmarks.py
@@ -28,6 +28,7 @@
 import ansys.health.heart.models as models
 from ansys.health.heart.utils.landmark_utils import (
     compute_aha17,
+    compute_aha17_landmarks,
     compute_anatomy_axis,
     compute_element_cs,
 )
@@ -91,6 +92,26 @@ def test_compute_aha17(model):
     assert np.sum(np.isnan(aha_ids)) == 361818
 
 
+def test_compute_aha17_landmarks(model):
+    l4cv = {
+        "center": np.array([30.05990578, 109.49603257, 375.8178529]),
+        "normal": np.array([-0.54847906, -0.10082087, -0.83006378]),
+    }
+    l4cv = {
+        "center": np.array([42.8065945, 105.236255, 367.91265619]),
+        "normal": np.array([0.57496125, 0.67530147, -0.46193883]),
+    }
+    short = {
+        "center": np.array([26.07305844, 125.37379059, 376.52369382]),
+        "normal": np.array([0.60711636, -0.73061828, -0.31242063]),
+    }
+    surf, coords = compute_aha17_landmarks(model, short, l4cv)
+
+    surf.save("surf.vtk")
+
+    np.savetxt("points.txt", coords)
+
+
 def test_compute_element_cs(model):
     l4cv = {
         "center": np.array([30.05990578, 109.49603257, 375.8178529]),

From 6248cf7513dc03b4c1865a98b71ed1c0a8088490 Mon Sep 17 00:00:00 2001
From: wenfengye <wenfeng.ye@ansys.com>
Date: Wed, 13 Aug 2025 13:24:02 +0200
Subject: [PATCH 02/14] clean

---
 src/ansys/health/heart/utils/landmark_utils.py | 17 ++++++-----------
 tests/heart/utils/test_landmarks.py            | 11 ++++++++---
 2 files changed, 14 insertions(+), 14 deletions(-)

diff --git a/src/ansys/health/heart/utils/landmark_utils.py b/src/ansys/health/heart/utils/landmark_utils.py
index 9229b1066..ffbbb7f43 100644
--- a/src/ansys/health/heart/utils/landmark_utils.py
+++ b/src/ansys/health/heart/utils/landmark_utils.py
@@ -122,14 +122,16 @@ def compute_aha17_landmarks(
     np.ndarray
         _description_.
     """
+    # Note: epicardium is not supported because it's incomplete for more than left ventricle model
+    surf: pv.PolyData = model.left_ventricle.endocardium
+
+    # get anatomical points
     p_basal, p_mid, p_apical, apex_endo, apex_epi = _calculate_longitudinal_points(
         model, short_axis
     )
     axe_60, axe_120, axe_180, axe_45, axe_135 = _calculate_rotation_axis(short_axis, long_axis)
 
-    # Note: epicardium is not supported because it's incomplete for more than left ventricle model
-    surf: pv.PolyData = model.left_ventricle.endocardium
-
+    # define points for the segments
     coords = []
 
     for center in [p_basal, p_mid, p_apical]:
@@ -172,9 +174,6 @@ def compute_aha17_landmarks(
         _project_line_segment(surf, apex_endo, coords[25], coords[22]),
     ]
 
-    hlines = pv.MultiBlock(hline)
-    hlines.save("hlines.vtm")
-
     vline = [
         _project_line_segment(surf, p_basal, coords[0], coords[6]),
         _project_line_segment(surf, p_basal, coords[1], coords[7]),
@@ -194,11 +193,7 @@ def compute_aha17_landmarks(
         _project_line_segment(surf, p_apical, coords[21], coords[25]),
     ]
 
-    vlines = pv.MultiBlock()
-    for h in vline:
-        vlines.append(h)
-    vlines.save("vlines.vtm")
-    return surf, coords
+    return surf, coords, hline, vline
 
 
 def _project_line_segment(
diff --git a/tests/heart/utils/test_landmarks.py b/tests/heart/utils/test_landmarks.py
index fc8688cf6..4cac50733 100644
--- a/tests/heart/utils/test_landmarks.py
+++ b/tests/heart/utils/test_landmarks.py
@@ -24,6 +24,7 @@
 
 import numpy as np
 import pytest
+import pyvista as pv
 
 import ansys.health.heart.models as models
 from ansys.health.heart.utils.landmark_utils import (
@@ -105,11 +106,15 @@ def test_compute_aha17_landmarks(model):
         "center": np.array([26.07305844, 125.37379059, 376.52369382]),
         "normal": np.array([0.60711636, -0.73061828, -0.31242063]),
     }
-    surf, coords = compute_aha17_landmarks(model, short, l4cv)
+    surf, coords, hline, vline = compute_aha17_landmarks(model, short, l4cv)
 
-    surf.save("surf.vtk")
+    # surf.save("surf.vtk")
 
-    np.savetxt("points.txt", coords)
+    # np.savetxt("points.txt", coords)
+    hlines = pv.MultiBlock(hline)
+    hlines.save("hlines.vtm")
+    vlines = pv.MultiBlock(vline)
+    vlines.save("vlines.vtm")
 
 
 def test_compute_element_cs(model):

From 690892d81cf9aaeee9b0887de782bb117898c00c Mon Sep 17 00:00:00 2001
From: wenfengye <wenfeng.ye@ansys.com>
Date: Wed, 13 Aug 2025 14:50:58 +0200
Subject: [PATCH 03/14] add test

---
 .../health/heart/utils/landmark_utils.py      | 149 ++++++++++--------
 tests/heart/utils/test_landmarks.py           |  36 ++---
 2 files changed, 98 insertions(+), 87 deletions(-)

diff --git a/src/ansys/health/heart/utils/landmark_utils.py b/src/ansys/health/heart/utils/landmark_utils.py
index ffbbb7f43..400ea8dc9 100644
--- a/src/ansys/health/heart/utils/landmark_utils.py
+++ b/src/ansys/health/heart/utils/landmark_utils.py
@@ -81,29 +81,34 @@ def compute_anatomy_axis(
 
 
 def compute_aha17_landmarks(
-    model: HeartModel,
-    short_axis: dict,
-    long_axis: dict,
-):
-    """_summary_.
+    model: HeartModel, short_axis: dict, long_axis: dict, surface: pv.PolyData = None
+) -> tuple[list[np.ndarray], list[pv.PolyData], list[pv.PolyData]]:
+    """Compute AHA17 landmarks for the left ventricle.
+
+    The AHA17 landmarks are defined by
+    - 26 points
+    - 26 horizontal lines
+    - 16 vertical lines
 
-    P1--H1----P2---H2---P3----H3--P4---H4---P5---H5---P6---H6---P1
+    They split left ventricle into 17 segments as follow:
+
+    P1---H1---P2---H2---P3----H3--P4---H4---P5---H5---P6---H6---P1
     |          |         |         |         |         |         |
     |          |         |         |         |         |         |
-    V1        V2        V3        V4        V5        V6        V1
+    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        V8        V9        V10       V11       V12       V7
+    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        V14         V15         V16         V17
+           V13  S15   V14   S16   V15   S13   V16   S14   V13
            |           |           |           |           |
            |           |           |           |           |
           P23---H23---P24---H24---P25---H25---P26---H26---P23
@@ -111,19 +116,27 @@ def compute_aha17_landmarks(
     Parameters
     ----------
     model : HeartModel
-        _description_.
+        Heart model.
     short_axis : dict
-        _description_.
-    l4cv_axis : dict
-        _description_.
+        Short axis.
+    long_axis : dict
+        Long axis.
+    surface : pv.PolyData, default: None
+        Surface to be evaluated, only endocardium is supported.
+        If not given, the endocardium from heart model is used.
 
     Returns
     -------
-    np.ndarray
-        _description_.
+    tuple[np.ndarray, list[pv.PolyData], list[pv.PolyData]]
+        List of coordinates of the landmarks
+
+        List of horizontal lines
+
+        List of vertical lines
     """
     # Note: epicardium is not supported because it's incomplete for more than left ventricle model
-    surf: pv.PolyData = model.left_ventricle.endocardium
+    if surface is None:
+        surface: pv.PolyData = model.left_ventricle.endocardium
 
     # get anatomical points
     p_basal, p_mid, p_apical, apex_endo, apex_epi = _calculate_longitudinal_points(
@@ -132,68 +145,66 @@ def compute_aha17_landmarks(
     axe_60, axe_120, axe_180, axe_45, axe_135 = _calculate_rotation_axis(short_axis, long_axis)
 
     # define points for the segments
-    coords = []
+    points = []
 
     for center in [p_basal, p_mid, p_apical]:
         for normal in [-axe_60, -axe_120, axe_180, axe_60, axe_120, -axe_180]:
-            point, _ = surf.ray_trace(center, center + 1e6 * normal, first_point=True)
-            coords.append(point)
+            point, _ = surface.ray_trace(center, center + 1e6 * normal, first_point=True)
+            points.append(point)
     for center in [p_apical, apex_endo]:
         for normal in [-axe_45, -axe_135, axe_45, axe_135]:
-            point, _ = surf.ray_trace(center, center + 1e6 * normal, first_point=True)
-            coords.append(point)
-
-    coords = np.array(coords)
+            point, _ = surface.ray_trace(center, center + 1e6 * normal, first_point=True)
+            points.append(point)
 
     hline = [
-        _project_line_segment(surf, p_basal, coords[0], coords[1]),
-        _project_line_segment(surf, p_basal, coords[1], coords[2]),
-        _project_line_segment(surf, p_basal, coords[2], coords[3]),
-        _project_line_segment(surf, p_basal, coords[3], coords[4]),
-        _project_line_segment(surf, p_basal, coords[4], coords[5]),
-        _project_line_segment(surf, p_basal, coords[5], coords[0]),
-        _project_line_segment(surf, p_mid, coords[6], coords[7]),
-        _project_line_segment(surf, p_mid, coords[7], coords[8]),
-        _project_line_segment(surf, p_mid, coords[8], coords[9]),
-        _project_line_segment(surf, p_mid, coords[9], coords[10]),
-        _project_line_segment(surf, p_mid, coords[10], coords[11]),
-        _project_line_segment(surf, p_mid, coords[11], coords[6]),
-        _project_line_segment(surf, p_apical, coords[12], coords[13]),
-        _project_line_segment(surf, p_apical, coords[13], coords[14]),
-        _project_line_segment(surf, p_apical, coords[14], coords[15]),
-        _project_line_segment(surf, p_apical, coords[15], coords[16]),
-        _project_line_segment(surf, p_apical, coords[16], coords[17]),
-        _project_line_segment(surf, p_apical, coords[17], coords[12]),
-        _project_line_segment(surf, p_apical, coords[18], coords[19]),
-        _project_line_segment(surf, p_apical, coords[19], coords[20]),
-        _project_line_segment(surf, p_apical, coords[20], coords[21]),
-        _project_line_segment(surf, p_apical, coords[21], coords[18]),
-        _project_line_segment(surf, apex_endo, coords[22], coords[23]),
-        _project_line_segment(surf, apex_endo, coords[23], coords[24]),
-        _project_line_segment(surf, apex_endo, coords[24], coords[25]),
-        _project_line_segment(surf, apex_endo, coords[25], coords[22]),
+        _project_line_segment(surface, p_basal, points[0], points[1]),
+        _project_line_segment(surface, p_basal, points[1], points[2]),
+        _project_line_segment(surface, p_basal, points[2], points[3]),
+        _project_line_segment(surface, p_basal, points[3], points[4]),
+        _project_line_segment(surface, p_basal, points[4], points[5]),
+        _project_line_segment(surface, p_basal, points[5], points[0]),
+        _project_line_segment(surface, p_mid, points[6], points[7]),
+        _project_line_segment(surface, p_mid, points[7], points[8]),
+        _project_line_segment(surface, p_mid, points[8], points[9]),
+        _project_line_segment(surface, p_mid, points[9], points[10]),
+        _project_line_segment(surface, p_mid, points[10], points[11]),
+        _project_line_segment(surface, p_mid, points[11], points[6]),
+        _project_line_segment(surface, p_apical, points[12], points[13]),
+        _project_line_segment(surface, p_apical, points[13], points[14]),
+        _project_line_segment(surface, p_apical, points[14], points[15]),
+        _project_line_segment(surface, p_apical, points[15], points[16]),
+        _project_line_segment(surface, p_apical, points[16], points[17]),
+        _project_line_segment(surface, p_apical, points[17], points[12]),
+        _project_line_segment(surface, p_apical, points[18], points[19]),
+        _project_line_segment(surface, p_apical, points[19], points[20]),
+        _project_line_segment(surface, p_apical, points[20], points[21]),
+        _project_line_segment(surface, p_apical, points[21], points[18]),
+        _project_line_segment(surface, apex_endo, points[22], points[23]),
+        _project_line_segment(surface, apex_endo, points[23], points[24]),
+        _project_line_segment(surface, apex_endo, points[24], points[25]),
+        _project_line_segment(surface, apex_endo, points[25], points[22]),
     ]
 
     vline = [
-        _project_line_segment(surf, p_basal, coords[0], coords[6]),
-        _project_line_segment(surf, p_basal, coords[1], coords[7]),
-        _project_line_segment(surf, p_basal, coords[2], coords[8]),
-        _project_line_segment(surf, p_basal, coords[3], coords[9]),
-        _project_line_segment(surf, p_basal, coords[4], coords[10]),
-        _project_line_segment(surf, p_basal, coords[5], coords[11]),
-        _project_line_segment(surf, p_mid, coords[6], coords[12]),
-        _project_line_segment(surf, p_mid, coords[7], coords[13]),
-        _project_line_segment(surf, p_mid, coords[8], coords[14]),
-        _project_line_segment(surf, p_mid, coords[9], coords[15]),
-        _project_line_segment(surf, p_mid, coords[10], coords[16]),
-        _project_line_segment(surf, p_mid, coords[11], coords[17]),
-        _project_line_segment(surf, p_apical, coords[18], coords[22]),
-        _project_line_segment(surf, p_apical, coords[19], coords[23]),
-        _project_line_segment(surf, p_apical, coords[20], coords[24]),
-        _project_line_segment(surf, p_apical, coords[21], coords[25]),
+        _project_line_segment(surface, p_basal, points[0], points[6]),
+        _project_line_segment(surface, p_basal, points[1], points[7]),
+        _project_line_segment(surface, p_basal, points[2], points[8]),
+        _project_line_segment(surface, p_basal, points[3], points[9]),
+        _project_line_segment(surface, p_basal, points[4], points[10]),
+        _project_line_segment(surface, p_basal, points[5], points[11]),
+        _project_line_segment(surface, p_mid, points[6], points[12]),
+        _project_line_segment(surface, p_mid, points[7], points[13]),
+        _project_line_segment(surface, p_mid, points[8], points[14]),
+        _project_line_segment(surface, p_mid, points[9], points[15]),
+        _project_line_segment(surface, p_mid, points[10], points[16]),
+        _project_line_segment(surface, p_mid, points[11], points[17]),
+        _project_line_segment(surface, p_apical, points[18], points[22]),
+        _project_line_segment(surface, p_apical, points[19], points[23]),
+        _project_line_segment(surface, p_apical, points[20], points[24]),
+        _project_line_segment(surface, p_apical, points[21], points[25]),
     ]
 
-    return surf, coords, hline, vline
+    return points, hline, vline
 
 
 def _project_line_segment(
@@ -414,12 +425,12 @@ def compute_aha17(
     return aha_ids
 
 
-def _calculate_rotation_axis(short_axis, l4cv_axis):
-    short_normal = short_axis["normal"]
+def _calculate_rotation_axis(short: dict, long: dict):
+    short_normal = short["normal"]
 
     # define reference cut plane
     # default: rotate 60 from long axis
-    long_axis = l4cv_axis["normal"]
+    long_axis = long["normal"]
     axe_60 = Rotation.from_rotvec(np.radians(60) * short_normal).apply(  # noqa:E501
         long_axis
     )
diff --git a/tests/heart/utils/test_landmarks.py b/tests/heart/utils/test_landmarks.py
index 4cac50733..714e5ffc9 100644
--- a/tests/heart/utils/test_landmarks.py
+++ b/tests/heart/utils/test_landmarks.py
@@ -24,7 +24,6 @@
 
 import numpy as np
 import pytest
-import pyvista as pv
 
 import ansys.health.heart.models as models
 from ansys.health.heart.utils.landmark_utils import (
@@ -83,9 +82,9 @@ def test_compute_aha17(model):
     }
 
     aha_ids = compute_aha17(model, short, l4cv)
-    # solid = model.mesh.extract_cells_by_type(10)
-    # solid.cell_data["aha"] = aha_ids
-    # solid.save("aha.vtk")
+    solid = model.mesh.extract_cells_by_type(10)
+    solid.cell_data["aha"] = aha_ids
+    solid.save("aha.vtk")
 
     assert np.sum(aha_ids == 10) == 10591
     assert np.sum(aha_ids == 12) == 11220
@@ -94,11 +93,8 @@ def test_compute_aha17(model):
 
 
 def test_compute_aha17_landmarks(model):
-    l4cv = {
-        "center": np.array([30.05990578, 109.49603257, 375.8178529]),
-        "normal": np.array([-0.54847906, -0.10082087, -0.83006378]),
-    }
-    l4cv = {
+    # note give l2cv
+    l2cv = {
         "center": np.array([42.8065945, 105.236255, 367.91265619]),
         "normal": np.array([0.57496125, 0.67530147, -0.46193883]),
     }
@@ -106,15 +102,19 @@ def test_compute_aha17_landmarks(model):
         "center": np.array([26.07305844, 125.37379059, 376.52369382]),
         "normal": np.array([0.60711636, -0.73061828, -0.31242063]),
     }
-    surf, coords, hline, vline = compute_aha17_landmarks(model, short, l4cv)
-
-    # surf.save("surf.vtk")
-
-    # np.savetxt("points.txt", coords)
-    hlines = pv.MultiBlock(hline)
-    hlines.save("hlines.vtm")
-    vlines = pv.MultiBlock(vline)
-    vlines.save("vlines.vtm")
+    coords, hline, vline = compute_aha17_landmarks(model, short, l2cv)
+
+    # np.savetxt("points.txt", np.array(coords))
+    # hlines = pv.MultiBlock(hline)
+    # hlines.save("hlines.vtm")
+    # vlines = pv.MultiBlock(vline)
+    # vlines.save("vlines.vtm")
+
+    assert np.allclose(coords[0], np.array([5.5856953, 113.95523, 363.41443]))
+    assert pytest.approx(hline[0].length, abs=0.1) == 31.6
+    assert pytest.approx(hline[10].length, abs=0.1) == 29.8
+    assert pytest.approx(vline[0].length, abs=0.1) == 26.9
+    assert pytest.approx(vline[10].length, abs=0.1) == 26.4
 
 
 def test_compute_element_cs(model):

From add0fe9761c1079dc8168b13153461630e53927c Mon Sep 17 00:00:00 2001
From: pyansys-ci-bot <92810346+pyansys-ci-bot@users.noreply.github.com>
Date: Wed, 13 Aug 2025 13:05:05 +0000
Subject: [PATCH 04/14] chore: adding changelog file 1212.added.md
 [dependabot-skip]

---
 doc/source/changelog/1212.added.md | 1 +
 1 file changed, 1 insertion(+)
 create mode 100644 doc/source/changelog/1212.added.md

diff --git a/doc/source/changelog/1212.added.md b/doc/source/changelog/1212.added.md
new file mode 100644
index 000000000..e96cc9dfd
--- /dev/null
+++ b/doc/source/changelog/1212.added.md
@@ -0,0 +1 @@
+New aha strain

From 89fbfcb7820fe49b37b1c3b812c51e048e27849e Mon Sep 17 00:00:00 2001
From: wenfengye <wenfeng.ye@ansys.com>
Date: Wed, 13 Aug 2025 16:27:16 +0200
Subject: [PATCH 05/14] add aha strain

---
 .../health/heart/post/strain_calculator.py    | 108 +++++++++++++++++-
 .../health/heart/utils/landmark_utils.py      |   3 +-
 tests/heart/post/test_postprocess.py          |  11 ++
 3 files changed, 120 insertions(+), 2 deletions(-)

diff --git a/src/ansys/health/heart/post/strain_calculator.py b/src/ansys/health/heart/post/strain_calculator.py
index f6c700ce3..a771d12c7 100644
--- a/src/ansys/health/heart/post/strain_calculator.py
+++ b/src/ansys/health/heart/post/strain_calculator.py
@@ -29,7 +29,11 @@
 
 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_landmarks,
+    compute_element_cs,
+)
 from ansys.health.heart.utils.vtk_utils import find_corresponding_points, generate_thickness_lines
 
 
@@ -54,6 +58,108 @@ def __init__(self, model: HeartModel, d3plot_file):
 
         self.d3plot = D3plotReader(d3plot_file)
 
+    def _compute_new_strain(self, time_array: np.ndarray | list = None):
+        save = True
+        if time_array is None:
+            time_array = self.d3plot.time
+
+        surface_endo = self.model.left_ventricle.endocardium.copy()
+        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"]]
+            _, hlines, vlines = compute_aha17_landmarks(
+                self.model, self.model.short_axis, self.model.l2cv_axis, surface=surface_endo
+            )
+
+            hlength_array.append([hl.length for hl in hlines])
+            vlength_array.append([vl.length for vl in vlines])
+
+            if save:
+                surface_endo.save(f"endo_{it}.vtp")
+                pv.MultiBlock(hlines).save(f"hlines_{it}.vtm")
+                pv.MultiBlock(vlines).save(f"vlines_{it}.vtm")
+
+        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.
 
diff --git a/src/ansys/health/heart/utils/landmark_utils.py b/src/ansys/health/heart/utils/landmark_utils.py
index 400ea8dc9..fb7102e08 100644
--- a/src/ansys/health/heart/utils/landmark_utils.py
+++ b/src/ansys/health/heart/utils/landmark_utils.py
@@ -263,7 +263,8 @@ def _calculate_longitudinal_points(model, short_axis):
     # to have a flat segment 17, project endocardical apex point on short axis
     x = apex_ed - apex_ep
     y = p_basal - apex_ep
-    apex_ed = y * np.dot(x, y) / np.dot(y, y) + apex_ep
+    # TODO: temporary solution, move the apex upper so it can be projected to surface
+    apex_ed = 1.5 * y * np.dot(x, y) / np.dot(y, y) + apex_ep
     return p_basal, p_mid, p_apical, apex_ed, apex_ep
 
 
diff --git a/tests/heart/post/test_postprocess.py b/tests/heart/post/test_postprocess.py
index a63070252..be07ccad6 100644
--- a/tests/heart/post/test_postprocess.py
+++ b/tests/heart/post/test_postprocess.py
@@ -62,6 +62,17 @@ def test_compute_thickness(get_left_ventricle):
     assert lines[0]["thickness"][0] == pytest.approx(6.75, abs=0.1)
 
 
+@pytest.mark.requires_dpf
+def test_compute_strain(get_left_ventricle):
+    test_dir = get_left_ventricle[0]
+    model = get_left_ventricle[1]
+    d3plot = os.path.join(os.path.join(test_dir, "main", "d3plot"))
+
+    s = AhaStrainCalculator(model, d3plot)
+    s._compute_new_strain()
+    assert True
+
+
 @pytest.mark.requires_dpf
 def test_compute_myocardial_strain(get_left_ventricle):
     test_dir = get_left_ventricle[0]

From 0186cab97141a7b388cb5da9ab613243e9b27c20 Mon Sep 17 00:00:00 2001
From: wenfengye <wenfeng.ye@ansys.com>
Date: Thu, 14 Aug 2025 10:34:38 +0200
Subject: [PATCH 06/14] raise error when ray_tracing fails

---
 .../health/heart/post/strain_calculator.py    | 11 +++++---
 .../health/heart/utils/landmark_utils.py      | 25 +++++++++++++------
 2 files changed, 25 insertions(+), 11 deletions(-)

diff --git a/src/ansys/health/heart/post/strain_calculator.py b/src/ansys/health/heart/post/strain_calculator.py
index a771d12c7..17f1db960 100644
--- a/src/ansys/health/heart/post/strain_calculator.py
+++ b/src/ansys/health/heart/post/strain_calculator.py
@@ -71,10 +71,13 @@ def _compute_new_strain(self, time_array: np.ndarray | list = None):
                 self.d3plot.model.results.coordinates.on_time_scoping(float(t)).eval()[0].data
             )
             surface_endo.points = coordinates[surface_endo["_global-point-ids"]]
-            _, hlines, vlines = compute_aha17_landmarks(
-                self.model, self.model.short_axis, self.model.l2cv_axis, surface=surface_endo
-            )
-
+            try:
+                _, hlines, vlines = compute_aha17_landmarks(
+                    self.model, self.model.short_axis, self.model.l2cv_axis, surface=surface_endo
+                )
+            except ValueError as e:
+                print(f"Error 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])
 
diff --git a/src/ansys/health/heart/utils/landmark_utils.py b/src/ansys/health/heart/utils/landmark_utils.py
index fb7102e08..53c27d0cf 100644
--- a/src/ansys/health/heart/utils/landmark_utils.py
+++ b/src/ansys/health/heart/utils/landmark_utils.py
@@ -149,12 +149,18 @@ def compute_aha17_landmarks(
 
     for center in [p_basal, p_mid, p_apical]:
         for normal in [-axe_60, -axe_120, axe_180, axe_60, axe_120, -axe_180]:
-            point, _ = surface.ray_trace(center, center + 1e6 * normal, first_point=True)
-            points.append(point)
+            point = surface.ray_trace(center, center + 1e6 * normal, first_point=True)[0]
+            if point.size == 0:
+                raise ValueError("Cannot find point on surface from basal, mid, or apical.")
+            else:
+                points.append(point)
     for center in [p_apical, apex_endo]:
         for normal in [-axe_45, -axe_135, axe_45, axe_135]:
-            point, _ = surface.ray_trace(center, center + 1e6 * normal, first_point=True)
-            points.append(point)
+            point = surface.ray_trace(center, center + 1e6 * normal, first_point=True)[0]
+            if point.size == 0:
+                raise ValueError("Cannot find point on surface from apical or apex.")
+            else:
+                points.append(point)
 
     hline = [
         _project_line_segment(surface, p_basal, points[0], points[1]),
@@ -234,9 +240,14 @@ def _project_line_segment(
     projected_points = []
     for pt in segment_points:
         start = center
-        end = center + 1e9 * (pt - center)
-        intersection = surf.ray_trace(start, end, first_point=True)[0]
-        if intersection.size > 0:
+        end = center + 5.0 * (pt - center)
+        intersection = surf.copy().ray_trace(start, end, first_point=True)[0]
+        if intersection.size == 0:
+            # pv.PolyData([p1, p2, center,pt]).save("debug_points.vtp")
+            # pv.Line(start,end).save("debug_points2.vtp")
+            # surf.save("debug_endo.vtp")
+            raise ValueError(f"Cannot find intersection for segment from {p1} to {p2} on surface.")
+        else:
             projected_points.append(intersection)
 
     # Convert to array

From 333923f6f90539e6f2773f9054fc08c9045cd7ad Mon Sep 17 00:00:00 2001
From: wenfengye <wenfeng.ye@ansys.com>
Date: Thu, 14 Aug 2025 12:01:11 +0200
Subject: [PATCH 07/14] minor

---
 src/ansys/health/heart/utils/landmark_utils.py | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/src/ansys/health/heart/utils/landmark_utils.py b/src/ansys/health/heart/utils/landmark_utils.py
index 53c27d0cf..65862e2a4 100644
--- a/src/ansys/health/heart/utils/landmark_utils.py
+++ b/src/ansys/health/heart/utils/landmark_utils.py
@@ -149,15 +149,17 @@ def compute_aha17_landmarks(
 
     for center in [p_basal, p_mid, p_apical]:
         for normal in [-axe_60, -axe_120, axe_180, axe_60, axe_120, -axe_180]:
-            point = surface.ray_trace(center, center + 1e6 * normal, first_point=True)[0]
+            point = surface.copy().ray_trace(center, center + 1e3 * normal, first_point=True)[0]
             if point.size == 0:
                 raise ValueError("Cannot find point on surface from basal, mid, or apical.")
             else:
                 points.append(point)
     for center in [p_apical, apex_endo]:
         for normal in [-axe_45, -axe_135, axe_45, axe_135]:
-            point = surface.ray_trace(center, center + 1e6 * normal, first_point=True)[0]
+            point = surface.copy().ray_trace(center, center + 1e3 * normal, first_point=True)[0]
             if point.size == 0:
+                surface.save("debug_endo.vtp")
+                pv.Line(center, center + 1e2 * normal).save("debug_points.vtp")
                 raise ValueError("Cannot find point on surface from apical or apex.")
             else:
                 points.append(point)
@@ -234,13 +236,14 @@ def _project_line_segment(
     pv.PolyData
         The projected line segment.
     """
-    segment_points = np.linspace(p1, p2, num=100)
+    segment_points = np.linspace(p1, p2, num=10)
 
     # Project each point
     projected_points = []
     for pt in segment_points:
         start = center
-        end = center + 5.0 * (pt - center)
+        end = center + 100.0 * (pt - center)
+        # Note: copy is needed, or ray_trace will fail after several times of calling
         intersection = surf.copy().ray_trace(start, end, first_point=True)[0]
         if intersection.size == 0:
             # pv.PolyData([p1, p2, center,pt]).save("debug_points.vtp")

From e38143642aadb9997f2cdca309b0f2bcd7c8f149 Mon Sep 17 00:00:00 2001
From: wenfengye <wenfeng.ye@ansys.com>
Date: Thu, 14 Aug 2025 16:16:01 +0200
Subject: [PATCH 08/14] split method

---
 .../health/heart/post/strain_calculator.py    |  6 ++--
 .../health/heart/utils/landmark_utils.py      | 30 +++++++++++++++++--
 tests/heart/utils/test_landmarks.py           |  9 ++++--
 3 files changed, 38 insertions(+), 7 deletions(-)

diff --git a/src/ansys/health/heart/post/strain_calculator.py b/src/ansys/health/heart/post/strain_calculator.py
index 17f1db960..d8817a43e 100644
--- a/src/ansys/health/heart/post/strain_calculator.py
+++ b/src/ansys/health/heart/post/strain_calculator.py
@@ -31,7 +31,8 @@
 from ansys.health.heart.post.dpf_utils import D3plotReader
 from ansys.health.heart.utils.landmark_utils import (
     compute_aha17,
-    compute_aha17_landmarks,
+    compute_aha17_lines,
+    compute_aha17_points,
     compute_element_cs,
 )
 from ansys.health.heart.utils.vtk_utils import find_corresponding_points, generate_thickness_lines
@@ -72,9 +73,10 @@ def _compute_new_strain(self, time_array: np.ndarray | list = None):
             )
             surface_endo.points = coordinates[surface_endo["_global-point-ids"]]
             try:
-                _, hlines, vlines = compute_aha17_landmarks(
+                points = compute_aha17_points(
                     self.model, self.model.short_axis, self.model.l2cv_axis, surface=surface_endo
                 )
+                hlines, vlines = compute_aha17_lines(surface_endo, points)
             except ValueError as e:
                 print(f"Error computing AHA strain at time {t}: {e}")
                 exit()
diff --git a/src/ansys/health/heart/utils/landmark_utils.py b/src/ansys/health/heart/utils/landmark_utils.py
index 65862e2a4..8ab1e6d6e 100644
--- a/src/ansys/health/heart/utils/landmark_utils.py
+++ b/src/ansys/health/heart/utils/landmark_utils.py
@@ -80,7 +80,7 @@ def compute_anatomy_axis(
     return (l4cv_axis, l2cv_axis, short_axis)
 
 
-def compute_aha17_landmarks(
+def compute_aha17_points(
     model: HeartModel, short_axis: dict, long_axis: dict, surface: pv.PolyData = None
 ) -> tuple[list[np.ndarray], list[pv.PolyData], list[pv.PolyData]]:
     """Compute AHA17 landmarks for the left ventricle.
@@ -164,6 +164,32 @@ def compute_aha17_landmarks(
             else:
                 points.append(point)
 
+    return points
+
+
+def compute_aha17_lines(
+    surface: pv.PolyData, points: list[np.ndarray]
+) -> tuple[list[pv.PolyData], list[pv.PolyData]]:
+    """Compute AHA 17 lines.
+
+    Parameters
+    ----------
+    surface : pv.PolyData
+        The surface to project the l
+        lines onto.
+    points : list[np.ndarray]
+        The list of points defining the lines.
+
+    Returns
+    -------
+    tuple[list[pv.PolyData], list[pv.PolyData]]
+        The horizontal and vertical lines.
+    """
+    p_basal = np.mean(np.array(points[0:6]), axis=0)
+    p_mid = np.mean(np.array(points[6:12]), axis=0)
+    p_apical = np.mean(np.array(points[12:18]), axis=0)
+    apex_endo = np.mean(np.array(points[18:24]), axis=0)
+
     hline = [
         _project_line_segment(surface, p_basal, points[0], points[1]),
         _project_line_segment(surface, p_basal, points[1], points[2]),
@@ -212,7 +238,7 @@ def compute_aha17_landmarks(
         _project_line_segment(surface, p_apical, points[21], points[25]),
     ]
 
-    return points, hline, vline
+    return hline, vline
 
 
 def _project_line_segment(
diff --git a/tests/heart/utils/test_landmarks.py b/tests/heart/utils/test_landmarks.py
index 714e5ffc9..e88c36148 100644
--- a/tests/heart/utils/test_landmarks.py
+++ b/tests/heart/utils/test_landmarks.py
@@ -28,7 +28,8 @@
 import ansys.health.heart.models as models
 from ansys.health.heart.utils.landmark_utils import (
     compute_aha17,
-    compute_aha17_landmarks,
+    compute_aha17_lines,
+    compute_aha17_points,
     compute_anatomy_axis,
     compute_element_cs,
 )
@@ -102,7 +103,8 @@ def test_compute_aha17_landmarks(model):
         "center": np.array([26.07305844, 125.37379059, 376.52369382]),
         "normal": np.array([0.60711636, -0.73061828, -0.31242063]),
     }
-    coords, hline, vline = compute_aha17_landmarks(model, short, l2cv)
+    coords = compute_aha17_points(model, short, l2cv)
+    assert np.allclose(coords[0], np.array([5.5856953, 113.95523, 363.41443]))
 
     # np.savetxt("points.txt", np.array(coords))
     # hlines = pv.MultiBlock(hline)
@@ -110,7 +112,8 @@ def test_compute_aha17_landmarks(model):
     # vlines = pv.MultiBlock(vline)
     # vlines.save("vlines.vtm")
 
-    assert np.allclose(coords[0], np.array([5.5856953, 113.95523, 363.41443]))
+    hline, vline = compute_aha17_lines(model.left_ventricle.endocardium, coords)
+
     assert pytest.approx(hline[0].length, abs=0.1) == 31.6
     assert pytest.approx(hline[10].length, abs=0.1) == 29.8
     assert pytest.approx(vline[0].length, abs=0.1) == 26.9

From b0f0cef8805f03569d4040a403a0f79f2a7d0c85 Mon Sep 17 00:00:00 2001
From: wenfengye <wenfeng.ye@ansys.com>
Date: Thu, 14 Aug 2025 16:27:30 +0200
Subject: [PATCH 09/14] trace node for AHA strain

---
 .../health/heart/post/strain_calculator.py     | 18 ++++++++++++++----
 1 file changed, 14 insertions(+), 4 deletions(-)

diff --git a/src/ansys/health/heart/post/strain_calculator.py b/src/ansys/health/heart/post/strain_calculator.py
index d8817a43e..a21bf93ca 100644
--- a/src/ansys/health/heart/post/strain_calculator.py
+++ b/src/ansys/health/heart/post/strain_calculator.py
@@ -64,7 +64,16 @@ def _compute_new_strain(self, time_array: np.ndarray | list = None):
         if time_array is None:
             time_array = self.d3plot.time
 
-        surface_endo = self.model.left_ventricle.endocardium.copy()
+        surface_endo: pv.PolyData = self.model.left_ventricle.endocardium.copy()
+
+        # AHA points are bounded to nodes of endocardium surface.
+        points = compute_aha17_points(
+            self.model, self.model.short_axis, self.model.l2cv_axis, surface=surface_endo
+        )
+        ids = []
+        for p in points:
+            ids.append(surface_endo.find_closest_point(p))
+
         hlength_array = []
         vlength_array = []
         for it, t in enumerate(time_array):
@@ -73,9 +82,10 @@ def _compute_new_strain(self, time_array: np.ndarray | list = None):
             )
             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 = 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:
                 print(f"Error computing AHA strain at time {t}: {e}")

From 7fd37760eb5a4ea3a0a2dbe6d1413fd44e312aea Mon Sep 17 00:00:00 2001
From: wenfengye <wenfeng.ye@ansys.com>
Date: Thu, 14 Aug 2025 17:20:48 +0200
Subject: [PATCH 10/14] minor fix

---
 .../health/heart/post/strain_calculator.py     |  6 +++---
 src/ansys/health/heart/utils/landmark_utils.py | 18 ++++++++++--------
 2 files changed, 13 insertions(+), 11 deletions(-)

diff --git a/src/ansys/health/heart/post/strain_calculator.py b/src/ansys/health/heart/post/strain_calculator.py
index a21bf93ca..7445ee784 100644
--- a/src/ansys/health/heart/post/strain_calculator.py
+++ b/src/ansys/health/heart/post/strain_calculator.py
@@ -70,9 +70,7 @@ def _compute_new_strain(self, time_array: np.ndarray | list = None):
         points = compute_aha17_points(
             self.model, self.model.short_axis, self.model.l2cv_axis, surface=surface_endo
         )
-        ids = []
-        for p in points:
-            ids.append(surface_endo.find_closest_point(p))
+        ids = [surface_endo.find_closest_point(p) for p in points]
 
         hlength_array = []
         vlength_array = []
@@ -87,9 +85,11 @@ def _compute_new_strain(self, time_array: np.ndarray | list = None):
                 # )
                 points = [surface_endo.points[id] for id in ids]
                 hlines, vlines = compute_aha17_lines(surface_endo, points)
+
             except ValueError as e:
                 print(f"Error 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])
 
diff --git a/src/ansys/health/heart/utils/landmark_utils.py b/src/ansys/health/heart/utils/landmark_utils.py
index 8ab1e6d6e..288d6baa7 100644
--- a/src/ansys/health/heart/utils/landmark_utils.py
+++ b/src/ansys/health/heart/utils/landmark_utils.py
@@ -82,7 +82,7 @@ def compute_anatomy_axis(
 
 def compute_aha17_points(
     model: HeartModel, short_axis: dict, long_axis: dict, surface: pv.PolyData = None
-) -> tuple[list[np.ndarray], list[pv.PolyData], list[pv.PolyData]]:
+) -> list[np.ndarray]:
     """Compute AHA17 landmarks for the left ventricle.
 
     The AHA17 landmarks are defined by
@@ -154,12 +154,13 @@ def compute_aha17_points(
                 raise ValueError("Cannot find point on surface from basal, mid, or apical.")
             else:
                 points.append(point)
+
     for center in [p_apical, apex_endo]:
         for normal in [-axe_45, -axe_135, axe_45, axe_135]:
             point = surface.copy().ray_trace(center, center + 1e3 * normal, first_point=True)[0]
             if point.size == 0:
-                surface.save("debug_endo.vtp")
-                pv.Line(center, center + 1e2 * normal).save("debug_points.vtp")
+                # surface.save("debug_endo.vtp")
+                # pv.Line(center, center + 1e2 * normal).save("debug_points.vtp")
                 raise ValueError("Cannot find point on surface from apical or apex.")
             else:
                 points.append(point)
@@ -187,8 +188,8 @@ def compute_aha17_lines(
     """
     p_basal = np.mean(np.array(points[0:6]), axis=0)
     p_mid = np.mean(np.array(points[6:12]), axis=0)
-    p_apical = np.mean(np.array(points[12:18]), axis=0)
-    apex_endo = np.mean(np.array(points[18:24]), axis=0)
+    p_apical = np.mean(np.array(points[12:22]), axis=0)
+    apex_endo = np.mean(np.array(points[22:25]), axis=0)
 
     hline = [
         _project_line_segment(surface, p_basal, points[0], points[1]),
@@ -266,20 +267,21 @@ def _project_line_segment(
 
     # Project each point
     projected_points = []
+
     for pt in segment_points:
         start = center
         end = center + 100.0 * (pt - center)
-        # Note: copy is needed, or ray_trace will fail after several times of calling
+        # NOTE: copy() is needed, or ray_trace will fail after several times of calling
         intersection = surf.copy().ray_trace(start, end, first_point=True)[0]
         if intersection.size == 0:
             # pv.PolyData([p1, p2, center,pt]).save("debug_points.vtp")
             # pv.Line(start,end).save("debug_points2.vtp")
             # surf.save("debug_endo.vtp")
-            raise ValueError(f"Cannot find intersection for segment from {p1} to {p2} on surface.")
+            # raise ValueError(f"Cannot find intersection for segment.")
+            print("Cannot find intersection on the surface.")
         else:
             projected_points.append(intersection)
 
-    # Convert to array
     projected_points = np.array(projected_points)
 
     # Create curve

From 76710b2ddd6e3f4395a8d30794258609f5059794 Mon Sep 17 00:00:00 2001
From: wenfengye <wenfeng.ye@ansys.com>
Date: Thu, 21 Aug 2025 15:12:42 +0200
Subject: [PATCH 11/14] clean up

---
 .../health/heart/post/strain_calculator.py    | 41 ++++++++++++++++---
 .../health/heart/utils/landmark_utils.py      |  3 +-
 tests/heart/post/test_postprocess.py          |  8 ++--
 3 files changed, 43 insertions(+), 9 deletions(-)

diff --git a/src/ansys/health/heart/post/strain_calculator.py b/src/ansys/health/heart/post/strain_calculator.py
index 7445ee784..4ffcd3bb7 100644
--- a/src/ansys/health/heart/post/strain_calculator.py
+++ b/src/ansys/health/heart/post/strain_calculator.py
@@ -24,9 +24,11 @@
 
 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 (
@@ -59,14 +61,30 @@ def __init__(self, model: HeartModel, d3plot_file):
 
         self.d3plot = D3plotReader(d3plot_file)
 
-    def _compute_new_strain(self, time_array: np.ndarray | list = None):
-        save = True
+    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()
 
-        # AHA points are bounded to nodes of endocardium surface.
+        # Find AHA points on the endocardium surface.
         points = compute_aha17_points(
             self.model, self.model.short_axis, self.model.l2cv_axis, surface=surface_endo
         )
@@ -87,13 +105,13 @@ def _compute_new_strain(self, time_array: np.ndarray | list = None):
                 hlines, vlines = compute_aha17_lines(surface_endo, points)
 
             except ValueError as e:
-                print(f"Error computing AHA strain at time {t}: {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 save:
+            if vtk_dir is not None:
                 surface_endo.save(f"endo_{it}.vtp")
                 pv.MultiBlock(hlines).save(f"hlines_{it}.vtm")
                 pv.MultiBlock(vlines).save(f"vlines_{it}.vtm")
@@ -105,6 +123,7 @@ def _compute_new_strain(self, time_array: np.ndarray | list = None):
         # 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:
 
@@ -265,6 +284,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:
@@ -316,6 +339,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.
@@ -353,6 +380,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]]:
diff --git a/src/ansys/health/heart/utils/landmark_utils.py b/src/ansys/health/heart/utils/landmark_utils.py
index 288d6baa7..b86a4e540 100644
--- a/src/ansys/health/heart/utils/landmark_utils.py
+++ b/src/ansys/health/heart/utils/landmark_utils.py
@@ -29,6 +29,7 @@
 import pyvista as pv
 from scipy.spatial.transform import Rotation
 
+from ansys.health.heart import LOG as LOGGER
 from ansys.health.heart.models import HeartModel
 from ansys.health.heart.objects import CapType
 
@@ -278,7 +279,7 @@ def _project_line_segment(
             # pv.Line(start,end).save("debug_points2.vtp")
             # surf.save("debug_endo.vtp")
             # raise ValueError(f"Cannot find intersection for segment.")
-            print("Cannot find intersection on the surface.")
+            LOGGER.warning("Cannot find intersection on the surface.")
         else:
             projected_points.append(intersection)
 
diff --git a/tests/heart/post/test_postprocess.py b/tests/heart/post/test_postprocess.py
index be07ccad6..79594943e 100644
--- a/tests/heart/post/test_postprocess.py
+++ b/tests/heart/post/test_postprocess.py
@@ -68,9 +68,11 @@ def test_compute_strain(get_left_ventricle):
     model = get_left_ventricle[1]
     d3plot = os.path.join(os.path.join(test_dir, "main", "d3plot"))
 
-    s = AhaStrainCalculator(model, d3plot)
-    s._compute_new_strain()
-    assert True
+    c = AhaStrainCalculator(model, d3plot)
+    l_strain, c_strain = c.compute_longitudinal_radial_strain()
+    assert l_strain.shape == (2, 17)
+    assert l_strain[1, 1] == pytest.approx(0.000386, abs=1e-3)
+    assert c_strain[1, 1] == pytest.approx(-0.00183, abs=1e-3)
 
 
 @pytest.mark.requires_dpf

From b399f31d96a0de03679901d48a66420c80966052 Mon Sep 17 00:00:00 2001
From: wenfengye <wenfeng.ye@ansys.com>
Date: Fri, 22 Aug 2025 15:12:45 +0200
Subject: [PATCH 12/14] enhance save

---
 src/ansys/health/heart/post/strain_calculator.py | 7 ++++---
 src/ansys/health/heart/utils/landmark_utils.py   | 6 ++++++
 2 files changed, 10 insertions(+), 3 deletions(-)

diff --git a/src/ansys/health/heart/post/strain_calculator.py b/src/ansys/health/heart/post/strain_calculator.py
index 4ffcd3bb7..0ea13260a 100644
--- a/src/ansys/health/heart/post/strain_calculator.py
+++ b/src/ansys/health/heart/post/strain_calculator.py
@@ -112,9 +112,10 @@ def compute_longitudinal_radial_strain(
             vlength_array.append([vl.length for vl in vlines])
 
             if vtk_dir is not None:
-                surface_endo.save(f"endo_{it}.vtp")
-                pv.MultiBlock(hlines).save(f"hlines_{it}.vtm")
-                pv.MultiBlock(vlines).save(f"vlines_{it}.vtm")
+                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)
diff --git a/src/ansys/health/heart/utils/landmark_utils.py b/src/ansys/health/heart/utils/landmark_utils.py
index b86a4e540..ea6293e31 100644
--- a/src/ansys/health/heart/utils/landmark_utils.py
+++ b/src/ansys/health/heart/utils/landmark_utils.py
@@ -240,6 +240,12 @@ def compute_aha17_lines(
         _project_line_segment(surface, p_apical, points[21], points[25]),
     ]
 
+    for id, line in enumerate(hline):
+        line.cell_data["id"] = id + 1
+
+    for id, line in enumerate(vline):
+        line.cell_data["id"] = id + 1
+
     return hline, vline
 
 

From eacbce6d4608b9eb11068f241a098773616346d8 Mon Sep 17 00:00:00 2001
From: wenfengye <wenfeng.ye@ansys.com>
Date: Fri, 22 Aug 2025 15:45:09 +0200
Subject: [PATCH 13/14] use actual coordinates of apex, as surface is deformed
 after zerop

---
 src/ansys/health/heart/utils/landmark_utils.py | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/src/ansys/health/heart/utils/landmark_utils.py b/src/ansys/health/heart/utils/landmark_utils.py
index ea6293e31..48f70f2a4 100644
--- a/src/ansys/health/heart/utils/landmark_utils.py
+++ b/src/ansys/health/heart/utils/landmark_utils.py
@@ -301,9 +301,11 @@ def _calculate_longitudinal_points(model, short_axis):
     """Define landmarks along the short axis."""
     for apex in model.left_ventricle.apex_points:
         if "endocardium" in apex.name:
-            apex_ed = apex.xyz
+            surface = model.left_ventricle.endocardium
+            apex_ed = surface.points[np.where(surface["_global-point-ids"] == apex.node_id)[0][0]]
         elif "epicardium" in apex.name:
-            apex_ep = apex.xyz
+            surface = model.left_ventricle.epicardium
+            apex_ep = surface.points[np.where(surface["_global-point-ids"] == apex.node_id)[0][0]]
 
     p_basal = short_axis["center"]
     p_mid = 1 / 3 * (apex_ep - p_basal) + p_basal
@@ -312,7 +314,7 @@ def _calculate_longitudinal_points(model, short_axis):
     # to have a flat segment 17, project endocardical apex point on short axis
     x = apex_ed - apex_ep
     y = p_basal - apex_ep
-    # TODO: temporary solution, move the apex upper so it can be projected to surface
+    # slightly move the apex upper so it can be projected to surface
     apex_ed = 1.5 * y * np.dot(x, y) / np.dot(y, y) + apex_ep
     return p_basal, p_mid, p_apical, apex_ed, apex_ep
 

From 513017f404a582b28ffc0cbfd29cdcdc9a084aa7 Mon Sep 17 00:00:00 2001
From: wenfengye <wenfeng.ye@ansys.com>
Date: Fri, 22 Aug 2025 16:28:52 +0200
Subject: [PATCH 14/14] set to default stran

---
 src/ansys/health/heart/post/auto_process.py | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/src/ansys/health/heart/post/auto_process.py b/src/ansys/health/heart/post/auto_process.py
index 8948f847f..bcc7f879b 100644
--- a/src/ansys/health/heart/post/auto_process.py
+++ b/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