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 diff --git a/src/ansys/health/heart/utils/compute_fiber_angles.py b/src/ansys/health/heart/utils/compute_fiber_angles.py new file mode 100644 index 000000000..d5eae28d9 --- /dev/null +++ b/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() diff --git a/src/ansys/health/heart/utils/misc.py b/src/ansys/health/heart/utils/misc.py index bba76f692..42ae7e3a0 100644 --- a/src/ansys/health/heart/utils/misc.py +++ b/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 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..cfcb18412 --- /dev/null +++ b/tests/heart/utils/test_compute_fiber_angles.py @@ -0,0 +1,99 @@ +# 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.misc import angle_between_vectors + + +def test_compute_signed_angles(): + """Test computing signed angles between vectors.""" + + # Local element coordinate system + # 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], + [0.0, -1.0, 1.0], + [0.0, -1.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] + -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] + ] + ) + + # 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 = 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)