BerGeo/h2/quick_hull.py

import random
from collections import namedtuple
from enum import Enum, auto
from math import sqrt
from typing import Set

import matplotlib.pyplot as plt

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)


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

display(points, hull)