diff options
Diffstat (limited to 'stsci')
| -rw-r--r-- | stsci/__init__.py | 8 | ||||
| -rw-r--r-- | stsci/sphere/__init__.py | 5 | ||||
| -rw-r--r-- | stsci/sphere/graph.py | 909 | ||||
| -rw-r--r-- | stsci/sphere/great_circle_arc.py | 368 | ||||
| -rw-r--r-- | stsci/sphere/polygon.py | 831 | ||||
| -rw-r--r-- | stsci/sphere/test/__init__.py | 1 | ||||
| -rw-r--r-- | stsci/sphere/test/benchmarks.py | 39 | ||||
| -rw-r--r-- | stsci/sphere/test/data/1904-66_TAN.fits.gz | bin | 0 -> 113196 bytes | |||
| -rw-r--r-- | stsci/sphere/test/data/2chipA.fits.gz | bin | 0 -> 227579 bytes | |||
| -rw-r--r-- | stsci/sphere/test/data/2chipB.fits.gz | bin | 0 -> 227493 bytes | |||
| -rw-r--r-- | stsci/sphere/test/data/2chipC.fits.gz | bin | 0 -> 183594 bytes | |||
| -rw-r--r-- | stsci/sphere/test/test.py | 277 | ||||
| -rw-r--r-- | stsci/sphere/test/test_intersection.py | 189 | ||||
| -rw-r--r-- | stsci/sphere/test/test_shared.py | 12 | ||||
| -rw-r--r-- | stsci/sphere/test/test_skyline.py | 225 | ||||
| -rw-r--r-- | stsci/sphere/test/test_union.py | 215 | ||||
| -rw-r--r-- | stsci/sphere/test/test_util.py | 14 | ||||
| -rw-r--r-- | stsci/sphere/vector.py | 243 |
18 files changed, 3336 insertions, 0 deletions
diff --git a/stsci/__init__.py b/stsci/__init__.py new file mode 100644 index 0000000..a5133cc --- /dev/null +++ b/stsci/__init__.py @@ -0,0 +1,8 @@ +# This is a special __init__.py required to place sample_package under the +# stsci namespace package. There should be no other code in this module. +try: + from pkg_resources import declare_namespace + declare_namespace(__name__) +except ImportError: + from pkgutil import extend_path + __path__ = extend_path(__path__, __name__) diff --git a/stsci/sphere/__init__.py b/stsci/sphere/__init__.py new file mode 100644 index 0000000..6a03273 --- /dev/null +++ b/stsci/sphere/__init__.py @@ -0,0 +1,5 @@ +import sys + +if sys.version_info[0] >= 3: + # Python 3 compatibility + __builtins__['xrange'] = range diff --git a/stsci/sphere/graph.py b/stsci/sphere/graph.py new file mode 100644 index 0000000..d8ccf75 --- /dev/null +++ b/stsci/sphere/graph.py @@ -0,0 +1,909 @@ +# -*- 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. + """ + point = vector.normalize_vector(*point) + + 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 + emphirical test cases. Relative threshold based on + the actual sizes of polygons is not implemented. + """ + # return np.array_equal(self._point, other._point) + return great_circle_arc.length(self._point, other._point, + degrees=False) < thres + + + 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._start_node = None + + 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]) + if self._start_node is None: + self._start_node = start_node + 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 + """ + new_node = self.Node(point, source_polygons) + + # Don't add nodes that already exist. Update the existing + # node's source_polygons list to include the new polygon. + for node in self._nodes: + if node.equals(new_node): + node._source_polygons.update(source_polygons) + return node + + 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) + 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 : `~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.equals(edge._nodes[0]) and + B.equals(edge._nodes[1])) or + (B.equals(edge._nodes[0]) and + A.equals(edge._nodes[1]))): + 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)): + from mpl_toolkits.basemap import Basemap + from matplotlib import pyplot as plt + fig = plt.figure() + m = Basemap() + + counts = {} + for node in self._nodes: + count = func(node) + counts.setdefault(count, []) + counts[count].append(list(node._point)) + + minx = np.inf + miny = np.inf + maxx = -np.inf + maxy = -np.inf + 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 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') + + 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_3ary_edges(large_first=True) + self._sanity_check("intersection - remove 3ary edges") + 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 + 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) + + # First, split all edges by any nodes that intersect them + 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)" + + distance = great_circle_arc.length(A, B) + for node in self._nodes: + if node not in AB._nodes: + distanceA = great_circle_arc.length(node._point, A) + distanceB = great_circle_arc.length(node._point, B) + if np.abs((distanceA + distanceB) - distance) < 1e-8: + newA, newB = self.split_edge(AB, node) + + new_edges = [ + edge for edge in (newA, newB) + if edge not in edges] + + edges = new_edges + edges + new_starts, new_ends = get_edge_points(new_edges) + starts = np.vstack( + (new_starts, starts)) + ends = np.vstack( + (new_ends, ends)) + break + + edges = list(self._edges) + starts, ends = get_edge_points(edges) + + # Next, 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: + C = starts[j] + D = ends[j] + 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 = 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_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 + for polygon in polygons: + if (not polygon in edge._source_polygons and + polygon.intersects_arc( + edge._nodes[0]._point, edge._nodes[1]._point)): + edge._count += 1 + + for edge in list(self._edges): + if edge._count >= 1: + self.remove_edge(edge) + + 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 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: + 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 + seen_nodes = set() + 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") diff --git a/stsci/sphere/great_circle_arc.py b/stsci/sphere/great_circle_arc.py new file mode 100644 index 0000000..d5e78f2 --- /dev/null +++ b/stsci/sphere/great_circle_arc.py @@ -0,0 +1,368 @@ +# -*- 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. + +""" +The `sphere.great_circle_arc` module contains functions for computing +the length, intersection, angle and midpoint of great circle arcs. + +Great circles are circles on the unit sphere whose center is +coincident with the center of the sphere. Great circle arcs are the +section of those circles between two points on the unit sphere. +""" + +from __future__ import with_statement, division, absolute_import, unicode_literals + +# THIRD-PARTY +import numpy as np + + +try: + from . import math_util + HAS_C_UFUNCS = True +except ImportError: + HAS_C_UFUNCS = False + +if HAS_C_UFUNCS: + inner1d = math_util.inner1d +else: + from numpy.core.umath_tests import inner1d + + +__all__ = ['angle', 'intersection', 'intersects', 'length', 'midpoint', + 'interpolate'] + + +def _fast_cross(a, b): + """ + This is a reimplementation of `numpy.cross` that only does 3D x + 3D, and is therefore faster since it doesn't need any + conditionals. + """ + if HAS_C_UFUNCS: + return math_util.cross(a, b) + + cp = np.empty(np.broadcast(a, b).shape) + aT = a.T + bT = b.T + cpT = cp.T + + cpT[0] = aT[1]*bT[2] - aT[2]*bT[1] + cpT[1] = aT[2]*bT[0] - aT[0]*bT[2] + cpT[2] = aT[0]*bT[1] - aT[1]*bT[0] + + return cp + +def _cross_and_normalize(A, B): + if HAS_C_UFUNCS: + return math_util.cross_and_norm(A, B) + + T = _fast_cross(A, B) + # Normalization + l = np.sqrt(np.sum(T ** 2, axis=-1)) + l = np.expand_dims(l, 2) + # Might get some divide-by-zeros, but we don't care + with np.errstate(invalid='ignore'): + TN = T / l + return TN + + +def intersection(A, B, C, D): + r""" + Returns the point of intersection between two great circle arcs. + The arcs are defined between the points *AB* and *CD*. Either *A* + and *B* or *C* and *D* may be arrays of points, but not both. + + Parameters + ---------- + A, B : (*x*, *y*, *z*) triples or Nx3 arrays of triples + Endpoints of the first great circle arc. + + C, D : (*x*, *y*, *z*) triples or Nx3 arrays of triples + Endpoints of the second great circle arc. + + Returns + ------- + T : (*x*, *y*, *z*) triples or Nx3 arrays of triples + If the given arcs intersect, the intersection is returned. If + the arcs do not intersect, the triple is set to all NaNs. + + Notes + ----- + The basic intersection is computed using linear algebra as follows + [1]_: + + .. math:: + + T = \lVert(A × B) × (C × D)\rVert + + To determine the correct sign (i.e. hemisphere) of the + intersection, the following four values are computed: + + .. math:: + + s_1 = ((A × B) × A) · T + + s_2 = (B × (A × B)) · T + + s_3 = ((C × D) × C) · T + + s_4 = (D × (C × D)) · T + + For :math:`s_n`, if all positive :math:`T` is returned as-is. If + all negative, :math:`T` is multiplied by :math:`-1`. Otherwise + the intersection does not exist and is undefined. + + References + ---------- + + .. [1] Method explained in an `e-mail + <http://www.mathworks.com/matlabcentral/newsreader/view_thread/276271>`_ + by Roger Stafford. + + http://www.mathworks.com/matlabcentral/newsreader/view_thread/276271 + """ + if HAS_C_UFUNCS: + return math_util.intersection(A, B, C, D) + + A = np.asanyarray(A) + B = np.asanyarray(B) + C = np.asanyarray(C) + D = np.asanyarray(D) + + A, B = np.broadcast_arrays(A, B) + C, D = np.broadcast_arrays(C, D) + + ABX = _fast_cross(A, B) + CDX = _fast_cross(C, D) + T = _cross_and_normalize(ABX, CDX) + T_ndim = len(T.shape) + + if T_ndim > 1: + s = np.zeros(T.shape[0]) + else: + s = np.zeros(1) + s += np.sign(inner1d(_fast_cross(ABX, A), T)) + s += np.sign(inner1d(_fast_cross(B, ABX), T)) + s += np.sign(inner1d(_fast_cross(CDX, C), T)) + s += np.sign(inner1d(_fast_cross(D, CDX), T)) + if T_ndim > 1: + s = np.expand_dims(s, 2) + + cross = np.where(s == -4, -T, np.where(s == 4, T, np.nan)) + + # If they share a common point, it's not an intersection. This + # gets around some rounding-error/numerical problems with the + # above. + equals = (np.all(A == C, axis = -1) | + np.all(A == D, axis = -1) | + np.all(B == C, axis = -1) | + np.all(B == D, axis = -1)) + + equals = np.expand_dims(equals, 2) + + return np.where(equals, np.nan, cross) + + +def length(A, B, degrees=True): + r""" + Returns the angular distance between two points (in vector space) + on the unit sphere. + + Parameters + ---------- + A, B : (*x*, *y*, *z*) triples or Nx3 arrays of triples + The endpoints of the great circle arc, in vector space. + + degrees : bool, optional + If `True` (default) the result is returned in decimal degrees, + otherwise radians. + + Returns + ------- + length : scalar or array of scalars + The angular length of the great circle arc. + + Notes + ----- + The length is computed using the following: + + .. math:: + + \Delta = \arccos(A \dot B) + """ + A = np.asanyarray(A) + B = np.asanyarray(B) + + dot = inner1d(A, B) + dot = np.clip(dot, -1.0, 1.0) + with np.errstate(invalid='ignore'): + result = np.arccos(dot) + + if degrees: + return np.rad2deg(result) + else: + return result + + +def intersects(A, B, C, D): + """ + Returns `True` if the great circle arcs between *AB* and *CD* + intersect. Either *A* and *B* or *C* and *D* may be arrays of + points, but not both. + + Parameters + ---------- + A, B : (*x*, *y*, *z*) triples or Nx3 arrays of triples + Endpoints of the first great circle arc. + + C, D : (*x*, *y*, *z*) triples or Nx3 arrays of triples + Endpoints of the second great circle arc. + + Returns + ------- + intersects : bool or array of bool + If the given arcs intersect, the intersection is returned as + `True`. + """ + with np.errstate(invalid='ignore'): + intersections = intersection(A, B, C, D) + return np.isfinite(intersections[..., 0]) + + +def angle(A, B, C, degrees=True): + """ + Returns the angle at *B* between *AB* and *BC*. + + This always returns the shortest angle < π. + + Parameters + ---------- + A, B, C : (*x*, *y*, *z*) triples or Nx3 arrays of triples + Points on sphere. + + degrees : bool, optional + If `True` (default) the result is returned in decimal degrees, + otherwise radians. + + Returns + ------- + angle : float or array of floats + The angle at *B* between *AB* and *BC*. + + References + ---------- + + .. [1] Miller, Robert D. Computing the area of a spherical + polygon. Graphics Gems IV. 1994. Academic Press. + """ + A = np.asanyarray(A) + B = np.asanyarray(B) + C = np.asanyarray(C) + + A, B, C = np.broadcast_arrays(A, B, C) + + ABX = _fast_cross(A, B) + ABX = _cross_and_normalize(B, ABX) + BCX = _fast_cross(C, B) + BCX = _cross_and_normalize(B, BCX) + with np.errstate(invalid='ignore'): + angle = np.arccos(inner1d(ABX, BCX)) + + if degrees: + angle = np.rad2deg(angle) + return angle + + +def midpoint(A, B): + """ + Returns the midpoint on the great circle arc between *A* and *B*. + + Parameters + ---------- + A, B : (*x*, *y*, *z*) triples or Nx3 arrays of triples + The endpoints of the great circle arc. It is assumed that + these points are already normalized. + + Returns + ------- + midpoint : (*x*, *y*, *z*) triple or Nx3 arrays of triples + The midpoint between *A* and *B*, normalized on the unit + sphere. + """ + P = (A + B) / 2.0 + # Now normalize... + l = np.sqrt(np.sum(P * P, axis=-1)) + l = np.expand_dims(l, 2) + return P / l + + +def interpolate(A, B, steps=50): + r""" + Interpolate along the great circle arc. + + Parameters + ---------- + A, B : (*x*, *y*, *z*) triples or Nx3 arrays of triples + The endpoints of the great circle arc. It is assumed thats + these points are already normalized. + + steps : int + The number of interpolation steps + + Returns + ------- + array : (*x*, *y*, *z*) triples + The points interpolated along the great circle arc + + Notes + ----- + + This uses Slerp interpolation where *Ω* is the angle subtended by + the arc, and *t* is the parameter 0 <= *t* <= 1. + + .. math:: + + \frac{\sin((1 - t)\Omega)}{\sin \Omega}A + \frac{\sin(t \Omega)}{\sin \Omega}B + """ + steps = max(steps, 2) + t = np.linspace(0.0, 1.0, steps, endpoint=True).reshape((steps, 1)) + + omega = length(A, B, degrees=False) + if omega == 0.0: + offsets = t + else: + sin_omega = np.sin(omega) + offsets = np.sin(t * omega) / sin_omega + + return offsets[::-1] * A + offsets * B diff --git a/stsci/sphere/polygon.py b/stsci/sphere/polygon.py new file mode 100644 index 0000000..e9ad379 --- /dev/null +++ b/stsci/sphere/polygon.py @@ -0,0 +1,831 @@ +# -*- 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. + +""" +The `polygon` module defines the `SphericalPolygon` class for managing +polygons on the unit sphere. +""" +from __future__ import division, print_function, unicode_literals, absolute_import + +# STDLIB +from copy import copy, deepcopy + +# THIRD-PARTY +import numpy as np + +# LOCAL +from . import great_circle_arc +from . import vector + +__all__ = ['SphericalPolygon'] + + +class SphericalPolygon(object): + r""" + Polygons are represented by both a set of points (in Cartesian + (*x*, *y*, *z*) normalized on the unit sphere), and an inside + point. The inside point is necessary, because both the inside and + outside of the polygon are finite areas on the great sphere, and + therefore we need a way of specifying which is which. + """ + + def __init__(self, points, inside=None): + r""" + Parameters + ---------- + points : An Nx3 array of (*x*, *y*, *z*) triples in vector space + These points define the boundary of the polygon. It must + be "closed", i.e., the last point is the same as the first. + + It may contain zero points, in which it defines the null + polygon. It may not contain one, two or three points. + Four points are needed to define a triangle, since the + polygon must be closed. + + inside : An (*x*, *y*, *z*) triple, optional + This point must be inside the polygon. If not provided, the + mean of the points will be used. + """ + if len(points) == 0: + # handle special case of initializing with an empty list of + # vertices (ticket #1079). + self._inside = np.zeros(3) + self._points = np.asanyarray(points) + return + elif len(points) < 3: + raise ValueError("Polygon made of too few points") + else: + assert np.array_equal(points[0], points[-1]), 'Polygon is not closed' + + self._points = np.asanyarray(points) + + if inside is None: + self._inside = np.mean(points[:-1], axis=0) + else: + self._inside = np.asanyarray(inside) + + # TODO: Detect self-intersection and fix + + def __copy__(self): + return deepcopy(self) + + def __repr__(self): + return '%s(%r, %r)' % (self.__class__.__name__, + self.points, self.inside) + + @property + def points(self): + """ + The points defining the polygon. It is an Nx3 array of + (*x*, *y*, *z*) vectors. The polygon will be explicitly + closed, i.e., the first and last points are the same. + """ + return self._points + + @property + def inside(self): + """ + Get the inside point of the polygon. + """ + return self._inside + + def to_radec(self): + """ + Convert `SphericalPolygon` footprint to RA and DEC. + + Returns + ------- + ra, dec : list of float + List of *ra* and *dec* in degrees corresponding + to `points`. + """ + return vector.vector_to_radec(self.points[:,0], self.points[:,1], + self.points[:,2], degrees=True) + + @classmethod + def from_radec(cls, ra, dec, center=None, degrees=True): + r""" + Create a new `SphericalPolygon` from a list of (*ra*, *dec*) + points. + + Parameters + ---------- + ra, dec : 1-D arrays of the same length + The vertices of the polygon in right-ascension and + declination. It must be \"closed\", i.e., that is, the + last point is the same as the first. + + center : (*ra*, *dec*) pair, optional + A point inside of the polygon to define its inside. If no + *center* point is provided, the mean of the polygon's + points in vector space will be used. That approach may + not work for concave polygons. + + degrees : bool, optional + If `True`, (default) *ra* and *dec* are in decimal degrees, + otherwise in radians. + + Returns + ------- + polygon : `SphericalPolygon` object + """ + # Convert to Cartesian + x, y, z = vector.radec_to_vector(ra, dec, degrees=degrees) + + if center is None: + xc = x.mean() + yc = y.mean() + zc = z.mean() + center = vector.normalize_vector(xc, yc, zc) + else: + center = vector.radec_to_vector(*center, degrees=degrees) + + return cls(np.dstack((x, y, z))[0], center) + + @classmethod + def from_cone(cls, ra, dec, radius, degrees=True, steps=16.0): + r""" + Create a new `SphericalPolygon` from a cone (otherwise known + as a "small circle") defined using (*ra*, *dec*, *radius*). + + The cone is not represented as an ideal circle on the sphere, + but as a series of great circle arcs. The resolution of this + conversion can be controlled using the *steps* parameter. + + Parameters + ---------- + ra, dec : float scalars + This defines the center of the cone + + radius : float scalar + The radius of the cone + + degrees : bool, optional + If `True`, (default) *ra*, *dec* and *radius* are in + decimal degrees, otherwise in radians. + + steps : int, optional + The number of steps to use when converting the small + circle to a polygon. + + Returns + ------- + polygon : `SphericalPolygon` object + """ + u, v, w = vector.radec_to_vector(ra, dec, degrees=degrees) + if degrees: + radius = np.deg2rad(radius) + + # Get an arbitrary perpendicular vector. This be be obtained + # by crossing (u, v, w) with any unit vector that is not itself. + which_min = np.argmin([u, v, w]) + if which_min == 0: + perp = np.cross([u, v, w], [1., 0., 0.]) + elif which_min == 1: + perp = np.cross([u, v, w], [0., 1., 0.]) + else: + perp = np.cross([u, v, w], [0., 0., 1.]) + + # Rotate by radius around the perpendicular vector to get the + # "pen" + x, y, z = vector.rotate_around( + u, v, w, perp[0], perp[1], perp[2], radius, degrees=False) + + # Then rotate the pen around the center point all 360 degrees + C = np.linspace(0, np.pi * 2.0, steps) + # Ensure that the first and last elements are exactly the + # same. 2π should equal 0, but with rounding error that isn't + # always the case. + C[-1] = 0 + C = C[::-1] + X, Y, Z = vector.rotate_around(x, y, z, u, v, w, C, degrees=False) + + return cls(np.dstack((X, Y, Z))[0], (u, v, w)) + + @classmethod + def from_wcs(cls, fitspath, steps=1, crval=None): + r""" + Create a new `SphericalPolygon` from the footprint of a FITS + WCS specification. + + This method requires having `pywcs` and `pyfits` installed. + + Parameters + ---------- + fitspath : path to a FITS file, `pyfits.Header`, or `pywcs.WCS` + Refers to a FITS header containing a WCS specification. + + steps : int, optional + The number of steps along each edge to convert into + polygon edges. + + Returns + ------- + polygon : `SphericalPolygon` object + """ + import pywcs + import pyfits + + if isinstance(fitspath, pyfits.Header): + header = fitspath + wcs = pywcs.WCS(header) + elif isinstance(fitspath, pywcs.WCS): + wcs = fitspath + else: + wcs = pywcs.WCS(fitspath) + if crval is not None: + wcs.wcs.crval = crval + xa, ya = [wcs.naxis1, wcs.naxis2] + + length = steps * 4 + 1 + X = np.empty(length) + Y = np.empty(length) + + # Now define each of the 4 edges of the quadrilateral + X[0 :steps ] = np.linspace(1, xa, steps, False) + Y[0 :steps ] = 1 + X[steps :steps*2] = xa + Y[steps :steps*2] = np.linspace(1, ya, steps, False) + X[steps*2:steps*3] = np.linspace(xa, 1, steps, False) + Y[steps*2:steps*3] = ya + X[steps*3:steps*4] = 1 + Y[steps*3:steps*4] = np.linspace(ya, 1, steps, False) + X[-1] = 1 + Y[-1] = 1 + + # Use wcslib to convert to (ra, dec) + ra, dec = wcs.all_pix2sky(X, Y, 1) + + # Convert to Cartesian + x, y, z = vector.radec_to_vector(ra, dec) + + # Calculate an inside point + ra, dec = wcs.all_pix2sky(xa / 2.0, ya / 2.0, 1) + xc, yc, zc = vector.radec_to_vector(ra, dec) + + return cls(np.dstack((x, y, z))[0], (xc[0], yc[0], zc[0])) + + def _unique_points(self): + """ + Return a copy of `points` with duplicates removed. + Order is preserved. + + .. note:: Output cannot be used to build a new + polygon. + """ + val = [] + for p in self.points: + v = tuple(p) + if v not in val: + val.append(v) + return np.array(val) + + def _sorted_points(self, preserve_order=True, unique=False): + """ + Return a copy of `points` sorted such that smallest + (*x*, *y*, *z*) is on top. + + .. note:: Output cannot be used to build a new + polygon. + + Parameters + ---------- + preserve_order : bool + Preserve original order? If `True`, polygon is + rotated around min point. If `False`, all points + are sorted in ascending order. + + unique : bool + Exclude duplicates. + """ + if len(self.points) == 0: + return [] + + if unique: + pts = self._unique_points() + else: + pts = self.points + + idx = np.lexsort((pts[:,0], pts[:,1], pts[:,2])) + + if preserve_order: + i_min = idx[0] + val = np.vstack([pts[i_min:], pts[:i_min]]) + else: + val = pts[idx] + + return val + + def same_points_as(self, other, do_sort=True, thres=0.01): + """ + Determines if this `SphericalPolygon` points are the same + as the other. Number of points and areas are also compared. + + When `do_sort` is `True`, even when *self* and *other* + have same points, they might not be equivalent because + the order of the points defines the polygon. + + Parameters + ---------- + other : `SphericalPolygon` + + do_sort : bool + Compare sorted unique points. + + thres : float + Fraction of area to use in equality decision. + + Returns + ------- + is_eq : bool + `True` or `False`. + """ + self_n = len(self.points) + + if self_n != len(other.points): + return False + + if self_n == 0: + return True + + self_a = self.area() + is_same_limit = thres * self_a + + if np.abs(self_a - other.area()) > is_same_limit: + return False + + if do_sort: + self_pts = self._sorted_points(preserve_order=False, unique=True) + other_pts = other._sorted_points(preserve_order=False, unique=True) + else: + self_pts = self.points + other_pts = other.points + + is_eq = True + + for self_p, other_p in zip(self_pts, other_pts): + x_sum = 0.0 + + for a,b in zip(self_p, other_p): + x_sum += (a - b) ** 2 + + if np.sqrt(x_sum) > is_same_limit: + is_eq = False + break + + return is_eq + + def contains_point(self, point): + r""" + Determines if this `SphericalPolygon` contains a given point. + + Parameters + ---------- + point : an (*x*, *y*, *z*) triple + The point to test. + + Returns + ------- + contains : bool + Returns `True` if the polygon contains the given *point*. + """ + P = self._points + r = self._inside + point = np.asanyarray(point) + + intersects = great_circle_arc.intersects(P[:-1], P[1:], r, point) + crossings = np.sum(intersects) + + return (crossings % 2) == 0 + + def intersects_poly(self, other): + r""" + Determines if this `SphericalPolygon` intersects another + `SphericalPolygon`. + + This method is much faster than actually computing the + intersection region between two polygons. + + Parameters + ---------- + other : `SphericalPolygon` + + Returns + ------- + intersects : bool + Returns `True` if this polygon intersects the *other* + polygon. + + Notes + ----- + + The algorithm proceeds as follows: + + 1. Determine if any single point of one polygon is contained + within the other. + + 2. Deal with the case where only the edges overlap as in:: + + : o---------o + : o----+---------+----o + : | | | | + : o----+---------+----o + : o---------o + + In this case, an edge from one polygon must cross an + edge from the other polygon. + """ + assert isinstance(other, SphericalPolygon) + + # The easy case is in which a point of one polygon is + # contained in the other polygon. + for point in other._points: + if self.contains_point(point): + return True + for point in self._points: + if other.contains_point(point): + return True + + # The hard case is when only the edges overlap, as in: + # + # o---------o + # o----+---------+----o + # | | | | + # o----+---------+----o + # o---------o + # + for i in range(len(self._points) - 1): + A = self._points[i] + B = self._points[i+1] + if np.any(great_circle_arc.intersects( + A, B, other._points[:-1], other._points[1:])): + return True + return False + + def intersects_arc(self, a, b): + """ + Determines if this `SphericalPolygon` intersects or contains + the given arc. + """ + return self.contains_point(great_circle_arc.midpoint(a, b)) + + if self.contains_point(a) and self.contains_point(b): + return True + + P = self._points + intersects = great_circle_arc.intersects(P[:-1], P[1:], a, b) + if np.any(intersects): + return True + return False + + def area(self): + r""" + Returns the area of the polygon on the unit sphere. + + The algorithm is not able to compute the area of polygons + that are larger than half of the sphere. Therefore, the + area will always be less than 2π. + + The area is computed by transforming the polygon to two + dimensions using the `Lambert azimuthal equal-area projection + <http://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection>`_ + + .. math:: + + X = \sqrt{\frac{2}{1-z}}x + + .. math:: + + Y = \sqrt{\frac{2}{1-z}}y + + The individual great arc circle segments are interpolated + before doing the transformation so that the curves are not + straightened in the process. + + It then uses a standard 2D algorithm to compute the area. + + .. math:: + + A = \left| \sum^n_{i=0} X_i Y_{i+1} - X_{i+1}Y_i \right| + """ + if len(self._points) < 3: + #return np.float64(0.0) + return np.array(0.0) + + points = self._points.copy() + + # Rotate polygon so that center of polygon is at north pole + centroid = np.mean(points[:-1], axis=0) + centroid = vector.normalize_vector(*centroid) + points = self._points - (centroid + np.array([0, 0, 1])) + vector.normalize_vector( + points[:, 0], points[:, 1], points[:, 2], inplace=True) + + X = [] + Y = [] + for A, B in zip(points[:-1], points[1:]): + length = great_circle_arc.length(A, B, degrees=True) + interp = great_circle_arc.interpolate(A, B, length * 4) + x, y, z = vector.normalize_vector( + interp[:, 0], interp[:, 1], interp[:, 2], inplace=True) + x, y = vector.equal_area_proj(x, y, z) + X.extend(x) + Y.extend(y) + + X = np.array(X) + Y = np.array(Y) + + return np.abs(np.sum(X[:-1] * Y[1:] - X[1:] * Y[:-1]) * 0.5 * np.pi) + + def union(self, other): + """ + Return a new `SphericalPolygon` that is the union of *self* + and *other*. + + If the polygons are disjoint, they result will be connected + using cut lines. For example:: + + : o---------o + : | | + : o---------o=====o----------o + : | | + : o----------o + + Parameters + ---------- + other : `SphericalPolygon` + + Returns + ------- + polygon : `SphericalPolygon` object + + See also + -------- + multi_union + + Notes + ----- + For implementation details, see the :mod:`~sphere.graph` + module. + """ + from . import graph + g = graph.Graph([self, other]) + + polygon = g.union() + + return self.__class__(polygon, self._inside) + + @classmethod + def multi_union(cls, polygons): + """ + Return a new `SphericalPolygon` that is the union of all of the + polygons in *polygons*. + + Parameters + ---------- + polygons : sequence of `SphericalPolygon` + + Returns + ------- + polygon : `SphericalPolygon` object + + See also + -------- + union + """ + assert len(polygons) + for polygon in polygons: + assert isinstance(polygon, SphericalPolygon) + + from . import graph + + g = graph.Graph(polygons) + polygon = g.union() + return cls(polygon, polygons[0]._inside) + + @staticmethod + def _find_new_inside(points, polygons): + """ + Finds an acceptable inside point inside of *points* that is + also inside of *polygons*. Used by the intersection + algorithm, and is really only useful in that context because + it requires existing polygons with known inside points. + """ + if len(points) < 4: + return np.array([0, 0, 0]) + + # Special case for a triangle + if len(points) == 4: + return np.sum(points[:3]) / 3.0 + + for i in range(len(points) - 1): + A = points[i] + # Skip the adjacent point, since it is by definition on + # the edge of the polygon, not potentially running through + # the middle. + for j in range(i + 2, len(points) - 1): + B = points[j] + C = great_circle_arc.midpoint(A, B) + in_all = True + for polygon in polygons: + if not polygon.contains_point(C): + in_all = False + break + if in_all: + return C + + raise RuntimeError("Suitable inside point could not be found") + + def intersection(self, other): + """ + Return a new `SphericalPolygon` that is the intersection of + *self* and *other*. + + If the intersection is empty, a `SphericalPolygon` with zero + points will be returned. + + If the result is disjoint, the pieces will be connected using + cut lines. For example:: + + : o---------o + : | | + : o---------o=====o----------o + : | | + : o----------o + + Parameters + ---------- + other : `SphericalPolygon` + + Returns + ------- + polygon : `SphericalPolygon` object + + Notes + ----- + For implementation details, see the :mod:`~sphere.graph` + module. + """ + # if not self.intersects_poly(other): + # return self.__class__([], [0, 0, 0]) + + from . import graph + if len(self._points) < 3 or len(other._points) < 3: + return self.__class__([], [0, 0, 0]) + + g = graph.Graph([self, other]) + + polygon = g.intersection() + + inside = self._find_new_inside(polygon, [self, other]) + + return self.__class__(polygon, inside) + + @classmethod + def multi_intersection(cls, polygons, method='parallel'): + """ + Return a new `SphericalPolygon` that is the intersection of + all of the polygons in *polygons*. + + Parameters + ---------- + polygons : sequence of `SphericalPolygon` + + method : 'parallel' or 'serial', optional + Specifies the method that is used to perform the + intersections: + + - 'parallel' (default): A graph is built using all of + the polygons, and the intersection operation is computed on + the entire thing globally. + + - 'serial': The polygon is built in steps by adding one + polygon at a time and computing the intersection at + each step. + + This option is provided because one may be faster than the + other depending on context, but it primarily exposed for + testing reasons. Both modes should theoretically provide + equivalent results. + + Returns + ------- + polygon : `SphericalPolygon` object + """ + assert len(polygons) + for polygon in polygons: + assert isinstance(polygon, SphericalPolygon) + + # for i in range(len(polygons)): + # polyA = polygons[i] + # for j in range(i + 1, len(polygons)): + # polyB = polygons[j] + # if not polyA.intersects_poly(polyB): + # return cls([], [0, 0, 0]) + + from . import graph + + if method.lower() == 'parallel': + g = graph.Graph(polygons) + polygon = g.intersection() + inside = cls._find_new_inside(polygon, polygons) + return cls(polygon, inside) + elif method.lower() == 'serial': + result = copy(polygons[0]) + for polygon in polygons[1:]: + result = result.intersection(polygon) + # If we have a null intersection already, we don't + # need to go any further. + if len(result._points) < 3: + return result + return result + else: + raise ValueError("method must be 'parallel' or 'serial'") + + def overlap(self, other): + r""" + Returns the fraction of *self* that is overlapped by *other*. + + Let *self* be *a* and *other* be *b*, then the overlap is + defined as: + + .. math:: + + \frac{S_a}{S_{a \cap b}} + + Parameters + ---------- + other : `SphericalPolygon` + + Returns + ------- + frac : float + The fraction of *self* that is overlapped by *other*. + """ + s1 = self.area() + intersection = self.intersection(other) + s2 = intersection.area() + return s2 / s1 + + def draw(self, m, **plot_args): + """ + Draws the polygon in a matplotlib.Basemap axes. + + Parameters + ---------- + m : Basemap axes object + + **plot_args : Any plot arguments to pass to basemap + """ + if not len(self._points): + return + if not len(plot_args): + plot_args = {'color': 'blue'} + points = self._points + + for A, B in zip(points[0:-1], points[1:]): + length = great_circle_arc.length(A, B, degrees=True) + if not np.isfinite(length): + length = 2 + interpolated = great_circle_arc.interpolate(A, B, length * 4) + ra, dec = vector.vector_to_radec( + interpolated[:, 0], interpolated[:, 1], interpolated[:, 2], + degrees=True) + for r0, d0, r1, d1 in zip(ra[0:-1], dec[0:-1], ra[1:], dec[1:]): + m.drawgreatcircle(r0, d0, r1, d1, **plot_args) + + ra, dec = vector.vector_to_radec( + *self._inside, degrees=True) + x, y = m(ra, dec) + m.scatter(x, y, 1, **plot_args) diff --git a/stsci/sphere/test/__init__.py b/stsci/sphere/test/__init__.py new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/stsci/sphere/test/__init__.py @@ -0,0 +1 @@ + diff --git a/stsci/sphere/test/benchmarks.py b/stsci/sphere/test/benchmarks.py new file mode 100644 index 0000000..91ef8db --- /dev/null +++ b/stsci/sphere/test/benchmarks.py @@ -0,0 +1,39 @@ +import os +import sys +import time + +import numpy as np +from sphere import * +from test_util import * +from test_shared import resolve_imagename + +def point_in_poly_lots(): + image_name = resolve_imagename(ROOT_DIR,'1904-66_TAN.fits') + + poly1 = SphericalPolygon.from_wcs(image_name, 64, crval=[0, 87]) + poly2 = SphericalPolygon.from_wcs(image_name, 64, crval=[20, 89]) + poly3 = SphericalPolygon.from_wcs(image_name, 64, crval=[180, 89]) + + points = get_point_set(density=25) + + count = 0 + for point in points: + if poly1.contains_point(point) or poly2.contains_point(point) or \ + poly3.contains_point(point): + count += 1 + + assert count == 5 + assert poly1.intersects_poly(poly2) + assert not poly1.intersects_poly(poly3) + assert not poly2.intersects_poly(poly3) + +if __name__ == '__main__': + for benchmark in [point_in_poly_lots]: + t = time.time() + sys.stdout.write(benchmark.__name__) + sys.stdout.write('...') + sys.stdout.flush() + benchmark() + sys.stdout.write(' %.03fs\n' % (time.time() - t)) + sys.stdout.flush() + diff --git a/stsci/sphere/test/data/1904-66_TAN.fits.gz b/stsci/sphere/test/data/1904-66_TAN.fits.gz Binary files differnew file mode 100644 index 0000000..34ff168 --- /dev/null +++ b/stsci/sphere/test/data/1904-66_TAN.fits.gz diff --git a/stsci/sphere/test/data/2chipA.fits.gz b/stsci/sphere/test/data/2chipA.fits.gz Binary files differnew file mode 100644 index 0000000..27a68ab --- /dev/null +++ b/stsci/sphere/test/data/2chipA.fits.gz diff --git a/stsci/sphere/test/data/2chipB.fits.gz b/stsci/sphere/test/data/2chipB.fits.gz Binary files differnew file mode 100644 index 0000000..209d5f9 --- /dev/null +++ b/stsci/sphere/test/data/2chipB.fits.gz diff --git a/stsci/sphere/test/data/2chipC.fits.gz b/stsci/sphere/test/data/2chipC.fits.gz Binary files differnew file mode 100644 index 0000000..10af627 --- /dev/null +++ b/stsci/sphere/test/data/2chipC.fits.gz diff --git a/stsci/sphere/test/test.py b/stsci/sphere/test/test.py new file mode 100644 index 0000000..b551090 --- /dev/null +++ b/stsci/sphere/test/test.py @@ -0,0 +1,277 @@ +from __future__ import absolute_import + +import os +import random + +import numpy as np +from numpy.testing import assert_almost_equal, assert_array_less + +from .. import graph +from .. import great_circle_arc +from .. import polygon +from .. import vector + +from .test_util import * +from .test_shared import resolve_imagename + +graph.DEBUG = True + + +def test_normalize_vector(): + x, y, z = np.ogrid[-100:100:11,-100:100:11,-100:100:11] + xn, yn, zn = vector.normalize_vector(x.flatten(), y.flatten(), z.flatten()) + l = np.sqrt(xn ** 2 + yn ** 2 + zn ** 2) + assert_almost_equal(l, 1.0) + + +def test_radec_to_vector(): + npx, npy, npz = vector.radec_to_vector(np.arange(-360, 360, 1), 90) + assert_almost_equal(npx, 0.0) + assert_almost_equal(npy, 0.0) + assert_almost_equal(npz, 1.0) + + spx, spy, spz = vector.radec_to_vector(np.arange(-360, 360, 1), -90) + assert_almost_equal(spx, 0.0) + assert_almost_equal(spy, 0.0) + assert_almost_equal(spz, -1.0) + + eqx, eqy, eqz = vector.radec_to_vector(np.arange(-360, 360, 1), 0) + assert_almost_equal(eqz, 0.0) + + +def test_vector_to_radec(): + ra, dec = vector.vector_to_radec(0, 0, 1) + assert_almost_equal(dec, 90) + + ra, dec = vector.vector_to_radec(0, 0, -1) + assert_almost_equal(dec, -90) + + ra, dec = vector.vector_to_radec(1, 1, 0) + assert_almost_equal(ra, 45.0) + assert_almost_equal(dec, 0.0) + + +def test_intersects_poly_simple(): + ra1 = [-10, 10, 10, -10, -10] + dec1 = [30, 30, 0, 0, 30] + poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1) + + ra2 = [-5, 15, 15, -5, -5] + dec2 = [20, 20, -10, -10, 20] + poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2) + + assert poly1.intersects_poly(poly2) + + # Make sure it isn't order-dependent + ra1 = ra1[::-1] + dec1 = dec1[::-1] + poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1) + + ra2 = ra2[::-1] + dec2 = dec2[::-1] + poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2) + + assert poly1.intersects_poly(poly2) + + +def test_intersects_poly_fully_contained(): + ra1 = [-10, 10, 10, -10, -10] + dec1 = [30, 30, 0, 0, 30] + poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1) + + ra2 = [-5, 5, 5, -5, -5] + dec2 = [20, 20, 10, 10, 20] + poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2) + + assert poly1.intersects_poly(poly2) + + # Make sure it isn't order-dependent + ra1 = ra1[::-1] + dec1 = dec1[::-1] + poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1) + + ra2 = ra2[::-1] + dec2 = dec2[::-1] + poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2) + + assert poly1.intersects_poly(poly2) + + +def test_hard_intersects_poly(): + ra1 = [-10, 10, 10, -10, -10] + dec1 = [30, 30, 0, 0, 30] + poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1) + + ra2 = [-20, 20, 20, -20, -20] + dec2 = [20, 20, 10, 10, 20] + poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2) + + assert poly1.intersects_poly(poly2) + + # Make sure it isn't order-dependent + ra1 = ra1[::-1] + dec1 = dec1[::-1] + poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1) + + ra2 = ra2[::-1] + dec2 = dec2[::-1] + poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2) + + assert poly1.intersects_poly(poly2) + + +def test_not_intersects_poly(): + ra1 = [-10, 10, 10, -10, -10] + dec1 = [30, 30, 5, 5, 30] + poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1) + + ra2 = [-20, 20, 20, -20, -20] + dec2 = [-20, -20, -10, -10, -20] + poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2) + + assert not poly1.intersects_poly(poly2) + + # Make sure it isn't order-dependent + ra1 = ra1[::-1] + dec1 = dec1[::-1] + poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1) + + ra2 = ra2[::-1] + dec2 = dec2[::-1] + poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2) + + assert not poly1.intersects_poly(poly2) + + +def test_point_in_poly(): + point = np.asarray([-0.27475449, 0.47588873, -0.83548781]) + points = np.asarray([[ 0.04821217, -0.29877206, 0.95310589], + [ 0.04451801, -0.47274119, 0.88007608], + [-0.14916503, -0.46369786, 0.87334649], + [-0.16101648, -0.29210164, 0.94273555], + [ 0.04821217, -0.29877206, 0.95310589]]) + inside = np.asarray([-0.03416009, -0.36858623, 0.9289657]) + poly = polygon.SphericalPolygon(points, inside) + assert not poly.contains_point(point) + + +def test_point_in_poly_lots(): + import pyfits + fits = pyfits.open(resolve_imagename(ROOT_DIR, '1904-66_TAN.fits')) + header = fits[0].header + + poly1 = polygon.SphericalPolygon.from_wcs( + header, 1, crval=[0, 87]) + poly2 = polygon.SphericalPolygon.from_wcs( + header, 1, crval=[20, 89]) + poly3 = polygon.SphericalPolygon.from_wcs( + header, 1, crval=[180, 89]) + + points = get_point_set() + count = 0 + for point in points: + if poly1.contains_point(point) or poly2.contains_point(point) or \ + poly3.contains_point(point): + count += 1 + + assert count == 5 + assert poly1.intersects_poly(poly2) + assert not poly1.intersects_poly(poly3) + assert not poly2.intersects_poly(poly3) + + +def test_great_circle_arc_intersection(): + A = [-10, -10] + B = [10, 10] + + C = [-25, 10] + D = [15, -10] + + E = [-20, 40] + F = [20, 40] + + correct = [0.99912414, -0.02936109, -0.02981403] + + A = vector.radec_to_vector(*A) + B = vector.radec_to_vector(*B) + C = vector.radec_to_vector(*C) + D = vector.radec_to_vector(*D) + E = vector.radec_to_vector(*E) + F = vector.radec_to_vector(*F) + + assert great_circle_arc.intersects(A, B, C, D) + r = great_circle_arc.intersection(A, B, C, D) + assert r.shape == (3,) + assert_almost_equal(r, correct) + + assert np.all(great_circle_arc.intersects([A], [B], [C], [D])) + r = great_circle_arc.intersection([A], [B], [C], [D]) + assert r.shape == (1, 3) + assert_almost_equal(r, [correct]) + + assert np.all(great_circle_arc.intersects([A], [B], C, D)) + r = great_circle_arc.intersection([A], [B], C, D) + assert r.shape == (1, 3) + assert_almost_equal(r, [correct]) + + assert not np.all(great_circle_arc.intersects([A, E], [B, F], [C], [D])) + r = great_circle_arc.intersection([A, E], [B, F], [C], [D]) + assert r.shape == (2, 3) + assert_almost_equal(r[0], correct) + assert np.all(np.isnan(r[1])) + + # Test parallel arcs + r = great_circle_arc.intersection(A, B, A, B) + assert np.all(np.isnan(r)) + + +def test_great_circle_arc_length(): + A = [90, 0] + B = [-90, 0] + A = vector.radec_to_vector(*A) + B = vector.radec_to_vector(*B) + assert great_circle_arc.length(A, B) == 180.0 + + A = [135, 0] + B = [-90, 0] + A = vector.radec_to_vector(*A) + B = vector.radec_to_vector(*B) + assert great_circle_arc.length(A, B) == 135.0 + + A = [0, 0] + B = [0, 90] + A = vector.radec_to_vector(*A) + B = vector.radec_to_vector(*B) + assert great_circle_arc.length(A, B) == 90.0 + + +def test_great_circle_arc_angle(): + A = [1, 0, 0] + B = [0, 1, 0] + C = [0, 0, 1] + assert great_circle_arc.angle(A, B, C) == 90.0 + + # TODO: More angle tests + + +def test_cone(): + random.seed(0) + for i in range(50): + ra = random.randrange(-180, 180) + dec = random.randrange(20, 90) + cone = polygon.SphericalPolygon.from_cone(ra, dec, 8, steps=64) + + +def test_area(): + triangles = [ + ([[90, 0], [0, -45], [0, 45], [90, 0]], np.pi * 0.5), + ([[90, 0], [0, -22.5], [0, 22.5], [90, 0]], np.pi * 0.25), + ([[90, 0], [0, -11.25], [0, 11.25], [90, 0]], np.pi * 0.125) + ] + + for tri, area in triangles: + tri = np.array(tri) + x, y, z = vector.radec_to_vector(tri[:,1], tri[:,0]) + points = np.dstack((x, y, z))[0] + poly = polygon.SphericalPolygon(points) + calc_area = poly.area() diff --git a/stsci/sphere/test/test_intersection.py b/stsci/sphere/test/test_intersection.py new file mode 100644 index 0000000..6dacb3d --- /dev/null +++ b/stsci/sphere/test/test_intersection.py @@ -0,0 +1,189 @@ +from __future__ import print_function, absolute_import + +# STDLIB +import functools +import itertools +import math +import os +import random +import sys + +# THIRD-PARTY +import numpy as np +from numpy.testing import assert_array_almost_equal + +# LOCAL +from .. import polygon +from .test_shared import resolve_imagename + +GRAPH_MODE = False +ROOT_DIR = os.path.join(os.path.dirname(__file__), 'data') + + +class intersection_test: + def __init__(self, lon_0, lat_0, proj='ortho'): + self._lon_0 = lon_0 + self._lat_0 = lat_0 + self._proj = proj + + def __call__(self, func): + @functools.wraps(func) + def run(*args, **kwargs): + if GRAPH_MODE: + from mpl_toolkits.basemap import Basemap + from matplotlib import pyplot as plt + + polys = func(*args, **kwargs) + + intersections = [] + num_permutations = math.factorial(len(polys)) + step_size = int(max(float(num_permutations) / 20.0, 1.0)) + + areas = [x.area() for x in polys] + + if GRAPH_MODE: + print("%d permutations" % num_permutations) + for method in ('parallel', 'serial'): + for i, permutation in enumerate( + itertools.islice( + itertools.permutations(polys), + None, None, step_size)): + filename = '%s_%s_intersection_%04d.svg' % ( + func.__name__, method, i) + print(filename) + + intersection = polygon.SphericalPolygon.multi_intersection( + permutation, method=method) + intersections.append(intersection) + intersection_area = intersection.area() + if GRAPH_MODE: + fig = plt.figure() + m = Basemap(projection=self._proj, + lon_0=self._lon_0, + lat_0=self._lat_0) + m.drawparallels(np.arange(-90., 90., 20.)) + m.drawmeridians(np.arange(0., 420., 20.)) + m.drawmapboundary(fill_color='white') + + intersection.draw(m, color='red', linewidth=3) + for poly in permutation: + poly.draw(m, color='blue', alpha=0.5) + plt.savefig(filename) + fig.clear() + + assert np.all(intersection_area * 0.9 <= areas) + + lengths = np.array([len(x._points) for x in intersections]) + assert np.all(lengths == [lengths[0]]) + areas = np.array([x.area() for x in intersections]) + assert_array_almost_equal(areas, areas[0], decimal=1) + + return run + + +@intersection_test(0, 90) +def test1(): + import pyfits + import os + + fits = pyfits.open(resolve_imagename(ROOT_DIR,'1904-66_TAN.fits')) + header = fits[0].header + + poly1 = polygon.SphericalPolygon.from_wcs( + header, 1, crval=[0, 87]) + poly2 = polygon.SphericalPolygon.from_wcs( + header, 1, crval=[20, 89]) + + return [poly1, poly2] + + +@intersection_test(0, 90) +def test2(): + poly1 = polygon.SphericalPolygon.from_cone(0, 60, 8, steps=16) + poly2 = polygon.SphericalPolygon.from_cone(0, 68, 8, steps=16) + poly3 = polygon.SphericalPolygon.from_cone(12, 66, 8, steps=16) + return [poly1, poly2, poly3] + + +@intersection_test(0, 90) +def test3(): + import pyfits + fits = pyfits.open(resolve_imagename(ROOT_DIR, '1904-66_TAN.fits')) + header = fits[0].header + + poly1 = polygon.SphericalPolygon.from_wcs( + header, 1, crval=[0, 87]) + poly3 = polygon.SphericalPolygon.from_wcs( + header, 1, crval=[175, 89]) + + return [poly1, poly3] + + +def test4(): + import pyfits + import pywcs + + A = pyfits.open(os.path.join(ROOT_DIR, '2chipA.fits.gz')) + B = pyfits.open(os.path.join(ROOT_DIR, '2chipB.fits.gz')) + + wcs = pywcs.WCS(A[1].header, fobj=A) + chipA1 = polygon.SphericalPolygon.from_wcs(wcs) + wcs = pywcs.WCS(A[4].header, fobj=A) + chipA2 = polygon.SphericalPolygon.from_wcs(wcs) + wcs = pywcs.WCS(B[1].header, fobj=B) + chipB1 = polygon.SphericalPolygon.from_wcs(wcs) + wcs = pywcs.WCS(B[4].header, fobj=B) + chipB2 = polygon.SphericalPolygon.from_wcs(wcs) + + Apoly = chipA1.union(chipA2) + Bpoly = chipB1.union(chipB2) + + X = Apoly.intersection(Bpoly) + + +def test5(): + import pyfits + import pywcs + + A = pyfits.open(os.path.join(ROOT_DIR, '2chipA.fits.gz')) + B = pyfits.open(os.path.join(ROOT_DIR, '2chipB.fits.gz')) + + wcs = pywcs.WCS(A[1].header, fobj=A) + chipA1 = polygon.SphericalPolygon.from_wcs(wcs) + wcs = pywcs.WCS(A[4].header, fobj=A) + chipA2 = polygon.SphericalPolygon.from_wcs(wcs) + wcs = pywcs.WCS(B[1].header, fobj=B) + chipB1 = polygon.SphericalPolygon.from_wcs(wcs) + wcs = pywcs.WCS(B[4].header, fobj=B) + chipB2 = polygon.SphericalPolygon.from_wcs(wcs) + + Apoly = chipA1.union(chipA2) + Bpoly = chipB1.union(chipB2) + + Apoly.overlap(chipB1) + + +@intersection_test(0, 90) +def test6(): + import pyfits + fits = pyfits.open(resolve_imagename(ROOT_DIR, '1904-66_TAN.fits')) + header = fits[0].header + + poly1 = polygon.SphericalPolygon.from_wcs( + header, 1) + poly2 = polygon.SphericalPolygon.from_wcs( + header, 1) + + return [poly1, poly2] + + +if __name__ == '__main__': + if '--profile' not in sys.argv: + GRAPH_MODE = True + from mpl_toolkits.basemap import Basemap + from matplotlib import pyplot as plt + + functions = [(k, v) for k, v in globals().items() if k.startswith('test')] + functions.sort() + for k, v in functions: + v() diff --git a/stsci/sphere/test/test_shared.py b/stsci/sphere/test/test_shared.py new file mode 100644 index 0000000..0853d3a --- /dev/null +++ b/stsci/sphere/test/test_shared.py @@ -0,0 +1,12 @@ +import os + +def resolve_imagename(ROOT_DIR, base_name): + """Resolve image name for tests.""" + + image_name = os.path.join(ROOT_DIR, base_name) + + # Is it zipped? + if not os.path.exists(image_name): + image_name = image_name.replace('.fits', '.fits.gz') + + return image_name diff --git a/stsci/sphere/test/test_skyline.py b/stsci/sphere/test/test_skyline.py new file mode 100644 index 0000000..929f0ba --- /dev/null +++ b/stsci/sphere/test/test_skyline.py @@ -0,0 +1,225 @@ +""" +SkyLine tests. + +:Author: Pey Lian Lim + +:Organization: Space Telescope Science Institute + +Examples +-------- +>>> cd path/to/project +>>> nosetests + +""" +from __future__ import absolute_import + +from copy import copy +from numpy.testing import assert_almost_equal + +from ..vector import vector_to_radec +from ..polygon import SphericalPolygon +from ..skyline import SkyLine + +from .test_util import ROOT_DIR +from .test_shared import resolve_imagename + + +#---------------------------------------# +# Load footprints used by all the tests # +#---------------------------------------# +f_2chipA = resolve_imagename(ROOT_DIR, '2chipA.fits') # ACS/WFC #1 +im_2chipA = SkyLine(f_2chipA) +f_2chipB = resolve_imagename(ROOT_DIR, '2chipB.fits') # ACS/WFC #2 +im_2chipB = SkyLine(f_2chipB) +f_2chipC = resolve_imagename(ROOT_DIR, '2chipC.fits') # WFC3/UVIS +im_2chipC = SkyLine(f_2chipC) +f_66_tan = resolve_imagename(ROOT_DIR, '1904-66_TAN.fits') +im_66_tan = SkyLine(f_66_tan, extname='primary') + + +#----- SHARED FUNCTIONS ----- + +def same_members(mem1, mem2): + assert len(mem1) == len(mem2) + + for m in mem1: + assert m in mem2 + + for m in mem2: + assert m in mem1 + +def subset_members(mem_child, mem_parent): + assert len(mem_parent) > len(mem_child) + + for m in mem_child: + assert m in mem_parent + +def subset_polygon(p_child, p_parent): + """Overlap not working. Do this instead until fixed.""" + assert p_parent.area() >= p_child.area() + assert p_parent.contains_point(p_child.inside) + +def no_polygon(p_child, p_parent): + """Overlap not working. Do this instead until fixed.""" + assert not p_parent.contains_point(p_child.inside) + + +#----- MEMBERSHIP ----- + +def do_member_overlap(im): + for m in im.members: + assert_almost_equal(m.polygon.overlap(im), 1.0) + +def test_membership(): + do_member_overlap(im_2chipA) + do_member_overlap(im_2chipB) + do_member_overlap(im_66_tan) + + assert len(im_2chipA.members) == 1 + assert im_2chipA.members[0].fname == f_2chipA + assert im_2chipA.members[0].ext == (1,4) + + +#----- COPY ----- + +def test_copy(): + a_copy = copy(im_2chipA) + assert a_copy is not im_2chipA + + +#----- SPHERICAL POLYGON RELATED ----- + +def test_sphericalpolygon(): + assert im_2chipA.contains_point(im_2chipA.inside) + + assert im_2chipA.intersects_poly(im_2chipB.polygon) + + assert im_2chipA.intersects_arc(im_2chipA.inside, im_2chipB.inside) + + assert im_2chipA.overlap(im_2chipB) < im_2chipA.overlap(im_2chipA) + + assert_almost_equal(im_2chipA.area(), im_2chipB.area()) + + ra_A, dec_A = im_2chipA.to_radec() + for i in xrange(len(im_2chipA.points)): + p = im_2chipA.points[i] + ra, dec = vector_to_radec(p[0], p[1], p[2], degrees=True) + assert_almost_equal(ra_A[i], ra) + assert_almost_equal(dec_A[i], dec) + + +#----- WCS ----- + +def test_wcs(): + wcs = im_2chipA.to_wcs() + new_p = SphericalPolygon.from_wcs(wcs) + subset_polygon(im_2chipA, new_p) + + +#----- UNION ----- + +def do_add_image(im1, im2): + u1 = im1.add_image(im2) + u2 = im2.add_image(im1) + + assert u1.same_points_as(u2) + same_members(u1.members, u2.members) + + all_mems = im1.members + im2.members + same_members(u1.members, all_mems) + + subset_polygon(im1, u1) + subset_polygon(im2, u1) + +def test_add_image(): + # Dithered + do_add_image(im_2chipA, im_2chipB) + + # Not related + do_add_image(im_2chipA, im_66_tan) + + +#----- INTERSECTION ----- + +def do_intersect_image(im1, im2): + i1 = im1.find_intersection(im2) + i2 = im2.find_intersection(im1) + + assert i1.same_points_as(i2) + same_members(i1.members, i2.members) + + if len(i1.points) > 0: + subset_members(im1.members, i1.members) + subset_members(im2.members, i1.members) + + subset_polygon(i1, im1) + subset_polygon(i1, im2) + +def test_find_intersection(): + # Dithered + do_intersect_image(im_2chipA, im_2chipB) + + # Not related + do_intersect_image(im_2chipA, im_66_tan) + + +#----- SKYLINE OVERLAP ----- + +def test_max_overlap(): + max_s, max_a = im_2chipA.find_max_overlap([im_2chipB, im_2chipC, im_66_tan]) + assert max_s is im_2chipB + assert_almost_equal(max_a, im_2chipA.intersection(im_2chipB).area()) + + max_s, max_a = im_2chipA.find_max_overlap([im_2chipB, im_2chipA]) + assert max_s is im_2chipA + assert_almost_equal(max_a, im_2chipA.area()) + +def test_max_overlap_pair(): + assert SkyLine.max_overlap_pair( + [im_2chipB, im_2chipC, im_2chipA, im_66_tan]) == (im_2chipB, im_2chipA) + + assert SkyLine.max_overlap_pair([im_2chipC, im_2chipA, im_66_tan]) is None + + +#----- INTENDED USE CASE ----- + +def test_science_1(): + mos, inc, exc = SkyLine.mosaic([im_2chipA, im_2chipB, im_2chipC, im_66_tan]) + + assert inc == [f_2chipA, f_2chipB] + assert exc == [f_2chipC, f_66_tan] + + subset_polygon(im_2chipA, mos) + subset_polygon(im_2chipB, mos) + + no_polygon(im_2chipC, mos) + no_polygon(im_66_tan, mos) + +def test_science_2(): + """Like `test_science_1` but different input order.""" + mos, inc, exc = SkyLine.mosaic([im_2chipB, im_66_tan, im_2chipC, im_2chipA]) + + assert inc == [f_2chipB, f_2chipA] + assert exc == [f_66_tan, f_2chipC] + + subset_polygon(im_2chipA, mos) + subset_polygon(im_2chipB, mos) + + no_polygon(im_2chipC, mos) + no_polygon(im_66_tan, mos) + + +#----- UNSTABLE ----- + +def DISABLED_unstable_overlap(): + i1 = im_2chipA.find_intersection(im_2chipB) + i2 = im_2chipB.find_intersection(im_2chipA) + + u1 = im_2chipA.add_image(im_2chipB) + u2 = im_2chipB.add_image(im_2chipA) + + # failed here before - known bug + # failure not always the same due to hash mapping + assert_almost_equal(i1.overlap(u1), 1.0) + assert_almost_equal(i1.overlap(i2), 1.0) + assert_almost_equal(u1.overlap(u2), 1.0) diff --git a/stsci/sphere/test/test_union.py b/stsci/sphere/test/test_union.py new file mode 100644 index 0000000..4b649f6 --- /dev/null +++ b/stsci/sphere/test/test_union.py @@ -0,0 +1,215 @@ +from __future__ import print_function, absolute_import + +# STDLIB +import functools +import itertools +import math +import os +import random +import sys + +# THIRD-PARTY +import numpy as np +from numpy.testing import assert_array_almost_equal + +# LOCAL +from .. import polygon +from .test_shared import resolve_imagename + +GRAPH_MODE = False +ROOT_DIR = os.path.join(os.path.dirname(__file__), 'data') + + +class union_test: + def __init__(self, lon_0, lat_0, proj='ortho'): + self._lon_0 = lon_0 + self._lat_0 = lat_0 + self._proj = proj + + def __call__(self, func): + @functools.wraps(func) + def run(*args, **kwargs): + if GRAPH_MODE: + from mpl_toolkits.basemap import Basemap + from matplotlib import pyplot as plt + + polys = func(*args, **kwargs) + + unions = [] + num_permutations = math.factorial(len(polys)) + step_size = int(max(float(num_permutations) / 20.0, 1.0)) + if GRAPH_MODE: + print("%d permutations" % num_permutations) + + areas = [x.area() for x in polys] + + for i, permutation in enumerate( + itertools.islice( + itertools.permutations(polys), + None, None, step_size)): + filename = '%s_union_%04d.svg' % ( + func.__name__, i) + print(filename) + + union = polygon.SphericalPolygon.multi_union( + permutation) + unions.append(union) + union_area = union.area() + print(union._points) + print(permutation[0]._points) + + if GRAPH_MODE: + fig = plt.figure() + m = Basemap(projection=self._proj, + lon_0=self._lon_0, + lat_0=self._lat_0) + m.drawmapboundary(fill_color='white') + m.drawparallels(np.arange(-90., 90., 20.)) + m.drawmeridians(np.arange(0., 420., 20.)) + + union.draw(m, color='red', linewidth=3) + for poly in permutation: + poly.draw(m, color='blue', alpha=0.5) + plt.savefig(filename) + fig.clear() + + print(union_area, areas) + assert np.all(union_area * 1.1 >= areas) + + lengths = np.array([len(x._points) for x in unions]) + assert np.all(lengths == [lengths[0]]) + areas = np.array([x.area() for x in unions]) + assert_array_almost_equal(areas, areas[0], 1) + + return run + + +@union_test(0, 90) +def test1(): + import pyfits + fits = pyfits.open(resolve_imagename(ROOT_DIR, '1904-66_TAN.fits')) + header = fits[0].header + + poly1 = polygon.SphericalPolygon.from_wcs( + header, 1, crval=[0, 87]) + poly2 = polygon.SphericalPolygon.from_wcs( + header, 1, crval=[20, 89]) + poly3 = polygon.SphericalPolygon.from_wcs( + header, 1, crval=[175, 89]) + poly4 = polygon.SphericalPolygon.from_cone( + 90, 70, 10, steps=50) + + return [poly1, poly2, poly3, poly4] + + +@union_test(0, 90) +def test2(): + poly1 = polygon.SphericalPolygon.from_cone(0, 60, 7, steps=16) + poly2 = polygon.SphericalPolygon.from_cone(0, 72, 7, steps=16) + poly3 = polygon.SphericalPolygon.from_cone(20, 60, 7, steps=16) + poly4 = polygon.SphericalPolygon.from_cone(20, 72, 7, steps=16) + poly5 = polygon.SphericalPolygon.from_cone(35, 55, 7, steps=16) + poly6 = polygon.SphericalPolygon.from_cone(60, 60, 3, steps=16) + return [poly1, poly2, poly3, poly4, poly5, poly6] + + +@union_test(0, 90) +def test3(): + random.seed(0) + polys = [] + for i in range(10): + polys.append(polygon.SphericalPolygon.from_cone( + random.randrange(-180, 180), + random.randrange(20, 90), + random.randrange(5, 16), + steps=16)) + return polys + + +@union_test(0, 15) +def test4(): + random.seed(64) + polys = [] + for i in range(10): + polys.append(polygon.SphericalPolygon.from_cone( + random.randrange(-30, 30), + random.randrange(-15, 60), + random.randrange(5, 16), + steps=16)) + return polys + + +def test5(): + import pyfits + import pywcs + + A = pyfits.open(os.path.join(ROOT_DIR, '2chipA.fits.gz')) + + wcs = pywcs.WCS(A[1].header, fobj=A) + chipA1 = polygon.SphericalPolygon.from_wcs(wcs) + wcs = pywcs.WCS(A[4].header, fobj=A) + chipA2 = polygon.SphericalPolygon.from_wcs(wcs) + + null_union = chipA1.union(chipA2) + + +def test6(): + import pyfits + import pywcs + + A = pyfits.open(os.path.join(ROOT_DIR, '2chipC.fits.gz')) + + wcs = pywcs.WCS(A[1].header, fobj=A) + chipA1 = polygon.SphericalPolygon.from_wcs(wcs) + wcs = pywcs.WCS(A[4].header, fobj=A) + chipA2 = polygon.SphericalPolygon.from_wcs(wcs) + + null_union = chipA1.union(chipA2) + + +@union_test(0, 90) +def test7(): + import pyfits + import pywcs + + A = pyfits.open(os.path.join(ROOT_DIR, '2chipA.fits.gz')) + + wcs = pywcs.WCS(A[1].header, fobj=A) + chipA1 = polygon.SphericalPolygon.from_wcs(wcs) + wcs = pywcs.WCS(A[4].header, fobj=A) + chipA2 = polygon.SphericalPolygon.from_wcs(wcs) + + B = pyfits.open(os.path.join(ROOT_DIR, '2chipB.fits.gz')) + + wcs = pywcs.WCS(B[1].header, fobj=B) + chipB1 = polygon.SphericalPolygon.from_wcs(wcs) + wcs = pywcs.WCS(B[4].header, fobj=B) + chipB2 = polygon.SphericalPolygon.from_wcs(wcs) + + return [chipA1, chipA2, chipB1, chipB2] + + +@union_test(0, 90) +def test8(): + import pyfits + fits = pyfits.open(resolve_imagename(ROOT_DIR, '1904-66_TAN.fits')) + header = fits[0].header + + poly1 = polygon.SphericalPolygon.from_wcs( + header, 1) + poly2 = polygon.SphericalPolygon.from_wcs( + header, 1) + + return [poly1, poly2] + + +if __name__ == '__main__': + if '--profile' not in sys.argv: + GRAPH_MODE = True + from mpl_toolkits.basemap import Basemap + from matplotlib import pyplot as plt + + functions = [(k, v) for k, v in globals().items() if k.startswith('test')] + functions.sort() + for k, v in functions: + v() diff --git a/stsci/sphere/test/test_util.py b/stsci/sphere/test/test_util.py new file mode 100644 index 0000000..6dc0393 --- /dev/null +++ b/stsci/sphere/test/test_util.py @@ -0,0 +1,14 @@ +import os + +import numpy as np +from .. import vector + +ROOT_DIR = os.path.join(os.path.dirname(__file__), 'data') + +def get_point_set(density=25): + points = [] + for i in np.linspace(-85, 85, density, True): + for j in np.linspace(-180, 180, np.cos(np.deg2rad(i)) * density): + points.append([j, i]) + points = np.asarray(points) + return np.dstack(vector.radec_to_vector(points[:,0], points[:,1]))[0] diff --git a/stsci/sphere/vector.py b/stsci/sphere/vector.py new file mode 100644 index 0000000..c6887c4 --- /dev/null +++ b/stsci/sphere/vector.py @@ -0,0 +1,243 @@ +# -*- 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. + +""" +The `sphere.vector` module contains the basic operations for handling +vectors and converting them to and from other representations. +""" + +from __future__ import unicode_literals + +# THIRD-PARTY +import numpy as np + + +__all__ = ['radec_to_vector', 'vector_to_radec', 'normalize_vector', + 'rotate_around'] + + +def radec_to_vector(ra, dec, degrees=True): + r""" + Converts a location on the unit sphere from right-ascension and + declination to an *x*, *y*, *z* vector. + + Parameters + ---------- + ra, dec : scalars or 1-D arrays + + degrees : bool, optional + + If `True`, (default) *ra* and *dec* are in decimal degrees, + otherwise in radians. + + Returns + ------- + x, y, z : tuple of scalars or 1-D arrays of the same length + + Notes + ----- + Where right-ascension is *α* and declination is *δ*: + + .. math:: + x = \cos\alpha \cos\delta + + y = \sin\alpha \cos\delta + + z = \sin\delta + """ + ra = np.asanyarray(ra) + dec = np.asanyarray(dec) + + if degrees: + ra_rad = np.deg2rad(ra) + dec_rad = np.deg2rad(dec) + else: + ra_rad = ra + dec_rad = dec + + cos_dec = np.cos(dec_rad) + + return ( + np.cos(ra_rad) * cos_dec, + np.sin(ra_rad) * cos_dec, + np.sin(dec_rad)) + + +def vector_to_radec(x, y, z, degrees=True): + r""" + Converts a vector to right-ascension and declination. + + Parameters + ---------- + x, y, z : scalars or 1-D arrays + The input vectors + + degrees : bool, optional + If `True` (default) the result is returned in decimal degrees, + otherwise radians. + + Returns + ------- + ra, dec : tuple of scalars or arrays of the same length + + Notes + ----- + Where right-ascension is *α* and declination is + *δ*: + + .. math:: + \alpha = \arctan2(y, x) + + \delta = \arctan2(z, \sqrt{x^2 + y^2}) + """ + x = np.asanyarray(x) + y = np.asanyarray(y) + z = np.asanyarray(z) + + result = ( + np.arctan2(y, x), + np.arctan2(z, np.sqrt(x ** 2 + y ** 2))) + + if degrees: + return np.rad2deg(result[0]), np.rad2deg(result[1]) + else: + return result + + +def normalize_vector(x, y, z, inplace=False): + r""" + Normalizes a vector so it falls on the unit sphere. + + *x*, *y*, *z* may be scalars or 1-D arrays + + Parameters + ---------- + x, y, z : scalars or 1-D arrays of the same length + The input vectors + + inplace : bool, optional + When `True`, the original arrays will be normalized in place, + otherwise a normalized copy is returned. + + Returns + ------- + X, Y, Z : scalars of 1-D arrays of the same length + The normalized output vectors + """ + x = np.asanyarray(x) + y = np.asanyarray(y) + z = np.asanyarray(z) + + l = np.sqrt(x ** 2 + y ** 2 + z ** 2) + + if inplace: + x /= l + y /= l + z /= l + return (x, y, z) + else: + return (x / l, y / l, z / l) + + +def rotate_around(x, y, z, u, v, w, theta, degrees=True): + r""" + Rotates the vector (*x*, *y*, *z*) around the arbitrary axis defined by + vector (*u*, *v*, *w*) by *theta*. + + It is assumed that both (*x*, *y*, *z*) and (*u*, *v*, *w*) are + already normalized. + + Parameters + ---------- + x, y, z : doubles + The normalized vector to rotate + + u, v, w : doubles + The normalized vector to rotate around + + theta : double, or array of doubles + The amount to rotate + + degrees : bool, optional + When `True`, *theta* is given in degrees, otherwise radians. + + Returns + ------- + X, Y, Z : doubles + The rotated vector + """ + if degrees: + theta = np.deg2rad(theta) + + costheta = np.cos(theta) + sintheta = np.sin(theta) + icostheta = 1.0 - costheta + + det = (-u*x - v*y - w*z) + X = (-u*det)*icostheta + x*costheta + (-w*y + v*z)*sintheta + Y = (-v*det)*icostheta + y*costheta + ( w*x - u*z)*sintheta + Z = (-w*det)*icostheta + z*costheta + (-v*x + u*y)*sintheta + + return X, Y, Z + + +def equal_area_proj(x, y, z): + """ + Transform the coordinates to a 2-dimensional plane using the + Lambert azimuthal equal-area projection. + + Parameters + ---------- + x, y, z : scalars or 1-D arrays + The input vectors + + Returns + ------- + X, Y : tuple of scalars or arrays of the same length + + Notes + ----- + + .. math:: + + X = \sqrt{\frac{2}{1-z}}x + + .. math:: + + Y = \sqrt{\frac{2}{1-z}}y + """ + scale = np.sqrt(2.0 / (1.0 - z)) + X = scale * x + Y = scale * y + return X, Y |