diff options
Diffstat (limited to 'lib')
| -rw-r--r-- | lib/graph.py | 58 | ||||
| -rw-r--r-- | lib/great_circle_arc.py | 262 | ||||
| -rw-r--r-- | lib/polygon.py | 13 | ||||
| -rw-r--r-- | lib/skyline.py | 108 |
4 files changed, 232 insertions, 209 deletions
diff --git a/lib/graph.py b/lib/graph.py index 491f2fa..a32c3be 100644 --- a/lib/graph.py +++ b/lib/graph.py @@ -60,14 +60,14 @@ class Graph: .. note:: This class is not meant to be used directly. Instead, use - `sphere.polygon.SphericalPolygon.union` and - `sphere.polygon.SphericalPolygon.intersection`. + `~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. + """ + A `~Graph.Node` represents a single point, connected by an arbitrary + number of `~Graph.Edge` objects to other `~Graph.Node` objects. """ def __init__(self, point, source_polygons=[]): """ @@ -91,12 +91,12 @@ class Graph: Parameters ---------- - edge : `Edge` instance + edge : `~Graph.Edge` instance The edge to follow away from. Returns ------- - other : `Edge` instance + other : `~Graph.Edge` instance The other edge. """ edges = list(self._edges) @@ -109,11 +109,14 @@ class Graph: def equals(self, other, thres=2e-8): """ Returns `True` if the location of this and the *other* - `Node` are the same. + `~Graph.Node` are the same. Parameters ---------- - thres: float + other : `~Graph.Node` instance + The other node. + + thres : float If difference is smaller than this, points are equal. The default value of 2e-8 radians is set based on emphirical test cases. Relative threshold based on @@ -126,14 +129,14 @@ class Graph: class Edge: """ - An `Edge` represents a connection between exactly two `Node` - objects. This `Edge` class has no direction. + An `~Graph.Edge` represents a connection between exactly two + `~Graph.Node` objects. This `~Graph.Edge` class has no direction. """ def __init__(self, A, B, source_polygons=[]): """ Parameters ---------- - A, B : `Node` instances + A, B : `~Graph.Node` instances source_polygon : `~sphere.polygon.SphericalPolygon` instance, optional The polygon this edge came from. Used for bookkeeping. @@ -154,11 +157,11 @@ class Graph: Parameters ---------- - node : `Node` instance + node : `~Graph.Node` instance Returns ------- - other : `Node` instance + other : `~Graph.Node` instance """ nodes = self._nodes try: @@ -168,11 +171,11 @@ class Graph: def equals(self, other): """ - Returns True if the other edge is between the same two nodes. + Returns `True` if the other edge is between the same two nodes. Parameters ---------- - other : `Edge` instance + other : `~Graph.Edge` instance Returns ------- @@ -261,7 +264,7 @@ class Graph: Returns ------- - node : `Node` instance + node : `~Graph.Node` instance The new node """ new_node = self.Node(point, source_polygons) @@ -285,7 +288,7 @@ class Graph: Parameters ---------- - node : `Node` instance + node : `~Graph.Node` instance """ assert node in self._nodes @@ -304,14 +307,14 @@ class Graph: Parameters ---------- - A, B : `Node` instances + A, B : `~Graph.Node` instances source_polygons : `~sphere.polygon.SphericalPolygon` instance, optional The polygon(s) this edge came from. Used for bookkeeping. Returns ------- - edge : `Edge` instance + edge : `~Graph.Edge` instance The new edge """ assert A in self._nodes @@ -338,7 +341,7 @@ class Graph: Parameters ---------- - edge : `Edge` instance + edge : `~Graph.Edge` instance """ assert edge in self._edges @@ -354,21 +357,22 @@ class Graph: 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. + Splits an `~Graph.Edge` *edge* at `~Graph.Node` *node*, removing + *edge* and replacing it with two new `~Graph.Edge` instances. + It is intended that *E* is along the original edge, but that is + not enforced. Parameters ---------- - edge : `Edge` instance + edge : `~Graph.Edge` instance The edge to split - node : `Node` instance + node : `~Graph.Node` instance The node to insert Returns ------- - edgeA, edgeB : `Edge` instances + edgeA, edgeB : `~Graph.Edge` instances The two new edges on either side of *node*. """ assert edge in self._edges diff --git a/lib/great_circle_arc.py b/lib/great_circle_arc.py index 399cb3d..d5e78f2 100644 --- a/lib/great_circle_arc.py +++ b/lib/great_circle_arc.py @@ -53,138 +53,145 @@ try: except ImportError: HAS_C_UFUNCS = False +if HAS_C_UFUNCS: + inner1d = math_util.inner1d +else: + from numpy.core.umath_tests import inner1d + __all__ = ['angle', 'intersection', 'intersects', 'length', 'midpoint', 'interpolate'] -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. + """ + if HAS_C_UFUNCS: + return math_util.cross(a, b) + + cp = np.empty(np.broadcast(a, b).shape) + aT = a.T + bT = b.T + cpT = cp.T + + cpT[0] = aT[1]*bT[2] - aT[2]*bT[1] + cpT[1] = aT[2]*bT[0] - aT[0]*bT[2] + cpT[2] = aT[0]*bT[1] - aT[1]*bT[0] + + return cp + +def _cross_and_normalize(A, B): + if HAS_C_UFUNCS: + return math_util.cross_and_norm(A, B) + + T = _fast_cross(A, B) + # Normalization + l = np.sqrt(np.sum(T ** 2, axis=-1)) + l = np.expand_dims(l, 2) + # Might get some divide-by-zeros, but we don't care + with np.errstate(invalid='ignore'): + TN = T / l + return TN - def _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): - r""" - Returns the point of intersection between two great circle arcs. - The arcs are defined between the points *AB* and *CD*. Either *A* - and *B* or *C* and *D* may be arrays of points, but not both. - - Parameters - ---------- - A, B : (*x*, *y*, *z*) triples or Nx3 arrays of triples - Endpoints of the first great circle arc. - - C, D : (*x*, *y*, *z*) triples or Nx3 arrays of triples - Endpoints of the second great circle arc. - - Returns - ------- - T : (*x*, *y*, *z*) triples or Nx3 arrays of triples - If the given arcs intersect, the intersection is returned. If - the arcs do not intersect, the triple is set to all NaNs. - - Notes - ----- - The basic intersection is computed using linear algebra as follows - [1]_: - - .. math:: - - T = \lVert(A × B) × (C × D)\rVert - - To determine the correct sign (i.e. hemisphere) of the - intersection, the following four values are computed: - - .. math:: - - s_1 = ((A × B) × A) · T - - s_2 = (B × (A × B)) · T - - s_3 = ((C × D) × C) · T - - s_4 = (D × (C × D)) · T - - For :math:`s_n`, if all positive :math:`T` is returned as-is. If - all negative, :math:`T` is multiplied by :math:`-1`. Otherwise - the intersection does not exist and is undefined. - - References - ---------- - - .. [1] Method explained in an `e-mail - <http://www.mathworks.com/matlabcentral/newsreader/view_thread/276271>`_ - by Roger Stafford. - - http://www.mathworks.com/matlabcentral/newsreader/view_thread/276271 - """ - 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 intersection(A, B, C, D): + r""" + Returns the point of intersection between two great circle arcs. + The arcs are defined between the points *AB* and *CD*. Either *A* + and *B* or *C* and *D* may be arrays of points, but not both. + + Parameters + ---------- + A, B : (*x*, *y*, *z*) triples or Nx3 arrays of triples + Endpoints of the first great circle arc. + + C, D : (*x*, *y*, *z*) triples or Nx3 arrays of triples + Endpoints of the second great circle arc. + + Returns + ------- + T : (*x*, *y*, *z*) triples or Nx3 arrays of triples + If the given arcs intersect, the intersection is returned. If + the arcs do not intersect, the triple is set to all NaNs. + + Notes + ----- + The basic intersection is computed using linear algebra as follows + [1]_: + + .. math:: + + T = \lVert(A × B) × (C × D)\rVert + + To determine the correct sign (i.e. hemisphere) of the + intersection, the following four values are computed: + + .. math:: + + s_1 = ((A × B) × A) · T + + s_2 = (B × (A × B)) · T + + s_3 = ((C × D) × C) · T + + s_4 = (D × (C × D)) · T + + For :math:`s_n`, if all positive :math:`T` is returned as-is. If + all negative, :math:`T` is multiplied by :math:`-1`. Otherwise + the intersection does not exist and is undefined. + + References + ---------- + + .. [1] Method explained in an `e-mail + <http://www.mathworks.com/matlabcentral/newsreader/view_thread/276271>`_ + by Roger Stafford. + + http://www.mathworks.com/matlabcentral/newsreader/view_thread/276271 + """ + if HAS_C_UFUNCS: + return math_util.intersection(A, B, C, D) + + A = np.asanyarray(A) + B = np.asanyarray(B) + C = np.asanyarray(C) + D = np.asanyarray(D) + + A, B = np.broadcast_arrays(A, B) + C, D = np.broadcast_arrays(C, D) + + ABX = _fast_cross(A, B) + CDX = _fast_cross(C, D) + T = _cross_and_normalize(ABX, CDX) + T_ndim = len(T.shape) + + if T_ndim > 1: + s = np.zeros(T.shape[0]) + else: + s = np.zeros(1) + s += np.sign(inner1d(_fast_cross(ABX, A), T)) + s += np.sign(inner1d(_fast_cross(B, ABX), T)) + s += np.sign(inner1d(_fast_cross(CDX, C), T)) + s += np.sign(inner1d(_fast_cross(D, CDX), T)) + if T_ndim > 1: + s = np.expand_dims(s, 2) + + cross = np.where(s == -4, -T, np.where(s == 4, T, np.nan)) + + # If they share a common point, it's not an intersection. This + # gets around some rounding-error/numerical problems with the + # above. + equals = (np.all(A == C, axis = -1) | + np.all(A == D, axis = -1) | + np.all(B == C, axis = -1) | + np.all(B == D, axis = -1)) + + equals = np.expand_dims(equals, 2) + + return np.where(equals, np.nan, cross) def length(A, B, degrees=True): @@ -261,7 +268,8 @@ def angle(A, B, C, degrees=True): Parameters ---------- - A, B, C : (*x*, *y*, *z*) triples or Nx3 arrays of triples. + A, B, C : (*x*, *y*, *z*) triples or Nx3 arrays of triples + Points on sphere. degrees : bool, optional If `True` (default) the result is returned in decimal degrees, diff --git a/lib/polygon.py b/lib/polygon.py index 8047ab1..fd67042 100644 --- a/lib/polygon.py +++ b/lib/polygon.py @@ -115,12 +115,14 @@ class SphericalPolygon(object): return self._inside def to_radec(self): - """Convert `SphericalPolygon` footprint to RA and DEC. + """ + Convert `SphericalPolygon` footprint to RA and DEC. Returns ------- - List of RA and DEC pairs (in degrees) corresponding to - `self.points`. + zip(ra,dec) : list of tuples + List of (*ra*, *dec*) pairs in degrees corresponding + to `points`. """ ra, dec = vector.vector_to_radec(self.points[:,0], self.points[:,1], self.points[:,2], degrees=True) @@ -314,7 +316,7 @@ class SphericalPolygon(object): Parameters ---------- - preserve_order: bool + preserve_order : bool Preserve original order? If `True`, polygon is rotated around min point. If `False`, all points are sorted in ascending order. @@ -361,7 +363,8 @@ class SphericalPolygon(object): Returns ------- - `True` or `False` + is_eq : bool + `True` or `False`. """ self_n = len(self.points) diff --git a/lib/skyline.py b/lib/skyline.py index 113192e..db4d6b8 100644 --- a/lib/skyline.py +++ b/lib/skyline.py @@ -33,11 +33,9 @@ # OF THE POSSIBILITY OF SUCH DAMAGE. """ -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:: +generalized steps: #. Initialize `SkyLine` objects for each input image. This object would be the union of all the input @@ -72,16 +70,6 @@ This process will work reasonably fast as most operations are performed using the `SkyLine` objects and WCS information solely, not image data itself. -.. note:: Requires Python 2.7 or later. - -:Authors: Pey Lian Lim, Warren Hack, Michael Droettboom - -:Organization: Space Telescope Science Institute - -:History: - * 2012-05-25 PLL started coding. Original class - structure by WJH. Parent class by MD. - """ from __future__ import division, print_function, absolute_import @@ -96,26 +84,28 @@ from stwcs.distortion.utils import output_wcs # LOCAL from .polygon import SphericalPolygon -__all__ = ['SkyLine'] +__all__ = ['SkyLineMember', 'SkyLine'] __version__ = '0.3a' __vdate__ = '20-Jun-2012' class SkyLineMember(object): + """ + Container for `SkyLine` members with these attributes: + + * `fname`: Image name (with path if given) + * `ext`: Extension read + * `wcs`: `HSTWCS` object the data data + * `polygon`: `~sphere.polygon.SphericalPolygon` object of the data + """ def __init__(self, fname, ext): """ - Container for `SkyLine` members. - - Given FITS image and extension, will store its - `SphericalPolygon` instance from WCS under - `polygon`. - Parameters ---------- - fname: str + fname : str FITS image. - ext: int + ext : int Image extension. """ @@ -145,17 +135,22 @@ class SkyLineMember(object): return self._polygon class SkyLine(object): + """ + Manage outlines on the sky. + + Each `SkyLine` has a list of `~SkyLine.members` and + a composite `~SkyLine.polygon` with all the + functionalities of `~sphere.polygon.SphericalPolygon`. + """ def __init__(self, fname, extname='SCI'): """ - Initialize `SkyLine` object instance. - Parameters ---------- - fname: str + fname : str FITS image. `None` to create empty `SkyLine`. - extname: str + extname : str EXTNAME to use. SCI is recommended for normal HST images. PRIMARY if image is single ext. @@ -179,7 +174,7 @@ class SkyLine(object): self.polygon = SphericalPolygon([]) def __getattr__(self, what): - """Control attribute access to `SphericalPolygon`.""" + """Control attribute access to `~sphere.polygon.SphericalPolygon`.""" if what in ('from_radec', 'from_cone', 'from_wcs', 'multi_union', 'multi_intersection', '_find_new_inside',): @@ -197,23 +192,31 @@ class SkyLine(object): @property def polygon(self): - """`SphericalPolygon` portion of `SkyLine`.""" + """ + `~sphere.polygon.SphericalPolygon` portion of `SkyLine` + that contains the composite skyline from `members` + belonging to *self*. + + """ return self._polygon @polygon.setter def polygon(self, value): - """Deep copy a `SphericalPolygon`.""" assert isinstance(value, SphericalPolygon) - self._polygon = copy(value) + self._polygon = copy(value) # Deep copy @property def members(self): - """List of `SkyLineMember` objects.""" + """ + List of `SkyLineMember` objects that belong to *self*. + Duplicate members are discarded. Members are kept in + the order of their additions to *self*. + + """ return self._members @members.setter def members(self, values): - """Make sure `SkyLineMember` entries are unique.""" self._members = [] # Not using set to preserve order @@ -226,8 +229,8 @@ class SkyLine(object): def to_wcs(self): """ - Combine WCS objects from all `members` and return a - new WCS object. If no `members`, return `None`. + Combine `HSTWCS` objects from all `members` and return + a new `HSTWCS` object. If no `members`, return `None`. """ wcs = None @@ -244,15 +247,15 @@ class SkyLine(object): Parameters ---------- - self: obj + self : obj `SkyLine` instance. - given_members: list + given_members : list List of `SkyLineMember` to consider. Returns ------- - new_members: list + new_members : list List of `SkyLineMember` belonging to *self*. """ @@ -310,16 +313,16 @@ class SkyLine(object): Parameters ---------- - skylines: list + skylines : list A list of `SkyLine` instances. Returns ------- - max_skyline: `SkyLine` instance or `None` + max_skyline : `SkyLine` instance or `None` `SkyLine` that overlaps the most or `None` if no overlap found. This is *not* a copy. - max_overlap_area: float + max_overlap_area : float Area of intersection. """ @@ -342,13 +345,13 @@ class SkyLine(object): Parameters ---------- - skylines: list + skylines : list A list of `SkyLine` instances. Returns ------- - max_pair: tuple - Pair of `SkyLine` instances with max overlap + max_pair : tuple + Pair of `SkyLine` objects with max overlap among given *skylines*. If no overlap found, return `None`. These are *not* copies. @@ -371,24 +374,29 @@ class SkyLine(object): """ Mosaic all overlapping *skylines*. + A pair of skylines with the most overlap is used as + a starting point. Then a skyline that overlaps the + most with the mosaic is used, and so forth until no + overlapping skyline is found. + Parameters ---------- - skylines: list - A list of `SkyLine` instances. + skylines : list + A list of `SkyLine` objects. Returns ------- - mosaic: `SkyLine` instance or `None` + mosaic : `SkyLine` instance or `None` Union of all overlapping *skylines*, or `None` if no overlap found. - included: list - List of image names added to `mosaic` in the order + included : list + List of image names added to mosaic in the order of addition. - excluded: list + excluded : list List of image names excluded because they do not - overlap with `mosaic`. + overlap with mosaic. """ out_order = [] |