summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/math_util.c90
-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
6 files changed, 124 insertions, 80 deletions
diff --git a/src/math_util.c b/src/math_util.c
index 9dbb32a..2aae11a 100644
--- a/src/math_util.c
+++ b/src/math_util.c
@@ -49,6 +49,9 @@ typedef npy_intp intp;
INIT_OUTER_LOOP_4 \
intp s4 = *steps++;
+#define BEGIN_OUTER_LOOP_2 \
+ for (N_ = 0; N_ < dN; N_++, args[0] += s0, args[1] += s1) {
+
#define BEGIN_OUTER_LOOP_3 \
for (N_ = 0; N_ < dN; N_++, args[0] += s0, args[1] += s1, args[2] += s2) {
@@ -97,13 +100,6 @@ save_pointl(const long double* in, char *out, const intp s) {
}
static inline void
-cross(const double *A, const double *B, double *C) {
- C[0] = A[1]*B[2] - A[2]*B[1];
- C[1] = A[2]*B[0] - A[0]*B[2];
- C[2] = A[0]*B[1] - A[1]*B[0];
-}
-
-static inline void
crossl(const long double *A, const long double *B, long double *C) {
C[0] = A[1]*B[2] - A[2]*B[1];
C[1] = A[2]*B[0] - A[0]*B[2];
@@ -111,13 +107,17 @@ crossl(const long double *A, const long double *B, long double *C) {
}
static inline void
-normalize(double *A) {
+normalize_output(long double *A, double *B) {
double l = A[0]*A[0] + A[1]*A[1] + A[2]*A[2];
if (l != 1.0) {
l = sqrt(l);
- A[0] /= l;
- A[1] /= l;
- A[2] /= l;
+ B[0] = A[0] / l;
+ B[1] = A[1] / l;
+ B[2] = A[2] / l;
+ } else {
+ B[0] = A[0];
+ B[1] = A[1];
+ B[2] = A[2];
}
}
@@ -207,6 +207,35 @@ static void * inner1d_data[] = { (void *)NULL };
static char inner1d_signatures[] = { PyArray_DOUBLE, PyArray_DOUBLE, PyArray_DOUBLE };
/*///////////////////////////////////////////////////////////////////////////
+ normalize
+*/
+
+char *normalize_signature = "(i)->(i)";
+
+static void
+DOUBLE_normalize(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+ long double IN[3];
+ double OUT[3];
+
+ INIT_OUTER_LOOP_2
+ intp is1=steps[0], is2=steps[1];
+ BEGIN_OUTER_LOOP_2
+ char *ip1=args[0], *op=args[1];
+
+ load_pointl(ip1, is1, IN);
+
+ normalize_output(IN, OUT);
+
+ save_point(OUT, op, is2);
+ END_OUTER_LOOP
+}
+
+static PyUFuncGenericFunction normalize_functions[] = { DOUBLE_normalize };
+static void * normalize_data[] = { (void *)NULL };
+static char normalize_signatures[] = { PyArray_DOUBLE, PyArray_DOUBLE };
+
+/*///////////////////////////////////////////////////////////////////////////
cross
*/
@@ -215,21 +244,21 @@ char *cross_signature = "(i),(i)->(i)";
static void
DOUBLE_cross(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
{
- double A[3];
- double B[3];
- double C[3];
+ long double A[3];
+ long double B[3];
+ long double C[3];
INIT_OUTER_LOOP_3
intp is1=steps[0], is2=steps[1], is3=steps[2];
BEGIN_OUTER_LOOP_3
char *ip1=args[0], *ip2=args[1], *op=args[2];
- load_point(ip1, is1, A);
- load_point(ip2, is2, B);
+ load_pointl(ip1, is1, A);
+ load_pointl(ip2, is2, B);
- cross(A, B, C);
+ crossl(A, B, C);
- save_point(C, op, is3);
+ save_pointl(C, op, is3);
END_OUTER_LOOP
}
@@ -250,22 +279,22 @@ char *cross_and_norm_signature = "(i),(i)->(i)";
static void
DOUBLE_cross_and_norm(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
{
- double A[3];
- double B[3];
- double C[3];
+ long double A[3];
+ long double B[3];
+ long double C[3];
INIT_OUTER_LOOP_3
intp is1=steps[0], is2=steps[1], is3=steps[2];
BEGIN_OUTER_LOOP_3
char *ip1=args[0], *ip2=args[1], *op=args[2];
- load_point(ip1, is1, A);
- load_point(ip2, is2, B);
+ load_pointl(ip1, is1, A);
+ load_pointl(ip2, is2, B);
- cross(A, B, C);
- normalize(C);
+ crossl(A, B, C);
+ normalizel(C);
- save_point(C, op, is3);
+ save_pointl(C, op, is3);
END_OUTER_LOOP
}
@@ -455,6 +484,15 @@ addUfuncs(PyObject *dictionary) {
Py_DECREF(f);
f = PyUFunc_FromFuncAndDataAndSignature(
+ normalize_functions, normalize_data, normalize_signatures, 2, 1, 1,
+ PyUFunc_None, "normalize",
+ "inner on the last dimension and broadcast on the rest \n" \
+ " \"(i)->(i)\" \n",
+ 0, normalize_signature);
+ PyDict_SetItemString(dictionary, "normalize", f);
+ Py_DECREF(f);
+
+ f = PyUFunc_FromFuncAndDataAndSignature(
cross_functions, cross_data, cross_signatures, 2, 2, 1,
PyUFunc_None, "cross",
"cross product of 3-vectors only \n" \
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-28 18:07:22 -0400