Merge pull request #332 from effigies/freesurfer4d

MRG: FreeSurfer nifti surfaces can have >3 dimensions

Quick workaround for https://mail.python.org/pipermail/neuroimaging/2015-July/000132.html.

This uses the FreeSurfer hacks on any shapes beginning with (-1, 1, 1) or (27307, 1, 6), which seems fairly safe, but we can restrict down to 4-dimensional, or possibly only allow (-1, 1, 1, 1) or (27307, 1, 6, 1) in the 4D case.

Author: matthew-brett

nibabel/nifti1.py

@@ -693,27 +693,23 @@ def get_data_shape(self):
 
         Notes
         -----
-        Allows for freesurfer hack for large vectors described in
-        https://github.com/nipy/nibabel/issues/100 and
-        https://code.google.com/p/fieldtrip/source/browse/trunk/external/freesurfer/save_nifti.m?spec=svn5022&r=5022#77
+        Applies freesurfer hack for large vectors described in `issue 100`_ and
+        `save_nifti.m <save77_>`_.
 
         Allows for freesurfer hack for 7th order icosahedron surface described
-        in
-        https://github.com/nipy/nibabel/issues/309
-        https://code.google.com/p/fieldtrip/source/browse/trunk/external/freesurfer/load_nifti.m?r=8776#86
-        https://code.google.com/p/fieldtrip/source/browse/trunk/external/freesurfer/save_nifti.m?r=8776#50
+        in `issue 309`_, load_nifti.m_, and `save_nifti.m <save50_>`_.
         '''
         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,31 +719,56 @@ 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 FreeSurfer-variant
+        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
            sequence of integers specifying data array shape
 
         Notes
         -----
-        Applies freesurfer hack for large vectors described in
-        https://github.com/nipy/nibabel/issues/100 and
-        https://code.google.com/p/fieldtrip/source/browse/trunk/external/freesurfer/save_nifti.m?spec=svn5022&r=5022#77
+        Applies freesurfer hack for large vectors described in `issue 100`_ and
+        `save_nifti.m <save77_>`_.
 
         Allows for freesurfer hack for 7th order icosahedron surface described
-        in
-        https://github.com/nipy/nibabel/issues/309
-        https://code.google.com/p/fieldtrip/source/browse/trunk/external/freesurfer/load_nifti.m?r=8776#86
-        https://code.google.com/p/fieldtrip/source/browse/trunk/external/freesurfer/save_nifti.m?r=8776#50
+        in `issue 309`_, load_nifti.m_, and `save_nifti.m <save50_>`_.
+
+        The Nifti1 `standard header`_ allows for the following "point set"
+        definition of a surface, not currently implemented in nibabel.
+
+        ::
+
+          To signify that the vector value at each voxel is really a
+          spatial coordinate (e.g., the vertices or nodes of a surface mesh):
+            - dataset must have a 5th dimension
+            - intent_code must be NIFTI_INTENT_POINTSET
+            - dim[0] = 5
+            - dim[1] = number of points
+            - dim[2] = dim[3] = dim[4] = 1
+            - dim[5] must be the dimensionality of space (e.g., 3 => 3D space).
+            - intent_name may describe the object these points come from
+              (e.g., "pial", "gray/white" , "EEG", "MEG").
+
+        .. _issue 100: https://github.com/nipy/nibabel/issues/100
+        .. _issue 309: https://github.com/nipy/nibabel/issues/309
+        .. _save77: https://code.google.com/p/fieldtrip/source/browse/trunk/external/freesurfer/save_nifti.m?spec=svn8776&r=8776#77
+        .. _save50: https://code.google.com/p/fieldtrip/source/browse/trunk/external/freesurfer/save_nifti.m?spec=svn8776&r=8776#50
+        .. _load_nifti.m: https://code.google.com/p/fieldtrip/source/browse/trunk/external/freesurfer/load_nifti.m?spec=svn8776&r=8776#86
+        .. _standard header: http://nifti.nimh.nih.gov/pub/dist/src/niftilib/nifti1.h
         '''
         hdr = self._structarr
         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 +781,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,29 @@ 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):
+            # First element in 'dim' field is number of dimensions
+            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 +305,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',