Fix open_edges, now checks all. Refactor edge_intersect

.gitignore

@@ -5,3 +5,5 @@ __pycache__
 poly.png
 .cache
 .idea
+runtime
+/examples

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()
 

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)

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]
-'''

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)