BerGeo/h2/quick_hull.py

from math import sqrt
from typing import Set

from util import Point, gen_point, display


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)


points = {gen_point() for _ in range(30)}
hull = quick_hull(points)

display(points, hull)