BF: FreeSurfer nifti surfaces can have >3 dimensions

Author: effigies

nibabel/nifti1.py

@@ -704,16 +704,16 @@ def get_data_shape(self):
         https://code.google.com/p/fieldtrip/source/browse/trunk/external/freesurfer/save_nifti.m?r=8776#50
         '''
         shape = super(Nifti1Header, self).get_data_shape()
-        # Apply freesurfer hack for vector
-        if shape == (-1, 1, 1):
+        # Apply freesurfer hack for large vectors
+        if shape[:3] == (-1, 1, 1):
             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)
+            return (vec_len, 1, 1) + shape[3:]
         # Apply freesurfer hack for ico7 surface
-        elif shape == (27307, 1, 6):
-            return (163842, 1, 1)
+        elif shape[:3] == (27307, 1, 6):
+            return (163842, 1, 1) + shape[3:]
         else:  # Normal case
             return shape
 
@@ -723,6 +723,11 @@ def set_data_shape(self, shape):
         If ``ndims == len(shape)`` then we set zooms for dimensions higher than
         ``ndims`` to 1.0
 
+        Nifti1 images can have up to seven dimensions. For Nifti surface files,
+        the first dimension is assumed to correspond to vertices/nodes on a
+        surface, and dimensions two and three are constrained to have depth of
+        1. Dimensions 4-7 are constrained only by type bounds.
+
         Parameters
         ----------
         shape : sequence
@@ -744,10 +749,10 @@ def set_data_shape(self, shape):
         shape = tuple(shape)
 
         # Apply freesurfer hack for ico7 surface
-        if shape == (163842, 1, 1):
-            shape = (27307, 1, 6)
-        # Apply freesurfer hack for vector
-        elif (len(shape) == 3 and shape[1:] == (1, 1) and
+        if shape[:3] == (163842, 1, 1):
+            shape = (27307, 1, 6) + shape[3:]
+        # Apply freesurfer hack for large vectors
+        elif (len(shape) >= 3 and shape[1:3] == (1, 1) and
                 shape[0] > np.iinfo(hdr['dim'].dtype.base).max):
             try:
                 hdr['glmin'] = shape[0]
@@ -760,7 +765,7 @@ def set_data_shape(self, shape):
                                       '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)
+            shape = (-1, 1, 1) + shape[3:]
         super(Nifti1Header, self).set_data_shape(shape)
 
     def get_qform_quaternion(self):

nibabel/tests/test_nifti1.py

@@ -261,19 +261,28 @@ def test_freesurfer_large_vector_hack(self):
         hdr.set_data_shape((too_big-1, 1, 1))
         assert_equal(hdr.get_data_shape(), (too_big-1, 1, 1))
         # The freesurfer case
+        full_shape = (too_big, 1, 1, 1, 1, 1, 1)
+        for dim in range(3, 8):
+            expected_dim = np.array([dim, -1, 1, 1, 1, 1, 1, 1])
+            with suppress_warnings():
+                hdr.set_data_shape(full_shape[:dim])
+            assert_equal(hdr.get_data_shape(), full_shape[:dim])
+            assert_array_equal(hdr['dim'], expected_dim)
+            assert_equal(hdr['glmin'], too_big)
+        # Allow the fourth dimension to vary
         with suppress_warnings():
-            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
+            hdr.set_data_shape((too_big, 1, 1, 4))
+        assert_equal(hdr.get_data_shape(), (too_big, 1, 1, 4))
+        assert_array_equal(hdr['dim'][:5], np.array([4, -1, 1, 1, 4]))
+        # This only works when the first 3 dimensions are -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))
+        assert_raises(HeaderDataError, hdr.set_data_shape, (1, 1, 1, too_big))
         # Outside range of glmin raises error
         far_too_big = int(np.iinfo(glmin).max) + 1
         with suppress_warnings():
@@ -295,9 +304,22 @@ def test_freesurfer_large_vector_hack(self):
     def test_freesurfer_ico7_hack(self):
         HC = self.header_class
         hdr = HC()
+        full_shape = (163842, 1, 1, 1, 1, 1, 1)
         # Test that using ico7 shape automatically uses factored dimensions
-        hdr.set_data_shape((163842, 1, 1))
-        assert_array_equal(hdr._structarr['dim'][1:4], np.array([27307, 1, 6]))
+        for dim in range(3, 8):
+            expected_dim = np.array([dim, 27307, 1, 6, 1, 1, 1, 1])
+            hdr.set_data_shape(full_shape[:dim])
+            assert_equal(hdr.get_data_shape(), full_shape[:dim])
+            assert_array_equal(hdr._structarr['dim'], expected_dim)
+        # Only works on dimensions >= 3
+        assert_raises(HeaderDataError, hdr.set_data_shape, full_shape[:1])
+        assert_raises(HeaderDataError, hdr.set_data_shape, full_shape[:2])
+        # Bad shapes
+        assert_raises(HeaderDataError, hdr.set_data_shape, (163842, 2, 1))
+        assert_raises(HeaderDataError, hdr.set_data_shape, (163842, 1, 2))
+        assert_raises(HeaderDataError, hdr.set_data_shape, (1, 163842, 1))
+        assert_raises(HeaderDataError, hdr.set_data_shape, (1, 1, 163842))
+        assert_raises(HeaderDataError, hdr.set_data_shape, (1, 1, 1, 163842))
         # Test consistency of data in .mgh and mri_convert produced .nii
         nitest_path = os.path.join(get_nibabel_data(), 'nitest-freesurfer')
         mgh = mghload(os.path.join(nitest_path, 'fsaverage', 'surf',