From dada4d0162050f075ca242157d5070dc2f463b33 Mon Sep 17 00:00:00 2001 From: Leo Singer Date: Tue, 30 Oct 2018 17:36:38 -0400 Subject: [PATCH] WIP: Rewrite core in C using Numpy ufuncs As a result, broadcasting across nside and ipix is supported automatically. Fixes #72, closes #75. --- astropy_healpix/_core.c | 301 +++++++++++++ astropy_healpix/core.py | 153 +++---- astropy_healpix/core_cython.pxd | 44 -- astropy_healpix/core_cython.pyx | 492 ---------------------- astropy_healpix/setup_package.py | 4 +- astropy_healpix/tests/test_core_cython.py | 129 ------ cextern/astrometry.net/healpix-utils.h | 2 +- 7 files changed, 362 insertions(+), 763 deletions(-) create mode 100644 astropy_healpix/_core.c delete mode 100644 astropy_healpix/core_cython.pxd delete mode 100644 astropy_healpix/core_cython.pyx delete mode 100644 astropy_healpix/tests/test_core_cython.py diff --git a/astropy_healpix/_core.c b/astropy_healpix/_core.c new file mode 100644 index 0000000..115ea32 --- /dev/null +++ b/astropy_healpix/_core.c @@ -0,0 +1,301 @@ +#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION + +#include +#include +#include +#include "healpix.h" +#include "healpix-utils.h" +#include "interpolation.h" + + +typedef struct { + int64_t (*order_to_xy)(int64_t, int); + int64_t (*xy_to_order)(int64_t, int); +} order_funcs; + + +static void healpix_to_lonlat_loop( + char **args, npy_intp *dimensions, npy_intp *steps, void *data) +{ + order_funcs *funcs = data; + + for (npy_intp i = 0; i < dimensions[0]; i ++) + { + int64_t pixel = *(int64_t *) &args[0][i * steps[0]]; + int nside = *(int *) &args[1][i * steps[1]]; + double dx = *(double *) &args[2][i * steps[2]]; + double dy = *(double *) &args[3][i * steps[3]]; + double *lon = (double *) &args[4][i * steps[4]]; + double *lat = (double *) &args[5][i * steps[5]]; + + int64_t xy = funcs->order_to_xy(pixel, nside); + if (xy >= 0) + healpixl_to_radec(xy, nside, dx, dy, lon, lat); + else + *lon = *lat = NPY_NAN; + } +} + + +static void lonlat_to_healpix_loop( + char **args, npy_intp *dimensions, npy_intp *steps, void *data) +{ + order_funcs *funcs = data; + + for (npy_intp i = 0; i < dimensions[0]; i ++) + { + double lon = *(double *) &args[0][i * steps[0]]; + double lat = *(double *) &args[1][i * steps[1]]; + int nside = *(int *) &args[2][i * steps[2]]; + int64_t *pixel = (int64_t *) &args[3][i * steps[3]]; + double *dx = (double *) &args[4][i * steps[4]]; + double *dy = (double *) &args[5][i * steps[5]]; + + int64_t xy = radec_to_healpixlf(lon, lat, nside, dx, dy); + if (xy >= 0) + *pixel = funcs->xy_to_order(xy, nside); + else { + *pixel = -1; + *dx = *dy = NPY_NAN; + } + } +} + + +static void nested_to_ring_loop( + char **args, npy_intp *dimensions, npy_intp *steps, void *data) +{ + for (npy_intp i = 0; i < dimensions[0]; i ++) + { + int64_t nested = *(int64_t *) &args[0][i * steps[0]]; + int nside = *(int *) &args[1][i * steps[1]]; + int64_t *ring = (int64_t *) &args[2][i * steps[2]]; + + int64_t xy = healpixl_nested_to_xy(nested, nside); + if (xy >= 0) + *ring = healpixl_xy_to_ring(xy, nside); + else + *ring = -1; + } +} + + +static void ring_to_nested_loop( + char **args, npy_intp *dimensions, npy_intp *steps, void *data) +{ + for (npy_intp i = 0; i < dimensions[0]; i ++) + { + int64_t ring = *(int64_t *) &args[0][i * steps[0]]; + int nside = *(int *) &args[1][i * steps[1]]; + int64_t *nested = (int64_t *) &args[2][i * steps[2]]; + + int64_t xy = healpixl_ring_to_xy(ring, nside); + if (xy >= 0) + *nested = healpixl_xy_to_nested(xy, nside); + else + *nested = -1; + } +} + + +static void bilinear_interpolation_weights_loop( + char **args, npy_intp *dimensions, npy_intp *steps, void *data) +{ + for (npy_intp i = 0; i < dimensions[0]; i ++) + { + double lon = *(double *) &args[0][i * steps[0]]; + double lat = *(double *) &args[1][i * steps[1]]; + int nside = *(int *) &args[2][i * steps[2]]; + int64_t indices[4]; + double weights[4]; + + interpolate_weights(lon, lat, indices, weights, nside); + for (int j = 0; j < 4; j ++) + { + *(int64_t *) &args[3 + j][i * steps[3 + j]] = indices[j]; + *(double *) &args[7 + j][i * steps[7 + j]] = weights[j]; + } + } +} + + +static void neighbours_loop( + char **args, npy_intp *dimensions, npy_intp *steps, void *data) +{ + order_funcs *funcs = data; + + for (npy_intp i = 0; i < dimensions[0]; i ++) + { + int64_t pixel = *(int64_t *) &args[0][i * steps[0]]; + int nside = *(int *) &args[1][i * steps[1]]; + int64_t neighbours[] = {-1, -1, -1, -1, -1, -1, -1, -1}; + + int64_t xy = funcs->order_to_xy(pixel, nside); + if (xy >= 0) + healpixl_get_neighbours(xy, neighbours, nside); + + for (int j = 0; j < 8; j ++) + { + xy = neighbours[(4 - j) % 8]; + if (xy >= 0) + pixel = funcs->xy_to_order(xy, nside); + else + pixel = -1; + *(int64_t *) &args[2 + j][i * steps[2 + j]] = pixel; + } + } +} + + +static PyObject *healpix_cone_search( + PyObject *self, PyObject *args, PyObject *kwargs) +{ + PyObject *result = NULL; + static const char *kws[] = {"lon", "lat", "radius", "nside", "order", NULL}; + double lon, lat, radius; + int nside; + char *order; + + if (!PyArg_ParseTupleAndKeywords( + args, kwargs, "dddis", kws, &lon, &lat, &radius, &nside, &order)) + return NULL; + + int64_t *indices; + int64_t n_indices = healpix_rangesearch_radec_simple( + lon, lat, radius, nside, 0, &indices); + if (!indices) + { + PyErr_SetString( + PyExc_RuntimeError, "healpix_rangesearch_radec_simple failed"); + goto done; + } + + npy_intp dims[] = {n_indices}; + result = PyArray_SimpleNew(1, dims, NPY_INT64); + if (!result) + goto done; + + int64_t *result_data = PyArray_DATA((PyArrayObject *) result); + + if (strcmp(order, "nested") == 0) + { + for (int i = 0; i < n_indices; i ++) + result_data[i] = healpixl_xy_to_nested(indices[i], nside); + } else { + for (int i = 0; i < n_indices; i ++) + result_data[i] = healpixl_xy_to_ring(indices[i], nside); + } + +done: + free(indices); + return result; +} + + +static PyMethodDef methods[] = { + {"healpix_cone_search", healpix_cone_search, METH_VARARGS | METH_KEYWORDS, NULL}, + {NULL, NULL, 0, NULL} +}; + + +static PyModuleDef moduledef = { + PyModuleDef_HEAD_INIT, + "_core", NULL, -1, methods +}; + +static const PyUFuncGenericFunction + healpix_to_lonlat_loops [] = {healpix_to_lonlat_loop}, + lonlat_to_healpix_loops [] = {lonlat_to_healpix_loop}, + nested_to_ring_loops [] = {nested_to_ring_loop}, + ring_to_nested_loops [] = {ring_to_nested_loop}, + bilinear_interpolation_weights_loops[] = {bilinear_interpolation_weights_loop}, + neighbours_loops [] = {neighbours_loop}; + +static const order_funcs + nested_order_funcs = {healpixl_nested_to_xy, healpixl_xy_to_nested}, + ring_order_funcs = {healpixl_ring_to_xy, healpixl_xy_to_ring}; + +static void *no_ufunc_data[] = {NULL}, + *nested_ufunc_data[] = {&nested_order_funcs}, + *ring_ufunc_data[] = {&ring_order_funcs}; + +static const char + healpix_to_lonlat_types[] = { + NPY_INT64, NPY_INT, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE}, + lonlat_to_healpix_types[] = { + NPY_DOUBLE, NPY_DOUBLE, NPY_INT, NPY_INT64, NPY_DOUBLE, NPY_DOUBLE}, + healpix_to_healpix_types[] = { + NPY_INT64, NPY_INT, NPY_INT64}, + bilinear_interpolation_weights_types[] = { + NPY_DOUBLE, NPY_DOUBLE, NPY_INT, + NPY_INT64, NPY_INT64, NPY_INT64, NPY_INT64, + NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE}, + neighbours_types[] = { + NPY_INT64, NPY_INT, + NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE, + NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE}; + + +PyMODINIT_FUNC PyInit__core(void) +{ + import_array(); + import_umath(); + + PyObject *module = PyModule_Create(&moduledef); + + PyModule_AddObject( + module, "healpix_nested_to_lonlat", PyUFunc_FromFuncAndData( + healpix_to_lonlat_loops, nested_ufunc_data, + healpix_to_lonlat_types, 1, 4, 2, PyUFunc_None, + "healpix_nested_to_lonlat", NULL, 0)); + + PyModule_AddObject( + module, "healpix_ring_to_lonlat", PyUFunc_FromFuncAndData( + healpix_to_lonlat_loops, ring_ufunc_data, + healpix_to_lonlat_types, 1, 4, 2, PyUFunc_None, + "healpix_ring_to_lonlat", NULL, 0)); + + PyModule_AddObject( + module, "lonlat_to_healpix_nested", PyUFunc_FromFuncAndData( + lonlat_to_healpix_loops, nested_ufunc_data, + lonlat_to_healpix_types, 1, 3, 3, PyUFunc_None, + "lonlat_to_healpix_nested", NULL, 0)); + + PyModule_AddObject( + module, "lonlat_to_healpix_ring", PyUFunc_FromFuncAndData( + lonlat_to_healpix_loops, ring_ufunc_data, + lonlat_to_healpix_types, 1, 3, 3, PyUFunc_None, + "lonlat_to_healpix_ring", NULL, 0)); + + PyModule_AddObject( + module, "nested_to_ring", PyUFunc_FromFuncAndData( + nested_to_ring_loops, no_ufunc_data, + healpix_to_healpix_types, 1, 2, 1, PyUFunc_None, + "nested_to_ring", NULL, 0)); + + PyModule_AddObject( + module, "ring_to_nested", PyUFunc_FromFuncAndData( + ring_to_nested_loops, no_ufunc_data, + healpix_to_healpix_types, 1, 2, 1, PyUFunc_None, + "ring_to_nested", NULL, 0)); + + PyModule_AddObject( + module, "bilinear_interpolation_weights", PyUFunc_FromFuncAndData( + bilinear_interpolation_weights_loops, no_ufunc_data, + bilinear_interpolation_weights_types, 1, 3, 8, PyUFunc_None, + "bilinear_interpolation_weights", NULL, 0)); + + PyModule_AddObject( + module, "neighbours_nested", PyUFunc_FromFuncAndData( + neighbours_loops, nested_ufunc_data, + neighbours_types, 1, 2, 8, PyUFunc_None, + "neighbours_nested", NULL, 0)); + + PyModule_AddObject( + module, "neighbours_ring", PyUFunc_FromFuncAndData( + neighbours_loops, ring_ufunc_data, + neighbours_types, 1, 2, 8, PyUFunc_None, + "neighbours_ring", NULL, 0)); + + return module; +} diff --git a/astropy_healpix/core.py b/astropy_healpix/core.py index 9cd0a42..ec0537f 100644 --- a/astropy_healpix/core.py +++ b/astropy_healpix/core.py @@ -9,8 +9,7 @@ from astropy import units as u from astropy.coordinates import Longitude, Latitude -from . import core_cython -from .core_cython import _validate_order +from . import _core __all__ = [ 'nside_to_pixel_area', @@ -44,6 +43,18 @@ def _restore_shape(*args, **kwargs): return np.asscalar(args[0]) +def _validate_order(order): + # We also support upper-case, to support directly the values + # ORDERING = {'RING', 'NESTED'} in FITS headers + # This is currently undocumented in the docstrings. + if order == 'nested' or order == 'NESTED': + return 'nested' + elif order == 'ring' or order == 'RING': + return 'ring' + else: + raise ValueError("order must be 'nested' or 'ring'") + + def _validate_healpix_index(label, healpix_index, nside): npix = nside_to_npix(nside) if np.any((healpix_index < 0) | (healpix_index > npix - 1)): @@ -341,7 +352,7 @@ def healpix_to_lonlat(healpix_index, nside, dx=None, dy=None, order='ring'): ---------- healpix_index : int or `~numpy.ndarray` HEALPix indices (as a scalar or array) - nside : int + nside : int or `~numpy.ndarray` Number of pixels along the side of each of the 12 top-level HEALPix tiles dx, dy : float or `~numpy.ndarray`, optional Offsets inside the HEALPix pixel, which must be in the range [0:1], @@ -357,42 +368,24 @@ def healpix_to_lonlat(healpix_index, nside, dx=None, dy=None, order='ring'): The latitude values """ - if (dx is None and dy is not None) or (dx is not None and dy is None): - raise ValueError('Either both or neither dx and dy must be specified') + if _validate_order(order) == 'ring': + func = _core.healpix_ring_to_lonlat + else: # _validate_order(order) == 'nested' + func = _core.healpix_nested_to_lonlat - healpix_index = np.asarray(healpix_index, dtype=np.int64) - - if dx is None and dy is not None: + if dx is None: dx = 0.5 - elif dx is not None and dy is None: + if dy is None: dy = 0.5 - if dx is not None: - dx = np.asarray(dx, dtype=np.float) - dy = np.asarray(dy, dtype=np.float) - _validate_offset('x', dx) - _validate_offset('y', dy) - healpix_index, dx, dy = np.broadcast_arrays(healpix_index, dx, dy) - dx = dx.ravel() - dy = dy.ravel() + nside = np.asarray(nside, dtype=np.intc) - shape = healpix_index.shape - healpix_index = healpix_index.ravel() - nside = int(nside) - - _validate_healpix_index('healpix_index', healpix_index, nside) - _validate_nside(nside) - order = _validate_order(order) - - if dx is None: - lon, lat = core_cython.healpix_to_lonlat(healpix_index, nside, order) - else: - lon, lat = core_cython.healpix_with_offset_to_lonlat(healpix_index, dx, dy, nside, order) + lon, lat = func(healpix_index, nside, dx, dy) lon = Longitude(lon, unit=u.rad, copy=False) lat = Latitude(lat, unit=u.rad, copy=False) - return _restore_shape(lon, lat, shape=shape) + return lon, lat def lonlat_to_healpix(lon, lat, nside, return_offsets=False, order='ring'): @@ -404,7 +397,7 @@ def lonlat_to_healpix(lon, lat, nside, return_offsets=False, order='ring'): lon, lat : :class:`~astropy.units.Quantity` The longitude and latitude values as :class:`~astropy.units.Quantity` instances with angle units. - nside : int + nside : int or `~numpy.ndarray` Number of pixels along the side of each of the 12 top-level HEALPix tiles order : { 'nested' | 'ring' } Order of HEALPix pixels @@ -422,26 +415,22 @@ def lonlat_to_healpix(lon, lat, nside, return_offsets=False, order='ring'): center of the HEALPix pixels """ - lon = lon.to(u.rad).value - lat = lat.to(u.rad).value + if _validate_order(order) == 'ring': + func = _core.lonlat_to_healpix_ring + else: # _validate_order(order) == 'nested' + func = _core.lonlat_to_healpix_nested - lon, lat = np.broadcast_arrays(lon, lat) + nside = np.asarray(nside, dtype=np.intc) - shape = np.shape(lon) + lon = lon.to_value(u.rad) + lat = lat.to_value(u.rad) - lon = lon.astype(float).ravel() - lat = lat.astype(float).ravel() - - nside = int(nside) - _validate_nside(nside) - order = _validate_order(order) + healpix_index, dx, dy = func(lon, lat, nside) if return_offsets: - healpix_index, dx, dy = core_cython.lonlat_to_healpix_with_offset(lon, lat, nside, order) - return _restore_shape(healpix_index, dx, dy, shape=shape) + return healpix_index, dx, dy else: - healpix_index = core_cython.lonlat_to_healpix(lon, lat, nside, order) - return _restore_shape(healpix_index, shape=shape) + return healpix_index def nested_to_ring(nested_index, nside): @@ -452,7 +441,7 @@ def nested_to_ring(nested_index, nside): ---------- nested_index : int or `~numpy.ndarray` Healpix index using the 'nested' ordering - nside : int + nside : int or `~numpy.ndarray` Number of pixels along the side of each of the 12 top-level HEALPix tiles Returns @@ -461,16 +450,9 @@ def nested_to_ring(nested_index, nside): Healpix index using the 'ring' ordering """ - nested_index = np.asarray(nested_index, dtype=np.int64) - shape = nested_index.shape - nested_index = nested_index.ravel() - nside = int(nside) - - _validate_healpix_index('nested_index', nested_index, nside) - _validate_nside(nside) + nside = np.asarray(nside, dtype=np.intc) - ring_index = core_cython.nested_to_ring(nested_index, nside) - return _restore_shape(ring_index, shape=shape) + return _core.nested_to_ring(nested_index, nside) def ring_to_nested(ring_index, nside): @@ -481,7 +463,7 @@ def ring_to_nested(ring_index, nside): ---------- ring_index : int or `~numpy.ndarray` Healpix index using the 'ring' ordering - nside : int + nside : int or `~numpy.ndarray` Number of pixels along the side of each of the 12 top-level HEALPix tiles Returns @@ -490,16 +472,9 @@ def ring_to_nested(ring_index, nside): Healpix index using the 'nested' ordering """ - ring_index = np.asarray(ring_index, dtype=np.int64) - shape = ring_index.shape - ring_index = ring_index.ravel() - nside = int(nside) - - _validate_healpix_index('ring_index', ring_index, nside) - _validate_nside(nside) + nside = np.asarray(nside, dtype=np.intc) - nested_index = core_cython.ring_to_nested(ring_index, nside) - return _restore_shape(nested_index, shape=shape) + return _core.ring_to_nested(ring_index, nside) def bilinear_interpolation_weights(lon, lat, nside, order='ring'): @@ -527,23 +502,17 @@ def bilinear_interpolation_weights(lon, lat, nside, order='ring'): interpolation """ - lon = lon.to(u.rad).value - lat = lat.to(u.rad).value - - lon, lat = np.broadcast_arrays(lon, lat) + lon = lon.to_value(u.rad) + lat = lat.to_value(u.rad) - shape = np.shape(lon) + nside = np.asarray(nside, dtype=np.intc) - lon = lon.astype(float).ravel() - lat = lat.astype(float).ravel() - nside = int(nside) + result = _core.bilinear_interpolation_weights(lon, lat, nside) + indices = np.stack(result[:4]) + weights = np.stack(result[4:]) - order = _validate_order(order) - _validate_nside(nside) - - indices, weights = core_cython.bilinear_interpolation_weights(lon, lat, nside, order) - indices = _restore_shape(indices, shape=(4,) + shape) - weights = _restore_shape(weights, shape=(4,) + shape) + if _validate_order(order) == 'nested': + indices = ring_to_nested(indices, nside) return indices, weights @@ -602,19 +571,14 @@ def neighbours(healpix_index, nside, order='ring'): ``neigh`` has shape (8, 2, 3). """ - healpix_index = np.asarray(healpix_index, dtype=np.int64) - nside = int(nside) - - shape = np.shape(healpix_index) - - healpix_index = healpix_index.ravel() + nside = np.asarray(nside, dtype=np.intc) - _validate_healpix_index('healpix_index', healpix_index, nside) - _validate_nside(nside) - order = _validate_order(order) + if _validate_order(order) == 'ring': + func = _core.neighbours_ring + else: # _validate_order(order) == 'nested' + func = _core.neighbours_nested - neigh = core_cython.neighbours(healpix_index, nside, order) - return _restore_shape(neigh, shape=(8,) + shape) + return np.stack(func(healpix_index, nside)) def healpix_cone_search(lon, lat, radius, nside, order='ring'): @@ -643,15 +607,14 @@ def healpix_cone_search(lon, lat, radius, nside, order='ring'): 1-D array with all the matching HEALPix pixel indices. """ - lon = float(lon.to(u.deg).value) - lat = float(lat.to(u.deg).value) - radius = float(radius.to(u.deg).value) - nside = int(nside) + lon = lon.to_value(u.deg) + lat = lat.to_value(u.deg) + radius = radius.to_value(u.deg) _validate_nside(nside) order = _validate_order(order) - return core_cython.healpix_cone_search(lon, lat, radius, nside, order) + return _core.healpix_cone_search(lon, lat, radius, nside, order) def boundaries_lonlat(healpix_index, step, nside, order='ring'): diff --git a/astropy_healpix/core_cython.pxd b/astropy_healpix/core_cython.pxd deleted file mode 100644 index bea16e5..0000000 --- a/astropy_healpix/core_cython.pxd +++ /dev/null @@ -1,44 +0,0 @@ -# Licensed under a 3-clause BSD style license - see LICENSE.rst - -# This file is needed in order to be able to cimport functions into other Cython -# files - -from numpy cimport int64_t - -cdef extern from "healpix.h": - - # Converts a healpix index from the XY scheme to the RING scheme. - int64_t healpixl_xy_to_ring(int64_t hp, int Nside) nogil - - # Converts a healpix index from the RING scheme to the XY scheme. - int64_t healpixl_ring_to_xy(int64_t ring_index, int Nside) nogil - - # Converts a healpix index from the XY scheme to the NESTED scheme. - int64_t healpixl_xy_to_nested(int64_t hp, int Nside) nogil - - # Converts a healpix index from the NESTED scheme to the XY scheme. - int64_t healpixl_nested_to_xy(int64_t nested_index, int Nside) nogil - - # Converts a healpix index, plus fractional offsets (dx,dy), into (x,y,z) - # coordinates on the unit sphere. (dx,dy) must be in [0, 1]. (0.5, 0.5) - # is the center of the healpix. (0,0) is the southernmost corner, (1,1) is - # the northernmost corner, (1,0) is the easternmost, and (0,1) the - # westernmost. - void healpixl_to_radec(int64_t hp, int Nside, double dx, double dy, - double* ra, double* dec) nogil - - # Converts (RA, DEC) coordinates (in radians) to healpix XY index - int64_t radec_to_healpixlf(double ra, double dec, int Nside, double* dx, double* dy) nogil - - # Finds the healpixes neighbouring the given healpix, placing them in the - # array "neighbour". Returns the number of neighbours. You must ensure - # that "neighbour" has 8 elements. - # - # Healpixes in the interior of a large healpix will have eight neighbours; - # pixels near the edges can have fewer. - int healpixl_get_neighbours(int64_t hp, int64_t* neighbours, int Nside) nogil - - int healpix_rangesearch_radec_simple(double ra, double dec, double radius, int Nside, int approx, int64_t** indices) nogil - - void interpolate_weights(double lon, double lat, int64_t *ring_indices, - double *weights, int Nside) nogil diff --git a/astropy_healpix/core_cython.pyx b/astropy_healpix/core_cython.pyx deleted file mode 100644 index 5c77b5f..0000000 --- a/astropy_healpix/core_cython.pyx +++ /dev/null @@ -1,492 +0,0 @@ -# Licensed under a 3-clause BSD style license - see LICENSE.rst -""" -This module contains vectorized Cython functions for common HEALPix operations. -Since they are written in Cython rather than Python, their input types are -strict and the functions will fail if the incorrect types are passed in. -""" - -import numpy as np -cimport numpy as np -import cython -from libc.stdlib cimport abort, malloc, free -from libc.math cimport sin, cos, sqrt - -ctypedef np.intp_t intp_t -ctypedef np.double_t double_t - -npy_double = np.double -npy_intp = np.intp -npy_int64 = np.int64 - - -def _validate_order(str order): - # We also support upper-case, to support directly the values - # ORDERING = {'RING', 'NESTED'} in FITS headers - # This is currently undocumented in the docstrings. - if order == 'nested' or order == 'NESTED': - return 'nested' - elif order == 'ring' or order == 'RING': - return 'ring' - else: - raise ValueError("order must be 'nested' or 'ring'") - - -@cython.boundscheck(False) -def healpix_to_lonlat(np.ndarray[int64_t, ndim=1, mode="c"] healpix_index, - int nside, str order): - """ - Convert HEALPix indices to longitudes/latitudes. - - This returns the longitudes/latitudes of the center of the HEALPix pixels. - If you also want to provide relative offsets inside the pixels, see - :func:`healpix_with_offset_to_lonlat`. - - Parameters - ---------- - healpix_index : `~numpy.ndarray` - 1-D array of HEALPix indices - nside : int - Number of pixels along the side of each of the 12 top-level HEALPix tiles - order : { 'nested' | 'ring' } - Order of HEALPix pixels - - Returns - ------- - lon, lat : `~numpy.ndarray` - 1-D arrays of longitude and latitude in radians - """ - - cdef intp_t n = healpix_index.shape[0] - cdef intp_t i - cdef int64_t xy_index - cdef double dx, dy; - cdef np.ndarray[double_t, ndim=1, mode="c"] lon = np.zeros(n, dtype=npy_double) - cdef np.ndarray[double_t, ndim=1, mode="c"] lat = np.zeros(n, dtype=npy_double) - - dx = 0.5 - dy = 0.5 - - order = _validate_order(order) - - if order == 'nested': - for i in range(n): - xy_index = healpixl_nested_to_xy(healpix_index[i], nside) - healpixl_to_radec(xy_index, nside, dx, dy, &lon[i], &lat[i]) - elif order == 'ring': - for i in range(n): - xy_index = healpixl_ring_to_xy(healpix_index[i], nside) - healpixl_to_radec(xy_index, nside, dx, dy, &lon[i], &lat[i]) - - return lon, lat - - -@cython.boundscheck(False) -def healpix_with_offset_to_lonlat(np.ndarray[int64_t, ndim=1, mode="c"] healpix_index, - np.ndarray[double_t, ndim=1, mode="c"] dx, - np.ndarray[double_t, ndim=1, mode="c"] dy, - int nside, str order): - """ - Convert HEALPix indices to longitudes/latitudes - - This function takes relative offsets in x and y inside the HEALPix pixels. - If you are only interested in the centers of the pixels, see - `healpix_to_lonlat`. - - Parameters - ---------- - healpix_index : `~numpy.ndarray` - 1-D array of HEALPix indices - dx, dy : `~numpy.ndarray` - 1-D arrays of offsets inside the HEALPix pixel, which must be in the - range [0:1] (0.5 is the center of the HEALPix pixels) - nside : int - Number of pixels along the side of each of the 12 top-level HEALPix tiles - order : { 'nested' | 'ring' } - Order of HEALPix pixels - - - Returns - ------- - lon, lat : `~numpy.ndarray` - 1-D arrays of longitude and latitude in radians - """ - - cdef intp_t n = healpix_index.shape[0] - cdef intp_t i - cdef int64_t xy_index - cdef np.ndarray[double_t, ndim=1, mode="c"] lon = np.zeros(n, dtype=npy_double) - cdef np.ndarray[double_t, ndim=1, mode="c"] lat = np.zeros(n, dtype=npy_double) - - order = _validate_order(order) - - if order == 'nested': - for i in range(n): - xy_index = healpixl_nested_to_xy(healpix_index[i], nside) - healpixl_to_radec(xy_index, nside, dx[i], dy[i], &lon[i], &lat[i]) - elif order == 'ring': - for i in range(n): - xy_index = healpixl_ring_to_xy(healpix_index[i], nside) - healpixl_to_radec(xy_index, nside, dx[i], dy[i], &lon[i], &lat[i]) - - return lon, lat - - -@cython.boundscheck(False) -def lonlat_to_healpix(np.ndarray[double_t, ndim=1, mode="c"] lon, - np.ndarray[double_t, ndim=1, mode="c"] lat, - int nside, str order): - """ - Convert longitudes/latitudes to HEALPix indices - - This returns only the HEALPix indices. If you also want to get relative - offsets inside the pixels, see :func:`lonlat_to_healpix_with_offset`. - - Parameters - ---------- - lon, lat : `~numpy.ndarray` - 1-D arrays of longitude and latitude in radians - nside : int - Number of pixels along the side of each of the 12 top-level HEALPix tiles - order : { 'nested' | 'ring' } - Order of HEALPix pixels - - - Returns - ------- - healpix_index : `~numpy.ndarray` - 1-D array of HEALPix indices - """ - - cdef intp_t n = lon.shape[0] - cdef intp_t i - cdef int64_t xy_index - cdef double dx, dy; - cdef np.ndarray[int64_t, ndim=1, mode="c"] healpix_index = np.zeros(n, dtype=npy_int64) - - order = _validate_order(order) - - if order == 'nested': - for i in range(n): - xy_index = radec_to_healpixlf(lon[i], lat[i], nside, &dx, &dy) - healpix_index[i] = healpixl_xy_to_nested(xy_index, nside) - elif order == 'ring': - for i in range(n): - xy_index = radec_to_healpixlf(lon[i], lat[i], nside, &dx, &dy) - healpix_index[i] = healpixl_xy_to_ring(xy_index, nside) - - return healpix_index - - -@cython.boundscheck(False) -def lonlat_to_healpix_with_offset(np.ndarray[double_t, ndim=1, mode="c"] lon, - np.ndarray[double_t, ndim=1, mode="c"] lat, - int nside, str order): - """ - Convert longitudes/latitudes to healpix indices - - This returns the HEALPix indices and relative offsets inside the pixels. If - you want only the HEALPix indices, see :func:`lonlat_to_healpix`. - - Parameters - ---------- - lon, lat : `~numpy.ndarray` - 1-D arrays of longitude and latitude in radians - nside : int - Number of pixels along the side of each of the 12 top-level HEALPix tiles - order : { 'nested' | 'ring' } - Order of HEALPix pixels - - Returns - ------- - healpix_index : `~numpy.ndarray` - 1-D array of HEALPix indices - dx, dy : `~numpy.ndarray` - 1-D arrays of offsets inside the HEALPix pixel in the range [0:1] (0.5 - is the center of the HEALPix pixels) - """ - - cdef intp_t n = lon.shape[0] - cdef intp_t i - cdef int64_t xy_index - cdef np.ndarray[int64_t, ndim=1, mode="c"] healpix_index = np.zeros(n, dtype=npy_int64) - cdef np.ndarray[double_t, ndim=1, mode="c"] dx = np.zeros(n, dtype=npy_double) - cdef np.ndarray[double_t, ndim=1, mode="c"] dy = np.zeros(n, dtype=npy_double) - - order = _validate_order(order) - - if order == 'nested': - for i in range(n): - xy_index = radec_to_healpixlf(lon[i], lat[i], nside, &dx[i], &dy[i]) - healpix_index[i] = healpixl_xy_to_nested(xy_index, nside) - elif order == 'ring': - for i in range(n): - xy_index = radec_to_healpixlf(lon[i], lat[i], nside, &dx[i], &dy[i]) - healpix_index[i] = healpixl_xy_to_ring(xy_index, nside) - - return healpix_index, dx, dy - - -@cython.boundscheck(False) -def nested_to_ring(np.ndarray[int64_t, ndim=1, mode="c"] nested_index, int nside): - """ - Convert a HEALPix 'nested' index to a HEALPix 'ring' index - - Parameters - ---------- - nested_index : `~numpy.ndarray` - Healpix index using the 'nested' ordering - nside : int - Number of pixels along the side of each of the 12 top-level HEALPix tiles - - Returns - ------- - ring_index : `~numpy.ndarray` - Healpix index using the 'ring' ordering - """ - - cdef intp_t n = nested_index.shape[0] - cdef intp_t i - cdef np.ndarray[int64_t, ndim=1, mode="c"] ring_index = np.zeros(n, dtype=npy_int64) - - for i in range(n): - ring_index[i] = healpixl_xy_to_ring(healpixl_nested_to_xy(nested_index[i], nside), nside) - - return ring_index - - -@cython.boundscheck(False) -def ring_to_nested(np.ndarray[int64_t, ndim=1, mode="c"] ring_index, int nside): - """ - Convert a HEALPix 'ring' index to a HEALPix 'nested' index - - Parameters - ---------- - ring_index : `~numpy.ndarray` - Healpix index using the 'ring' ordering - nside : int - Number of pixels along the side of each of the 12 top-level HEALPix tiles - - Returns - ------- - nested_index : `~numpy.ndarray` - Healpix index using the 'nested' ordering - """ - - cdef intp_t n = ring_index.shape[0] - cdef intp_t i - cdef np.ndarray[int64_t, ndim=1, mode="c"] nested_index = np.zeros(n, dtype=npy_int64) - - for i in range(n): - nested_index[i] = healpixl_xy_to_nested(healpixl_ring_to_xy(ring_index[i], nside), nside) - - return nested_index - - - -@cython.boundscheck(False) -def bilinear_interpolation_weights(np.ndarray[double_t, ndim=1, mode="c"] lon, - np.ndarray[double_t, ndim=1, mode="c"] lat, - int nside, str order): - """ - Get the four neighbours for each (lon, lat) position and the weight - associated with each one for bilinear interpolation. - - Parameters - ---------- - lon, lat : `~numpy.ndarray` - 1-D arrays of longitude and latitude in radians - nside : int - Number of pixels along the side of each of the 12 top-level HEALPix tiles - order : { 'nested' | 'ring' } - Order of HEALPix pixels - - Returns - ------- - indices : `~numpy.ndarray` - 2-D array with shape (4, N) giving the four indices to use for the - interpolation - weights : `~numpy.ndarray` - 2-D array with shape (4, N) giving the four weights to use for the - interpolation - """ - - cdef intp_t n = lon.shape[0] - cdef intp_t i, j - cdef np.ndarray[int64_t, ndim=2, mode="c"] indices = np.zeros((4, n), dtype=npy_int64) - cdef np.ndarray[double_t, ndim=2, mode="c"] weights = np.zeros((4, n), dtype=npy_double) - cdef int order_int - - cdef double *weights_indiv - cdef int64_t *indices_indiv - - order = _validate_order(order) - - if order == 'nested': - order_int = 0 - elif order == 'ring': - order_int = 1 - - with nogil: - - indices_indiv = malloc(sizeof(int64_t) * 4) - if indices_indiv == NULL: - abort() - - weights_indiv = malloc(sizeof(double) * 4) - if weights_indiv == NULL: - abort() - - for i in range(n): - interpolate_weights(lon[i], lat[i], indices_indiv, weights_indiv, nside) - for j in range(4): - if order_int == 0: - indices[j, i] = healpixl_xy_to_nested(healpixl_ring_to_xy(indices_indiv[j], nside), nside) - else: - indices[j, i] = indices_indiv[j] - weights[j, i] = weights_indiv[j] - - return indices, weights - - -@cython.boundscheck(False) -def neighbours(np.ndarray[int64_t, ndim=1, mode="c"] healpix_index, - int nside, str order): - """ - Find all the HEALPix pixels that are the neighbours of a HEALPix pixel - - Parameters - ---------- - healpix_index : `~numpy.ndarray` - 1-D array of HEALPix indices - nside : int - Number of pixels along the side of each of the 12 top-level HEALPix tiles - order : { 'nested' | 'ring' } - Order of HEALPix pixels - - Returns - ------- - neighbours : `~numpy.ndarray` - 2-D array with shape (8, N) giving the neighbours starting SW and - rotating clockwise. - """ - - cdef intp_t n = healpix_index.shape[0] - cdef intp_t i - cdef int64_t xy_index - cdef int j, k - cdef np.ndarray[int64_t, ndim=2, mode="c"] neighbours = np.zeros((8, n), dtype=npy_int64) - cdef int64_t * neighbours_indiv - cdef int order_int - - order = _validate_order(order) - - if order == 'nested': - order_int = 0 - elif order == 'ring': - order_int = 1 - - with nogil: - - neighbours_indiv = malloc(sizeof(int64_t) * 8) - if neighbours_indiv == NULL: - abort() - - # The neighbours above are ordered as follows: - # - # 3 2 1 - # 4 X 0 - # 5 6 7 - # - # but we want: - # - # 2 3 4 - # 1 X 5 - # 0 7 6 - # - # so we reorder these on-the-fly - - if order_int == 0: - - for i in range(n): - - xy_index = healpixl_nested_to_xy(healpix_index[i], nside) - healpixl_get_neighbours(xy_index, neighbours_indiv, nside) - - for j in range(8): - k = 4 - j - if k < 0: - k = k + 8 - if neighbours_indiv[k] < 0: - neighbours[j, i] = -1 - else: - neighbours[j, i] = healpixl_xy_to_nested(neighbours_indiv[k], nside) - - elif order_int == 1: - - for i in range(n): - - xy_index = healpixl_ring_to_xy(healpix_index[i], nside) - - healpixl_get_neighbours(xy_index, neighbours_indiv, nside) - - for j in range(8): - k = 4 - j - if k < 0: - k = k + 8 - if neighbours_indiv[k] < 0: - neighbours[j, i] = -1 - else: - neighbours[j, i] = healpixl_xy_to_ring(neighbours_indiv[k], nside) - - free(neighbours_indiv) - - return neighbours - - -@cython.boundscheck(False) -def healpix_cone_search(double lon, double lat, double radius, int nside, str order): - """ - Find all the HEALPix pixels within a given radius of a longitude/latitude. - - Note that this returns all pixels that overlap, including partially, with - the search cone. This function can only be used for a single lon/lat pair at - a time, since different calls to the function may result in a different - number of matches. - - Parameters - ---------- - lon, lat : float - The longitude and latitude to search around, in degrees - radius : float - The search radius, in degrees - nside : int - Number of pixels along the side of each of the 12 top-level HEALPix tiles - order : { 'nested' | 'ring' } - Order of HEALPix pixels - - Returns - ------- - healpix_index : `~numpy.ndarray` - 1-D array with all the matching HEALPix pixel indices. - """ - cdef intp_t i - cdef int64_t *indices - cdef int64_t n_indices - cdef int64_t index - - n_indices = healpix_rangesearch_radec_simple(lon, lat, radius, nside, 0, &indices) - - cdef np.ndarray[int64_t, ndim=1, mode="c"] result = np.zeros(n_indices, dtype=npy_int64) - - order = _validate_order(order) - - if order == 'nested': - for i in range(n_indices): - index = indices[i] - result[i] = healpixl_xy_to_nested(index, nside) - elif order == 'ring': - for i in range(n_indices): - index = indices[i] - result[i] = healpixl_xy_to_ring(index, nside) - - return result diff --git a/astropy_healpix/setup_package.py b/astropy_healpix/setup_package.py index cb73fb3..eb66d89 100644 --- a/astropy_healpix/setup_package.py +++ b/astropy_healpix/setup_package.py @@ -24,12 +24,12 @@ def get_extensions(): sources = [os.path.join(C_DIR, filename) for filename in C_FILES] sources.append(os.path.join(HEALPIX_ROOT, 'interpolation.c')) - sources.append(os.path.join(HEALPIX_ROOT, 'core_cython.pyx')) + sources.append(os.path.join(HEALPIX_ROOT, '_core.c')) include_dirs = ['numpy', C_DIR, HEALPIX_ROOT] extension = Extension( - name="astropy_healpix.core_cython", + name="astropy_healpix._core", sources=sources, include_dirs=include_dirs, libraries=libraries, diff --git a/astropy_healpix/tests/test_core_cython.py b/astropy_healpix/tests/test_core_cython.py deleted file mode 100644 index 9977794..0000000 --- a/astropy_healpix/tests/test_core_cython.py +++ /dev/null @@ -1,129 +0,0 @@ -# Licensed under a 3-clause BSD style license - see LICENSE.rst - -from __future__ import absolute_import, print_function, division - -from itertools import product - -import pytest - -import numpy as np -from numpy.testing import assert_equal, assert_allclose - -from astropy import units as u -from astropy.coordinates.angle_utilities import angular_separation - -from ..core import nside_to_pixel_resolution -from .. import core_cython - -NSIDE_POWERS = range(0, 21) -ORDERS = ('nested', 'ring') - - -def get_test_indices(nside): - # For large number of pixels, we only compute a random subset of points - if nside > 2 ** 8: - try: - return np.random.randint(0, 12 * nside ** 2, 12 * 8 ** 2, dtype=np.int64) - except TypeError: # Numpy 1.9 and 1.10 - return (np.random.random(12 * 8 ** 2) * (12 * float(nside) ** 2)).astype(np.int64, copy=False) - else: - return np.arange(12 * nside ** 2, dtype=np.int64) - - -# NOTE: we use capfd in all tests here to make sure no errors/warnings are being -# raised by the C code. - -@pytest.mark.parametrize(('order', 'nside_power'), product(ORDERS, NSIDE_POWERS)) -def test_roundtrip_healpix_no_offsets(order, nside_power, capfd): - nside = 2 ** nside_power - index = get_test_indices(nside) - lon, lat = core_cython.healpix_to_lonlat(index, nside, order) - index_new = core_cython.lonlat_to_healpix(lon, lat, nside, order) - assert_equal(index, index_new) - out, err = capfd.readouterr() - assert out == "" and err == "" - - -@pytest.mark.parametrize(('order', 'nside_power'), product(ORDERS, NSIDE_POWERS)) -def test_roundtrip_healpix_with_offsets(order, nside_power, capfd): - nside = 2 ** nside_power - index = get_test_indices(nside) - dx = np.random.random(index.shape) - dy = np.random.random(index.shape) - lon, lat = core_cython.healpix_with_offset_to_lonlat(index, dx, dy, nside, order) - index_new, dx_new, dy_new = core_cython.lonlat_to_healpix_with_offset(lon, lat, nside, order) - assert_equal(index, index_new) - assert_allclose(dx, dx_new, atol=1e-8) - assert_allclose(dy, dy_new, atol=1e-8) - out, err = capfd.readouterr() - assert out == "" and err == "" - - -@pytest.mark.parametrize('nside_power', NSIDE_POWERS) -def test_roundtrip_nested_ring(nside_power, capfd): - nside = 2 ** nside_power - nested_index = get_test_indices(nside) - ring_index = core_cython.nested_to_ring(nested_index, nside) - nested_index_new = core_cython.ring_to_nested(ring_index, nside) - assert_equal(nested_index, nested_index_new) - if nside == 1: - assert np.all(nested_index == ring_index) - else: - assert not np.all(nested_index == ring_index) - out, err = capfd.readouterr() - assert out == "" and err == "" - - -@pytest.mark.parametrize(('order', 'nside_power'), product(ORDERS, NSIDE_POWERS)) -def test_neighbours(order, nside_power, capfd): - # This just makes sure things run, but doesn't check the validity of result - nside = 2 ** nside_power - index = get_test_indices(nside) - neighbours = core_cython.neighbours(index, nside, order) - assert np.all(neighbours >= -1) and np.all(neighbours < 12 * nside ** 2) - out, err = capfd.readouterr() - assert out == "" and err == "" - - -@pytest.mark.parametrize(('order', 'nside_power'), product(ORDERS, NSIDE_POWERS)) -def test_healpix_cone_search(order, nside_power, capfd): - nside = 2 ** nside_power - lon0, lat0 = 12., 40. - radius = nside_to_pixel_resolution(nside).to(u.degree).value * 10 - index_inside = core_cython.healpix_cone_search(lon0, lat0, radius, nside, order) - n_inside = len(index_inside) - dx = np.array([[0.0, 0.0, 1.0, 1.0]]) - dy = np.array([[0.0, 1.0, 1.0, 0.0]]) - dx = np.repeat(dx, n_inside, axis=0).ravel() - dy = np.repeat(dy, n_inside, axis=0).ravel() - index_inside = np.repeat(index_inside, 4).ravel() - lon, lat = core_cython.healpix_with_offset_to_lonlat(index_inside, dx, dy, nside, order) - lon, lat = lon.reshape((n_inside, 4)), lat.reshape((n_inside, 4)) - sep = angular_separation(lon0 * u.deg, lat0 * u.deg, lon * u.rad, lat * u.rad) - sep = sep.min(axis=1) - assert np.all(sep.to(u.degree).value < radius * 1.05) - out, err = capfd.readouterr() - assert out == "" and err == "" - - -# Regression tests for fixed bugs - -def test_regression_healpix_to_lonlat_sqrt(): - - # Regression test for a bug that caused the ring index decomposition to fail - # and return a negative longitude. - - index = np.array([9007199120523263], dtype=np.int64) - lon, lat = core_cython.healpix_to_lonlat(index, 2**26, order='ring') - assert_allclose(lon, 6.283185295476241, rtol=1e-14) - assert_allclose(lat, 0.729727669554970, rtol=1e-14) - - index = np.array([720575940916150240], dtype=np.int64) - lon, lat = core_cython.healpix_to_lonlat(index, 2**28, order='ring') - assert_allclose(lon, 6.283185122851909, rtol=1e-14) - assert_allclose(lat, -0.729727656226966, rtol=1e-14) - - index = np.array([180143985363255292], dtype=np.int64) - lon, lat = core_cython.healpix_to_lonlat(index, 2**27, order='ring') - assert_allclose(lon, 6.283185266217880, rtol=1e-14) - assert_allclose(lat, -0.729727656226966, rtol=1e-14) diff --git a/cextern/astrometry.net/healpix-utils.h b/cextern/astrometry.net/healpix-utils.h index a59aa71..f1ad43e 100644 --- a/cextern/astrometry.net/healpix-utils.h +++ b/cextern/astrometry.net/healpix-utils.h @@ -15,7 +15,7 @@ il* healpix_rangesearch_xyz(const double* xyz, double radius, int Nside, il* hps il* healpix_rangesearch_xyz_approx(const double* xyz, double radius, int Nside, il* hps); il* healpix_rangesearch_radec_approx(double ra, double dec, double radius, int Nside, il* hps); il* healpix_rangesearch_radec(double ra, double dec, double radius, int Nside, il* hps); -int healpix_rangesearch_radec_simple(double ra, double dec, double radius, int Nside, int **indices); +int64_t healpix_rangesearch_radec_simple(double ra, double dec, double radius, int Nside, int approx, int64_t **indices); /** Starting from a "seed" or list of "seeds" healpixes, grows a region