#################################################################################################### # # Patro - A Python library to make patterns for fashion design # Copyright (C) 2017 Fabrice Salvaire # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program. If not, see . # #################################################################################################### """Module to implement polygon. """ #################################################################################################### __all__ = ['Polygon2D'] #################################################################################################### import math import numpy as np from .Primitive import PrimitiveNP, ClosedPrimitiveMixin, PathMixin, Primitive2DMixin from .Segment import Segment2D from .Triangle import Triangle2D from Patro.Common.Math.Functions import sign #################################################################################################### # Fixme: PrimitiveNP last ??? class Polygon2D(Primitive2DMixin, ClosedPrimitiveMixin, PathMixin, PrimitiveNP): """Class to implements 2D Polygon.""" ############################################## # def __new__(cls, *points): # # remove consecutive duplicates # no_duplicate = [] # for point in points: # if no_duplicate and point == no_duplicate[-1]: # continue # no_duplicate.append(point) # if len(no_duplicate) > 1 and no_duplicate[-1] == no_duplicate[0]: # no_duplicate.pop() # last point was same as first # # remove collinear points # i = -3 # while i < len(no_duplicate) - 3 and len(no_duplicate) > 2: # a, b, c = no_duplicate[i], no_duplicate[i + 1], no_duplicate[i + 2] # if Point.is_collinear(a, b, c): # no_duplicate.pop(i + 1) # if a == c: # no_duplicate.pop(i) # else: # i += 1 # if len(vertices) > 3: # return GeometryEntity.__new__(cls, *vertices, **kwargs) # elif len(vertices) == 3: # return Triangle(*vertices, **kwargs) # elif len(vertices) == 2: # return Segment(*vertices, **kwargs) # else: # return Point(*vertices, **kwargs) ############################################## def __init__(self, *points): if len(points) < 3: raise ValueError('Polygon require at least 3 vertexes') PrimitiveNP.__init__(self, points) self._edges = None self._is_simple = None self._is_convex = None self._area = None # self._cross = None # self._barycenter = None # self._major_axis_angle = None self._major_axis = None # self._minor_axis = None # self._axis_ratio = None ############################################## @property def is_triangle(self): return self.number_of_points == 3 def to_triangle(self): if self.is_triangle: return Triangle2D(*self.points) else: raise ValueError('Polygon is not a triangle') ############################################## @property def edges(self): if self._edges is None: edges = [] N = self.number_of_points for i in range(N): j = (i+1) % N edge = Segment2D(self._points[i], self._points[j]) edges.append(edge) self._edges = edges return iter(self._edges) ############################################## def _test_is_simple(self): edges = list(self.edges) # intersections = [] # Test for edge intersection for edge1 in edges: for edge2 in edges: if edge1 != edge2: # Fixme: recompute line for edge intersection, intersect = edge1.intersection(edge2) if intersect: common_vertex = edge1.share_vertex_with(edge2) if common_vertex is not None: if common_vertex == intersection: continue else: # degenerated case where a vertex lie on an edge return False else: # two edge intersect # intersections.append(intersection) return False return True ############################################## def _test_is_convex(self): # https://en.wikipedia.org/wiki/Convex_polygon # http://mathworld.wolfram.com/ConvexPolygon.html if not self.is_simple: return False edges = list(self.edges) # a polygon is convex if all turns from one edge vector to the next have the same sense # sign = edges[-1].perp_dot(edges[0]) sign0 = sign(edges[-1].cross(edges[0])) for i in range(len(edges)): if sign(edges[i].cross(edges[i+1])) != sign0: return False return True ############################################## @property def is_simple(self): """Test if the polygon is simple, i.e. if it doesn't self-intersect.""" if self._is_simple is None: self._is_simple = self._test_is_simple() return self._is_simple ############################################## @property def is_convex(self): if self._is_convex is None: self._is_convex = self._test_is_convex() return self._is_convex @property def is_concave(self): return not self.is_convex ############################################## @property def perimeter(self): return sum([edge.length for edge in self.edges]) ############################################## @property def point_barycenter(self): center = self.start_point for point in self.iter_from_second_point(): center += point return center / self.number_of_points ############################################## def _compute_area_barycenter(self): r"""Compute polygon area and barycenter. Polygon area is determined by .. math:: \begin{align} \mathbf{A} &= \frac{1}{2} \sum_{i=0}^{n-1} P_i \otimes P_{i+1} \\ &= \frac{1}{2} \sum_{i=0}^{n-1} \begin{vmatrix} x_i & x_{i+1} \\ y_i & y_{i+1} \end{vmatrix} \\ &= \frac{1}{2} \sum_{i=0}^{n-1} x_i y_{i+1} - x_{i+1} y_i \end{align} where :math:`x_n = x_0` Polygon barycenter is determined by .. math:: \begin{align} \mathbf{C} &= \frac{1}{6\mathbf{A}} \sum_{i=0}^{n-1} (P_i + P_{i+1}) \times (P_i \otimes P_{i+1}) \\ &= \frac{1}{6\mathbf{A}} \sum_{i=0}^{n-1} \begin{pmatrix} (x_i + x_{i+1}) (x_i y_{i+1} - x_{i+1} y_i) \\ (y_i + y_{i+1}) (x_i y_{i+1} - x_{i+1} y_i) \end{pmatrix} \end{align} References * On the Calculation of Arbitrary Moments of Polygons, Carsten Steger, Technical Report FGBV–96–05, October 1996 * http://mathworld.wolfram.com/PolygonArea.html * https://en.wikipedia.org/wiki/Polygon#Area_and_centroid """ if not self.is_simple: return None # area = self._points[-1].cross(self._points[0]) # for i in range(self.number_of_points): # area *= self._points[i].cross(self._points[i+1]) # P0, P1, Pn-1, P0 points = self.closed_point_array # from 0 to n-1 : P0, ..., Pn-1 xi = points[0,:-1] yi = points[1,:-1] # from 1 to n : P1, ..., Pn-1, P0 xi1 = points[0,1:] yi1 = points[1,1:] # Fixme: np.cross ??? cross = xi * yi1 - xi1 * yi self._cross = cross area = .5 * np.sum(cross) if area == 0: # print('Null area') self._area = 0 self._barycenter = self.start_point else: factor = 1 / (6*area) x = factor * np.sum((xi + xi1) * cross) y = factor * np.sum((yi + yi1) * cross) # area of a convex polygon is defined to be positive if the points are arranged in a # counterclockwise order, and negative if they are in clockwise order (Beyer 1987). self._area = abs(area) self._barycenter = self.__vector_cls__(x, y) ############################################## def _compute_inertia_moment(self): r"""Compute inertia moment on vertices. .. warning:: untrusted formulae .. math:: \begin{align} I_x &= \frac{1}{12} \sum (y_i^2 + y_i y_{i+1} + y_{i+1}^2) (x_i y_{i+1} - x_{i+1} y_i) \\ I_y &= \frac{1}{12} \sum (x_i^2 + x_i x_{i+1} + x_{i+1}^2) (x_i y_{i+1} - x_{i+1} y_i) \\ I_{xy} &= \frac{1}{24} \sum (x_i y_{i+1} + 2 x_i y_i + 2 x_{i+1} y_{i+1} + x_{i+1} y_i) (x_i y_{i+1} - x_{i+1} y_i) \end{align} Reference * https://en.wikipedia.org/wiki/Second_moment_of_area#Any_cross_section_defined_as_polygon """ # self.recenter() # Fixme: duplicated code # P0, P1, Pn-1, P0 points = self.closed_point_array # from 0 to n-1 : P0, ..., Pn-1 xi = points[0,:-1] yi = points[1,:-1] # from 1 to n : P1, ..., Pn-1, P0 xi1 = points[0,1:] yi1 = points[1,1:] # computation on vertices number_of_points = self.number_of_points Ix = np.sum(yi**2) / number_of_points Iy = np.sum(xi**2) / number_of_points Ixy = - np.sum(xi*yi) / number_of_points # cross = xi * yi1 - xi1 * yi # cross = self._cross # Ix = 1/(12*self._area) * np.sum((yi**2 + yi*yi1 + yi1**2) * cross) # Iy = 1/(12*self._area) * np.sum((xi**2 + xi*xi1 + xi1**2) * cross) # Ixy = 1/(24*self._area) * np.sum((xi*yi1 + 2*(xi*yi + xi1*yi1) + xi1*yi) * cross) # cx, cy = self._barycenter # Ix -= cy**2 # Iy -= cx**2 # Ixy -= cx*cy # Ix = -Ix # Iy = -Iy # print(Ix, Iy, Ixy) if Ixy == 0: if Iy >= Ix: self._major_axis_angle = 0 lambda1 = Iy lambda2 = Ix vx = 0 v1y = 1 v2y = 0 else: self._major_axis_angle = 90 lambda1 = Ix lambda2 = Iy vx = 1 v1y = 0 v2y = 1 else: Is = Iy + Ix Id = Ix - Iy sqrt0 = math.sqrt(Id*Id + 4*Ixy*Ixy) lambda1 = (Is + sqrt0) / 2 lambda2 = (Is - sqrt0) / 2 vx = Ixy v1y = (Id + sqrt0) / 2 v2y = (Id - sqrt0) / 2 if lambda1 < lambda2: v1y, v2y = v2y, v1y lambda1, lambda2 = lambda2, lambda1 self._major_axis_angle = - math.degrees(math.atan(v1y/vx)) self._major_axis = 4 * math.sqrt(math.fabs(lambda1)) self._minor_axis = 4 * math.sqrt(math.fabs(lambda2)) if self._minor_axis != 0: self._axis_ratio = self._major_axis / self._minor_axis else: self._axis_ratio = 0 ############################################## def _check_area(self): if self.is_simple and self._area is None: self._compute_area_barycenter() ############################################## @property def area(self): """Return polygon area.""" self._check_area() return self._area ############################################## @property def barycenter(self): """Return polygon barycenter.""" self._check_area() return self._barycenter ############################################## def recenter(self): """Recenter the polygon to the barycenter.""" # if self._centred: # return barycenter = self._barycenter for point in self._points: point -= barycenter # self._centred = True ############################################## def _check_moment(self): if self.is_simple and self._major_axis is None: self._compute_inertia_moment() ############################################## @property def major_axis_angle(self): self._check_moment() return self._major_axis_angle @property def major_axis(self): self._check_moment() return self._major_axis @property def minor_axis(self): self._check_moment() return self._minor_axis @property def axis_ratio(self): self._check_moment() return self._axis_ratio ############################################## def _crossing_number_test(self, point): """Crossing number test for a point in a polygon.""" # Wm. Randolph Franklin, "PNPOLY - Point Inclusion in Polygon Test" Web Page (2000) # https://www.ecse.rpi.edu/Homepages/wrf/research/geom/pnpoly.html crossing_number = 0 x = point.x y = point.y for edge in self.edges: if ((edge.p0.y <= y < edge.p1.y) or # upward crossing (edge.p1.y <= y < edge.p0.y)): # downward crossing xi = edge.p0.x + (y - edge.p0.y) / edge.vector.slope if x < xi: crossing_number += 1 # Fixme: even/odd func return (crossing_number & 1) == 1 # odd => in ############################################## def _winding_number_test(self, point): """Winding number test for a point in a polygon.""" # more accurate than crossing number test # http://geomalgorithms.com/a03-_inclusion.html#wn_PnPoly() winding_number = 0 y = point.y for edge in self.edges: if edge.p0.y <= y: if edge.p1.y > y: # upward crossing if edge.is_left(point): winding_number += 1 else: if edge.p1.y <= y: # downward crossing if edge.is_right(point): winding_number -= 1 return winding_number > 0 ############################################## def is_point_inside(self, point): # http://geomalgorithms.com/a03-_inclusion.html # http://paulbourke.net/geometry/polygonmesh/#insidepoly # Fixme: bounding box test return self._winding_number_test(point)