This commit is contained in:
Casper 2018-10-11 15:54:52 +02:00
parent 9fee0dac99
commit c6bfb7d30e
No known key found for this signature in database
GPG Key ID: B1156723DB3BDDA8
4 changed files with 41 additions and 36 deletions

View File

@ -1,7 +1,7 @@
# Use atan2 instead of acos to calc angle; atan2(x,y) of the point we potentially want to add # Use atan2 instead of acos to calc angle; atan2(x,y) of the point we potentially want to add
from math import acos, sqrt from math import acos, sqrt
from util import Vector, Point, gen_point, display from util import Vector, Point, gen_point
def calc_angle(v1: Vector, v2: Vector) -> float: def calc_angle(v1: Vector, v2: Vector) -> float:
@ -10,8 +10,7 @@ def calc_angle(v1: Vector, v2: Vector) -> float:
len_2 = sqrt(v2.x**2 + v2.y**2) len_2 = sqrt(v2.x**2 + v2.y**2)
tmp = dot / (len_1 * len_2) tmp = dot / (len_1 * len_2)
return acos(max(min(tmp, 1), -1)) # acos is only defined in [-1,1]
return acos(round(tmp, 6))
def calc_vector(p1: Point, p2: Point) -> Vector: def calc_vector(p1: Point, p2: Point) -> Vector:

View File

@ -1,13 +1,17 @@
import random
import statistics import statistics
from math import inf from math import inf
from typing import Set, List, Tuple from typing import Set, List, Tuple
from scipy.optimize import linprog
import util
from util import Side, Point, gen_point, display from util import Side, Point, gen_point, display
def sidedness(slope: float, intersection: float, p3: Point, linprog_flipper: callable, eps=0.0000001) -> Side: def sidedness(slope: float, intersection: float, p3: Point, flipper: callable, eps=0.0000001) -> Side:
# finds where a point is in regards to a line # finds where a point is in regards to a line
if linprog_flipper(p3.y) - eps <= linprog_flipper(slope * p3.x + intersection) <= linprog_flipper(p3.y) + eps: if flipper(p3.y) - eps <= flipper(slope * p3.x + intersection) <= flipper(p3.y) + eps:
return Side.ON return Side.ON
elif p3.y > slope * p3.x + intersection: elif p3.y > slope * p3.x + intersection:
return Side.ABOVE return Side.ABOVE
@ -45,17 +49,16 @@ def solve_2dlp(c: Tuple[float, float], constraints: List[Tuple[Tuple[float, floa
x1, x2 = (-inf, -inf) if c1 > 0 else (inf, inf) x1, x2 = (-inf, -inf) if c1 > 0 else (inf, inf)
#random.shuffle(constraints) #random.shuffle(constraints)
for i, ((a1, a2), b) in enumerate(constraints[1:], start=1): for i, ((a1, a2), b) in enumerate(constraints[1:], start=0):
if a1*x1 + a2*x2 <= b: if not a1*x1 + a2*x2 <= b:
continue x1 = solve_1dlp(c1 - c2*a1/a2,
[(ai1 - ai2*a1 / a2, bi - ai2*b / a2) for (ai1, ai2), bi in constraints[:i]])
x1 = solve_1dlp(c1 - c2*a1/a2, [(ai1 - ai2*a1 / a2, bi - ai2*b / a2) for (ai1, ai2), bi in constraints[:i]])
x2 = (b - a1*x1) / a2 x2 = (b - a1*x1) / a2
return x1, x2 return x1, x2
def mbc_ch(points: Set[Point], linprog_flipper: callable) -> Set[Point]: def mbc_ch(points: Set[Point], flipper: callable) -> Set[Point]:
if len(points) <= 2: if len(points) <= 2:
return points return points
@ -67,26 +70,19 @@ def mbc_ch(points: Set[Point], linprog_flipper: callable) -> Set[Point]:
pr = {p for p in points if p.x >= med_x} pr = {p for p in points if p.x >= med_x}
# Find the bridge over the vertical line in pm # Find the bridge over the vertical line in pm
c = (linprog_flipper(med_x), linprog_flipper(1)) slope, intercept = solve_2dlp((flipper(med_x), flipper(1)),
A = [(linprog_flipper(-p.x), linprog_flipper(-1)) for p in points] [((flipper(-p.x), flipper(-1)), flipper(-p.y)) for p in points])
b = [linprog_flipper(-p.y) for p in points]
#result = linprog(c, A, b, bounds=[(None, None), (None, None)], options={"tol": 0.00001})
#slope, intercept = result.x[0], result.x[1]
slope, intercept = solve_2dlp(c, list(zip(A, b)))
#display_line_only(points, slope, intercept, [])
# Find the two points which are on the line, should work # Find the two points which are on the line, should work
left_point = next(p for p in pl if sidedness(slope, intercept, p, linprog_flipper) == Side.ON) on = {p for p in points if sidedness(slope, intercept, p, flipper) == Side.ON}
right_point = next(p for p in pr if sidedness(slope, intercept, p, linprog_flipper) == Side.ON) left_point = min(on)
right_point = max(on)
# Prune the points between the two line points # Prune the points between the two line points
pl = {p for p in pl if p.x <= left_point.x} pl = {p for p in points if p.x <= left_point.x}
pr = {p for p in pr if p.x >= right_point.x} pr = {p for p in points if p.x >= right_point.x}
return set.union(mbc_ch(pl, linprog_flipper), {left_point, right_point}, mbc_ch(pr, linprog_flipper)) return set.union(mbc_ch(pl, flipper), {left_point, right_point}, mbc_ch(pr, flipper))
def mbc(points: Set[Point]) -> Set[Point]: def mbc(points: Set[Point]) -> Set[Point]:

View File

@ -7,16 +7,26 @@ from graham import graham_scan
from mbc import mbc from mbc import mbc
from quick_hull import quick_hull from quick_hull import quick_hull
random.seed(1337_420) #random.seed(1337_420)
points = {util.gen_point(-100, 100) for i in range(5_000)} while True:
print("="*100)
points = {util.gen_point(-100, 100) for i in range(100)}
print("Points:", sorted(points))
# Sanity check # Sanity check
graham = set(graham_scan(points)) graham = set(graham_scan(points))
gift = set(rapper(points)) gift = set(rapper(points))
quick = quick_hull(points) quick = quick_hull(points)
mbch = set.union(mbc(points)) mbch = set.union(mbc(points))
assert gift == graham == quick == mbch
print("graham:", sorted(graham))
print("gift :", sorted(gift))
print("quick :", sorted(quick))
print("mbch :", sorted(mbch))
if not gift == graham == quick == mbch:
util.display(points=points, hull=mbch)
def time_it(f: callable, args: tuple = (), iterations=20): def time_it(f: callable, args: tuple = (), iterations=20):

View File

@ -32,7 +32,7 @@ def display(points: Set[Point], hull: Set[Point]):
plt.show() plt.show()
def display_line_only(points: Set[Point], slope: int, intercept: int, line_points: Set[Point]): def display_line_only(points: Set[Point], slope: float, intercept: float, line_points: Set[Point]):
x = [point.x for point in points] x = [point.x for point in points]
y = [point.y for point in points] y = [point.y for point in points]