BerGeo/h2/quick_hull.py
2018-09-20 16:54:47 +02:00

97 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):
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
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}
points = points - hull
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 not is_left(right, left, 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: abs(distance(p, q, point)))
hull.add(farthest)
points.remove(farthest)
find_hull({po for po in points if not is_left(p, farthest, po)},
p,
farthest,
hull)
find_hull({po for po in points if not is_left(farthest, q, po)},
farthest,
q,
hull)
peepees = {gen_point() for _ in range(30)}
hulliees = quick_hull(peepees)
display(peepees, hulliees)