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
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
diff --git a/src/ansys/health/heart/post/strain_calculator.py b/src/ansys/health/heart/post/strain_calculator.py
index f6c700ce3..0ea13260a 100644
--- a/src/ansys/health/heart/post/strain_calculator.py
+++ b/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]]:
diff --git a/src/ansys/health/heart/utils/landmark_utils.py b/src/ansys/health/heart/utils/landmark_utils.py
index 31ad2844c..48f70f2a4 100644
--- a/src/ansys/health/heart/utils/landmark_utils.py
+++ b/src/ansys/health/heart/utils/landmark_utils.py
@@ -26,8 +26,10 @@
 
 from deprecated import deprecated
 import numpy as np
+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
 
@@ -79,6 +81,244 @@ def compute_anatomy_axis(
     return (l4cv_axis, l2cv_axis, short_axis)
 
 
+def compute_aha17_points(
+    model: HeartModel, short_axis: dict, long_axis: dict, surface: pv.PolyData = None
+) -> list[np.ndarray]:
+    """Compute AHA17 landmarks for the left ventricle.
+
+    The AHA17 landmarks are defined by
+    - 26 points
+    - 26 horizontal lines
+    - 16 vertical lines
+
+    They split left ventricle into 17 segments as follow:
+
+    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
+
+    Parameters
+    ----------
+    model : HeartModel
+        Heart model.
+    short_axis : dict
+        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
+    -------
+    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
+    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(
+        model, short_axis
+    )
+    axe_60, axe_120, axe_180, axe_45, axe_135 = _calculate_rotation_axis(short_axis, long_axis)
+
+    # define points for the segments
+    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 = 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.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)
+
+    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: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]),
+        _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(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]),
+    ]
+
+    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
+
+
+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=10)
+
+    # 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
+        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.")
+            LOGGER.warning("Cannot find intersection on the surface.")
+        else:
+            projected_points.append(intersection)
+
+    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:
+            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:
+            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
+    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
+    # 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
+
+
 def compute_aha17(
     model: HeartModel,
     short_axis: dict,
@@ -237,6 +477,24 @@ def compute_aha17(
     return aha_ids
 
 
+def _calculate_rotation_axis(short: dict, long: dict):
+    short_normal = short["normal"]
+
+    # define reference cut plane
+    # default: rotate 60 from long axis
+    long_axis = long["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/post/test_postprocess.py b/tests/heart/post/test_postprocess.py
index a63070252..79594943e 100644
--- a/tests/heart/post/test_postprocess.py
+++ b/tests/heart/post/test_postprocess.py
@@ -62,6 +62,19 @@ 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"))
+
+    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
 def test_compute_myocardial_strain(get_left_ventricle):
     test_dir = get_left_ventricle[0]
diff --git a/tests/heart/utils/test_landmarks.py b/tests/heart/utils/test_landmarks.py
index d4f5fadc0..e88c36148 100644
--- a/tests/heart/utils/test_landmarks.py
+++ b/tests/heart/utils/test_landmarks.py
@@ -28,6 +28,8 @@
 import ansys.health.heart.models as models
 from ansys.health.heart.utils.landmark_utils import (
     compute_aha17,
+    compute_aha17_lines,
+    compute_aha17_points,
     compute_anatomy_axis,
     compute_element_cs,
 )
@@ -81,9 +83,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
@@ -91,6 +93,33 @@ def test_compute_aha17(model):
     assert np.sum(np.isnan(aha_ids)) == 361818
 
 
+def test_compute_aha17_landmarks(model):
+    # note give l2cv
+    l2cv = {
+        "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]),
+    }
+    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)
+    # hlines.save("hlines.vtm")
+    # vlines = pv.MultiBlock(vline)
+    # vlines.save("vlines.vtm")
+
+    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
+    assert pytest.approx(vline[10].length, abs=0.1) == 26.4
+
+
 def test_compute_element_cs(model):
     l4cv = {
         "center": np.array([30.05990578, 109.49603257, 375.8178529]),