From af99439c057871ea8974f5e4cdd674b8c6110dfe Mon Sep 17 00:00:00 2001 From: lim Date: Thu, 29 Mar 2012 15:31:13 +0000 Subject: lim added files to sphere branch git-svn-id: http://svn.stsci.edu/svn/ssb/stsci_python/stsci_python/branches/sphere@15878 fe389314-cf27-0410-b35b-8c050e845b92 Former-commit-id: b51412b5b545598ae101152eeaed470b08616176 --- lib/__init__.py | 0 lib/graph.py | 792 +++++++++++++++++++++++++++++++ lib/great_circle_arc.py | 318 +++++++++++++ lib/polygon.py | 679 ++++++++++++++++++++++++++ lib/skyline.py | 48 ++ lib/test/__init__.py | 1 + lib/test/benchmarks.py | 39 ++ lib/test/data/1904-66_TAN.fits | 375 +++++++++++++++ lib/test/data/2chipA.fits.REMOVED.git-id | 1 + lib/test/data/2chipB.fits.REMOVED.git-id | 1 + lib/test/data/2chipC.fits.REMOVED.git-id | 1 + lib/test/test.py | 260 ++++++++++ lib/test/test_intersection.py | 167 +++++++ lib/test/test_union.py | 169 +++++++ lib/test/test_util.py | 14 + lib/vector.py | 213 +++++++++ 16 files changed, 3078 insertions(+) create mode 100644 lib/__init__.py create mode 100644 lib/graph.py create mode 100644 lib/great_circle_arc.py create mode 100644 lib/polygon.py create mode 100644 lib/skyline.py create mode 100644 lib/test/__init__.py create mode 100644 lib/test/benchmarks.py create mode 100644 lib/test/data/1904-66_TAN.fits create mode 100644 lib/test/data/2chipA.fits.REMOVED.git-id create mode 100644 lib/test/data/2chipB.fits.REMOVED.git-id create mode 100644 lib/test/data/2chipC.fits.REMOVED.git-id create mode 100644 lib/test/test.py create mode 100644 lib/test/test_intersection.py create mode 100644 lib/test/test_union.py create mode 100644 lib/test/test_util.py create mode 100644 lib/vector.py (limited to 'lib') diff --git a/lib/__init__.py b/lib/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/lib/graph.py b/lib/graph.py new file mode 100644 index 0000000..e9f03ba --- /dev/null +++ b/lib/graph.py @@ -0,0 +1,792 @@ +# -*- coding: utf-8 -*- + +# Copyright (C) 2011 Association of Universities for Research in +# Astronomy (AURA) +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions +# are met: +# +# 1. Redistributions of source code must retain the above +# copyright notice, this list of conditions and the following +# disclaimer. +# +# 2. Redistributions in binary form must reproduce the above +# copyright notice, this list of conditions and the following +# disclaimer in the documentation and/or other materials +# provided with the distribution. +# +# 3. The name of AURA and its representatives may not be used to +# endorse or promote products derived from this software without +# specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY AURA ``AS IS'' AND ANY EXPRESS OR +# IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL AURA BE LIABLE FOR ANY DIRECT, +# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, +# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED +# OF THE POSSIBILITY OF SUCH DAMAGE. + +""" +This contains the code that does the actual unioning of regions. +""" +# TODO: Weak references for memory management problems? + +# STDLIB +import itertools +import weakref + +# THIRD-PARTY +import numpy as np + +# LOCAL +#from . import great_circle_arc +#from . import vector +import great_circle_arc +import vector + +# Set to True to enable some sanity checks +DEBUG = False + + +class Graph: + """ + A graph of nodes connected by edges. The graph is used to build + unions between polygons. + + .. note:: + This class is not meant to be used directly. Instead, use + `sphere.polygon.SphericalPolygon.union` and + `sphere.polygon.SphericalPolygon.intersection`. + """ + + class Node: + """ + A `Node` represents a single point, connected by an arbitrary + number of `Edge` objects to other `Node` objects. + """ + def __init__(self, point, source_polygons=[]): + """ + Parameters + ---------- + point : 3-sequence (*x*, *y*, *z*) coordinate + + source_polygon : `~sphere.polygon.SphericalPolygon` instance, optional + The polygon(s) this node came from. Used for bookkeeping. + """ + self._point = np.asanyarray(point) + self._source_polygons = set(source_polygons) + self._edges = weakref.WeakSet() + + def __repr__(self): + return "Node(%s %d)" % (str(self._point), len(self._edges)) + + def follow(self, edge): + """ + Follows from one edge to another across this node. + + Parameters + ---------- + edge : `Edge` instance + The edge to follow away from. + + Returns + ------- + other : `Edge` instance + The other edge. + """ + edges = list(self._edges) + assert len(edges) == 2 + try: + return edges[not edges.index(edge)] + except IndexError: + raise ValueError("Following from disconnected edge") + + def equals(self, other): + """ + Returns `True` if the location of this and the *other* + `Node` are the same. + """ + return np.array_equal(self._point, other._point) + + + class Edge: + """ + An `Edge` represents a connection between exactly two `Node` + objects. This `Edge` class has no direction. + """ + def __init__(self, A, B, source_polygons=[]): + """ + Parameters + ---------- + A, B : `Node` instances + + source_polygon : `~sphere.polygon.SphericalPolygon` instance, optional + The polygon this edge came from. Used for bookkeeping. + """ + self._nodes = [A, B] + for node in self._nodes: + node._edges.add(self) + self._source_polygons = set(source_polygons) + + def __repr__(self): + nodes = self._nodes + return "Edge(%s -> %s)" % (nodes[0]._point, nodes[1]._point) + + def follow(self, node): + """ + Follow along the edge from the given *node* to the other + node. + + Parameters + ---------- + node : `Node` instance + + Returns + ------- + other : `Node` instance + """ + nodes = self._nodes + try: + return nodes[not nodes.index(node)] + except IndexError: + raise ValueError("Following from disconnected node") + + + def __init__(self, polygons): + """ + Parameters + ---------- + polygons : sequence of `~sphere.polygon.SphericalPolygon` instances + Build a graph from this initial set of polygons. + """ + self._nodes = set() + self._edges = set() + self._source_polygons = set() + + self.add_polygons(polygons) + + def add_polygons(self, polygons): + """ + Add more polygons to the graph. + + .. note:: + Must be called before `union` or `intersection`. + + Parameters + ---------- + polygons : sequence of `~sphere.polygon.SphericalPolygon` instances + Set of polygons to add to the graph + """ + for polygon in polygons: + self.add_polygon(polygon) + + def add_polygon(self, polygon): + """ + Add a single polygon to the graph. + + .. note:: + Must be called before `union` or `intersection`. + + Parameters + ---------- + polygon : `~sphere.polygon.SphericalPolygon` instance + Polygon to add to the graph + """ + points = polygon._points + + if len(points) < 3: + raise ValueError("Too few points in polygon") + + self._source_polygons.add(polygon) + + start_node = nodeA = self.add_node(points[0], [polygon]) + for i in range(1, len(points) - 1): + # Don't create self-pointing edges + if not np.array_equal(points[i], nodeA._point): + nodeB = self.add_node(points[i], [polygon]) + self.add_edge(nodeA, nodeB, [polygon]) + nodeA = nodeB + # Close the polygon + self.add_edge(nodeA, start_node, [polygon]) + + def add_node(self, point, source_polygons=[]): + """ + Add a node to the graph. It will be disconnected until used + in a call to `add_edge`. + + Parameters + ---------- + point : 3-sequence (*x*, *y*, *z*) coordinate + + source_polygon : `~sphere.polygon.SphericalPolygon` instance, optional + The polygon this node came from. Used for bookkeeping. + + Returns + ------- + node : `Node` instance + The new node + """ + node = self.Node(point, source_polygons) + self._nodes.add(node) + return node + + def remove_node(self, node): + """ + Removes a node and all of the edges that touch it. + + .. note:: + It is assumed that *Node* is already a part of the graph. + + Parameters + ---------- + node : `Node` instance + """ + assert node in self._nodes + + for edge in list(node._edges): + nodeB = edge.follow(node) + nodeB._edges.remove(edge) + self._edges.remove(edge) + self._nodes.remove(node) + + def add_edge(self, A, B, source_polygons=[]): + """ + Add an edge between two nodes. + + .. note:: + It is assumed both nodes already belong to the graph. + + Parameters + ---------- + A, B : `Node` instances + + source_polygons : `~sphere.polygon.SphericalPolygon` instance, optional + The polygon(s) this edge came from. Used for bookkeeping. + + Returns + ------- + edge : `Edge` instance + The new edge + """ + assert A in self._nodes + assert B in self._nodes + + edge = self.Edge(A, B, source_polygons) + self._edges.add(edge) + return edge + + def remove_edge(self, edge): + """ + Remove an edge from the graph. The nodes it points to remain intact. + + .. note:: + It is assumed that *edge* is already a part of the graph. + + Parameters + ---------- + edge : `Edge` instance + """ + assert edge in self._edges + + A, B = edge._nodes + A._edges.remove(edge) + if len(A._edges) == 0: + self.remove_node(A) + if A is not B: + B._edges.remove(edge) + if len(B._edges) == 0: + self.remove_node(B) + self._edges.remove(edge) + + def split_edge(self, edge, node): + """ + Splits an `Edge` *edge* at `Node` *node*, removing *edge* and + replacing it with two new `Edge` instances. It is intended + that *E* is along the original edge, but that is not enforced. + + Parameters + ---------- + edge : `Edge` instance + The edge to split + + node : `Node` instance + The node to insert + + Returns + ------- + edgeA, edgeB : `Edge` instances + The two new edges on either side of *node*. + """ + assert edge in self._edges + assert node in self._nodes + + A, B = edge._nodes + edgeA = self.add_edge(A, node, edge._source_polygons) + edgeB = self.add_edge(node, B, edge._source_polygons) + self.remove_edge(edge) + return [edgeA, edgeB] + + def _sanity_check(self, node_is_2=False): + """ + For debugging purposes: assert that edges and nodes are + connected to each other correctly and there are no orphaned + edges or nodes. + """ + if not DEBUG: + return + + try: + unique_edges = set() + for edge in self._edges: + for node in edge._nodes: + assert edge in node._edges + assert node in self._nodes + edge_repr = [tuple(x._point) for x in edge._nodes] + edge_repr.sort() + edge_repr = tuple(edge_repr) + # assert edge_repr not in unique_edges + unique_edges.add(edge_repr) + + for node in self._nodes: + if node_is_2: + assert len(node._edges) == 2 + else: + assert len(node._edges) >= 2 + for edge in node._edges: + assert node in edge._nodes + assert edge in self._edges + except AssertionError as e: + import traceback + traceback.print_exc() + self._dump_graph() + raise + + def _dump_graph(self, title=None, lon_0=0, lat_0=90, + projection='ortho', func=lambda x: len(x._edges)): + from mpl_toolkits.basemap import Basemap + from matplotlib import pyplot as plt + fig = plt.figure() + m = Basemap(projection="ortho", + lon_0=lon_0, + lat_0=lat_0) + m.drawparallels(np.arange(-90., 90., 20.)) + m.drawmeridians(np.arange(0., 420., 20.)) + m.drawmapboundary(fill_color='white') + + counts = {} + for node in self._nodes: + count = func(node) + counts.setdefault(count, []) + counts[count].append(list(node._point)) + + for k, v in counts.items(): + v = np.array(v) + ra, dec = vector.vector_to_radec(v[:, 0], v[:, 1], v[:, 2]) + x, y = m(ra, dec) + m.plot(x, y, 'o', label=str(k)) + + for edge in list(self._edges): + A, B = [x._point for x in edge._nodes] + r0, d0 = vector.vector_to_radec(A[0], A[1], A[2]) + r1, d1 = vector.vector_to_radec(B[0], B[1], B[2]) + m.drawgreatcircle(r0, d0, r1, d1, color='blue') + + if title: + plt.title(title) + plt.legend() + plt.show() + + def union(self): + """ + Once all of the polygons have been added to the graph, + join the polygons together. + + Returns + ------- + points : Nx3 array of (*x*, *y*, *z*) points + This is a list of points outlining the union of the + polygons that were given to the constructor. If the + original polygons are disjunct or contain holes, cut lines + will be included in the output. + """ + self._remove_cut_lines() + self._sanity_check() + self._find_all_intersections() + self._sanity_check() + self._remove_interior_nodes() + self._sanity_check() + self._remove_3ary_edges() + self._sanity_check(True) + return self._trace() + + def intersection(self): + """ + Once all of the polygons have been added to the graph, + calculate the intersection. + + Returns + ------- + points : Nx3 array of (*x*, *y*, *z*) points + This is a list of points outlining the intersection of the + polygons that were given to the constructor. If the + resulting polygons are disjunct or contain holes, cut lines + will be included in the output. + """ + self._remove_cut_lines() + self._sanity_check() + self._find_all_intersections() + self._sanity_check() + self._remove_exterior_edges() + self._sanity_check() + self._remove_3ary_edges(large_first=True) + self._sanity_check(True) + return self._trace() + + def _remove_cut_lines(self): + """ + Removes any cutlines that may already have existed in the + input polygons. This is so any cutlines in the final result + will be optimized to be as short as possible and won't + intersect each other. + + This works by finding coincident edges that are reverse to + each other, and then splicing around them. + """ + # As this proceeds, edges are removed from the graph. It + # iterates over a static list of all edges that exist at the + # start, so each time one is selected, we need to ensure it + # still exists as part of the graph. + + # This transforms the following (where = is the cut line) + # + # \ / + # A' + + B' + # | | + # A +====================+ B + # + # D +====================+ C + # | | + # D' + + C' + # / \ + # + # to this: + # + # \ / + # A' + + B' + # | | + # A + + C + # | | + # D' + + C' + # / \ + # + + edges = list(self._edges) + for i in xrange(len(edges)): + AB = edges[i] + if AB not in self._edges: + continue + A, B = AB._nodes + for j in xrange(i + 1, len(edges)): + CD = edges[j] + if CD not in self._edges: + continue + C, D = CD._nodes + # To be a cutline, the candidate edges need to run in + # the opposite direction, hence A == D and B == C, not + # A == C and B == D. + if (A.equals(D) and B.equals(C)): + # Create new edges A -> D' and C -> B' + self.add_edge( + A, D.follow(CD).follow(D), + AB._source_polygons | CD._source_polygons) + self.add_edge( + C, B.follow(AB).follow(B), + AB._source_polygons | CD._source_polygons) + + # Remove B and D which are identical to C and A + # respectively. We do not need to remove AB and + # CD because this will remove it for us. + self.remove_node(D) + self.remove_node(B) + + break + + def _find_all_intersections(self): + """ + Find all the intersecting edges in the graph. For each + intersecting pair, four new edges are created around the + intersection point. + """ + def get_edge_points(edges): + return (np.array([x._nodes[0]._point for x in edges]), + np.array([x._nodes[1]._point for x in edges])) + + # For speed, we want to vectorize all of the intersection + # calculations. Therefore, there is a list of edges, and two + # arrays containing the end points of those edges. They all + # need to have things added and removed from them at the same + # time to keep them in sync, but of course the interface for + # doing so is different between Python lists and numpy arrays. + + edges = list(self._edges) + starts, ends = get_edge_points(edges) + + while len(edges) > 1: + AB = edges.pop(0) + A = starts[0]; starts = starts[1:] # numpy equiv of "pop(0)" + B = ends[0]; ends = ends[1:] # numpy equiv of "pop(0)" + + # Calculate the intersection points between AB and all + # other remaining edges + with np.errstate(invalid='ignore'): + intersections = great_circle_arc.intersection( + A, B, starts, ends) + # intersects is `True` everywhere intersections has an + # actual intersection + intersects = np.isfinite(intersections[..., 0]) + + intersection_indices = np.nonzero(intersects)[0] + + # Iterate through the candidate intersections, if any -- + # we want to eliminate intersections that only intersect + # at the end points + for j in intersection_indices: + C = starts[j] + D = ends[j] + if (np.array_equal(C, A) or np.array_equal(C, B) or + np.array_equal(D, A) or np.array_equal(D, B)): + continue + + CD = edges[j] + E = intersections[j] + + # This is a bona-fide intersection, and E is the + # point at which the two lines intersect. Make a + # new node for it -- this must belong to the all + # of the source polygons of both of the edges that + # crossed. + + # A + # | + # C--E--D + # | + # B + + E = self.add_node( + E, AB._source_polygons | CD._source_polygons) + newA, newB = self.split_edge(AB, E) + newC, newD = self.split_edge(CD, E) + + new_edges = [newA, newB, newC, newD] + + # Delete CD, and push the new edges to the + # front so they will be tested for intersection + # against all remaining edges. + del edges[j] # CD + edges = new_edges + edges + new_starts, new_ends = get_edge_points(new_edges) + starts = np.vstack( + (new_starts, starts[:j], starts[j+1:])) + ends = np.vstack( + (new_ends, ends[:j], ends[j+1:])) + break + + def _remove_interior_nodes(self): + """ + Removes any nodes that are contained inside other polygons. + What's left is the (possibly disjunct) outline. + """ + polygons = self._source_polygons + + # node._count is the number of polygons that the node is + # inside, other than the polygon that it came from + for node in self._nodes: + node._count = 0 + for polygon in polygons: + if polygon in node._source_polygons: + continue + if polygon.contains_point(node._point): + node._count += 1 + + for node in list(self._nodes): + # Nodes with a count >= 1 are completely contained in the + # outer polygon, so they should be removed. What's left + # are only outer nodes. + if node._count >= 1: + self.remove_node(node) + + def _remove_exterior_edges(self): + """ + Removes any edges that are not contained in all of the source + polygons. What's left is the (possibly disjunct) outline. + """ + polygons = self._source_polygons + + for edge in self._edges: + edge._count = 0 + for polygon in polygons: + if (polygon in edge._source_polygons or + polygon.intersects_arc( + edge._nodes[0]._point, edge._nodes[1]._point)): + edge._count += 1 + + for edge in list(self._edges): + if edge._count < len(polygons): + self.remove_edge(edge) + + def _remove_3ary_edges(self, large_first=False): + """ + Remove edges between pairs of nodes that have more than 3 + edges. This removes triangles that can't be traced. + """ + if large_first: + max_ary = 0 + for node in self._nodes: + max_ary = max(len(node._edges), max_ary) + order = range(max_ary + 1, 2, -1) + else: + order = [3] + + for i in order: + removals = [] + for edge in list(self._edges): + if (len(edge._nodes[0]._edges) >= i and + len(edge._nodes[1]._edges) >= i): + removals.append(edge) + + for edge in removals: + if edge in self._edges: + self.remove_edge(edge) + + def _trace(self): + """ + Given a graph that has had cutlines removed and all + intersections found, traces it to find a resulting single + polygon. + """ + polygons = [] + edges = set(self._edges) # copy + seen_nodes = set() + while len(edges): + polygon = [] + edge = edges.pop() + start_node = node = edge._nodes[0] + while True: + # TODO: Do we need this if clause any more? + if len(polygon): + if not np.array_equal(polygon[-1], node._point): + polygon.append(node._point) + else: + polygon.append(node._point) + edge = node.follow(edge) + edges.discard(edge) + node = edge.follow(node) + if node is start_node: + polygon.append(node._point) + break + + polygons.append(np.asarray(polygon)) + + if len(polygons) == 1: + return polygons[0] + elif len(polygons) == 0: + return [] + else: + return self._join(polygons) + + def _join(self, polygons): + """ + If the graph is disjunct, joins the parts with cutlines. + + The closest nodes between each pair that don't intersect + any other edges are used as cutlines. + + TODO: This is not optimal, because the closest distance + between two polygons may not in fact be between two vertices, + but may be somewhere along an edge. + """ + def do_join(polygons): + all_polygons = polygons[:] + + skipped = 0 + + polyA = polygons.pop(0) + while len(polygons): + polyB = polygons.pop(0) + + # If fewer than 3 edges, it's not a polygon, + # just throw it out + if len(polyB) < 4: + continue + + # Find the closest set of vertices between polyA and + # polyB that don't cross any of the edges in *any* of + # the polygons + closest = np.inf + closest_pair_idx = (None, None) + for a in xrange(len(polyA) - 1): + A = polyA[a] + distances = great_circle_arc.length(A, polyB[:-1]) + b = np.argmin(distances) + distance = distances[b] + if distance < closest: + B = polyB[b] + # Does this candidate line cross other edges? + crosses = False + for poly in all_polygons: + if np.any( + great_circle_arc.intersects( + A, B, poly[:-1], poly[1:])): + crosses = True + break + if not crosses: + closest = distance + closest_pair_idx = (a, b) + + if not np.isfinite(closest): + # We didn't find a pair of points that don't cross + # something else, so we want to try to join another + # polygon. Defer the current polygon until later. + if len(polygons) in (0, skipped): + return None + polygons.append(polyB) + skipped += 1 + else: + # Splice the two polygons together using a cut + # line + a, b = closest_pair_idx + new_poly = np.vstack(( + # polyA up to and including the cut point + polyA[:a+1], + + # polyB starting with the cut point and + # wrapping around back to the cut point. + # Ignore the last point in polyB, because it + # is the same as the first + np.roll(polyB[:-1], -b, 0), + + # The cut point on polyB + polyB[b:b+1], + + # the rest of polyA, starting with the cut + # point + polyA[a:] + )) + + skipped = 0 + polyA = new_poly + + return polyA + + for permutation in itertools.permutations(polygons): + poly = do_join(list(permutation)) + if poly is not None: + return poly + + raise RuntimeError("Could not find cut points") diff --git a/lib/great_circle_arc.py b/lib/great_circle_arc.py new file mode 100644 index 0000000..0fff174 --- /dev/null +++ b/lib/great_circle_arc.py @@ -0,0 +1,318 @@ +# -*- 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 + +# THIRD-PARTY +import numpy as np + + +try: + #from . import math_util + import math_util + HAS_C_UFUNCS = True +except ImportError: + HAS_C_UFUNCS = False + + +__all__ = ['angle', 'intersection', 'intersects', 'length', 'midpoint'] + + +if not HAS_C_UFUNCS: + from numpy.core.umath_tests import inner1d + + 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. + """ + 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): + 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): + ur""" + 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 + `_ + by Roger Stafford. + + http://www.mathworks.com/matlabcentral/newsreader/view_thread/276271 + """ + 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) +else: + inner1d = math_util.inner1d + _fast_cross = math_util.cross + _cross_and_normalize = math_util.cross_and_norm + intersection = math_util.intersection + + +def length(A, B, degrees=True): + ur""" + 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 · B) + """ + A = np.asanyarray(A) + B = np.asanyarray(B) + + dot = inner1d(A, B) + 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. + + 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 diff --git a/lib/polygon.py b/lib/polygon.py new file mode 100644 index 0000000..b9b5d8c --- /dev/null +++ b/lib/polygon.py @@ -0,0 +1,679 @@ +# -*- 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 + +# STDLIB +from copy import copy + +# THIRD-PARTY +import numpy as np + +# LOCAL +#from . import great_circle_arc +#from . import vector +import great_circle_arc +import vector + +__all__ = ['SphericalPolygon'] + + +class SphericalPolygon(object): + ur""" + 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): + ur""" + 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 + This point must be inside the polygon. + """ + if len(points) == 0: + pass + 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) + self._inside = np.asanyarray(inside) + + # TODO: Detect self-intersection and fix + + def __copy__(self): + return self.__class__(np.copy(self._points), np.copy(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 + + @classmethod + def from_radec(cls, ra, dec, center=None, degrees=True): + ur""" + 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): + ur""" + 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 + 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): + ur""" + 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 contains_point(self, point): + ur""" + 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): + ur""" + 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): + ur""" + 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 based on the sum of the radian angles: + + .. math:: + + S = \sum^n_{i=0} \theta_i - (n - 2) \pi + """ + if len(self._points) < 3: + return np.array(0.0) + + A = self._points[:-1] + B = np.roll(A, 1, 0) + C = np.roll(B, 1, 0) + + angles = great_circle_arc.angle(A, B, C, degrees=False) + + midpoints = great_circle_arc.midpoint(A, C) + interior = [self.contains_point(x) for x in midpoints] + angles = np.where(interior, angles, 2.*np.pi - angles) + + sum = np.sum(angles) + area = sum - float(len(angles) - 2) * np.pi + + return area + + 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. + """ + import graph + g = graph.Graph([self, other]) + + polygon = g.union() + + return self.__class__(polygon, self._inside) + + @classmethod + def multi_union(cls, polygons, method='parallel'): + """ + Return a new `SphericalPolygon` that is the union 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 unions: + + - 'parallel' (default): A graph is built using all of + the polygons, and the union operation is computed on + the entire thing globally. + + - 'serial': The polygon is built in steps by adding one + polygon at a time and unioning 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 + + See also + -------- + union + """ + assert len(polygons) + for polygon in polygons: + assert isinstance(polygon, SphericalPolygon) + + import graph + + if method.lower() == 'parallel': + g = graph.Graph(polygons) + polygon = g.union() + return cls(polygon, polygons[0]._inside) + elif method.lower() == 'serial': + result = copy(polygons[0]) + for polygon in polygons[1:]: + result = result.union(polygon) + return result + else: + raise ValueError("method must be 'parallel' or 'serial'") + + @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]) + + import graph + 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]) + + 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) + return result + else: + raise ValueError("method must be 'parallel' or 'serial'") + + def overlap(self, other): + ur""" + 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 + ra, dec = vector.vector_to_radec( + points[:, 0], points[:, 1], points[:, 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/skyline.py b/lib/skyline.py new file mode 100644 index 0000000..05fd279 --- /dev/null +++ b/lib/skyline.py @@ -0,0 +1,48 @@ +""" +skyline -- Manage outlines on the sky + +This module provides support for working with footprints on the sky. +Primary use case would use the following generalized steps:: + +1. Initialize SkyLine objects for each input image. This object would be the + union of all the input image's individual chips WCS footprints. +2. Determine overlap between all images. The determination would employ a + recursive operation to return the extended list of all overlap values + computed as [img1 vs [img2,img3,...,imgN],img2 vs [img3,...,imgN],...] +3. Select the pair with the largest overlap, or the pair which produces the + largest overlap with the first input image. This defines the initial + reference SkyLine object. +4. Perform some operation on the 2 images: for example, match sky in intersecting + regions, or aligning second image with the first (reference) image. +5. Update the second image, either apply the sky value or correct the WCS, then + generate a new SkyLine object for that image. +6. Create a new reference SkyLine object as the union of the initial reference + object and the newly updated SkyLine object. +7. Repeat Steps 2-6 for all remaining input images. + +This process will work reasonably fast as most operations are performed using +the SkyLine objects and WCS information solely, not image data itself. +""" +import pyfits + +import sphere + +class SkyLine(object): + def __init__(self,fname): + pass + + def add_image(self,skyline): + pass + + def compute_overlap(self, skyline): + pass + + def find_intersection(self, skyline): + pass + + def within(self, pos): + pass + + def create_wcs(self): + pass + diff --git a/lib/test/__init__.py b/lib/test/__init__.py new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/lib/test/__init__.py @@ -0,0 +1 @@ + diff --git a/lib/test/benchmarks.py b/lib/test/benchmarks.py new file mode 100644 index 0000000..9711d7b --- /dev/null +++ b/lib/test/benchmarks.py @@ -0,0 +1,39 @@ +import os +import sys +import time + +import numpy as np +from sphere import * +from test_util import * + +def point_in_poly_lots(): + poly1 = SphericalPolygon.from_wcs( + os.path.join(ROOT_DIR, '1904-66_TAN.fits'), 64, crval=[0, 87]) + poly2 = SphericalPolygon.from_wcs( + os.path.join(ROOT_DIR, '1904-66_TAN.fits'), 64, crval=[20, 89]) + poly3 = SphericalPolygon.from_wcs( + os.path.join(ROOT_DIR, '1904-66_TAN.fits'), 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 b/lib/test/data/1904-66_TAN.fits new file mode 100644 index 0000000..09c3929 --- /dev/null +++ b/lib/test/data/1904-66_TAN.fits @@ -0,0 +1,375 @@ +SIMPLE = T BITPIX = -32 / IEEE (big-endian) 32-bit floating point data NAXIS = 2 NAXIS1 = 192 NAXIS2 = 192 BUNIT = 'JY/BEAM ' CTYPE1 = 'RA---TAN' CRPIX1 = -2.680658087122E+02 CDELT1 = -6.666666666667E-02 CRVAL1 = 0.000000000000E+00 CTYPE2 = 'DEC--TAN' CRPIX2 = -5.630437201085E-01 CDELT2 = 6.666666666667E-02 CRVAL2 = -9.000000000000E+01 LONPOLE = 1.800000000000E+02 / Native longitude of celestial pole LATPOLE = -9.000000000000E+01 / Native latitude of celestial pole EQUINOX = 2.000000000000E+03 / Equinox of equatorial coordinates BMAJ = 2.399999936422E-01 / Beam major axis in degrees BMIN = 2.399999936422E-01 / Beam minor axis in degrees BPA = 0.000000000000E+00 / Beam position angle in degrees RESTFRQ = 1.420405750000E+09 / Line rest frequency, Hz HISTORY Parkes Multibeam continuum map HISTORY Formed on Mon 2004/02/09 01:23:37 GMT by "pksgridzilla" which was HISTORY compiled on Feb 9 2004 12:08:02 (local time) within HISTORY AIPS++ version 19.405.00 dated . HISTORY Polarization mode: A and B aggregated HISTORY Gridding parameters: HISTORY Method: WGTMED HISTORY Clip fraction: 0.000 HISTORY Tsys weighting: applied HISTORY Beam weight order: 1 HISTORY Beam FWHM: 14.4 arcmin HISTORY Beam normalization: applied HISTORY Smoothing kernel type: TOP-HAT HISTORY Kernel FWHM: 12.0 arcmin HISTORY Cutoff radius: 6.0 arcmin HISTORY Beam RSS cutoff: 0.0 HISTORY Input data sets: HISTORY 97-10-09_0356_193558-66_206a.sdfits HISTORY 97-10-12_0142_182123-66_193a.sdfits HISTORY 97-10-12_0151_182707-66_194a.sdfits HISTORY 97-10-12_0200_183252-66_195a.sdfits HISTORY 97-11-07_0510_183836-66_196a.sdfits HISTORY 97-11-07_0519_184420-66_197a.sdfits HISTORY 97-11-07_0528_185004-66_198a.sdfits HISTORY 97-11-07_0537_185548-66_199a.sdfits HISTORY 97-11-07_0546_190132-66_200a.sdfits HISTORY 97-11-07_0556_190717-66_201a.sdfits HISTORY 97-11-07_0645_191301-66_202a.sdfits HISTORY 97-11-07_0654_191845-66_203a.sdfits HISTORY 97-11-07_0703_192429-66_204a.sdfits HISTORY 97-11-07_0712_193013-66_205a.sdfits HISTORY 97-11-07_0724_194142-66_207a.sdfits HISTORY 97-11-18_0256_193815-66_206c.sdfits HISTORY 97-11-18_0306_194359-66_207c.sdfits HISTORY 97-11-19_0447_182341-66_193c.sdfits HISTORY 97-11-19_0456_182925-66_194c.sdfits HISTORY 97-11-19_0507_190350-66_200c.sdfits HISTORY 97-11-19_0516_190934-66_201c.sdfits HISTORY 97-11-19_0525_191519-66_202c.sdfits HISTORY 97-11-19_0534_192103-66_203c.sdfits HISTORY 97-11-19_0544_192647-66_204c.sdfits HISTORY 97-11-19_0553_193231-66_205c.sdfits HISTORY 97-11-19_0602_183509-66_195c.sdfits HISTORY 97-11-19_0612_184053-66_196c.sdfits HISTORY 97-11-19_0622_184638-66_197c.sdfits HISTORY 97-11-19_0631_185222-66_198c.sdfits HISTORY 97-11-19_0640_185806-66_199c.sdfits HISTORY 98-03-24_2107_193706-66_206b.sdfits HISTORY 98-03-24_2116_194251-66_207b.sdfits HISTORY 98-03-25_2020_190826-66_201b.sdfits HISTORY 98-03-25_2029_191410-66_202b.sdfits HISTORY 98-03-25_2038_191954-66_203b.sdfits HISTORY 98-03-25_2047_192538-66_204b.sdfits HISTORY 98-03-25_2056_193122-66_205b.sdfits HISTORY 98-03-26_2048_190459-66_200d.sdfits HISTORY 98-03-27_2034_191627-66_202d.sdfits HISTORY 98-03-27_2043_192212-66_203d.sdfits HISTORY 98-03-27_2052_192756-66_204d.sdfits HISTORY 98-03-27_2102_193340-66_205d.sdfits HISTORY 98-03-27_2111_193924-66_206d.sdfits HISTORY 98-03-27_2120_194508-66_207d.sdfits HISTORY 98-03-27_2130_191043-66_201d.sdfits HISTORY 98-05-10_2123_182232-66_193b.sdfits HISTORY 98-05-10_2133_182816-66_194b.sdfits HISTORY 98-05-10_2142_183400-66_195b.sdfits HISTORY 98-05-10_2151_183945-66_196b.sdfits HISTORY 98-05-10_2200_184529-66_197b.sdfits HISTORY 98-05-10_2209_185113-66_198b.sdfits HISTORY 98-05-10_2219_185657-66_199b.sdfits HISTORY 98-05-10_2228_190241-66_200b.sdfits HISTORY 98-05-13_2132_182450-66_193d.sdfits HISTORY 98-05-13_2151_183034-66_194d.sdfits HISTORY 98-05-13_2200_183618-66_195d.sdfits HISTORY 98-05-13_2210_184202-66_196d.sdfits HISTORY 98-05-13_2219_184746-66_197d.sdfits HISTORY 98-05-13_2228_185331-66_198d.sdfits HISTORY 98-05-13_2237_185915-66_199d.sdfits HISTORY 98-05-25_1711_182559-66_193e.sdfits HISTORY 98-05-25_1720_183143-66_194e.sdfits HISTORY 98-05-25_1729_183727-66_195e.sdfits HISTORY 98-05-25_1738_184311-66_196e.sdfits HISTORY 98-05-25_1747_184855-66_197e.sdfits HISTORY 98-05-25_1756_185439-66_198e.sdfits HISTORY 98-05-25_1806_190024-66_199e.sdfits HISTORY 98-05-25_1815_190608-66_200e.sdfits HISTORY 98-05-25_1824_191152-66_201e.sdfits HISTORY 98-05-25_1833_191736-66_202e.sdfits HISTORY 98-05-25_1842_192320-66_203e.sdfits HISTORY 98-05-25_1851_192905-66_204e.sdfits HISTORY 98-05-25_1901_193449-66_205e.sdfits HISTORY 98-05-25_1910_194033-66_206e.sdfits HISTORY 98-05-25_1919_194617-66_207e.sdfits HISTORY Original FITS filename "1904-66_TAN.continuum.fits". HISTORY Noise level of continuum map: 59 mJy (RMS) END sPwľ'O$ `zģz@l-T$H(&Dƫdҽ,.޾V67rH|fJ'jվe up>'w24$U)?R-پSQ$yG:(b žos?7'M<.6]zg{p~2Ƚ}@<?c+( +żF*p,E־F/[ k>V詤Q"UXҨHp^K+F5;*nҽD`:HtCʍmaO<߽Xk>68½@ѽj1wXLϽ -`n˔Yۋм߽gμj<С< <#= +̂;5===u>/N>H>?y?.z)?,,>_>>W[$T>?18={48=1N=y< <$=]==T> +5>5>GuC>W#?S?v?K>>,=Ӽw2rT| <\uFGk=!5=zO=Jͺ=Vq=HK{|k+4~7:=`d=Q=.P=<=><=t;el>&t=@<ޠ=M=%>y> >o0>==g3==&=Gi<>j>l=Ħ=17}I= =g; *J=f~=m>}6>>9>p>??S?=> m=\t='󿽑n䱽`A;^=S: =&i=/C<ܽo7m#c oL@9&ɽ؃OV<ټdh=kb=zoi0{y6"=OJ= 9==&&-<=<;r==<Һm <=؈=>Cv>X%>c s>j>> x==>!s>ɝ>r>Zko=-&on75*$%y=N=n>=l=>90`>>>û-f=qS&P!ƽ kT弨ۚ< &=t='=S-΅>H>>b->:<=>>;>0)T>2>,j=V5<Ԭ`<'G1e1Y罖7ԽRAfv l=%XE>(>P=m=K$Ž{Qo佩ýo8+ =W_ <:AϽq{ F%ɽ?=11]=!=>Λ>,=x= == 9=p6= J +6l>>> =&H=v +=> w>E=<䗣XHą`=7=='A=>k==B酐𫽟޾"xoỼ3: ʰTｘhҫ/ֽ?Jս+:y=3=sV>p>T +D>,v> A=R9w;Jl<=]='i%zv +ȽV9pɽQAU~𽃭4rUY>R~e}:9gzO϶D;"bwμ.= W=ކ==z=+=.<D̽Fi|ݽ ==kG==s=ыL==[|<}`;g=J=C=x=t=W1=^W[=w=XZ=N=!=1:>>=l==p<<;9{:4F,LkKˈdⲽý@c0$$<\=f5<Ŝ=Vr>H>n<> =G> l==ޏ8 l>5f +jHOZ>Je> m>n==Lt==X='k=׽F_rH|EI=+p=9b<> G v$MLžDL _y,DF3G⽚`< :====D?( F=^=^<>d>Ά=R=$F<٩ ==8E=S=؅=B=^w='=@]==tZ=-=Q==Ѱ=ژ==pf= ~k;P^;p>ּg"E_Wｩ3ﺚ0]$;S>I[>kY>>Fk>>F=<,Rͽnpcm˾%+#$q<4;P9XhY@N,ɒLOE<=O?<L:f;=b=8=z <,-zbUq"> +M^>I=<ϼ&R\Fw/C ܾZƾ(IHF$WF`[6crG@kyAἢGvA7==h=ʙ= =y=lɻɻJ;r'۾JQ =}<=oɼASitzn0[pؽ>߼HSg>UN%>u>zA>a>Z?>b=> '2>&L>2^==M$>3+k%Z]<=7===놡=Y>!F=ҕ=w=r=X={K>/>Ω>eE>=l<@pLú,#Pl; )9t}C<4.v=1N- a*Vq(:ⴽN^ܽnR=a==w=јuQh,ٸü===ٝ{ +iX)&0.4a +t4ۼh=="C;};⋯9wG߫,?z-cd=c<=!뼚\=>>Bj>UJ>]>>8>1 >.jH=ę=鐏=먷=Z="=7=/$=VY<?=sn7>z=x$F=>Q<DCeY=9 F>a>=5>ngD>YG!> +[=){=>l=o:<_<-<Ic>vx>vc>3=V2=c_<,<&h{F<[|QnU&p: +/A=Y=3r=4 <'O{c{-;ASNAʾ@.je +~.m˜<=.= H=y=8;0<=v=9=@=5Y= =ϛz>x֯>>>>b->T>Q>N>A6=== +N=Gi=Y=2=hD=2 :==&3Rp>i?#z>h>U>;g>"ql=П=0= +=Ze9=>M?>Yn>y.>s2>Z>"]v=؊(<}弇:_ۮH$u _;Fcͼe\`F R{p<ѭ'=y=]> |>e>I>=ƌ; ׀Dz /AyƋ +ڽn]rS9Ge[ѽ\~Erټ61弒*T~:H#Cv ƾT"ݼ7Ǣ-= +$=L7=<66O5d>]z>p@>Ƿ>Օ?#t='>l>*> Q=>'f>>c>+>f>>Q=m>a>+=b=ԉ=>=!==͝=O;켶===IA<}=T~<==ɮ=Ӆ=Ɵ=:J="P=|<E|S>=>I͵>ǒ>>Ɩ?y>>gX>q>H+>3=j}=eX>>>>f>2>!>b>J=0옼 (7ֽ|PMG)ܼu/μF!@"Ǽ~;Z=MF(==4>>=!='m<~" +Gݽ#b'ΑP꽃u7뢽N{ƽ9ý#!,rǽ(,.ν'@f'"Q|U)aKZcC%m;Š_<:;;=5=v=D==!?<;G޽Bw@Uғ; 災&k'.$ r齰$_g|оX 'ɽSbey9q<-F<6K;Քts9 >@%8];۲n;ȫ:-{>Wu>?yw?%?%?0?1?r?z>'>&=k">=2>g>u>f9> ,=^==lq<=b=S>6=ͺ>/=]ؽW4JpA"bh%H=B=.k򜼼6=؂<. 箑>ò>V/>>>t-?b>~>r=>>N;>f=Q\=$>4>=>>>>I= +(ҽjKl)XKz8vD⽋%3u^pg;WXusykKU;=O=B> 9=M==Fy*9z7k@˽Umýq%ݽ/:(# T; rz;N&ڸ>,>t>>?+? 0>> S>k>>>9`>)> %E>0> >!%=PZ=1-e;F'}>oOGᐽF^*+]w<( TQ [{ĽDŽxF O>"a>D">= >n>>g>c+>i >|/>2>J,=~K=>>>RZ>6l>lo>JF=ؽ,H$ͮd +w ė½ޓݽ@~[G0wd:FgmCU8ǽtW0{:aJsżȼħ}5`ܽ|$.^4ǽJxk}d>(*)㟽QC7f7*)"e1 <+;:`H_EϽٽ;[1y4ܽ:-kvlN*z <6>@,>0>">To@>c\>'=1=^= +>%I>>cC>c>;=͐ +=n;=ཟ=4=#== K;EK=#'{3>Ll\=<ϻ>,t>6V>2[>1v>1=@f=aA<(=;><>v}>@0>A=?[ly +놽Ŗ6˅ƎhɃнF( īx=[NM(Z!Sý'XrGT/;>_G>r>? 91?Vl>_>[>y==>.e>>jF=]:=*=#f<-=FP]=$$U<ԧa3нѽk)]> 4>>aO>.=/Pq=>N>-t><={>==C3nj><>)?Y@5@x@@U)@/?a>ٰ>\=!;>D0>G=y=+_=b=jE=1=>b>2`}>2b2> O=/>=d=Q=m=\ =^==;F5<=P& <=E~=u==);OգmYƽ܇~!?aýS4l +?p6ս֑ } jdEYyӱ:xͼ'S8}=0*Lj-ƻXB! /#Xs hVϽ gX̼{ȼ +ͻ =4:= %=E=vO<}(?Hd0˽x0^5<@>z;B<=#*==e;dϽ扽!cڼ n<Tƽ_{h<9mO;қ<=$>k>?~@?r@^,@Ͱ@_@' @ +@6W?/>>& >Hi>J>19=Q=|AW=&=z/<;yۼUmi2gټdx;Yo_9T=c=;m>"|~<=<* X!VxJ.y'鵽p@򻸥޻mּ⽊V7JNnOһＭ5/xd<485`%;RdB<1=n==t=?>>G>q>s>6b==O1iꍑxٽI`no`\mC!^D{׽]A<= v=T=Xt&gNk +>^R>?T@Y(@!A^A7/A&Ғ@_@?>+>C<>#>9>#>*Y>$Y>0B=.=:eC,ɇνfm x +t5ht"V:DR:)8 &7S`> <>O=R=$R=P]=9;z==S=K=JE{>m?x?T?1?\*>x=jM4C˺ٽVYϽ}ltos,꽞2Dۦ<Ƽ*q\ dm\ofyK:;n9>>Ņ>f=bw<ٽLGEK轀NF,Iz>g8;rX=4=Ķ>"">j>(?yVz@'@@͞AfA-q@~W@ز?`>f>b>,=ya> s>*> t>/ܳ>*=]:nĽr૽ھ u߾2BX=>p/'eJP"پc,^>L CJؽ[`Rtm ^<n=#=Z2=?!=*=jD=SEsi+ 9.%Bݽ5#L5#Fw.:G=^e>I=9=́=8==l=ņ>>?p???p?q?"K`>W9]˼A+ƽˍ̝XzIT䢽Aʽ+@Ҡp}>6qNAPK1l#rażϠ<<)UO>֙>I>>>I?>;fB E+q)"3;,+̽L'b?7t;I<[%=K>m>?.޺?d@<@ckM@@Q@X@<@E?|>=P="ʼ;g-=[Y>5>C>.>N>*< 彅UѧF'.<[KF+9p6 uƆr_Aa<|+y6v^׽dn3~{zRg=P:=i=)=3=ȹ= p(;Y޽ eRyyry竽@QX ;Z369[1DCgw96lf~6v<B=z==ev=3=[==>>Uh?n$?3@N@?X?a;>."wLSG G W +<ɒ6!PqٽĽL[ z:|| ).荽<:J <ɤ=j +>  >+?J?X?Q?)}>?Q>?p=1<.+1Ent9󮽸 ý\7!h[q29a=^&>ׄ?d ??R@'@/A@O}@@@Q?>T=Zl=18!9ý(\1<<=sT>4N>:&> =`O &.ăwdlǠF9I3 7NqоR1v9vm +W"h?-*+Mǽ(WJU4>o>0>o)>ӌ=Ja< 7̽7>Fo;P<{G<(ҽJ;GVF +SϼA?wyir;i\Y%Ii<=F7=;A< s?[d???̟??dY>qû]佾 #}⹾)k DʽҨ3>:. Ui׽9vmI ὧzu-=h>X%>?&?ch?Zy?U>>#s=Ѿ=<*11ؽ⽡}彅k@G +Dc;}>~fO?4?ǒ(?E@@,@ Z?Z?ːr???Z>}'K$m*l=Nel====ݽtսxMC!]-=XX.!>&a><.^q3j1;|b7=8=v|=u=ۄƼ.3S>׽3b<ް/ +i`.w\ۧ A4 8y +Rl]={6>=?)?V?+?B>}>6.>R??_?1Ua?F?#D?>ˇ>]=\=.sӼ ՝F;<7޼"սǽPl;<$8>p#?;ag?U@ +@&{@#?{?r?c>?W2!>>H=(_䎣#C,d:OZ='8=ƕ=O<CD\Js1H ; p佊#q)ʾ1*]۽ 䨽ɝTne%&' #M>@s=c=^Vǫ=70x=h=&`S=Zq<һa;G;[Hvnqn,қZHz;=&=Nht4ּF<+%)J>/>ρt>>Ns>7=]>!? t?5? +۹>>>>N=h=GSսhǽg5,SjL<kB>#?.V?)~@܁@@ ??U?X> >'=8VTԼ|,ܗ½.o'Mft׼ ;d<ח8T{==(V=k=3=->u=s>&=^= .<=9Ι= ɺN_/N¼s%B#|˾Ӿ +,)aUE;=l%>>>^F>>Z>6L><ۼHм{nqGG +JC&}3]4ʽ>½S?cf1 )k}f޽jNԼf{!}<{=?>l'>@>&>Z>a->|>>>Lp>9| ߺr2~ؽHf"$Ƃ`< >mz? l)?GG?e??4%>=`=.x+q)սmG]{󇽨gH+\i@/ :z;=,<7C8b>>">-O==>=> $>tE>XR>SH>#>s==^,E=qK=n=;C:@;ťx;7 S3x۠B!pNɮs1\2guTսlRJRk؈}t8F@&ƅ52=@+>>e>>?>[9>%*=f;gm^ԽF|;b;LD( " t*: +B~d gZp|h+:;;.r;I;;+ͼjF?8<=N$>&{>Y>}>2L>>>U>r>Yz>5H=x,>$? >3>~> =W^)?3k*<ˏ~ھ&;߽Ű۽4*Ľ;/%XؼC$[p== =5>yO>U>*V>>*>> "> )2>r8=2=<==,:E|.|;C >?U>S>U0c>4' +>'>,>+>#<>>3R>,ݸ=R4=ʭ<Q Fzf 좽Cy[~_nY潿(F no'=s_>#>`= B&zlͽHȯJ; \M2*/{T8>?[:б;kJ,>!`=⛂>3>Ǝ>>g>Dj>*>>>@,=2=!<;9<,9Tnk%n[%Խp`߽5jP 9@LJӽw<<'=m=&>-G>2>=<ΒD#ZN.:ZpcWS7=ls-dzм{' <=>{(=-=:=9=Z;>====Ng==mz=Xw>==F=="4>E>>}F>7a>F>"~==Q=<}5<<w=B^>!>F> =њ==ڽ=,m>'e=?=DŽ=>?={=q=S;]7>f~hhW/ZI R 伿gQ$Ӽ0!˟ +ȼB(dt턼26.lV/Q?YֽE<|ֽ!RB]H+<D =1>=Ԙ=:;KTȼ~ݽ C)^Vn*{-g=7t +=d>(>0>=`=j=-"=mo=9=Jwe==Ax=>1>8>"=q>t=hq=#z<\ =MZ O=;*~EǶ𼍹뽥!LϾ1 @u|P࿻)l!<=p?y?@W> 9>"f>ط>p%>%=C=+==o=m2<=x +=XS_=`=A3=?= +6=¼>#♼=7\;:t];<%E<{cw" _\ɹϭ_'c;[EOP.8>a.̷MCHEˉ~X +f _bN6)1>[>8>w>yz>$=f===v=(=>D> >,H)>D>;E>3K)=;:=j< n<<]="= Ҥ<<+t]:YGb:8/YzKzoz;= +λ2<]<4<|;s_a<7Z½zνm[DApR{~⼢)]νg@r+7z<<=<8<]D$s%G.=H<ӷ?wz1J3%mV;(ʽ3d车/8)`O|7'˼qaJ<= c>(>W>NA>P>D>!>O>]m==}b="m<<": >>>b>J>t>%@3=Ҩ=Z>> >` >PJ>M,>G >):> +$Z=@ +<.;; ;= vZ˽$3HF21Y+Ƚٷ:ȼ}> >};)>y>U.>B">uJZSUܼ/W!|6)ǽG)gM/ ּ|DihӼq;=7<<;jvf;]4佑11sS>м:=W>>uX>f>*J>#>>V +>2֥>!w>8>~>/> >>?>`=Ѹ=;;][($$ʽ8)ؽ.`9/; zbzrCyO:@{T.5]1D;2^ܼQ:sm;K;%:|8;<=NNo=ANY@w>:XF=,>MrP>CF->x>2d=tn=T=C@="HNA 85/M{QoDj*onV*.4\>w>j>>1>{X>D>1>YR>U>>~>>W>8[=qSFt"tnaeX6zeĽCqOYP(z߼ѡϽ%si^@#ez 6uý6KoGD>̝'6ZvhtwBۼTۨu +%t;tEk<ɲ<#=U<$<=EO q+q)uQ<[ +Ŝ>ف==fSݸ +==]̼1 +[TA*3y7*-N<H9<=5>ej> y=x:'XhYCy ü"-2n /Q̪(KԛRpS=(oiݿiDr&; 9ۼǃUZʽU9 3°q;ͼyM;Ƹ<= =Q>!n>X>DG>>u>MS>"h>R8>) +>Q>??>~n>]=t0<}Fz49Dշ#w;@b㽼FyQ彤Ľn_"Yսms 㨽y^;!D<#<];WJ4m7輧μz11;s;]-t$6R.Pڽ7=v ( սm&>8>R>51=̧<dfrC27 =-A=a=U8=9=H|/&>a>>;>h%>M>3= +>tG>H>L>>U>;=z{;J/ D2Ahü3½IUݽ}Ž@?/Zʠa[ʽ+FI/R'iǽ4btɽ(ֈDKԼwq;8 <<6;4|S;R؁;Uw< <ijs,25yݢEs喏K*,>罤>~D^v}H'm=@>>)l?uW>{>Q=gJ>6甾#ž%*=-=M0==? =i=(<hv\"`Ӽkҽ<ý%ؽ ׽mmcYʽh!gx3B=(>>BK>>}rs>;'=N=> >Hwx>P^4>!<dE P@/DNljzĆm$ý gvuF(9޽(Ƚ):wsv;C:~l׼Pt;)c<4<<<.1= ;~ڼr޼ +vmĆ??? +>Y=C~νh&D߼6 +`>g+'Iہ9X<~= +==Q<΍;a81n +MAXG0'-;d==0b>^*>bdy>c=I=q`κ ڼ<#JG=17=ͫ>8|>v>Q>2o=y====.==yf5>>(>͢>S=*To~Qܽtg<:8$;MG<<Ҵ= J<;K:ÓI⼞'Ͻ.ݽ0,naoW=>z>>h?B>B>A+=m.G;Rt- +R;m=>>B=<3;я:=.=y³==zh<=C=36hU>1^=#"~@ێڽh FNiZ. +߻M뷿;<|!=#o?=M=&<⹙+<0(~?;?v??-P?>>T=u`T·ʽlr#2j:t 2`;=y=K=RV=a=p(nRRrZ2`> =e==].< xj(~ƌ…-Ye\ziݽ%XI,"oܔB\,<~<<ߕ<ߤ=b=OE=?P2=&=)= jq=,=TCH6|&˽4D#Ɗе8"SŘL.5_fa;:"԰X&5o۞P-,Q7ñ ' $=p/"ͻC:<}?J=Ao =W=, =&=P=:=Z'=Ɨ=I=O<{S"iJ/_a^5h;< <=;Th8Q˼iL93;-ԭ;gysB⽘S`ï<;{g>t?a?t???Z? =h4LȽ*xω?6>nl.9=g>P#>d>>H>=xDܽJcVٺ0;;9P<=7CH. %a}T6ʽ$gk\k 7!ս+m#pܺ~d;DžxN$?Q%?G?ֶ?ћ?v>?= 6 ukнٽٽϽp:W+}`=<>=> >%>>=[R'GY|)< =5n=qC=jԸT;;go9۹ʽ߼蔧tr(EY4r{5］4~f<%^F>6??1 ?}o<E,_+; -C]WJzPd﬽^K,~g<:>e>>?Y>S>=48tFl:ndg=aЃ===B#=ޗ=OU=d_=>b=g<8d;1GOŽT ۽yV D SƼY;sf; ;S@;7 >W=q=d:3&P eٺ;!};q;~;Z>:>E>>>9>N= m%nCߊR?>EA.znnqޞ=0>ar>g>Ε8>>6,7k~#EJ~F~bjR:<^v,矼"+W;"<=?RhdS Ӵ-ɾI¾˕ "轮UyǙ Bp1CbbOh+,,mҾ ޾noK1M>Czx>5 >== l}=<P<#<6[xDۻ;<R&Rʹ1ۻ=8B==l=:=#.ͽ}tGz$,ok!lKkF:f/ =0=$>RG >>> =\9h; 7==3=d=V=)U=븠=J< <qFE<RԽ>Ľ.xԫ &amJf;46@;-0<[[<Α<=Ks=W=m=>g==u;j!ϰGK>ǜk x;g;q;Hu _UXACtkS4nލLvͦ8l1һ} ;%7=+<"fO&彺&ɾ9 v [ͽ-utg@/}(ڷN:,X!g¾ϓ P ;λc=P=+}d=4=IQ=*?_=˺]_<njB:>E2>̌>>_p1>K=====_=ƽ&vԽG,='=a=e= +='=K=,V:Ps˻̻`h1'JR. +Rν"$lX>I\^q嶼5_$t<<̜d+<֡i KsC1TdOt+j弚suOr>7> +c>"$>ua>> +J>U`>\P=?=H<|8<=="r<(x{l"ׯ#;y;w!u=Ӻ=/>->a/>_9>"=$v{(>$> =R<Ӝ +BidD5;(<0*#T>r=ʯ=<8:M|яSu ӽF[1ν牽FȔr"ҽKHh @ɬP-=z%=>>C>*> έ>Q=u=W=8a=%=hA=y=sv==^=5(==w=H%=9===)> =>c>M5>]1>|>4>kV>1g=j~ }>>>>y>0=H;4< MH< =B_>>R5:>b>K=;=1Ȼơ==ENC=P=6E<¼¼NtKS`{GvaUXPL<*s:<8?w:>+]>.>zi>M>!=i={=@=a}=q֓=O=>:><=ۨ]=E===ۢ===m==|=׾=uK> z>;OS>7>+F>p=K\=*^s> ?>rD>va=7=4==n=|y>,>vJb>j>G=+<݄߼d2Of<w=iK'=`)=~;Ҙ'x( j sU^&սgcc҇V[:'6<H!!$~UJug 5 ~1Lֽ &&O_Gb,<<"ûl\h;<[;}:> +h>@%=(=s=A}C>aȟ=]=݄=䠒===cKD=J=rV>9^>>6>g>;=I=þ=qj==ٟ&= =4]=gGM==^==!=%===7=#}Z.5Oͽsܽq\˽,|rp!@:PFh:P;<}[btϼ=P >>u~>ce>ͥ>~7>Ȟ=<¨=Ӎ=l.i> +R>\>Vׂ>- +=E<~(gLZlnO2]&sv5s{E8Q齆;H ͼ7@P芼k^SIGT\d9UJmD!Q8wݼƽ0cw̲(";:; 5:ț;0<=8\= 1,g>h&8>zZ>Sl> 4!=1V=<7<6y˼dP>d>{S?> >>j= ===Β== =E= +%<ɎU=^=2h=e.=X=E=<;tQQA:;; > >>>Ey==@R<e>S==<raݶ@tv!;"MfKRSu̽aCHi꽚Z +_lv/9½9E%[S7kSM ;׻kkSCʗ;U䭝*mBP_;i<>b>u>b)e>>!Kp=>=SՇNցƼ3ýuvt̽4mqCӽ@ڼ<<*<[<]>+?7?HM?&y>[>\&=M=O=== s<2<0EݽAc<PHs&]r;/kP ,Ž6w;+/24-Te踽 CӼ s_/h +R;#: -6.lg3:zJ[XN;L;3Y|= eC=+W=>>ET3>?8>>9==+g얽6 ?0ɓ5Rὔ罖D%!½HO{A9h= +=b=HI<ͪ/\??N?Wg?:> >R%=>==6W=<`y:9;;B5ػ@ UZwTyv[мF1^M¹HvCGuozk}ڽi7׽~C!1P}[goJ@:9@"Hbb;iԺ8;=Ggda ߼^)ŝ,m1TG + D¼>7=Ub=> +L>:==>=k<@;/= ~LoE1ow[)IJ5bӽ p"B>x>?o)??>ծ>>:=B={Q#8&<;=_=<<;<Y;M~2t==_"ʽ]>N6$ӽ\Zґ+-ʼس<$V=!=^=ɺ=d=-v==^S3g6MGOCG伯@gX>7F>=>>6>*>o +>#=|<<\*P;7VF橼:39l8!27VD~L] @~PNfY",ڽ|'#kVؽϽ?[P;M@@1Vݽ^~O?ڽI/f'Cn{;ڽ<}==<=$=r<;TY$1FW:X>=B=Ȥ=i<<=c>J>Y9>$>W(>\'R> n=f;,:ݼ>ƽB4 v7IǠ4lKFZ)_]n9׽amQ o<-<9=*:vAt6˽~Iv`½"zWZ<"<2%>=^>=5>9R> |t=:*==~=~E,->h{>uk>U> Z=fv:;e)U:\KQѦM(|wUk8l,e<<<9T; ;VZHqL֪;UctȔJp +.UYҽ==>==<=\w;琽5w n c>&sD^iּܔ*racEFzTJȞڽZ彷ƥk Ⱦ|1}jr51:9<Ѳ>Yb>^h>5[==$?=m=O=#"= :á޼xx}u7#w3H3Mѽ}WՒ9ngxr?P9A?5У߽GջC==A;O=Yev==$M=/.<{:< ==œ<ܕmAp<W+h^5>>|>? =HJ&@$ x)NQk`>?J"P5Cxu9ڀp<&<=b=P=G=(5vyqeoŶ]{Q4Ǫ6 ^n}M`j_[ӽǽЊ|%3AM 5Ej<m>"M>: >4>"=Ƞ=ì|=gd<+=1 G= =:=>Ў=L 4Z>gگ>n!>#= r~K#a!!M6)G,K[kT>5 S<:̜>!jB> =-=jl<,riq뼔׼BpRJPU9!rzc^T/"@5yѽـپ,**w1^=E=x_=7b> =~><>=T>>&1> +=^="=> a=.==ݩ=S@1==焎=Fa=Խ=t='4<ޘdDi_I;HI;B<c==M< =1v==S=t9=;;V;`~<1<2<<:[tj0;<=E>$|>* O= !=q%<}ۓ1~Rp_wYTwԽwh4彋ˆN< +O%s;B$=5{==7=(= +>#}">,G>==g۽Ru=A~<f=WΈ>̏>SD>E<>94>=s>> >v>p===&[<^<~yg5 &>O>[.>Oi>>`>v=ܳ=(^=a=Z:8^Zؼ \6'-`vY7ⲽzVӼ$q;Y=$"=UM="==bʒt=ذ==J=ϧ>V>= D=<&ߓ>E|>1~>=K6<3&>=< =C!j=]=8? m4ڵqB +JYL(ֽ&⽇AԼ%<=='vX=0S=# +Q]wZKU?H4'SѮ󼑘3 P{R)0.X=Ʋ=[> >;.>2L>9=-=1==Y=.7=§=_=3L==G=Z=J]= E=kx fٶU6n#<==_==3=4b>>>=>D> N=Q=Z=U<[};ل<-^>Z >>s> P>=-T,BWڻ;8*;h;5㛼lｨ ~U4Ҿ*?B̾%(A$ҽ⮗RT1.<34=.U=0v= y$U=Aku½ iT ;̻0z<=h=R=c=ƶ=J==h=(><`^j24B ?M\5+J\ʼBkݽ0%a3%伊ٺ>\>b=.&=O|====Z>!jy>*`齛&qo><>==%=^X; ;t=#5>>ZY>~>i>> +"=K%:\Xc50Z@> *=O&=l^=h=Jy=J}=7= [>P= (>Ma2>~=>>P> u=;YIC26ƽ(=%TIH ޚ Z홄y:ؽ˚x60<=>Hg>*D>vN>X>=S=[= s=HI==$=0=]T<.;ȷ9 G)ɽ-<*_+/V¼ӗ";;;ܻM:=R>NQ<s=M=[;=\<;<"뻧;CC Yt{ѽaHEU Q +O\<ʟ̄ƻsmQt*n<~;ֻ {"ڔaa' >yCq1 O\&vKp^8 h`\-̼olq5 ,>۩>9-= ==AN&R>,>>i>>n֥> <=Ń=xٹ==Ys=w=e="=XK<;Tȼzr69&_軁!E `MO½% ռj.:i8;*;ϼ;!;< =Cn=^W=OK_]z<'u7,X;LY]N;1;p)<^'=*=~4==Y<=/7===+?=w10[Y':f> >>? t|? m >>46>4=aC=7D====Q=0M=v@?= j;*Ҽk6|ݽPPfC 5&}2@ʼ)Cc.[0:^;f<-.==F =.E=4t=cB===5=%L3<~w=\=8=ir======i<2-<<7:'Yl  𦽟'P|oq HU/@/A< <}ͻp+( ̼ٴB Z>~ 5*:̼HE/Xy䞽ϼ<vq>9>>ǚ>>K>>=#_>>>u>>O9=>=H<#;Z ;x.!̈́˼#Z\|=\Q="<_ LqIo:=T;}:̻.:>K>>N>.>>-=ݝ<=R +zǼvXӢ'UkaiĽՊB޽4o=:=J>E>.8>>>k==5uֻ:]l7:} S"+І^ PU$z9)jSh.-\`==[ȼp+|Nn[Wz̽qsf(=>1>u>מ>̃>u>W=88;ʒQ]I X{7H^yN弦I:O|˕k9ټju;m<\ !>=E#==4]==un> U<| =D>58>>R>'>>[j=m<򹲺!0?k !Ğ!rսԢؽMӜ{Bw/{ 0.^ȼI@!Ij#(ݽ!V ཐ,(M~ c" eS҆:\=u7=w>N>+Xw=s=Rd=YLM>? l? +7>c>/ =^;==/t`[8krh`-'мKD*o9mx)`<M*=F=(<< ڛ<!b<:A Gdo5;QIetR켆*<:ǹ׺*2:<>>;=26@>#>>>0>=ﲕ<)>1̑¼b1F)>ha⽈;sc&\ϼ_:C;1<<@<[>9w=e<(PfT|޼Q) &^A42R4 ]Վr?Lļ1&[KiJ'`(rFY3/FJ׼_zܽ\jT][5ҽhսh>o~E罹E1ZKڽri֢5W:~=. >1>_4>j>>>J=iݣFVDٻf1aO%ɼښPZj**ȽͭPD/?:@<4 =\=Bб<;/;|<"=U;i=>)V>Ma2>׾=Q==;z;h΋X6VνC.Bv;G@(8BlֽDJR{Jf_p^"5P]rԽ!ZA8 Sɵ".}zU!$Jw`=)r:8qOUZ]Nǽi"ݽ'*𼤷[M-$m;j==&=O<";m⇼޽ZC`CCoԽ .;F[ά12ޛQJ:N:%.7 (o&3#IQ/2.<1 +>;>AEf>5=Bm<>ERcC٢--:~޾V2(팽N/iǽ1;V=u =40U=5i=<9r<<]{]Q!ʼbt8xwW8Hפ"Zļa'C9*>|H>?>W}>R,M>Z->$w_>3Y>b>1=r6=ic=⽡9Ͻ񽢥Q2l'н3Z F<'A=us=Ԕ>g>> >=<\!$\~-C| 8NuI;"==?<;Qg>~Ƚ} atsk&C`;N<>=f=wo=[=D=U+=j=6b<I'i>Bs>z=ķ =|?> >),=>#=.=>!=o<<<=3=0=Ma2K6O;*<{===V=L<*;H_4Blp3DTpWMyżRԽ qRA(!f]8cQ᳽7WԽֽCaw ĽIQw)f( $O<`e>J=ǭ=H2<<ϻʽ;#iIzﾽ\(.½5`.'=[t=V=>c==e,i[76<4=,0{=oy>!Uj>:U>;L= =o +=zzr=fz9>*><ԕ>]=6=Hr=(=<=L<,I "v G#7࠻O;;H,^H8 :lٻ&>h>O>}@= =b=h==x.XԼt+q)~޽XkO;0heIa+DԽXps1Ż𻆴q<*=>>=q=S}+>gh]8;F<=n>> dA>N>B> a=u =;ώ;%2svn? +[幷*_?c\Dt==>1HF>>>g>R}>W:p>TQ>l>j>e>=U=0eꩡCK*»oǻT6Aqj0xeh sC +>Zu>>Z=G==_={ =B<<<lzk ټ⛽1u１u g>=>g >{UH>1G=n<;Z$Z7ҽMBf;i-:@#ߴ(WhͼWQO8K9d<$ռw[J꽰- #:ə<5L<=*-=>|>R>>e7>.==^<W;*fk+<D-e a={d>S?3?>b?\`?*>H>t =謌?)\A=)x=,!>h>O>>Y>oq>/=ЋN=!=x<;<ҟ;Q >.?= E>>==8z<ޱ<f;뻊41e߼)-\ͽ*"2V3~tʭFYHv<=>>/>#.=/===;==l#=e_x=/=/"T<䀳;SJR3`IMQB;lN<ħ=V>Fi>E2>V)>V>f&=n^=1;I1^{vȻ/im~r!Or3tn;ٸu) +>2?F?Z?.>b>X=C1<< +LikW,T;(;<{j'<$߻Q༟9μm6ؼżxȼgҼdkZ(+}v{"<<Ȳ=n +=ʞX>J>`>s~>>v6><01==t=V0=1<Ƶ:*[j:!l^~Ԣ;Z;Rb<>'d>o?<}?h[?"??:>=W;_|3 2>OZ>a>,>ik>3j> =:==5<><=7$2==T^Ud=ʢ)=~==*=#[<N>Wa>r4>s>: >>N>* >r=LU=k@=Bˏ44>==$'= ;#O5]L:o<`?$!?P?`?}?2>e=Zn=Rn>=և>x>Ю>B,>=GZ=.1=og3=9R#<<<@:e%<,W}>y?!?u??ږ?:,>(>D ;6lT;&^ļo/2>8>=V>====M@\;˭*@<0,=)=w={=E}A>]ҽQ䨻 ϻ<:@=/==Ͳn=°;𓽉2fP2h #^B=-%>9C>{>>~>>>K>2>A +>N4]>YӜ>7==AJ+<<>U?o3\?b?X?]?c>[>7&=Q<Ƕ<>޺ *ѻ;>S;P<E<< X+\w|71ht_whr&l|H;Fi:`LEJ'[;<"==Z=_>w>>5>8Sx>!> +== = <<2<1mߑ;f: ;C8; +c==>\MO>.?)?;K?1 +>?>j=:=阽H,xὧ>CyifmE==6Np =!|?di@PCr՟3<8=f>%$>,>=|CJ`󂽶Y}x6# =4>XP>0H>>T>X>f>Z>V&1>]Ŝ>P>![==7<|^;Ph;z;,hW=gԕ2,PUm.V<и=̽=&U=B?1?#?w;??w +> >!``<\<-)<<`jd;Ph<D<U +;~J;(1(jZݐ}r k+Źbӽڗy<%]<$<<;< <[=R=A7= ==|\=ڼ=R}=o==~Z=6*<<@<% +<w=2x=%F<]E< {=.&>|>[>8>G> f=_ ;[rN۽~$v}&4yI< 3yk>BO>'$=|x<&;(J)~ۅǽ]=xkP=_<ϹS">B>2U= +! glR6>;9=Ts> D>{>\>z>>9>zG>G >7AE>5=0=a=!z;s_ ȦT7eWCD|5+4Bc 4QJ-HQ)災}q-n\H>_??!IY? >>=y<E<;4<5\0❼' :(~߼ļ0W} 轖Iet v#59&I[[>#}6<ax>j>op>I=8"t;e>t=y=Jm=AV<9OH >/>6> =!~{oPYZNƃ =_*=Ӡ>hp>>)>>%>ӻ>CP>5= S==-=D eMk,Bi<齅]ʽi^Xh½3pg3"?&ȼn! =<ɂ=.==>I>>3>>:޵=s(ݘ+ɐ2t88n5JP^SzI۽L1ٻSh=6as==|6=f<佽 wZ㼫ZDp;-;p<<;'= @i=6L=} =s==d=={==Ϧ=d-=fa<6 ڽ4hcT5|fA<'hM49GER #]sir`v4=Vkt@C *E<J2̼|?<-G= =ne=ji>"Ϊ>Q*>W>;T=nn*i/AA +ʾv8M|=Wh=>R>;>) ) \B>=פ=Z). ؼ=ʼ?XN=S=>3=>oFO>5>:>z7>N><ј=R=5<;%<ũ<'9;:`hM^8mt9[% OY$k9`o$<=M=]=^L=N=v 6=~M=f=A?=5 ++=X[=\Cf<%:ʼn1 +^XĽ1 +kAteI;E>x=dz<DݪBG"۶M8ݾu hB=?=h2=2=<;}Zy#x_0GoA஽+RLcBA~:(<4=Z=̗/=y=ǜ=bՙŢN 弎2g>;>>>MW>Eŭ>k=<7AoKDk`chOF< 8=c=%< +;< <`;J! <<#6;5{䘄:;52B;XcT>!;AdҽgDF +F))~ʃ輟T8:iS幝Ch6=[=f=[]=p<\<@8=!{ӽ8ξ\<8=.|2=w=[==L.=l<;;YKѼ}wC+=(+ԅ;λ# 9;;<<$;qr:|;n<9ݷ?6񱭼`r.;q"<)zw<=j`==x>s> f=w=bEq#o*,:} +%=Me=:=8^>X=i=ef(>0> rc=g=;-㺕g#ܼbw:SJ˼ѽ!0S+% ܐW}s1;뒶;<(ݴ>>:7>A>y=Ԉ=J%`Ž=Y.j =j=&i<ʯ<έ==Q=hw<"%*P(x'.Ž+`{˼L, +`ՅL/U!;<}#ۑj:|23<pJ=ˬr= 佼S9p&＂#+C3nH8>?!3v?D>ŚI>sPb=1=< 8ͅZ߂vtJu#;#W8񃻒ܼ8*׻鏮tT:н$ ݽPEoh,Qh̽wT>XRb7y۽Oƽw{#W< +=:f=N=T@rm(*}ɽhǾV^Z5ȈFTj; ʂs}wܽqd6Hk꽤EJY+33<1-zֻּ;s>&2>(G!> @=AAeżn_&Ѽ#8;1k>f ?)M}?OO?F6?$s>1]>/K= }GS,ORݼ5l1軆嬼RgN5rMguH潐By}<Ӯ;"i' ?5)7c][_<-T:9Y<4=(T=<ߴ<'ZK鼍nd< +Uُ+|ؼ;HwHPF~pMהļV&d ޼ ּ5=og;֊<~c<$<ߛ<x>D>C+>g=l;61IllѼ N+lU< s=T=;:<;ml=u=ю=n==ax=q<:<\ < ;ݼ=hU1EH3%;F !1)~h'0D90?1c?ɼ<&H:XF$o>3? ?5?A i?$>>CWO=| +=$Vk9BF+7k5P|Lk;\;$YA9Zo41"Q]G70ýb[EϽ:tEj;;BӼ.o% B5Y)O;|C1]0ՔKμѶ0<<7Q==$=s=;Mϛg(u<;Rka_7IYz͒<-< ~mױfN<8=8=Y8=3h<<;ļBҽx눐9漳Zl<^ +3>==Cu: +HڗK(!DT*mϿ*׽IOüN2<=^4>'>>>>]? L>&8>]>=G'\<+ZlShxʼmݼ91h0; f;FV=L=(=R=G2=f=$i<`7h%;D<rۼyn9(˻I;<\7<~;߼;| 5i> }>=c=U;LI;r; +L OhR);v<Ր=CS=,<;ӲNOv2뼎J9K@z;<^+R=`>(ϭ>z>_ >f>|}G>G1=舼o؛dm0?5 1z =D=`?u +/iMc$<=sF>j>@>>_D:>O={=c=;h{ͪμaTC:DM/&;::N(;S<*c;Mƺ:!弹͓VrF+pǼZ4Tǘly[>=q4=ەM<=t=WzEgԽF"ͼ~y<JD%*Y# ;v:\I;<;3រRcּ8 o⻻=;T:~;ݒ=;5K=SJ=*6=f=MN> +>iE>hu>o>,>>MA>3==h=v=԰(s+eRXb;/K$O"ev:S;ovjr$Xӽ&J!m𼟔߼Xl(9K;;L; qBetba+(L^%c:OL3I_:9x;a|L[!=6ĎϽ"NAX5 5=t4&fG\y/u>:g>(t;==>ˤ>@v>e==5l=^=/ ;ٽҼ{JeBg BֶݼXҽ-4GPJm ^GWڼN/ =Di=e[=uY4===:>>B>he>N>_>eIo>#/=_=<;P&$q$p%igEtA7o(!eټVfXB +Gj{8m::O;є;!l;)v$3%/#z+*B8LPaQz) RYI：k½2伅?:]/<*z~t߽_ڽ `ƒ„ ;t=-> C>y&s>a'a>o>K=,Y>{>t>0p>*e>#M== [=4iTռX zϼnV0[纽F2ֽʽ GS_p6oؽeV߽*мf䟼wa;uj=Ko!=m==>nE> 1>0>JfS==1=q=V=k8ļDN;@{x'PI :Rbdz10RhRe\jāoZ(mŲü;e(l7$Ӽz\Z t_/F11ng5ysR"?G~2 o"hR;S;Β8:1?"սh?3k<#'P{AJCr;e:ho<<҂<<Ҕ<7]<ծ]c&ɽMltsB6Ľo;z)B=N>@>>h+'>#|=>H>pu>Hc>7=Q=G@;}C+\}/FآSa2 =Z3p,½miZ-޽u{&1U5ۓ,:<=8=#T=98=*I =a=fL=-H=?= }=a=T<;K<&<z". v̽5+T^@ٱՌý>2dj Gc_sqoK(V4ɽ߼!rp;'h/ {@~-(5A,2d{9V7bڽuMPSM==U?Sf +ݽ!۽:3;Ͻb-RBe'C%:Q; 4<~F%o>A>(Z=*<&Vh{BL-&꽆WNuS~i Q/ 2>: TWQ +N(\r.I 9μ+J՞.:-;AM<Y&?* 쇽w3➼"r;8:i YRM=@)@)]M<x&ͼvJƽ0(:f 6:Q(̱; +e|CT|[Pc%#Pv<}=.==kz=4bz<=>==뛐>'>"2=)M=9:ݽ~H IŽi;IyBǽ&ډŦFbfSP/rս&Ԝ˼7{/ +:=G< ;;ek;;JAh"fwJ)I=)+L琽^P|2Kt.:z-6~R]2쵽pD0b(6>rǼr>} zr 7Z*-EBI)jE4necO+15(%C*WVo)O8E7&t|=)K s[܌v=ļ̽7q蘼*"?nHX;K%<^0= ?~=&n(= #2=6=1=<<-,ʼv&=~=s=< *>~@=-k==u=`<0<8<3<{&󽠓OӒ轻2?pbf9W趼W-l`T!ཱི̢'WAj:6Ƽǻlr%^;:<$<8;Ǿ7̼R_8F|k)q5\ٽ?D%e"; zN_33JJ:{Nz!q IOs5mV +4H +Q^'aY:E +l;0 C>.>E7_>.> B=#=oV:J +;q3F>~sX=>==r=X2==o&=#=@`=.S Gquh*f80z7Hœ<-uX:t:ek ZqzBo}콆ͽD%Sλ5<5< +~-ϼ9sͽ<^ rYH<(OK>}>p==F<8AFڽEmfMAw=ZV#_e޻mѩ :FYL/¼޻3B ּܽPDgV&<}e> +1>?>j>;ai>e=ۋD=R:wz]q::襻^PLIcmQ0]˽WTż}򒼂Sn *f~bcּ)/K7>[6>]q>>=ŦE==l> =ڙ={l=8 >>1>>yX>?@=w=:<<;&!:`F+/$◼?꼑uax <<մ<;qBFOc&T<==>D5>o=;==C|du'ǯ;Jw;yL6#[L0t~)̽`ji9ýa_ㇽet) + NȽYq;˽u2|'aO Sfa~/p꼝RuT?eC?'?>Y>W>>=g<;=Mɼ4O<#o=8on}K<8l=@=v?=Զ=="}>*=ۈ>AO>>^)>0I>" =H=(>T>>4>W8==Se=_=o>=3=;ӻe*<>St<<*2AeUJnؽ jŽ>½h t/?'&?89?򽝯:KX?սyL1 ۻ":; B<=>F1?_?'Ms?1R>>==[J= +|O(ag%7A*7O'4ͽm.wu63C98~&;͘<%=O= &==m=;==h==X=h=y{==l<81&>D==+=>5Β==l=V=<<|;f3-=#=>I>>T>>.> $=bvWY?˧=)>>6>>KA>>< =e=?S= 9\伹漗M nK2&Uh6Ƚ|Xa޼k;BB^=*3==䝌=!9=y=v< = < p>|.=ޥ`=&=AN=pW={=j= = +< ;3>#<<:o=D5=0< =2e.¼ mr:5$)e\ƾsp@"ES;f<ٚ=g>p>,b>}>O><>=mр<^V >U+>] >,=Q=m_d='=<ܻmGtΕ \m߾9(3@;z==H===Ek=J< >1h>=#<\wM~;Յ<<:2BJbc3ROt;n;TZ;|36A;;P<2<=C=s>!R>XA>N9>)d== N=%%= +u<1=&=4$`uF.jMѼ>cͺnN^{;Z߼H<4;ʈ=+̽<<<1=v>,Vۼƻo;NMbVo].Z6GمTf^dzڻ= =H=>"D>m>>>">Z=a <&|mIط7ٻ"Q=v=]0=;M<=*e=D)j=R/;`Ո߼;- <=SQ=n0<umSC`ې>>=H>WMz>}t>hr>mq>7=tR<迊<.<=We=Q/7=<<5$7'ȼ𝼄 hSחxS]<m4x>F4>=Im=JP=?:YX0'R;c + ȼYkK s>D==n?p>c>>>;>]@L>q5= +ҏ<0;$< +<:͵! >l +>>j>N>M>$=tM;WbJ92ځt O$y3=;&_ +yU;u<<=?<ٿ=Q=-?==1=c<2<<,cL<]ǽ8 +jRApRv C2=>J>3>}E>V[>>{9> @=(;d^5xvᅟJ?yzU? +nO:<7:u~/= =ȼ<F|== :`lP kҊ$ Ƽ<։.hq; #<I'<]4\<9;o>'Tg+mۻr0u֩;ά:<><=h==?=#S<c;mD;;Ӥ8Hթ; +:z;<-x=1=&2O>@>|7>w>H$V>M`==ǐ6 -暽 ټy o&;I7L E@[ +1)~(FΫ;} +l=2=L=7>; 4=1>k> *=6=>=1S<=D=#==v=aB$OR-HAiq-=^Խؽaʽ6+4F~</= H<[<ƃl?ڵ"ͽ[;+3;C~<7M%<2N9#ʼ񋍽H=J&_%v<:!h;9yoL ;?SS"rȼ9m8#Cdg$ȼDjcܻ6FX=& =;,p<[=Wj=W=;=zl(%==A=2.:x9jzAV+1+jm;c6#5<5S;;ҼϼHһs:GK~;{D;P<^C<<8<ޕ=Ar=s5o=X2= +sn4>u>J%>*#;> ={=<3\< <5029\wyNOepż8;#;|;+`7q/y==z>=={=S=z==s<މQ h=R=Yئ2.46*aؽe?孽򂖽#5nO +Z}3K(x/,<-<<<::9{#r";!g +‚4;<JU` +:ӼD@ۅIA ߼ ػMӻμF'V_6>0>'S==@~;XTI(Fͽ:h=]= /=x8=</4a>83> +={=b;׼0˼/ +Ye"\;J<~<d=Zc=z=w=uA=/7= `m~s} ONڻ(=Wff=*= +r=q==A=r<<<fz;@ 4"$ DXw60]39=(zA˖w欏奻0q; k$*>Y +E>K7>=sA=s~1>>5W>1>y={0Y>Kg>)>ӳ>PH=Y=;k*<&3c<F<=?]t=a=G5=#ͼ5id ǽ'뫽h9OSQ4߼nڏỸ6ٺ ƺl99t:E;Ѳ<*<62>t1>zf>==  =$&==04<Քcduۼ':d4o"_<=VH=> ~>=§=<;(;Žd\˽ ⼐'&g;zg;tA*4; +<=l=6P=Z<;F7 N:{;Z:l%;k ?D??G?(3>T=E+<>3>>+T>h=MV=^=4=}-=r =.=.;V + *(l7p̼:⻥Ӽ-#¼0Oڨ| 4R +0< (8t6$מo3zq +ｵb=]j<:J0;׻<w=ļ;@3y<;mΜ+*ʽmA=ѻ}7K4Ȼׄ>'?)i?m5?R?2?g>> =J}=@2=,=f13=Aڃ=W=!=Hc= = <7C"k:Gxvbt+ܽe -i;m-㍼!;ݼĀs׽vAɻY}!3>K?(??G?s>؂>!C=}=\n=y=z=j=)$=k=m=aq=^=Qת~yCTsOֽu ȽBJ$(½[ܽ<&B*mRUE>Ҽ7󪼂9Ķ*< k .>'=H===NM?&0?x I?e w?3 >gw>==q>==wK=xX=c=bA=b=<{p6JɽZr#DC9Ͻ X[ýMg'V%yEu$ +Pټ/ռ;W<Q< $<{]f;滏<-b0.=JRpN%?-;ὢ8WؽTڼ*elrD[[+BqM, N漼Ω*켎RhxLo\E F;b>Z>}>2>K(==n;I;)b<@?=@p>1=l=/;M`ں:e<<2ڼPuO ýhy"ܡ<;QhϽMʼp>ڼhs=?<=)>>x>?>~>4><==#=>_'> y==kN=Z =s===:=G-[u;uK;6S彍au,ha:/!T=i)t;};W 2+;s<'<{^]J(T0>Z:LZ/:P}:vEA뻨?f`4V8vûӿY2K<=fp=UL=<;p';+==4>$f(>P_>pj>~>[>}>>Z =c<;׹<<;'=s=Y=t==}>/:>=[!U=%=ݷ=B|==SZ=h= 5#uĽ> HüʼdpAl o왽El@ļE%'.lO=:=i= =In >mi>E>Dg>>f>-=V=7=[=4=<=R=xS.;ļٹeHqiӽaN1꽭&.HIrv&8 :)3/h<>>=K=d}=W[;<_I< +l7F;/Hx`/beOμ:X<<}h<6= ,=<9v,;=J=>S>i>g>M],>;==כU=<[L<+>> =,=+6<R==Ϛ=K0N]8Vƽ[*@OXJ TJhѝ<':<];,>>o>cW=v=7<2m<=w0c=~>o> => +U="=Bi<(Xս0BY~&;k;&0[͜Xl]^gXsR.j~nݽ/L}T@"<s>>_z%>&= =bn,wQȽ޼ +=*=<˸ɽ3D'z)648j؟|(mZ׼Ycp;z=`$ ?(1hlt 8#,şqt-9૽ cר/dT潛p O H㼄a{R#Ӽ +MD";h=z> >Bw8>1>e>8 =]&;;Cw˼!='===7=Wl=,=o,N~:e `̽ ڽYk27e=@yqzΰd,JXM:LP녽7.(N4Vq>d>=i =wB=`O#yt>Fo>z>z>1$=\;%cHI=J=n=v=ʖ=}Q=J<l;;C)q\~+ нia0٫{{N`|jcO7>?>q==W>#P=J;[vI84:~( +U6g7`4~,r տ*߼tC. e)\Sj^Gڽ"2:Wa;R';!ٹ,A4aW 7( @Ǒ}k: حW<\=6<< +;#p |؀K +$z>"B=e<ɽoJ{<== ==%k\{< }K:ս^5qvYSS+~$ +hz4Ѽ^SC=@=-D@ZpR1J= !>&>]>,j>]=>hy>!0oT兽IqK:U/Yz"w^a𖽋󧽘uI'"쓽-;mq< <9%i .R* čT?=AN튽m<߼3^7%Sn' Ak";<^-<<x- ͽ ༟ܶh`u}jT%gռۼv*K +*B'<=˽6~,7T? j１]ɯ<ܻ +4=C==*"<ݼt:1s-_9=u4=\ҼD\6hb0==Ec>,4>>X>^M>?zIǽY_e#ņ}9;9Pd(j F9䝼s.N꾽a{l:zyѽj a&`uzUz{R^FݽIǽ<:}9样)<=B >p8><==U=O|Hu,<0U +ǼV;V< ;3Z;G69C~,L"kVEX!x+f2$Հ; +St +Qνѽכ/`-2VDV:<3>#>&2=3==W%:(*b{ҽȸCm;(⼐I22Cn'ws4 +Ⱥ;)#{;<;8KARX& +g~ebd!ż)s-@<.9!dV,)=;O==J<4%>!FJ=#=;Z ̛^ӽ6zӽ {i"FEo:r3Ta4F誽BٺW<𝕼 #R=楽W%a|M6<򼼼AOc# #v3fDe/62kP=t>,K>P>>h>J.=TU=8k;gJ )cѴv-9sZ< K<+f;;t\[l͍Bzׄl|'B> w+ +3l}F 4 ý#)uĽ9hgBrdof$'9<^H=EN=I=9=h8=Kb;_;!:'.:t; 3>?h>_e{>g==5K=ú>DvqEoƾG}̼Žy1";H;/t:[/$A; < <[FmF!= +"ڙ=nf ;]`F{?g$t"7Ts.D~w=D&t>V>N>w>>>=S+溝7MBcM#@("?|k;8<<<|>>/?>|>wۜ>=o;ѼR'N!ǽ҂ ڇ?ɲj2_LSJ0_. e T;U;h:;UN~9Mm KEӽmսCݼVp듽iQaɼ?qb< m'>lQ>/><0>k1=v=5h^>=9==,X<:*;  <X<|[>>n>־>q>>86<=7~\$=m='A=y=65A=*2=<a<<=W=w<%<&=%=l=v2=qe]=#IGK +p_}ԼkϽ!r=V[}z1kؕC +Zu򎽎 x"ˌ<3i]NHoK ^r%k@qUpE;l<mt(:+;Wz QAUfѼ猼ծT<|V==R>-H>nG>s3>S=a=73Pa,;ypV_C? <\;E0~oƽDoؽ HZڽ)c~;w.c>nn>fpQ>@:>q=ύR=0>>->>L>$> =>%Z=Kn==v<<G=Rp<nj:#x2 u\kp%M6-Zͽ4Zi4=G= /=o#==SZ=H`=='s=M=t=x=};F* v8>==On>f>)7>b>>9=g=;/SNI J;P<S;;_7CJ=qE3_ JLm_${l+M쥻N!NKŽQ_T6'ӏ<ܘk=>a4>>6>>H=1= <jɽM28]l,@&bTK=Z8>/`>Zs>mc>a~>Jg>=\ == D<۪'>> >>>T>>=9r+4q ߼'O( ;r["2߼/Wݽ-|o5}Ӵm`1e-z'7P )<8[=ŧ{>$9>/Bt>Ks+=]=7=Vj}>G>8*>>x@>*v==<Ҽ1c'YJ; eY{=|;EL>i;L<5L=>=~ =~==A[=J\;=cٻ!P<#=r= =R=@X=*==38='u4AbǽԼH7<͈L>I\>p>->Ն>S>%`4=5ȜER @q+}Ž1\Vwܶ+e@!:GD=*3=^=wc->q>> >>U\> F=2;=1;p ]\ủm>1=* ]ŻN;E=Я===2>k7>z=ڶamJc!cbkqSXˈSq#sxEC>I >{C>>c>@o=;τf`>}>xJ>]C>m>g;=C=8<kPgDov0[=<@=E@i= Z==[<ֲ<<:al=n>+>"b=]==K:k;q=.=P=v>==Ǭ=Fo"\?:(a9';==n$={j=Ȫ===D= v==[-=t=°o=f=L3=xOK5c*s>FR>>s=<= M=q<-S磼iW缳ۭ:87"J< +=c=^=<;Ѽ*"<֞=Z=F>> \ܥ{=9f=8H~<4<;ҽz‾_{.* #;?==#=v====¶=3=[f<,:, [8%<}Ѽ;<7B= m<먎&>B>-!>n=u==)z=44>JJa=<;h"R⎽yELM8haryV߃ϗW=5缸\w<<<<=P +>={=?b=ѧ>>MO>N><8˼V#0`.'<^<._">Taq>@> v>V==d=rG> >>`> !?MT>=3>$"==բNO0k{hNν\өbqF>>|>-=> >a>p>MJ>IpQ=kݻҸ;uۈO!">-> +=P=ӔS=F=h;#L7-&,*񦼣 ̼@0*wGa<<=s=Fs==x8==|=j=Z:j>>&>$>6a=ӻkT_[\\o¼ݽ5罽#Q̽O\<_=jz>>$>E>K>J>=< a[L~&&c_@<۶=&=_"=3=4;=߻}e/~ʽ?:du7½_ܽY%ݦK,~';i ,> e>>E>TT=={;Ǧ<\P=i< X=Y=:y ;=]M=;>_=O>zB=L=v(=A;3Y'b9$4W糽j?$u%$+ռ۽h)ބ'! +;;9q >6>mi>|>b =CА?FL^|nrk~M7x.C=9Y>%97>?4?|?oɸ?,>>=e=}nA>I-B=== F}ĬI%^Ѽ-59lLnFAKPM/>4+pN$:ھ$4=D^==2=5< \ɔ;n<\f=pK====Ê=S=<:8< d;)>/w> l=2Y#rнѼTٻq3= t=R>=-=tJ= =S=V; =s+;Qq`Lv%_mw6(06J(;1K= +Հ=gy==a=Ա=<<;'H7߻ua;L 7<=K>>$\> =ڏ<ҽZ\D0Q0FlѭD+;t=>9X?OC??"??b̫?<>O==GLɼ; +=0>;">>">v>6u'=6ɐ* zu5K꼽V MCz=w}p}S4!U,M?ɾ*XVyP9t<ż#:T =t=G=zD;!;6S=_=S= +=}=-o!>">L>n=*Z ۽)|}/ : 2=UZ=/>>ؾ=s=KL){GU+:}E<@ =-=;==\3<;;lȼع7< A??h?y1?c?ˢ)?)?=> D==ƈ6;u =UO>0e>X>>>kП= <*<- Ih"?:+E\^lbO݁ǽ ƻ5s">9r5> $>=8=2 +dPӼ"5&ֽ+aɽB+߼$[w?gSϼSa;d"|2u{B6WnQq B꽃ڽH}ν)뿜=H>{??;s= /)>#>V>u A>S=Xp>#>=w==<9T`[ѽo8b+[B-ǽ<0Xa*;;!~=%|v>_?? ?*?Ӊ? >(:>=V;bbݼ_d,;b<4=8>>)]>!=ՅM=p<;?"ﻦ^;Q;D=|S`fBQs,u/ɾ*{YMgQ<= tB;ۈIR>?!+?#X?>-_=K<\ϼνY!&fnI;j4 <|=b=LV=uN=^=ϵ<>>az>'K#=s<y/;le#-5iٽYjJ /3F͊B08w$/W!Y뮼GU0@pPCFw?kEXsL~!찼޶?x׼,ꜽ:ȖCU<==3>,>=ʊ;H@Ν*|bfhɼH=<{Q< +Ľh奮>=!j==N=<;ν0 kֽ%;DS;:< ) OÚ꽑B|b8ɽ(ؽ;|ӼKһ9Hj罳DR޽^d}}ɼ%M04DaGy'<~<*7}K 1;4:Y%ʽ nR?Լ&<*s>x^cVҽx~ýRǼDِfvܽ .}G[Ipp;ʽD)W54\@7C <{>6>QRQ>ep> +=,=T<޺ +{3Y1==X=;"<]b IнQuo=58=={> ==m5W9's7 ywR@נ; 92ܻ%;r@pt>&6>(0=U=8N=&;A$x~?<3=q$"=a=F"7~ƯTcr$]N'' Z.2ȼ2Q){9@y:4W<;?B7=Dv=g=]>()>@M>-S>=i==Y?i=x*3Z>+W>%>I>>Bk>>k>>N==s +s;c:=%=˕=C=Ԙ=;1ͽ'YxxMXsҽ|ͽфƮ=D(_uF;^6n# <*0ݽ@z*FD߽)߽SrIyq1A@׼a':5;_]<+<9lPL], kOؼ ߼=ӼW~kQɺ5._;3+<$B;x; ;y((~){ +6rSWІ!87B`+$<)>.CH>;b>Wf?I?yc0?^J?A$? i>޳>[~>>m=;utmMjsyV C1G @֦|ڢ=c^PW +4bDͶ]K8 ݱt'/Z* ;-<ںTt[rmS*h>:?î@-@EV@@[@6+?wv?+k>im>>-=D;+V ¼'= x==.=𼯂2ƻ(wՂm܉n?fb3 +=hs>?@5@@@@"]@??>c>^3=<Q$Nؽ@M<'<@>ͼ8Yz8%fQu txOb4ؼ毽<\gQRsWO?8@jx@A3aA6cAf@@G?fK>#>2=:[fss mJHmϽ-8܎`L1jX`/f#Ǽ#g@;-F+o̽W_PEdFEeKқ-M5>KŁq_o }QW<=Bp.=oX=U'<~p;(JȽoN|;EҽKK⪒=o7=ëe6!d0+5 Uk<:0=3;č;\;l:>&yHaWR\_ZziA󽂣4*ͽMHxoܼ vw4 ;.o#?@b8@RA<AVABjA^@?2>>=r@UŽf' ]k>,ƽ&ܽGh.Ѵv$4|BcFмjoDs 4[Sz̽ ~Ƽ#ۢ:w%ս6<\4;%<ᮻ~(xPҽH콤d;> ~=(;./ ]ҽ56ĽÑvld Iy_:[6Ͻ,w8촚@E;=?a=(=$=4սVս2" |y^ZD +i=Xm=(6f@¢n2$<:&uf==\= t&,%ݡ<^|=ļU:o$9O<= =;.Cm`ݽ#[2ս6 OݎP:Y;^AUs?@UM@ї@'@g@~s??<+>Kb[=<=;ռ%ZڽK|YDwWgSs4]K]ֽ&<%<`k"jq;<=2= +R=6=#]<` 9h'漗u9dHۼz* I;;4p_>G>8>N$𥼓w A`a0=u:üu񗝽Eн4=Jr +43*9$W=ڼ=XY<#ED0e; +X== T>Q>?ܜ?ɾ@7?9?>}>[== =O< +˽S;~ٽV\[wG½ĵȽd걼ë"{= =<)<@U>t>xb>':%MOF,$־!@<ߦ?a]ZֻEQ=sj=&e=_;b0 o==eQ=m==<) >]>q +>x>|>v>><>= =T<>j|P>N=>-V.A6ݖKV򽋄=7>#=Z}!;٬<:D<=U,=V7<=,c= +O= =m`?=:=-y=2=}=eM=w3%=>n4>R==$<:*q}oX;O}煽lͽU0}>AVZy5rLL< < y<<-],e'5;>^>P:>=J= 7=l=ezp6Q'E= +=Us>^:_e=._<=ѷ>)A.>>Z>9#>ƵB=EG3/*\DZ0%>= N="=W<;<\=A=Ϡ=>!%>l>f>=X0غ<=Ct{==Z<t<,K[;q%=. =:=*l=e_=aZt='i=C <<%;!;|&<<6^ټ &H'Gwdud5 RN_\mh袘#,@ ah3(>0;$ib&4潗Ž3fU9FbŽvj+ϽL]Ǿڽq ·w:̺A&> W'>'>= Q>Lޙ=SP"4TW B~DtƼ#3ż]EO|f2¾NI=lj iUrޮuŽ'뼑Yʻ٥;<= ==_{=\= ]<<{v`9L;An;A> g>8>&>۪>! |Dz.Š=]ڻ[Ͻ2g- +:j;;=_=2h=1 "/^SuW~Pyw=6=x :>>#K>J>u>#J>fg sؽ0:=L =F==-Y2+6ʽ`g/;u<<}a_>5>!#>bUr>%>U>֧8='="=ɭ{=8=4=:;9J==p@=!=<Қ4,xݍPM^-:z<,UHJ\ ?3< wgf<];x ̻…{ +ҽ!~` +ֽv!x}n<> +ߺ>B=ؽ$!eN(<);X>2>l>:>Q>>H>C/>L1X==j= =ۃU=E>u>7JX> .=?<.=L=o<;.`A`^*U#oC;O_<>C >>ҸA>=߼ٺ!?½½W<~@==hK= =Z=Pi=ż=yW=4(<(;06#<=;@$<_=O=ˎ==U>J>:>D>B]>>;HE>9=悴>26>Q>2==Wm;X;<2; 3: &!up߽.Y=?F=#MtW2׿:|0ー;f:<fʼ숻<7=WP7= =V=-2<=x<ڃ< }Uļ*<W=b +=/=d=9r|=oV==-=kv>3>w?&?Yf?Oa?.>H= (:x4>;@>>C \ No newline at end of file diff --git a/lib/test/data/2chipA.fits.REMOVED.git-id b/lib/test/data/2chipA.fits.REMOVED.git-id new file mode 100644 index 0000000..50242ca --- /dev/null +++ b/lib/test/data/2chipA.fits.REMOVED.git-id @@ -0,0 +1 @@ +fe2343477504847a56ac369c96fbe2d1b9d835db \ No newline at end of file diff --git a/lib/test/data/2chipB.fits.REMOVED.git-id b/lib/test/data/2chipB.fits.REMOVED.git-id new file mode 100644 index 0000000..1d70532 --- /dev/null +++ b/lib/test/data/2chipB.fits.REMOVED.git-id @@ -0,0 +1 @@ +8a32e201693484546f631c52e3774dfbf278151b \ No newline at end of file diff --git a/lib/test/data/2chipC.fits.REMOVED.git-id b/lib/test/data/2chipC.fits.REMOVED.git-id new file mode 100644 index 0000000..7afd241 --- /dev/null +++ b/lib/test/data/2chipC.fits.REMOVED.git-id @@ -0,0 +1 @@ +167e2502832be1609cdde29c52e9dc8c571cc3ae \ No newline at end of file diff --git a/lib/test/test.py b/lib/test/test.py new file mode 100644 index 0000000..c6c1fcd --- /dev/null +++ b/lib/test/test.py @@ -0,0 +1,260 @@ +import os +import random + +import numpy as np +from numpy.testing import assert_almost_equal, assert_array_less + +from sphere import graph +from sphere import great_circle_arc +from sphere import polygon +from sphere import vector + +from .test_util import * + + +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(os.path.join(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) diff --git a/lib/test/test_intersection.py b/lib/test/test_intersection.py new file mode 100644 index 0000000..81512d1 --- /dev/null +++ b/lib/test/test_intersection.py @@ -0,0 +1,167 @@ +from __future__ import print_function + +# 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 sphere import polygon + +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): + polys = func(*args, **kwargs) + + intersections = [] + 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) + 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) + areas = [x.area() for x in permutation] + intersection_area = intersection.area() + assert np.all(intersection_area < areas) + + 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() + + 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]) + + return run + + +@intersection_test(0, 90) +def test1(): + import pyfits + fits = pyfits.open(os.path.join(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(os.path.join(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) + + + +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_union.py b/lib/test/test_union.py new file mode 100644 index 0000000..66abda8 --- /dev/null +++ b/lib/test/test_union.py @@ -0,0 +1,169 @@ +from __future__ import print_function + +# 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 sphere import polygon + +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): + 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) + for method in ('parallel', 'serial'): + for i, permutation in enumerate( + itertools.islice( + itertools.permutations(polys), + None, None, step_size)): + filename = '%s_%s_union_%04d.svg' % ( + func.__name__, method, i) + print(filename) + + union = polygon.SphericalPolygon.multi_union( + permutation, method=method) + unions.append(union) + areas = [x.area() for x in permutation] + union_area = union.area() + assert np.all(union_area >= areas) + + 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() + + 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(os.path.join(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) + + +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 new file mode 100644 index 0000000..47a3444 --- /dev/null +++ b/lib/test/test_util.py @@ -0,0 +1,14 @@ +import os + +import numpy as np +from sphere 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 new file mode 100644 index 0000000..7323c92 --- /dev/null +++ b/lib/vector.py @@ -0,0 +1,213 @@ +# -*- 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. +""" + +# 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): + ur""" + 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): + ur""" + 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): + ur""" + 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) + + u2 = u**2 + v2 = v**2 + w2 = w**2 + 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 -- cgit