"""Implementation of a 2D closed polygon class."""
import warnings
from functools import cached_property
import numpy as np
from skimage.measure import points_in_poly
from .line_geometry_m import first_self_intersection
from torx.decorators import AutoDocClass
[docs]
class Polygon2D(AutoDocClass):
"""2D polygon class on a poloidal plane."""
[docs]
def __init__(self, x_points: np.ndarray, y_points: np.ndarray,
invert_polygon: bool=False):
"""
Initialize the polygon with the vertices.
The polygon is treated as closed, though the first and last points do
not have to be the same.
If invert_polygon is true, points outside the polygon are considered
"inside," i.e. will return True from the routine points_inside.
"""
self.x_points = np.array(x_points).flatten()
self.y_points = np.array(y_points).flatten()
self.invert_polygon = invert_polygon
assert self.x_points.shape == self.y_points.shape
assert len(self.x_points) == self.y_points.size
[docs]
def __getitem__(self, index):
"""
Allow array indexing.
Like:
polygon[i] -> (x, y) at index i
polygon[i, 0] -> x at index i
polygon[i, 1] -> y at index i
"""
if isinstance(index, tuple):
i, j = index
if j == 0:
return self.x_points[i]
elif j == 1:
return self.y_points[i]
else:
raise IndexError("Index must be 0 or 1 for second dimension")
else:
# Handle single index
x = self.x_points[index]
y = self.y_points[index]
return np.column_stack([x, y])
def flip_Z(self):
"""
Flips the vertical direction.
To interface with NUMERICAL equilibria which use this for field
reversal. Also reverses the vertex arrays to
ensure that the points remain defined in the same counter-/clockwise
orientation as before
"""
self.x_points = self.x_points[::-1]
self.y_points = -self.y_points[::-1]
self.clear_cached_properties()
def clear_cached_properties(self):
"""
Clear the computed values of all cached properties.
Quantities decorated with @cached_property (for example the
centroid of an OrderedPolygon2D) are computed once from the
vertices and stored on the instance. Every method that changes
the vertices must call this, so that they are recomputed from
the new points on next access. The whole class hierarchy is
searched, so properties defined on subclasses are cleared too.
"""
for cls in type(self).__mro__:
for name, value in vars(cls).items():
if isinstance(value, cached_property):
self.__dict__.pop(name, None)
@property
def verts(self):
"""
Return the point vertices in a (M, 2) array.
Compatible with skimage.measure.points_in_poly
"""
return np.column_stack((self.x_points, self.y_points))
def points_inside(self, x_tests, y_tests):
"""Check whether points are inside the polygon."""
assert x_tests.shape == y_tests.shape
original_shape = x_tests.shape
points = np.column_stack((x_tests.flatten(), y_tests.flatten()))
mask = points_in_poly(points=points,
verts=self.verts).reshape(original_shape)
if not (self.invert_polygon):
return mask
else:
return np.logical_not(mask)
def point_inside(self, x_test, y_test):
"""Return whether a test point is within the polygon."""
return self.points_inside(np.atleast_1d(x_test),
np.atleast_1d(y_test))[0]
def signed_area(self):
"""
Return the signed area of a non-intersecting polygon.
If the sign is positive, points are counterclockwise. Otherwise, points
are clockwise.
See https://mathworld.wolfram.com/PolygonArea.html
"""
area = np.sum( self.x_points * np.roll(self.y_points, -1) \
- self.y_points * np.roll(self.x_points, -1) \
) / 2.0
return area
def is_simple(self) -> bool:
"""
Return whether the polygon has no self-intersections.
Returns
-------
bool
True if no two non-adjacent edges cross.
"""
return first_self_intersection(self.x_points,
self.y_points) is None
def split_self_intersections(self,
max_iter: int=500) -> list["Polygon2D"]:
"""
Split the polygon into simple (non-self-intersecting) loops.
Repeatedly splits each loop at its first self-crossing into the
loop between the two crossing edges and the remainder, both
closed through the crossing point, until no loop crosses itself.
A polygon that is already simple yields a single unchanged copy.
The loops are built with _new_from_loop, so they have the same
class and invert_polygon setting as this polygon.
Parameters
----------
max_iter : int, optional
Safety bound on the total number of splits. Loops still
waiting to be split when it is reached are returned as they
are, so they may not be simple, and a warning is issued.
Returns
-------
list[Polygon2D]
The loops, sorted by decreasing absolute area.
"""
pending = [(self.x_points.astype(float).copy(),
self.y_points.astype(float).copy())]
loops = []
n_splits = 0
n_unsplit = 0
while pending:
x, y = pending.pop()
hit = first_self_intersection(x, y)
if hit is None or n_splits >= max_iter:
loops.append(self._new_from_loop(x, y))
n_unsplit += hit is not None
continue
n_splits += 1
i, j, point = hit
# Loop between the two crossings, closed by the crossing point.
pending.append((np.concatenate([[point[0]], x[i + 1:j + 1]]),
np.concatenate([[point[1]], y[i + 1:j + 1]])))
# Everything else, rejoined through the crossing point.
pending.append((
np.concatenate([x[:i + 1], [point[0]], x[j + 1:]]),
np.concatenate([y[:i + 1], [point[1]], y[j + 1:]])))
if n_unsplit:
warnings.warn(
f"Reached max_iter={max_iter} in split_self_intersections, "
f"leaving {n_unsplit} loop(s) still self-intersecting.")
loops.sort(key=lambda loop: abs(loop.signed_area()), reverse=True)
return loops
def return_as_simple(self, keep: str="largest",
max_iter: int=500) -> "Polygon2D":
"""
Return a simple (non-self-intersecting) copy of the polygon.
Splits the polygon into simple loops with
split_self_intersections and keeps one of them. The default
keeps the largest-area loop, dropping small finger excursions
(for example the doubling-back of outer scrape-off layer flux
surfaces). A polygon that is already simple is returned as an
unchanged copy. The result has the same class as this polygon.
Parameters
----------
keep : {'largest', 'smallest'}, optional
Which loop to keep, by absolute area. Use
split_self_intersections to get all of them.
max_iter : int, optional
Safety bound on the number of splits.
Returns
-------
Polygon2D
A polygon with the self-intersection loops removed.
"""
if keep not in ("largest", "smallest"):
raise ValueError(
f"Expected keep to be 'largest' or 'smallest', got {keep!r}.")
loops = self.split_self_intersections(max_iter)
return loops[0] if keep == "largest" else loops[-1]
def _new_from_loop(self, x_points: np.ndarray,
y_points: np.ndarray) -> "Polygon2D":
"""
Build a polygon of the same kind from one loop of this polygon.
Used by split_self_intersections so that the loops keep the
settings of the polygon they were split from. Subclasses with
extra invariants override it to restore them on the loop.
Parameters
----------
x_points, y_points : numpy.ndarray
Vertices of the loop.
Returns
-------
Polygon2D
The loop, with invert_polygon copied from this polygon.
"""
return Polygon2D(x_points, y_points,
invert_polygon=self.invert_polygon)
def repeat_first_point(self):
"""
Conditionally repeat the first point at the end.
If the polygon is not already closed, adds the final edge to link the
first and last points. Useful for plotting.
"""
if not np.isclose(self.x_points[0], self.x_points[-1]) or \
not np.isclose(self.y_points[0], self.y_points[-1]):
self.x_points = np.append(self.x_points, self.x_points[0])
self.y_points = np.append(self.y_points, self.y_points[0])
self.clear_cached_properties()
def to_ordered(self, start_idx: int=0, k: int=20,
n_interp: int=0, resample_method: str="spline",
make_simple: bool=False):
"""
Return an OrderedPolygon2D version of this polygon.
Orders the points counterclockwise starting from the point with
index start_idx using a nearest-neighbor algorithm with up to k
nearest neighbors. See OrderedPolygon2D for the resample_method
and make_simple options. The invert_polygon setting is kept.
"""
from .ordered_polygon_2d_m import OrderedPolygon2D
return OrderedPolygon2D(self.x_points, self.y_points,
start_idx=start_idx, k=k,
n_interp=n_interp,
resample_method=resample_method,
make_simple=make_simple,
invert_polygon=self.invert_polygon)
[docs]
@classmethod
def from_2d_array(cls, data_slice: np.ndarray):
"""Create an instance from a single 2D array slice (npoints, xy)."""
if data_slice.shape[-1] != 2:
raise ValueError("The last dimension (dim_xy) must have length 2.")
x_points = data_slice[:, 0]
y_points = data_slice[:, 1]
return cls(x_points, y_points)
