path: root/stsci/sphere/graph.py

# -*- 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?
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)

        # TODO: Vectorize this
        for node in self._nodes:
            if np.all(np.abs(node._point - point) < 2 ** -32):
                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

        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(title=title)
            raise

    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 xrange(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.
                del edges[j]  # CD
                edges = new_edges + edges
                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 = 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.
        """
        polygons = []
        edges = set(self._edges)  # copy
        while len(edges):
            polygon = []
            # Carefully pick out an "original" edge first.  Synthetic
            # edges may not be pointing in the right direction to
            # properly calculate the area.
            for edge in edges:
                if len(edge._source_polygons) == 1:
                    break
            edges.remove(edge)
            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")