From 69037ec7714d8fa5dcc8853dd4254f7bcda397f3 Mon Sep 17 00:00:00 2001 From: Alexander Munch-Hansen Date: Thu, 20 Sep 2018 14:50:16 +0200 Subject: [PATCH] hmmmm --- h2/gift_wrapper.py | 28 ++++++++---- h2/quick_hull.py | 106 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 125 insertions(+), 9 deletions(-) create mode 100644 h2/quick_hull.py diff --git a/h2/gift_wrapper.py b/h2/gift_wrapper.py index b2a4e54..edc28dd 100644 --- a/h2/gift_wrapper.py +++ b/h2/gift_wrapper.py @@ -2,7 +2,7 @@ import random from collections import namedtuple import matplotlib.pyplot as plt -from math import atan2, degrees, tau, pi, acos, sqrt +from math import acos, sqrt Point = namedtuple('Point', 'x y') Vector = namedtuple('Vector', 'x y') @@ -18,8 +18,19 @@ def gen_point(): return p_i -def calc_angle(v: Vector, upside_down=False) -> float: - return (90 - degrees(atan2(v.y, v.x)) - 180*upside_down) % 360 +def calc_angle(v1: Vector, v2: Vector) -> float: + dot = (v1.x * v2.x) + (v1.y * v2.y) + len_1 = sqrt(v1.x**2 + v1.y**2) + len_2 = sqrt(v2.x**2 + v2.y**2) + + tmp = dot / (len_1 * len_2) + + # pls + hmmmmmmm = round(tmp, 6) + + return acos(hmmmmmmm) + + def calc_vector(p1: Point, p2: Point) -> Vector: @@ -41,16 +52,15 @@ def display(points, hull): def rapper(points: set): min_pt = min(points) - max_pt = max(points) hull = [min_pt] - upside_down = False - + comp_vec = Vector(0, 1) while True: hull.append(min(points - {hull[-1]}, - key=lambda p: calc_angle(calc_vector(hull[-1], p), upside_down))) - if hull[-1].x == max_pt.x: - upside_down = True + key=lambda p: calc_angle(comp_vec, calc_vector(hull[-1], p)))) + comp_vec = calc_vector(hull[-2], hull[-1]) + display(points, hull) + if hull[-1] == min_pt: break diff --git a/h2/quick_hull.py b/h2/quick_hull.py new file mode 100644 index 0000000..d84db6f --- /dev/null +++ b/h2/quick_hull.py @@ -0,0 +1,106 @@ +import random +from collections import namedtuple +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 sidedness(p1, p2, p3, eps=0.0000001): + y = p1.y * (p3.x - p2.x) + x = p1.x + a = (p3.y - p2.y) + b = p3.y * (p3.x - p2.x) - a * p3.x + + line = a*x+b + + if y - eps < line < y + eps: + return [Side.ON, y - line] + elif y > line: + return [Side.ABOVE, y - line] + return [Side.BELOW, y - line] + + +def prune_below_find_highest(points: set, line_segment, side = Side.BELOW, lel=max): + min_pt, max_pt = line_segment + + above_pts = [] + + for point in points: + sidedness_res = sidedness(min_pt, max_pt, point) + if sidedness_res[0] != side: + above_pts.append(point) + + highest_above = lel(above_pts, key=lambda pt: pt[1]) + + return above_pts, highest_above + + +def find_line_segment(points): + min_pt = min(points, key=lambda pt: pt.x) + max_pt = max(points, key=lambda pt: pt.x) + return min_pt, max_pt + + +def quick_hull(points, side=Side.BELOW, lel=max, hull=None): + if hull is None: + hull = set() + + print(len(points)) + if len(points) <= 2: + return hull + + min_pt, max_pt = find_line_segment(points) + + above_points, max_dist_pt = prune_below_find_highest(points, (min_pt, max_pt), side, lel) + + left_side = prune_below_find_highest(points, (min_pt, max_dist_pt), side, lel)[0] + right_side = prune_below_find_highest(points, (max_dist_pt, max_pt), side, lel)[0] + + hull.add(max_dist_pt) + + quick_hull(left_side, side, lel, hull=hull).union(quick_hull(right_side, side, lel, hull=hull)) + + return hull + +points = {Point(2, 5), Point(3, 1), Point(3, 4), Point(9, 4), Point(1, 5), Point(5, 7)} +#points = {gen_point() for _ in range(10)} + +line_segment = find_line_segment(points) +uh_hull = quick_hull(points) +lh_hull = quick_hull(points, Side.ABOVE, min) + +hull = uh_hull.union(lh_hull) + +display(points, hull)