BerGeo/h2/sidedness.py

import random
from collections import namedtuple
from enum import Enum, auto
import matplotlib.pyplot as plt

Point = namedtuple('Point', ['x', 'y'])


class Side(Enum):
    ON = auto()
    ABOVE = auto()
    BELOW = auto()


def sidedness(p1, p2, p3, eps=0.0000001):

    # Find line from p1 to p2, ask where p3 is in relation to this

    y = p3.y * (p2.x - p1.x)
    x = p3.x

    a = (p2.y - p1.y)
    b = p2.y * (p2.x - p1.x) - a * p2.x

    if y - eps < a * x + b < y + eps:
        return Side.ON
    elif y > a * x + b:
        return Side.ABOVE
    return Side.BELOW


# test


p1 = Point(4, 4)
p2 = Point(0, 0)
p3 = Point(5, 2)


# print(sidedness(p1, p2, p3))


def genPoint():
    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 graham_scan(points):

    # A funky issue where both a and b become negative in the sidedness test causes us to have to use
    # Side.ABOVE for both tests, regardless of UH or LH.

    sorted_points = sorted(points)

    UH = sorted_points[:2]
    #del sorted_points[0]

    for s in sorted_points[2:]:
        while len(UH) > 1 and (sidedness(UH[-2], UH[-1], s) != Side.ABOVE):
            del UH[-1]
        UH.append(s)

    reversed_list = list(reversed(sorted_points))
    reversed_list.append(UH[0])
    LH = reversed_list[:2]
    #del reversed_list[0]

    for s in reversed_list[2:]:
        while len(LH) > 1 and (sidedness(LH[-2], LH[-1], s) != Side.ABOVE):
            del LH[-1]
        LH.append(s)

    return UH, LH


p = [genPoint() for _ in range(30)]
UH, LH = graham_scan(p)

x = [point.x for point in p]
y = [point.y for point in p]

UH_x = [point.x for point in UH]
UH_y = [point.y for point in UH]

LH_x = [point.x for point in LH]
LH_y = [point.y for point in LH]

plt.plot(UH_x, UH_y, 'ro-')
plt.plot(LH_x, LH_y, 'ro-')

plt.scatter(x,y)
plt.show()