"""Implementation of a 2D ordered closed polygon class."""
from functools import cached_property
import numpy as np
from scipy.interpolate import splprep, splev, PchipInterpolator
from .polygon_2d_m import Polygon2D
from .line_geometry_m import order_points, compute_arc_length
[docs]
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,
resample_method: str="spline", make_simple: bool=False,
order: bool=True, invert_polygon: bool=False):
"""
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. With order=False the points are taken to
be in order already and only the orientation is made counterclockwise.
If n_interp > 0, the polygon points are resampled to n_interp points,
evenly spaced along the curve, with a periodic cubic spline
(resample_method='spline', default) or a pchip that does not
overshoot (resample_method='pchip', see resample). Resampling never
smooths; use resample with smooth > 0, or
PoincarePolygonSet.smooth_surfaces, for that.
If make_simple is True, any self-intersection loop (for example an outer
scrape-off layer finger doubling back) is removed after resampling and
the polygon is re-resampled to n_interp points.
Parameters
----------
x_points, y_points : numpy.ndarray
Polygon vertices.
start_idx : int, optional
Index of the ordering start point.
k : int, optional
Number of nearest neighbors used for ordering.
n_interp : int, optional
Number of resampled points; 0 disables resampling.
resample_method : {'spline', 'pchip'}, optional
Resampling scheme applied when n_interp > 0.
make_simple : bool, optional
Remove self-intersection loops after resampling.
order : bool, optional
Order the points with the nearest-neighbor algorithm (default
True). Set to False for points that are already in sequence,
such as a loop split off this polygon, where reordering could
scramble points at sharp tips. start_idx and k are then
unused.
invert_polygon : bool, optional
Treat points outside the polygon as inside (see Polygon2D).
"""
super().__init__(x_points, y_points, invert_polygon=invert_polygon)
if order:
# Order the polygon points using a nearest neighbor algorithm
self.order(start_idx, k)
else:
self.ordered = True
if n_interp > 0:
self.resample(n_interp, resample_method)
if make_simple and not self.is_simple():
simple = self.return_as_simple()
self.x_points = simple.x_points
self.y_points = simple.y_points
if n_interp > 0:
self.resample(n_interp, resample_method)
# Make sure the polygon points are ordered counterclockwise
self._ensure_ccw()
def resample(self, n_interp: int, method: str="spline",
smooth: float=0.0):
"""
Resample the polygon in place to n_interp points.
Both schemes space the points evenly along the curve and keep
the orientation of the polygon. method='spline' evaluates a
periodic cubic spline (see interpolate), optionally smoothed.
It can overshoot at corners and, on coarse thin fingers, cross
the gap and make the polygon self-intersect. method='pchip' does
not overshoot, so straight edges stay straight and no crossings
are added. Cached quantities derived from the old points
(arclength, centroid) are dropped.
Parameters
----------
n_interp : int
Number of resampled points.
method : {'spline', 'pchip'}, optional
Resampling scheme.
smooth : float, optional
Dimensionless spline smoothing strength (see interpolate),
only for method='spline'. The default of 0.0 makes the
spline pass through every point.
Raises
------
ValueError
If method is not 'spline' or 'pchip', or if smooth > 0 is
combined with method='pchip'.
"""
if method == "pchip" and smooth > 0.0:
raise ValueError(
"Smoothing is only available with method='spline', got "
f"smooth={smooth} with method='pchip'.")
if method == "pchip":
x_points, y_points = self._resample_pchip(n_interp)
elif method == "spline":
x_points, y_points = self.interpolate(n_interp, smooth)
else:
raise ValueError(
f"Unknown resample method '{method}'; choose 'spline' or "
"'pchip'.")
self.x_points, self.y_points = x_points, y_points
self.clear_cached_properties()
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
self.clear_cached_properties()
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:
self.order()
if self.signed_area() < 0.0:
# Reverse the point order 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.clear_cached_properties()
self.ccw = True
def _new_from_loop(self, x_points: np.ndarray,
y_points: np.ndarray) -> "OrderedPolygon2D":
"""
Build an ordered polygon from one loop of this polygon.
A loop split off at a self-crossing is a contiguous run of the
already ordered vertices, so it is built with order=False:
reordering could scramble points at the sharp tip left by the
split and bring a crossing back. The constructor still restores
the counterclockwise orientation, since the two loops of a
crossing run in opposite directions.
Parameters
----------
x_points, y_points : numpy.ndarray
Vertices of the loop, in order.
Returns
-------
OrderedPolygon2D
The loop, ordered counterclockwise.
"""
return type(self)(x_points, y_points, order=False,
invert_polygon=self.invert_polygon)
def interpolate(self, n_interp: int=100, smooth: float=0.0):
"""
Interpolate the polygon points using a periodic spline representation.
Returns arrays of n_interp interpolated x and y points. splprep
parametrizes the curve by chord length, so the points are spaced
evenly along it.
The smoothing strength is dimensionless, so one value behaves the
same for every polygon size and point count. It is the allowed
RMS distance of the spline from the input points, as a fraction
of the equivalent radius sqrt(|A| / pi) of the polygon. It is
converted to the splprep residual bound
s = m * (smooth * sqrt(|A| / pi))^2 for m input points.
Parameters
----------
n_interp : int, optional
Number of interpolated points to return.
smooth : float, optional
Dimensionless smoothing strength. The default of 0.0 makes the
spline pass through every input point. A value near the noise
amplitude relative to the equivalent radius (typically
0.005 to 0.05) removes the noise; much larger values flatten
real features.
Returns
-------
tuple of numpy.ndarray
The interpolated x and y point arrays.
"""
# splprep requires the spline degree to be strictly less than the
# number of points, so clamp it for small polygons (cubic otherwise).
spline_degree = min(3, len(self.x_points) - 1)
equivalent_radius = np.sqrt(abs(self.signed_area()) / np.pi)
s = len(self.x_points) * (smooth * equivalent_radius)**2
x, y = self.x_points, self.y_points
is_closed = np.isclose(x[0], x[-1]) and np.isclose(y[0], y[-1])
# With per=True splprep ignores the last point and closes the
# curve back to the first, so an open polygon must be closed
# explicitly or its last vertex is lost.
if not is_closed:
x, y = np.append(x, x[0]), np.append(y, y[0])
tck, u = splprep([x, y], s=s, per=True, k=spline_degree)
self.tck = tck
u_new = np.linspace(0, 1, n_interp, endpoint=is_closed)
interpolated_x, interpolated_y = splev(u_new, tck)
return np.asarray(interpolated_x), np.asarray(interpolated_y)
def _resample_pchip(self, n_interp: int=100):
"""
Resample the polygon evenly in arc length with a pchip.
x and y are each fitted with a pchip in the cumulative arc
length. Unlike a cubic spline, a pchip does not overshoot, so
no self-intersections are added at corners or thin fingers.
The polygon must already be ordered.
Parameters
----------
n_interp : int, optional
Number of resampled points.
Returns
-------
tuple of numpy.ndarray
The resampled x and y point arrays.
Raises
------
RuntimeError
If the polygon has not been ordered yet (arc-length
parametrization is meaningless on unordered points).
"""
if not getattr(self, "ordered", False):
raise RuntimeError(
"Pchip resampling requires an ordered polygon; "
"call order() first.")
x = np.asarray(self.x_points, float)
y = np.asarray(self.y_points, float)
# Close the loop for a continuous arc-length parametrization.
if not np.isclose(x[0], x[-1]) or not np.isclose(y[0], y[-1]):
x = np.append(x, x[0])
y = np.append(y, y[0])
s = compute_arc_length(x, y)
keep = np.concatenate([[True], np.diff(s) > 1e-12])
s, x, y = s[keep], x[keep], y[keep]
length = s[-1]
# Pchip uses one-sided slopes at the ends of its data, so we pad one
# point at each end to avoid a corner at the seam.
s_padded = np.concatenate([[s[-2] - length], s, [s[1] + length]])
x_padded = np.concatenate([[x[-2]], x, [x[1]]])
y_padded = np.concatenate([[y[-2]], y, [y[1]]])
s_req = np.linspace(0.0, length, n_interp, endpoint=False)
return (PchipInterpolator(s_padded, x_padded)(s_req),
PchipInterpolator(s_padded, y_padded)(s_req))
@cached_property
def arclength(self):
"""Return the cumulative arclength of the polygon."""
self.repeat_first_point()
return compute_arc_length(self.x_points, self.y_points)
@cached_property
def centroid(self):
"""
Return the centroid (center of mass) of the polygon.
Shape: (2,)
Uses the polygon centroid formula. See
https://en.wikipedia.org/wiki/Centroid#Of_a_polygon
"""
x = self.x_points
xp1 = np.roll(x, -1)
y = self.y_points
yp1 = np.roll(y, -1)
cross = x * yp1 - xp1 * y
factor = 6 * self.signed_area()
cx = np.sum((x + xp1) * cross) / factor
cy = np.sum((yp1 + y) * cross) / factor
return np.array([cx, cy])
[docs]
@classmethod
def from_2d_array(cls, data_slice: np.ndarray, start_idx: int=0,
k: int=20, n_interp: int=0,
resample_method: str="spline", make_simple: bool=False):
"""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, resample_method, make_simple)
