diff --git a/nibabel/affines.py b/nibabel/affines.py index 5ed07e5e39..325e5e8b6c 100644 --- a/nibabel/affines.py +++ b/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, diff --git a/nibabel/processing.py b/nibabel/processing.py new file mode 100644 index 0000000000..b1cc867496 --- /dev/null +++ b/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) diff --git a/nibabel/testing/__init__.py b/nibabel/testing/__init__.py index a2fa2a83f1..c28a54f6a3 100644 --- a/nibabel/testing/__init__.py +++ b/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): diff --git a/nibabel/tests/data/.gitignore b/nibabel/tests/data/.gitignore new file mode 100644 index 0000000000..215e61ce01 --- /dev/null +++ b/nibabel/tests/data/.gitignore @@ -0,0 +1,2 @@ +anat_moved.nii +resampled_functional.nii diff --git a/nibabel/tests/data/anatomical.nii b/nibabel/tests/data/anatomical.nii new file mode 100644 index 0000000000..2d48e4770d Binary files /dev/null and b/nibabel/tests/data/anatomical.nii differ diff --git a/nibabel/tests/data/functional.nii b/nibabel/tests/data/functional.nii new file mode 100644 index 0000000000..2768d4d391 Binary files /dev/null and b/nibabel/tests/data/functional.nii differ diff --git a/nibabel/tests/data/make_moved_anat.py b/nibabel/tests/data/make_moved_anat.py new file mode 100644 index 0000000000..6fba2d0902 --- /dev/null +++ b/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') diff --git a/nibabel/tests/data/reoriented_anat_moved.nii b/nibabel/tests/data/reoriented_anat_moved.nii new file mode 100644 index 0000000000..2f2411d115 Binary files /dev/null and b/nibabel/tests/data/reoriented_anat_moved.nii differ diff --git a/nibabel/tests/data/resample_using_spm.m b/nibabel/tests/data/resample_using_spm.m new file mode 100644 index 0000000000..350fcb8382 --- /dev/null +++ b/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); diff --git a/nibabel/tests/data/resampled_anat_moved.nii b/nibabel/tests/data/resampled_anat_moved.nii new file mode 100644 index 0000000000..c6b4549c21 Binary files /dev/null and b/nibabel/tests/data/resampled_anat_moved.nii differ diff --git a/nibabel/tests/data/to_canonical.m b/nibabel/tests/data/to_canonical.m new file mode 100644 index 0000000000..08d5abd327 --- /dev/null +++ b/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 diff --git a/nibabel/tests/test_affines.py b/nibabel/tests/test_affines.py index c950dfefac..a45b344470 100644 --- a/nibabel/tests/test_affines.py +++ b/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(): diff --git a/nibabel/tests/test_processing.py b/nibabel/tests/test_processing.py new file mode 100644 index 0000000000..72076918a4 --- /dev/null +++ b/nibabel/tests/test_processing.py @@ -0,0 +1,398 @@ +# 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. +# +### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## +""" Testing processing module +""" +from __future__ import division, print_function + +from os.path import dirname, join as pjoin + +import numpy as np +import numpy.linalg as npl + +from nibabel.optpkg import optional_package +spnd, have_scipy, _ = optional_package('scipy.ndimage') + +import nibabel as nib +from nibabel.processing import (sigma2fwhm, fwhm2sigma, adapt_affine, + resample_from_to, resample_to_output, smooth_image) +from nibabel.nifti1 import Nifti1Image +from nibabel.nifti2 import Nifti2Image +from nibabel.orientations import flip_axis, inv_ornt_aff +from nibabel.affines import AffineError, from_matvec, to_matvec, apply_affine +from nibabel.eulerangles import euler2mat + +from numpy.testing import (assert_almost_equal, + assert_array_equal) +from numpy.testing.decorators import skipif + +from nose.tools import (assert_true, assert_false, assert_raises, + assert_equal, assert_not_equal) + +from nibabel.tests.test_spaces import assert_all_in, get_outspace_params +from nibabel.testing import assert_allclose_safely + +needs_scipy = skipif(not have_scipy, 'These tests need scipy') + +DATA_DIR = pjoin(dirname(__file__), 'data') + +def test_sigma2fwhm(): + # Test from constant + assert_almost_equal(sigma2fwhm(1), 2.3548200) + assert_almost_equal(sigma2fwhm([1, 2, 3]), np.arange(1, 4) * 2.3548200) + assert_almost_equal(fwhm2sigma(2.3548200), 1) + assert_almost_equal(fwhm2sigma(np.arange(1, 4) * 2.3548200), [1, 2, 3]) + # direct test fwhm2sigma and sigma2fwhm are inverses of each other + fwhm = np.arange(1.0, 5.0, 0.1) + sigma = np.arange(1.0, 5.0, 0.1) + assert_true(np.allclose(sigma2fwhm(fwhm2sigma(fwhm)), fwhm)) + assert_true(np.allclose(fwhm2sigma(sigma2fwhm(sigma)), sigma)) + + +def test_adapt_affine(): + # Adapt affine to missing or extra input dimensions + aff_3d = from_matvec(np.arange(9).reshape((3, 3)), [11, 12, 13]) + # For 4x4 affine, 3D image, no-op + assert_array_equal(adapt_affine(aff_3d, 3), aff_3d) + # For 4x4 affine, 4D image, add extra identity dimension + assert_array_equal(adapt_affine(aff_3d, 4), + [[ 0, 1, 2, 0, 11], + [ 3, 4, 5, 0, 12], + [ 6, 7, 8, 0, 13], + [ 0, 0, 0, 1, 0], + [ 0, 0, 0, 0, 1]]) + # For 5x5 affine, 4D image, identity + aff_4d = from_matvec(np.arange(16).reshape((4, 4)), [11, 12, 13, 14]) + assert_array_equal(adapt_affine(aff_4d, 4), aff_4d) + # For 4x4 affine, 2D image, dropped column + assert_array_equal(adapt_affine(aff_3d, 2), + [[ 0, 1, 11], + [ 3, 4, 12], + [ 6, 7, 13], + [ 0, 0, 1]]) + # For 4x4 affine, 1D image, 2 dropped columnn + assert_array_equal(adapt_affine(aff_3d, 1), + [[ 0, 11], + [ 3, 12], + [ 6, 13], + [ 0, 1]]) + # For 3x3 affine, 2D image, identity + aff_2d = from_matvec(np.arange(4).reshape((2, 2)), [11, 12]) + assert_array_equal(adapt_affine(aff_2d, 2), aff_2d) + + +@needs_scipy +def test_resample_from_to(): + # Test resampling from image to image / image space + data = np.arange(24).reshape((2, 3, 4)) + affine = np.diag([-4, 5, 6, 1]) + img = Nifti1Image(data, affine) + img.header['descrip'] = 'red shirt image' + out = resample_from_to(img, img) + assert_almost_equal(img.dataobj, out.dataobj) + assert_array_equal(img.affine, out.affine) + # Check resampling reverses effect of flipping axes + # This will also test translations + flip_ornt = np.array([[0, 1], [1, 1], [2, 1]]) + for axis in (0, 1, 2): + ax_flip_ornt = flip_ornt.copy() + ax_flip_ornt[axis, 1] = -1 + aff_flip_i = inv_ornt_aff(ax_flip_ornt, (2, 3, 4)) + flipped_img = Nifti1Image(flip_axis(data, axis), + np.dot(affine, aff_flip_i)) + out = resample_from_to(flipped_img, ((2, 3, 4), affine)) + assert_almost_equal(img.dataobj, out.dataobj) + assert_array_equal(img.affine, out.affine) + # A translation of one voxel on each axis + trans_aff = from_matvec(np.diag([-4, 5, 6]), [4, -5, -6]) + trans_img = Nifti1Image(data, trans_aff) + out = resample_from_to(trans_img, img) + exp_out = np.zeros_like(data) + exp_out[:-1, :-1, :-1] = data[1:, 1:, 1:] + assert_almost_equal(out.dataobj, exp_out) + out = resample_from_to(img, trans_img) + trans_exp_out = np.zeros_like(data) + trans_exp_out[1:, 1:, 1:] = data[:-1, :-1, :-1] + assert_almost_equal(out.dataobj, trans_exp_out) + # Test mode with translation of first axis only + # Default 'constant' mode first + trans1_aff = from_matvec(np.diag([-4, 5, 6]), [4, 0, 0]) + trans1_img = Nifti1Image(data, trans1_aff) + out = resample_from_to(img, trans1_img) + exp_out = np.zeros_like(data) + exp_out[1:, :, :] = data[:-1, :, :] + assert_almost_equal(out.dataobj, exp_out) + # Then 'nearest' mode + out = resample_from_to(img, trans1_img, mode='nearest') + exp_out[0, :, :] = exp_out[1, :, :] + assert_almost_equal(out.dataobj, exp_out) + # Test order + trans_p_25_aff = from_matvec(np.diag([-4, 5, 6]), [1, 0, 0]) + trans_p_25_img = Nifti1Image(data, trans_p_25_aff) + # Suprising to me, but all points outside are set to 0, even with NN + out = resample_from_to(img, trans_p_25_img, order=0) + exp_out = np.zeros_like(data) + exp_out[1:, :, :] = data[1, :, :] + assert_almost_equal(out.dataobj, exp_out) + out = resample_from_to(img, trans_p_25_img) + exp_out = spnd.affine_transform(data, [1, 1, 1], [-0.25, 0, 0], order=3) + assert_almost_equal(out.dataobj, exp_out) + # Test cval + out = resample_from_to(img, trans_img, cval=99) + exp_out = np.zeros_like(data) + 99 + exp_out[1:, 1:, 1:] = data[:-1, :-1, :-1] + assert_almost_equal(out.dataobj, exp_out) + # Out class + out = resample_from_to(img, trans_img) + assert_equal(out.__class__, Nifti1Image) + # By default, type of from_img makes no difference + n1_img = Nifti2Image(data, affine) + out = resample_from_to(n1_img, trans_img) + assert_equal(out.__class__, Nifti1Image) + # Passed as keyword arg + out = resample_from_to(img, trans_img, out_class=Nifti2Image) + assert_equal(out.__class__, Nifti2Image) + # If keyword arg is None, use type of from_img + out = resample_from_to(n1_img, trans_img, out_class=None) + assert_equal(out.__class__, Nifti2Image) + # to_img type irrelevant in all cases + n1_trans_img = Nifti2Image(data, trans_aff) + out = resample_from_to(img, n1_trans_img, out_class=None) + assert_equal(out.__class__, Nifti1Image) + # From 2D to 3D, error, the fixed affine is not invertible + img_2d = Nifti1Image(data[:, :, 0], affine) + assert_raises(AffineError, resample_from_to, img_2d, img) + # 3D to 2D, we don't need to invert the fixed matrix + out = resample_from_to(img, img_2d) + assert_array_equal(out.dataobj, data[:, :, 0]) + # Same for tuple as to_img imput + out = resample_from_to(img, (img_2d.shape, img_2d.affine)) + assert_array_equal(out.dataobj, data[:, :, 0]) + # 4D input and output also OK + data_4d = np.arange(24 * 5).reshape((2, 3, 4, 5)) + img_4d = Nifti1Image(data_4d, affine) + out = resample_from_to(img_4d, img_4d) + assert_almost_equal(data_4d, out.dataobj) + assert_array_equal(img_4d.affine, out.affine) + # Errors trying to match 3D to 4D + assert_raises(ValueError, resample_from_to, img_4d, img) + assert_raises(ValueError, resample_from_to, img, img_4d) + + +@needs_scipy +def test_resample_to_output(): + # Test routine to sample iamges to output space + # Image aligned to output axes - no-op + data = np.arange(24).reshape((2, 3, 4)) + img = Nifti1Image(data, np.eye(4)) + # Check default resampling + img2 = resample_to_output(img) + assert_array_equal(img2.shape, (2, 3, 4)) + assert_array_equal(img2.affine, np.eye(4)) + assert_array_equal(img2.dataobj, data) + # Check resampling with different voxel size specifications + for vox_sizes in (None, 1, [1, 1, 1]): + img2 = resample_to_output(img, vox_sizes) + assert_array_equal(img2.shape, (2, 3, 4)) + assert_array_equal(img2.affine, np.eye(4)) + assert_array_equal(img2.dataobj, data) + img2 = resample_to_output(img, vox_sizes) + # Check 2D works + img_2d = Nifti1Image(data[0], np.eye(4)) + img3 = resample_to_output(img_2d) + assert_array_equal(img3.shape, (3, 4, 1)) + assert_array_equal(img3.affine, np.eye(4)) + assert_array_equal(img3.dataobj, data[0][..., None]) + # Even 1D + img_1d = Nifti1Image(data[0, 0], np.eye(4)) + img3 = resample_to_output(img_1d) + assert_array_equal(img3.shape, (4, 1, 1)) + assert_array_equal(img3.affine, np.eye(4)) + assert_array_equal(img3.dataobj, data[0, 0][..., None, None]) + # But 4D does not + img_4d = Nifti1Image(data.reshape(2, 3, 2, 2), np.eye(4)) + assert_raises(ValueError, resample_to_output, img_4d) + # Run vox2vox_out tests, checking output shape, coordinate transform + for in_shape, in_aff, vox, out_shape, out_aff in get_outspace_params(): + # Allow for expansion of image shape from < 3D + in_n_dim = len(in_shape) + if in_n_dim < 3: + in_shape = in_shape + (1,) * (3 - in_n_dim) + if not vox is None: + vox = vox + (1,) * (3 - in_n_dim) + assert len(out_shape) == in_n_dim + out_shape = out_shape + (1,) * (3 - in_n_dim) + img = Nifti1Image(np.ones(in_shape), in_aff) + out_img = resample_to_output(img, vox) + assert_all_in(in_shape, in_aff, out_img.shape, out_img.affine) + assert_equal(out_img.shape, out_shape) + assert_almost_equal(out_img.affine, out_aff) + # Check data is as expected with some transforms + # Flip first axis + out_img = resample_to_output(Nifti1Image(data, np.diag([-1, 1, 1, 1]))) + assert_array_equal(out_img.dataobj, np.flipud(data)) + # Subsample voxels + out_img = resample_to_output(Nifti1Image(data, np.diag([4, 5, 6, 1]))) + exp_out = spnd.affine_transform(data, + [1/4, 1/5, 1/6], + output_shape = (5, 11, 19)) + assert_array_equal(out_img.dataobj, exp_out) + # Unsubsample with voxel sizes + out_img = resample_to_output(Nifti1Image(data, np.diag([4, 5, 6, 1])), + [4, 5, 6]) + assert_array_equal(out_img.dataobj, data) + # A rotation to test nearest, order, cval + rot_3 = from_matvec(euler2mat(np.pi / 4), [0, 0, 0]) + rot_3_img = Nifti1Image(data, rot_3) + out_img = resample_to_output(rot_3_img) + exp_shape = (4, 4, 4) + assert_equal(out_img.shape, exp_shape) + exp_aff = np.array([[1, 0, 0, -2 * np.cos(np.pi / 4)], + [0, 1, 0, 0], + [0, 0, 1, 0], + [0, 0, 0, 1]]) + assert_almost_equal(out_img.affine, exp_aff) + rzs, trans = to_matvec(np.dot(npl.inv(rot_3), exp_aff)) + exp_out = spnd.affine_transform(data, rzs, trans, exp_shape) + assert_almost_equal(out_img.dataobj, exp_out) + # Order + assert_almost_equal( + resample_to_output(rot_3_img, order=0).dataobj, + spnd.affine_transform(data, rzs, trans, exp_shape, order=0)) + # Cval + assert_almost_equal( + resample_to_output(rot_3_img, cval=99).dataobj, + spnd.affine_transform(data, rzs, trans, exp_shape, cval=99)) + # Mode + assert_almost_equal( + resample_to_output(rot_3_img, mode='nearest').dataobj, + spnd.affine_transform(data, rzs, trans, exp_shape, mode='nearest')) + # out_class + img_ni1 = Nifti2Image(data, np.eye(4)) + img_ni2 = Nifti2Image(data, np.eye(4)) + # Default is Nifti1Image + assert_equal( + resample_to_output(img_ni2).__class__, + Nifti1Image) + # Can be overriden + assert_equal( + resample_to_output(img_ni1, out_class=Nifti2Image).__class__, + Nifti2Image) + # None specifies out_class from input + assert_equal( + resample_to_output(img_ni2, out_class=None).__class__, + Nifti2Image) + + +@needs_scipy +def test_smooth_image(): + # Test image smoothing + data = np.arange(24).reshape((2, 3, 4)) + aff = np.diag([-4, 5, 6, 1]) + img = Nifti1Image(data, aff) + # Zero smoothing is no-op + out_img = smooth_image(img, 0) + assert_array_equal(out_img.affine, img.affine) + assert_array_equal(out_img.shape, img.shape) + assert_array_equal(out_img.dataobj, data) + # Isotropic smoothing + sd = fwhm2sigma(np.true_divide(8, [4, 5, 6])) + exp_out = spnd.gaussian_filter(data, sd, mode='nearest') + assert_array_equal(smooth_image(img, 8).dataobj, exp_out) + assert_array_equal(smooth_image(img, [8, 8, 8]).dataobj, exp_out) + assert_raises(ValueError, smooth_image, img, [8, 8]) + # Not isotropic + mixed_sd = fwhm2sigma(np.true_divide([8, 7, 6], [4, 5, 6])) + exp_out = spnd.gaussian_filter(data, mixed_sd, mode='nearest') + assert_array_equal(smooth_image(img, [8, 7, 6]).dataobj, exp_out) + # In 2D + img_2d = Nifti1Image(data[0], aff) + exp_out = spnd.gaussian_filter(data[0], sd[:2], mode='nearest') + assert_array_equal(smooth_image(img_2d, 8).dataobj, exp_out) + assert_array_equal(smooth_image(img_2d, [8, 8]).dataobj, exp_out) + assert_raises(ValueError, smooth_image, img_2d, [8, 8, 8]) + # Isotropic in 4D has zero for last dimension in scalar case + data_4d = np.arange(24 * 5).reshape((2, 3, 4, 5)) + img_4d = Nifti1Image(data_4d, aff) + exp_out = spnd.gaussian_filter(data_4d, list(sd) + [0], mode='nearest') + assert_array_equal(smooth_image(img_4d, 8).dataobj, exp_out) + # But raises error for vector case + assert_raises(ValueError, smooth_image, img_4d, [8, 8, 8]) + # mode, cval + exp_out = spnd.gaussian_filter(data, sd, mode='constant') + assert_array_equal(smooth_image(img, 8, mode='constant').dataobj, exp_out) + exp_out = spnd.gaussian_filter(data, sd, mode='constant', cval=99) + assert_array_equal(smooth_image(img, 8, mode='constant', cval=99).dataobj, + exp_out) + # out_class + img_ni1 = Nifti2Image(data, np.eye(4)) + img_ni2 = Nifti2Image(data, np.eye(4)) + # Default is Nifti1Image + assert_equal( + smooth_image(img_ni2, 0).__class__, + Nifti1Image) + # Can be overriden + assert_equal( + smooth_image(img_ni1, 0, out_class=Nifti2Image).__class__, + Nifti2Image) + # None specifies out_class from input + assert_equal( + smooth_image(img_ni2, 0, out_class=None).__class__, + Nifti2Image) + + +def assert_spm_resampling_close(from_img, our_resampled, spm_resampled): + """ Assert our resampling is close to SPM's, allowing for edge effects + """ + # To allow for differences in the way SPM and scipy.ndimage handle off-edge + # interpolation, mask out voxels off edge + to_img_shape = spm_resampled.shape + to_img_affine = spm_resampled.affine + to_vox_coords = np.indices(to_img_shape).transpose((1, 2, 3, 0)) + # Coordinates of to_img mapped to from_img + to_to_from = npl.inv(from_img.affine).dot(to_img_affine) + resamp_coords = apply_affine(to_to_from, to_vox_coords) + # Places where SPM may not return default value but scipy.ndimage will (SPM + # does not return zeros <0.05 from image edges). + # See: https://github.com/nipy/nibabel/pull/255#issuecomment-186774173 + outside_vol = np.any((resamp_coords < 0) | + (np.subtract(resamp_coords, from_img.shape) > -1), + axis=-1) + spm_res = np.where(outside_vol, np.nan, np.array(spm_resampled.dataobj)) + assert_allclose_safely(our_resampled.dataobj, spm_res) + assert_almost_equal(our_resampled.affine, spm_resampled.affine, 5) + + +@needs_scipy +def test_against_spm_resample(): + # Test resampling against images resampled with SPM12 + # anatomical.nii has a diagonal -2, 2 2 affine; + # functional.nii has a diagonal -4, 4 4 affine; + # These are a bit boring, so first add some rotations and translations to + # the anatomical image affine, and then resample to the first volume in the + # functional, and compare to the same thing in SPM. + # See ``make_moved_anat.py`` script in this directory for input to SPM. + anat = nib.load(pjoin(DATA_DIR, 'anatomical.nii')) + func = nib.load(pjoin(DATA_DIR, 'functional.nii')) + some_rotations = euler2mat(0.1, 0.2, 0.3) + extra_affine = from_matvec(some_rotations, [3, 4, 5]) + moved_anat = nib.Nifti1Image(anat.get_data().astype(float), + extra_affine.dot(anat.affine), + anat.header) + one_func = nib.Nifti1Image(func.dataobj[..., 0], + func.affine, + func.header) + moved2func = resample_from_to(moved_anat, one_func, order=1, cval=np.nan) + spm_moved = nib.load(pjoin(DATA_DIR, 'resampled_anat_moved.nii')) + assert_spm_resampling_close(moved_anat, moved2func, spm_moved) + # Next we resample the rotated anatomical image to output space, and compare + # to the same operation done with SPM (our own version of 'reorient.m' by + # John Ashburner). + moved2output = resample_to_output(moved_anat, 4, order=1, cval=np.nan) + spm2output = nib.load(pjoin(DATA_DIR, 'reoriented_anat_moved.nii')) + assert_spm_resampling_close(moved_anat, moved2output, spm2output); diff --git a/nibabel/tests/test_spaces.py b/nibabel/tests/test_spaces.py index 667059333c..2520d32225 100644 --- a/nibabel/tests/test_spaces.py +++ b/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)