BF+TST - SPM analyze had 111 voxel offset in mats

The .mat files for SPM analyze images were not being interpreted right
because the affines therein have 111 voxel origins not 000 as for nifti.

Author: matthew-brett

nibabel/spm99analyze.py

@@ -261,7 +261,6 @@ def from_file_map(klass, file_map):
                               'using first')
                 mat = mat[:, :, 0]
             ret._affine = mat
-            return ret
         elif 'M' in mats: # the 'M' matrix does not include flips
             hdr = ret._header
             if hdr.default_x_flip:
@@ -270,6 +269,10 @@ def from_file_map(klass, file_map):
                 ret._affine = mats['M']
         else:
             raise ValueError('mat file found but no "mat" or "M" in it')
+        # Adjust for matlab 1,1,1 voxel origin
+        to_111 = np.eye(4)
+        to_111[:3,3] = 1
+        ret._affine = np.dot(ret._affine, to_111)
         return ret
 
     def to_file_map(self, file_map=None):
@@ -295,6 +298,11 @@ def to_file_map(self, file_map=None):
             M = np.dot(np.diag([-1, 1, 1, 1]), mat)
         else:
             M = mat
+        # Adjust for matlab 1,1,1 voxel origin
+        from_111 = np.eye(4)
+        from_111[:3,3] = -1
+        M = np.dot(M, from_111)
+        mat = np.dot(mat, from_111)
         # use matlab 4 format to allow gzipped write without error
         mfobj = file_map['mat'].get_prepare_fileobj(mode='wb')
         sio.savemat(mfobj, {'M': M, 'mat': mat}, format='4')

nibabel/tests/test_spm99analyze.py

@@ -105,6 +105,59 @@ class TestSpm99AnalyzeImage(test_analyze.TestAnalyzeImage):
         test_analyze.TestAnalyzeImage.test_data_hdr_cache
     ))
 
+    @scipy_skip
+    def test_mat_read(self):
+        # Test mat file reading and writing for the SPM analyze types
+        img_klass = self.image_class
+        arr = np.arange(24).reshape((2,3,4))
+        aff = np.diag([2,3,4,1]) # no LR flip in affine
+        img = img_klass(arr, aff)
+        fm = img.file_map
+        for key, value in fm.items():
+            value.fileobj = BytesIO()
+        # Test round trip
+        img.to_file_map()
+        r_img = img_klass.from_file_map(fm)
+        assert_array_equal(r_img.get_data(), arr)
+        assert_array_equal(r_img.get_affine(), aff)
+        # mat files are for matlab and have 111 voxel origins.  We need to
+        # adjust for that, when loading and saving.  Check for signs of that in
+        # the saved mat file
+        mat_fileobj = img.file_map['mat'].fileobj
+        from scipy.io import loadmat, savemat
+        mat_fileobj.seek(0)
+        mats = loadmat(mat_fileobj)
+        assert_true('M' in mats and 'mat' in mats)
+        from_111 = np.eye(4)
+        from_111[:3,3] = -1
+        to_111 = np.eye(4)
+        to_111[:3,3] = 1
+        assert_array_equal(mats['mat'], np.dot(aff, from_111))
+        # The M matrix does not include flips, so if we only
+        # have the M matrix in the mat file, and we have default flipping, the
+        # mat resulting should have a flip.  The 'mat' matrix does include flips
+        # and so should be unaffected by the flipping.  If both are present we
+        # prefer the the 'mat' matrix.
+        assert_true(img.get_header().default_x_flip) # check the default
+        flipper = np.diag([-1,1,1,1])
+        assert_array_equal(mats['M'], np.dot(aff, np.dot(flipper, from_111)))
+        mat_fileobj.seek(0)
+        savemat(mat_fileobj, dict(M=np.diag([3,4,5,1]), mat=np.diag([6,7,8,1])))
+        # Check we are preferring the 'mat' matrix
+        r_img = img_klass.from_file_map(fm)
+        assert_array_equal(r_img.get_data(), arr)
+        assert_array_equal(r_img.get_affine(),
+                           np.dot(np.diag([6,7,8,1]), to_111))
+        # But will use M if present
+        mat_fileobj.seek(0)
+        mat_fileobj.truncate(0)
+        savemat(mat_fileobj, dict(M=np.diag([3,4,5,1])))
+        r_img = img_klass.from_file_map(fm)
+        assert_array_equal(r_img.get_data(), arr)
+        assert_array_equal(r_img.get_affine(),
+                           np.dot(np.diag([3,4,5,1]), np.dot(flipper, to_111)))
+
+
 
 def test_origin_affine():
     hdr = Spm99AnalyzeHeader()