Document edge case behavior and add tests for smoothing small polygons

lapy/__init__.py

@@ -1,5 +1,5 @@
 from ._version import __version__  # noqa: F401
-from .polygon import Polygon
+from .polygon import Polygon  # noqa: F401
 from .solver import Solver  # noqa: F401
 from .tet_mesh import TetMesh  # noqa: F401
 from .tria_mesh import TriaMesh  # noqa: F401

lapy/polygon.py

@@ -70,14 +70,27 @@ def __init__(self, points: np.ndarray, closed: bool = False):
         if self.points.size == 0:
             raise ValueError("Polygon has no points (empty)")
 
-        # Transpose if necessary
-        if self.points.shape[0] < self.points.shape[1]:
+        # Ensure points array is 2-dimensional
+        if self.points.ndim != 2:
+            raise ValueError("Points array must be 2-dimensional")
+
+        n_rows, n_cols = self.points.shape
+
+        # Support both (n_points, dim) and (dim, n_points) where dim is 2 or 3.
+        # Only transpose when it is unambiguous that the first dimension is dim.
+        if n_cols not in (2, 3) and n_rows in (2, 3):
+            logger.warning(
+                "Transposing points array from shape %s to %s; expected shape (n_points, dim).",
+                self.points.shape,
+                self.points.T.shape,
+            )
             self.points = self.points.T
+            n_rows, n_cols = self.points.shape
 
         # Support both 2D and 3D points
-        if self.points.shape[1] == 2:
+        if n_cols == 2:
             self._is_2d = True
-        elif self.points.shape[1] == 3:
+        elif n_cols == 3:
             self._is_2d = False
         else:
             raise ValueError("Points should have 2 or 3 coordinates")
@@ -177,14 +190,14 @@ def centroid(self) -> np.ndarray:
         y_closed = np.append(y, y[0])
         # Shoelace formula components
         cross = x_closed[:-1] * y_closed[1:] - x_closed[1:] * y_closed[:-1]
-        area = 0.5 * np.abs(cross.sum())
+        signed_area = 0.5 * cross.sum()
 
-        if area < sys.float_info.epsilon:
-            # Degenerate case: zero area
+        if abs(signed_area) < sys.float_info.epsilon:
+            # Degenerate case: zero or near-zero area
             return np.mean(self.points, axis=0)
 
-        cx = np.sum((x_closed[:-1] + x_closed[1:]) * cross) / (6.0 * area)
-        cy = np.sum((y_closed[:-1] + y_closed[1:]) * cross) / (6.0 * area)
+        cx = np.sum((x_closed[:-1] + x_closed[1:]) * cross) / (6.0 * signed_area)
+        cy = np.sum((y_closed[:-1] + y_closed[1:]) * cross) / (6.0 * signed_area)
         return np.array([cx, cy])
 
     def area(self) -> float:
@@ -241,6 +254,10 @@ def resample(
         Polygon
             Resampled polygon. Returns self if inplace=True, new instance otherwise.
         """
+        if n_points < 2:
+            raise ValueError("n_points must be at least 2")
+        if n_iter < 1:
+            raise ValueError("n_iter must be at least 1")
         def _resample_once(p: np.ndarray, n: int, is_closed: bool) -> np.ndarray:
             """Single resampling pass."""
             if is_closed:
@@ -277,6 +294,17 @@ def _construct_smoothing_matrix(self) -> sparse.csc_matrix:
         its neighbors (previous and next point). For open polygons, boundary
         points (first and last) are kept fixed.
 
+        The method handles polygons of any size:
+
+        - For open polygons with 2 points: Both boundary points remain fixed
+          (identity matrix), so smoothing has no effect.
+        - For open polygons with 3+ points: Boundary points are fixed, interior
+          points are averaged with their neighbors.
+        - For closed polygons with 2 points: Each point is averaged with its
+          neighbor, causing them to converge to their midpoint.
+        - For closed polygons with 3+ points: All points are averaged with
+          their neighbors in a circular manner.
+
         Returns
         -------
         scipy.sparse.csc_matrix
@@ -344,6 +372,11 @@ def smooth_laplace(
         Polygon
             Smoothed polygon. Returns self if inplace=True, new instance otherwise.
         """
+        # Input validation to enforce documented parameter ranges
+        if not isinstance(n, int) or n <= 0:
+            raise ValueError(f"n must be a positive integer, got {n!r}")
+        if not (0.0 <= lambda_ <= 1.0):
+            raise ValueError(f"lambda_ must be in the range [0, 1], got {lambda_!r}")
         mat = self._construct_smoothing_matrix()
         points_smooth = self.points.copy()
 
@@ -386,6 +419,12 @@ def smooth_taubin(
         Polygon
             Smoothed polygon. Returns self if inplace=True, new instance otherwise.
         """
+        if n <= 0:
+            raise ValueError("n must be a positive integer")
+        if lambda_ <= 0:
+            raise ValueError("lambda_ must be positive")
+        if mu >= 0:
+            raise ValueError("mu must be negative")
         mat = self._construct_smoothing_matrix()
         points_smooth = self.points.copy()
 

lapy/utils/tests/test_polygon.py

@@ -35,7 +35,7 @@ def test_init_empty_raises(self):
     def test_init_invalid_dimensions_raises(self):
         """Test that invalid dimensions raise ValueError."""
         with pytest.raises(ValueError, match="2 or 3 coordinates"):
-            Polygon(np.array([[0.0], [1.0]]))
+            Polygon(np.array([[0.0, 1.0, 2.0, 3.0]]))
 
     def test_length_open(self):
         """Test length computation for open polygon."""
@@ -173,9 +173,9 @@ def test_smooth_laplace_open(self):
             "Should preserve number of points"
         # First and last points should remain unchanged for open polygon
         assert np.allclose(smoothed.get_points()[0], points[0]), \
-            "First point should not change much"
+            "First point should remain unchanged (up to numerical precision)"
         assert np.allclose(smoothed.get_points()[-1], points[-1]), \
-            "Last point should not change much"
+            "Last point should remain unchanged (up to numerical precision)"
 
     def test_smooth_laplace_closed(self):
         """Test Laplace smoothing on closed polygon."""
@@ -221,6 +221,42 @@ def test_smooth_taubin_inplace(self):
 
         assert result is poly, "Should return self when inplace=True"
 
+    def test_smooth_two_point_open(self):
+        """Test smoothing on two-point open polygon."""
+        points = np.array([[0.0, 0.0], [1.0, 1.0]])
+        poly = Polygon(points, closed=False)
+        smoothed = poly.smooth_laplace(n=5, lambda_=0.5, inplace=False)
+
+        # Both points should remain unchanged (boundary points are fixed)
+        assert np.allclose(smoothed.get_points(), points), \
+            "Two-point open polygon should not change when smoothed"
+
+    def test_smooth_two_point_closed(self):
+        """Test smoothing on two-point closed polygon."""
+        points = np.array([[0.0, 0.0], [2.0, 2.0]])
+        poly = Polygon(points, closed=True)
+        smoothed = poly.smooth_laplace(n=5, lambda_=0.5, inplace=False)
+
+        # Points should converge to their midpoint
+        midpoint = np.mean(points, axis=0)
+        assert np.allclose(smoothed.get_points(), midpoint, atol=1e-10), \
+            "Two-point closed polygon should converge to midpoint"
+
+    def test_smooth_three_point_open(self):
+        """Test smoothing on three-point open polygon."""
+        points = np.array([[0.0, 0.0], [0.5, 2.0], [1.0, 0.0]])
+        poly = Polygon(points, closed=False)
+        smoothed = poly.smooth_laplace(n=5, lambda_=0.5, inplace=False)
+
+        # First and last points should remain fixed
+        assert np.allclose(smoothed.get_points()[0], points[0]), \
+            "First point should remain fixed in open polygon"
+        assert np.allclose(smoothed.get_points()[-1], points[-1]), \
+            "Last point should remain fixed in open polygon"
+        # Middle point should be smoothed (moved toward average of neighbors)
+        assert smoothed.get_points()[1, 1] < points[1, 1], \
+            "Middle point should be smoothed downward"
+
 
 if __name__ == "__main__":
     pytest.main([__file__, "-v"])