BerGeo/h2/quick_hull.py
2018-09-20 15:44:39 +02:00

102 lines
2.0 KiB
Python

import random
from collections import namedtuple
from typing import Set
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
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
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
def quick_hull(points: Set[Point]):
assert len(points) >= 2
left = min(points)
right = max(points)
hull = {left, right}
find_hull({p for p in points if not is_left(left, right, p)},
left,
right,
hull)
find_hull({p for p in points if is_left(left, right, p)},
right,
left,
hull)
return hull
def find_hull(points: Set[Point], p: Point, q: Point, hull: Set[Point]):
if not points:
return
farthest = max(points, key=lambda point: distance(p, q, point))
hull.add(farthest)
s1 = {point
for point in points
if not is_left(p, farthest, point)}
print("--")
print(s1)
s2 = {point
for point in points
if not is_left(farthest, q, point)}
print(s2)
find_hull(s1, p, farthest, hull)
find_hull(s2, farthest, q, hull)
points = {gen_point() for _ in range(10)}
hull = quick_hull(points)
display(points, hull)