BerGeo/h2/quick_hull.py

from math import sqrt
from typing import Set

from profile import Profiler
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


@Profiler("quick_hull")
def quick_hull(points: Set[Point]):
    left = min(points)
    right = max(points)

    hull = {left, right}
    points = points - hull

    with Profiler("partitioning set"):
        partition = {p for p in points if not is_left(left, right, p)}
    find_hull(partition,
              left,
              right,
              hull)

    with Profiler("partitioning set"):
        partition = {p for p in points if not is_left(right, left, p)}
    find_hull(partition,
              right,
              left,
              hull)

    return hull


def find_hull(points: Set[Point], p: Point, q: Point, hull: Set[Point]):
    if not points:
        return

    with Profiler("finding farthest point from line"):
        farthest = max(points, key=lambda point: abs(distance(p, q, point)))
        hull.add(farthest)
        points.remove(farthest)

    with Profiler("partitioning set"):
        partition = {po for po in points if not is_left(p, farthest, po)}
    find_hull(partition,
              p,
              farthest,
              hull)

    with Profiler("partitioning set"):
        partition = {po for po in points if not is_left(farthest, q, po)}
    find_hull(partition,
              farthest,
              q,
              hull)


if __name__ == '__main__':
    points = {gen_point() for _ in range(30)}
    hull = quick_hull(points)

    display(points, hull)