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

1 files changed, 318 insertions, 0 deletions

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
+           <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 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