Use Cython to wrap C overlap area code

.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

.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

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

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

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) 

reproject/overlapArea.c

@@ -6,6 +6,7 @@
 #include <stdio.h>
 #include <math.h>
 #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;
 }
 
 /*

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

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():

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)