diff --git a/doc/source/changelog/1147.added.md b/doc/source/changelog/1147.added.md new file mode 100644 index 000000000..280802fdc --- /dev/null +++ b/doc/source/changelog/1147.added.md @@ -0,0 +1 @@ +Update conductionpath diff --git a/src/ansys/health/heart/models_utils.py b/src/ansys/health/heart/models_utils.py index 7799247cd..ee1c97bee 100644 --- a/src/ansys/health/heart/models_utils.py +++ b/src/ansys/health/heart/models_utils.py @@ -31,7 +31,7 @@ from ansys.health.heart import LOG as LOGGER from ansys.health.heart.landmarks import LandMarks import ansys.health.heart.models as models -from ansys.health.heart.objects import CapType, Point +from ansys.health.heart.objects import CapType, Point, SurfaceMesh from ansys.health.heart.pre.conduction_path import ConductionPath, ConductionPathType @@ -372,8 +372,11 @@ def define_full_conduction_system( left_bundle.up_path = his_left left_bundle.down_path = left_purkinje - surface_ids = [model.right_ventricle.endocardium.id, model.right_ventricle.septum.id] - endo_surface = model.mesh.get_surface(surface_ids) + # create complete endocardial surface for the right ventricle + endo_surface = SurfaceMesh( + model.right_ventricle.endocardium.merge(model.right_ventricle.septum), + name="right_ventricle_endo_septum", + ) right_bundle = ConductionPath.create_from_keypoints( name=ConductionPathType.RIGHT_BUNDLE_BRANCH, diff --git a/src/ansys/health/heart/pre/conduction_path.py b/src/ansys/health/heart/pre/conduction_path.py index 467d86b06..7d6c11197 100644 --- a/src/ansys/health/heart/pre/conduction_path.py +++ b/src/ansys/health/heart/pre/conduction_path.py @@ -83,7 +83,7 @@ def __init__( mesh: Mesh, id: int, is_connected: np.ndarray, - relying_surface: pv.PolyData, + relying_surface: pv.PolyData | SurfaceMesh, material: EPMaterial = EPMaterial.DummyMaterial(), up_path: ConductionPath | None = None, down_path: ConductionPath | None = None, @@ -101,7 +101,7 @@ def __init__( ID of the conduction path. is_connected : np.ndarray Mask array of points connected to the solid mesh. - relying_surface : pv.PolyData + relying_surface : pv.PolyData|SurfaceMesh Surface mesh that the conduction path relies on. material : EPMaterial, default: EPMaterial.DummyMaterial() EP Material property. @@ -311,6 +311,73 @@ def add_pmj_path( return self + def _compute_relative_coord(self) -> list[tuple[int, np.ndarray]]: + """ + Compute barycentric coordinates of the mesh points relative to the relying surface. + + Returns + ------- + list[tuple[int, np.ndarray]] + List of tuples containing cell ID and barycentric coordinates + for each point in the mesh. + """ + bary_info = [] + for pt in self.mesh.points: + cell_id = self.relying_surface.find_containing_cell(pt) + if cell_id < 0: + LOGGER.warning( + f"Point {pt} is not contained in the relying surface. Use closest cell search." + ) + cell_id = self.relying_surface.find_closest_cell(pt) + + cell = self.relying_surface.get_cell(cell_id) + verts = self.relying_surface.points[cell.point_ids] + # Compute barycentric coordinates + v0, v1, v2 = verts + a = pt - v0 + # Solve for barycentric coordinates (u, v) + u, v = np.linalg.lstsq(np.column_stack((v1 - v0, v2 - v0)), a, rcond=None)[0] + w = 1 - u - v + bary_info.append((cell_id, np.array([w, u, v]))) + return bary_info + + def deform_to_surface(self, surface: pv.PolyData) -> ConductionPath | None: + """Deform the conduction path to a new surface. + + Parameters + ---------- + surface : pv.PolyData + New surface to deform the conduction path to. + + Returns + ------- + ConductionPath | None + Updated conduction path if successful, None otherwise. + """ + if not np.array_equal(surface.faces, self.relying_surface.faces): + LOGGER.error(f"New surface does not match the relying surface of {self.name}.") + return + + try: + bary_info = self._compute_relative_coord() + except ValueError as e: + LOGGER.error( + f"Failed to compute barycentric coordinates for {self.name} on the surface: {e}" + ) + return + + # Map each mesh point to the new surface using barycentric coordinates + new_points = self.mesh.points.copy() + for i, (cell_id, bary) in enumerate(bary_info): + verts = surface.get_cell(cell_id).points + new_pt = bary[0] * verts[0] + bary[1] * verts[1] + bary[2] * verts[2] + new_points[i] = new_pt + + self.mesh.points = new_points + self.relying_surface = surface + + return self + @staticmethod def create_from_keypoints( name: ConductionPathType, @@ -737,21 +804,26 @@ def _create_path_in_solid( tetras[:, [0, 2, 3]], tetras[:, [1, 2, 3]], ) - ) # TODO: replace by pv extract_surface() + ) + segment = [] for i, j in zip(ids[0:-1], ids[1:]): for tri in triangles: if i in tri and j in tri: segment.append(tri) - break + break # keeping only one triangle per line is enough segment = np.array(segment) - surf = SurfaceMesh( - name="his_bundle_segment", # NOTE - triangles=segment, - nodes=sub_mesh.points, + surf = pv.PolyData( + sub_mesh.points, + np.hstack([np.insert(segment[i], 0, 3) for i in range(segment.shape[0])]), ) - return path, surf + + # keep global point IDs as a necessary property for the surface mesh + surf.point_data["_global-point-ids"] = sub_mesh.point_data["_global-point-ids"] + surf = surf.clean() # remove unused nodes + + return path, SurfaceMesh(surf) def _mesh_to_nx_graph(mesh: pv.UnstructuredGrid) -> nx.Graph: diff --git a/src/ansys/health/heart/simulator.py b/src/ansys/health/heart/simulator.py index abb21640a..78cc1a05e 100644 --- a/src/ansys/health/heart/simulator.py +++ b/src/ansys/health/heart/simulator.py @@ -808,6 +808,41 @@ def _write_main_simulation_files( return + def update_conduction_paths(self): + """Update conduction paths in the model after stress-free computation.""" + multi_surface = pv.MultiBlock() + multi_path = pv.MultiBlock() + new_paths = [] + for path in self.model.conduction_paths: + if isinstance(path.relying_surface, SurfaceMesh): + ids = path.relying_surface.global_node_ids_triangles + elif isinstance(path.relying_surface, pv.PolyData): + try: + ids = path.relying_surface["_global-point-ids"] + except KeyError: + msg = f"Conduction path {path.name} relying surface does not have" + "'_global-point-ids' array." + + LOGGER.error(msg) + raise KeyError(msg) + + # Update relying surface coordinates + new_coords = self.model.mesh.points[ids] + new_surface = path.relying_surface.copy() + new_surface.points = new_coords + + # deform path to the new surface + path2 = path.deform_to_surface(new_surface) + new_paths.append(path2) + + multi_surface.append(path2.relying_surface) + multi_path.append(path2.mesh) + + multi_surface.save(os.path.join(self.root_directory, "update_conduction_surfaces.vtm")) + multi_path.save(os.path.join(self.root_directory, "update_conduction_paths.vtm")) + + self.model.assign_conduction_paths(new_paths) + def _kill_all_ansyscl(): """Kill all Ansys license clients.""" diff --git a/tests/heart/pre/test_conduction_path.py b/tests/heart/pre/test_conduction_path.py index 65d19e9bb..63dca1d80 100644 --- a/tests/heart/pre/test_conduction_path.py +++ b/tests/heart/pre/test_conduction_path.py @@ -238,7 +238,7 @@ def simple_conduction_path(): surface = pv.PolyData(surface_points, faces=surface_faces) # Create a simple line mesh (2 points, 1 line) - line_points = np.array([[1, 1, 0], [0.3, 0.3, 0]]) + line_points = np.array([[0, 1, 0], [0.3, 0.3, 0]]) line_lines = np.array([2, 0, 1]) mesh = pv.PolyData(line_points, lines=line_lines) is_connected = np.array([0, 0]) @@ -282,8 +282,21 @@ def test_add_pmj_path_node(simple_conduction_path): assert np.array_equal(cp.is_connected, np.array([0, 1, 0])) +def test_update_to_surface(simple_conduction_path): + """Test update to surface.""" + cp: ConductionPath = simple_conduction_path + + surface_points = np.array([[0, 0, 0], [2, 0, 0], [0, 2, 0]]) + surface_faces = np.hstack([[3, 0, 1, 2]]) + surface = pv.PolyData(surface_points, faces=surface_faces) + + cp.deform_to_surface(surface) + assert np.allclose(cp.mesh.points[0], np.array([0, 2, 0])) + assert np.allclose(cp.mesh.points[1], np.array([0.6, 0.6, 0])) + + def test_create_path_on_surface_center(): - """Test conductionbeams can be initialized correctly on a pyvista sphere with 3 keypoints.""" + """Test that conduction beams can be initialized correctly on a PyVista sphere with three keypoints.""" # Create a sphere mesh sphere = pv.Sphere(radius=1.0, theta_resolution=30, phi_resolution=30) # Define 3 keypoints on the sphere surface diff --git a/tests/heart/simulator/test_simulator.py b/tests/heart/simulator/test_simulator.py index d5ffa4833..67300b8b4 100644 --- a/tests/heart/simulator/test_simulator.py +++ b/tests/heart/simulator/test_simulator.py @@ -32,6 +32,7 @@ from ansys.health.heart.exceptions import LSDYNATerminationError import ansys.health.heart.models as models +from ansys.health.heart.objects import SurfaceMesh from ansys.health.heart.settings.settings import DynaSettings # import after mocking. @@ -377,3 +378,50 @@ def test_find_dynain_file(dynain_files, expected, expected_error, mechanics_simu assert mechanics_simulator._find_dynain_file(zerop_folder) == expected mock_glob.assert_called_once() + + +@mock.patch("ansys.health.heart.simulator.pv.MultiBlock") +@pytest.mark.parametrize("surface_type", ["SurfaceMesh", "PolyData"]) +def test_update_conduction_paths(mock_multiblock, surface_type): + # Mocks + mock_model = mock.MagicMock() + mock_mesh = mock.MagicMock() + mock_model.mesh = mock_mesh + mock_model.assign_conduction_paths = mock.Mock() + + # Prepare conduction path mock + mock_path = mock.MagicMock() + mock_path.name = "test_path" + mock_path.deform_to_surface = mock.Mock(return_value=mock_path) + mock_path.mesh = pv.PolyData() + + if surface_type == "SurfaceMesh": + # Mock SurfaceMesh type and attributes + + mock_surface = mock.MagicMock(spec=SurfaceMesh) + mock_surface.__class__ = SurfaceMesh + mock_surface.global_node_ids_triangles = [0, 1] + mock_path.relying_surface = mock_surface + else: + # PolyData type and attributes + mock_surface = mock.MagicMock(spec=pv.PolyData) + mock_surface.__class__ = pv.PolyData + mock_surface.__getitem__.side_effect = ( + lambda key: [0, 1] if key == "_global-point-ids" else None + ) + mock_path.relying_surface = mock_surface + + mock_mesh.points = np.array([[1, 2, 3], [4, 5, 6]]) + mock_surface.copy.return_value = mock_surface + + mock_model.conduction_paths = [mock_path] + + with mock.patch.object(shutil, "which", return_value=1): + simulator = simulators.EPMechanicsSimulator(mock_model, dyna_settings=None) + simulator.model = mock_model + simulator.update_conduction_paths() + + # Check deform_to_surface called with updated surface + mock_path.deform_to_surface.assert_called_once_with(mock_surface) + # Check assign_conduction_paths called with new_paths + mock_model.assign_conduction_paths.assert_called_once_with([mock_path])