import math import sys from collections import namedtuple Point = namedtuple('Point', ['x', 'y']) EPSILON = math.sqrt(sys.float_info.epsilon) def earclip(polygon): """ Simple earclipping algorithm for a given polygon p. polygon is expected to be an array of 2-tuples of the cartesian points of the polygon For a polygon with n points it will return n-2 triangles. The triangles are returned as an array of 3-tuples where each item in the tuple is a 2-tuple of the cartesian point. e.g >>> polygon = [(0,1), (-1, 0), (0, -1), (1, 0)] >>> triangles = tripy.earclip(polygon) >>> triangles [((1, 0), (0, 1), (-1, 0)), ((1, 0), (-1, 0), (0, -1))] Implementation Reference: - https://www.geometrictools.com/Documentation/TriangulationByEarClipping.pdf """ ear_vertex = [] triangles = [] polygon = [Point(*point) for point in polygon] if _is_clockwise(polygon): polygon.reverse() point_count = len(polygon) for i in range(point_count): prev_index = i - 1 prev_point = polygon[prev_index] point = polygon[i] next_index = (i + 1) % point_count next_point = polygon[next_index] if _is_ear(prev_point, point, next_point, polygon): ear_vertex.append(point) while ear_vertex and point_count >= 3: ear = ear_vertex.pop(0) i = polygon.index(ear) prev_index = i - 1 prev_point = polygon[prev_index] next_index = (i + 1) % point_count next_point = polygon[next_index] polygon.remove(ear) point_count -= 1 triangles.append(((prev_point.x, prev_point.y), (ear.x, ear.y), (next_point.x, next_point.y))) if point_count > 3: prev_prev_point = polygon[prev_index - 1] next_next_index = (i + 1) % point_count next_next_point = polygon[next_next_index] groups = [ (prev_prev_point, prev_point, next_point, polygon), (prev_point, next_point, next_next_point, polygon), ] for group in groups: p = group[1] if _is_ear(*group): if p not in ear_vertex: ear_vertex.append(p) elif p in ear_vertex: ear_vertex.remove(p) return triangles def _is_clockwise(polygon): s = 0 polygon_count = len(polygon) for i in range(polygon_count): point = polygon[i] point2 = polygon[(i + 1) % polygon_count] s += (point2.x - point.x) * (point2.y + point.y) return s > 0 def _is_convex(prev, point, next): return _triangle_sum(prev.x, prev.y, point.x, point.y, next.x, next.y) < 0 def _is_ear(p1, p2, p3, polygon): ear = _contains_no_points(p1, p2, p3, polygon) and \ _is_convex(p1, p2, p3) and \ _triangle_area(p1.x, p1.y, p2.x, p2.y, p3.x, p3.y) > 0 return ear def _contains_no_points(p1, p2, p3, polygon): for pn in polygon: if pn in (p1, p2, p3): continue elif _is_point_inside(pn, p1, p2, p3): return False return True def _is_point_inside(p, a, b, c): area = _triangle_area(a.x, a.y, b.x, b.y, c.x, c.y) area1 = _triangle_area(p.x, p.y, b.x, b.y, c.x, c.y) area2 = _triangle_area(p.x, p.y, a.x, a.y, c.x, c.y) area3 = _triangle_area(p.x, p.y, a.x, a.y, b.x, b.y) areadiff = abs(area - sum([area1, area2, area3])) < EPSILON return areadiff def _triangle_area(x1, y1, x2, y2, x3, y3): return abs((x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)) / 2.0) def _triangle_sum(x1, y1, x2, y2, x3, y3): return x1 * (y3 - y2) + x2 * (y1 - y3) + x3 * (y2 - y1) def calculate_total_area(triangles): result = [] for triangle in triangles: sides = [] for i in range(3): next_index = (i + 1) % 3 pt = triangle[i] pt2 = triangle[next_index] # Distance between two points side = math.sqrt(math.pow(pt2[0] - pt[0], 2) + math.pow(pt2[1] - pt[1], 2)) sides.append(side) # Heron's numerically stable forumla for area of a triangle: # https://en.wikipedia.org/wiki/Heron%27s_formula # However, for line-like triangles of zero area this formula can produce an infinitesimally negative value # as an input to sqrt() due to the cumulative arithmetic errors inherent to floating point calculations: # https://people.eecs.berkeley.edu/~wkahan/Triangle.pdf # For this purpose, abs() is used as a reasonable guard against this condition. c, b, a = sorted(sides) area = .25 * math.sqrt(abs((a + (b + c)) * (c - (a - b)) * (c + (a - b)) * (a + (b - c)))) result.append((area, a, b, c)) triangle_area = sum(tri[0] for tri in result) return triangle_area