From 95e6b75794ca3026fca7d4c3a4b8f1a1da0f666f Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Fri, 22 Mar 2013 15:14:47 -0700 Subject: [PATCH 1/9] ENH: Adding vox2ras_tkr --- nibabel/freesurfer/mghformat.py | 11 +++++++++++ nibabel/freesurfer/tests/test_mghformat.py | 6 +++++- 2 files changed, 16 insertions(+), 1 deletion(-) diff --git a/nibabel/freesurfer/mghformat.py b/nibabel/freesurfer/mghformat.py index 645a995427..f502ec9d5b 100644 --- a/nibabel/freesurfer/mghformat.py +++ b/nibabel/freesurfer/mghformat.py @@ -234,6 +234,17 @@ def get_vox2ras(self): ''' 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 + ''' + pcrs_c = self._header_data['dims'][:3] / 2.0 + v2rtkr = np.array([[-1, 0, 0, pcrs_c[0]], + [0, 0, 1, -pcrs_c[2]], + [0, -1, 0, pcrs_c[1]], + [0, 0, 0, 1]], dtype=np.float32) + return v2rtkr + def get_ras2vox(self): '''return the inverse get_affine() ''' diff --git a/nibabel/freesurfer/tests/test_mghformat.py b/nibabel/freesurfer/tests/test_mghformat.py index 903408dbd5..77d6956f64 100644 --- a/nibabel/freesurfer/tests/test_mghformat.py +++ b/nibabel/freesurfer/tests/test_mghformat.py @@ -17,9 +17,12 @@ 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, 1.5,], [0, 0, 1, -2.5], + [0, -1, 0, 2], [0, 0, 0, 1]], dtype=np.float32) def test_read_mgh(): @@ -39,6 +42,7 @@ def test_read_mgh(): 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_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() From d5f8c079c1c7fa749472a958ea4ba1947804a1b2 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Wed, 27 Mar 2013 11:17:10 -0700 Subject: [PATCH 2/9] FIX(?): Add zooms test --- nibabel/freesurfer/mghformat.py | 12 +++++++----- nibabel/freesurfer/tests/test_mghformat.py | 5 ++++- 2 files changed, 11 insertions(+), 6 deletions(-) diff --git a/nibabel/freesurfer/mghformat.py b/nibabel/freesurfer/mghformat.py index f502ec9d5b..16d16515b8 100644 --- a/nibabel/freesurfer/mghformat.py +++ b/nibabel/freesurfer/mghformat.py @@ -238,10 +238,11 @@ def get_vox2ras_tkr(self): ''' Get the vox2ras-tkr transform. See "Torig" here: http://surfer.nmr.mgh.harvard.edu/fswiki/CoordinateSystems ''' - pcrs_c = self._header_data['dims'][:3] / 2.0 - v2rtkr = np.array([[-1, 0, 0, pcrs_c[0]], - [0, 0, 1, -pcrs_c[2]], - [0, -1, 0, pcrs_c[1]], + 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 @@ -557,7 +558,8 @@ def update_header(self): 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)) + delta = hdr['delta'] + #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], diff --git a/nibabel/freesurfer/tests/test_mghformat.py b/nibabel/freesurfer/tests/test_mghformat.py index 77d6956f64..04cdf3976d 100644 --- a/nibabel/freesurfer/tests/test_mghformat.py +++ b/nibabel/freesurfer/tests/test_mghformat.py @@ -21,8 +21,10 @@ 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, 1.5,], [0, 0, 1, -2.5], +v2rtkr = np.array([[-1, 0, 0, 1.5], [0, 0, 1, -2.5], [0, -1, 0, 2], [0, 0, 0, 1]], dtype=np.float32) +# sample voxel sizes (mri_info --res) +zooms = np.array([1.000, 1.000, 1.000]) def test_read_mgh(): @@ -41,6 +43,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_vox2ras(), v2r) assert_array_almost_equal(h.get_vox2ras_tkr(), v2rtkr) From 39ab09f8e698ff53d5c62ef20c57025ae1abd3a9 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Wed, 27 Mar 2013 14:12:29 -0700 Subject: [PATCH 3/9] FIX: Deal with borked test dataset --- nibabel/freesurfer/mghformat.py | 3 +-- nibabel/freesurfer/tests/test_mghformat.py | 14 +++++++++----- 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/nibabel/freesurfer/mghformat.py b/nibabel/freesurfer/mghformat.py index 16d16515b8..0cd789fc3e 100644 --- a/nibabel/freesurfer/mghformat.py +++ b/nibabel/freesurfer/mghformat.py @@ -558,8 +558,7 @@ def update_header(self): if not self._affine is None: # for more information, go through save_mgh.m in FreeSurfer dist MdcD = self._affine[:3, :3] - delta = hdr['delta'] - #delta = np.sqrt(np.sum(MdcD * MdcD, axis=0)) + 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], diff --git a/nibabel/freesurfer/tests/test_mghformat.py b/nibabel/freesurfer/tests/test_mghformat.py index 04cdf3976d..b020ddd06c 100644 --- a/nibabel/freesurfer/tests/test_mghformat.py +++ b/nibabel/freesurfer/tests/test_mghformat.py @@ -20,11 +20,15 @@ # 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, 1.5], [0, 0, 1, -2.5], - [0, -1, 0, 2], [0, 0, 0, 1]], dtype=np.float32) -# sample voxel sizes (mri_info --res) -zooms = np.array([1.000, 1.000, 1.000]) +# sample voxel sizes (mri_info --res gives different result because the +# direction cosines in our test file are borked) +zooms = np.array([3.7416575, 3.7416575, 3.7416575]) +# sample voxel to ras - tkr matrix (mri_info --vox2ras-tkr also gives +# different results because of the direction cosines, again) +v2rtkr = np.array([[-3.7416575, 0.0, 0.0, 5.61248636], + [0.0, 0.0, 3.7416575, -9.3541441], + [0.0, -3.7416575, 0.0, 7.48331499], + [0.0, 0.0, 0.0, 1.0]], dtype=np.float32) def test_read_mgh(): From 57f108933f2e7535d5515b73dc7a0a305a556c45 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Wed, 27 Mar 2013 15:47:02 -0700 Subject: [PATCH 4/9] FIX: Match MRIRead.m --- nibabel/freesurfer/mghformat.py | 4 ++-- nibabel/freesurfer/tests/test_mghformat.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/nibabel/freesurfer/mghformat.py b/nibabel/freesurfer/mghformat.py index 0cd789fc3e..1bade28105 100644 --- a/nibabel/freesurfer/mghformat.py +++ b/nibabel/freesurfer/mghformat.py @@ -240,9 +240,9 @@ def get_vox2ras_tkr(self): ''' 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]], + v2rtkr = np.array([[-ds[1], 0, 0, ns[1]], [0, 0, ds[2], -ns[2]], - [0, -ds[1], 0, ns[1]], + [0, -ds[0], 0, ns[0]], [0, 0, 0, 1]], dtype=np.float32) return v2rtkr diff --git a/nibabel/freesurfer/tests/test_mghformat.py b/nibabel/freesurfer/tests/test_mghformat.py index b020ddd06c..9b3ec26860 100644 --- a/nibabel/freesurfer/tests/test_mghformat.py +++ b/nibabel/freesurfer/tests/test_mghformat.py @@ -25,9 +25,9 @@ zooms = np.array([3.7416575, 3.7416575, 3.7416575]) # sample voxel to ras - tkr matrix (mri_info --vox2ras-tkr also gives # different results because of the direction cosines, again) -v2rtkr = np.array([[-3.7416575, 0.0, 0.0, 5.61248636], +v2rtkr = np.array([[-3.7416575, 0.0, 0.0, 7.48331499], [0.0, 0.0, 3.7416575, -9.3541441], - [0.0, -3.7416575, 0.0, 7.48331499], + [0.0, -3.7416575, 0.0, 5.61248636], [0.0, 0.0, 0.0, 1.0]], dtype=np.float32) From 6724c781fbeec63908d02a810eeb74d4f0491dfa Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Wed, 27 Mar 2013 15:52:32 -0700 Subject: [PATCH 5/9] FIX?: Order --- nibabel/freesurfer/mghformat.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/nibabel/freesurfer/mghformat.py b/nibabel/freesurfer/mghformat.py index 1bade28105..0cd789fc3e 100644 --- a/nibabel/freesurfer/mghformat.py +++ b/nibabel/freesurfer/mghformat.py @@ -240,9 +240,9 @@ def get_vox2ras_tkr(self): ''' ds = np.array(self._header_data['delta']) ns = (np.array(self._header_data['dims'][:3]) * ds) / 2.0 - v2rtkr = np.array([[-ds[1], 0, 0, ns[1]], + v2rtkr = np.array([[-ds[0], 0, 0, ns[0]], [0, 0, ds[2], -ns[2]], - [0, -ds[0], 0, ns[0]], + [0, -ds[1], 0, ns[1]], [0, 0, 0, 1]], dtype=np.float32) return v2rtkr From 4e0d3ead0985c00b08fda2490a0263292dfa1155 Mon Sep 17 00:00:00 2001 From: Matthew Brett Date: Wed, 27 Mar 2013 16:34:44 -0700 Subject: [PATCH 6/9] RF: make header update only if changed affine We were having trouble with the header updating on load of an image, and changing the underlying header information. https://github.com/nipy/nibabel/pull/164 Use the Analyze trick to only update the header if the affine has changed. --- nibabel/freesurfer/mghformat.py | 34 +++++++++------ nibabel/freesurfer/tests/test_mghformat.py | 48 ++++++++++++++++++++++ 2 files changed, 69 insertions(+), 13 deletions(-) diff --git a/nibabel/freesurfer/mghformat.py b/nibabel/freesurfer/mghformat.py index 0cd789fc3e..2047dc83b7 100644 --- a/nibabel/freesurfer/mghformat.py +++ b/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 diff --git a/nibabel/freesurfer/tests/test_mghformat.py b/nibabel/freesurfer/tests/test_mghformat.py index 9b3ec26860..fb7aafd857 100644 --- a/nibabel/freesurfer/tests/test_mghformat.py +++ b/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 @@ -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) From 13e87c0b0e7ea60548931d12e81c6e4f1b2dca1d Mon Sep 17 00:00:00 2001 From: Matthew Brett Date: Wed, 27 Mar 2013 16:38:56 -0700 Subject: [PATCH 7/9] BF: with new save affine algorithm, zooms=1 No we are not resetting the deltas on loading the image, we should get zooms of 1 from the test image. --- nibabel/freesurfer/tests/test_mghformat.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nibabel/freesurfer/tests/test_mghformat.py b/nibabel/freesurfer/tests/test_mghformat.py index fb7aafd857..18baa3f450 100644 --- a/nibabel/freesurfer/tests/test_mghformat.py +++ b/nibabel/freesurfer/tests/test_mghformat.py @@ -53,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) From d5588bfbc99c4475b5400e444750e59bced6a737 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Wed, 27 Mar 2013 16:51:00 -0700 Subject: [PATCH 8/9] FIX: Reverting to mri_info test compat --- nibabel/freesurfer/tests/test_mghformat.py | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/nibabel/freesurfer/tests/test_mghformat.py b/nibabel/freesurfer/tests/test_mghformat.py index 18baa3f450..9e735d777b 100644 --- a/nibabel/freesurfer/tests/test_mghformat.py +++ b/nibabel/freesurfer/tests/test_mghformat.py @@ -26,14 +26,10 @@ # 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 sizes (mri_info --res gives different result because the -# direction cosines in our test file are borked) -zooms = np.array([3.7416575, 3.7416575, 3.7416575]) -# sample voxel to ras - tkr matrix (mri_info --vox2ras-tkr also gives -# different results because of the direction cosines, again) -v2rtkr = np.array([[-3.7416575, 0.0, 0.0, 7.48331499], - [0.0, 0.0, 3.7416575, -9.3541441], - [0.0, -3.7416575, 0.0, 5.61248636], +# 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) From 2b0ecbe35ac96becc58f4d8cc11b2f0bc6405df6 Mon Sep 17 00:00:00 2001 From: Matthew Brett Date: Wed, 27 Mar 2013 16:56:46 -0700 Subject: [PATCH 9/9] TST: add test for cosine vector order We weren't previously testing the transpose / not-transpose of the cosine vectors because the test dataset had a symmetical 3x3 in the affine. --- nibabel/freesurfer/tests/test_mghformat.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/nibabel/freesurfer/tests/test_mghformat.py b/nibabel/freesurfer/tests/test_mghformat.py index 9e735d777b..ed46f60053 100644 --- a/nibabel/freesurfer/tests/test_mghformat.py +++ b/nibabel/freesurfer/tests/test_mghformat.py @@ -194,3 +194,19 @@ def test_header_updating(): 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)