feat: update conductionpath

doc/source/changelog/1147.added.md

@@ -0,0 +1 @@
+Update conductionpath

src/ansys/health/heart/models_utils.py

@@ -31,7 +31,7 @@
 
 from ansys.health.heart import LOG as LOGGER
 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,
@@ -384,6 +387,7 @@ def define_full_conduction_system(
         )
         right_bundle.up_path = his_right
         right_bundle.down_path = right_purkinje
+
         return [
             left_purkinje,
             right_purkinje,

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
+            The 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  # keep 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:

src/ansys/health/heart/simulator.py

@@ -805,6 +805,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."""

tests/heart/pre/test_conduction_path.py

@@ -268,7 +268,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])
@@ -312,6 +312,19 @@ 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."""
     # Create a sphere mesh

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])