diff options
Diffstat (limited to 'stsci/sphere/polygon.py')
| -rw-r--r-- | stsci/sphere/polygon.py | 831 |
1 files changed, 831 insertions, 0 deletions
diff --git a/stsci/sphere/polygon.py b/stsci/sphere/polygon.py new file mode 100644 index 0000000..e9ad379 --- /dev/null +++ b/stsci/sphere/polygon.py @@ -0,0 +1,831 @@ +# -*- coding: utf-8 -*- + +# Copyright (C) 2011 Association of Universities for Research in +# Astronomy (AURA) +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions +# are met: +# +# 1. Redistributions of source code must retain the above +# copyright notice, this list of conditions and the following +# disclaimer. +# +# 2. Redistributions in binary form must reproduce the above +# copyright notice, this list of conditions and the following +# disclaimer in the documentation and/or other materials +# provided with the distribution. +# +# 3. The name of AURA and its representatives may not be used to +# endorse or promote products derived from this software without +# specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY AURA ``AS IS'' AND ANY EXPRESS OR +# IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL AURA BE LIABLE FOR ANY DIRECT, +# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, +# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED +# OF THE POSSIBILITY OF SUCH DAMAGE. + +""" +The `polygon` module defines the `SphericalPolygon` class for managing +polygons on the unit sphere. +""" +from __future__ import division, print_function, unicode_literals, absolute_import + +# STDLIB +from copy import copy, deepcopy + +# THIRD-PARTY +import numpy as np + +# LOCAL +from . import great_circle_arc +from . import vector + +__all__ = ['SphericalPolygon'] + + +class SphericalPolygon(object): + r""" + Polygons are represented by both a set of points (in Cartesian + (*x*, *y*, *z*) normalized on the unit sphere), and an inside + point. The inside point is necessary, because both the inside and + outside of the polygon are finite areas on the great sphere, and + therefore we need a way of specifying which is which. + """ + + def __init__(self, points, inside=None): + r""" + Parameters + ---------- + points : An Nx3 array of (*x*, *y*, *z*) triples in vector space + These points define the boundary of the polygon. It must + be "closed", i.e., the last point is the same as the first. + + It may contain zero points, in which it defines the null + polygon. It may not contain one, two or three points. + Four points are needed to define a triangle, since the + polygon must be closed. + + inside : An (*x*, *y*, *z*) triple, optional + This point must be inside the polygon. If not provided, the + mean of the points will be used. + """ + if len(points) == 0: + # handle special case of initializing with an empty list of + # vertices (ticket #1079). + self._inside = np.zeros(3) + self._points = np.asanyarray(points) + return + elif len(points) < 3: + raise ValueError("Polygon made of too few points") + else: + assert np.array_equal(points[0], points[-1]), 'Polygon is not closed' + + self._points = np.asanyarray(points) + + if inside is None: + self._inside = np.mean(points[:-1], axis=0) + else: + self._inside = np.asanyarray(inside) + + # TODO: Detect self-intersection and fix + + def __copy__(self): + return deepcopy(self) + + def __repr__(self): + return '%s(%r, %r)' % (self.__class__.__name__, + self.points, self.inside) + + @property + def points(self): + """ + The points defining the polygon. It is an Nx3 array of + (*x*, *y*, *z*) vectors. The polygon will be explicitly + closed, i.e., the first and last points are the same. + """ + return self._points + + @property + def inside(self): + """ + Get the inside point of the polygon. + """ + return self._inside + + def to_radec(self): + """ + Convert `SphericalPolygon` footprint to RA and DEC. + + Returns + ------- + ra, dec : list of float + List of *ra* and *dec* in degrees corresponding + to `points`. + """ + return vector.vector_to_radec(self.points[:,0], self.points[:,1], + self.points[:,2], degrees=True) + + @classmethod + def from_radec(cls, ra, dec, center=None, degrees=True): + r""" + Create a new `SphericalPolygon` from a list of (*ra*, *dec*) + points. + + Parameters + ---------- + ra, dec : 1-D arrays of the same length + The vertices of the polygon in right-ascension and + declination. It must be \"closed\", i.e., that is, the + last point is the same as the first. + + center : (*ra*, *dec*) pair, optional + A point inside of the polygon to define its inside. If no + *center* point is provided, the mean of the polygon's + points in vector space will be used. That approach may + not work for concave polygons. + + degrees : bool, optional + If `True`, (default) *ra* and *dec* are in decimal degrees, + otherwise in radians. + + Returns + ------- + polygon : `SphericalPolygon` object + """ + # Convert to Cartesian + x, y, z = vector.radec_to_vector(ra, dec, degrees=degrees) + + if center is None: + xc = x.mean() + yc = y.mean() + zc = z.mean() + center = vector.normalize_vector(xc, yc, zc) + else: + center = vector.radec_to_vector(*center, degrees=degrees) + + return cls(np.dstack((x, y, z))[0], center) + + @classmethod + def from_cone(cls, ra, dec, radius, degrees=True, steps=16.0): + r""" + Create a new `SphericalPolygon` from a cone (otherwise known + as a "small circle") defined using (*ra*, *dec*, *radius*). + + The cone is not represented as an ideal circle on the sphere, + but as a series of great circle arcs. The resolution of this + conversion can be controlled using the *steps* parameter. + + Parameters + ---------- + ra, dec : float scalars + This defines the center of the cone + + radius : float scalar + The radius of the cone + + degrees : bool, optional + If `True`, (default) *ra*, *dec* and *radius* are in + decimal degrees, otherwise in radians. + + steps : int, optional + The number of steps to use when converting the small + circle to a polygon. + + Returns + ------- + polygon : `SphericalPolygon` object + """ + u, v, w = vector.radec_to_vector(ra, dec, degrees=degrees) + if degrees: + radius = np.deg2rad(radius) + + # Get an arbitrary perpendicular vector. This be be obtained + # by crossing (u, v, w) with any unit vector that is not itself. + which_min = np.argmin([u, v, w]) + if which_min == 0: + perp = np.cross([u, v, w], [1., 0., 0.]) + elif which_min == 1: + perp = np.cross([u, v, w], [0., 1., 0.]) + else: + perp = np.cross([u, v, w], [0., 0., 1.]) + + # Rotate by radius around the perpendicular vector to get the + # "pen" + x, y, z = vector.rotate_around( + u, v, w, perp[0], perp[1], perp[2], radius, degrees=False) + + # Then rotate the pen around the center point all 360 degrees + C = np.linspace(0, np.pi * 2.0, steps) + # Ensure that the first and last elements are exactly the + # same. 2π should equal 0, but with rounding error that isn't + # always the case. + C[-1] = 0 + C = C[::-1] + X, Y, Z = vector.rotate_around(x, y, z, u, v, w, C, degrees=False) + + return cls(np.dstack((X, Y, Z))[0], (u, v, w)) + + @classmethod + def from_wcs(cls, fitspath, steps=1, crval=None): + r""" + Create a new `SphericalPolygon` from the footprint of a FITS + WCS specification. + + This method requires having `pywcs` and `pyfits` installed. + + Parameters + ---------- + fitspath : path to a FITS file, `pyfits.Header`, or `pywcs.WCS` + Refers to a FITS header containing a WCS specification. + + steps : int, optional + The number of steps along each edge to convert into + polygon edges. + + Returns + ------- + polygon : `SphericalPolygon` object + """ + import pywcs + import pyfits + + if isinstance(fitspath, pyfits.Header): + header = fitspath + wcs = pywcs.WCS(header) + elif isinstance(fitspath, pywcs.WCS): + wcs = fitspath + else: + wcs = pywcs.WCS(fitspath) + if crval is not None: + wcs.wcs.crval = crval + xa, ya = [wcs.naxis1, wcs.naxis2] + + length = steps * 4 + 1 + X = np.empty(length) + Y = np.empty(length) + + # Now define each of the 4 edges of the quadrilateral + X[0 :steps ] = np.linspace(1, xa, steps, False) + Y[0 :steps ] = 1 + X[steps :steps*2] = xa + Y[steps :steps*2] = np.linspace(1, ya, steps, False) + X[steps*2:steps*3] = np.linspace(xa, 1, steps, False) + Y[steps*2:steps*3] = ya + X[steps*3:steps*4] = 1 + Y[steps*3:steps*4] = np.linspace(ya, 1, steps, False) + X[-1] = 1 + Y[-1] = 1 + + # Use wcslib to convert to (ra, dec) + ra, dec = wcs.all_pix2sky(X, Y, 1) + + # Convert to Cartesian + x, y, z = vector.radec_to_vector(ra, dec) + + # Calculate an inside point + ra, dec = wcs.all_pix2sky(xa / 2.0, ya / 2.0, 1) + xc, yc, zc = vector.radec_to_vector(ra, dec) + + return cls(np.dstack((x, y, z))[0], (xc[0], yc[0], zc[0])) + + def _unique_points(self): + """ + Return a copy of `points` with duplicates removed. + Order is preserved. + + .. note:: Output cannot be used to build a new + polygon. + """ + val = [] + for p in self.points: + v = tuple(p) + if v not in val: + val.append(v) + return np.array(val) + + def _sorted_points(self, preserve_order=True, unique=False): + """ + Return a copy of `points` sorted such that smallest + (*x*, *y*, *z*) is on top. + + .. note:: Output cannot be used to build a new + polygon. + + Parameters + ---------- + preserve_order : bool + Preserve original order? If `True`, polygon is + rotated around min point. If `False`, all points + are sorted in ascending order. + + unique : bool + Exclude duplicates. + """ + if len(self.points) == 0: + return [] + + if unique: + pts = self._unique_points() + else: + pts = self.points + + idx = np.lexsort((pts[:,0], pts[:,1], pts[:,2])) + + if preserve_order: + i_min = idx[0] + val = np.vstack([pts[i_min:], pts[:i_min]]) + else: + val = pts[idx] + + return val + + def same_points_as(self, other, do_sort=True, thres=0.01): + """ + Determines if this `SphericalPolygon` points are the same + as the other. Number of points and areas are also compared. + + When `do_sort` is `True`, even when *self* and *other* + have same points, they might not be equivalent because + the order of the points defines the polygon. + + Parameters + ---------- + other : `SphericalPolygon` + + do_sort : bool + Compare sorted unique points. + + thres : float + Fraction of area to use in equality decision. + + Returns + ------- + is_eq : bool + `True` or `False`. + """ + self_n = len(self.points) + + if self_n != len(other.points): + return False + + if self_n == 0: + return True + + self_a = self.area() + is_same_limit = thres * self_a + + if np.abs(self_a - other.area()) > is_same_limit: + return False + + if do_sort: + self_pts = self._sorted_points(preserve_order=False, unique=True) + other_pts = other._sorted_points(preserve_order=False, unique=True) + else: + self_pts = self.points + other_pts = other.points + + is_eq = True + + for self_p, other_p in zip(self_pts, other_pts): + x_sum = 0.0 + + for a,b in zip(self_p, other_p): + x_sum += (a - b) ** 2 + + if np.sqrt(x_sum) > is_same_limit: + is_eq = False + break + + return is_eq + + def contains_point(self, point): + r""" + Determines if this `SphericalPolygon` contains a given point. + + Parameters + ---------- + point : an (*x*, *y*, *z*) triple + The point to test. + + Returns + ------- + contains : bool + Returns `True` if the polygon contains the given *point*. + """ + P = self._points + r = self._inside + point = np.asanyarray(point) + + intersects = great_circle_arc.intersects(P[:-1], P[1:], r, point) + crossings = np.sum(intersects) + + return (crossings % 2) == 0 + + def intersects_poly(self, other): + r""" + Determines if this `SphericalPolygon` intersects another + `SphericalPolygon`. + + This method is much faster than actually computing the + intersection region between two polygons. + + Parameters + ---------- + other : `SphericalPolygon` + + Returns + ------- + intersects : bool + Returns `True` if this polygon intersects the *other* + polygon. + + Notes + ----- + + The algorithm proceeds as follows: + + 1. Determine if any single point of one polygon is contained + within the other. + + 2. Deal with the case where only the edges overlap as in:: + + : o---------o + : o----+---------+----o + : | | | | + : o----+---------+----o + : o---------o + + In this case, an edge from one polygon must cross an + edge from the other polygon. + """ + assert isinstance(other, SphericalPolygon) + + # The easy case is in which a point of one polygon is + # contained in the other polygon. + for point in other._points: + if self.contains_point(point): + return True + for point in self._points: + if other.contains_point(point): + return True + + # The hard case is when only the edges overlap, as in: + # + # o---------o + # o----+---------+----o + # | | | | + # o----+---------+----o + # o---------o + # + for i in range(len(self._points) - 1): + A = self._points[i] + B = self._points[i+1] + if np.any(great_circle_arc.intersects( + A, B, other._points[:-1], other._points[1:])): + return True + return False + + def intersects_arc(self, a, b): + """ + Determines if this `SphericalPolygon` intersects or contains + the given arc. + """ + return self.contains_point(great_circle_arc.midpoint(a, b)) + + if self.contains_point(a) and self.contains_point(b): + return True + + P = self._points + intersects = great_circle_arc.intersects(P[:-1], P[1:], a, b) + if np.any(intersects): + return True + return False + + def area(self): + r""" + Returns the area of the polygon on the unit sphere. + + The algorithm is not able to compute the area of polygons + that are larger than half of the sphere. Therefore, the + area will always be less than 2π. + + The area is computed by transforming the polygon to two + dimensions using the `Lambert azimuthal equal-area projection + <http://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection>`_ + + .. math:: + + X = \sqrt{\frac{2}{1-z}}x + + .. math:: + + Y = \sqrt{\frac{2}{1-z}}y + + The individual great arc circle segments are interpolated + before doing the transformation so that the curves are not + straightened in the process. + + It then uses a standard 2D algorithm to compute the area. + + .. math:: + + A = \left| \sum^n_{i=0} X_i Y_{i+1} - X_{i+1}Y_i \right| + """ + if len(self._points) < 3: + #return np.float64(0.0) + return np.array(0.0) + + points = self._points.copy() + + # Rotate polygon so that center of polygon is at north pole + centroid = np.mean(points[:-1], axis=0) + centroid = vector.normalize_vector(*centroid) + points = self._points - (centroid + np.array([0, 0, 1])) + vector.normalize_vector( + points[:, 0], points[:, 1], points[:, 2], inplace=True) + + X = [] + Y = [] + for A, B in zip(points[:-1], points[1:]): + length = great_circle_arc.length(A, B, degrees=True) + interp = great_circle_arc.interpolate(A, B, length * 4) + x, y, z = vector.normalize_vector( + interp[:, 0], interp[:, 1], interp[:, 2], inplace=True) + x, y = vector.equal_area_proj(x, y, z) + X.extend(x) + Y.extend(y) + + X = np.array(X) + Y = np.array(Y) + + return np.abs(np.sum(X[:-1] * Y[1:] - X[1:] * Y[:-1]) * 0.5 * np.pi) + + def union(self, other): + """ + Return a new `SphericalPolygon` that is the union of *self* + and *other*. + + If the polygons are disjoint, they result will be connected + using cut lines. For example:: + + : o---------o + : | | + : o---------o=====o----------o + : | | + : o----------o + + Parameters + ---------- + other : `SphericalPolygon` + + Returns + ------- + polygon : `SphericalPolygon` object + + See also + -------- + multi_union + + Notes + ----- + For implementation details, see the :mod:`~sphere.graph` + module. + """ + from . import graph + g = graph.Graph([self, other]) + + polygon = g.union() + + return self.__class__(polygon, self._inside) + + @classmethod + def multi_union(cls, polygons): + """ + Return a new `SphericalPolygon` that is the union of all of the + polygons in *polygons*. + + Parameters + ---------- + polygons : sequence of `SphericalPolygon` + + Returns + ------- + polygon : `SphericalPolygon` object + + See also + -------- + union + """ + assert len(polygons) + for polygon in polygons: + assert isinstance(polygon, SphericalPolygon) + + from . import graph + + g = graph.Graph(polygons) + polygon = g.union() + return cls(polygon, polygons[0]._inside) + + @staticmethod + def _find_new_inside(points, polygons): + """ + Finds an acceptable inside point inside of *points* that is + also inside of *polygons*. Used by the intersection + algorithm, and is really only useful in that context because + it requires existing polygons with known inside points. + """ + if len(points) < 4: + return np.array([0, 0, 0]) + + # Special case for a triangle + if len(points) == 4: + return np.sum(points[:3]) / 3.0 + + for i in range(len(points) - 1): + A = points[i] + # Skip the adjacent point, since it is by definition on + # the edge of the polygon, not potentially running through + # the middle. + for j in range(i + 2, len(points) - 1): + B = points[j] + C = great_circle_arc.midpoint(A, B) + in_all = True + for polygon in polygons: + if not polygon.contains_point(C): + in_all = False + break + if in_all: + return C + + raise RuntimeError("Suitable inside point could not be found") + + def intersection(self, other): + """ + Return a new `SphericalPolygon` that is the intersection of + *self* and *other*. + + If the intersection is empty, a `SphericalPolygon` with zero + points will be returned. + + If the result is disjoint, the pieces will be connected using + cut lines. For example:: + + : o---------o + : | | + : o---------o=====o----------o + : | | + : o----------o + + Parameters + ---------- + other : `SphericalPolygon` + + Returns + ------- + polygon : `SphericalPolygon` object + + Notes + ----- + For implementation details, see the :mod:`~sphere.graph` + module. + """ + # if not self.intersects_poly(other): + # return self.__class__([], [0, 0, 0]) + + from . import graph + if len(self._points) < 3 or len(other._points) < 3: + return self.__class__([], [0, 0, 0]) + + g = graph.Graph([self, other]) + + polygon = g.intersection() + + inside = self._find_new_inside(polygon, [self, other]) + + return self.__class__(polygon, inside) + + @classmethod + def multi_intersection(cls, polygons, method='parallel'): + """ + Return a new `SphericalPolygon` that is the intersection of + all of the polygons in *polygons*. + + Parameters + ---------- + polygons : sequence of `SphericalPolygon` + + method : 'parallel' or 'serial', optional + Specifies the method that is used to perform the + intersections: + + - 'parallel' (default): A graph is built using all of + the polygons, and the intersection operation is computed on + the entire thing globally. + + - 'serial': The polygon is built in steps by adding one + polygon at a time and computing the intersection at + each step. + + This option is provided because one may be faster than the + other depending on context, but it primarily exposed for + testing reasons. Both modes should theoretically provide + equivalent results. + + Returns + ------- + polygon : `SphericalPolygon` object + """ + assert len(polygons) + for polygon in polygons: + assert isinstance(polygon, SphericalPolygon) + + # for i in range(len(polygons)): + # polyA = polygons[i] + # for j in range(i + 1, len(polygons)): + # polyB = polygons[j] + # if not polyA.intersects_poly(polyB): + # return cls([], [0, 0, 0]) + + from . import graph + + if method.lower() == 'parallel': + g = graph.Graph(polygons) + polygon = g.intersection() + inside = cls._find_new_inside(polygon, polygons) + return cls(polygon, inside) + elif method.lower() == 'serial': + result = copy(polygons[0]) + for polygon in polygons[1:]: + result = result.intersection(polygon) + # If we have a null intersection already, we don't + # need to go any further. + if len(result._points) < 3: + return result + return result + else: + raise ValueError("method must be 'parallel' or 'serial'") + + def overlap(self, other): + r""" + Returns the fraction of *self* that is overlapped by *other*. + + Let *self* be *a* and *other* be *b*, then the overlap is + defined as: + + .. math:: + + \frac{S_a}{S_{a \cap b}} + + Parameters + ---------- + other : `SphericalPolygon` + + Returns + ------- + frac : float + The fraction of *self* that is overlapped by *other*. + """ + s1 = self.area() + intersection = self.intersection(other) + s2 = intersection.area() + return s2 / s1 + + def draw(self, m, **plot_args): + """ + Draws the polygon in a matplotlib.Basemap axes. + + Parameters + ---------- + m : Basemap axes object + + **plot_args : Any plot arguments to pass to basemap + """ + if not len(self._points): + return + if not len(plot_args): + plot_args = {'color': 'blue'} + points = self._points + + for A, B in zip(points[0:-1], points[1:]): + length = great_circle_arc.length(A, B, degrees=True) + if not np.isfinite(length): + length = 2 + interpolated = great_circle_arc.interpolate(A, B, length * 4) + ra, dec = vector.vector_to_radec( + interpolated[:, 0], interpolated[:, 1], interpolated[:, 2], + degrees=True) + for r0, d0, r1, d1 in zip(ra[0:-1], dec[0:-1], ra[1:], dec[1:]): + m.drawgreatcircle(r0, d0, r1, d1, **plot_args) + + ra, dec = vector.vector_to_radec( + *self._inside, degrees=True) + x, y = m(ra, dec) + m.scatter(x, y, 1, **plot_args) |