# -*- 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. """ from __future__ import absolute_import, division, unicode_literals, print_function # STDLIB import itertools import weakref # THIRD-PARTY import numpy as np # LOCAL from . import great_circle_arc from . import vector # Set to True to enable some sanity checks DEBUG = True 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 `~Graph.Node` represents a single point, connected by an arbitrary number of `~Graph.Edge` objects to other `~Graph.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 : `~Graph.Edge` instance The edge to follow away from. Returns ------- other : `~Graph.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, thres=2e-8): """ Returns `True` if the location of this and the *other* `~Graph.Node` are the same. Parameters ---------- other : `~Graph.Node` instance The other node. thres : float If difference is smaller than this, points are equal. The default value of 2e-8 radians is set based on empirical test cases. Relative threshold based on the actual sizes of polygons is not implemented. """ return np.array_equal(self._point, other._point) class Edge: """ An `~Graph.Edge` represents a connection between exactly two `~Graph.Node` objects. This `~Graph.Edge` class has no direction. """ def __init__(self, A, B, source_polygons=[]): """ Parameters ---------- A, B : `~Graph.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 : `~Graph.Node` instance Returns ------- other : `~Graph.Node` instance """ nodes = self._nodes try: return nodes[not nodes.index(node)] except IndexError: raise ValueError("Following from disconnected node") def equals(self, other): """ Returns `True` if the other edge is between the same two nodes. Parameters ---------- other : `~Graph.Edge` instance Returns ------- equals : bool """ if (self._nodes[0].equals(other._nodes[0]) and self._nodes[1].equals(other._nodes[1])): return True if (self._nodes[1].equals(other._nodes[0]) and self._nodes[0].equals(other._nodes[1])): return True return False 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: return self._source_polygons.add(polygon) start_node = nodeA = self.add_node(points[0], [polygon]) for i in range(1, len(points) - 1): nodeB = self.add_node(points[i], [polygon]) # Don't create self-pointing edges if nodeB is not nodeA: 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 : `~Graph.Node` instance The new node """ # Any nodes whose Cartesian coordinates are closer together # than 2 ** -32 will cause numerical problems in the # intersection calculations, so we merge any nodes that # are closer together than that. # Don't add nodes that already exist. Update the existing # node's source_polygons list to include the new polygon. point = vector.normalize_vector(point) if len(self._nodes): nodes = list(self._nodes) node_array = np.array([node._point for node in nodes]) diff = np.all(np.abs(point - node_array) < 2 ** -32, axis=-1) indices = np.nonzero(diff)[0] if len(indices): node = nodes[indices[0]] node._source_polygons.update(source_polygons) return node new_node = self.Node(point, source_polygons) self._nodes.add(new_node) return new_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 : `~Graph.Node` instance """ assert node in self._nodes for edge in list(node._edges): nodeB = edge.follow(node) nodeB._edges.remove(edge) if len(nodeB._edges) == 0: self._nodes.remove(nodeB) self._edges.remove(edge) if node in self._nodes: 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 : `~Graph.Node` instances source_polygons : `~sphere.polygon.SphericalPolygon` instance, optional The polygon(s) this edge came from. Used for bookkeeping. Returns ------- edge : `~Graph.Edge` instance The new edge """ assert A in self._nodes assert B in self._nodes # Don't add any edges that already exist. Update the edge's # source polygons list to include the new polygon. Care needs # to be taken here to not create an Edge until we know we need # one, otherwise the Edge will get hooked up to the nodes but # be orphaned. for edge in self._edges: if ((A is edge._nodes[0] and B is edge._nodes[1]) or (A is edge._nodes[1] and B is edge._nodes[0])): edge._source_polygons.update(source_polygons) return edge new_edge = self.Edge(A, B, source_polygons) self._edges.add(new_edge) return new_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 : `~Graph.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 `~Graph.Edge` *edge* at `~Graph.Node` *node*, removing *edge* and replacing it with two new `~Graph.Edge` instances. It is intended that *E* is along the original edge, but that is not enforced. Parameters ---------- edge : `~Graph.Edge` instance The edge to split node : `~Graph.Node` instance The node to insert Returns ------- edgeA, edgeB : `~Graph.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) if edge not in (edgeA, edgeB): self.remove_edge(edge) return [edgeA, edgeB] def _sanity_check(self, title, 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 log = open("/tmp/sanity.log", 'a') log.write("\n" + title + "\n") try: unique_edges = set() log.write("list of edges\n") for edge in self._edges: log.write("edge: " + str(edge) + "\n") for node in edge._nodes: log.write("node: " + str(node) + "\n") 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) log.write("list of nodes\n") for node in self._nodes: log.write("node: " + str(node) + "\n") if node_is_2: assert len(node._edges) % 2 == 0 else: assert len(node._edges) >= 2 for edge in node._edges: log.write("edge: " + str(edge) + "\n") assert node in edge._nodes assert edge in self._edges except AssertionError as e: log.close() import traceback traceback.print_exc() self._dump_graph(title=title) raise log.close() def _dump_graph(self, title=None, lon_0=0, lat_0=90, projection='vandg', func=lambda x: len(x._edges), highlight_edge=None, intersections=[]): from mpl_toolkits.basemap import Basemap from matplotlib import pyplot as plt fig = plt.figure() m = Basemap() minx = np.inf miny = np.inf maxx = -np.inf maxy = -np.inf for polygon in self._source_polygons: polygon.draw(m, lw=10, alpha=0.1, color="black") v = polygon._points ra, dec = vector.vector_to_radec(v[:, 0], v[:, 1], v[:, 2]) x, y = m(ra, dec) for x0 in x: minx = min(x0, minx) maxx = max(x0, maxx) for y0 in y: miny = min(y0, miny) maxy = max(y0, maxy) counts = {} for node in self._nodes: count = func(node) counts.setdefault(count, []) counts[count].append(list(node._point)) 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]) if edge is highlight_edge: color = 'red' lw = 2 else: color = 'blue' lw = 0.5 m.drawgreatcircle(r0, d0, r1, d1, color=color, lw=lw) 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 x0 in x: minx = min(x0, minx) maxx = max(x0, maxx) for y0 in y: miny = min(y0, miny) maxy = max(y0, maxy) for v in intersections: if np.all(np.isfinite(v)): ra, dec = vector.vector_to_radec(v[0], v[1], v[2]) x, y = m(ra, dec) m.plot(x, y, 'x', markersize=20) plt.xlim(minx, maxx) plt.ylim(miny, maxy) if title: plt.title("%s, %d v, %d e" % ( title, len(self._nodes), len(self._edges))) 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("union - remove cut lines") self._find_all_intersections() self._sanity_check("union - find all intersections") self._remove_interior_edges() self._sanity_check("union - remove interior edges") self._remove_3ary_edges() self._sanity_check("union - remove 3ary edges") self._remove_orphaned_nodes() self._sanity_check("union - remove orphan nodes", 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("intersection - remove cut lines") self._find_all_intersections() self._sanity_check("intersection - find all intersections") self._remove_exterior_edges() self._sanity_check("intersection - remove exterior edges") self._remove_cut_lines() self._sanity_check("intersection - remove cut lines [2]") self._remove_orphaned_nodes() self._sanity_check("intersection - remove orphan nodes", 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 range(len(edges)): AB = edges[i] if AB not in self._edges: continue A, B = AB._nodes if len(A._edges) == 3 and len(B._edges) == 3: self.remove_edge(AB) def _get_edge_points(self, edges): return (np.array([x._nodes[0]._point for x in edges]), np.array([x._nodes[1]._point for x in edges])) def _find_point_to_arc_intersections(self): # For speed, we want to vectorize all of the intersection # calculations. Therefore, there is a list of edges, and an # array of points for all of the nodes. Then calculating the # intersection between an edge and all other nodes becomes a # fast, vectorized operation. edges = list(self._edges) starts, ends = self._get_edge_points(edges) nodes = list(self._nodes) nodes_array = np.array([x._point for x in nodes]) # Split all edges by any nodes that intersect them while len(edges) > 1: AB = edges.pop(0) A, B = list(AB._nodes) intersects = great_circle_arc.intersects_point( A._point, B._point, nodes_array) intersection_indices = np.nonzero(intersects)[0] for index in intersection_indices: node = nodes[index] if node not in AB._nodes: newA, newB = self.split_edge(AB, node) new_edges = [ edge for edge in (newA, newB) if edge not in edges] for end_point in AB._nodes: node._source_polygons.update( end_point._source_polygons) edges = new_edges + edges break def _find_arc_to_arc_intersections(self): # 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 = self._get_edge_points(edges) # Calculate edge-to-edge intersections and break # edges on the intersection point. 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: 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 = [ edge for edge in (newA, newB, newC, newD) if edge not in edges] # Delete CD, and push the new edges to the # front so they will be tested for intersection # against all remaining edges. edges = new_edges + edges[:j] + edges[j+1:] new_starts, new_ends = self._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 _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. """ self._find_point_to_arc_intersections() self._find_arc_to_arc_intersections() def _remove_interior_edges(self): """ Removes any nodes that are contained inside other polygons. What's left is the (possibly disjunct) outline. """ polygons = self._source_polygons for edge in self._edges: edge._count = 0 A, B = edge._nodes for polygon in polygons: if (not polygon in edge._source_polygons and ((polygon in A._source_polygons or polygon.contains_point(A._point)) and (polygon in B._source_polygons or polygon.contains_point(B._point))) and polygon.contains_point( great_circle_arc.midpoint(A._point, B._point))): edge._count += 1 for edge in list(self._edges): if edge._count >= 1: self.remove_edge(edge) self._remove_orphaned_nodes() 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 A, B = edge._nodes for polygon in polygons: if polygon in edge._source_polygons: edge._count += 1 elif ((polygon in A._source_polygons or polygon.contains_point(A._point)) and (polygon in B._source_polygons or polygon.contains_point(B._point)) and polygon.contains_point( great_circle_arc.midpoint(A._point, B._point))): edge._count += 1 for edge in list(self._edges): if edge._count < len(polygons): self.remove_edge(edge) self._remove_orphaned_nodes() def _remove_3ary_edges(self, large_first=False): """ Remove edges between pairs of nodes that have 3 or more 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 = list(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 _remove_orphaned_nodes(self): """ Remove nodes with fewer than 2 edges. """ while True: removes = [] for node in list(self._nodes): if len(node._edges) < 2: removes.append(node) if len(removes): for node in removes: if node in self._nodes: self.remove_node(node) else: break def _trace(self): """ Given a graph that has had cutlines removed and all intersections found, traces it to find a resulting single polygon. """ from . import polygon as mpolygon polygons = [] edges = set(self._edges) # copy for edge in self._edges: edge._followed = False while len(edges): points = [] edge = edges.pop() start_node = node = edge._nodes[0] while True: # TODO: Do we need this if clause any more? if len(points): if not np.array_equal(points[-1], node._point): points.append(node._point) else: points.append(node._point) for edge in node._edges: if edge._followed is False: break else: raise ValueError("No more edges to follow") edge._followed = True edges.discard(edge) node = edge.follow(node) if node is start_node: points.append(node._point) break polygon = mpolygon._SingleSphericalPolygon(points) area = polygon.area() reverse_polygon = mpolygon._SingleSphericalPolygon(points[::-1]) reverse_area = reverse_polygon.area() if reverse_area < area: polygon = reverse_polygon polygons.append(polygon) return mpolygon.SphericalPolygon(polygons)