BerGeo/h2/quick_hull.py

103 lines
2.1 KiB
Python
Raw Normal View History

2018-09-20 12:50:16 +00:00
import random
from collections import namedtuple
2018-09-20 13:44:39 +00:00
from typing import Set
2018-09-20 12:50:16 +00:00
import matplotlib.pyplot as plt
from enum import Enum, auto
from math import atan2, degrees, tau, pi, acos, sqrt
Point = namedtuple('Point', 'x y')
class Side(Enum):
ON = auto()
ABOVE = auto()
BELOW = auto()
def display(points, hull):
x = [point.x for point in points]
y = [point.y for point in points]
h_x = [point.x for point in hull]
h_y = [point.y for point in hull]
plt.plot(h_x, h_y, 'ro-')
plt.scatter(x, y)
plt.show()
def gen_point():
a = random.uniform(1, 5)
b = random.uniform(1, 5)
x_i = random.uniform(1, 5)
p_i = Point(x_i, a * x_i + b)
return p_i
2018-09-20 13:44:39 +00:00
def distance(a, b, c):
try:
nom = abs((b.y - a.y) * c.x - (b.x - a.x) * c.y + b.x * a.y - b.y * a.x)
den = sqrt((b.y - a.y) ** 2 + (b.x - a.x) ** 2)
return nom / den
except ZeroDivisionError:
return 0
2018-09-20 12:50:16 +00:00
2018-09-20 13:44:39 +00:00
def is_left(a: Point, b: Point, c: Point):
return ((b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x)) > 0
2018-09-20 12:50:16 +00:00
2018-09-20 13:44:39 +00:00
def quick_hull(points: Set[Point]):
assert len(points) >= 2
2018-09-20 12:50:16 +00:00
2018-09-20 13:44:39 +00:00
left = min(points)
right = max(points)
2018-09-20 12:50:16 +00:00
2018-09-20 13:44:39 +00:00
hull = {left, right}
2018-09-20 13:58:02 +00:00
points = points - hull
2018-09-20 12:50:16 +00:00
2018-09-20 13:44:39 +00:00
find_hull({p for p in points if not is_left(left, right, p)},
left,
right,
hull)
2018-09-20 12:50:16 +00:00
2018-09-20 13:58:02 +00:00
find_hull({p for p in points if not is_left(right, left, p)},
2018-09-20 13:44:39 +00:00
right,
left,
hull)
2018-09-20 12:50:16 +00:00
2018-09-20 13:44:39 +00:00
return hull
2018-09-20 12:50:16 +00:00
2018-09-20 13:44:39 +00:00
def find_hull(points: Set[Point], p: Point, q: Point, hull: Set[Point]):
if not points:
return
2018-09-20 12:50:16 +00:00
2018-09-20 13:44:39 +00:00
farthest = max(points, key=lambda point: distance(p, q, point))
hull.add(farthest)
2018-09-20 12:50:16 +00:00
2018-09-20 13:44:39 +00:00
s1 = {point
for point in points
if not is_left(p, farthest, point)}
print("--")
print(s1)
2018-09-20 12:50:16 +00:00
2018-09-20 13:44:39 +00:00
s2 = {point
for point in points
if not is_left(farthest, q, point)}
print(s2)
2018-09-20 12:50:16 +00:00
2018-09-20 13:44:39 +00:00
find_hull(s1, p, farthest, hull)
find_hull(s2, farthest, q, hull)
2018-09-20 12:50:16 +00:00
2018-09-20 13:58:02 +00:00
peepees = {gen_point() for _ in range(11)}
hulliees = quick_hull(peepees)
2018-09-20 12:50:16 +00:00
2018-09-20 13:58:02 +00:00
display(peepees, hulliees)