add bspline cython extension

nipy · oesteban · Sep 20, 2019 · Sep 20, 2019 · Sep 22, 2019 · Sep 23, 2019
commit 2acd9f39c4500181ef5c52c238bcd880d0cd1f29
diff --git a/.gitignore b/.gitignore
@@ -84,3 +84,4 @@ Thumbs.db
 doc/source/reference
 venv/
 .buildbot.patch
+nibabel/maths/bspline.c
diff --git a/nibabel/maths/__init__.py b/nibabel/maths/__init__.py
@@ -0,0 +1,17 @@
+# emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
+#
+#   See COPYING file distributed along with the NiBabel package for the
+#   copyright and license terms.
+#
+### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
+"""Geometric transforms
+
+.. currentmodule:: nibabel.maths
+
+.. autosummary::
+   :toctree: ../generated
+
+   maths
+"""
diff --git a/nibabel/maths/bspline.pyx b/nibabel/maths/bspline.pyx
@@ -0,0 +1,43 @@
+# emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
+#
+#   See COPYING file distributed along with the NiBabel package for the
+#   copyright and license terms.
+#
+### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
+''' Cython math extension '''
+from __future__ import division
+import numpy as np
+cimport numpy as np
+cimport cython
+from libc.math cimport fabs
+
+
+cdef double c_cubic(double x) nogil:
+    cdef:
+        double x_t = fabs(x)
+
+    if x_t >= 2.0:
+        return(0.0)
+    if x_t <= 1.0:
+        return(2.0 / 3.0 - x_t**2 + 0.5 * x_t**3)
+    elif x_t <= 2.0:
+        return((2 - x_t)**3 / 6.0)
+
+
+def cubic(double x):
+    """
+    Evaluate the univariate cubic bspline at x
+
+    Pure python implementation: ::
+
+        def bspl(x):
+            if x >= 2.0:
+                return 0.0
+            if x <= 1.0:
+                return 2.0 / 3.0 - x**2 + 0.5 * x**3
+            elif x <= 2.0:
+                return (2 - x)**3 / 6.0
+    """
+    return(c_cubic(x))
diff --git a/nibabel/transform/nonlinear.py b/nibabel/transform/nonlinear.py
@@ -105,3 +105,31 @@ def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True,
 
     def map_voxel(self, index, moving=None):
         return tuple(self._moving[index + self.__s])
+
+    def map_coordinates(self, coordinates, order=3, mode='constant', cval=0.0,
+                        prefilter=True):
+        coordinates = np.array(coordinates)
+        # Extract shapes and dimensions, then flatten
+        ndim = coordinates.shape[-1]
+        output_shape = coordinates.shape[:-1]
+        flatcoord = np.moveaxis(coordinates, -1, 0).reshape(ndim, -1)
+
+        # Convert coordinates to voxel indices
+        ijk = np.tensordot(
+            np.linalg.inv(self.reference.affine),
+            np.vstack((flatcoord, np.ones((1, flatcoord.shape[1])))),
+            axes=1)
+        deltas = ndi.map_coordinates(
+            self._field,
+            ijk,
+            order=order,
+            mode=mode,
+            cval=cval,
+            prefilter=prefilter)
+
+        print(deltas)
+
+        deltas = np.moveaxis(deltas[0:3, :].reshape((ndim, ) + output_shape),
+                             0, -1)
+
+        return coordinates + deltas