145 lines
4.8 KiB
Python
145 lines
4.8 KiB
Python
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
|