diff options
Diffstat (limited to 'lib/graph.py')
| -rw-r--r-- | lib/graph.py | 792 |
1 files changed, 792 insertions, 0 deletions
diff --git a/lib/graph.py b/lib/graph.py new file mode 100644 index 0000000..e9f03ba --- /dev/null +++ b/lib/graph.py @@ -0,0 +1,792 @@ +# -*- coding: utf-8 -*- + +# Copyright (C) 2011 Association of Universities for Research in +# Astronomy (AURA) +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions +# are met: +# +# 1. Redistributions of source code must retain the above +# copyright notice, this list of conditions and the following +# disclaimer. +# +# 2. Redistributions in binary form must reproduce the above +# copyright notice, this list of conditions and the following +# disclaimer in the documentation and/or other materials +# provided with the distribution. +# +# 3. The name of AURA and its representatives may not be used to +# endorse or promote products derived from this software without +# specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY AURA ``AS IS'' AND ANY EXPRESS OR +# IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL AURA BE LIABLE FOR ANY DIRECT, +# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, +# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED +# OF THE POSSIBILITY OF SUCH DAMAGE. + +""" +This contains the code that does the actual unioning of regions. +""" +# TODO: Weak references for memory management problems? + +# STDLIB +import itertools +import weakref + +# THIRD-PARTY +import numpy as np + +# LOCAL +#from . import great_circle_arc +#from . import vector +import great_circle_arc +import vector + +# Set to True to enable some sanity checks +DEBUG = False + + +class Graph: + """ + A graph of nodes connected by edges. The graph is used to build + unions between polygons. + + .. note:: + This class is not meant to be used directly. Instead, use + `sphere.polygon.SphericalPolygon.union` and + `sphere.polygon.SphericalPolygon.intersection`. + """ + + class Node: + """ + A `Node` represents a single point, connected by an arbitrary + number of `Edge` objects to other `Node` objects. + """ + def __init__(self, point, source_polygons=[]): + """ + Parameters + ---------- + point : 3-sequence (*x*, *y*, *z*) coordinate + + source_polygon : `~sphere.polygon.SphericalPolygon` instance, optional + The polygon(s) this node came from. Used for bookkeeping. + """ + self._point = np.asanyarray(point) + self._source_polygons = set(source_polygons) + self._edges = weakref.WeakSet() + + def __repr__(self): + return "Node(%s %d)" % (str(self._point), len(self._edges)) + + def follow(self, edge): + """ + Follows from one edge to another across this node. + + Parameters + ---------- + edge : `Edge` instance + The edge to follow away from. + + Returns + ------- + other : `Edge` instance + The other edge. + """ + edges = list(self._edges) + assert len(edges) == 2 + try: + return edges[not edges.index(edge)] + except IndexError: + raise ValueError("Following from disconnected edge") + + def equals(self, other): + """ + Returns `True` if the location of this and the *other* + `Node` are the same. + """ + return np.array_equal(self._point, other._point) + + + class Edge: + """ + An `Edge` represents a connection between exactly two `Node` + objects. This `Edge` class has no direction. + """ + def __init__(self, A, B, source_polygons=[]): + """ + Parameters + ---------- + A, B : `Node` instances + + source_polygon : `~sphere.polygon.SphericalPolygon` instance, optional + The polygon this edge came from. Used for bookkeeping. + """ + self._nodes = [A, B] + for node in self._nodes: + node._edges.add(self) + self._source_polygons = set(source_polygons) + + def __repr__(self): + nodes = self._nodes + return "Edge(%s -> %s)" % (nodes[0]._point, nodes[1]._point) + + def follow(self, node): + """ + Follow along the edge from the given *node* to the other + node. + + Parameters + ---------- + node : `Node` instance + + Returns + ------- + other : `Node` instance + """ + nodes = self._nodes + try: + return nodes[not nodes.index(node)] + except IndexError: + raise ValueError("Following from disconnected node") + + + def __init__(self, polygons): + """ + Parameters + ---------- + polygons : sequence of `~sphere.polygon.SphericalPolygon` instances + Build a graph from this initial set of polygons. + """ + self._nodes = set() + self._edges = set() + self._source_polygons = set() + + self.add_polygons(polygons) + + def add_polygons(self, polygons): + """ + Add more polygons to the graph. + + .. note:: + Must be called before `union` or `intersection`. + + Parameters + ---------- + polygons : sequence of `~sphere.polygon.SphericalPolygon` instances + Set of polygons to add to the graph + """ + for polygon in polygons: + self.add_polygon(polygon) + + def add_polygon(self, polygon): + """ + Add a single polygon to the graph. + + .. note:: + Must be called before `union` or `intersection`. + + Parameters + ---------- + polygon : `~sphere.polygon.SphericalPolygon` instance + Polygon to add to the graph + """ + points = polygon._points + + if len(points) < 3: + raise ValueError("Too few points in polygon") + + self._source_polygons.add(polygon) + + start_node = nodeA = self.add_node(points[0], [polygon]) + for i in range(1, len(points) - 1): + # Don't create self-pointing edges + if not np.array_equal(points[i], nodeA._point): + nodeB = self.add_node(points[i], [polygon]) + self.add_edge(nodeA, nodeB, [polygon]) + nodeA = nodeB + # Close the polygon + self.add_edge(nodeA, start_node, [polygon]) + + def add_node(self, point, source_polygons=[]): + """ + Add a node to the graph. It will be disconnected until used + in a call to `add_edge`. + + Parameters + ---------- + point : 3-sequence (*x*, *y*, *z*) coordinate + + source_polygon : `~sphere.polygon.SphericalPolygon` instance, optional + The polygon this node came from. Used for bookkeeping. + + Returns + ------- + node : `Node` instance + The new node + """ + node = self.Node(point, source_polygons) + self._nodes.add(node) + return node + + def remove_node(self, node): + """ + Removes a node and all of the edges that touch it. + + .. note:: + It is assumed that *Node* is already a part of the graph. + + Parameters + ---------- + node : `Node` instance + """ + assert node in self._nodes + + for edge in list(node._edges): + nodeB = edge.follow(node) + nodeB._edges.remove(edge) + self._edges.remove(edge) + self._nodes.remove(node) + + def add_edge(self, A, B, source_polygons=[]): + """ + Add an edge between two nodes. + + .. note:: + It is assumed both nodes already belong to the graph. + + Parameters + ---------- + A, B : `Node` instances + + source_polygons : `~sphere.polygon.SphericalPolygon` instance, optional + The polygon(s) this edge came from. Used for bookkeeping. + + Returns + ------- + edge : `Edge` instance + The new edge + """ + assert A in self._nodes + assert B in self._nodes + + edge = self.Edge(A, B, source_polygons) + self._edges.add(edge) + return edge + + def remove_edge(self, edge): + """ + Remove an edge from the graph. The nodes it points to remain intact. + + .. note:: + It is assumed that *edge* is already a part of the graph. + + Parameters + ---------- + edge : `Edge` instance + """ + assert edge in self._edges + + A, B = edge._nodes + A._edges.remove(edge) + if len(A._edges) == 0: + self.remove_node(A) + if A is not B: + B._edges.remove(edge) + if len(B._edges) == 0: + self.remove_node(B) + self._edges.remove(edge) + + def split_edge(self, edge, node): + """ + Splits an `Edge` *edge* at `Node` *node*, removing *edge* and + replacing it with two new `Edge` instances. It is intended + that *E* is along the original edge, but that is not enforced. + + Parameters + ---------- + edge : `Edge` instance + The edge to split + + node : `Node` instance + The node to insert + + Returns + ------- + edgeA, edgeB : `Edge` instances + The two new edges on either side of *node*. + """ + assert edge in self._edges + assert node in self._nodes + + A, B = edge._nodes + edgeA = self.add_edge(A, node, edge._source_polygons) + edgeB = self.add_edge(node, B, edge._source_polygons) + self.remove_edge(edge) + return [edgeA, edgeB] + + def _sanity_check(self, node_is_2=False): + """ + For debugging purposes: assert that edges and nodes are + connected to each other correctly and there are no orphaned + edges or nodes. + """ + if not DEBUG: + return + + try: + unique_edges = set() + for edge in self._edges: + for node in edge._nodes: + assert edge in node._edges + assert node in self._nodes + edge_repr = [tuple(x._point) for x in edge._nodes] + edge_repr.sort() + edge_repr = tuple(edge_repr) + # assert edge_repr not in unique_edges + unique_edges.add(edge_repr) + + for node in self._nodes: + if node_is_2: + assert len(node._edges) == 2 + else: + assert len(node._edges) >= 2 + for edge in node._edges: + assert node in edge._nodes + assert edge in self._edges + except AssertionError as e: + import traceback + traceback.print_exc() + self._dump_graph() + raise + + def _dump_graph(self, title=None, lon_0=0, lat_0=90, + projection='ortho', func=lambda x: len(x._edges)): + from mpl_toolkits.basemap import Basemap + from matplotlib import pyplot as plt + fig = plt.figure() + m = Basemap(projection="ortho", + lon_0=lon_0, + lat_0=lat_0) + m.drawparallels(np.arange(-90., 90., 20.)) + m.drawmeridians(np.arange(0., 420., 20.)) + m.drawmapboundary(fill_color='white') + + counts = {} + for node in self._nodes: + count = func(node) + counts.setdefault(count, []) + counts[count].append(list(node._point)) + + for k, v in counts.items(): + v = np.array(v) + ra, dec = vector.vector_to_radec(v[:, 0], v[:, 1], v[:, 2]) + x, y = m(ra, dec) + m.plot(x, y, 'o', label=str(k)) + + for edge in list(self._edges): + A, B = [x._point for x in edge._nodes] + r0, d0 = vector.vector_to_radec(A[0], A[1], A[2]) + r1, d1 = vector.vector_to_radec(B[0], B[1], B[2]) + m.drawgreatcircle(r0, d0, r1, d1, color='blue') + + if title: + plt.title(title) + plt.legend() + plt.show() + + def union(self): + """ + Once all of the polygons have been added to the graph, + join the polygons together. + + Returns + ------- + points : Nx3 array of (*x*, *y*, *z*) points + This is a list of points outlining the union of the + polygons that were given to the constructor. If the + original polygons are disjunct or contain holes, cut lines + will be included in the output. + """ + self._remove_cut_lines() + self._sanity_check() + self._find_all_intersections() + self._sanity_check() + self._remove_interior_nodes() + self._sanity_check() + self._remove_3ary_edges() + self._sanity_check(True) + return self._trace() + + def intersection(self): + """ + Once all of the polygons have been added to the graph, + calculate the intersection. + + Returns + ------- + points : Nx3 array of (*x*, *y*, *z*) points + This is a list of points outlining the intersection of the + polygons that were given to the constructor. If the + resulting polygons are disjunct or contain holes, cut lines + will be included in the output. + """ + self._remove_cut_lines() + self._sanity_check() + self._find_all_intersections() + self._sanity_check() + self._remove_exterior_edges() + self._sanity_check() + self._remove_3ary_edges(large_first=True) + self._sanity_check(True) + return self._trace() + + def _remove_cut_lines(self): + """ + Removes any cutlines that may already have existed in the + input polygons. This is so any cutlines in the final result + will be optimized to be as short as possible and won't + intersect each other. + + This works by finding coincident edges that are reverse to + each other, and then splicing around them. + """ + # As this proceeds, edges are removed from the graph. It + # iterates over a static list of all edges that exist at the + # start, so each time one is selected, we need to ensure it + # still exists as part of the graph. + + # This transforms the following (where = is the cut line) + # + # \ / + # A' + + B' + # | | + # A +====================+ B + # + # D +====================+ C + # | | + # D' + + C' + # / \ + # + # to this: + # + # \ / + # A' + + B' + # | | + # A + + C + # | | + # D' + + C' + # / \ + # + + edges = list(self._edges) + for i in xrange(len(edges)): + AB = edges[i] + if AB not in self._edges: + continue + A, B = AB._nodes + for j in xrange(i + 1, len(edges)): + CD = edges[j] + if CD not in self._edges: + continue + C, D = CD._nodes + # To be a cutline, the candidate edges need to run in + # the opposite direction, hence A == D and B == C, not + # A == C and B == D. + if (A.equals(D) and B.equals(C)): + # Create new edges A -> D' and C -> B' + self.add_edge( + A, D.follow(CD).follow(D), + AB._source_polygons | CD._source_polygons) + self.add_edge( + C, B.follow(AB).follow(B), + AB._source_polygons | CD._source_polygons) + + # Remove B and D which are identical to C and A + # respectively. We do not need to remove AB and + # CD because this will remove it for us. + self.remove_node(D) + self.remove_node(B) + + break + + def _find_all_intersections(self): + """ + Find all the intersecting edges in the graph. For each + intersecting pair, four new edges are created around the + intersection point. + """ + def get_edge_points(edges): + return (np.array([x._nodes[0]._point for x in edges]), + np.array([x._nodes[1]._point for x in edges])) + + # For speed, we want to vectorize all of the intersection + # calculations. Therefore, there is a list of edges, and two + # arrays containing the end points of those edges. They all + # need to have things added and removed from them at the same + # time to keep them in sync, but of course the interface for + # doing so is different between Python lists and numpy arrays. + + edges = list(self._edges) + starts, ends = get_edge_points(edges) + + while len(edges) > 1: + AB = edges.pop(0) + A = starts[0]; starts = starts[1:] # numpy equiv of "pop(0)" + B = ends[0]; ends = ends[1:] # numpy equiv of "pop(0)" + + # Calculate the intersection points between AB and all + # other remaining edges + with np.errstate(invalid='ignore'): + intersections = great_circle_arc.intersection( + A, B, starts, ends) + # intersects is `True` everywhere intersections has an + # actual intersection + intersects = np.isfinite(intersections[..., 0]) + + intersection_indices = np.nonzero(intersects)[0] + + # Iterate through the candidate intersections, if any -- + # we want to eliminate intersections that only intersect + # at the end points + for j in intersection_indices: + C = starts[j] + D = ends[j] + if (np.array_equal(C, A) or np.array_equal(C, B) or + np.array_equal(D, A) or np.array_equal(D, B)): + continue + + CD = edges[j] + E = intersections[j] + + # This is a bona-fide intersection, and E is the + # point at which the two lines intersect. Make a + # new node for it -- this must belong to the all + # of the source polygons of both of the edges that + # crossed. + + # A + # | + # C--E--D + # | + # B + + E = self.add_node( + E, AB._source_polygons | CD._source_polygons) + newA, newB = self.split_edge(AB, E) + newC, newD = self.split_edge(CD, E) + + new_edges = [newA, newB, newC, newD] + + # Delete CD, and push the new edges to the + # front so they will be tested for intersection + # against all remaining edges. + del edges[j] # CD + edges = new_edges + edges + new_starts, new_ends = get_edge_points(new_edges) + starts = np.vstack( + (new_starts, starts[:j], starts[j+1:])) + ends = np.vstack( + (new_ends, ends[:j], ends[j+1:])) + break + + def _remove_interior_nodes(self): + """ + Removes any nodes that are contained inside other polygons. + What's left is the (possibly disjunct) outline. + """ + polygons = self._source_polygons + + # node._count is the number of polygons that the node is + # inside, other than the polygon that it came from + for node in self._nodes: + node._count = 0 + for polygon in polygons: + if polygon in node._source_polygons: + continue + if polygon.contains_point(node._point): + node._count += 1 + + for node in list(self._nodes): + # Nodes with a count >= 1 are completely contained in the + # outer polygon, so they should be removed. What's left + # are only outer nodes. + if node._count >= 1: + self.remove_node(node) + + def _remove_exterior_edges(self): + """ + Removes any edges that are not contained in all of the source + polygons. What's left is the (possibly disjunct) outline. + """ + polygons = self._source_polygons + + for edge in self._edges: + edge._count = 0 + for polygon in polygons: + if (polygon in edge._source_polygons or + polygon.intersects_arc( + edge._nodes[0]._point, edge._nodes[1]._point)): + edge._count += 1 + + for edge in list(self._edges): + if edge._count < len(polygons): + self.remove_edge(edge) + + def _remove_3ary_edges(self, large_first=False): + """ + Remove edges between pairs of nodes that have more than 3 + edges. This removes triangles that can't be traced. + """ + if large_first: + max_ary = 0 + for node in self._nodes: + max_ary = max(len(node._edges), max_ary) + order = range(max_ary + 1, 2, -1) + else: + order = [3] + + for i in order: + removals = [] + for edge in list(self._edges): + if (len(edge._nodes[0]._edges) >= i and + len(edge._nodes[1]._edges) >= i): + removals.append(edge) + + for edge in removals: + if edge in self._edges: + self.remove_edge(edge) + + def _trace(self): + """ + Given a graph that has had cutlines removed and all + intersections found, traces it to find a resulting single + polygon. + """ + polygons = [] + edges = set(self._edges) # copy + seen_nodes = set() + while len(edges): + polygon = [] + edge = edges.pop() + start_node = node = edge._nodes[0] + while True: + # TODO: Do we need this if clause any more? + if len(polygon): + if not np.array_equal(polygon[-1], node._point): + polygon.append(node._point) + else: + polygon.append(node._point) + edge = node.follow(edge) + edges.discard(edge) + node = edge.follow(node) + if node is start_node: + polygon.append(node._point) + break + + polygons.append(np.asarray(polygon)) + + if len(polygons) == 1: + return polygons[0] + elif len(polygons) == 0: + return [] + else: + return self._join(polygons) + + def _join(self, polygons): + """ + If the graph is disjunct, joins the parts with cutlines. + + The closest nodes between each pair that don't intersect + any other edges are used as cutlines. + + TODO: This is not optimal, because the closest distance + between two polygons may not in fact be between two vertices, + but may be somewhere along an edge. + """ + def do_join(polygons): + all_polygons = polygons[:] + + skipped = 0 + + polyA = polygons.pop(0) + while len(polygons): + polyB = polygons.pop(0) + + # If fewer than 3 edges, it's not a polygon, + # just throw it out + if len(polyB) < 4: + continue + + # Find the closest set of vertices between polyA and + # polyB that don't cross any of the edges in *any* of + # the polygons + closest = np.inf + closest_pair_idx = (None, None) + for a in xrange(len(polyA) - 1): + A = polyA[a] + distances = great_circle_arc.length(A, polyB[:-1]) + b = np.argmin(distances) + distance = distances[b] + if distance < closest: + B = polyB[b] + # Does this candidate line cross other edges? + crosses = False + for poly in all_polygons: + if np.any( + great_circle_arc.intersects( + A, B, poly[:-1], poly[1:])): + crosses = True + break + if not crosses: + closest = distance + closest_pair_idx = (a, b) + + if not np.isfinite(closest): + # We didn't find a pair of points that don't cross + # something else, so we want to try to join another + # polygon. Defer the current polygon until later. + if len(polygons) in (0, skipped): + return None + polygons.append(polyB) + skipped += 1 + else: + # Splice the two polygons together using a cut + # line + a, b = closest_pair_idx + new_poly = np.vstack(( + # polyA up to and including the cut point + polyA[:a+1], + + # polyB starting with the cut point and + # wrapping around back to the cut point. + # Ignore the last point in polyB, because it + # is the same as the first + np.roll(polyB[:-1], -b, 0), + + # The cut point on polyB + polyB[b:b+1], + + # the rest of polyA, starting with the cut + # point + polyA[a:] + )) + + skipped = 0 + polyA = new_poly + + return polyA + + for permutation in itertools.permutations(polygons): + poly = do_join(list(permutation)) + if poly is not None: + return poly + + raise RuntimeError("Could not find cut points") |