diff --git a/.gitignore b/.gitignore index e3a3035..07e6b81 100644 --- a/.gitignore +++ b/.gitignore @@ -5,3 +5,5 @@ __pycache__ poly.png .cache .idea +runtime +/examples diff --git a/graph.py b/graph.py index ced0e91..5c79c0b 100644 --- a/graph.py +++ b/graph.py @@ -43,6 +43,7 @@ def get_adjacent(self, point): def __contains__(self, point): return point in self.points + # TODO: This runs slow, self.points should be a set in __init__ def __eq__(self, edge): return set(self.points) == set(edge.points) @@ -86,6 +87,7 @@ def __init__(self, polygons): def get_adjacent_points(self, point): return [edge.get_adjacent(point) for edge in self.graph[point]] + ''' should return a generator/iterator ''' def get_points(self): return self.graph.keys() diff --git a/test_country.py b/test_country.py new file mode 100644 index 0000000..d3a61fe --- /dev/null +++ b/test_country.py @@ -0,0 +1,42 @@ +import cProfile +import pstats +import shapefile +import matplotlib.pyplot as plt +import matplotlib.patches as patches +from matplotlib.patches import Polygon +from matplotlib.collections import PatchCollection +from graph import Graph, Point, Edge +from visible_vertices import (edge_intersect, point_edge_distance, + visible_vertices, angle) + + +def test_country(): + sf = shapefile.Reader("examples/GSHHS_c_L1.dbf") + shapes = sf.shapes() + + poly1 = [Point(a[0], a[1]) for a in shapes[4].points] + graph = Graph([poly1]) + ship = Point(-80, -30) + + fig = plt.figure(1, figsize=(10, 10), dpi=90) + ax = fig.add_subplot(111) + + x = [a.x for a in poly1] + y = [a.y for a in poly1] + ax.plot(x, y, color='gray', alpha=0.7, linewidth=1, + solid_capstyle='round', zorder=2) + + for pp in graph.get_points(): + visible = visible_vertices(pp, graph, ship, None) + for point in visible: + x = [pp.x, point.x] + y = [pp.y, point.y] + ax.plot(x, y, color='red', alpha=0.7, linewidth=1) + + ax.set_title("Python visibility graph") + fig.savefig("poly.png", bbox_inches='tight') + +cProfile.run('test_country()', 'runtime') +p = pstats.Stats('runtime') +p.sort_stats('cumulative').print_stats(10) +p.sort_stats('time').print_stats(10) diff --git a/tests.py b/tests.py index a4545e8..e4ce501 100644 --- a/tests.py +++ b/tests.py @@ -168,159 +168,3 @@ def test_visible_vertices_3(self): ship = Point(2.5, 3.5) visible = visible_vertices(ship, graph, ship, None) assert visible == [point_c, point_b, point_d, point_f] - -''' -class TestShortestPaths: - - def setup_method(self, method): - # test obstacle A - self.point_a = Point(1.0, 2.0) - self.point_b = Point(4.0, 2.5) - self.point_c = Point(5.0, 3.0) - self.point_d = Point(4.5, 5.0) - self.point_e = Point(3.0, 5.0) - self.point_f = Point(3.0, 7.0) - self.point_g = Point(4.0, 7.0) - self.point_h = Point(4.0, 8.0) - self.point_i = Point(1.0, 8.0) - self.points = [self.point_a, self.point_b, self.point_c, self.point_d, - self.point_e, self.point_f, self.point_g, self.point_h, - self.point_i] - - self.edges = [] - self.edges.append(Edge(self.point_a, self.point_b)) - self.edges.append(Edge(self.point_a, self.point_i)) - self.edges.append(Edge(self.point_b, self.point_a)) - self.edges.append(Edge(self.point_b, self.point_c)) - self.edges.append(Edge(self.point_c, self.point_b)) - self.edges.append(Edge(self.point_c, self.point_d)) - self.edges.append(Edge(self.point_d, self.point_c)) - self.edges.append(Edge(self.point_d, self.point_e)) - self.edges.append(Edge(self.point_e, self.point_d)) - self.edges.append(Edge(self.point_e, self.point_f)) - self.edges.append(Edge(self.point_f, self.point_e)) - self.edges.append(Edge(self.point_f, self.point_g)) - self.edges.append(Edge(self.point_g, self.point_f)) - self.edges.append(Edge(self.point_g, self.point_h)) - self.edges.append(Edge(self.point_h, self.point_g)) - self.edges.append(Edge(self.point_h, self.point_i)) - self.edges.append(Edge(self.point_i, self.point_a)) - self.edges.append(Edge(self.point_i, self.point_h)) - - self.polygon_a = Polygon(self.points, self.edges) - - # test obstacle B - self.point_a = Point(6.0, 1.0) - self.point_b = Point(7.0, 2.0) - self.point_c = Point(8.0, 6.0) - self.point_d = Point(10.0, 1.5) - self.points = [self.point_a, self.point_b, self.point_c, self.point_d] - - self.edges = [] - self.edges.append(Edge(self.point_a, self.point_b)) - self.edges.append(Edge(self.point_a, self.point_d)) - self.edges.append(Edge(self.point_b, self.point_a)) - self.edges.append(Edge(self.point_b, self.point_c)) - self.edges.append(Edge(self.point_c, self.point_b)) - self.edges.append(Edge(self.point_c, self.point_d)) - self.edges.append(Edge(self.point_d, self.point_a)) - self.edges.append(Edge(self.point_d, self.point_c)) - - self.polygon_b = Polygon(self.points, self.edges) - - # test operating network - self.point_a = Point(1.0, 2.0) - self.point_b = Point(4.0, 2.5) - self.point_c = Point(5.0, 3.0) - self.point_d = Point(4.5, 5.0) - self.point_e = Point(4.0, 8.0) - self.point_f = Point(1.0, 8.0) - self.point_g = Point(6.0, 1.0) - self.point_h = Point(8.0, 6.0) - self.point_i = Point(10.0, 1.5) - self.point_j = Point(4.0, 7.0) - self.points = [self.point_a, self.point_b, self.point_c, self.point_d, - self.point_e, self.point_f, self.point_g, self.point_h, - self.point_i, self.point_j] - - self.edges = [] - self.edges.append(Edge(self.point_a, self.point_b)) - self.edges.append(Edge(self.point_a, self.point_f)) - self.edges.append(Edge(self.point_a, self.point_g)) - self.edges.append(Edge(self.point_b, self.point_a)) - self.edges.append(Edge(self.point_b, self.point_c)) - self.edges.append(Edge(self.point_c, self.point_b)) - self.edges.append(Edge(self.point_c, self.point_g)) - self.edges.append(Edge(self.point_c, self.point_h)) - self.edges.append(Edge(self.point_c, self.point_d)) - self.edges.append(Edge(self.point_d, self.point_c)) - self.edges.append(Edge(self.point_d, self.point_e)) - self.edges.append(Edge(self.point_d, self.point_g)) - self.edges.append(Edge(self.point_e, self.point_d)) - self.edges.append(Edge(self.point_e, self.point_f)) - self.edges.append(Edge(self.point_e, self.point_j)) - self.edges.append(Edge(self.point_e, self.point_h)) - self.edges.append(Edge(self.point_f, self.point_a)) - self.edges.append(Edge(self.point_f, self.point_e)) - self.edges.append(Edge(self.point_g, self.point_a)) - self.edges.append(Edge(self.point_g, self.point_c)) - self.edges.append(Edge(self.point_g, self.point_d)) - self.edges.append(Edge(self.point_g, self.point_h)) - self.edges.append(Edge(self.point_g, self.point_i)) - self.edges.append(Edge(self.point_h, self.point_c)) - self.edges.append(Edge(self.point_h, self.point_e)) - self.edges.append(Edge(self.point_h, self.point_g)) - self.edges.append(Edge(self.point_h, self.point_i)) - self.edges.append(Edge(self.point_i, self.point_h)) - self.edges.append(Edge(self.point_i, self.point_g)) - self.edges.append(Edge(self.point_j, self.point_e)) - - self.op_network = Polygon(self.points, self.edges) - self.graph = Graph([self.polygon_a, self.polygon_b]) - self.op_graph = Graph([self.op_network]) - - def test_shortest_path_1(self): - ship = Point(10.0, 5.5) - port = Point(1.0, 2.0) - shortest = shortest_path(self.op_graph, ship, port) - edge_a = Edge(ship, Point(8.0, 6.0)) - edge_b = Edge(Point(8.0, 6.0), Point(5.0, 3.0)) - edge_c = Edge(Point(5.0, 3.0), Point(4.0, 2.5)) - edge_d = Edge(Point(4.0, 2.5), Point(1.0, 2.0)) - - # Draw the obstacles and operating network - fig = plt.figure(1, figsize=(5, 5), dpi=90) - ax = fig.add_subplot(111) - - # Draw the obstacles - for polygon in self.graph.polygons: - for edge in polygon.edges: - x = [edge.points[0].x, edge.points[1].x] - y = [edge.points[0].y, edge.points[1].y] - ax.plot(x, y, color='gray', alpha=0.7, linewidth=5, - solid_capstyle='round', zorder=2) - - # Draw the operating network - for edge in self.op_graph.get_edges(): - x = [edge.points[0].x, edge.points[1].x] - y = [edge.points[0].y, edge.points[1].y] - ax.plot(x, y, color='red', alpha=0.7, linewidth=1) - - # draw shortest path - for edge in shortest: - x = [edge.points[0].x, edge.points[1].x] - y = [edge.points[0].y, edge.points[1].y] - ax.plot(x, y, color='green', alpha=0.7, linewidth=2) - - ax.set_title("Python visibility graph and shortest path") - xrange = [0, 11] - yrange = [0, 9] - ax.set_xlim(*xrange) - ax.set_xticks(range(*xrange) + [xrange[-1]]) - ax.set_ylim(*yrange) - ax.set_yticks(range(*yrange) + [yrange[-1]]) - ax.set_aspect(1) - fig.savefig("poly.png", bbox_inches='tight') - - assert shortest == [edge_a, edge_b, edge_c, edge_d] -''' diff --git a/visible_vertices.py b/visible_vertices.py index 11ffe9f..3f05095 100644 --- a/visible_vertices.py +++ b/visible_vertices.py @@ -30,9 +30,16 @@ def visible_vertices(point, graph, ship, port): open_edges.remove(edge) except ValueError: pass - is_visible = False - if not open_edges or edge_distance(point, p) <= point_edge_distance(point, p, open_edges[0]): + ''' Not ideal, checking each open_edge, but should be able to add new + open edges in sorted order, so you only check the closest edge. + Is likely due to concave polygons.''' + smallest_edge = float('inf') + for e in open_edges: + dist = point_edge_distance(point, p, e) + if dist < smallest_edge: + smallest_edge = dist + if not open_edges or edge_distance(point, p) <= smallest_edge: if previous_point is not None and angle(point, p) == angle(point, previous_point): if edge_distance(point, p) < edge_distance(point, previous_point): is_visible = True @@ -44,6 +51,7 @@ def visible_vertices(point, graph, ship, port): if is_visible: visible.append(p) + ''' This should only be needed if I fix open_edges order ''' edge_order = [] for edge in graph[p]: if (point not in edge) and counterclockwise(point, edge, p): @@ -64,6 +72,7 @@ def angle2(point_a, point_b, point_c): def counterclockwise(point, edge, endpoint): + ''' TODO: merge with ccw() ''' edge_point1, edge_point2 = edge.points if edge_point1 == endpoint: angle_diff = angle(point, edge_point2) - angle(point, endpoint) @@ -112,49 +121,33 @@ def angle(center, point): """ dx = point.x - center.x dy = point.y - center.y - if dx == 0: if dy < 0: return pi*3 / 2 - else: - return pi / 2 + return pi / 2 if dy == 0: if dx < 0: return pi - else: - return 0 - + return 0 if dx < 0: return pi + atan(dy/dx) if dy < 0: return 2*pi + atan(dy/dx) - return atan(dy/dx) -def edge_intersect(point1, point2, edge): - """ - Return True if 'edge' is interesected by the line going - through 'point1' and 'point2', False otherwise. - If edge contains either 'point1' or 'point2', return False. - """ - edge_point1, edge_point2 = edge.points - if point1 in edge or point2 in edge: - return False - - if point1.x == point2.x: - x1_left = edge_point1.x < point1.x - x2_left = edge_point2.x < point1.x - return not (x1_left == x2_left) - - slope = (point1.y - point2.y) / (point1.x - point2.x) +def ccw(A, B, C): + return (C.y-A.y) * (B.x-A.x) > (B.y-A.y) * (C.x-A.x) - y1_ex = slope * (edge_point1.x - point1.x) + point1.y - y2_ex = slope * (edge_point2.x - point1.x) + point1.y - if y1_ex == edge_point1.y or y2_ex == edge_point2.y: +def edge_intersect(A, B, edge): + """ + Return True if 'edge' is interesected by the line going through 'A' and + 'B', False otherwise. If edge contains either 'A' or 'B', return False. + TODO: May be an issue with colinerity here. See: + http://www.geeksforgeeks.org/check-if-two-given-line-segments-intersect/ + """ + C, D = edge.points + if A in edge or B in edge: return False - - y1_below = (y1_ex > edge_point1.y) - y2_below = (y2_ex > edge_point2.y) - return not (y1_below == y2_below) + return ccw(A, C, D) != ccw(B, C, D) and ccw(A, B, C) != ccw(A, B, D)