TST: add tests against SPM resampling

Build some test images by resampling with SPM, and test our resampling
against SPM's resampling.

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