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_processing.py b/nibabel/tests/test_processing.py index 20eb218054..60ad23eb13 100644 --- a/nibabel/tests/test_processing.py +++ b/nibabel/tests/test_processing.py @@ -10,19 +10,22 @@ """ from __future__ import division, print_function +from os.path import dirname, join as pjoin + import numpy as np import numpy.linalg as npl -from ..optpkg import optional_package +from nibabel.optpkg import optional_package spnd, have_scipy, _ = optional_package('scipy.ndimage') -from ..processing import (sigma2fwhm, fwhm2sigma, adapt_affine, - resample_from_to, resample_to_output, smooth_image) -from ..nifti1 import Nifti1Image -from ..nifti2 import Nifti2Image -from ..orientations import flip_axis, inv_ornt_aff -from ..affines import AffineError, from_matvec, to_matvec -from ..eulerangles import euler2mat +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) @@ -31,10 +34,12 @@ from nose.tools import (assert_true, assert_false, assert_raises, assert_equal, assert_not_equal) -from .test_spaces import assert_all_in, get_outspace_params +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 @@ -339,3 +344,54 @@ def test_smooth_image(): 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, rtol=1e-4, atol=1e-5) + assert_almost_equal(our_resampled.affine, spm_resampled.affine, 5) + + +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);