diff options
Diffstat (limited to 'lib')
| -rw-r--r-- | lib/__init__.py | 5 | ||||
| -rw-r--r-- | lib/graph.py | 909 | ||||
| -rw-r--r-- | lib/great_circle_arc.py | 368 | ||||
| -rw-r--r-- | lib/polygon.py | 831 | ||||
| -rw-r--r-- | lib/test/__init__.py | 1 | ||||
| -rw-r--r-- | lib/test/benchmarks.py | 39 | ||||
| -rw-r--r-- | lib/test/data/1904-66_TAN.fits.gz | bin | 113196 -> 0 bytes | |||
| -rw-r--r-- | lib/test/data/2chipA.fits.gz | bin | 227579 -> 0 bytes | |||
| -rw-r--r-- | lib/test/data/2chipB.fits.gz | bin | 227493 -> 0 bytes | |||
| -rw-r--r-- | lib/test/data/2chipC.fits.gz | bin | 183594 -> 0 bytes | |||
| -rw-r--r-- | lib/test/test.py | 277 | ||||
| -rw-r--r-- | lib/test/test_intersection.py | 189 | ||||
| -rw-r--r-- | lib/test/test_shared.py | 12 | ||||
| -rw-r--r-- | lib/test/test_skyline.py | 225 | ||||
| -rw-r--r-- | lib/test/test_union.py | 215 | ||||
| -rw-r--r-- | lib/test/test_util.py | 14 | ||||
| -rw-r--r-- | lib/vector.py | 243 |
17 files changed, 0 insertions, 3328 deletions
diff --git a/lib/__init__.py b/lib/__init__.py deleted file mode 100644 index 6a03273..0000000 --- a/lib/__init__.py +++ /dev/null @@ -1,5 +0,0 @@ -import sys - -if sys.version_info[0] >= 3: - # Python 3 compatibility - __builtins__['xrange'] = range diff --git a/lib/graph.py b/lib/graph.py deleted file mode 100644 index d8ccf75..0000000 --- a/lib/graph.py +++ /dev/null @@ -1,909 +0,0 @@ -# -*- 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/lib/great_circle_arc.py b/lib/great_circle_arc.py deleted file mode 100644 index d5e78f2..0000000 --- a/lib/great_circle_arc.py +++ /dev/null @@ -1,368 +0,0 @@ -# -*- 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/lib/polygon.py b/lib/polygon.py deleted file mode 100644 index e9ad379..0000000 --- a/lib/polygon.py +++ /dev/null @@ -1,831 +0,0 @@ -# -*- 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/lib/test/__init__.py b/lib/test/__init__.py deleted file mode 100644 index 8b13789..0000000 --- a/lib/test/__init__.py +++ /dev/null @@ -1 +0,0 @@ - diff --git a/lib/test/benchmarks.py b/lib/test/benchmarks.py deleted file mode 100644 index 91ef8db..0000000 --- a/lib/test/benchmarks.py +++ /dev/null @@ -1,39 +0,0 @@ -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/lib/test/data/1904-66_TAN.fits.gz b/lib/test/data/1904-66_TAN.fits.gz Binary files differdeleted file mode 100644 index 34ff168..0000000 --- a/lib/test/data/1904-66_TAN.fits.gz +++ /dev/null diff --git a/lib/test/data/2chipA.fits.gz b/lib/test/data/2chipA.fits.gz Binary files differdeleted file mode 100644 index 27a68ab..0000000 --- a/lib/test/data/2chipA.fits.gz +++ /dev/null diff --git a/lib/test/data/2chipB.fits.gz b/lib/test/data/2chipB.fits.gz Binary files differdeleted file mode 100644 index 209d5f9..0000000 --- a/lib/test/data/2chipB.fits.gz +++ /dev/null diff --git a/lib/test/data/2chipC.fits.gz b/lib/test/data/2chipC.fits.gz Binary files differdeleted file mode 100644 index 10af627..0000000 --- a/lib/test/data/2chipC.fits.gz +++ /dev/null diff --git a/lib/test/test.py b/lib/test/test.py deleted file mode 100644 index b551090..0000000 --- a/lib/test/test.py +++ /dev/null @@ -1,277 +0,0 @@ -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/lib/test/test_intersection.py b/lib/test/test_intersection.py deleted file mode 100644 index 6dacb3d..0000000 --- a/lib/test/test_intersection.py +++ /dev/null @@ -1,189 +0,0 @@ -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/lib/test/test_shared.py b/lib/test/test_shared.py deleted file mode 100644 index 0853d3a..0000000 --- a/lib/test/test_shared.py +++ /dev/null @@ -1,12 +0,0 @@ -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/lib/test/test_skyline.py b/lib/test/test_skyline.py deleted file mode 100644 index 929f0ba..0000000 --- a/lib/test/test_skyline.py +++ /dev/null @@ -1,225 +0,0 @@ -""" -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/lib/test/test_union.py b/lib/test/test_union.py deleted file mode 100644 index 4b649f6..0000000 --- a/lib/test/test_union.py +++ /dev/null @@ -1,215 +0,0 @@ -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/lib/test/test_util.py b/lib/test/test_util.py deleted file mode 100644 index 6dc0393..0000000 --- a/lib/test/test_util.py +++ /dev/null @@ -1,14 +0,0 @@ -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/lib/vector.py b/lib/vector.py deleted file mode 100644 index c6887c4..0000000 --- a/lib/vector.py +++ /dev/null @@ -1,243 +0,0 @@ -# -*- 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 |