diff --git a/.gitignore b/.gitignore index 2ed2b75ea..8739f0df0 100644 --- a/.gitignore +++ b/.gitignore @@ -2,6 +2,12 @@ .pydevproject .settings build +dist *.pyc *.so __pycache__ +reproject/cython_version.py +reproject/reproject.cfg +reproject/version.py +reproject.egg-info +reproject/_overlap.c diff --git a/.travis.yml b/.travis.yml index c397bfabf..e97daf623 100644 --- a/.travis.yml +++ b/.travis.yml @@ -27,8 +27,9 @@ before_install: install: - - $PIP_WHEEL_COMMAND "numpy==1.8.0" - - $PIP_WHEEL_COMMAND "astropy==0.3" + - $PIP_WHEEL_COMMAND "numpy==1.8.0" + - $PIP_WHEEL_COMMAND "astropy==0.3" + - $PIP_WHEEL_COMMAND "Cython==0.19" script: - if [[ $SETUP_CMD == 'test' ]]; then python setup.py test; fi diff --git a/reproject/__init__.py b/reproject/__init__.py index 6e2649899..debf43233 100644 --- a/reproject/__init__.py +++ b/reproject/__init__.py @@ -114,4 +114,4 @@ def test(package=None, test_path=None, args=None, plugins=None, remote_data=remote_data, pep8=pep8, pdb=pdb, coverage=coverage, open_files=open_files, **kwargs) -from .overlap import * +from .overlap import compute_overlap diff --git a/reproject/_overlap.pyx b/reproject/_overlap.pyx new file mode 100644 index 000000000..3cf5c2245 --- /dev/null +++ b/reproject/_overlap.pyx @@ -0,0 +1,28 @@ +import numpy as np +cimport numpy as np +import cython + +ctypedef np.double_t DOUBLE_T + +cdef extern from "overlapArea.h": + double computeOverlap(double * ilon, double * ilat, double * olon, double * olat, + int energyMode, double refArea, double * areaRatio) + + +# @cython.wraparound(False) +# @cython.boundscheck(False) +def _compute_overlap(np.ndarray[double, ndim=2] ilon, + np.ndarray[double, ndim=2] ilat, + np.ndarray[double, ndim=2] olon, + np.ndarray[double, ndim=2] olat): + cdef int i + cdef int n = ilon.shape[0] + + cdef np.ndarray[double, ndim = 1] overlap = np.empty(n, dtype=np.double) + cdef np.ndarray[double, ndim = 1] area_ratio = np.empty(n, dtype=np.double) + + for i in range(n): + overlap[i] = computeOverlap(& ilon[i, 0], & ilat[i, 0], & olon[i, 0], & olat[i, 0], + 0, 1, & area_ratio[i]) + + return overlap, area_ratio diff --git a/reproject/_overlap_wrapper.c b/reproject/_overlap_wrapper.c deleted file mode 100644 index 73c39af31..000000000 --- a/reproject/_overlap_wrapper.c +++ /dev/null @@ -1,97 +0,0 @@ -#include -#include -#include - -/* Define docstrings */ -static char module_docstring[] = "Wrap Montage overlap routines"; -static char computeOverlap_docstring[] = "Compute spherical polygon overlap"; - -/* Declare the C functions here. */ -static PyObject *_computeOverlap(PyObject *self, PyObject *args); - -/* Define the methods that will be available on the module. */ -static PyMethodDef module_methods[] = { - {"_computeOverlap", _computeOverlap, METH_VARARGS, computeOverlap_docstring}, - {NULL, NULL, 0, NULL} -}; - -/* This is the function that is called on import. */ - -#if PY_MAJOR_VERSION >= 3 - #define MOD_ERROR_VAL NULL - #define MOD_SUCCESS_VAL(val) val - #define MOD_INIT(name) PyMODINIT_FUNC PyInit_##name(void) - #define MOD_DEF(ob, name, doc, methods) \ - static struct PyModuleDef moduledef = { \ - PyModuleDef_HEAD_INIT, name, doc, -1, methods, }; \ - ob = PyModule_Create(&moduledef); -#else - #define MOD_ERROR_VAL - #define MOD_SUCCESS_VAL(val) - #define MOD_INIT(name) void init##name(void) - #define MOD_DEF(ob, name, doc, methods) \ - ob = Py_InitModule3(name, methods, doc); -#endif - -MOD_INIT(_overlap_wrapper) -{ - PyObject *m; - MOD_DEF(m, "_overlap_wrapper", module_docstring, module_methods); - if (m == NULL) - return MOD_ERROR_VAL; - import_array(); - return MOD_SUCCESS_VAL(m); -} - -/* Define interfaces */ - -double computeOverlap(double *, double *, double *, double *, int, double, double *); - -/* Do the heavy lifting here */ - -static PyObject *_computeOverlap(PyObject *self, PyObject *args) -{ - - PyObject *ilon_obj, *ilat_obj, *olon_obj, *olat_obj; - double refArea, areaRatio; - int energyMode; - - /* Parse the input tuple */ - if (!PyArg_ParseTuple(args, "OOOOid", &ilon_obj, &ilat_obj, &olon_obj, &olat_obj, &energyMode, &refArea)) - return NULL; - - /* Interpret the input objects as `numpy` arrays. */ - PyObject *ilon_array = PyArray_FROM_OTF(ilon_obj, NPY_DOUBLE, NPY_IN_ARRAY); - PyObject *ilat_array = PyArray_FROM_OTF(ilat_obj, NPY_DOUBLE, NPY_IN_ARRAY); - PyObject *olon_array = PyArray_FROM_OTF(olon_obj, NPY_DOUBLE, NPY_IN_ARRAY); - PyObject *olat_array = PyArray_FROM_OTF(olat_obj, NPY_DOUBLE, NPY_IN_ARRAY); - - /* If that didn't work, throw an `Exception`. */ - if (ilon_array == NULL || ilat_array == NULL || olon_array == NULL || olat_array == NULL) { - PyErr_SetString(PyExc_TypeError, "Couldn't parse the input arrays."); - Py_XDECREF(ilon_array); - Py_XDECREF(ilat_array); - Py_XDECREF(olon_array); - Py_XDECREF(olat_array); - return NULL; - } - - - /* Get pointers to the data as C-types. */ - double *ilon = (double*)PyArray_DATA(ilon_array); - double *ilat = (double*)PyArray_DATA(ilat_array); - double *olon = (double*)PyArray_DATA(olon_array); - double *olat = (double*)PyArray_DATA(olat_array); - - /* Compute overlap using Montage routine */ - double overlap = computeOverlap(ilon, ilat, olon, olat, energyMode, refArea, &areaRatio); - - /* Build the output tuple */ - PyObject *ret = Py_BuildValue("dd", overlap, areaRatio); - if (ret == NULL) { - PyErr_SetString(PyExc_RuntimeError, "Couldn't build output."); - return NULL; - } - - return ret; -} \ No newline at end of file diff --git a/reproject/overlap.py b/reproject/overlap.py index bd169d3f0..84e27517c 100644 --- a/reproject/overlap.py +++ b/reproject/overlap.py @@ -1,4 +1,5 @@ -from ._overlap_wrapper import _computeOverlap +import numpy as np +from ._overlap import _compute_overlap __all__ = ['compute_overlap'] @@ -8,20 +9,25 @@ def compute_overlap(ilon, ilat, olon, olat): Parameters ---------- - ilon : np.ndarray + ilon : np.ndarray with shape (N, 4) The longitudes (in radians) defining the four corners of the input pixel - ilat : np.ndarray + ilat : np.ndarray with shape (N, 4) The latitudes (in radians) defining the four corners of the input pixel - olon : np.ndarray + olon : np.ndarray with shape (N, 4) The longitudes (in radians) defining the four corners of the output pixel - olat : np.ndarray + olat : np.ndarray with shape (N, 4) The latitudes (in radians) defining the four corners of the output pixel Returns ------- - overlap : np.ndarray + overlap : np.ndarray of length N Pixel overlap solid angle in steradians - area_ratio : np.ndarray + area_ratio : np.ndarray of length N TODO """ - return _computeOverlap(ilon, ilat, olon, olat, 0, 1.) + ilon = np.asarray(ilon, dtype=np.float64) + ilat = np.asarray(ilat, dtype=np.float64) + olon = np.asarray(olon, dtype=np.float64) + olat = np.asarray(olat, dtype=np.float64) + + return _compute_overlap(ilon, ilat, olon, olat) diff --git a/reproject/overlapArea.c b/reproject/overlapArea.c index 1c3aa9ff2..56c798778 100644 --- a/reproject/overlapArea.c +++ b/reproject/overlapArea.c @@ -6,6 +6,7 @@ #include #include #include "mNaN.h" +#include "overlapArea.h" // Constants @@ -41,9 +42,6 @@ typedef struct vec { // Function prototypes -double computeOverlap(double *ilon, double *ilat, double *olon, double *olat, - int energyMode, double refArea, double *areaRatio); - int DirectionCalculator(Vec *a, Vec *b, Vec *c); int SegSegIntersect(Vec *a, Vec *b, Vec *c, Vec *d, Vec *e, Vec *f, Vec *p, Vec *q); @@ -819,8 +817,6 @@ double Dot(Vec *a, Vec *b) { double Normalize(Vec *v) { double len; - len = 0.; - len = sqrt(v->x * v->x + v->y * v->y + v->z * v->z); if (len == 0.) @@ -858,10 +854,10 @@ double Girard() { double sumang, cosAng, sinAng; - sumang = 0.; + sumang = 0; if (nv < 3) - return (0.); + return 0; if (DEBUG >= 4) { for (i = 0; i < nv; ++i) { @@ -878,7 +874,7 @@ double Girard() { for (i = 0; i < nv; ++i) { Cross(&V[i], &V[(i + 1) % nv], &side[i]); - (void) Normalize(&side[i]); + Normalize(&side[i]); } for (i = 0; i < nv; ++i) { @@ -935,7 +931,7 @@ double Girard() { fflush(stdout); } - return (area); + return area; } /* diff --git a/reproject/overlapArea.h b/reproject/overlapArea.h new file mode 100644 index 000000000..9d86d7b57 --- /dev/null +++ b/reproject/overlapArea.h @@ -0,0 +1,7 @@ +#ifndef OVERLAP_AREA +#define OVERLAP_AREA + +double computeOverlap(double *ilon, double *ilat, double *olon, double *olat, + int energyMode, double refArea, double *areaRatio); + +#endif diff --git a/reproject/setup_package.py b/reproject/setup_package.py index 26fc31729..7f91248da 100644 --- a/reproject/setup_package.py +++ b/reproject/setup_package.py @@ -1,17 +1,24 @@ +import os from distutils.core import Extension +REPROJECT_ROOT = os.path.relpath(os.path.dirname(__file__)) def get_extensions(): + sources = [os.path.join(REPROJECT_ROOT, "_overlap.pyx")] + include_dirs = ['numpy'] + libraries = [] - from numpy import get_include as get_numpy_include + sources.append("reproject/overlapArea.c") + include_dirs.append('reproject') - numpy_includes = get_numpy_include() + extension = Extension( + name="reproject._overlap", + sources=sources, + include_dirs=include_dirs, + libraries=libraries, + language="c",) - ext_modules = [Extension("reproject._overlap_wrapper", - ['reproject/_overlap_wrapper.c', 'reproject/overlapArea.c'], - include_dirs=[numpy_includes])] - - return ext_modules + return [extension] def requires_2to3(): diff --git a/reproject/tests/test_overlap.py b/reproject/tests/test_overlap.py index 93a47b31e..86ec6e4ac 100644 --- a/reproject/tests/test_overlap.py +++ b/reproject/tests/test_overlap.py @@ -3,15 +3,15 @@ def test_overlap(): EPS = np.radians(1e-2) - lon, lat = [0, EPS, EPS, 0], [0, 0, EPS, EPS] + lon, lat = np.array([[0, EPS, EPS, 0]]), np.array([[0, 0, EPS, EPS]]) overlap, area_ratio = compute_overlap(lon, lat, lon, lat) np.testing.assert_allclose(overlap, EPS ** 2, rtol=1e-6) np.testing.assert_allclose(area_ratio, 1, rtol=1e-6) def test_overlap2(): EPS = np.radians(1e-2) - ilon, ilat = [0, EPS, EPS, 0], [0, 0, EPS, EPS] - olon, olat = [0.5 * EPS, 1.5 * EPS, 1.5 * EPS, 0.5 * EPS], [0, 0, EPS, EPS] + ilon, ilat = np.array([[0, EPS, EPS, 0]]), np.array([[0, 0, EPS, EPS]]) + olon, olat = np.array([[0.5 * EPS, 1.5 * EPS, 1.5 * EPS, 0.5 * EPS]]), np.array([[0, 0, EPS, EPS]]) overlap, area_ratio = compute_overlap(ilon, ilat, olon, olat) np.testing.assert_allclose(overlap, 0.5 * EPS ** 2, rtol=1e-6)