Merge pull request #464 from matthew-brett/minc-defense

MRG: check for MINC using processing routines

The processing routines assume that image have spatial axes before time etc
axes, but this may well not be true for MINC images.

Barf on MINC images > 3D for now. Later maybe we can investigate the image
axis names to work out which axes are spatial.

Author: matthew-brett

nibabel/imageclasses.py

@@ -104,3 +104,29 @@ def __getitem__(self, *args, **kwargs):
     ('mgz', '.mgz'),
     ('par', '.par'),
 ))
+
+# Image classes known to require spatial axes to be first in index ordering.
+# When adding an image class, consider whether the new class should be listed
+# here.
+KNOWN_SPATIAL_FIRST = (Nifti1Pair, Nifti1Image, Nifti2Pair, Nifti2Image,
+                       Spm2AnalyzeImage, Spm99AnalyzeImage, AnalyzeImage,
+                       MGHImage, PARRECImage)
+
+
+def spatial_axes_first(img):
+    """ True if spatial image axes for `img` always preceed other axes
+
+    Parameters
+    ----------
+    img : object
+        Image object implementing at least ``shape`` attribute.
+
+    Returns
+    -------
+    spatial_axes_first : bool
+        True if image only has spatial axes (number of axes < 4) or image type
+        known to have spatial axes preceeding other axes.
+    """
+    if len(img.shape) < 4:
+        return True
+    return type(img) in KNOWN_SPATIAL_FIRST

nibabel/processing.py

@@ -25,6 +25,7 @@
 from .affines import AffineError, to_matvec, from_matvec, append_diag
 from .spaces import vox2out_vox
 from .nifti1 import Nifti1Image
+from .imageclasses import spatial_axes_first
 
 SIGMA2FWHM = np.sqrt(8 * np.log(2))
 
@@ -124,9 +125,9 @@ def resample_from_to(from_img,
     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.
+        Object having attributes ``dataobj``, ``affine``, ``header`` and
+        ``shape``. 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
@@ -153,6 +154,10 @@ def resample_from_to(from_img,
         resampling `from_img` into axes aligned to the output space of
         ``from_img.affine``
     """
+    # This check requires `shape` attribute of image
+    if not spatial_axes_first(from_img):
+        raise ValueError('Cannot predict position of spatial axes for Image '
+                         'type ' + str(type(from_img)))
     try:
         to_shape, to_affine = to_vox_map.shape, to_vox_map.affine
     except AttributeError:
@@ -248,9 +253,9 @@ def smooth_image(img,
     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.
+        Object having attributes ``dataobj``, ``affine``, ``header`` and
+        ``shape``. 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
@@ -280,6 +285,10 @@ def smooth_image(img,
         Image of instance specified by `out_class`, containing data output from
         smoothing `img` data by given FWHM kernel.
     """
+    # This check requires `shape` attribute of image
+    if not spatial_axes_first(img):
+        raise ValueError('Cannot predict position of spatial axes for Image '
+                         'type ' + str(type(img)))
     if out_class is None:
         out_class = img.__class__
     n_dim = len(img.shape)

nibabel/tests/test_imageclasses.py

@@ -0,0 +1,48 @@
+""" Testing imageclasses module
+"""
+
+from os.path import dirname, join as pjoin
+
+import numpy as np
+
+from nibabel.optpkg import optional_package
+
+import nibabel as nib
+from nibabel.analyze import AnalyzeImage
+from nibabel.nifti1 import Nifti1Image
+from nibabel.nifti2 import Nifti2Image
+
+from nibabel.imageclasses import spatial_axes_first
+
+from nose.tools import (assert_true, assert_false)
+
+DATA_DIR = pjoin(dirname(__file__), 'data')
+
+have_h5py = optional_package('h5py')[1]
+
+MINC_3DS = ('minc1_1_scale.mnc',)
+MINC_4DS = ('minc1_4d.mnc',)
+if have_h5py:
+    MINC_3DS = MINC_3DS + ('minc2_1_scale.mnc',)
+    MINC_4DS = MINC_4DS + ('minc2_4d.mnc',)
+
+
+def test_spatial_axes_first():
+    # Function tests is spatial axes are first three axes in image
+    # Always True for Nifti and friends
+    affine = np.eye(4)
+    for shape in ((2, 3), (4, 3, 2), (5, 4, 1, 2), (2, 3, 5, 2, 1)):
+        for img_class in (AnalyzeImage, Nifti1Image, Nifti2Image):
+            data = np.zeros(shape)
+            img = img_class(data, affine)
+            assert_true(spatial_axes_first(img))
+    # True for MINC images < 4D
+    for fname in MINC_3DS:
+        img = nib.load(pjoin(DATA_DIR, fname))
+        assert_true(len(img.shape) == 3)
+        assert_true(spatial_axes_first(img))
+    # False for MINC images < 4D
+    for fname in MINC_4DS:
+        img = nib.load(pjoin(DATA_DIR, fname))
+        assert_true(len(img.shape) == 4)
+        assert_false(spatial_axes_first(img))

nibabel/tests/test_processing.py

@@ -24,7 +24,8 @@
 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.affines import (AffineError, from_matvec, to_matvec, apply_affine,
+                             voxel_sizes)
 from nibabel.eulerangles import euler2mat
 
 from numpy.testing import (assert_almost_equal,
@@ -41,6 +42,15 @@
 
 DATA_DIR = pjoin(dirname(__file__), 'data')
 
+# 3D MINC work correctly with processing, but not 4D MINC
+from .test_imageclasses import MINC_3DS, MINC_4DS
+
+# Filenames of other images that should work correctly with processing
+OTHER_IMGS = ('anatomical.nii', 'functional.nii',
+              'example4d.nii.gz', 'example_nifti2.nii.gz',
+              'phantom_EPI_asc_CLEAR_2_1.PAR')
+
+
 def test_sigma2fwhm():
     # Test from constant
     assert_almost_equal(sigma2fwhm(1), 2.3548200)
@@ -346,6 +356,26 @@ def test_smooth_image():
         Nifti2Image)
 
 
+@needs_scipy
+def test_spatial_axes_check():
+    for fname in MINC_3DS + OTHER_IMGS:
+        img = nib.load(pjoin(DATA_DIR, fname))
+        s_img = smooth_image(img, 0)
+        assert_array_equal(img.dataobj, s_img.dataobj)
+        out = resample_from_to(img, img, mode='nearest')
+        assert_almost_equal(img.dataobj, out.dataobj)
+        if len(img.shape) > 3:
+            continue
+        # Resample to output does not raise an error
+        out = resample_to_output(img, voxel_sizes(img.affine))
+    for fname in MINC_4DS:
+        img = nib.load(pjoin(DATA_DIR, fname))
+        assert_raises(ValueError, smooth_image, img, 0)
+        assert_raises(ValueError, resample_from_to, img, img, mode='nearest')
+        assert_raises(ValueError,
+                      resample_to_output, img, voxel_sizes(img.affine))
+
+
 def assert_spm_resampling_close(from_img, our_resampled, spm_resampled):
     """ Assert our resampling is close to SPM's, allowing for edge effects
     """