PoinCloud deformation improvements

Author: gabe-rbo

CITATION.cff

@@ -5,7 +5,7 @@ authors:
     given-names: "Gabriel"
     email: "gabriel.ribeiro@dcc.ufmg.br"
 title: "pySurgery: Computational Surgery Theory in Python"
-version: 2.2.6
+version: 2.2.7
 date-released: 01-05-2026
 url: "https://github.com/gabe-rbo/pysurgery"
 repository-code: "https://github.com/gabe-rbo/pysurgery"

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "pysurgery"
-version = "2.2.6"
+version = "2.2.7"
 description = "A high-performance Python library for Computational Surgery Theory."
 dependencies = [
     "numpy>=2.4.4",
@@ -83,7 +83,7 @@ python_files = ["test_*.py"]
 "*" = ["*.jl", "*.yaml", "*.json"]
 
 [tool.bumpversion]
-current_version = "2.2.6"
+current_version = "2.2.7"
 commit = true
 tag = true
 

pysurgery/__init__.py

@@ -316,7 +316,7 @@
 
 from . import integrations
 
-__version__ = "2.2.6"
+__version__ = "2.2.7"
 
 def __getattr__(name):
     if name == "JuliaBridge":

pysurgery/geometry/point_cloud.py

@@ -2315,6 +2315,7 @@ class SimplicialComplex(ChainComplex):
     _simplices_table: Dict[int, List[Tuple[int, ...]]] = PrivateAttr(default_factory=dict)
     _point_cloud_to_simplices: Dict[int, List[Tuple[int, ...]]] = PrivateAttr(default_factory=dict)
     _simplices_to_point_cloud: Dict[Tuple[int, ...], np.ndarray] = PrivateAttr(default_factory=dict)
+    _point_cloud_cache: Optional[Any] = PrivateAttr(default=None)
     coefficient_ring: str = "Z"
     filtration: Dict[Tuple[int, ...], float] = Field(default_factory=dict)
 
@@ -3149,8 +3150,13 @@ def simplices_dict(self) -> Dict[int, List[Tuple[int, ...]]]:
     def point_cloud(self) -> Optional["PointCloud"]:
         """Return the PointCloud wrapper for coordinates, if coordinates exist."""
         if hasattr(self, "_coordinates") and self._coordinates is not None:
+            if getattr(self, "_point_cloud_cache", None) is not None:
+                if self._point_cloud_cache.points is self._coordinates:
+                    return self._point_cloud_cache
             from ..geometry.point_cloud import PointCloud
-            return PointCloud(self._coordinates, parent=self)
+            pc = PointCloud(self._coordinates, parent=self)
+            self._point_cloud_cache = pc
+            return pc
         return None
 
     @point_cloud.setter
@@ -3160,15 +3166,18 @@ def point_cloud(self, pc: Optional[Union["PointCloud", np.ndarray]]) -> None:
             self._coordinates = None
             self._point_cloud_to_simplices = {}
             self._simplices_to_point_cloud = {}
+            self._point_cloud_cache = None
         else:
             from ..geometry.point_cloud import PointCloud
             if isinstance(pc, PointCloud):
                 pts = pc.points
+                self._point_cloud_cache = pc
             else:
                 pts = np.asarray(pc, dtype=np.float64)
+                self._point_cloud_cache = PointCloud(pts, parent=self)
             self._coordinates = pts
             self._generate_point_cloud_mappings(pts)
-            self._link_point_cloud(pc)
+            self._link_point_cloud(self._point_cloud_cache)
 
     @property
     def point_cloud_to_simplices(self) -> Dict[int, List[Tuple[int, ...]]]:
@@ -3202,7 +3211,102 @@ def _link_point_cloud(self, points: Any) -> None:
         if isinstance(points, PointCloud):
             points._parent = self
 
+    def verify_transformation_collision(self, tol: float = 1e-8) -> List[Dict[str, Any]]:
+        """Verifies if any self-intersections (collisions) occurred during the transformation history.
+
+        This method walks through the sequential history of geometric deformations applied
+        to the linked PointCloud. At each step (including the initial undeformed state), it
+        recreates the coordinate realization, constructs a piecewise-linear map (PLMap),
+        and runs broad-phase/narrow-phase intersection detection to identify any overlaps
+        between non-adjacent simplices.
+
+        Algorithm & Mathematical Foundations:
+            1. Initializes a temporary, parent-less `PointCloud` using the `original_points`
+               coordinates of the linked point cloud.
+            2. For each state (Step 0 up to Step N):
+               - Constructs a PLMap f: |K| -> R^D using the active coordinates.
+               - Checks for self-intersections by running `detect_self_intersections(pl_map)`.
+                 An intersection is detected if the realizations of two non-adjacent simplices
+                 intersect in ambient space within the given numerical tolerance.
+               - If intersections are found, records details about the colliding simplices
+                 and the ambient coordinates of their vertices.
+               - Applies the next transformation in the history sequence to the coordinates.
+            3. Returns the log of all detected collisions.
+
+        Args:
+            tol (float): Numerical tolerance for intersection tests (distance below which
+                simplices are considered to overlap). Defaults to 1e-8.
+
+        Returns:
+            List[Dict[str, Any]]: A list of dictionaries representing steps where collisions
+            were detected. Each dictionary contains:
+                - "step" (int): The step index (0 for initial state, 1..N for transformations).
+                - "method" (str): The name of the transformation applied (or "initial_state").
+                - "args" (Dict[str, Any]): The arguments supplied to the transformation.
+                - "witnesses" (List[Dict[str, Any]]): Details of each collision witness, including:
+                    - "simplex_a" (Tuple[int, ...]): Vertices of the first simplex.
+                    - "simplex_b" (Tuple[int, ...]): Vertices of the second simplex.
+                    - "simplex_a_coordinates" (List[List[float]]): Ambient coordinates of simplex_a vertices.
+                    - "simplex_b_coordinates" (List[List[float]]): Ambient coordinates of simplex_b vertices.
+                    - "kind" (str): Intersection type (e.g. 'segment_segment').
+                    - "distance" (float): Minimum distance between the simplices.
+                    - "overlap_dimension" (int): Dimension of the intersection set.
+                    - "notes" (List[str]): Diagnostic notes.
+
+        Raises:
+            ValueError: If the simplicial complex does not have a linked PointCloud/coordinates.
+        """
+        pc = self.point_cloud
+        if pc is None:
+            raise ValueError("SimplicialComplex has no coordinates/PointCloud linked.")
+
+        history = pc._history
+        original_pts = pc._original_points
+
+        from ..geometry.point_cloud import PointCloud
+        from ..geometry.embedding import PLMap, detect_self_intersections
+
+        temp_pc = PointCloud(original_pts.copy(), original_points=original_pts)
+        collisions = []
+
+        def check_collisions(step_idx: int, method_name: str, args: Dict[str, Any]) -> None:
+            pl_map = PLMap.from_source(self, coordinates=temp_pc.points)
+            report = detect_self_intersections(pl_map, tol=tol)
+            if report.has_intersections:
+                witness_list = []
+                for w in report.witnesses:
+                    witness_list.append({
+                        "simplex_a": w.simplex_a,
+                        "simplex_b": w.simplex_b,
+                        "simplex_a_coordinates": temp_pc.points[list(w.simplex_a)].tolist(),
+                        "simplex_b_coordinates": temp_pc.points[list(w.simplex_b)].tolist(),
+                        "kind": w.kind,
+                        "distance": w.distance,
+                        "overlap_dimension": w.overlap_dimension,
+                        "notes": w.notes,
+                    })
+                collisions.append({
+                    "step": step_idx,
+                    "method": method_name,
+                    "args": args,
+                    "witnesses": witness_list,
+                })
+
+        # Check initial state (Step 0)
+        check_collisions(step_idx=0, method_name="initial_state", args={})
+
+        # Apply transformations sequentially and check each step
+        for i, entry in enumerate(history):
+            method_name = entry["method"]
+            args = entry["args"]
+            method = getattr(temp_pc, method_name)
+            method(**args)
+            check_collisions(step_idx=i + 1, method_name=method_name, args=args)
+
+        return collisions
+
     @classmethod
+
     def concatenate(cls, complexes: Iterable["SimplicialComplex"]) -> "SimplicialComplex":
         """Concatenate multiple simplicial complexes (simplex trees) into a single larger complex.
 

pysurgery/topology/complexes.py

@@ -667,3 +667,49 @@ def test_simplicial_complex_from_distance_matrix():
     if julia_engine.available:
         sc_jl = SimplicialComplex.from_distance_matrix(dist_mat, epsilon=1.1, max_dimension=2, backend="julia")
         assert sc_jl.count_simplices(2) == 1
+
+
+def test_simplicial_complex_verify_transformation_collision():
+    from pysurgery.geometry import SpaceBlock
+    from pysurgery.topology.complexes import SimplicialComplex
+
+    # 1. Create a simplicial complex representing two parallel, disjoint edges
+    pts = np.array([
+        [0.0, 0.0],  # vertex 0
+        [2.0, 0.0],  # vertex 1
+        [0.0, 1.0],  # vertex 2
+        [2.0, 1.0]   # vertex 3
+    ])
+    sc = SimplicialComplex.from_maximal_simplices([(0, 1), (2, 3)])
+    sc.point_cloud = pts
+
+    # Initial state (Step 0) has no collisions
+    collisions_init = sc.verify_transformation_collision()
+    assert len(collisions_init) == 0
+
+    # 2. Deform: shift vertex 2 down by Y=-2.0 so that the edges cross
+    # Define a SpaceBlock containing vertex 2 only
+    sb = SpaceBlock([-0.1, 0.9], [0.1, 1.1])
+    sc.point_cloud.translate([0.0, -2.0], movable_blocks=sb)
+
+    # Now verify collisions
+    collisions = sc.verify_transformation_collision()
+    
+    # Step 1 should contain a collision
+    assert len(collisions) == 1
+    report = collisions[0]
+    assert report["step"] == 1
+    assert report["method"] == "translate"
+    assert len(report["witnesses"]) == 1
+    
+    witness = report["witnesses"][0]
+    # The two edges (0, 1) and (2, 3) collided
+    assert set(witness["simplex_a"]) == {0, 1}
+    assert set(witness["simplex_b"]) == {2, 3}
+    assert witness["kind"] == "segment_segment"
+    
+    # Verify coordinate values of simplex A
+    np.testing.assert_allclose(witness["simplex_a_coordinates"], [[0.0, 0.0], [2.0, 0.0]])
+    # Vertex 2 coordinates were shifted to [0.0, -1.0], vertex 3 remained at [2.0, 1.0]
+    np.testing.assert_allclose(witness["simplex_b_coordinates"], [[0.0, -1.0], [2.0, 1.0]])
+

tests/test_complexes.py

@@ -632,17 +632,86 @@ def test_point_cloud_space_blocks():
     )
     pc_reverted = pc_revert_test.revert()
     np.testing.assert_allclose(pc_reverted.points, pts_def)
-
-
     pc_undone = pc_trans.undo()
     np.testing.assert_allclose(pc_undone.points, pts_def)
 
-
     # Test static_blocks
     pc_static = pc_def.translate([10.0, 10.0], static_blocks=top_right_block)
     np.testing.assert_allclose(pc_static.points, [[8.0, 8.0], [2.0, 2.0]])
     np.testing.assert_allclose(pc_static.revert().points, pts_def)
 
 
+def test_point_cloud_seam_alignment():
+    # 8 points on a circle of radius 1 centered at (0, 1)
+    # The center of curvature for curvature=1.0 at anchor=(0, 0) is (0, 1)
+    angles = np.linspace(0, 2 * np.pi, 8, endpoint=False)
+    pts = np.column_stack([np.sin(angles), 1.0 - np.cos(angles)])
+    pc = PointCloud(pts)
+
+    # Let's open at index 2 (angle = pi/2, coordinates approx [1.0, 1.0])
+    # The seam is at angle = pi/2
+    unbent = pc.unbend(curvature=1.0, axis=0, control_axis=1, anchor=(0, 0), open_at=2)
+    
+    # After unbending, the seam vertex (index 2) should land at the boundary (-pi or pi)
+    x_seam = unbent.points[2, 0]
+    # Check if x_seam is close to -pi or pi
+    assert np.isclose(np.abs(x_seam), np.pi)
+
+    # Revert should return to the original circle coordinates
+    reverted = unbent.revert()
+    np.testing.assert_allclose(reverted.points, pts, atol=1e-8)
+
+    # Undo should also return to original coordinates
+    unbent2 = pc.unbend(curvature=1.0, axis=0, control_axis=1, anchor=(0, 0), open_at=2)
+    undone = unbent2.undo()
+    np.testing.assert_allclose(undone.points, pts, atol=1e-8)
+
+    # Test open_at with a coordinate point
+    seam_pt = np.array([1.0, 1.0])
+    unbent_pt = pc.unbend(curvature=1.0, axis=0, control_axis=1, anchor=(0, 0), open_at=seam_pt)
+    assert np.isclose(np.abs(unbent_pt.points[2, 0]), np.pi)
+    np.testing.assert_allclose(unbent_pt.revert().points, pts, atol=1e-8)
 
 
+def test_point_cloud_frames():
+    pts = np.array([[1.0, 2.0]])
+    pc = PointCloud(pts)
+
+    # Translate frames
+    frames_list = list(pc.frames("translate", steps=4, translation_vector=[4.0, 8.0]))
+    assert len(frames_list) == 5
+    
+    # Check intermediate steps
+    expected_t = [0.0, 0.25, 0.5, 0.75, 1.0]
+    for idx, (t, frame) in enumerate(frames_list):
+        assert np.isclose(t, expected_t[idx])
+        expected_pts = pts + t * np.array([4.0, 8.0])
+        np.testing.assert_allclose(frame.points, expected_pts)
+
+
+def test_point_cloud_collision_and_clearance():
+    pc1 = PointCloud(np.array([[0.0, 0.0]]))
+    pc2 = PointCloud(np.array([[3.0, 4.0]]))
+
+    # KD-Tree distance check
+    dist = pc1.min_distance_to(pc2)
+    assert np.isclose(dist, 5.0)
+
+    # clearance check with sequence of actions
+    actions = [
+        ("translate", {"translation_vector": [4.0, 0.0]}),
+        ("translate", {"translation_vector": [0.0, 3.0]})
+    ]
+    obstacle = PointCloud(np.array([[2.5, 0.0]]))
+
+    # Under safety margin 1.0, translating from 0.0 to 4.0 along X will cross
+    # the obstacle at X=2.5. Distance to obstacle is |X - 2.5|.
+    # At t=0.0 (X=0.0): dist = 2.5
+    # At t=0.5 (X=2.0): dist = 0.5 < 1.0 -> VIOLATION!
+    # At t=0.75 (X=3.0): dist = 0.5 < 1.0 -> VIOLATION!
+    # At t=1.0 (X=4.0): dist = 1.5
+    violations = pc1.verify_isotopy_clearance(obstacle, actions, steps_per_action=4, safety_margin=1.0)
+    assert len(violations) > 0
+    # The first violation should be during the first translate
+    assert violations[0][1] == "translate"
+    assert violations[0][2] < 1.0