From bcd063a03e7fe3b2bb9a692266a78f57fb4971eb Mon Sep 17 00:00:00 2001
From: Martin Reuter <reuter@mit.edu>
Date: Fri, 19 Dec 2025 02:14:17 +0100
Subject: [PATCH 01/14] Create Polygon class with a few functions and tests

---
 lapy/__init__.py                 |   1 +
 lapy/polygon.py                  | 469 ++++++++++++++++++++++++-------
 lapy/tria_mesh.py                |   7 +-
 lapy/utils/tests/test_polygon.py | 249 ++++++++++++----
 4 files changed, 579 insertions(+), 147 deletions(-)

diff --git a/lapy/__init__.py b/lapy/__init__.py
index 456d8916..30fb8830 100644
--- a/lapy/__init__.py
+++ b/lapy/__init__.py
@@ -1,4 +1,5 @@
 from ._version import __version__  # noqa: F401
+from .polygon import Polygon
 from .solver import Solver  # noqa: F401
 from .tet_mesh import TetMesh  # noqa: F401
 from .tria_mesh import TriaMesh  # noqa: F401
diff --git a/lapy/polygon.py b/lapy/polygon.py
index 291c9296..273e9747 100644
--- a/lapy/polygon.py
+++ b/lapy/polygon.py
@@ -1,115 +1,394 @@
-"""Functions for open and closed polygon paths.
+"""Polygon class for open and closed polygon paths.
 
-This module provides utilities for resampling 2D and 3D polygon paths with
-equidistant spacing. These functions are useful for processing level set
-paths extracted from triangle meshes or any other polyline data.
-
-Functions
----------
-resample
-    Resample a 2D or 3D polygon path to have a specified number of equidistant points.
+This module provides a Polygon class for processing 2D and 3D polygon paths with
+various geometric operations including resampling, smoothing, and metric computations.
 """
 
+import logging
+import sys
+
 import numpy as np
+from scipy import sparse
+
+logger = logging.getLogger(__name__)
 
 
-def resample(path: np.ndarray, n_points: int = 100, n_iter: int = 1, closed: bool = False) -> np.ndarray:
-    """Resample a 2D or 3D polygon path to have equidistant points.
+class Polygon:
+    """Class representing a polygon path (open or closed).
 
-    This function resamples a polygon path (open or closed) by creating
-    n_points that are approximately equidistantly spaced along the cumulative
-    Euclidean distance of the path. The resampling is performed using linear
-    interpolation independently for each coordinate. Optionally, the resampling
-    can be performed iteratively to achieve better numerical stability and more
-    accurate equidistant spacing.
+    This class handles 2D and 3D polygon paths with operations for resampling,
+    smoothing, and computing geometric properties like length, centroid, and area.
 
     Parameters
     ----------
-    path : np.ndarray
-        Array of shape (n, d) containing coordinates of the polygon vertices
+    points : np.ndarray
+        Array of shape (n, d) containing coordinates of polygon vertices
         in order, where d is 2 or 3 for 2D (x, y) or 3D (x, y, z) paths.
         For closed polygons, the last point should not duplicate the first point.
-    n_points : int, default=100
-        Number of points in the resampled path. Must be at least 2.
-    n_iter : int, default=1
-        Number of resampling iterations to perform. The default value of 1
-        performs a single resampling pass. Higher values (e.g., 3-5) provide
-        better equidistant spacing but increase computation time. Must be at
-        least 1. Iterative resampling is particularly useful for paths with
-        highly variable point density or strong curvature.
     closed : bool, default=False
-        If True, treats the path as a closed polygon and includes the segment
-        from the last point back to the first point in the resampling. If False,
-        treats the path as an open polyline. For closed paths, the output will
-        not duplicate the first point at the end.
-
-    Returns
-    -------
-    np.ndarray
-        Array of shape (n_points, d) containing the resampled coordinates
-        with approximately equidistant spacing along the path, where d
-        matches the input dimensionality.
-
-    Notes
-    -----
-    The function computes cumulative Euclidean distances between successive
-    points and uses linear interpolation to place new points at equally spaced
-    distance values.
-
-    When n_iter > 1, the resampling is applied iteratively. Each iteration
-    refines the point distribution by resampling the result from the previous
-    iteration. The iterative process converges to an approximately equidistant
-    point distribution, with each iteration making the spacing more uniform.
-
-    For closed paths (closed=True), the total path length includes the distance
-    from the last point back to the first point, and resampled points are
-    distributed along this closed loop. The returned points form a closed path
-    without duplicating the first point at the end.
+        If True, treats the path as a closed polygon. If False, treats it as
+        an open polyline.
+
+    Attributes
+    ----------
+    points : np.ndarray
+        Polygon vertex coordinates, shape (n_points, d).
+    closed : bool
+        Whether the polygon is closed or open.
+    _is_2d : bool
+        Internal flag indicating if polygon is 2D (True) or 3D (False).
+
+    Raises
+    ------
+    ValueError
+        If points array is empty.
+        If points don't have 2 or 3 coordinates.
 
     Examples
     --------
     >>> import numpy as np
-    >>> from lapy import polygon
-    >>> # Create a simple 3D open path and resample once
-    >>> path_3d = np.array([[0, 0, 0], [1, 0, 0], [1, 1, 0], [1, 1, 1]])
-    >>> resampled_3d = polygon.resample(path_3d, n_points=10)
-    >>> resampled_3d.shape
-    (10, 3)
-    >>> # Create a 2D closed path (e.g., a square without duplicating first point)
+    >>> from lapy.polygon import Polygon
+    >>> # Create a 2D closed polygon (square)
     >>> square = np.array([[0, 0], [1, 0], [1, 1], [0, 1]])
-    >>> resampled_square = polygon.resample(square, n_points=40, closed=True)
-    >>> resampled_square.shape
-    (40, 2)
-    >>> # Iterative resampling with closed path
-    >>> path_2d = np.array([[0, 0], [0.1, 0], [1, 0], [1, 1]])
-    >>> resampled_2d = polygon.resample(path_2d, n_points=20, n_iter=5, closed=True)
-    >>> resampled_2d.shape
-    (20, 2)
+    >>> poly = Polygon(square, closed=True)
+    >>> poly.is_2d()
+    True
+    >>> poly.length()
+    4.0
+    >>> # Create a 3D open path
+    >>> path_3d = np.array([[0, 0, 0], [1, 0, 0], [1, 1, 0], [1, 1, 1]])
+    >>> poly3d = Polygon(path_3d, closed=False)
+    >>> poly3d.is_2d()
+    False
     """
-    def _resample_once(p: np.ndarray, n: int, is_closed: bool) -> np.ndarray:
-        """Single resampling pass."""
-        if is_closed:
-            # For closed paths, append the first point to close the loop
-            p_closed = np.vstack([p, p[0]])
-            # Cumulative Euclidean distance between successive polygon points
-            d = np.cumsum(np.r_[0, np.sqrt((np.diff(p_closed, axis=0) ** 2).sum(axis=1))])
-            # Get linearly spaced points along the cumulative Euclidean distance
-            # Exclude the endpoint (d.max()) to avoid duplicating the first point
-            d_sampled = np.linspace(0, d.max(), n + 1)[:-1]
+
+    def __init__(self, points: np.ndarray, closed: bool = False):
+        self.points = np.array(points)
+        self.closed = closed
+
+        # Validate non-empty polygon
+        if self.points.size == 0:
+            raise ValueError("Polygon has no points (empty)")
+
+        # Transpose if necessary
+        if self.points.shape[0] < self.points.shape[1]:
+            self.points = self.points.T
+
+        # Support both 2D and 3D points
+        if self.points.shape[1] == 2:
+            self._is_2d = True
+        elif self.points.shape[1] == 3:
+            self._is_2d = False
+        else:
+            raise ValueError("Points should have 2 or 3 coordinates")
+
+    def is_2d(self) -> bool:
+        """Check if the polygon is 2D.
+
+        Returns
+        -------
+        bool
+            True if polygon is 2D, False if 3D.
+        """
+        return self._is_2d
+
+    def is_closed(self) -> bool:
+        """Check if the polygon is closed.
+
+        Returns
+        -------
+        bool
+            True if polygon is closed, False if open.
+        """
+        return self.closed
+
+    def n_points(self) -> int:
+        """Get number of points in polygon.
+
+        Returns
+        -------
+        int
+            Number of points.
+        """
+        return self.points.shape[0]
+
+    def get_points(self) -> np.ndarray:
+        """Get polygon points.
+
+        Returns
+        -------
+        np.ndarray
+            Point array of shape (n, 2) or (n, 3).
+        """
+        return self.points
+
+    def length(self) -> float:
+        """Compute total length of polygon path.
+
+        For closed polygons, includes the segment from last to first point.
+
+        Returns
+        -------
+        float
+            Total path length.
+        """
+        if self.closed:
+            points_closed = np.vstack([self.points, self.points[0]])
+            edge_vecs = np.diff(points_closed, axis=0)
+        else:
+            edge_vecs = np.diff(self.points, axis=0)
+
+        edge_lens = np.sqrt((edge_vecs**2).sum(axis=1))
+        return edge_lens.sum()
+
+    def centroid(self) -> np.ndarray:
+        """Compute centroid of polygon.
+
+        For open polygons or closed 3D polygons, returns the simple arithmetic mean
+        of all vertex coordinates.
+
+        For closed 2D polygons, returns the area-weighted centroid (geometric center
+        of mass). The area weighting accounts for the shape's geometry, ensuring the
+        centroid lies at the balance point of the polygon as if it were a uniform
+        plate. This differs from the simple average of vertices, which would not
+        account for how vertices are distributed around the polygon's boundary.
+
+        Returns
+        -------
+        np.ndarray
+            Centroid coordinates, shape (2,) or (3,).
+
+        Notes
+        -----
+        For closed 2D polygons, uses the standard formula:
+        C_x = (1 / (6*A)) * sum((x_i + x_{i+1}) * (x_i * y_{i+1} - x_{i+1} * y_i))
+        C_y = (1 / (6*A)) * sum((y_i + y_{i+1}) * (x_i * y_{i+1} - x_{i+1} * y_i))
+        where A is the polygon area.
+        """
+        if not self.closed or not self._is_2d:
+            # Simple average for open polygons or 3D closed polygons
+            return np.mean(self.points, axis=0)
+
+        # Area-weighted centroid for closed 2D polygons
+        x = self.points[:, 0]
+        y = self.points[:, 1]
+        # Append first point to close polygon
+        x_closed = np.append(x, x[0])
+        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())
+
+        if area < sys.float_info.epsilon:
+            # Degenerate case: 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)
+        return np.array([cx, cy])
+
+    def area(self) -> float:
+        """Compute area enclosed by closed 2D polygon.
+
+        Uses the shoelace formula. Only valid for closed 2D polygons.
+
+        Returns
+        -------
+        float
+            Enclosed area (always positive).
+
+        Raises
+        ------
+        ValueError
+            If polygon is not closed or not 2D.
+        """
+        if not self.closed:
+            raise ValueError("Area computation requires closed polygon.")
+        if not self._is_2d:
+            raise ValueError("Area computation only valid for 2D polygons.")
+
+        x = self.points[:, 0]
+        y = self.points[:, 1]
+        # Append first point to close polygon
+        x_closed = np.append(x, x[0])
+        y_closed = np.append(y, y[0])
+        # Shoelace formula
+        area = 0.5 * np.abs(
+            np.sum(x_closed[:-1] * y_closed[1:] - x_closed[1:] * y_closed[:-1])
+        )
+        return area
+
+    def resample(
+        self, n_points: int = 100, n_iter: int = 1, inplace: bool = False
+    ) -> "Polygon":
+        """Resample polygon to have equidistant points.
+
+        Creates n_points that are approximately equidistantly spaced along
+        the cumulative Euclidean distance. Uses linear interpolation.
+
+        Parameters
+        ----------
+        n_points : int, default=100
+            Number of points in resampled polygon. Must be at least 2.
+        n_iter : int, default=1
+            Number of resampling iterations. Higher values (e.g., 3-5) provide
+            better equidistant spacing. Must be at least 1.
+        inplace : bool, default=False
+            If True, modify this polygon in-place. If False, return new polygon.
+
+        Returns
+        -------
+        Polygon
+            Resampled polygon. Returns self if inplace=True, new instance otherwise.
+        """
+        def _resample_once(p: np.ndarray, n: int, is_closed: bool) -> np.ndarray:
+            """Single resampling pass."""
+            if is_closed:
+                p_closed = np.vstack([p, p[0]])
+                d = np.cumsum(
+                    np.r_[0, np.sqrt((np.diff(p_closed, axis=0) ** 2).sum(axis=1))]
+                )
+                d_sampled = np.linspace(0, d.max(), n + 1)[:-1]
+            else:
+                d = np.cumsum(np.r_[0, np.sqrt((np.diff(p, axis=0) ** 2).sum(axis=1))])
+                d_sampled = np.linspace(0, d.max(), n)
+                p_closed = p
+
+            n_dims = p.shape[1]
+            return np.column_stack(
+                [np.interp(d_sampled, d, p_closed[:, i]) for i in range(n_dims)]
+            )
+
+        # Perform resampling n_iter times
+        points_resampled = _resample_once(self.points, n_points, self.closed)
+        for _ in range(n_iter - 1):
+            points_resampled = _resample_once(points_resampled, n_points, self.closed)
+
+        if inplace:
+            self.points = points_resampled
+            return self
         else:
-            # For open paths, use the original behavior
-            d = np.cumsum(np.r_[0, np.sqrt((np.diff(p, axis=0) ** 2).sum(axis=1))])
-            d_sampled = np.linspace(0, d.max(), n)
-            p_closed = p
-
-        # Interpolate each coordinate dimension
-        n_dims = p.shape[1]
-        return np.column_stack([
-            np.interp(d_sampled, d, p_closed[:, i]) for i in range(n_dims)
-        ])
-
-    # Perform resampling n_iter times
-    path_resampled = _resample_once(path, n_points, closed)
-    for _ in range(n_iter - 1):
-        path_resampled = _resample_once(path_resampled, n_points, closed)
-    return path_resampled
+            return Polygon(points_resampled, closed=self.closed)
+
+    def _construct_smoothing_matrix(self) -> sparse.csc_matrix:
+        """Construct smoothing matrix for Laplacian smoothing.
+
+        Creates a row-stochastic matrix where each point is connected to
+        its neighbors (previous and next point).
+
+        Returns
+        -------
+        scipy.sparse.csc_matrix
+            Sparse smoothing matrix.
+        """
+        n = self.points.shape[0]
+
+        if self.closed:
+            # For closed polygons, connect last to first
+            i = np.arange(n)
+            j_prev = np.roll(np.arange(n), 1)
+            j_next = np.roll(np.arange(n), -1)
+        else:
+            # For open polygons, first and last points have only one neighbor
+            i = np.arange(n)
+            j_prev = np.clip(np.arange(n) - 1, 0, n - 1)
+            j_next = np.clip(np.arange(n) + 1, 0, n - 1)
+
+        # Create adjacency with neighbors
+        i_all = np.concatenate([i, i])
+        j_all = np.concatenate([j_prev, j_next])
+        data = np.ones(len(i_all))
+
+        adj = sparse.csc_matrix((data, (i_all, j_all)), shape=(n, n))
+
+        # Normalize rows to create stochastic matrix
+        row_sum = np.array(adj.sum(axis=1)).ravel()
+        row_sum[row_sum == 0] = 1.0  # Avoid division by zero
+        adj = adj.multiply(1.0 / row_sum[:, np.newaxis])
+
+        return adj
+
+    def smooth_laplace(
+        self,
+        n: int = 1,
+        lambda_: float = 0.5,
+        inplace: bool = False,
+    ) -> "Polygon":
+        """Smooth polygon using Laplace smoothing.
+
+        Applies iterative smoothing: p_new = (1-lambda)*p + lambda * M*p
+        where M is the neighbor-averaging matrix.
+
+        Parameters
+        ----------
+        n : int, default=1
+            Number of smoothing iterations.
+        lambda_ : float, default=0.5
+            Diffusion speed parameter in range [0, 1].
+        inplace : bool, default=False
+            If True, modify this polygon in-place. If False, return new polygon.
+
+        Returns
+        -------
+        Polygon
+            Smoothed polygon. Returns self if inplace=True, new instance otherwise.
+        """
+        mat = self._construct_smoothing_matrix()
+        points_smooth = self.points.copy()
+
+        for _ in range(n):
+            points_smooth = (1.0 - lambda_) * points_smooth + lambda_ * mat.dot(
+                points_smooth
+            )
+
+        if inplace:
+            self.points = points_smooth
+            return self
+        else:
+            return Polygon(points_smooth, closed=self.closed)
+
+    def smooth_taubin(
+        self,
+        n: int = 1,
+        lambda_: float = 0.5,
+        mu: float = -0.53,
+        inplace: bool = False,
+    ) -> "Polygon":
+        """Smooth polygon using Taubin smoothing.
+
+        Alternates between shrinking (positive lambda) and expanding (negative mu)
+        steps to reduce shrinkage while smoothing.
+
+        Parameters
+        ----------
+        n : int, default=1
+            Number of smoothing iterations.
+        lambda_ : float, default=0.5
+            Positive diffusion parameter for shrinking step.
+        mu : float, default=-0.53
+            Negative diffusion parameter for expanding step.
+        inplace : bool, default=False
+            If True, modify this polygon in-place. If False, return new polygon.
+
+        Returns
+        -------
+        Polygon
+            Smoothed polygon. Returns self if inplace=True, new instance otherwise.
+        """
+        mat = self._construct_smoothing_matrix()
+        points_smooth = self.points.copy()
+
+        for _ in range(n):
+            # Lambda step (shrinking)
+            points_smooth = (1.0 - lambda_) * points_smooth + lambda_ * mat.dot(
+                points_smooth
+            )
+            # Mu step (expanding)
+            points_smooth = (1.0 - mu) * points_smooth + mu * mat.dot(points_smooth)
+
+        if inplace:
+            self.points = points_smooth
+            return self
+        else:
+            return Polygon(points_smooth, closed=self.closed)
+
diff --git a/lapy/tria_mesh.py b/lapy/tria_mesh.py
index de5464cd..a40e5845 100644
--- a/lapy/tria_mesh.py
+++ b/lapy/tria_mesh.py
@@ -1801,8 +1801,11 @@ def level_path(
             if n_points:
                 if get_edges:
                     raise ValueError("n_points cannot be combined with get_edges=True.")
-                path3d = polygon.resample(path3d, n_points, n_iter=3)
+                poly = polygon.Polygon(path3d, closed=False)
+                path3d = poly.resample(n_points=n_points, n_iter=3, inplace=False)
+                path3d = path3d.get_points()
             if get_edges:
                 return path3d, llength, edges_vidxs, edges_relpos
             else:
-                return path3d, llength
\ No newline at end of file
+                return path3d, llength
+
diff --git a/lapy/utils/tests/test_polygon.py b/lapy/utils/tests/test_polygon.py
index 60a6a13b..fc29c56b 100644
--- a/lapy/utils/tests/test_polygon.py
+++ b/lapy/utils/tests/test_polygon.py
@@ -3,75 +3,224 @@
 import numpy as np
 import pytest
 
-from ... import polygon
+from ... import Polygon
 
 
-class TestResample:
-    """Test cases for the polygon.resample function."""
+class TestPolygonClass:
+    """Test cases for the Polygon class."""
 
-    def test_resample_open_path(self):
-        """Test resampling a 2D open path."""
-        # Create a simple L-shaped path
-        path = np.array([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0]])
-        resampled = polygon.resample(path, n_points=10, closed=False)
+    def test_init_2d(self):
+        """Test initialization with 2D points."""
+        points = np.array([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
+        poly = Polygon(points, closed=True)
 
-        assert resampled.shape == (10, 2), "Output shape should be (10, 2)"
-        assert np.allclose(resampled[0], [0.0, 0.0]), "First point should match"
-        assert np.allclose(resampled[-1], [1.0, 1.0]), "Last point should match"
+        assert poly.is_2d(), "Should be 2D polygon"
+        assert poly.is_closed(), "Should be closed"
+        assert poly.n_points() == 4, "Should have 4 points"
 
-    def test_resample_closed_path(self):
-        """Test resampling a 2D closed path (square)."""
-        # Create a unit square (4 vertices, no duplication)
-        square = np.array([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
-        resampled = polygon.resample(square, n_points=12, closed=True)
+    def test_init_3d(self):
+        """Test initialization with 3D points."""
+        points = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [1.0, 1.0, 0.0]])
+        poly = Polygon(points, closed=False)
 
-        assert resampled.shape == (12, 2), "Output shape should be (12, 2)"
+        assert not poly.is_2d(), "Should be 3D polygon"
+        assert not poly.is_closed(), "Should be open"
+        assert poly.n_points() == 3, "Should have 3 points"
 
-        # Check that path covers the full perimeter
-        resampled_with_wrap = np.vstack([resampled, resampled[0]])
-        dists = np.sqrt(np.sum(np.diff(resampled_with_wrap, axis=0)**2, axis=1))
-        total_length = np.sum(dists)
+    def test_init_empty_raises(self):
+        """Test that empty points raise ValueError."""
+        with pytest.raises(ValueError, match="empty"):
+            Polygon(np.array([]))
 
-        assert np.allclose(total_length, 4.0, atol=1e-10), \
-            f"Closed square should have perimeter 4.0, got {total_length}"
+    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]]))
 
-    def test_resample_closed_vs_open(self):
-        """Test that closed parameter makes a difference."""
-        square = np.array([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
+    def test_length_open(self):
+        """Test length computation for open polygon."""
+        points = np.array([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0]])
+        poly = Polygon(points, closed=False)
+        length = poly.length()
 
-        resampled_open = polygon.resample(square, n_points=12, closed=False)
-        resampled_closed = polygon.resample(square, n_points=12, closed=True)
+        expected = 2.0  # 1.0 + 1.0
+        assert np.isclose(length, expected), f"Expected {expected}, got {length}"
 
-        # Calculate path lengths
-        dists_open = np.sqrt(np.sum(np.diff(resampled_open, axis=0)**2, axis=1))
-        total_open = np.sum(dists_open)
+    def test_length_closed(self):
+        """Test length computation for closed polygon (square)."""
+        square = np.array([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
+        poly = Polygon(square, closed=True)
+        length = poly.length()
+
+        expected = 4.0
+        assert np.isclose(length, expected), f"Expected {expected}, got {length}"
+
+    def test_centroid_open(self):
+        """Test centroid for open polygon."""
+        points = np.array([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
+        poly = Polygon(points, closed=False)
+        centroid = poly.centroid()
+
+        expected = np.array([0.5, 0.5])
+        assert np.allclose(centroid, expected), f"Expected {expected}, got {centroid}"
+
+    def test_centroid_closed_2d(self):
+        """Test area-weighted centroid for closed 2D polygon."""
+        # Unit square centered at origin
+        square = np.array([[-0.5, -0.5], [0.5, -0.5], [0.5, 0.5], [-0.5, 0.5]])
+        poly = Polygon(square, closed=True)
+        centroid = poly.centroid()
+
+        expected = np.array([0.0, 0.0])
+        assert np.allclose(centroid, expected, atol=1e-10), \
+            f"Expected {expected}, got {centroid}"
+
+    def test_centroid_closed_3d(self):
+        """Test centroid for closed 3D polygon (simple average)."""
+        points = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [1.0, 1.0, 0.0]])
+        poly = Polygon(points, closed=True)
+        centroid = poly.centroid()
+
+        expected = np.mean(points, axis=0)
+        assert np.allclose(centroid, expected), f"Expected {expected}, got {centroid}"
+
+    def test_area_closed_2d(self):
+        """Test area computation for closed 2D polygon."""
+        square = np.array([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
+        poly = Polygon(square, closed=True)
+        area = poly.area()
+
+        expected = 1.0
+        assert np.isclose(area, expected), f"Expected {expected}, got {area}"
+
+    def test_area_open_raises(self):
+        """Test that area computation raises for open polygon."""
+        points = np.array([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0]])
+        poly = Polygon(points, closed=False)
+
+        with pytest.raises(ValueError, match="closed polygon"):
+            poly.area()
+
+    def test_area_3d_raises(self):
+        """Test that area computation raises for 3D polygon."""
+        points = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [1.0, 1.0, 0.0]])
+        poly = Polygon(points, closed=True)
+
+        with pytest.raises(ValueError, match="2D polygons"):
+            poly.area()
+
+    def test_resample_open(self):
+        """Test resampling an open polygon."""
+        points = np.array([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0]])
+        poly = Polygon(points, closed=False)
+        resampled = poly.resample(n_points=10, inplace=False)
+
+        assert resampled.n_points() == 10, "Should have 10 points"
+        assert not resampled.is_closed(), "Should remain open"
+        assert np.allclose(resampled.get_points()[0], [0.0, 0.0]), \
+            "First point should match"
+        assert np.allclose(resampled.get_points()[-1], [1.0, 1.0]), \
+            "Last point should match"
+
+    def test_resample_closed(self):
+        """Test resampling a closed polygon."""
+        square = np.array([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
+        poly = Polygon(square, closed=True)
+        resampled = poly.resample(n_points=12, inplace=False)
 
-        resampled_closed_wrap = np.vstack([resampled_closed, resampled_closed[0]])
-        dists_closed = np.sqrt(np.sum(np.diff(resampled_closed_wrap, axis=0)**2, axis=1))
-        total_closed = np.sum(dists_closed)
+        assert resampled.n_points() == 12, "Should have 12 points"
+        assert resampled.is_closed(), "Should remain closed"
+        # Check perimeter
+        length = resampled.length()
+        assert np.isclose(length, 4.0, atol=1e-10), \
+            f"Perimeter should be 4.0, got {length}"
 
-        # Closed path should cover full perimeter, open path should not
-        assert total_open < total_closed, "Open path should be shorter"
-        assert np.allclose(total_closed, 4.0, atol=1e-10), "Closed perimeter should be 4.0"
+    def test_resample_inplace(self):
+        """Test in-place resampling."""
+        points = np.array([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0]])
+        poly = Polygon(points, closed=False)
+        result = poly.resample(n_points=10, inplace=True)
+
+        assert result is poly, "Should return self when inplace=True"
+        assert poly.n_points() == 10, "Should have 10 points after in-place resample"
 
     def test_resample_iterative(self):
-        """Test that iterative resampling improves spacing uniformity."""
-        # Create a path with non-uniform spacing
-        path = np.array([[0.0, 0.0], [0.1, 0.0], [1.0, 0.0], [1.0, 1.0]])
+        """Test that iterative resampling improves uniformity."""
+        # Create path with non-uniform spacing
+        points = np.array([[0.0, 0.0], [0.1, 0.0], [1.0, 0.0], [1.0, 1.0]])
+        poly = Polygon(points, closed=False)
 
-        result_single = polygon.resample(path, n_points=20, n_iter=1, closed=False)
-        result_multiple = polygon.resample(path, n_points=20, n_iter=5, closed=False)
+        result1 = poly.resample(n_points=20, n_iter=1, inplace=False)
+        result5 = poly.resample(n_points=20, n_iter=5, inplace=False)
 
         # Calculate spacing uniformity
-        dists_single = np.sqrt(np.sum(np.diff(result_single, axis=0)**2, axis=1))
-        dists_multiple = np.sqrt(np.sum(np.diff(result_multiple, axis=0)**2, axis=1))
-
-        # Multiple iterations should have more uniform spacing (lower std dev)
-        assert np.std(dists_multiple) <= np.std(dists_single), \
-            "Iterative resampling should improve uniformity"
-
+        pts1 = result1.get_points()
+        pts5 = result5.get_points()
+        dists1 = np.sqrt(np.sum(np.diff(pts1, axis=0)**2, axis=1))
+        dists5 = np.sqrt(np.sum(np.diff(pts5, axis=0)**2, axis=1))
+
+        assert np.std(dists5) <= np.std(dists1), \
+            "More iterations should improve uniformity"
+
+    def test_smooth_laplace_open(self):
+        """Test Laplace smoothing on open polygon."""
+        # Create a jagged path
+        points = np.array([[0.0, 0.0], [0.5, 0.5], [1.0, 0.0], [1.5, 0.5], [2.0, 0.0]])
+        poly = Polygon(points, closed=False)
+        smoothed = poly.smooth_laplace(n=3, lambda_=0.5, inplace=False)
+
+        assert smoothed.n_points() == poly.n_points(), \
+            "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"
+        assert np.allclose(smoothed.get_points()[-1], points[-1]), \
+            "Last point should not change much"
+
+    def test_smooth_laplace_closed(self):
+        """Test Laplace smoothing on closed polygon."""
+        # Create a slightly irregular square
+        square = np.array([
+            [0.0, 0.0], [1.0, 0.1], [1.0, 1.0], [0.1, 1.0]
+        ])
+        poly = Polygon(square, closed=True)
+        smoothed = poly.smooth_laplace(n=5, lambda_=0.5, inplace=False)
+
+        assert smoothed.n_points() == 4, "Should preserve number of points"
+        assert smoothed.is_closed(), "Should remain closed"
+
+    def test_smooth_laplace_inplace(self):
+        """Test in-place Laplace smoothing."""
+        points = np.array([[0.0, 0.0], [0.5, 0.5], [1.0, 0.0]])
+        poly = Polygon(points, closed=False)
+        original_points = poly.get_points().copy()
+        result = poly.smooth_laplace(n=1, lambda_=0.5, inplace=True)
+
+        assert result is poly, "Should return self when inplace=True"
+        assert not np.allclose(poly.get_points(), original_points), \
+            "Points should be modified in-place"
+
+    def test_smooth_taubin(self):
+        """Test Taubin smoothing on polygon."""
+        # Create a jagged path
+        points = np.array([[0.0, 0.0], [0.5, 0.5], [1.0, 0.0], [1.5, 0.5], [2.0, 0.0]])
+        poly = Polygon(points, closed=False)
+        smoothed = poly.smooth_taubin(n=3, lambda_=0.5, mu=-0.53, inplace=False)
+
+        assert smoothed.n_points() == poly.n_points(), \
+            "Should preserve number of points"
+        # Taubin should preserve overall shape better than pure Laplace
+        assert not np.allclose(smoothed.get_points(), points), \
+            "Points should be smoothed"
+
+    def test_smooth_taubin_inplace(self):
+        """Test in-place Taubin smoothing."""
+        points = np.array([[0.0, 0.0], [0.5, 0.5], [1.0, 0.0]])
+        poly = Polygon(points, closed=False)
+        result = poly.smooth_taubin(n=1, inplace=True)
+
+        assert result is poly, "Should return self when inplace=True"
 
 
 if __name__ == "__main__":
     pytest.main([__file__, "-v"])
-

From 9e066de1c6dc6167f8ca3fc11fa297a982239b46 Mon Sep 17 00:00:00 2001
From: Martin Reuter <reuter@mit.edu>
Date: Fri, 19 Dec 2025 02:21:19 +0100
Subject: [PATCH 02/14] Fix boundary point in open polygons when smoothing

---
 lapy/polygon.py | 30 +++++++++++++++++++++---------
 1 file changed, 21 insertions(+), 9 deletions(-)

diff --git a/lapy/polygon.py b/lapy/polygon.py
index 273e9747..e3b4b672 100644
--- a/lapy/polygon.py
+++ b/lapy/polygon.py
@@ -274,7 +274,8 @@ def _construct_smoothing_matrix(self) -> sparse.csc_matrix:
         """Construct smoothing matrix for Laplacian smoothing.
 
         Creates a row-stochastic matrix where each point is connected to
-        its neighbors (previous and next point).
+        its neighbors (previous and next point). For open polygons, boundary
+        points (first and last) are kept fixed.
 
         Returns
         -------
@@ -288,16 +289,22 @@ def _construct_smoothing_matrix(self) -> sparse.csc_matrix:
             i = np.arange(n)
             j_prev = np.roll(np.arange(n), 1)
             j_next = np.roll(np.arange(n), -1)
+
+            # Create adjacency with neighbors
+            i_all = np.concatenate([i, i])
+            j_all = np.concatenate([j_prev, j_next])
+            data = np.ones(len(i_all))
         else:
-            # For open polygons, first and last points have only one neighbor
-            i = np.arange(n)
-            j_prev = np.clip(np.arange(n) - 1, 0, n - 1)
-            j_next = np.clip(np.arange(n) + 1, 0, n - 1)
+            # For open polygons, first and last points stay fixed
+            # Interior points are connected to their neighbors
+            i_interior = np.arange(1, n - 1)
+            j_prev_interior = i_interior - 1
+            j_next_interior = i_interior + 1
 
-        # Create adjacency with neighbors
-        i_all = np.concatenate([i, i])
-        j_all = np.concatenate([j_prev, j_next])
-        data = np.ones(len(i_all))
+            # Create adjacency for interior points
+            i_all = np.concatenate([i_interior, i_interior])
+            j_all = np.concatenate([j_prev_interior, j_next_interior])
+            data = np.ones(len(i_all))
 
         adj = sparse.csc_matrix((data, (i_all, j_all)), shape=(n, n))
 
@@ -306,6 +313,11 @@ def _construct_smoothing_matrix(self) -> sparse.csc_matrix:
         row_sum[row_sum == 0] = 1.0  # Avoid division by zero
         adj = adj.multiply(1.0 / row_sum[:, np.newaxis])
 
+        # For open polygons, add identity for boundary points
+        if not self.closed:
+            adj[0, 0] = 1.0
+            adj[n - 1, n - 1] = 1.0
+
         return adj
 
     def smooth_laplace(

From 44032aadb1f09d960062abe93afc2c456a79c51c Mon Sep 17 00:00:00 2001
From: Martin Reuter <reuter@mit.edu>
Date: Fri, 19 Dec 2025 02:26:54 +0100
Subject: [PATCH 03/14] Fix sparse matrix construction

---
 lapy/polygon.py | 37 ++++++++++++++++++-------------------
 1 file changed, 18 insertions(+), 19 deletions(-)

diff --git a/lapy/polygon.py b/lapy/polygon.py
index e3b4b672..8c5ee505 100644
--- a/lapy/polygon.py
+++ b/lapy/polygon.py
@@ -271,7 +271,7 @@ def _resample_once(p: np.ndarray, n: int, is_closed: bool) -> np.ndarray:
             return Polygon(points_resampled, closed=self.closed)
 
     def _construct_smoothing_matrix(self) -> sparse.csc_matrix:
-        """Construct smoothing matrix for Laplacian smoothing.
+        """Construct smoothing matrix for Laplace smoothing.
 
         Creates a row-stochastic matrix where each point is connected to
         its neighbors (previous and next point). For open polygons, boundary
@@ -294,30 +294,29 @@ def _construct_smoothing_matrix(self) -> sparse.csc_matrix:
             i_all = np.concatenate([i, i])
             j_all = np.concatenate([j_prev, j_next])
             data = np.ones(len(i_all))
-        else:
-            # For open polygons, first and last points stay fixed
-            # Interior points are connected to their neighbors
-            i_interior = np.arange(1, n - 1)
-            j_prev_interior = i_interior - 1
-            j_next_interior = i_interior + 1
-
-            # Create adjacency for interior points
-            i_all = np.concatenate([i_interior, i_interior])
-            j_all = np.concatenate([j_prev_interior, j_next_interior])
-            data = np.ones(len(i_all))
 
-        adj = sparse.csc_matrix((data, (i_all, j_all)), shape=(n, n))
+            adj = sparse.csc_matrix((data, (i_all, j_all)), shape=(n, n))
 
-        # Normalize rows to create stochastic matrix
-        row_sum = np.array(adj.sum(axis=1)).ravel()
-        row_sum[row_sum == 0] = 1.0  # Avoid division by zero
-        adj = adj.multiply(1.0 / row_sum[:, np.newaxis])
+            # Normalize rows to create stochastic matrix
+            row_sum = np.array(adj.sum(axis=1)).ravel()
+            row_sum[row_sum == 0] = 1.0  # Avoid division by zero
+            adj = adj.multiply(1.0 / row_sum[:, np.newaxis])
+        else:
+            # For open polygons, use LIL format for easier construction
+            adj = sparse.lil_matrix((n, n))
 
-        # For open polygons, add identity for boundary points
-        if not self.closed:
+            # Set identity for boundary points (they stay fixed)
             adj[0, 0] = 1.0
             adj[n - 1, n - 1] = 1.0
 
+            # For interior points, connect to neighbors
+            for i in range(1, n - 1):
+                adj[i, i - 1] = 0.5
+                adj[i, i + 1] = 0.5
+
+            # Convert to CSC for efficient operations
+            adj = adj.tocsc()
+
         return adj
 
     def smooth_laplace(

From 1495ffc80ab88c7451e64c5673525a8d823c1811 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Fri, 19 Dec 2025 06:47:29 +0000
Subject: [PATCH 04/14] Initial plan


From c7e331856794483ea59d2dd110d3f4e94a19897c Mon Sep 17 00:00:00 2001
From: Martin Reuter <m-reuter@users.noreply.github.com>
Date: Fri, 19 Dec 2025 07:51:28 +0100
Subject: [PATCH 05/14] Apply suggestions from code review

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
---
 lapy/__init__.py                 |  2 +-
 lapy/polygon.py                  | 30 ++++++++++++++++++++----------
 lapy/utils/tests/test_polygon.py |  4 ++--
 3 files changed, 23 insertions(+), 13 deletions(-)

diff --git a/lapy/__init__.py b/lapy/__init__.py
index 30fb8830..0fccfc2c 100644
--- a/lapy/__init__.py
+++ b/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
diff --git a/lapy/polygon.py b/lapy/polygon.py
index 8c5ee505..1ffd80d1 100644
--- a/lapy/polygon.py
+++ b/lapy/polygon.py
@@ -4,15 +4,10 @@
 various geometric operations including resampling, smoothing, and metric computations.
 """
 
-import logging
 import sys
 
 import numpy as np
 from scipy import sparse
-
-logger = logging.getLogger(__name__)
-
-
 class Polygon:
     """Class representing a polygon path (open or closed).
 
@@ -177,14 +172,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 +236,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:
@@ -344,6 +343,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 +390,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()
 
diff --git a/lapy/utils/tests/test_polygon.py b/lapy/utils/tests/test_polygon.py
index fc29c56b..657b23bf 100644
--- a/lapy/utils/tests/test_polygon.py
+++ b/lapy/utils/tests/test_polygon.py
@@ -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."""

From 60e25c6839a003d3327703c94da69bc20ff68709 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Fri, 19 Dec 2025 06:52:30 +0000
Subject: [PATCH 06/14] Add documentation and tests for edge cases in smoothing
 matrix

Co-authored-by: m-reuter <8526484+m-reuter@users.noreply.github.com>
---
 lapy/polygon.py                  | 11 ++++++++++
 lapy/utils/tests/test_polygon.py | 36 ++++++++++++++++++++++++++++++++
 2 files changed, 47 insertions(+)

diff --git a/lapy/polygon.py b/lapy/polygon.py
index 8c5ee505..699fa122 100644
--- a/lapy/polygon.py
+++ b/lapy/polygon.py
@@ -277,6 +277,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
diff --git a/lapy/utils/tests/test_polygon.py b/lapy/utils/tests/test_polygon.py
index fc29c56b..5181723a 100644
--- a/lapy/utils/tests/test_polygon.py
+++ b/lapy/utils/tests/test_polygon.py
@@ -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"])

From ebaccad3bf0b43f0572c721900276e33fc1297f6 Mon Sep 17 00:00:00 2001
From: Martin Reuter <reuter@mit.edu>
Date: Fri, 19 Dec 2025 08:01:09 +0100
Subject: [PATCH 07/14] fix include order

---
 lapy/polygon.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/lapy/polygon.py b/lapy/polygon.py
index 1ffd80d1..50c2fe38 100644
--- a/lapy/polygon.py
+++ b/lapy/polygon.py
@@ -8,6 +8,8 @@
 
 import numpy as np
 from scipy import sparse
+
+
 class Polygon:
     """Class representing a polygon path (open or closed).
 

From 3154eace067e5346632b929ad6d89bb9faa42f34 Mon Sep 17 00:00:00 2001
From: Martin Reuter <m-reuter@users.noreply.github.com>
Date: Fri, 19 Dec 2025 08:04:41 +0100
Subject: [PATCH 08/14] Update lapy/polygon.py

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
---
 lapy/polygon.py | 21 +++++++++++++++++----
 1 file changed, 17 insertions(+), 4 deletions(-)

diff --git a/lapy/polygon.py b/lapy/polygon.py
index 50c2fe38..c43e1613 100644
--- a/lapy/polygon.py
+++ b/lapy/polygon.py
@@ -67,14 +67,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")

From e1845cd5bcd3d5182a4d13993b97bc4307d2a540 Mon Sep 17 00:00:00 2001
From: Martin Reuter <m-reuter@users.noreply.github.com>
Date: Fri, 19 Dec 2025 07:51:28 +0100
Subject: [PATCH 09/14] Apply suggestions from code review

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
---
 lapy/__init__.py                 |  2 +-
 lapy/polygon.py                  | 30 ++++++++++++++++++++----------
 lapy/utils/tests/test_polygon.py |  4 ++--
 3 files changed, 23 insertions(+), 13 deletions(-)

diff --git a/lapy/__init__.py b/lapy/__init__.py
index 30fb8830..0fccfc2c 100644
--- a/lapy/__init__.py
+++ b/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
diff --git a/lapy/polygon.py b/lapy/polygon.py
index 699fa122..8ceead88 100644
--- a/lapy/polygon.py
+++ b/lapy/polygon.py
@@ -4,15 +4,10 @@
 various geometric operations including resampling, smoothing, and metric computations.
 """
 
-import logging
 import sys
 
 import numpy as np
 from scipy import sparse
-
-logger = logging.getLogger(__name__)
-
-
 class Polygon:
     """Class representing a polygon path (open or closed).
 
@@ -177,14 +172,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 +236,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:
@@ -355,6 +354,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()
 
@@ -397,6 +401,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()
 
diff --git a/lapy/utils/tests/test_polygon.py b/lapy/utils/tests/test_polygon.py
index 5181723a..4258931c 100644
--- a/lapy/utils/tests/test_polygon.py
+++ b/lapy/utils/tests/test_polygon.py
@@ -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."""

From 4e453e5078a1f58ec24c6ebc4470d3c9b9789c14 Mon Sep 17 00:00:00 2001
From: Martin Reuter <reuter@mit.edu>
Date: Fri, 19 Dec 2025 08:01:09 +0100
Subject: [PATCH 10/14] fix include order

---
 lapy/polygon.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/lapy/polygon.py b/lapy/polygon.py
index 8ceead88..c0989e0c 100644
--- a/lapy/polygon.py
+++ b/lapy/polygon.py
@@ -8,6 +8,8 @@
 
 import numpy as np
 from scipy import sparse
+
+
 class Polygon:
     """Class representing a polygon path (open or closed).
 

From 28c2468f0660324c867a57276f3f35b44a5778fc Mon Sep 17 00:00:00 2001
From: Martin Reuter <m-reuter@users.noreply.github.com>
Date: Fri, 19 Dec 2025 08:04:41 +0100
Subject: [PATCH 11/14] Update lapy/polygon.py

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
---
 lapy/polygon.py | 21 +++++++++++++++++----
 1 file changed, 17 insertions(+), 4 deletions(-)

diff --git a/lapy/polygon.py b/lapy/polygon.py
index c0989e0c..7724699f 100644
--- a/lapy/polygon.py
+++ b/lapy/polygon.py
@@ -67,14 +67,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")

From 83dbfb060f3620432b1175e32fad3286bbf7c69d Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Fri, 19 Dec 2025 06:47:29 +0000
Subject: [PATCH 12/14] Initial plan


From 8a19940ae126584cdf4524793e981ff61ac2bdad Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Fri, 19 Dec 2025 07:12:14 +0000
Subject: [PATCH 13/14] Fix logger import and test after rebase

Co-authored-by: m-reuter <8526484+m-reuter@users.noreply.github.com>
---
 lapy/polygon.py                  | 3 +++
 lapy/utils/tests/test_polygon.py | 2 +-
 2 files changed, 4 insertions(+), 1 deletion(-)

diff --git a/lapy/polygon.py b/lapy/polygon.py
index 7724699f..568b8a86 100644
--- a/lapy/polygon.py
+++ b/lapy/polygon.py
@@ -4,11 +4,14 @@
 various geometric operations including resampling, smoothing, and metric computations.
 """
 
+import logging
 import sys
 
 import numpy as np
 from scipy import sparse
 
+logger = logging.getLogger(__name__)
+
 
 class Polygon:
     """Class representing a polygon path (open or closed).
diff --git a/lapy/utils/tests/test_polygon.py b/lapy/utils/tests/test_polygon.py
index 4258931c..869c7a4b 100644
--- a/lapy/utils/tests/test_polygon.py
+++ b/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."""

From c9a525c575eb5e5558134f29624ba619c82afc1f Mon Sep 17 00:00:00 2001
From: Martin Reuter <reuter@mit.edu>
Date: Fri, 19 Dec 2025 08:16:40 +0100
Subject: [PATCH 14/14] add logger

---
 lapy/polygon.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/lapy/polygon.py b/lapy/polygon.py
index c43e1613..62488084 100644
--- a/lapy/polygon.py
+++ b/lapy/polygon.py
@@ -3,12 +3,13 @@
 This module provides a Polygon class for processing 2D and 3D polygon paths with
 various geometric operations including resampling, smoothing, and metric computations.
 """
-
+import logging
 import sys
 
 import numpy as np
 from scipy import sparse
 
+logger = logging.getLogger(__name__)
 
 class Polygon:
     """Class representing a polygon path (open or closed).