From c11ecbb5186305a8b2de6910adeeddfacdac5c5d Mon Sep 17 00:00:00 2001
From: wenfengye <wenfeng.ye@ansys.com>
Date: Thu, 19 Dec 2024 14:40:39 +0100
Subject: [PATCH 01/21] add David's script

---
 .../heart/misc/compute_fiber_aha_angles.py    | 282 ++++++++++++++++++
 1 file changed, 282 insertions(+)
 create mode 100644 src/ansys/heart/misc/compute_fiber_aha_angles.py

diff --git a/src/ansys/heart/misc/compute_fiber_aha_angles.py b/src/ansys/heart/misc/compute_fiber_aha_angles.py
new file mode 100644
index 000000000..2602edfa6
--- /dev/null
+++ b/src/ansys/heart/misc/compute_fiber_aha_angles.py
@@ -0,0 +1,282 @@
+# Copyright (C) 2023 - 2024 ANSYS, Inc. and/or its affiliates.
+# SPDX-License-Identifier: MIT
+#
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in all
+# copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.
+
+"""
+Compute average fiber orientations with respect to UHCs in each AHA region in the LV.
+
+from david.oks@medtronic.com
+"""
+
+import os
+import pathlib
+
+import matplotlib.pyplot as plt
+import numpy as np
+from numpy.linalg import norm
+import pandas as pd
+import pyvista as pv
+
+import ansys.heart.core.models as models
+from ansys.heart.simulator.simulator import DynaSettings, EPSimulator
+
+
+def angle_btw_vectors(x, y):
+    """Computes angles between N number of M-dimensional vectors
+    Input:  shape(x): (N,M)
+            shape(y): (N,M)
+    Output: shape(): (N,)
+    """  # noqa
+    return (
+        np.arccos(np.clip(np.sum(x * y, axis=1) / (norm(x, axis=1) * norm(y, axis=1)), -1, 1))
+        * 180
+        / np.pi
+    )
+
+
+def signed_angle_btw_vectors(x, y, n):
+    """Computes signed angles between N number of M-dimensional vectors
+        Input:  shape(x): (N,M)
+                shape(y): (N,M)
+                shape(n): (N,M)
+        Output: shape(): (N,)
+
+    Va . Vb == |Va| * |Vb| * cos(alpha)    (by definition)
+            == |Va| * |Vb| * cos(beta)     (cos(alpha) == cos(-alpha) == cos(360° - alpha)
+
+
+    Va x Vb == |Va| * |Vb| * sin(alpha) * n1
+        (by definition; n1 is a unit vector perpendicular to Va and Vb with
+         orientation matching the right-hand rule)
+
+    Therefore (again assuming Vn is normalized):
+       n1 . Vn == 1 when beta < 180
+       n1 . Vn == -1 when beta > 180
+
+    ==>  (Va x Vb) . Vn == |Va| * |Vb| * sin(beta)
+
+    """  # noqa
+    n /= norm(n, axis=1)[:, None]  # normalize normal vectors
+
+    return np.arctan2(np.sum(np.cross(x, y) * n, axis=1), np.sum(x * y, axis=1)) * 180 / np.pi
+
+
+def compute_aha_fiber_angles(mesh: pv.UnstructuredGrid, out_dir):
+    """Computes average fiber inclination in each AHA region
+    as a function of the transmural depth
+    """  # noqa
+    assert "aha17" in mesh.cell_data
+    assert "fiber" in mesh.cell_data
+    assert "sheet" in mesh.cell_data
+    assert "transmural" in mesh.point_data
+    assert "rotational" in mesh.point_data
+
+    aha_ids = mesh["aha17"]
+    aha_elements = np.where(~np.isnan(aha_ids))[0]
+    aha_model = mesh.extract_cells(aha_elements)
+    aha_model.cell_data["AHA"] = aha_ids[aha_elements]
+
+    # print("Nmbr of points:\t{:d}".format(mesh.n_points))
+    # print("Nmbr of cells:\t{:d}".format(mesh.n_cells))
+    # print("Nmbr of AHA cells:\t{:d}".format(len(aha_elements)))
+
+    # load fibers and sheets at cells
+    el_fibers = mesh.cell_data["fiber"]
+    el_sheets = mesh.cell_data["sheet"]
+    el_fibers = el_fibers[aha_elements]
+    el_sheets = el_sheets[aha_elements]
+    el_fibers /= np.linalg.norm(el_fibers, axis=1)[:, None]  # make sure fibers are normalized
+
+    # interpolate transmural depth from points to cells
+    el_depths = mesh.point_data_to_cell_data()["transmural"][aha_elements]
+    el_depths = 2.0 * el_depths - 1.0  # map from [0,1] -> [-1,1]
+    # interpolate rotational coordinates from points to cells
+    # el_rotat = mesh.point_data_to_cell_data()["rotational"][aha_elements]
+
+    # compute transmural vector from derivative of transmural depth
+    # TODO : interpolate t instead of grad_t
+    pt_grad_t = mesh.compute_derivative(scalars="transmural", preference="cell")["gradient"]
+    mesh.point_data.set_scalars(name="grad_t", scalars=pt_grad_t)
+    el_grad_t = mesh.point_data_to_cell_data()["grad_t"][aha_elements]
+    el_grad_t = el_grad_t / np.linalg.norm(el_grad_t, axis=1)[:, None]
+
+    # compute transmural vector from derivative of rotational coordinate
+    pt_grad_r = mesh.compute_derivative(scalars="rotational", preference="cell")["gradient"]
+    mesh.point_data.set_scalars(name="grad_r", scalars=pt_grad_r)
+    el_grad_r = mesh.point_data_to_cell_data()["grad_r"][aha_elements]
+    # set elements at the rotational coordinate discontinuity to NaN
+    # TODO: prescribe to average gradient from nearest neighbors
+    norm_grad_r = norm(el_grad_r, axis=1)
+    id_discont = np.where(norm_grad_r > np.average(norm_grad_r) + 2 * np.std(norm_grad_r))[0]
+    nans = np.empty((3,))
+    nans[:] = np.nan
+    el_grad_r[id_discont, :] = nans
+    el_grad_r = el_grad_r / np.linalg.norm(el_grad_r, axis=1)[:, None]
+
+    # compute angle between fibers and transmural vectors
+    el_angles_t = signed_angle_btw_vectors(el_fibers, el_grad_t, el_grad_r)
+
+    # compute angle between fibers and rotational vectors
+    el_angles_r = signed_angle_btw_vectors(el_grad_r, el_fibers, el_grad_t)
+
+    # get aha17 label for left ventricle elements
+    aha17_label = aha_ids[aha_elements]
+
+    # average fiber angles wrt to transmural and rotational coords in each AHA and depth
+    ndepths = 9
+    depth_bin_edges = np.linspace(-1.0, 1.0, ndepths + 1)
+    depth_bin_ctrs = 0.5 * (depth_bin_edges[:-1] + depth_bin_edges[1:])
+    el_angles_t_avg = np.zeros((len(aha_elements)))
+    aha_angles_t = np.zeros((17, ndepths))
+    el_angles_r_avg = np.zeros((len(aha_elements)))
+    aha_angles_r = np.zeros((17, ndepths))
+    for i in range(1, 18):
+        idx_aha = aha17_label == i
+        for j in range(ndepths):
+            depth_min = depth_bin_edges[j]
+            depth_max = depth_bin_edges[j + 1]
+            idx_depth = (el_depths >= depth_min) & (el_depths < depth_max)
+            idx = np.where(idx_aha & idx_depth)[0]
+            aha_angles_t[i - 1, j] = np.nanmean(el_angles_t[idx])
+            el_angles_t_avg[idx] = aha_angles_t[i - 1, j]
+            aha_angles_r[i - 1, j] = np.nanmean(el_angles_r[idx])
+            el_angles_r_avg[idx] = aha_angles_r[i - 1, j]
+
+    # save to vtk
+    aha_model.cell_data["transmural_angles"] = el_angles_t
+    aha_model.cell_data["transmural_angles_aha"] = el_angles_t_avg
+    aha_model.cell_data["rotational_angles"] = el_angles_r
+    aha_model.cell_data["rotational_angles_aha"] = el_angles_r_avg
+    aha_model.cell_data["transmural_depth"] = el_depths
+    aha_model.cell_data.set_vectors(el_fibers, "fibers", deep_copy=True)
+    aha_model.cell_data.set_vectors(el_sheets, "sheets", deep_copy=True)
+    aha_model.cell_data.set_vectors(el_grad_t, "grad_t", deep_copy=True)
+    aha_model.cell_data.set_vectors(el_grad_r, "grad_r", deep_copy=True)
+    aha_model.save(os.path.join(pathlib.Path(out_dir), "aha_averaged_angles.vtk"))
+
+    # save to csv
+    cols = ["{:1.2f}".format(x) for x in depth_bin_ctrs]
+    rows = ["{:d}".format(x) for x in range(1, 18)]
+    df_r = pd.DataFrame(data=aha_angles_r, index=rows, columns=cols)
+    df_r.to_csv(os.path.join(out_dir, "AHA_fiber_angles_r.csv"), index=True)
+    df_t = pd.DataFrame(data=aha_angles_t, index=rows, columns=cols)
+    df_t.to_csv(os.path.join(out_dir, "AHA_fiber_angles_t.csv"), index=True)
+
+    # print(df_r)
+    # print(df_t)
+
+    return df_r, df_t
+
+
+def plot_fiber_aha_angles(data: pd.DataFrame | str):
+    """
+    Plot average fiber helical orientation in each AHA as a function of transmural depth
+
+    """  # noqa
+    if isinstance(data, str):
+        df = pd.read_csv(os.path.join(data, "AHA_fiber_angles_t.csv"))
+    elif isinstance(data, pd.DataFrame):
+        df = data
+
+    aha_names = [
+        "Basal anterior",
+        "Basal anteroseptal",
+        "Basal inferoseptal",
+        "Basal inferior",
+        "Basal inferolateral",
+        "Basal anterolateral",
+        "Mid anterior",
+        "Mid anteroseptal",
+        "Mid inferoseptal",
+        "Mid inferior",
+        "Mid inferolateral",
+        "Mid anterolateral",
+        "Apical anterior",
+        "Apical inferior",
+        "Apical septal",
+        "Apical lateral",
+        "Apex",
+    ]
+
+    fig, axs = plt.subplots(3, 6)
+    fig.set_size_inches(31, 18)
+    # print(df.columns.tolist()[1:])
+    depths = [float(x) for x in df.columns.tolist()[1:]]
+    for iaha in range(17):
+        i = iaha // 6
+        j = iaha % 6
+        alphas = df.iloc[iaha].tolist()[1:]
+        print(alphas)
+        axs[i, j].plot(depths, alphas, "o-b")
+        axs[i, j].plot([-1, 1], [-60, -60], "--k")
+        axs[i, j].plot([-1, 1], [60, 60], "--k")
+        axs[i, j].set_title(aha_names[iaha])
+        axs[i, j].set_xlim(xmin=-1, xmax=1)
+        axs[i, j].set_ylim(ymin=-100, ymax=100)
+
+    for ax in axs.flat:
+        ax.set(xlabel="Transmural Depth", ylabel="$\\alpha_h$")
+
+    for ax in axs.flat:
+        ax.label_outer()
+
+    axs[2, 5].remove()
+
+    # plt.savefig("fiber_helical_angles.png", bbox_inches="tight")
+    plt.show()
+
+
+if __name__ == "__main__":
+    """Demo
+    """
+    # assumed LS-DYNA
+    lsdyna_path = r"D:\ansysdev\lsdyna\ls-dyna_smp_d_DEV_112901-gcbb8e36701_winx64_ifort190.exe"
+    dyna_settings = DynaSettings(
+        lsdyna_path=lsdyna_path, dynatype="smp", num_cpus=4, platform="windows"
+    )
+
+    # assumed case
+    workdir = os.path.abspath(os.path.join("downloads", "Rodero2021", "01", "BV"))
+    path_to_model = os.path.join(workdir, "heart_model.vtu")
+
+    model: models.BiVentricle = models.BiVentricle(workdir)
+    model.load_model_from_mesh(path_to_model, path_to_model.replace(".vtu", ".partinfo.json"))
+
+    # get fields data
+    simulator = EPSimulator(
+        model=model,
+        dyna_settings=dyna_settings,
+        simulation_directory=os.path.join(workdir, "simulation-EP"),
+    )
+    simulator.settings.load_defaults()
+    simulator.compute_uhc()
+    simulator.compute_fibers()
+
+    # get AHA labels
+    from ansys.heart.core.helpers.landmarks import compute_aha17
+
+    aha_ids = compute_aha17(simulator.model, simulator.model.short_axis, simulator.model.l4cv_axis)
+    grid = simulator.model.mesh.extract_cells_by_type(10)
+    grid.cell_data["aha17"] = aha_ids
+
+    # compute angles
+    df_r, df_t = compute_aha_fiber_angles(grid, workdir)
+    plot_fiber_aha_angles(df_t)

From eab5477a80ee85cf461c6fabe999a9e41ae6ca3c Mon Sep 17 00:00:00 2001
From: wenfengye <wenfeng.ye@ansys.com>
Date: Thu, 19 Dec 2024 15:06:32 +0100
Subject: [PATCH 02/21] minor

---
 src/ansys/heart/misc/compute_fiber_aha_angles.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/ansys/heart/misc/compute_fiber_aha_angles.py b/src/ansys/heart/misc/compute_fiber_aha_angles.py
index 2602edfa6..7ee1a7e37 100644
--- a/src/ansys/heart/misc/compute_fiber_aha_angles.py
+++ b/src/ansys/heart/misc/compute_fiber_aha_angles.py
@@ -23,7 +23,7 @@
 """
 Compute average fiber orientations with respect to UHCs in each AHA region in the LV.
 
-from david.oks@medtronic.com
+from davoks
 """
 
 import os

From ff57fc8f267ee7c0a8407517fb0de1c4ddb1c7fc Mon Sep 17 00:00:00 2001
From: wenfengye <wenfeng.ye@ansys.com>
Date: Thu, 19 Dec 2024 15:06:32 +0100
Subject: [PATCH 03/21] minor

---
 src/ansys/heart/misc/compute_fiber_aha_angles.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/ansys/heart/misc/compute_fiber_aha_angles.py b/src/ansys/heart/misc/compute_fiber_aha_angles.py
index 2602edfa6..e362c0758 100644
--- a/src/ansys/heart/misc/compute_fiber_aha_angles.py
+++ b/src/ansys/heart/misc/compute_fiber_aha_angles.py
@@ -23,7 +23,7 @@
 """
 Compute average fiber orientations with respect to UHCs in each AHA region in the LV.
 
-from david.oks@medtronic.com
+from davoks
 """
 
 import os
@@ -224,7 +224,7 @@ def plot_fiber_aha_angles(data: pd.DataFrame | str):
         i = iaha // 6
         j = iaha % 6
         alphas = df.iloc[iaha].tolist()[1:]
-        print(alphas)
+        # print(alphas)
         axs[i, j].plot(depths, alphas, "o-b")
         axs[i, j].plot([-1, 1], [-60, -60], "--k")
         axs[i, j].plot([-1, 1], [60, 60], "--k")

From 661a74b00ea37ed5aaa0e98cbc58b1cf07c2e92f Mon Sep 17 00:00:00 2001
From: mhoeijm <102799582+mhoeijm@users.noreply.github.com>
Date: Fri, 22 Aug 2025 10:32:07 +0200
Subject: [PATCH 04/21] fix: imports

---
 src/ansys/heart/misc/compute_fiber_aha_angles.py | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/ansys/heart/misc/compute_fiber_aha_angles.py b/src/ansys/heart/misc/compute_fiber_aha_angles.py
index e362c0758..4cdc10eba 100644
--- a/src/ansys/heart/misc/compute_fiber_aha_angles.py
+++ b/src/ansys/heart/misc/compute_fiber_aha_angles.py
@@ -1,4 +1,4 @@
-# Copyright (C) 2023 - 2024 ANSYS, Inc. and/or its affiliates.
+# Copyright (C) 2023 - 2025 ANSYS, Inc. and/or its affiliates.
 # SPDX-License-Identifier: MIT
 #
 #
@@ -35,11 +35,11 @@
 import pandas as pd
 import pyvista as pv
 
-import ansys.heart.core.models as models
-from ansys.heart.simulator.simulator import DynaSettings, EPSimulator
+import ansys.health.heart.models as models
+from ansys.health.heart.simulator import DynaSettings, EPSimulator
 
 
-def angle_btw_vectors(x, y):
+def angle_btw_vectors(x: np.ndarray, y: np.ndarray) -> np.ndarray:
     """Computes angles between N number of M-dimensional vectors
     Input:  shape(x): (N,M)
             shape(y): (N,M)
@@ -52,7 +52,7 @@ def angle_btw_vectors(x, y):
     )
 
 
-def signed_angle_btw_vectors(x, y, n):
+def signed_angle_between_vectors(x, y, n):
     """Computes signed angles between N number of M-dimensional vectors
         Input:  shape(x): (N,M)
                 shape(y): (N,M)
@@ -132,10 +132,10 @@ def compute_aha_fiber_angles(mesh: pv.UnstructuredGrid, out_dir):
     el_grad_r = el_grad_r / np.linalg.norm(el_grad_r, axis=1)[:, None]
 
     # compute angle between fibers and transmural vectors
-    el_angles_t = signed_angle_btw_vectors(el_fibers, el_grad_t, el_grad_r)
+    el_angles_t = signed_angle_between_vectors(el_fibers, el_grad_t, el_grad_r)
 
     # compute angle between fibers and rotational vectors
-    el_angles_r = signed_angle_btw_vectors(el_grad_r, el_fibers, el_grad_t)
+    el_angles_r = signed_angle_between_vectors(el_grad_r, el_fibers, el_grad_t)
 
     # get aha17 label for left ventricle elements
     aha17_label = aha_ids[aha_elements]
@@ -271,7 +271,7 @@ def plot_fiber_aha_angles(data: pd.DataFrame | str):
     simulator.compute_fibers()
 
     # get AHA labels
-    from ansys.heart.core.helpers.landmarks import compute_aha17
+    from ansys.health.heart.utils.landmark_utils import compute_aha17
 
     aha_ids = compute_aha17(simulator.model, simulator.model.short_axis, simulator.model.l4cv_axis)
     grid = simulator.model.mesh.extract_cells_by_type(10)

From d1084883d515a77e3b8d09cc1a57705b3e13ddd6 Mon Sep 17 00:00:00 2001
From: mhoeijm <102799582+mhoeijm@users.noreply.github.com>
Date: Fri, 22 Aug 2025 11:06:11 +0200
Subject: [PATCH 05/21] fix: remove asserts

---
 src/ansys/heart/misc/compute_fiber_aha_angles.py | 11 +++++------
 1 file changed, 5 insertions(+), 6 deletions(-)

diff --git a/src/ansys/heart/misc/compute_fiber_aha_angles.py b/src/ansys/heart/misc/compute_fiber_aha_angles.py
index 4cdc10eba..41967793b 100644
--- a/src/ansys/heart/misc/compute_fiber_aha_angles.py
+++ b/src/ansys/heart/misc/compute_fiber_aha_angles.py
@@ -39,7 +39,7 @@
 from ansys.health.heart.simulator import DynaSettings, EPSimulator
 
 
-def angle_btw_vectors(x: np.ndarray, y: np.ndarray) -> np.ndarray:
+def angle_between_vectors(x: np.ndarray, y: np.ndarray) -> np.ndarray:
     """Computes angles between N number of M-dimensional vectors
     Input:  shape(x): (N,M)
             shape(y): (N,M)
@@ -83,11 +83,10 @@ def compute_aha_fiber_angles(mesh: pv.UnstructuredGrid, out_dir):
     """Computes average fiber inclination in each AHA region
     as a function of the transmural depth
     """  # noqa
-    assert "aha17" in mesh.cell_data
-    assert "fiber" in mesh.cell_data
-    assert "sheet" in mesh.cell_data
-    assert "transmural" in mesh.point_data
-    assert "rotational" in mesh.point_data
+    expected_arrays = ["aha17", "fiber", "sheet", "transmural", "rotational"]
+    for arr in expected_arrays:
+        if arr not in mesh.cell_data and arr not in mesh.point_data:
+            raise ValueError(f"Expected array '{arr}' not found in mesh data.")
 
     aha_ids = mesh["aha17"]
     aha_elements = np.where(~np.isnan(aha_ids))[0]

From a0000b8cac2cfd332803225d7bea4383560daf12 Mon Sep 17 00:00:00 2001
From: mhoeijm <102799582+mhoeijm@users.noreply.github.com>
Date: Fri, 22 Aug 2025 13:28:46 +0200
Subject: [PATCH 06/21] refactor: move to utils

---
 .../heart/utils}/compute_fiber_aha_angles.py  | 37 ++++++++++++++++---
 1 file changed, 32 insertions(+), 5 deletions(-)
 rename src/ansys/{heart/misc => health/heart/utils}/compute_fiber_aha_angles.py (91%)

diff --git a/src/ansys/heart/misc/compute_fiber_aha_angles.py b/src/ansys/health/heart/utils/compute_fiber_aha_angles.py
similarity index 91%
rename from src/ansys/heart/misc/compute_fiber_aha_angles.py
rename to src/ansys/health/heart/utils/compute_fiber_aha_angles.py
index 41967793b..92b9fc727 100644
--- a/src/ansys/heart/misc/compute_fiber_aha_angles.py
+++ b/src/ansys/health/heart/utils/compute_fiber_aha_angles.py
@@ -52,8 +52,26 @@ def angle_between_vectors(x: np.ndarray, y: np.ndarray) -> np.ndarray:
     )
 
 
-def signed_angle_between_vectors(x, y, n):
-    """Computes signed angles between N number of M-dimensional vectors
+def signed_angle_between_vectors(x: np.ndarray, y: np.ndarray, n: np.ndarray):
+    """Compute signed angles between N number of M-dimensional vectors.
+
+    Parameters
+    ----------
+    x : np.ndarray
+        (N,M) array for which to compute the angles.
+    y : np.ndarray
+        (N,M) array with reference vectors.
+    n : np.ndarray
+        (N,M) array with normal vectors of the reference plane.
+
+    Returns
+    -------
+    np.ndarray
+        Array of signed angles in degrees of size (N,).
+
+    Notes
+    -----
+    Computes signed angles between N number of M-dimensional vectors
         Input:  shape(x): (N,M)
                 shape(y): (N,M)
                 shape(n): (N,M)
@@ -73,10 +91,19 @@ def signed_angle_between_vectors(x, y, n):
 
     ==>  (Va x Vb) . Vn == |Va| * |Vb| * sin(beta)
 
-    """  # noqa
-    n /= norm(n, axis=1)[:, None]  # normalize normal vectors
+    """
+    # normalize vectors
+    n /= norm(n, axis=1)[:, None]
+    x /= norm(x, axis=1)[:, None]
+    y /= norm(y, axis=1)[:, None]
+
+    nan_mask = np.isclose(np.linalg.norm(np.cross(x, n), axis=1), 0)
+
+    angles = np.arctan2(np.sum(np.cross(x, y) * n, axis=1), np.sum(x * y, axis=1))
+
+    angles[nan_mask] = np.nan  # Set angles where x is parallel to n to NaN
 
-    return np.arctan2(np.sum(np.cross(x, y) * n, axis=1), np.sum(x * y, axis=1)) * 180 / np.pi
+    return angles * 180 / np.pi
 
 
 def compute_aha_fiber_angles(mesh: pv.UnstructuredGrid, out_dir):

From b12e5234578f5d2b4330da5a4f6d933f404e4fdb Mon Sep 17 00:00:00 2001
From: mhoeijm <102799582+mhoeijm@users.noreply.github.com>
Date: Fri, 22 Aug 2025 13:31:00 +0200
Subject: [PATCH 07/21] refactor: clean up docstring

---
 src/ansys/health/heart/utils/compute_fiber_aha_angles.py | 6 +-----
 1 file changed, 1 insertion(+), 5 deletions(-)

diff --git a/src/ansys/health/heart/utils/compute_fiber_aha_angles.py b/src/ansys/health/heart/utils/compute_fiber_aha_angles.py
index 92b9fc727..2a446b4ed 100644
--- a/src/ansys/health/heart/utils/compute_fiber_aha_angles.py
+++ b/src/ansys/health/heart/utils/compute_fiber_aha_angles.py
@@ -20,11 +20,7 @@
 # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 # SOFTWARE.
 
-"""
-Compute average fiber orientations with respect to UHCs in each AHA region in the LV.
-
-from davoks
-"""
+"""Compute average fiber orientations with respect to UHCs in each AHA region in the LV."""
 
 import os
 import pathlib

From 153cc692ad9021a79f5ce53030a4d66f63d6d422 Mon Sep 17 00:00:00 2001
From: mhoeijm <102799582+mhoeijm@users.noreply.github.com>
Date: Fri, 22 Aug 2025 13:33:54 +0200
Subject: [PATCH 08/21] feat: add tests

---
 .../heart/utils/test_compute_fiber_angles.py  | 60 +++++++++++++++++++
 1 file changed, 60 insertions(+)
 create mode 100644 tests/heart/utils/test_compute_fiber_angles.py

diff --git a/tests/heart/utils/test_compute_fiber_angles.py b/tests/heart/utils/test_compute_fiber_angles.py
new file mode 100644
index 000000000..3e7005ea8
--- /dev/null
+++ b/tests/heart/utils/test_compute_fiber_angles.py
@@ -0,0 +1,60 @@
+# Copyright (C) 2023 - 2025 ANSYS, Inc. and/or its affiliates.
+# SPDX-License-Identifier: MIT
+#
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in all
+# copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.
+
+"""Test to verify the computation of fiber inclination angles."""
+
+import numpy as np
+
+from ansys.health.heart.utils.compute_fiber_aha_angles import signed_angle_between_vectors
+
+
+def test_compute_signed_angle():
+    """Test computing signed angles between vectors."""
+
+    # Local coordinate system
+    e_c = np.array([1.0, 0.0, 0.0])  # circumferential direction
+    e_t = np.array([0.0, 1.0, 0.0])  # transverse direction
+    e_l = np.array([0.0, 0.0, 1.0])  # longitudinal direction
+
+    # Test fiber directions
+    fiber = np.array([e_c, -e_c, e_l, e_t, -e_t, [1.0, 1.0, 1.0], [1.0, 0.0, -1.0]])
+    expected_angles = np.array(
+        [
+            0.0,  # angle between two equal vectors is 0
+            180.0,  # angle between opposite vectors is 180
+            90.0,  # angle between circumferential and longitudinal is 90 degrees
+            np.nan,  # undefined due to parallel with normal vector
+            np.nan,  # undefined due to parallel with normal vector
+            45.0,  # angle should be 45 degrees for the vector [1, 1, 1]
+            -45.0,  # angle should be -45 degrees for the vector [1, 0, -1]
+        ]
+    )
+
+    # Extend input vectors to right shape
+    num_fibers = fiber.shape[0]
+    e_c = np.tile(e_c, (num_fibers, 1))
+    e_t = np.tile(e_t, (num_fibers, 1))
+    e_l = np.tile(e_l, (num_fibers, 1))
+
+    computed_angles = signed_angle_between_vectors(x=fiber, y=e_c, n=e_t)
+
+    assert np.allclose(computed_angles, expected_angles, equal_nan=True)

From 01f3034d6e730b6b07f013853eafc38f50425e7c Mon Sep 17 00:00:00 2001
From: mhoeijm <102799582+mhoeijm@users.noreply.github.com>
Date: Fri, 22 Aug 2025 13:34:38 +0200
Subject: [PATCH 09/21] refactor: remove unused method

---
 .../health/heart/utils/compute_fiber_aha_angles.py  | 13 -------------
 1 file changed, 13 deletions(-)

diff --git a/src/ansys/health/heart/utils/compute_fiber_aha_angles.py b/src/ansys/health/heart/utils/compute_fiber_aha_angles.py
index 2a446b4ed..24a56d741 100644
--- a/src/ansys/health/heart/utils/compute_fiber_aha_angles.py
+++ b/src/ansys/health/heart/utils/compute_fiber_aha_angles.py
@@ -35,19 +35,6 @@
 from ansys.health.heart.simulator import DynaSettings, EPSimulator
 
 
-def angle_between_vectors(x: np.ndarray, y: np.ndarray) -> np.ndarray:
-    """Computes angles between N number of M-dimensional vectors
-    Input:  shape(x): (N,M)
-            shape(y): (N,M)
-    Output: shape(): (N,)
-    """  # noqa
-    return (
-        np.arccos(np.clip(np.sum(x * y, axis=1) / (norm(x, axis=1) * norm(y, axis=1)), -1, 1))
-        * 180
-        / np.pi
-    )
-
-
 def signed_angle_between_vectors(x: np.ndarray, y: np.ndarray, n: np.ndarray):
     """Compute signed angles between N number of M-dimensional vectors.
 

From 0c55914e716dfbe432f2b49247dce2dd8b67bcba Mon Sep 17 00:00:00 2001
From: mhoeijm <102799582+mhoeijm@users.noreply.github.com>
Date: Fri, 22 Aug 2025 13:35:13 +0200
Subject: [PATCH 10/21] refactor: move method to misc

---
 .../heart/utils/compute_fiber_aha_angles.py   | 55 +------------------
 src/ansys/health/heart/utils/misc.py          | 55 +++++++++++++++++++
 .../heart/utils/test_compute_fiber_angles.py  |  2 +-
 3 files changed, 57 insertions(+), 55 deletions(-)

diff --git a/src/ansys/health/heart/utils/compute_fiber_aha_angles.py b/src/ansys/health/heart/utils/compute_fiber_aha_angles.py
index 24a56d741..bb7282cb8 100644
--- a/src/ansys/health/heart/utils/compute_fiber_aha_angles.py
+++ b/src/ansys/health/heart/utils/compute_fiber_aha_angles.py
@@ -33,60 +33,7 @@
 
 import ansys.health.heart.models as models
 from ansys.health.heart.simulator import DynaSettings, EPSimulator
-
-
-def signed_angle_between_vectors(x: np.ndarray, y: np.ndarray, n: np.ndarray):
-    """Compute signed angles between N number of M-dimensional vectors.
-
-    Parameters
-    ----------
-    x : np.ndarray
-        (N,M) array for which to compute the angles.
-    y : np.ndarray
-        (N,M) array with reference vectors.
-    n : np.ndarray
-        (N,M) array with normal vectors of the reference plane.
-
-    Returns
-    -------
-    np.ndarray
-        Array of signed angles in degrees of size (N,).
-
-    Notes
-    -----
-    Computes signed angles between N number of M-dimensional vectors
-        Input:  shape(x): (N,M)
-                shape(y): (N,M)
-                shape(n): (N,M)
-        Output: shape(): (N,)
-
-    Va . Vb == |Va| * |Vb| * cos(alpha)    (by definition)
-            == |Va| * |Vb| * cos(beta)     (cos(alpha) == cos(-alpha) == cos(360° - alpha)
-
-
-    Va x Vb == |Va| * |Vb| * sin(alpha) * n1
-        (by definition; n1 is a unit vector perpendicular to Va and Vb with
-         orientation matching the right-hand rule)
-
-    Therefore (again assuming Vn is normalized):
-       n1 . Vn == 1 when beta < 180
-       n1 . Vn == -1 when beta > 180
-
-    ==>  (Va x Vb) . Vn == |Va| * |Vb| * sin(beta)
-
-    """
-    # normalize vectors
-    n /= norm(n, axis=1)[:, None]
-    x /= norm(x, axis=1)[:, None]
-    y /= norm(y, axis=1)[:, None]
-
-    nan_mask = np.isclose(np.linalg.norm(np.cross(x, n), axis=1), 0)
-
-    angles = np.arctan2(np.sum(np.cross(x, y) * n, axis=1), np.sum(x * y, axis=1))
-
-    angles[nan_mask] = np.nan  # Set angles where x is parallel to n to NaN
-
-    return angles * 180 / np.pi
+from ansys.health.heart.utils.misc import signed_angle_between_vectors
 
 
 def compute_aha_fiber_angles(mesh: pv.UnstructuredGrid, out_dir):
diff --git a/src/ansys/health/heart/utils/misc.py b/src/ansys/health/heart/utils/misc.py
index bba76f692..0e9726ce9 100644
--- a/src/ansys/health/heart/utils/misc.py
+++ b/src/ansys/health/heart/utils/misc.py
@@ -25,6 +25,7 @@
 import os
 
 import numpy as np
+from numpy.linalg import norm
 from scipy.spatial import cKDTree
 
 from ansys.health.heart import LOG as LOGGER
@@ -279,3 +280,57 @@ def interpolate_with_k_nearest(query_point: np.ndarray, k: int = 4) -> np.ndarra
         result[i] = interpolate_with_k_nearest(target_pos[i])
 
     return result
+
+
+def signed_angle_between_vectors(x: np.ndarray, y: np.ndarray, n: np.ndarray):
+    """Compute signed angles between N number of M-dimensional vectors.
+
+    Parameters
+    ----------
+    x : np.ndarray
+        (N,M) array for which to compute the angles.
+    y : np.ndarray
+        (N,M) array with reference vectors.
+    n : np.ndarray
+        (N,M) array with normal vectors of the reference plane.
+
+    Returns
+    -------
+    np.ndarray
+        Array of signed angles in degrees of size (N,).
+
+    Notes
+    -----
+    Computes signed angles between N number of M-dimensional vectors
+        Input:  shape(x): (N,M)
+                shape(y): (N,M)
+                shape(n): (N,M)
+        Output: shape(): (N,)
+
+    Va . Vb == |Va| * |Vb| * cos(alpha)    (by definition)
+            == |Va| * |Vb| * cos(beta)     (cos(alpha) == cos(-alpha) == cos(360° - alpha)
+
+
+    Va x Vb == |Va| * |Vb| * sin(alpha) * n1
+        (by definition; n1 is a unit vector perpendicular to Va and Vb with
+         orientation matching the right-hand rule)
+
+    Therefore (again assuming Vn is normalized):
+       n1 . Vn == 1 when beta < 180
+       n1 . Vn == -1 when beta > 180
+
+    ==>  (Va x Vb) . Vn == |Va| * |Vb| * sin(beta)
+
+    """
+    # normalize vectors
+    n /= norm(n, axis=1)[:, None]
+    x /= norm(x, axis=1)[:, None]
+    y /= norm(y, axis=1)[:, None]
+
+    nan_mask = np.isclose(np.linalg.norm(np.cross(x, n), axis=1), 0)
+
+    angles = np.arctan2(np.sum(np.cross(x, y) * n, axis=1), np.sum(x * y, axis=1))
+
+    angles[nan_mask] = np.nan  # Set angles where x is parallel to n to NaN
+
+    return angles * 180 / np.pi
diff --git a/tests/heart/utils/test_compute_fiber_angles.py b/tests/heart/utils/test_compute_fiber_angles.py
index 3e7005ea8..9eaac73be 100644
--- a/tests/heart/utils/test_compute_fiber_angles.py
+++ b/tests/heart/utils/test_compute_fiber_angles.py
@@ -24,7 +24,7 @@
 
 import numpy as np
 
-from ansys.health.heart.utils.compute_fiber_aha_angles import signed_angle_between_vectors
+from ansys.health.heart.utils.misc import signed_angle_between_vectors
 
 
 def test_compute_signed_angle():

From af962b142613618ca62b9fa40810fe686c5a603e Mon Sep 17 00:00:00 2001
From: mhoeijm <102799582+mhoeijm@users.noreply.github.com>
Date: Fri, 22 Aug 2025 13:48:43 +0200
Subject: [PATCH 11/21] doc: update docstring

---
 .../health/heart/utils/compute_fiber_aha_angles.py  | 13 ++++++++++---
 1 file changed, 10 insertions(+), 3 deletions(-)

diff --git a/src/ansys/health/heart/utils/compute_fiber_aha_angles.py b/src/ansys/health/heart/utils/compute_fiber_aha_angles.py
index bb7282cb8..1037085bb 100644
--- a/src/ansys/health/heart/utils/compute_fiber_aha_angles.py
+++ b/src/ansys/health/heart/utils/compute_fiber_aha_angles.py
@@ -37,9 +37,16 @@
 
 
 def compute_aha_fiber_angles(mesh: pv.UnstructuredGrid, out_dir):
-    """Computes average fiber inclination in each AHA region
-    as a function of the transmural depth
-    """  # noqa
+    """Compute the average fiber helix (inclination) and transverse angles.
+
+    Notes
+    -----
+    The helix (or inclination) angle is computed as the signed angle between the fiber direction
+    and the circumferential direction, with the transmural direction as the normal vector.
+
+    The transverse angle is computed as the signed angle between the fiber direction and the
+    transmural direction, with the longitudinal direction as the normal vector.
+    """
     expected_arrays = ["aha17", "fiber", "sheet", "transmural", "rotational"]
     for arr in expected_arrays:
         if arr not in mesh.cell_data and arr not in mesh.point_data:

From eed65c2546f6b82d7a9f0fc70ecf50b36105f854 Mon Sep 17 00:00:00 2001
From: pyansys-ci-bot <92810346+pyansys-ci-bot@users.noreply.github.com>
Date: Fri, 22 Aug 2025 11:49:15 +0000
Subject: [PATCH 12/21] chore: adding changelog file 827.added.md
 [dependabot-skip]

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

diff --git a/doc/source/changelog/827.added.md b/doc/source/changelog/827.added.md
new file mode 100644
index 000000000..3222de168
--- /dev/null
+++ b/doc/source/changelog/827.added.md
@@ -0,0 +1 @@
+Compute fiber helical angle

From 2c722214c129ae4401679400b3fa5e2eda80a32a Mon Sep 17 00:00:00 2001
From: mhoeijm <102799582+mhoeijm@users.noreply.github.com>
Date: Fri, 22 Aug 2025 15:02:36 +0200
Subject: [PATCH 13/21] refactor: rename module

---
 .../{compute_fiber_aha_angles.py => compute_fiber_angles.py}      | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename src/ansys/health/heart/utils/{compute_fiber_aha_angles.py => compute_fiber_angles.py} (100%)

diff --git a/src/ansys/health/heart/utils/compute_fiber_aha_angles.py b/src/ansys/health/heart/utils/compute_fiber_angles.py
similarity index 100%
rename from src/ansys/health/heart/utils/compute_fiber_aha_angles.py
rename to src/ansys/health/heart/utils/compute_fiber_angles.py

From eaec30741c5da31a0d0aed14c9bb3c9a44816096 Mon Sep 17 00:00:00 2001
From: mhoeijm <102799582+mhoeijm@users.noreply.github.com>
Date: Fri, 22 Aug 2025 15:03:23 +0200
Subject: [PATCH 14/21] refactor: rename test

---
 tests/heart/utils/test_compute_fiber_angles.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/heart/utils/test_compute_fiber_angles.py b/tests/heart/utils/test_compute_fiber_angles.py
index 9eaac73be..bcd0dc4af 100644
--- a/tests/heart/utils/test_compute_fiber_angles.py
+++ b/tests/heart/utils/test_compute_fiber_angles.py
@@ -27,7 +27,7 @@
 from ansys.health.heart.utils.misc import signed_angle_between_vectors
 
 
-def test_compute_signed_angle():
+def test_compute_signed_angles():
     """Test computing signed angles between vectors."""
 
     # Local coordinate system

From df0e395d8f45479ac83a5dd36f5ec13c01510669 Mon Sep 17 00:00:00 2001
From: mhoeijm <102799582+mhoeijm@users.noreply.github.com>
Date: Fri, 22 Aug 2025 15:04:09 +0200
Subject: [PATCH 15/21] refactor: extend edge case

---
 tests/heart/utils/test_compute_fiber_angles.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/tests/heart/utils/test_compute_fiber_angles.py b/tests/heart/utils/test_compute_fiber_angles.py
index bcd0dc4af..a2784801a 100644
--- a/tests/heart/utils/test_compute_fiber_angles.py
+++ b/tests/heart/utils/test_compute_fiber_angles.py
@@ -36,12 +36,13 @@ def test_compute_signed_angles():
     e_l = np.array([0.0, 0.0, 1.0])  # longitudinal direction
 
     # Test fiber directions
-    fiber = np.array([e_c, -e_c, e_l, e_t, -e_t, [1.0, 1.0, 1.0], [1.0, 0.0, -1.0]])
+    fiber = np.array([e_c, -e_c, e_l, 2 * e_l, e_t, -e_t, [1.0, 1.0, 1.0], [1.0, 0.0, -1.0]])
     expected_angles = np.array(
         [
             0.0,  # angle between two equal vectors is 0
             180.0,  # angle between opposite vectors is 180
             90.0,  # angle between circumferential and longitudinal is 90 degrees
+            90.0,  # angle between circumferential and longitudinal is 90 degrees
             np.nan,  # undefined due to parallel with normal vector
             np.nan,  # undefined due to parallel with normal vector
             45.0,  # angle should be 45 degrees for the vector [1, 1, 1]

From 8f77774ef3f7bdcf5a58aa22175dad2757d5b76f Mon Sep 17 00:00:00 2001
From: mhoeijm <102799582+mhoeijm@users.noreply.github.com>
Date: Mon, 25 Aug 2025 14:39:01 +0200
Subject: [PATCH 16/21] tests: add 4th quadrant angle

---
 tests/heart/utils/test_compute_fiber_angles.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/tests/heart/utils/test_compute_fiber_angles.py b/tests/heart/utils/test_compute_fiber_angles.py
index a2784801a..637f5113a 100644
--- a/tests/heart/utils/test_compute_fiber_angles.py
+++ b/tests/heart/utils/test_compute_fiber_angles.py
@@ -36,7 +36,9 @@ def test_compute_signed_angles():
     e_l = np.array([0.0, 0.0, 1.0])  # longitudinal direction
 
     # Test fiber directions
-    fiber = np.array([e_c, -e_c, e_l, 2 * e_l, e_t, -e_t, [1.0, 1.0, 1.0], [1.0, 0.0, -1.0]])
+    fiber = np.array(
+        [e_c, -e_c, e_l, 2 * e_l, e_t, -e_t, [1.0, 1.0, 1.0], [1.0, 0.0, -1.0], [-1.0, 0.0, -1.0]]
+    )
     expected_angles = np.array(
         [
             0.0,  # angle between two equal vectors is 0
@@ -47,6 +49,7 @@ def test_compute_signed_angles():
             np.nan,  # undefined due to parallel with normal vector
             45.0,  # angle should be 45 degrees for the vector [1, 1, 1]
             -45.0,  # angle should be -45 degrees for the vector [1, 0, -1]
+            -135.0,  # angle should be -135 degrees for the vector [1, 0, -1]
         ]
     )
 

From 5258bd2c63088b439e9b3fcb4cb42a7bdd3a370a Mon Sep 17 00:00:00 2001
From: mhoeijm <102799582+mhoeijm@users.noreply.github.com>
Date: Tue, 26 Aug 2025 15:32:11 +0200
Subject: [PATCH 17/21] feat: add projection to angle computation

---
 src/ansys/health/heart/utils/misc.py          | 57 ++++++++++---------
 .../heart/utils/test_compute_fiber_angles.py  | 53 ++++++++++++++---
 2 files changed, 74 insertions(+), 36 deletions(-)

diff --git a/src/ansys/health/heart/utils/misc.py b/src/ansys/health/heart/utils/misc.py
index 0e9726ce9..59778768d 100644
--- a/src/ansys/health/heart/utils/misc.py
+++ b/src/ansys/health/heart/utils/misc.py
@@ -23,6 +23,7 @@
 """Module containing miscellaneous methods."""
 
 import os
+from typing import Literal
 
 import numpy as np
 from numpy.linalg import norm
@@ -282,17 +283,26 @@ def interpolate_with_k_nearest(query_point: np.ndarray, k: int = 4) -> np.ndarra
     return result
 
 
-def signed_angle_between_vectors(x: np.ndarray, y: np.ndarray, n: np.ndarray):
-    """Compute signed angles between N number of M-dimensional vectors.
+def angle_between_vectors(
+    x: np.ndarray,
+    y: np.ndarray,
+    n: np.ndarray,
+    representation: Literal["2-quadrant-signed", "4-quadrant"] = "2-quadrant-signed",
+):
+    """Compute signed angles from between N number of M-dimensional vectors.
 
     Parameters
     ----------
     x : np.ndarray
-        (N,M) array for which to compute the angles.
-    y : np.ndarray
         (N,M) array with reference vectors.
+    y : np.ndarray
+        (N,M) array with vectors for which to compute the angle.
     n : np.ndarray
-        (N,M) array with normal vectors of the reference plane.
+        (N,M) array with normal vectors of the plane on which to project y.
+    representation : Literal["2-quadrant-signed", "4-quadrant"], default: "2-quadrant-signed"
+        Representation of the angle:
+            - "2-quadrant-signed": angles in the range [-180, 180)
+            - "4-quadrant": angles in the range [0, 360)
 
     Returns
     -------
@@ -301,25 +311,8 @@ def signed_angle_between_vectors(x: np.ndarray, y: np.ndarray, n: np.ndarray):
 
     Notes
     -----
-    Computes signed angles between N number of M-dimensional vectors
-        Input:  shape(x): (N,M)
-                shape(y): (N,M)
-                shape(n): (N,M)
-        Output: shape(): (N,)
-
-    Va . Vb == |Va| * |Vb| * cos(alpha)    (by definition)
-            == |Va| * |Vb| * cos(beta)     (cos(alpha) == cos(-alpha) == cos(360° - alpha)
-
-
-    Va x Vb == |Va| * |Vb| * sin(alpha) * n1
-        (by definition; n1 is a unit vector perpendicular to Va and Vb with
-         orientation matching the right-hand rule)
-
-    Therefore (again assuming Vn is normalized):
-       n1 . Vn == 1 when beta < 180
-       n1 . Vn == -1 when beta > 180
-
-    ==>  (Va x Vb) . Vn == |Va| * |Vb| * sin(beta)
+    Projects y vectors onto the planes defined by the n vectors and then computes the angle
+    between the x and the projected y vectors.
 
     """
     # normalize vectors
@@ -329,8 +322,18 @@ def signed_angle_between_vectors(x: np.ndarray, y: np.ndarray, n: np.ndarray):
 
     nan_mask = np.isclose(np.linalg.norm(np.cross(x, n), axis=1), 0)
 
-    angles = np.arctan2(np.sum(np.cross(x, y) * n, axis=1), np.sum(x * y, axis=1))
+    # project y onto n
+    y_proj = y - np.sum(y * n, axis=1, keepdims=True) * n
+    y_proj /= norm(y_proj, axis=1)[:, None]
+
+    # compute angle between projected y and x, positive angle following right-hand rule around n
+    angles = np.arctan2(np.sum(np.cross(x, y_proj) * n, axis=1), np.sum(x * y_proj, axis=1))
+
+    angles[nan_mask] = np.nan  # Set angles where y is parallel to n to NaN
 
-    angles[nan_mask] = np.nan  # Set angles where x is parallel to n to NaN
+    angles = np.degrees(angles)
 
-    return angles * 180 / np.pi
+    if representation == "4-quadrant":
+        return np.mod(angles + 360, 360)
+    else:
+        return angles
diff --git a/tests/heart/utils/test_compute_fiber_angles.py b/tests/heart/utils/test_compute_fiber_angles.py
index 637f5113a..cfcb18412 100644
--- a/tests/heart/utils/test_compute_fiber_angles.py
+++ b/tests/heart/utils/test_compute_fiber_angles.py
@@ -24,20 +24,33 @@
 
 import numpy as np
 
-from ansys.health.heart.utils.misc import signed_angle_between_vectors
+from ansys.health.heart.utils.misc import angle_between_vectors
 
 
 def test_compute_signed_angles():
     """Test computing signed angles between vectors."""
 
-    # Local coordinate system
-    e_c = np.array([1.0, 0.0, 0.0])  # circumferential direction
-    e_t = np.array([0.0, 1.0, 0.0])  # transverse direction
-    e_l = np.array([0.0, 0.0, 1.0])  # longitudinal direction
+    # Local element coordinate system
+    # <e_t, e_c, e_l> transmural, circumferential, longitudinal directions
+    e_t = np.array([1.0, 0.0, 0.0])  # transmural direction
+    e_c = np.array([0.0, 1.0, 0.0])  # circumferential direction
+    e_l = np.array([0.0, 0.0, 1.0])  # longitudinal direction e_l = e_t X e_c
 
     # Test fiber directions
     fiber = np.array(
-        [e_c, -e_c, e_l, 2 * e_l, e_t, -e_t, [1.0, 1.0, 1.0], [1.0, 0.0, -1.0], [-1.0, 0.0, -1.0]]
+        [
+            e_c,
+            -e_c,
+            e_l,
+            2 * e_l,
+            e_t,
+            -e_t,
+            [1.0, 1.0, 1.0],
+            [1.0, 0.0, -1.0],
+            [-1.0, 0.0, -1.0],
+            [0.0, -1.0, 1.0],
+            [0.0, -1.0, -1.0],
+        ]
     )
     expected_angles = np.array(
         [
@@ -48,8 +61,10 @@ def test_compute_signed_angles():
             np.nan,  # undefined due to parallel with normal vector
             np.nan,  # undefined due to parallel with normal vector
             45.0,  # angle should be 45 degrees for the vector [1, 1, 1]
-            -45.0,  # angle should be -45 degrees for the vector [1, 0, -1]
-            -135.0,  # angle should be -135 degrees for the vector [1, 0, -1]
+            -90.0,  # angle should be -90 degrees for the vector [1, 0, -1]
+            -90.0,  # angle should be -90 degrees for the vector [1, 0, -1]
+            135.0,  # angle should be -45 degrees for the vector [0, -1, 1]
+            -135.0,  # angle should be -135 degrees for the vector [0, -1, -1]
         ]
     )
 
@@ -59,6 +74,26 @@ def test_compute_signed_angles():
     e_t = np.tile(e_t, (num_fibers, 1))
     e_l = np.tile(e_l, (num_fibers, 1))
 
-    computed_angles = signed_angle_between_vectors(x=fiber, y=e_c, n=e_t)
+    computed_angles = angle_between_vectors(
+        x=e_c, y=fiber, n=e_t, representation="2-quadrant-signed"
+    )
+
+    assert np.allclose(computed_angles, expected_angles, equal_nan=True)
 
+    computed_angles = angle_between_vectors(x=e_c, y=fiber, n=e_t, representation="4-quadrant")
+    expected_angles = np.array(
+        [
+            0.0,  # angle between two equal vectors is 0
+            180.0,  # angle between opposite vectors is 180
+            90.0,  # angle between circumferential and longitudinal is 90 degrees
+            90.0,  # angle between circumferential and longitudinal is 90 degrees
+            np.nan,  # undefined due to parallel with normal vector
+            np.nan,  # undefined due to parallel with normal vector
+            45.0,  # angle should be 45 degrees for the vector [1, 1, 1]
+            270,
+            270,
+            135.0,  # angle should be 315 degrees for the vector [1, 0, -1]
+            225.0,  # angle should be 225 degrees for the vector [1, 0, -1]
+        ]
+    )
     assert np.allclose(computed_angles, expected_angles, equal_nan=True)

From 2cb58cb5eec00f287addac0018a5f328fa0808cc Mon Sep 17 00:00:00 2001
From: mhoeijm <102799582+mhoeijm@users.noreply.github.com>
Date: Wed, 27 Aug 2025 14:04:38 +0200
Subject: [PATCH 18/21] docs: fix docstring

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

diff --git a/src/ansys/health/heart/utils/misc.py b/src/ansys/health/heart/utils/misc.py
index 59778768d..42ae7e3a0 100644
--- a/src/ansys/health/heart/utils/misc.py
+++ b/src/ansys/health/heart/utils/misc.py
@@ -289,14 +289,14 @@ def angle_between_vectors(
     n: np.ndarray,
     representation: Literal["2-quadrant-signed", "4-quadrant"] = "2-quadrant-signed",
 ):
-    """Compute signed angles from between N number of M-dimensional vectors.
+    """Compute signed angles between N number of M-dimensional vectors.
 
     Parameters
     ----------
     x : np.ndarray
-        (N,M) array with reference vectors.
+        (N,M) array with first (reference) vectors.
     y : np.ndarray
-        (N,M) array with vectors for which to compute the angle.
+        (N,M) array with second vectors. The angle is computed from vectors x to vectors y.
     n : np.ndarray
         (N,M) array with normal vectors of the plane on which to project y.
     representation : Literal["2-quadrant-signed", "4-quadrant"], default: "2-quadrant-signed"

From f3f66409f314cf80b5bf978d778f79a4a55405d2 Mon Sep 17 00:00:00 2001
From: mhoeijm <102799582+mhoeijm@users.noreply.github.com>
Date: Wed, 27 Aug 2025 14:59:31 +0200
Subject: [PATCH 19/21] feat: improve and cleanup fiber angle

---
 .../heart/utils/compute_fiber_angles.py       | 280 ++++++++++--------
 1 file changed, 158 insertions(+), 122 deletions(-)

diff --git a/src/ansys/health/heart/utils/compute_fiber_angles.py b/src/ansys/health/heart/utils/compute_fiber_angles.py
index 1037085bb..f8c8714d3 100644
--- a/src/ansys/health/heart/utils/compute_fiber_angles.py
+++ b/src/ansys/health/heart/utils/compute_fiber_angles.py
@@ -23,7 +23,6 @@
 """Compute average fiber orientations with respect to UHCs in each AHA region in the LV."""
 
 import os
-import pathlib
 
 import matplotlib.pyplot as plt
 import numpy as np
@@ -31,22 +30,71 @@
 import pandas as pd
 import pyvista as pv
 
-import ansys.health.heart.models as models
-from ansys.health.heart.simulator import DynaSettings, EPSimulator
-from ansys.health.heart.utils.misc import signed_angle_between_vectors
+from ansys.health.heart.utils.misc import angle_between_vectors
 
 
-def compute_aha_fiber_angles(mesh: pv.UnstructuredGrid, out_dir):
-    """Compute the average fiber helix (inclination) and transverse angles.
+def compute_gradient(mesh: pv.UnstructuredGrid, scalar: str, interpolate_to_cell: bool = False):
+    """Compute gradient of a scalar field on the mesh.
 
-    Notes
-    -----
-    The helix (or inclination) angle is computed as the signed angle between the fiber direction
-    and the circumferential direction, with the transmural direction as the normal vector.
+    Parameters
+    ----------
+    mesh : pv.UnstructuredGrid
+        Input mesh.
+    scalar : str
+        Name of the scalar field.
 
-    The transverse angle is computed as the signed angle between the fiber direction and the
-    transmural direction, with the longitudinal direction as the normal vector.
+    Returns
+    -------
+    gradient : np.ndarray
+        Gradient of the scalar field at each point.
     """
+    if scalar not in mesh.point_data:
+        raise ValueError(f"Scalar '{scalar}' not found in mesh point data.")
+
+    # compute gradient at points
+    mesh_with_grad = mesh.compute_derivative(scalars=scalar, preference="point")
+
+    if interpolate_to_cell:
+        grad = mesh_with_grad.point_data_to_cell_data().cell_data["gradient"]
+    else:
+        grad = mesh_with_grad.point_data["gradient"]
+    return grad
+
+
+def _get_angles_from_mesh(grid: pv.UnstructuredGrid) -> pd.DataFrame:
+    """Retrieve the fiber angles from the mesh and store in pandas dataframe."""
+    expected_arrays = ["aha17", "helix_angle", "transverse_angle", "depth_bin"]
+    for arr in expected_arrays:
+        if arr not in grid.cell_data and arr not in grid.point_data:
+            raise ValueError(f"Expected array '{arr}' not found in mesh data.")
+
+    # Retrieve angles for each AHA segment and depth bin
+    ndepths = np.unique(grid["depth_bin"]).shape[0]
+    n_ahas = np.unique(grid["aha17"]).shape[0]
+    aha_ids = np.unique(grid["aha17"])
+
+    data_helix = np.full((ndepths, n_ahas), np.nan)
+    data_transverse = np.full((ndepths, n_ahas), np.nan)
+
+    # Retrieve helix and transverse angles for each AHA segment and depth bin
+    for ii in range(1, 18):
+        for jj in range(1, ndepths + 1):
+            mask = grid["aha17"] == ii
+            mask &= grid["depth_bin"] == jj
+
+            helix_angle = np.nanmean(grid.extract_cells(mask)["helix_angle"])
+            transverse_angle = np.nanmean(grid.extract_cells(mask)["transverse_angle"])
+
+            data_helix[jj - 1, ii - 1] = helix_angle  # depth in row, AHA in column
+            data_transverse[jj - 1, ii - 1] = transverse_angle  # depth in row, AHA in column
+
+    legend_entries = [f"AHA_{i}" for i in aha_ids]
+
+    return data_helix, data_transverse, legend_entries
+
+
+def compute_aha_fiber_angles1(mesh: pv.UnstructuredGrid, out_dir: str):
+    """Compute the average fiber helix and transverse angles."""
     expected_arrays = ["aha17", "fiber", "sheet", "transmural", "rotational"]
     for arr in expected_arrays:
         if arr not in mesh.cell_data and arr not in mesh.point_data:
@@ -54,37 +102,29 @@ def compute_aha_fiber_angles(mesh: pv.UnstructuredGrid, out_dir):
 
     aha_ids = mesh["aha17"]
     aha_elements = np.where(~np.isnan(aha_ids))[0]
-    aha_model = mesh.extract_cells(aha_elements)
-    aha_model.cell_data["AHA"] = aha_ids[aha_elements]
+    aha_model: pv.UnstructuredGrid = mesh.extract_cells(aha_elements)
+    # aha_model.cell_data["fiber"] = aha_model.cell_data["fiber"] * -1  # flip fiber direction
+    aha_ids = aha_model.cell_data["aha17"]
 
-    # print("Nmbr of points:\t{:d}".format(mesh.n_points))
-    # print("Nmbr of cells:\t{:d}".format(mesh.n_cells))
-    # print("Nmbr of AHA cells:\t{:d}".format(len(aha_elements)))
+    # flip fibers
+    # aha_model.cell_data["fiber"] = 1 * aha_model.cell_data["fiber"]
 
     # load fibers and sheets at cells
-    el_fibers = mesh.cell_data["fiber"]
-    el_sheets = mesh.cell_data["sheet"]
-    el_fibers = el_fibers[aha_elements]
-    el_sheets = el_sheets[aha_elements]
-    el_fibers /= np.linalg.norm(el_fibers, axis=1)[:, None]  # make sure fibers are normalized
-
-    # interpolate transmural depth from points to cells
-    el_depths = mesh.point_data_to_cell_data()["transmural"][aha_elements]
-    el_depths = 2.0 * el_depths - 1.0  # map from [0,1] -> [-1,1]
-    # interpolate rotational coordinates from points to cells
-    # el_rotat = mesh.point_data_to_cell_data()["rotational"][aha_elements]
-
-    # compute transmural vector from derivative of transmural depth
+    el_fibers = aha_model.cell_data["fiber"]
+    el_sheets = aha_model.cell_data["sheet"]  # noqa N841
+
     # TODO : interpolate t instead of grad_t
-    pt_grad_t = mesh.compute_derivative(scalars="transmural", preference="cell")["gradient"]
-    mesh.point_data.set_scalars(name="grad_t", scalars=pt_grad_t)
-    el_grad_t = mesh.point_data_to_cell_data()["grad_t"][aha_elements]
-    el_grad_t = el_grad_t / np.linalg.norm(el_grad_t, axis=1)[:, None]
-
-    # compute transmural vector from derivative of rotational coordinate
-    pt_grad_r = mesh.compute_derivative(scalars="rotational", preference="cell")["gradient"]
-    mesh.point_data.set_scalars(name="grad_r", scalars=pt_grad_r)
-    el_grad_r = mesh.point_data_to_cell_data()["grad_r"][aha_elements]
+    aha_model.cell_data["grad_transmural"] = compute_gradient(
+        aha_model, "transmural", interpolate_to_cell=True
+    )
+    el_grad_t = aha_model.cell_data["grad_transmural"]
+
+    # compute rotational vector from derivative of rotational coordinate
+    aha_model.cell_data["grad_rotational"] = compute_gradient(
+        aha_model, "rotational", interpolate_to_cell=True
+    )
+    el_grad_r = aha_model.cell_data["grad_rotational"]
+
     # set elements at the rotational coordinate discontinuity to NaN
     # TODO: prescribe to average gradient from nearest neighbors
     norm_grad_r = norm(el_grad_r, axis=1)
@@ -94,59 +134,93 @@ def compute_aha_fiber_angles(mesh: pv.UnstructuredGrid, out_dir):
     el_grad_r[id_discont, :] = nans
     el_grad_r = el_grad_r / np.linalg.norm(el_grad_r, axis=1)[:, None]
 
-    # compute angle between fibers and transmural vectors
-    el_angles_t = signed_angle_between_vectors(el_fibers, el_grad_t, el_grad_r)
+    aha_model.cell_data["grad_rotational"] = el_grad_r
+
+    # visualize: rotational vector
+    # plotter = pv.Plotter()
+    # aha_model.set_active_scalars("aha17")
+    # glyph = aha_model.glyph("grad_transmural", scale=False)
+    # plotter.add_mesh(aha_model, opacity=0.1, color="white")
+    # plotter.add_mesh(glyph, color="blue")
+    # plotter.show()
 
-    # compute angle between fibers and rotational vectors
-    el_angles_r = signed_angle_between_vectors(el_grad_r, el_fibers, el_grad_t)
+    # compute longitudinal vector as cross product of transmural and rotational vectors
+    el_grad_l = np.cross(el_grad_t, el_grad_r)
+    aha_model.cell_data["grad_longitudinal"] = el_grad_l
 
-    # get aha17 label for left ventricle elements
-    aha17_label = aha_ids[aha_elements]
+    # compute angle between rotational and fiber vectors in plane defined
+    # by longitudinal vector (transverse angle)
 
-    # average fiber angles wrt to transmural and rotational coords in each AHA and depth
+    el_angles_t = angle_between_vectors(el_grad_r, el_fibers, el_grad_l, "2-quadrant-signed")
+    aha_model.cell_data["transverse_angle"] = el_angles_t
+
+    # compute angle between fibers and rotational vectors in plane defined
+    # by transmural vector (helix angle)
+    el_angles_r = angle_between_vectors(el_grad_r, el_fibers, el_grad_t, "2-quadrant-signed")
+    aha_model.cell_data["helix_angle"] = el_angles_r
+
+    # Add depth bins to the mesh
     ndepths = 9
-    depth_bin_edges = np.linspace(-1.0, 1.0, ndepths + 1)
-    depth_bin_ctrs = 0.5 * (depth_bin_edges[:-1] + depth_bin_edges[1:])
-    el_angles_t_avg = np.zeros((len(aha_elements)))
-    aha_angles_t = np.zeros((17, ndepths))
-    el_angles_r_avg = np.zeros((len(aha_elements)))
-    aha_angles_r = np.zeros((17, ndepths))
-    for i in range(1, 18):
-        idx_aha = aha17_label == i
-        for j in range(ndepths):
-            depth_min = depth_bin_edges[j]
-            depth_max = depth_bin_edges[j + 1]
-            idx_depth = (el_depths >= depth_min) & (el_depths < depth_max)
-            idx = np.where(idx_aha & idx_depth)[0]
-            aha_angles_t[i - 1, j] = np.nanmean(el_angles_t[idx])
-            el_angles_t_avg[idx] = aha_angles_t[i - 1, j]
-            aha_angles_r[i - 1, j] = np.nanmean(el_angles_r[idx])
-            el_angles_r_avg[idx] = aha_angles_r[i - 1, j]
-
-    # save to vtk
-    aha_model.cell_data["transmural_angles"] = el_angles_t
-    aha_model.cell_data["transmural_angles_aha"] = el_angles_t_avg
-    aha_model.cell_data["rotational_angles"] = el_angles_r
-    aha_model.cell_data["rotational_angles_aha"] = el_angles_r_avg
-    aha_model.cell_data["transmural_depth"] = el_depths
-    aha_model.cell_data.set_vectors(el_fibers, "fibers", deep_copy=True)
-    aha_model.cell_data.set_vectors(el_sheets, "sheets", deep_copy=True)
-    aha_model.cell_data.set_vectors(el_grad_t, "grad_t", deep_copy=True)
-    aha_model.cell_data.set_vectors(el_grad_r, "grad_r", deep_copy=True)
-    aha_model.save(os.path.join(pathlib.Path(out_dir), "aha_averaged_angles.vtk"))
+    depth_lower_limit = np.linspace(-1.01, 1.0, ndepths + 1)[0:-1]
+    depth_upper_limit = np.append(depth_lower_limit[1:], 1.01)
+    depth_bin_ctrs = np.mean(np.vstack([depth_lower_limit, depth_upper_limit]), axis=0)
+
+    # Compute the depth of each cell in the mesh, normalized to [-1 (endocardium), 1 (epicardium)]
+    depth_cells = aha_model.point_data_to_cell_data()["transmural"]
+    aha_model.cell_data["depth"] = 2 * depth_cells - 1.0
+
+    # Split in 9 transmural bins and store bin number in cell data
+    aha_model.cell_data["depth_bin"] = 0.0
+    depth_bin = 1
+    for lower, upper in zip(depth_lower_limit, depth_upper_limit):
+        mask = aha_model["depth"] > lower
+        mask &= aha_model["depth"] < upper
+
+        aha_model.cell_data["depth_bin"][mask] = float(depth_bin)
+        depth_bin += 1
+
+    if np.any(aha_model.cell_data["depth_bin"] == 0):
+        raise ValueError("Some cells were not assigned to a depth bin - check tolerances.")
+
+    data_helix, data_transverse, legend_entries = _get_angles_from_mesh(aha_model)
 
     # save to csv
     cols = ["{:1.2f}".format(x) for x in depth_bin_ctrs]
     rows = ["{:d}".format(x) for x in range(1, 18)]
-    df_r = pd.DataFrame(data=aha_angles_r, index=rows, columns=cols)
-    df_r.to_csv(os.path.join(out_dir, "AHA_fiber_angles_r.csv"), index=True)
-    df_t = pd.DataFrame(data=aha_angles_t, index=rows, columns=cols)
-    df_t.to_csv(os.path.join(out_dir, "AHA_fiber_angles_t.csv"), index=True)
-
-    # print(df_r)
-    # print(df_t)
-
-    return df_r, df_t
+    df_helix_angle = pd.DataFrame(data=data_helix.T, index=rows, columns=cols)
+    df_helix_angle.to_csv(os.path.join(out_dir, "AHA_fiber_angles_r.csv"), index=True)
+    df_transverse_angle = pd.DataFrame(data=data_transverse.T, index=rows, columns=cols)
+    df_transverse_angle.to_csv(os.path.join(out_dir, "AHA_fiber_angles_t.csv"), index=True)
+
+    # fig, axs = plt.subplots(1, 2)
+    # axs[0].plot(depth_bin_ctrs, data_helix, "o-")
+    # axs[0].set_title(r"Helix Angle alpha_h")
+    # axs[0].set_xlim([-1, 1])
+    # axs[0].set_xlabel("transmural depth")
+    # axs[0].set_ylabel("angle [deg]")
+    # axs[1].plot(depth_bin_ctrs, data_transverse, "o-")
+    # axs[1].set_title(r"Transverse Angle alpha_t")
+    # axs[1].set_xlabel("transmural depth [-]")
+    # axs[1].set_xlim([-1, 1])
+    # axs[0].legend(legend_entries)
+    # fig.show()
+
+    # aha_model.save(os.path.join(out_dir, "aha_model.vtu"))
+
+    # sub_model: pv.UnstructuredGrid = aha_model.threshold([1, 1], "aha17").threshold(
+    #     [1, 1], "depth_bin"
+    # )
+
+    # # plot components: transmural, rotational and longitudinal vectors
+    # plotter = pv.Plotter()
+    # plotter.add_mesh(sub_model.glyph("grad_transmural", scale=False), color="b")
+    # plotter.add_mesh(sub_model.glyph("grad_rotational", scale=False), color="r")
+    # plotter.add_mesh(sub_model.glyph("grad_longitudinal", scale=False), color="g")
+    # plotter.add_mesh(sub_model.glyph("fiber", scale=False), color="white")
+    # plotter.add_mesh(aha_model, opacity=0.5, color="white")
+    # plotter.show()
+
+    return df_helix_angle, df_transverse_angle
 
 
 def plot_fiber_aha_angles(data: pd.DataFrame | str):
@@ -205,41 +279,3 @@ def plot_fiber_aha_angles(data: pd.DataFrame | str):
 
     # plt.savefig("fiber_helical_angles.png", bbox_inches="tight")
     plt.show()
-
-
-if __name__ == "__main__":
-    """Demo
-    """
-    # assumed LS-DYNA
-    lsdyna_path = r"D:\ansysdev\lsdyna\ls-dyna_smp_d_DEV_112901-gcbb8e36701_winx64_ifort190.exe"
-    dyna_settings = DynaSettings(
-        lsdyna_path=lsdyna_path, dynatype="smp", num_cpus=4, platform="windows"
-    )
-
-    # assumed case
-    workdir = os.path.abspath(os.path.join("downloads", "Rodero2021", "01", "BV"))
-    path_to_model = os.path.join(workdir, "heart_model.vtu")
-
-    model: models.BiVentricle = models.BiVentricle(workdir)
-    model.load_model_from_mesh(path_to_model, path_to_model.replace(".vtu", ".partinfo.json"))
-
-    # get fields data
-    simulator = EPSimulator(
-        model=model,
-        dyna_settings=dyna_settings,
-        simulation_directory=os.path.join(workdir, "simulation-EP"),
-    )
-    simulator.settings.load_defaults()
-    simulator.compute_uhc()
-    simulator.compute_fibers()
-
-    # get AHA labels
-    from ansys.health.heart.utils.landmark_utils import compute_aha17
-
-    aha_ids = compute_aha17(simulator.model, simulator.model.short_axis, simulator.model.l4cv_axis)
-    grid = simulator.model.mesh.extract_cells_by_type(10)
-    grid.cell_data["aha17"] = aha_ids
-
-    # compute angles
-    df_r, df_t = compute_aha_fiber_angles(grid, workdir)
-    plot_fiber_aha_angles(df_t)

From ef337679d82fe68e06d8c17ef7046788082c062a Mon Sep 17 00:00:00 2001
From: mhoeijm <102799582+mhoeijm@users.noreply.github.com>
Date: Wed, 27 Aug 2025 15:25:28 +0200
Subject: [PATCH 20/21] refactor: cleanup

---
 .../heart/utils/compute_fiber_angles.py       | 52 ++++---------------
 1 file changed, 9 insertions(+), 43 deletions(-)

diff --git a/src/ansys/health/heart/utils/compute_fiber_angles.py b/src/ansys/health/heart/utils/compute_fiber_angles.py
index f8c8714d3..ea9177d3c 100644
--- a/src/ansys/health/heart/utils/compute_fiber_angles.py
+++ b/src/ansys/health/heart/utils/compute_fiber_angles.py
@@ -33,7 +33,7 @@
 from ansys.health.heart.utils.misc import angle_between_vectors
 
 
-def compute_gradient(mesh: pv.UnstructuredGrid, scalar: str, interpolate_to_cell: bool = False):
+def _compute_gradient(mesh: pv.UnstructuredGrid, scalar: str, interpolate_to_cell: bool = False):
     """Compute gradient of a scalar field on the mesh.
 
     Parameters
@@ -93,7 +93,7 @@ def _get_angles_from_mesh(grid: pv.UnstructuredGrid) -> pd.DataFrame:
     return data_helix, data_transverse, legend_entries
 
 
-def compute_aha_fiber_angles1(mesh: pv.UnstructuredGrid, out_dir: str):
+def compute_aha_fiber_angles(mesh: pv.UnstructuredGrid, out_dir: str):
     """Compute the average fiber helix and transverse angles."""
     expected_arrays = ["aha17", "fiber", "sheet", "transmural", "rotational"]
     for arr in expected_arrays:
@@ -114,13 +114,13 @@ def compute_aha_fiber_angles1(mesh: pv.UnstructuredGrid, out_dir: str):
     el_sheets = aha_model.cell_data["sheet"]  # noqa N841
 
     # TODO : interpolate t instead of grad_t
-    aha_model.cell_data["grad_transmural"] = compute_gradient(
+    aha_model.cell_data["grad_transmural"] = _compute_gradient(
         aha_model, "transmural", interpolate_to_cell=True
     )
     el_grad_t = aha_model.cell_data["grad_transmural"]
 
     # compute rotational vector from derivative of rotational coordinate
-    aha_model.cell_data["grad_rotational"] = compute_gradient(
+    aha_model.cell_data["grad_rotational"] = _compute_gradient(
         aha_model, "rotational", interpolate_to_cell=True
     )
     el_grad_r = aha_model.cell_data["grad_rotational"]
@@ -136,14 +136,6 @@ def compute_aha_fiber_angles1(mesh: pv.UnstructuredGrid, out_dir: str):
 
     aha_model.cell_data["grad_rotational"] = el_grad_r
 
-    # visualize: rotational vector
-    # plotter = pv.Plotter()
-    # aha_model.set_active_scalars("aha17")
-    # glyph = aha_model.glyph("grad_transmural", scale=False)
-    # plotter.add_mesh(aha_model, opacity=0.1, color="white")
-    # plotter.add_mesh(glyph, color="blue")
-    # plotter.show()
-
     # compute longitudinal vector as cross product of transmural and rotational vectors
     el_grad_l = np.cross(el_grad_t, el_grad_r)
     aha_model.cell_data["grad_longitudinal"] = el_grad_l
@@ -173,8 +165,8 @@ def compute_aha_fiber_angles1(mesh: pv.UnstructuredGrid, out_dir: str):
     aha_model.cell_data["depth_bin"] = 0.0
     depth_bin = 1
     for lower, upper in zip(depth_lower_limit, depth_upper_limit):
-        mask = aha_model["depth"] > lower
-        mask &= aha_model["depth"] < upper
+        mask = aha_model.cell_data["depth"] > lower
+        mask &= aha_model.cell_data["depth"] < upper
 
         aha_model.cell_data["depth_bin"][mask] = float(depth_bin)
         depth_bin += 1
@@ -192,35 +184,9 @@ def compute_aha_fiber_angles1(mesh: pv.UnstructuredGrid, out_dir: str):
     df_transverse_angle = pd.DataFrame(data=data_transverse.T, index=rows, columns=cols)
     df_transverse_angle.to_csv(os.path.join(out_dir, "AHA_fiber_angles_t.csv"), index=True)
 
-    # fig, axs = plt.subplots(1, 2)
-    # axs[0].plot(depth_bin_ctrs, data_helix, "o-")
-    # axs[0].set_title(r"Helix Angle alpha_h")
-    # axs[0].set_xlim([-1, 1])
-    # axs[0].set_xlabel("transmural depth")
-    # axs[0].set_ylabel("angle [deg]")
-    # axs[1].plot(depth_bin_ctrs, data_transverse, "o-")
-    # axs[1].set_title(r"Transverse Angle alpha_t")
-    # axs[1].set_xlabel("transmural depth [-]")
-    # axs[1].set_xlim([-1, 1])
-    # axs[0].legend(legend_entries)
-    # fig.show()
-
-    # aha_model.save(os.path.join(out_dir, "aha_model.vtu"))
-
-    # sub_model: pv.UnstructuredGrid = aha_model.threshold([1, 1], "aha17").threshold(
-    #     [1, 1], "depth_bin"
-    # )
-
-    # # plot components: transmural, rotational and longitudinal vectors
-    # plotter = pv.Plotter()
-    # plotter.add_mesh(sub_model.glyph("grad_transmural", scale=False), color="b")
-    # plotter.add_mesh(sub_model.glyph("grad_rotational", scale=False), color="r")
-    # plotter.add_mesh(sub_model.glyph("grad_longitudinal", scale=False), color="g")
-    # plotter.add_mesh(sub_model.glyph("fiber", scale=False), color="white")
-    # plotter.add_mesh(aha_model, opacity=0.5, color="white")
-    # plotter.show()
-
-    return df_helix_angle, df_transverse_angle
+    aha_model.save(os.path.join(out_dir, "aha_model.vtu"))
+
+    return df_helix_angle, df_transverse_angle, aha_model
 
 
 def plot_fiber_aha_angles(data: pd.DataFrame | str):

From 0d85743f3a5604ff680e169a3471299ca56655eb Mon Sep 17 00:00:00 2001
From: mhoeijm <102799582+mhoeijm@users.noreply.github.com>
Date: Thu, 28 Aug 2025 12:46:08 +0200
Subject: [PATCH 21/21] refactor: move to private method and add proper ylabel

---
 .../heart/utils/compute_fiber_angles.py       | 46 ++++++++++++++++---
 1 file changed, 39 insertions(+), 7 deletions(-)

diff --git a/src/ansys/health/heart/utils/compute_fiber_angles.py b/src/ansys/health/heart/utils/compute_fiber_angles.py
index ea9177d3c..d5eae28d9 100644
--- a/src/ansys/health/heart/utils/compute_fiber_angles.py
+++ b/src/ansys/health/heart/utils/compute_fiber_angles.py
@@ -23,6 +23,7 @@
 """Compute average fiber orientations with respect to UHCs in each AHA region in the LV."""
 
 import os
+from typing import Literal
 
 import matplotlib.pyplot as plt
 import numpy as np
@@ -93,22 +94,44 @@ def _get_angles_from_mesh(grid: pv.UnstructuredGrid) -> pd.DataFrame:
     return data_helix, data_transverse, legend_entries
 
 
-def compute_aha_fiber_angles(mesh: pv.UnstructuredGrid, out_dir: str):
-    """Compute the average fiber helix and transverse angles."""
+def _compute_aha_fiber_angles(
+    mesh: pv.UnstructuredGrid, out_dir: str, num_layers: int = 9
+) -> tuple[pd.DataFrame, pd.DataFrame, pv.UnstructuredGrid]:
+    """Compute the average helix and transverse angles over N number of layers.
+
+    Parameters
+    ----------
+    mesh : pv.UnstructuredGrid
+        Mesh containing the fiber orientations, transmural, and rotational data.
+    out_dir : str
+        Output directory for saving results.
+    num_layers : int, default: 9
+        Number of layers to compute average angles over.
+
+    Returns
+    -------
+    tuple[pd.DataFrame, pd.DataFrame, pv.UnstructuredGrid]
+        DataFrames containing the average helix and transverse angles, and the updated mesh.
+
+    Raises
+    ------
+    ValueError
+        If the expected arrays are not found in the mesh data.
+    """
     expected_arrays = ["aha17", "fiber", "sheet", "transmural", "rotational"]
     for arr in expected_arrays:
         if arr not in mesh.cell_data and arr not in mesh.point_data:
             raise ValueError(f"Expected array '{arr}' not found in mesh data.")
 
+    if num_layers < 1 or not isinstance(num_layers, int):
+        raise ValueError("num_layers must be a positive integer.")
+
     aha_ids = mesh["aha17"]
     aha_elements = np.where(~np.isnan(aha_ids))[0]
     aha_model: pv.UnstructuredGrid = mesh.extract_cells(aha_elements)
     # aha_model.cell_data["fiber"] = aha_model.cell_data["fiber"] * -1  # flip fiber direction
     aha_ids = aha_model.cell_data["aha17"]
 
-    # flip fibers
-    # aha_model.cell_data["fiber"] = 1 * aha_model.cell_data["fiber"]
-
     # load fibers and sheets at cells
     el_fibers = aha_model.cell_data["fiber"]
     el_sheets = aha_model.cell_data["sheet"]  # noqa N841
@@ -189,7 +212,9 @@ def compute_aha_fiber_angles(mesh: pv.UnstructuredGrid, out_dir: str):
     return df_helix_angle, df_transverse_angle, aha_model
 
 
-def plot_fiber_aha_angles(data: pd.DataFrame | str):
+def plot_fiber_aha_angles(
+    data: pd.DataFrame | str, angle_type: Literal["helix", "transverse"] = "helix"
+):
     """
     Plot average fiber helical orientation in each AHA as a function of transmural depth
 
@@ -235,8 +260,15 @@ def plot_fiber_aha_angles(data: pd.DataFrame | str):
         axs[i, j].set_xlim(xmin=-1, xmax=1)
         axs[i, j].set_ylim(ymin=-100, ymax=100)
 
+    if angle_type == "transverse":
+        y_label = "$\\alpha_t$"
+    elif angle_type == "helix":
+        y_label = "$\\alpha_h$"
+    else:
+        y_label = "$\\alpha$"
+
     for ax in axs.flat:
-        ax.set(xlabel="Transmural Depth", ylabel="$\\alpha_h$")
+        ax.set(xlabel="Transmural Depth", ylabel=y_label)
 
     for ax in axs.flat:
         ax.label_outer()