MRG: add resampling and smoothing functions

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);

nibabel/tests/data/to_canonical.m

@@ -0,0 +1,61 @@
+function to_canonical(imgs, vox_sizes, prefix, hold)
+% Resample images to canonical (transverse) orientation with given voxel sizes
+%
+% Inspired by ``reorient.m`` by John Ashburner:
+% http://blogs.warwick.ac.uk/files/nichols/reorient.m
+%
+% Parameters
+% ----------
+% imgs : char or cell array or struct array
+%     Images to resample to canonical orientation.
+% vox_sizes : vector (3, 1), optional
+%     Voxel sizes for output image.
+% prefix : char, optional
+%     Prefix for output resampled images, default = 'r'
+% hold : float, optional
+%     Hold (resampling method) value, default = 3.
+
+if ~isstruct(imgs)
+    imgs = spm_vol(imgs);
+end
+if nargin < 2
+    vox_sizes = [1 1 1];
+elseif numel(vox_sizes) == 1
+    vox_sizes = [vox_sizes vox_sizes vox_sizes];
+end
+vox_sizes = vox_sizes(:);
+if nargin < 3
+    prefix = 'r';
+end
+if nargin < 4
+    hold = 3;
+end
+
+for vol_no = 1:numel(imgs)
+    vol = imgs{vol_no}(1);
+    % From:
+    % http://stackoverflow.com/questions/4165859/generate-all-possible-combinations-of-the-elements-of-some-vectors-cartesian-pr
+    sets = {[1, vol.dim(1)], [1, vol.dim(2)], [1, vol.dim(3)]};
+    [x y z] = ndgrid(sets{:});
+    corners = [x(:) y(:) z(:)];
+    corner_coords = [corners ones(length(corners), 1)]';
+    corner_mm = vol.mat * corner_coords;
+    min_xyz = min(corner_mm(1:3, :), [], 2);
+    max_xyz = max(corner_mm(1:3, :), [], 2);
+    % Make output volume
+    out_vol = vol;
+    out_vol.private = [];
+    out_vol.mat = diag([vox_sizes' 1]);
+    out_vol.mat(1:3, 4) = min_xyz - vox_sizes;
+    out_vol.dim(1:3) = ceil((max_xyz - min_xyz) ./ vox_sizes) + 1;
+    [dpath, froot, ext] = fileparts(vol.fname);
+    out_vol.fname = fullfile(dpath, [prefix froot ext]);
+    out_vol = spm_create_vol(out_vol);
+    % Resample original volume at output volume grid
+    plane_size = out_vol.dim(1:2);
+    for slice_no = 1:out_vol.dim(3)
+        resamp_affine = inv(spm_matrix([0 0 -slice_no]) * inv(out_vol.mat) * vol.mat);
+        slice_vals = spm_slice_vol(vol, resamp_affine, plane_size, hold);
+        out_vol = spm_write_plane(out_vol, slice_vals, slice_no);
+    end
+end

nibabel/tests/test_affines.py

@@ -3,8 +3,8 @@
 
 import numpy as np
 
-from ..affines import (apply_affine, append_diag, to_matvec, from_matvec,
-                       dot_reduce)
+from ..affines import (AffineError, apply_affine, append_diag, to_matvec,
+                       from_matvec, dot_reduce)
 
 
 from nose.tools import assert_equal, assert_raises
@@ -123,7 +123,7 @@ def test_append_diag():
                         [0, 0, 0, 5, 9],
                         [0, 0, 0, 0, 1]])
     # Length of starts has to match length of steps
-    assert_raises(ValueError, append_diag, aff, [5, 6], [9])
+    assert_raises(AffineError, append_diag, aff, [5, 6], [9])
 
 
 def test_dot_reduce():

nibabel/tests/test_spaces.py

@@ -34,17 +34,14 @@ def assert_all_in(in_shape, in_affine, out_shape, out_affine):
     assert_true(np.all(out_grid < np.array(out_shape) + TINY))
 
 
-def test_vox2out_vox():
-    # Test world space bounding box
-    # Test basic case, identity, no voxel sizes passed
-    shape, aff = vox2out_vox(((2, 3, 4), np.eye(4)))
-    assert_array_equal(shape, (2, 3, 4))
-    assert_array_equal(aff, np.eye(4))
+def get_outspace_params():
+    # Return in_shape, in_aff, vox, out_shape, out_aff for output space tests
+    # Put in function to use also for resample_to_output tests
     # Some affines as input to the tests
     trans_123 = [[1, 0, 0, 1], [0, 1, 0, 2], [0, 0, 1, 3], [0, 0, 0, 1]]
     trans_m123 = [[1, 0, 0, -1], [0, 1, 0, -2], [0, 0, 1, -3], [0, 0, 0, 1]]
     rot_3 = from_matvec(euler2mat(np.pi / 4), [0, 0, 0])
-    for in_shape, in_aff, vox, out_shape, out_aff in (
+    return ( # in_shape, in_aff, vox, out_shape, out_aff
         # Identity
         ((2, 3, 4), np.eye(4), None, (2, 3, 4), np.eye(4)),
         # Flip first axis
@@ -80,7 +77,16 @@ def test_vox2out_vox():
         # Number of voxel sizes matches length
         ((2, 3), np.diag([4, 5, 6, 1]), (4, 5),
          (2, 3), np.diag([4, 5, 1, 1])),
-    ):
+    )
+
+
+def test_vox2out_vox():
+    # Test world space bounding box
+    # Test basic case, identity, no voxel sizes passed
+    shape, aff = vox2out_vox(((2, 3, 4), np.eye(4)))
+    assert_array_equal(shape, (2, 3, 4))
+    assert_array_equal(aff, np.eye(4))
+    for in_shape, in_aff, vox, out_shape, out_aff in get_outspace_params():
         img = Nifti1Image(np.ones(in_shape), in_aff)
         for input in ((in_shape, in_aff), img):
             shape, aff = vox2out_vox(input, vox)