Skip to content

ENH: Adding vox2ras_tkr #164

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 11 commits into from
Mar 28, 2013
Merged

ENH: Adding vox2ras_tkr #164

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
95e6b75
ENH: Adding vox2ras_tkr
larsoner Mar 22, 2013
d5f8c07
FIX(?): Add zooms test
larsoner Mar 27, 2013
39ab09f
FIX: Deal with borked test dataset
larsoner Mar 27, 2013
57f1089
FIX: Match MRIRead.m
larsoner Mar 27, 2013
6724c78
FIX?: Order
larsoner Mar 27, 2013
4e0d3ea
RF: make header update only if changed affine
matthew-brett Mar 27, 2013
13e87c0
BF: with new save affine algorithm, zooms=1
matthew-brett Mar 27, 2013
1ab91f0
Merge pull request #1 from matthew-brett/setting-mgh-affines
larsoner Mar 27, 2013
d5588bf
FIX: Reverting to mri_info test compat
larsoner Mar 27, 2013
2b0ecbe
TST: add test for cosine vector order
matthew-brett Mar 27, 2013
06a0a23
Merge pull request #2 from matthew-brett/mgh-cosine-vector-order
larsoner Mar 28, 2013
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
15 changes: 14 additions & 1 deletion nibabel/freesurfer/mghformat.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -234,6 +234,18 @@ 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
'''
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)
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These are the lines I'm unclear about (not that there are so many to choose from). I expected the 1's in the upper left 3x3 of this matrix to be given by the values in get_zooms(), but using those values makes the test script fail (since the zooms are not unity, but mri_info --vox2ras-tkr on the test .mgz file has only +/-1's in the upper left 3x3...

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Using mri_info to get the resolution in each direction gives:

$ mri_info --res test.mgz
1.000 1.000 1.000 2.000

However, in nibabel getting the deltas I get:

>>> nib.load('test.mgz').get_header().get_zooms()
(3.7416575, 3.7416575, 3.7416575)

So, clearly, get_zooms() is doing something else that I mis-interpreted. I'm not sure where to get the equivalent information in nibabel, though...

return v2rtkr

def get_ras2vox(self):
'''return the inverse get_affine()
'''
Expand Down Expand Up @@ -546,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))
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This was the potentially offending line. Not sure why it's necessary to "re-compute" the deltas if we have them stored, but this re-computation was breaking vox2ras-tkr equivalence.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@satra any idea why the delta values were being re-calculated here? Why is it necessary? They replace the original (correct?) values of hdr['delta'] a couple lines later...

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],
Expand Down
9 changes: 8 additions & 1 deletion nibabel/freesurfer/tests/test_mghformat.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,9 +17,14 @@
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)
# sample voxel sizes (mri_info --res)
zooms = np.array([1.000, 1.000, 1.000])


def test_read_mgh():
Expand All @@ -38,7 +43,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(), zooms)
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note that this test now passes, whereas before it didn't.

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()
Expand Down