"""Implementation of a 2D ordered closed polygon class."""
import numpy as np
from torx.autodoc_decorators_m import autodoc_class
from scipy.interpolate import splprep, splev
from torx.geometry.polygon_2d_m import Polygon2D
from torx.geometry.line_geometry_m import order_points
[docs]
@autodoc_class
class OrderedPolygon2D(Polygon2D):
"""2D polygon which ensures that points are ordered counterclockwise."""
[docs]
def __init__(self, x_points: np.ndarray, y_points: np.ndarray,
start_idx: int=0, k: int=20, n_interp: int=0, smooth: float=None):
"""
Initialize the polygon with the vertices given.
The points are ordered counterclockwise using a nearest-neighbor
algorithm starting from the point with index start_idx and looking for
up to k nearest neighbors.
If n_interp > 0, the polygon points are interpolated to n_interp points
using a periodic spline representation with optional smoothing parameter
'smooth'. Note: smoothing can be useful but high values may lead to very
distorted polygons.
"""
super().__init__(x_points, y_points)
# Order the polygon points using a nearest neighbor algorithm
self.order(start_idx, k)
# Optional: interpolate if n_interp > 0
if n_interp > 0:
x_interp, y_interp = self.interpolate(n_interp, smooth)
self.x_points = x_interp
self.y_points = y_interp
# Make sure the polygon points are ordered counterclockwise
self._ensure_ccw()
[docs]
def order(self, start_idx: int=0, k: int=20):
"""Order the polygon points using a nearest-neighbor algorithm."""
self.x_points, self.y_points = order_points(self.x_points, self.y_points, start_idx, k)
self.ordered = True
def _ensure_ccw(self):
"""Ensure that the polygon points are ordered counterclockwise."""
if not hasattr(self, 'ccw') or not self.ccw:
if not self.ordered:
order(self)
if self.signed_area() < 0.0:
# Reverse the order of the points but keep the first point fixed
self.x_points = np.concatenate(([self.x_points[0]], self.x_points[:0:-1]))
self.y_points = np.concatenate(([self.y_points[0]], self.y_points[:0:-1]))
self.ccw = True
[docs]
def interpolate(self, n_interp: int=100, smooth: float=None):
"""
Interpolate the polygon points using a periodic spline representation.
Returns arrays of n_interp interpolated x and y points.
"""
tck, u = splprep([self.x_points, self.y_points], s=smooth, per=True)
self.tck = tck
# Check if polygon is closed
if not np.isclose(self.x_points[0], self.x_points[-1]) or \
not np.isclose(self.y_points[0], self.y_points[-1]):
u_new = np.linspace(0, 1, n_interp, endpoint=False)
else:
u_new = np.linspace(0, 1, n_interp, endpoint=True)
interpolated_x, interpolated_y = splev(u_new, tck)
return np.asarray(interpolated_x), np.asarray(interpolated_y)
@property
def arclength(self):
"""Return the cumulative arclength of the polygon."""
self.repeat_first_point()
if not hasattr(self, '_arclength'):
ds = np.sqrt(np.diff(self.x_points)**2 + np.diff(self.y_points)**2)
self._arclength = np.concatenate([[0], np.cumsum(ds)])
return self._arclength
@property
def centroid(self):
"""
Return the centroid (center of mass) of the polygon.
Shape: (2,)
Uses the polygon centroid formula https://en.wikipedia.org/wiki/Centroid#Of_a_polygon
"""
if not hasattr(self, '_centroid'):
x = self.x_points
xp1 = np.roll(x, -1)
y = self.y_points
yp1 = np.roll(y, -1)
n = len(x)
cx = np.sum((x + xp1) * (x * yp1 - xp1 * y)) / (6 * self.signed_area())
cy = np.sum((yp1 + y) * (x * yp1 - xp1 * y)) / (6 * self.signed_area())
self._centroid = np.array([cx, cy])
return self._centroid
[docs]
@classmethod
def from_2d_array(cls, data_slice: np.ndarray,
start_idx: int=0, k: int=20, n_interp: int=0, smooth: float=None):
"""Create an instance from a single 2D array slice (npoints, xy)."""
return Polygon2D.from_2d_array(data_slice).to_ordered(start_idx, k, n_interp, smooth)
