Merge pull request nipy#255 from nipy/add-processing

MRG: add resampling and smoothing functions

This pull request adds some simple processing functions when scipy
is available. The functions are:

* resample_from_to - resample first image into space of second
* resample_to_output - resample image into output (world) space
* smooth_image - smooth image over voxel axes

I also put in the utility functions converting between FWHM and sigma.

nibabel/affines.py

@@ -8,6 +8,13 @@
 from .externals.six.moves import reduce
 
 
+class AffineError(ValueError):
+    """ Errors in calculating or using affines """
+    # Inherits from ValueError to keep compatibility with ValueError previously
+    # raised in append_diag
+    pass
+
+
 def apply_affine(aff, pts):
     """ Apply affine matrix `aff` to points `pts`
 
@@ -213,7 +220,7 @@ def append_diag(aff, steps, starts=()):
     if len(starts) == 0:
         starts = np.zeros(n_steps, dtype=steps.dtype)
     elif len(starts) != n_steps:
-        raise ValueError('Steps should have same length as starts')
+        raise AffineError('Steps should have same length as starts')
     old_n_out, old_n_in = aff.shape[0] - 1, aff.shape[1] - 1
     # make new affine
     aff_plus = np.zeros((old_n_out + n_steps + 1,

nibabel/processing.py

@@ -0,0 +1,304 @@
+# 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.
+#
+### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
+""" Image processing functions for:
+
+* smoothing
+* resampling
+* converting sd to and from FWHM
+
+Smoothing and resampling routines need scipy
+"""
+from __future__ import print_function, division, absolute_import
+
+import numpy as np
+import numpy.linalg as npl
+
+from .optpkg import optional_package
+spnd, _, _ = optional_package('scipy.ndimage')
+
+from .affines import AffineError, to_matvec, from_matvec, append_diag
+from .spaces import vox2out_vox
+from .nifti1 import Nifti1Image
+
+SIGMA2FWHM = np.sqrt(8 * np.log(2))
+
+
+def fwhm2sigma(fwhm):
+    """ Convert a FWHM value to sigma in a Gaussian kernel.
+
+    Parameters
+    ----------
+    fwhm : array-like
+       FWHM value or values
+
+    Returns
+    -------
+    sigma : array or float
+       sigma values corresponding to `fwhm` values
+
+    Examples
+    --------
+    >>> sigma = fwhm2sigma(6)
+    >>> sigmae = fwhm2sigma([6, 7, 8])
+    >>> sigma == sigmae[0]
+    True
+    """
+    return np.asarray(fwhm) / SIGMA2FWHM
+
+
+def sigma2fwhm(sigma):
+    """ Convert a sigma in a Gaussian kernel to a FWHM value
+
+    Parameters
+    ----------
+    sigma : array-like
+       sigma value or values
+
+    Returns
+    -------
+    fwhm : array or float
+       fwhm values corresponding to `sigma` values
+
+    Examples
+    --------
+    >>> fwhm = sigma2fwhm(3)
+    >>> fwhms = sigma2fwhm([3, 4, 5])
+    >>> fwhm == fwhms[0]
+    True
+    """
+    return np.asarray(sigma) * SIGMA2FWHM
+
+
+def adapt_affine(affine, n_dim):
+    """ Adapt input / output dimensions of spatial `affine` for `n_dims`
+
+    Adapts a spatial (4, 4) affine that is being applied to an image with fewer
+    than 3 spatial dimensions, or more than 3 dimensions.  If there are more
+    than three dimensions, assume an identity transformation for these
+    dimensions.
+
+    Parameters
+    ----------
+    affine : array-like
+        affine transform. Usually shape (4, 4).  For what follows ``N, M =
+        affine.shape``
+    n_dims : int
+        Number of dimensions of underlying array, and therefore number of input
+        dimensions for affine.
+
+    Returns
+    -------
+    adapted : shape (M, n_dims+1) array
+        Affine array adapted to number of input dimensions.  Columns of the
+        affine corresponding to missing input dimensions have been dropped,
+        columns corresponding to extra input dimensions have an extra identity
+        column added
+    """
+    affine = np.asarray(affine)
+    rzs, trans = to_matvec(affine)
+    # For missing input dimensions, drop columns in rzs
+    rzs = rzs[:, :n_dim]
+    adapted = from_matvec(rzs, trans)
+    n_extra_columns = n_dim - adapted.shape[1] + 1
+    if n_extra_columns > 0:
+        adapted = append_diag(adapted, np.ones((n_extra_columns,)))
+    return adapted
+
+
+def resample_from_to(from_img,
+                     to_vox_map,
+                     order=3,
+                     mode='constant',
+                     cval=0.,
+                     out_class=Nifti1Image):
+    """ Resample image `from_img` to mapped voxel space `to_vox_map`
+
+    Resample using N-d spline interpolation.
+
+    Parameters
+    ----------
+    from_img : object
+        Object having attributes ``dataobj``, ``affine``, ``header``. If
+        `out_class` is not None, ``img.__class__`` should be able to construct
+        an image from data, affine and header.
+    to_vox_map : image object or length 2 sequence
+        If object, has attributes ``shape`` giving input voxel shape, and
+        ``affine`` giving mapping of input voxels to output space. If length 2
+        sequence, elements are (shape, affine) with same meaning as above. The
+        affine is a (4, 4) array-like.
+    order : int, optional
+        The order of the spline interpolation, default is 3.  The order has to
+        be in the range 0-5 (see ``scipy.ndimage.affine_transform``)
+    mode : str, optional
+        Points outside the boundaries of the input are filled according
+        to the given mode ('constant', 'nearest', 'reflect' or 'wrap').
+        Default is 'constant' (see ``scipy.ndimage.affine_transform``)
+    cval : scalar, optional
+        Value used for points outside the boundaries of the input if
+        ``mode='constant'``. Default is 0.0 (see
+        ``scipy.ndimage.affine_transform``)
+    out_class : None or SpatialImage class, optional
+        Class of output image.  If None, use ``from_img.__class__``.
+
+    Returns
+    -------
+    out_img : object
+        Image of instance specified by `out_class`, containing data output from
+        resampling `from_img` into axes aligned to the output space of
+        ``from_img.affine``
+    """
+    try:
+        to_shape, to_affine = to_vox_map.shape, to_vox_map.affine
+    except AttributeError:
+        to_shape, to_affine = to_vox_map
+    a_to_affine = adapt_affine(to_affine, len(to_shape))
+    if out_class is None:
+        out_class = from_img.__class__
+    from_n_dim = len(from_img.shape)
+    if from_n_dim < 3:
+        raise AffineError('from_img must be at least 3D')
+    a_from_affine = adapt_affine(from_img.affine, from_n_dim)
+    to_vox2from_vox = npl.inv(a_from_affine).dot(a_to_affine)
+    rzs, trans = to_matvec(to_vox2from_vox)
+    data = spnd.affine_transform(from_img.dataobj,
+                                 rzs,
+                                 trans,
+                                 to_shape,
+                                 order=order,
+                                 mode=mode,
+                                 cval=cval)
+    return out_class(data, to_affine, from_img.header)
+
+
+def resample_to_output(in_img,
+                       voxel_sizes=None,
+                       order=3,
+                       mode='constant',
+                       cval=0.,
+                       out_class=Nifti1Image):
+    """ Resample image `in_img` to output voxel axes (world space)
+
+    Parameters
+    ----------
+    in_img : object
+        Object having attributes ``dataobj``, ``affine``, ``header``. If
+        `out_class` is not None, ``img.__class__`` should be able to construct
+        an image from data, affine and header.
+    voxel_sizes : None or sequence
+        Gives the diagonal entries of ``out_img.affine` (except the trailing 1
+        for the homogenous coordinates) (``out_img.affine ==
+        np.diag(voxel_sizes + [1])``). If None, return identity
+        `out_img.affine`.  If scalar, interpret as vector ``[voxel_sizes] *
+        len(in_img.shape)``.
+    order : int, optional
+        The order of the spline interpolation, default is 3.  The order has to
+        be in the range 0-5 (see ``scipy.ndimage.affine_transform``).
+    mode : str, optional
+        Points outside the boundaries of the input are filled according to the
+        given mode ('constant', 'nearest', 'reflect' or 'wrap').  Default is
+        'constant' (see ``scipy.ndimage.affine_transform``).
+    cval : scalar, optional
+        Value used for points outside the boundaries of the input if
+        ``mode='constant'``. Default is 0.0 (see
+        ``scipy.ndimage.affine_transform``).
+    out_class : None or SpatialImage class, optional
+        Class of output image.  If None, use ``in_img.__class__``.
+
+    Returns
+    -------
+    out_img : object
+        Image of instance specified by `out_class`, containing data output from
+        resampling `in_img` into axes aligned to the output space of
+        ``in_img.affine``
+    """
+    if out_class is None:
+        out_class = in_img.__class__
+    in_shape = in_img.shape
+    n_dim = len(in_shape)
+    if voxel_sizes is not None:
+        voxel_sizes = np.asarray(voxel_sizes)
+        if voxel_sizes.ndim == 0:  # Scalar
+            voxel_sizes = np.repeat(voxel_sizes, n_dim)
+    # Allow 2D images by promoting to 3D.  We might want to see what a slice
+    # looks like when resampled into world coordinates
+    if n_dim < 3:  # Expand image to 3D, make voxel sizes match
+        new_shape = in_shape + (1,) * (3 - n_dim)
+        data = in_img.get_data().reshape(new_shape)  # 2D data should be small
+        in_img = out_class(data, in_img.affine, in_img.header)
+        if voxel_sizes is not None and len(voxel_sizes) == n_dim:
+            # Need to pad out voxel sizes to match new image dimensions
+            voxel_sizes = tuple(voxel_sizes) + (1,) * (3 - n_dim)
+    out_vox_map = vox2out_vox((in_img.shape, in_img.affine), voxel_sizes)
+    return resample_from_to(in_img, out_vox_map, order, mode, cval, out_class)
+
+
+def smooth_image(img,
+                 fwhm,
+                 mode='nearest',
+                 cval=0.,
+                 out_class=Nifti1Image):
+    """ Smooth image `img` along voxel axes by FWHM `fwhm` millimeters
+
+    Parameters
+    ----------
+    img : object
+        Object having attributes ``dataobj``, ``affine``, ``header``. If
+        `out_class` is not None, ``img.__class__`` should be able to construct
+        an image from data, affine and header.
+    fwhm : scalar or length 3 sequence
+        FWHM *in mm* over which to smooth.  The smoothing applies to the voxel
+        axes, not to the output axes, but is in millimeters.  The function
+        adjusts the FWHM to voxels using the voxel sizes calculated from the
+        affine. A scalar implies the same smoothing across the spatial
+        dimensions of the image, but 0 smoothing over any further dimensions
+        such as time.  A vector should be the same length as the number of
+        image dimensions.
+    mode : str, optional
+        Points outside the boundaries of the input are filled according
+        to the given mode ('constant', 'nearest', 'reflect' or 'wrap').
+        Default is 'nearest'. This is different from the default for
+        ``scipy.ndimage.affine_transform``, which is 'constant'. 'nearest'
+        might be a better choice when smoothing to the edge of an image where
+        there is still strong brain signal, otherwise this signal will get
+        blurred towards zero.
+    cval : scalar, optional
+        Value used for points outside the boundaries of the input if
+        ``mode='constant'``. Default is 0.0 (see
+        ``scipy.ndimage.affine_transform``).
+    out_class : None or SpatialImage class, optional
+        Class of output image.  If None, use ``img.__class__``.
+
+    Returns
+    -------
+    smoothed_img : object
+        Image of instance specified by `out_class`, containing data output from
+        smoothing `img` data by given FWHM kernel.
+    """
+    if out_class is None:
+        out_class = img.__class__
+    n_dim = len(img.shape)
+    # TODO: make sure time axis is last
+    # Pad out fwhm from scalar, adding 0 for fourth etc (time etc) dimensions
+    fwhm = np.asarray(fwhm)
+    if fwhm.size == 1:
+        fwhm_scalar = fwhm
+        fwhm = np.zeros((n_dim,))
+        fwhm[:3] = fwhm_scalar
+    # Voxel sizes
+    RZS = img.affine[:-1, :n_dim]
+    vox = np.sqrt(np.sum(RZS ** 2, 0))
+    # Smoothing in terms of voxels
+    vox_fwhm = fwhm / vox
+    vox_sd = fwhm2sigma(vox_fwhm)
+    # Do the smoothing
+    sm_data = spnd.gaussian_filter(img.dataobj,
+                                   vox_sd,
+                                   mode=mode,
+                                   cval=cval)
+    return out_class(sm_data, img.affine, img.header)

nibabel/testing/__init__.py

@@ -37,7 +37,7 @@ def assert_dt_equal(a, b):
     assert_equal(np.dtype(a).str, np.dtype(b).str)
 
 
-def assert_allclose_safely(a, b, match_nans=True):
+def assert_allclose_safely(a, b, match_nans=True, rtol=1e-5, atol=1e-8):
     """ Allclose in integers go all wrong for large integers
     """
     a = np.atleast_1d(a)  # 0d arrays cannot be indexed
@@ -57,7 +57,7 @@ def assert_allclose_safely(a, b, match_nans=True):
         a = a.astype(float)
     if b.dtype.kind in 'ui':
         b = b.astype(float)
-    assert_true(np.allclose(a, b))
+    assert_true(np.allclose(a, b, rtol=rtol, atol=atol))
 
 
 def assert_re_in(regex, c, flags=0):

nibabel/tests/data/.gitignore

@@ -0,0 +1,2 @@
+anat_moved.nii
+resampled_functional.nii

nibabel/tests/data/make_moved_anat.py

@@ -0,0 +1,22 @@
+""" Make anatomical image with altered affine
+
+* Add some rotations and translations to affine;
+* Save as ``.nii`` file so SPM can read it.
+
+See ``resample_using_spm.m`` for processing of this generated image by SPM.
+"""
+
+import numpy as np
+
+import nibabel as nib
+from nibabel.eulerangles import euler2mat
+from nibabel.affines import from_matvec
+
+img = nib.load('anatomical.nii')
+some_rotations = euler2mat(0.1, 0.2, 0.3)
+extra_affine = from_matvec(some_rotations, [3, 4, 5])
+moved_anat = nib.Nifti1Image(img.dataobj,
+                             extra_affine.dot(img.affine),
+                             img.header)
+moved_anat.set_data_dtype(np.float32)
+nib.save(moved_anat, 'anat_moved.nii')

nibabel/tests/data/resample_using_spm.m

@@ -0,0 +1,15 @@
+% Script uses SPM to resample moved anatomical image.
+%
+% Run `python make_moved_anat.py` to generate file to work on.
+%
+% Run from the directory containing this file.
+% Works with Octave or MATLAB.
+% Needs SPM (5, 8 or 12) on the MATLAB path.
+P = {'functional.nii', 'anat_moved.nii'};
+% Resample without masking
+flags = struct('mask', false, 'mean', false, ...
+               'interp', 1, 'which', 1, ...
+               'prefix', 'resampled_');
+spm_reslice(P, flags);
+% Reorient to canonical orientation at 4mm resolution, polynomial interpolation
+to_canonical({'anat_moved.nii'}, 4, 'reoriented_', 1);