summaryrefslogtreecommitdiff
path: root/stsci
diff options
context:
space:
mode:
Diffstat (limited to 'stsci')
-rw-r--r--stsci/sphere/graph.py14
-rw-r--r--stsci/sphere/polygon.py31
-rw-r--r--stsci/sphere/test/test.py5
-rw-r--r--stsci/sphere/test/test_union.py2
-rw-r--r--stsci/sphere/vector.py62
5 files changed, 60 insertions, 54 deletions
diff --git a/stsci/sphere/graph.py b/stsci/sphere/graph.py
index 0fbc9a7..bae1c60 100644
--- a/stsci/sphere/graph.py
+++ b/stsci/sphere/graph.py
@@ -270,11 +270,17 @@ class Graph:
# Don't add nodes that already exist. Update the existing
# node's source_polygons list to include the new polygon.
- point = vector.normalize_vector(*point)
+ point = vector.normalize_vector(point)
- # TODO: Vectorize this
- for node in self._nodes:
- if np.all(np.abs(node._point - point) < 2 ** -32):
+ if len(self._nodes):
+ nodes = list(self._nodes)
+ node_array = np.array([node._point for node in nodes])
+
+ diff = np.all(np.abs(point - node_array) < 2 ** -32, axis=-1)
+
+ indices = np.nonzero(diff)[0]
+ if len(indices):
+ node = nodes[indices[0]]
node._source_polygons.update(source_polygons)
return node
diff --git a/stsci/sphere/polygon.py b/stsci/sphere/polygon.py
index f7c44be..5167171 100644
--- a/stsci/sphere/polygon.py
+++ b/stsci/sphere/polygon.py
@@ -168,11 +168,11 @@ class SphericalPolygon(object):
# Convert to Cartesian
x, y, z = vector.radec_to_vector(ra, dec, degrees=degrees)
+ points = np.dstack((x, y, z))[0]
+
if center is None:
- xc = x.mean()
- yc = y.mean()
- zc = z.mean()
- center = vector.normalize_vector(xc, yc, zc)
+ center = points.mean(axis=0)
+ vector.normalize_vector(center, output=center)
else:
center = vector.radec_to_vector(*center, degrees=degrees)
@@ -556,24 +556,21 @@ class SphericalPolygon(object):
# Rotate polygon so that center of polygon is at north pole
centroid = np.mean(points[:-1], axis=0)
- centroid = vector.normalize_vector(*centroid)
+ 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)
+ vector.normalize_vector(points, output=points)
- X = []
- Y = []
+ XYs = []
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)
+ vector.normalize_vector(interp, output=interp)
+ XY = vector.equal_area_proj(interp)
+ XYs.append(XY)
+
+ XY = np.vstack(XYs)
+ X = XY[..., 0]
+ Y = XY[..., 1]
return np.abs(np.sum(X[:-1] * Y[1:] - X[1:] * Y[:-1]) * 0.5 * np.pi)
diff --git a/stsci/sphere/test/test.py b/stsci/sphere/test/test.py
index 8e6e270..5ab7bbb 100644
--- a/stsci/sphere/test/test.py
+++ b/stsci/sphere/test/test.py
@@ -19,8 +19,9 @@ 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)
+ xyz = np.dstack((x.flatten(), y.flatten(), z.flatten()))[0]
+ xyzn = vector.normalize_vector(xyz)
+ l = np.sqrt(np.sum(xyzn * xyzn, axis=-1))
assert_almost_equal(l, 1.0)
diff --git a/stsci/sphere/test/test_union.py b/stsci/sphere/test/test_union.py
index 60ea572..a6e5793 100644
--- a/stsci/sphere/test/test_union.py
+++ b/stsci/sphere/test/test_union.py
@@ -209,7 +209,7 @@ def test_difficult_unions():
poly = polygon.SphericalPolygon(points, inside)
polys.append(poly)
- polygon.SphericalPolygon.multi_union(polys)
+ polygon.SphericalPolygon.multi_union(polys[:len(polys)/2])
if __name__ == '__main__':
diff --git a/stsci/sphere/vector.py b/stsci/sphere/vector.py
index c6887c4..bb444e6 100644
--- a/stsci/sphere/vector.py
+++ b/stsci/sphere/vector.py
@@ -42,6 +42,12 @@ from __future__ import unicode_literals
# THIRD-PARTY
import numpy as np
+try:
+ from . import math_util
+ HAS_C_UFUNCS = True
+except ImportError:
+ HAS_C_UFUNCS = False
+
__all__ = ['radec_to_vector', 'vector_to_radec', 'normalize_vector',
'rotate_around']
@@ -121,9 +127,9 @@ def vector_to_radec(x, y, z, degrees=True):
\delta = \arctan2(z, \sqrt{x^2 + y^2})
"""
- x = np.asanyarray(x)
- y = np.asanyarray(y)
- z = np.asanyarray(z)
+ x = np.asanyarray(x, dtype=np.float64)
+ y = np.asanyarray(y, dtype=np.float64)
+ z = np.asanyarray(z, dtype=np.float64)
result = (
np.arctan2(y, x),
@@ -135,39 +141,37 @@ def vector_to_radec(x, y, z, degrees=True):
return result
-def normalize_vector(x, y, z, inplace=False):
+def normalize_vector(xyz, output=None):
r"""
Normalizes a vector so it falls on the unit sphere.
- *x*, *y*, *z* may be scalars or 1-D arrays
-
Parameters
----------
- x, y, z : scalars or 1-D arrays of the same length
+ xyz : Nx3 array of vectors
The input vectors
- inplace : bool, optional
- When `True`, the original arrays will be normalized in place,
- otherwise a normalized copy is returned.
+ output : Nx3 array of vectors, optional
+ The array to store the results in. If `None`, a new array
+ will be created and returned.
Returns
-------
- X, Y, Z : scalars of 1-D arrays of the same length
- The normalized output vectors
+ output : Nx3 array of vectors
"""
- x = np.asanyarray(x)
- y = np.asanyarray(y)
- z = np.asanyarray(z)
+ xyz = np.asanyarray(xyz, dtype=np.float64)
- l = np.sqrt(x ** 2 + y ** 2 + z ** 2)
+ if output is None:
+ output = np.empty(xyz.shape, dtype=np.float64)
- if inplace:
- x /= l
- y /= l
- z /= l
- return (x, y, z)
- else:
- return (x / l, y / l, z / l)
+ if HAS_C_UFUNCS:
+ math_util.normalize(xyz, output)
+ return output
+
+ l = np.sqrt(np.sum(xyz * xyz, axis=-1))
+
+ output = xyz / np.expand_dims(l, 2)
+
+ return output
def rotate_around(x, y, z, u, v, w, theta, degrees=True):
@@ -212,19 +216,19 @@ def rotate_around(x, y, z, u, v, w, theta, degrees=True):
return X, Y, Z
-def equal_area_proj(x, y, z):
+def equal_area_proj(points):
"""
Transform the coordinates to a 2-dimensional plane using the
Lambert azimuthal equal-area projection.
Parameters
----------
- x, y, z : scalars or 1-D arrays
+ points : Nx3 array of vectors
The input vectors
Returns
-------
- X, Y : tuple of scalars or arrays of the same length
+ output : Nx2 array of points
Notes
-----
@@ -237,7 +241,5 @@ def equal_area_proj(x, y, z):
Y = \sqrt{\frac{2}{1-z}}y
"""
- scale = np.sqrt(2.0 / (1.0 - z))
- X = scale * x
- Y = scale * y
- return X, Y
+ scale = np.sqrt(2.0 / (1.0 - points[..., 2]))
+ return np.expand_dims(scale, 2) * points[:, 0:2]
generated by cgit (git 2.34.1) at 2025-04-29 00:41:00 -0400