Merge pull request #24 from oesteban/enh/22

ENH: Uber-refactor of code style, method names, etc.

Author: oesteban

nitransforms/base.py

@@ -15,13 +15,14 @@
 EQUALITY_TOL = 1e-5
 
 
-class ImageSpace(object):
+class ImageGrid(object):
     """Class to represent spaces of gridded data (images)."""
 
     __slots__ = ['_affine', '_shape', '_ndim', '_ndindex', '_coords', '_nvox',
                  '_inverse']
 
     def __init__(self, image):
+        """Create a gridded sampling reference."""
         self._affine = image.affine
         self._shape = image.shape
         self._ndim = len(image.shape)
@@ -34,26 +35,32 @@ def __init__(self, image):
 
     @property
     def affine(self):
+        """Access the indexes-to-RAS affine."""
         return self._affine
 
     @property
     def inverse(self):
+        """Access the RAS-to-indexes affine."""
         return self._inverse
 
     @property
     def shape(self):
+        """Access the space's size of each dimension."""
         return self._shape
 
     @property
     def ndim(self):
+        """Access the number of dimensions."""
         return self._ndim
 
     @property
     def nvox(self):
+        """Access the total number of voxels."""
         return self._nvox
 
     @property
     def ndindex(self):
+        """List the indexes corresponding to the space grid."""
         if self._ndindex is None:
             indexes = tuple([np.arange(s) for s in self._shape])
             self._ndindex = np.array(np.meshgrid(
@@ -62,6 +69,7 @@ def ndindex(self):
 
     @property
     def ndcoords(self):
+        """List the physical coordinates of this gridded space samples."""
         if self._coords is None:
             self._coords = np.tensordot(
                 self._affine,
@@ -70,14 +78,15 @@ def ndcoords(self):
             )[:3, ...]
         return self._coords
 
-    def map_voxels(self, coordinates):
-        coordinates = np.array(coordinates)
-        ncoords = coordinates.shape[-1]
-        coordinates = np.vstack((coordinates, np.ones((1, ncoords))))
+    def ras(self, ijk):
+        """Get RAS+ coordinates from input indexes."""
+        ras = self._affine.dot(_as_homogeneous(ijk).T)[:3, ...]
+        return ras.T
 
-        # Back to grid coordinates
-        return np.tensordot(np.linalg.inv(self._affine),
-                            coordinates, axes=1)[:3, ...]
+    def index(self, x):
+        """Get the image array's indexes corresponding to coordinates."""
+        ijk = self._inverse.dot(_as_homogeneous(x).T)[:3, ...]
+        return ijk.T
 
     def _to_hdf5(self, group):
         group.attrs['Type'] = 'image'
@@ -86,14 +95,9 @@ def _to_hdf5(self, group):
         group.create_dataset('shape', data=self.shape)
 
     def __eq__(self, other):
-        try:
-            return (
-                np.allclose(self.affine, other.affine, rtol=EQUALITY_TOL)
-                and self.shape == other.shape
-            )
-        except AttributeError:
-            pass
-        return False
+        """Overload equals operator."""
+        return (np.allclose(self.affine, other.affine, rtol=EQUALITY_TOL) and
+                self.shape == other.shape)
 
 
 class TransformBase(object):
@@ -102,6 +106,7 @@ class TransformBase(object):
     __slots__ = ['_reference']
 
     def __init__(self):
+        """Instantiate a transform."""
         self._reference = None
 
     def __eq__(self, other):
@@ -110,25 +115,31 @@ def __eq__(self, other):
             return False
         return np.allclose(self.matrix, other.matrix, rtol=EQUALITY_TOL)
 
+    def __call__(self, x, inverse=False, index=0):
+        """Apply y = f(x)."""
+        return self.map(x, inverse=inverse, index=index)
+
     @property
     def reference(self):
-        '''A reference space where data will be resampled onto'''
+        """Access a reference space where data will be resampled onto."""
         if self._reference is None:
             raise ValueError('Reference space not set')
         return self._reference
 
     @reference.setter
     def reference(self, image):
-        self._reference = ImageSpace(image)
+        self._reference = ImageGrid(image)
 
     @property
     def ndim(self):
+        """Access the dimensions of the reference space."""
         return self.reference.ndim
 
     def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True,
                  output_dtype=None):
         """
         Resample the moving image in reference space.
+
         Parameters
         ----------
         moving : `spatialimage`
@@ -150,43 +161,58 @@ def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True,
             slightly blurred if *order > 1*, unless the input is prefiltered,
             i.e. it is the result of calling the spline filter on the original
             input.
+
         Returns
         -------
         moved_image : `spatialimage`
             The moving imaged after resampling to reference space.
+
         """
         moving_data = np.asanyarray(moving.dataobj)
         if output_dtype is None:
             output_dtype = moving_data.dtype
 
         moved = ndi.geometric_transform(
             moving_data,
-            mapping=self.map_voxel,
+            mapping=self._map_index,
             output_shape=self.reference.shape,
             output=output_dtype,
             order=order,
             mode=mode,
             cval=cval,
             prefilter=prefilter,
-            extra_keywords={'moving': moving},
+            extra_keywords={'moving': ImageGrid(moving)},
         )
 
         moved_image = moving.__class__(moved, self.reference.affine, moving.header)
         moved_image.header.set_data_dtype(output_dtype)
         return moved_image
 
-    def map_point(self, coords):
-        """Apply y = f(x), where x is the argument `coords`."""
-        raise NotImplementedError
+    def map(self, x, inverse=False, index=0):
+        r"""
+        Apply :math:`y = f(x)`.
 
-    def map_voxel(self, index, moving=None):
-        """Apply ijk' = f_ijk((i, j, k)), equivalent to the above with indexes."""
-        raise NotImplementedError
+        Parameters
+        ----------
+        x : N x D numpy.ndarray
+            Input RAS+ coordinates (i.e., physical coordinates).
+        inverse : bool
+            If ``True``, apply the inverse transform :math:`x = f^{-1}(y)`.
+        index : int, optional
+            Transformation index
 
-    def _to_hdf5(self, x5_root):
-        """Serialize this object into the x5 file format."""
+        Returns
+        -------
+        y : N x D numpy.ndarray
+            Transformed (mapped) RAS+ coordinates (i.e., physical coordinates).
+
+        """
         raise NotImplementedError
 
+    def _map_index(self, ijk, moving):
+        x = self.reference.ras(_as_homogeneous(ijk))
+        return moving.index(self.map(x))
+
     def to_filename(self, filename, fmt='X5'):
         """Store the transform in BIDS-Transforms HDF5 file format (.x5)."""
         with h5py.File(filename, 'w') as out_file:
@@ -196,3 +222,28 @@ def to_filename(self, filename, fmt='X5'):
             self._to_hdf5(root)
 
         return filename
+
+    def _to_hdf5(self, x5_root):
+        """Serialize this object into the x5 file format."""
+        raise NotImplementedError
+
+
+def _as_homogeneous(xyz, dtype='float32'):
+    """
+    Convert 2D and 3D coordinates into homogeneous coordinates.
+
+    Examples
+    --------
+    >>> _as_homogeneous((4, 5), dtype='int8').tolist()
+    [[4, 5, 1]]
+
+    >>> _as_homogeneous((4, 5, 6),dtype='int8').tolist()
+    [[4, 5, 6, 1]]
+
+    >>> _as_homogeneous([(1, 2, 3), (4, 5, 6)]).tolist()
+    [[1.0, 2.0, 3.0, 1.0], [4.0, 5.0, 6.0, 1.0]]
+
+
+    """
+    xyz = np.atleast_2d(np.array(xyz, dtype=dtype))
+    return np.hstack((xyz, np.ones((xyz.shape[0], 1), dtype=dtype)))

nitransforms/conftest.py

@@ -39,7 +39,7 @@ def data_path():
 
 
 @pytest.fixture
-def get_data():
+def get_testdata():
     """Generate data in the requested orientation."""
     global _data
 

nitransforms/linear.py

@@ -143,7 +143,7 @@ def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True,
             nvols = moving.shape[3]
 
         movaff = moving.affine
-        movingdata = moving.get_data()
+        movingdata = np.asanyarray(moving.dataobj)
         if nvols == 1:
             movingdata = movingdata[..., np.newaxis]
 
@@ -178,42 +178,45 @@ def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True,
 
         return moved_image
 
-    def map_point(self, coords, index=0, forward=True):
-        """
-        Apply y = f(x), where x is the argument `coords`.
+    def map(self, x, inverse=False, index=0):
+        r"""
+        Apply :math:`y = f(x)`.
 
         Parameters
         ----------
-        coords : array_like
-            RAS coordinates to map
+        x : N x D numpy.ndarray
+            Input RAS+ coordinates (i.e., physical coordinates).
+        inverse : bool
+            If ``True``, apply the inverse transform :math:`x = f^{-1}(y)`.
         index : int, optional
             Transformation index
-        forward: bool, optional
-            Direction of mapping. Default is set to ``True``. If ``False``,
-            the inverse transformation is applied.
 
         Returns
         -------
-        out: ndarray
-            Transformed coordinates
+        y : N x D numpy.ndarray
+            Transformed (mapped) RAS+ coordinates (i.e., physical coordinates).
 
         Examples
         --------
         >>> xfm = Affine([[1, 0, 0, 1], [0, 1, 0, 2], [0, 0, 1, 3], [0, 0, 0, 1]])
-        >>> xfm.map_point((0,0,0))
+        >>> xfm.map((0,0,0))
         array([1, 2, 3])
 
-        >>> xfm.map_point((0,0,0), forward=False)
+        >>> xfm.map((0,0,0), inverse=True)
         array([-1., -2., -3.])
 
         """
-        coords = np.array(coords)
+        coords = np.array(x)
         if coords.shape[0] == self._matrix[index].shape[0] - 1:
             coords = np.append(coords, [1])
-        affine = self._matrix[index] if forward else np.linalg.inv(self._matrix[index])
+        affine = self._matrix[index]
+
+        if inverse is True:
+            affine = np.linalg.inv(self._matrix[index])
+
         return affine.dot(coords)[:-1]
 
-    def map_voxel(self, index, nindex=0, moving=None):
+    def _map_voxel(self, index, nindex=0, moving=None):
         """Apply ijk' = f_ijk((i, j, k)), equivalent to the above with indexes."""
         try:
             reference = self.reference