Merge pull request #1 from matthew-brett/setting-mgh-affines

Setting mgh affines

Author: larsoner

nibabel/freesurfer/mghformat.py

@@ -229,6 +229,9 @@ def get_affine(self):
         M[0:3, 3] = pxyz_0.T
         return M
 
+    # For compatibility with nifti (multiple affines)
+    get_best_affine = get_affine
+
     def get_vox2ras(self):
         '''return the get_affine()
         '''
@@ -554,20 +557,25 @@ def update_header(self):
         hdr = self._header
         if not self._data is None:
             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
+        if self._affine is None:
+            return
+        if np.allclose(self._affine, hdr.get_best_affine()):
+            return
+        # for more information, go through save_mgh.m in FreeSurfer dist
+        MdcD = self._affine[:3, :3]
+        delta = np.sqrt(np.sum(MdcD * MdcD, axis=0))
+        Mdc = MdcD / np.tile(delta, (3, 1))
+        Pcrs_c = np.array([0, 0, 0, 1], dtype=np.float)
+        Pcrs_c[:3] = np.array([self._data.shape[0], self._data.shape[1],
+                                self._data.shape[2]], dtype=np.float) / 2.0
+        Pxyz_c = np.dot(self._affine, Pcrs_c)
+
+        hdr['delta'][:] = delta
+        hdr['Mdc'][:, :] = Mdc.T
+        hdr['Pxyz_c'][:] = Pxyz_c[:3]
 
-        if not self._affine is None:
-            # for more information, go through save_mgh.m in FreeSurfer dist
-            MdcD = self._affine[:3, :3]
-            delta = np.sqrt(np.sum(MdcD * MdcD, axis=0))
-            Mdc = MdcD / np.tile(delta, (3, 1))
-            Pcrs_c = np.array([0, 0, 0, 1], dtype=np.float)
-            Pcrs_c[:3] = np.array([self._data.shape[0], self._data.shape[1],
-                                   self._data.shape[2]], dtype=np.float) / 2.0
-            Pxyz_c = np.dot(self._affine, Pcrs_c)
-
-            hdr['delta'][:] = delta
-            hdr['Mdc'][:, :] = Mdc.T
-            hdr['Pxyz_c'][:] = Pxyz_c[:3]
 
 load = MGHImage.load
 save = MGHImage.instance_to_filename

nibabel/freesurfer/tests/test_mghformat.py

@@ -8,11 +8,17 @@
 ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
 '''Tests for mghformat reading writing'''
 from __future__ import with_statement
+
 import os
+from io import StringIO, BytesIO
+
 import numpy as np
 from .. import load, save, MGHImage
 from ..mghformat import MGHError
 from ...tmpdirs import InTemporaryDirectory
+from ...py3k import unicode
+from ...fileholders import FileHolder
+
 from ...testing import data_path
 from numpy.testing import assert_equal, assert_array_equal, \
     assert_array_almost_equal, assert_almost_equal, assert_raises
@@ -47,7 +53,7 @@ def test_read_mgh():
     assert_equal(h['goodRASFlag'], 1)
     assert_array_equal(h['dims'], [3, 4, 5, 2])
     assert_array_almost_equal(h['mrparms'], [2.0, 0.0, 0.0, 0.0])
-    assert_array_almost_equal(h.get_zooms(), zooms)
+    assert_array_almost_equal(h.get_zooms(), 1)
     assert_array_almost_equal(h.get_vox2ras(), v2r)
     assert_array_almost_equal(h.get_vox2ras_tkr(), v2rtkr)
 
@@ -150,3 +156,45 @@ def test_filename_exts():
             img_back = load(fname)
             assert_array_equal(img_back.get_data(), v)
             del img_back
+
+
+def _mgh_rt(img, fobj):
+    file_map = {'image': FileHolder(fileobj=fobj)}
+    img.to_file_map(file_map)
+    return MGHImage.from_file_map(file_map)
+
+
+def test_header_updating():
+    # Don't update the header information if the affine doesn't change.
+    # Luckily the test.mgz dataset had a bad set of cosine vectors, so these
+    # will be changed if the affine gets updated
+    mgz_path = os.path.join(data_path, 'test.mgz')
+    mgz = load(mgz_path)
+    hdr = mgz.get_header()
+    # Test against mri_info output
+    exp_aff = np.loadtxt(StringIO(unicode("""
+    1.0000   2.0000   3.0000   -13.0000
+    2.0000   3.0000   1.0000   -11.5000
+    3.0000   1.0000   2.0000   -11.5000
+    0.0000   0.0000   0.0000     1.0000""")))
+    assert_almost_equal(mgz.get_affine(), exp_aff, 6)
+    assert_almost_equal(hdr.get_affine(), exp_aff, 6)
+    # Test that initial wonky header elements have not changed
+    assert_equal(hdr['delta'], 1)
+    assert_almost_equal(hdr['Mdc'], exp_aff[:3, :3].T)
+    # Save, reload, same thing
+    img_fobj = BytesIO()
+    mgz2 = _mgh_rt(mgz, img_fobj)
+    hdr2 = mgz2.get_header()
+    assert_almost_equal(hdr2.get_affine(), exp_aff, 6)
+    assert_equal(hdr2['delta'], 1)
+    # Change affine, change underlying header info
+    exp_aff_d = exp_aff.copy()
+    exp_aff_d[0, -1] = -14
+    # This will (probably) become part of the official API
+    mgz2._affine[:] = exp_aff_d
+    mgz2.update_header()
+    assert_almost_equal(hdr2.get_affine(), exp_aff_d, 6)
+    RZS = exp_aff_d[:3, :3]
+    assert_almost_equal(hdr2['delta'], np.sqrt(np.sum(RZS ** 2, axis=0)))
+    assert_almost_equal(hdr2['Mdc'], (RZS / hdr2['delta']).T)