feat: compute fiber helical angle

doc/source/changelog/827.added.md

@@ -0,0 +1 @@
+Compute fiber helical angle

src/ansys/health/heart/utils/compute_fiber_angles.py

@@ -0,0 +1,279 @@
+# 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.
+
+"""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
+from numpy.linalg import norm
+import pandas as pd
+import pyvista as pv
+
+from ansys.health.heart.utils.misc import angle_between_vectors
+
+
+def _compute_gradient(mesh: pv.UnstructuredGrid, scalar: str, interpolate_to_cell: bool = False):
+    """Compute gradient of a scalar field on the mesh.
+
+    Parameters
+    ----------
+    mesh : pv.UnstructuredGrid
+        Input mesh.
+    scalar : str
+        Name of the scalar field.
+
+    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_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"]
+
+    # load fibers and sheets at cells
+    el_fibers = aha_model.cell_data["fiber"]
+    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, "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)
+    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]
+
+    aha_model.cell_data["grad_rotational"] = el_grad_r
+
+    # 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
+
+    # compute angle between rotational and fiber vectors in plane defined
+    # by longitudinal vector (transverse angle)
+
+    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_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.cell_data["depth"] > lower
+        mask &= aha_model.cell_data["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_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)
+
+    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, angle_type: Literal["helix", "transverse"] = "helix"
+):
+    """
+    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)
+
+    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=y_label)
+
+    for ax in axs.flat:
+        ax.label_outer()
+
+    axs[2, 5].remove()
+
+    # plt.savefig("fiber_helical_angles.png", bbox_inches="tight")
+    plt.show()

src/ansys/health/heart/utils/misc.py

@@ -23,8 +23,10 @@
 """Module containing miscellaneous methods."""
 
 import os
+from typing import Literal
 
 import numpy as np
+from numpy.linalg import norm
 from scipy.spatial import cKDTree
 
 from ansys.health.heart import LOG as LOGGER
@@ -279,3 +281,59 @@ 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 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 between N number of M-dimensional vectors.
+
+    Parameters
+    ----------
+    x : np.ndarray
+        (N,M) array with first (reference) vectors.
+    y : np.ndarray
+        (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"
+        Representation of the angle:
+            - "2-quadrant-signed": angles in the range [-180, 180)
+            - "4-quadrant": angles in the range [0, 360)
+
+    Returns
+    -------
+    np.ndarray
+        Array of signed angles in degrees of size (N,).
+
+    Notes
+    -----
+    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
+    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)
+
+    # 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 = np.degrees(angles)
+
+    if representation == "4-quadrant":
+        return np.mod(angles + 360, 360)
+    else:
+        return angles