Merge pull request #101 from matthew-brett/freesurfer-hack

Freesurfer large vector hack

Add hack to allow Freesurfer files that put large vector sizes in glmin.  Alexandre Gramfort spotted files using this hack.

Author: matthew-brett

nibabel/analyze.py

@@ -626,7 +626,7 @@ def set_data_shape(self, shape):
         # Check that dimensions fit
         if not np.all(dims[1:ndims+1] == shape):
             raise HeaderDataError('shape %s does not fit in dim datatype' %
-                                   shape)
+                                   (shape,))
         self._structarr['pixdim'][ndims+1:] = 1.0
 
     def get_base_affine(self):
@@ -974,7 +974,7 @@ def update_header(self):
         # We need to update the header if the data shape has changed.  It's a
         # bit difficult to change the data shape using the standard API, but
         # maybe it happened
-        if not self._data is None:
+        if not self._data is None and hdr.get_data_shape() != self._data.shape:
             hdr.set_data_shape(self._data.shape)
         # If the affine is not None, and it is different from the main affine in
         # the header, update the heaader

nibabel/nifti1.py

@@ -8,6 +8,8 @@
 ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
 ''' Header reading / writing functions for nifti1 image format
 '''
+import warnings
+
 import numpy as np
 import numpy.linalg as npl
 
@@ -620,6 +622,75 @@ def default_structarr(klass, endianness=None):
             hdr_data['vox_offset'] = 0
         return hdr_data
 
+    def get_data_shape(self):
+        ''' Get shape of data
+
+        Examples
+        --------
+        >>> hdr = Nifti1Header()
+        >>> hdr.get_data_shape()
+        (0,)
+        >>> hdr.set_data_shape((1,2,3))
+        >>> hdr.get_data_shape()
+        (1, 2, 3)
+
+        Expanding number of dimensions gets default zooms
+
+        >>> hdr.get_zooms()
+        (1.0, 1.0, 1.0)
+
+        Notes
+        -----
+        Allows for freesurfer hack for large vectors described in
+        https://github.com/nipy/nibabel/issues/100 and
+        http://code.google.com/p/fieldtrip/source/browse/trunk/external/freesurfer/save_nifti.m?spec=svn5022&r=5022#77
+        '''
+        shape = super(Nifti1Header, self).get_data_shape()
+        # Apply freesurfer hack for vector
+        if shape != (-1, 1, 1): # Normal case
+            return shape
+        vec_len = int(self._structarr['glmin'])
+        if vec_len == 0:
+            raise HeaderDataError('-1 in dim[1] but 0 in glmin; inconsistent '
+                                  'freesurfer type header?')
+        return (vec_len, 1, 1)
+
+    def set_data_shape(self, shape):
+        ''' Set shape of data
+
+        If ``ndims == len(shape)`` then we set zooms for dimensions higher than
+        ``ndims`` to 1.0
+
+        Parameters
+        ----------
+        shape : sequence
+           sequence of integers specifying data array shape
+
+        Notes
+        -----
+        Applies freesurfer hack for large vectors described in
+        https://github.com/nipy/nibabel/issues/100 and
+        http://code.google.com/p/fieldtrip/source/browse/trunk/external/freesurfer/save_nifti.m?spec=svn5022&r=5022#77
+        '''
+        # Apply freesurfer hack for vector
+        hdr = self._structarr
+        shape = tuple(shape)
+        if (len(shape) == 3 and shape[1:] == (1, 1) and
+            shape[0] > np.iinfo(hdr['dim'].dtype.base).max): # Freesurfer case
+            try:
+                hdr['glmin'] = shape[0]
+            except OverflowError:
+                overflow = True
+            else:
+                overflow = hdr['glmin'] != shape[0]
+            if overflow:
+                raise HeaderDataError('shape[0] %s does not fit in glmax datatype' %
+                                      shape[0])
+            warnings.warn('Using large vector Freesurfer hack; header will '
+                          'not be compatible with SPM or FSL', stacklevel=2)
+            shape = (-1, 1, 1)
+        super(Nifti1Header, self).set_data_shape(shape)
+
     def get_qform_quaternion(self):
         ''' Compute quaternion from b, c, d of quaternion
 

nibabel/tests/test_analyze.py

@@ -234,6 +234,11 @@ def test_shapes(self):
         assert_equal(hdr.get_data_shape(), shape)
         shape = (mx+1,)
         assert_raises(HeaderDataError, hdr.set_data_shape, shape)
+        # Lists or tuples or arrays will work for setting shape
+        shape = (2, 3, 4)
+        for constructor in (list, tuple, np.array):
+            hdr.set_data_shape(constructor(shape))
+            assert_equal(hdr.get_data_shape(), shape)
 
     def test_read_write_data(self):
         # Check reading and writing of data

nibabel/tests/test_nifti1.py

@@ -169,6 +169,49 @@ def test_nifti_log_checks(self):
         assert_equal(message, 'sform_code -1 not valid; '
                            'setting to 0')
 
+    def test_freesurfer_hack(self):
+        # For large vector images, Freesurfer appears to set dim[1] to -1 and
+        # then use glmin for the vector length (an i4)
+        HC = self.header_class
+        # The standard case
+        hdr = HC()
+        hdr.set_data_shape((2, 3, 4))
+        assert_equal(hdr.get_data_shape(), (2, 3, 4))
+        assert_equal(hdr['glmin'], 0)
+        # Just left of the freesurfer case
+        dim_type = hdr.template_dtype['dim'].base
+        glmin = hdr.template_dtype['glmin'].base
+        too_big = int(np.iinfo(dim_type).max) + 1
+        hdr.set_data_shape((too_big-1, 1, 1))
+        assert_equal(hdr.get_data_shape(), (too_big-1, 1, 1))
+        # The freesurfer case
+        hdr.set_data_shape((too_big, 1, 1))
+        assert_equal(hdr.get_data_shape(), (too_big, 1, 1))
+        assert_array_equal(hdr['dim'][:4], [3, -1, 1, 1])
+        assert_equal(hdr['glmin'], too_big)
+        # This only works for the case of a 3D with -1, 1, 1
+        assert_raises(HeaderDataError, hdr.set_data_shape, (too_big,))
+        assert_raises(HeaderDataError, hdr.set_data_shape, (too_big,1))
+        assert_raises(HeaderDataError, hdr.set_data_shape, (too_big,1,2))
+        assert_raises(HeaderDataError, hdr.set_data_shape, (too_big,2,1))
+        assert_raises(HeaderDataError, hdr.set_data_shape, (1, too_big))
+        assert_raises(HeaderDataError, hdr.set_data_shape, (1, too_big, 1))
+        assert_raises(HeaderDataError, hdr.set_data_shape, (1, 1, too_big))
+        # Outside range of glmin raises error
+        far_too_big = int(np.iinfo(glmin).max) + 1
+        hdr.set_data_shape((far_too_big-1, 1, 1))
+        assert_equal(hdr.get_data_shape(), (far_too_big-1, 1, 1))
+        assert_raises(HeaderDataError, hdr.set_data_shape, (far_too_big,1,1))
+        # glmin of zero raises error (implausible vector length)
+        hdr.set_data_shape((-1,1,1))
+        hdr['glmin'] = 0
+        assert_raises(HeaderDataError, hdr.get_data_shape)
+        # Lists or tuples or arrays will work for setting shape
+        for shape in ((too_big-1, 1, 1), (too_big, 1, 1)):
+            for constructor in (list, tuple, np.array):
+                hdr.set_data_shape(constructor(shape))
+                assert_equal(hdr.get_data_shape(), shape)
+
 
 class TestNifti1SingleHeader(TestNifti1PairHeader):