lim added files to sphere branch

git-svn-id: http://svn.stsci.edu/svn/ssb/stsci_python/stsci_python/branches/sphere@15878 fe389314-cf27-0410-b35b-8c050e845b92

Former-commit-id: b51412b5b545598ae101152eeaed470b08616176

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