diff --git a/nibabel/freesurfer/mghformat.py b/nibabel/freesurfer/mghformat.py index 645a995427..2047dc83b7 100644 --- a/nibabel/freesurfer/mghformat.py +++ b/nibabel/freesurfer/mghformat.py @@ -229,11 +229,26 @@ 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() ''' return self.get_affine() + def get_vox2ras_tkr(self): + ''' Get the vox2ras-tkr transform. See "Torig" here: + http://surfer.nmr.mgh.harvard.edu/fswiki/CoordinateSystems + ''' + ds = np.array(self._header_data['delta']) + ns = (np.array(self._header_data['dims'][:3]) * ds) / 2.0 + v2rtkr = np.array([[-ds[0], 0, 0, ns[0]], + [0, 0, ds[2], -ns[2]], + [0, -ds[1], 0, ns[1]], + [0, 0, 0, 1]], dtype=np.float32) + return v2rtkr + def get_ras2vox(self): '''return the inverse get_affine() ''' @@ -542,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 diff --git a/nibabel/freesurfer/tests/test_mghformat.py b/nibabel/freesurfer/tests/test_mghformat.py index 903408dbd5..ed46f60053 100644 --- a/nibabel/freesurfer/tests/test_mghformat.py +++ b/nibabel/freesurfer/tests/test_mghformat.py @@ -8,18 +8,29 @@ ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## '''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 -# sample voxel to ras matrix +# sample voxel to ras matrix (mri_info --vox2ras) v2r = np.array([[1, 2, 3, -13], [2, 3, 1, -11.5], [3, 1, 2, -11.5], [0, 0, 0, 1]], dtype=np.float32) +# sample voxel to ras - tkr matrix (mri_info --vox2ras-tkr) +v2rtkr = np.array([[-1.0, 0.0, 0.0, 1.5], + [0.0, 0.0, 1.0, -2.5], + [0.0, -1.0, 0.0, 2.0], + [0.0, 0.0, 0.0, 1.0]], dtype=np.float32) def test_read_mgh(): @@ -38,7 +49,9 @@ 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(), 1) assert_array_almost_equal(h.get_vox2ras(), v2r) + assert_array_almost_equal(h.get_vox2ras_tkr(), v2rtkr) # data. will be different for your own mri_volsynth invocation v = mgz.get_data() @@ -139,3 +152,61 @@ 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) + + +def test_cosine_order(): + # Test we are interpreting the cosine order right + data = np.arange(60).reshape((3, 4, 5)).astype(np.int32) + aff = np.diag([2., 3, 4, 1]) + aff[0] = [2, 1, 0, 10] + img = MGHImage(data, aff) + assert_almost_equal(img.get_affine(), aff, 6) + img_fobj = BytesIO() + img2 = _mgh_rt(img, img_fobj) + hdr2 = img2.get_header() + RZS = aff[:3, :3] + zooms = np.sqrt(np.sum(RZS ** 2, axis=0)) + assert_almost_equal(hdr2['Mdc'], (RZS / zooms).T) + assert_almost_equal(hdr2['delta'], zooms)