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.

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

merged 11 commits into from
Mar 28, 2013

Conversation

larsoner
Copy link
Contributor

This adds the support requested in #163.

Note that I am not 100% sure about this formulation. Based on the powerpoint linked from this page (http://surfer.nmr.mgh.harvard.edu/fswiki/CoordinateSystems) I thought that the [-1, -1, 1] triplet that is used to fill in the upper left corner of the matrix would need to be [-delta[0], -delta[2], delta[1]], but that wasn't what the output of mri_info --vox2ras-tkr gave me. Can someone confirm before merging?

@agramfort
Copy link
Contributor

tried with 2 MRIs. Results with mri_info match.

+1 for merge

@larsoner
Copy link
Contributor Author

@agramfort do you know what's up with the difference between the freesurfer definition of vox2ras-tkr (which seems to use dC, dR, and dS) and the output of mri_info --vox2ras-tkr? I can't quite make sense of it, even after reading through the powerpoint.

@agramfort
Copy link
Contributor

Nope sorry

@agramfort
Copy link
Contributor

can this be merged or @Eric89GXL you think it's not ready?

@matthew-brett any opinion?

@larsoner
Copy link
Contributor Author

@agramfort I'm just concerned about the mismatch between what Freesurfer documented, and the results from mri_info --vox2ras-tkr that I've replicated in this PR. It would be good to get someone more familiar with the transforms to double-check it before merging, if possible.

@agramfort
Copy link
Contributor

ping @satra any idea?

or you can ask doug on the FS mailing list...

@larsoner
Copy link
Contributor Author

I asked about the transforms at one point, and I was directed to the powerpoint. Part of the issue is that I'm not 100% sure what the dR, dC, dS terms are in the powerpoint. I assumed this was the size of the voxels in each dimension, and assumed this information was stored in img.get_header().get_zooms(), but that must not be the case since then there should be non-1's in the upper 3x3 of the vox2ras-tkr...

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)
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...

@larsoner
Copy link
Contributor Author

I also shot an email to Doug from the Freesurfer list to ask about how to get the "xsize", "ysize", and "zsize" parameters (which can be obtained from mri_info --res), since I'll probably need to add support for that in this PR, too.

@agramfort
Copy link
Contributor

can't you obtain the voxels sizes for the get_affine(). You should.

@larsoner
Copy link
Contributor Author

@agramfort the get_affine() command returns the same thing as mri_info --vox2ras command. Here are --vox2ras and --vox2ras-tkr for one subject:

larsoner@bunk:~/Documents/structurals/AKCLEE_101/mri$ mri_info --vox2ras-tkr orig.mgz 
  -1.00000    0.00000    0.00000  128.00000 
   0.00000    0.00000    1.00000 -128.00000 
   0.00000   -1.00000    0.00000  128.00000 
   0.00000    0.00000    0.00000    1.00000 
larsoner@bunk:~/Documents/structurals/AKCLEE_101/mri$ mri_info --vox2ras orig.mgz 
  -1.00000    0.00000   -0.00000  126.47404 
  -0.00000    0.00000    1.00000 -116.68102 
   0.00000   -1.00000   -0.00000  121.04381 
   0.00000    0.00000    0.00000    1.00000

And for the test dataset:

larsoner@bunk:~/custombuilds/nibabel/nibabel/tests/data$ mri_info test.mgz --vox2ras-tkr
  -1.00000    0.00000    0.00000    1.50000 
   0.00000    0.00000    1.00000   -2.50000 
   0.00000   -1.00000    0.00000    2.00000 
   0.00000    0.00000    0.00000    1.00000 
larsoner@bunk:~/custombuilds/nibabel/nibabel/tests/data$ mri_info test.mgz --vox2ras
   1.00000    2.00000    3.00000  -13.00000 
   2.00000    3.00000    1.00000  -11.50000 
   3.00000    1.00000    2.00000  -11.50000 
   0.00000    0.00000    0.00000    1.00000

I don't see an obvious mapping between the two, and looking at the code and powerpoint, I'm not sure what it would be, either...

@satra
Copy link
Member

satra commented Mar 27, 2013

i believe mri_info --vox2ras-tkr will always return the observed rotation matrix and the fov/2 for translation. is the size of test.mgz 3, 4, 5?

@larsoner
Copy link
Contributor Author

@satra yes, but the issue is whether or not it is just a pure rotation (all 1's), or a rotation with scaling (non-1 values). In the Freesurfer coordinate system .ppt, it says these values must be the voxel sizes dR/dC/dS, which happen to be 1 for these datasets, but I don't think have to in general...? Thus I'm trying to figure out where to get the csize/rsize/ssize-equivalent information from the header.

@satra
Copy link
Member

satra commented Mar 27, 2013

(devpype)[satra@warpspeed temp]$ mri_info out.mgz 
Volume information for out.mgz
          type: MGH
    dimensions: 64 x 64 x 35
   voxel sizes: 3.4375, 3.4375, 4.0000
          type: FLOAT (3)
           fov: 220.000
           dof: 0
        xstart: -110.0, xend: 110.0
        ystart: -110.0, yend: 110.0
        zstart: -70.0, zend: 70.0
            TR: 3000.00 msec, TE: 0.00 msec, TI: 0.00 msec, flip angle: 0.00 degrees
       nframes: 1
       PhEncDir: UNKNOWN
ras xform present
    xform info: x_r =  -0.9973, y_r =   0.0637, z_r =   0.0379, c_r =     1.7908
              : x_a =  -0.0720, y_a =  -0.9532, z_a =  -0.2937, c_a =     1.5082
              : x_s =   0.0174, y_s =  -0.2957, z_s =   0.9551, c_s =     5.3327

talairach xfm : 
Orientation   : LPS
Primary Slice Direction: axial

voxel to ras transform:
               -3.4281   0.2189   0.1514   101.8337
               -0.2476  -3.2765  -1.1749   134.8401
                0.0597  -1.0163   3.8206   -30.9163
                0.0000   0.0000   0.0000     1.0000

voxel-to-ras determinant 47.2656

ras to voxel transform:
               -0.2901  -0.0210   0.0051    32.5242
                0.0185  -0.2773  -0.0860    32.8437
                0.0095  -0.0734   0.2388    16.3204
                0.0000   0.0000   0.0000     1.0000
(devpype)[satra@warpspeed temp]$ mri_info --vox2ras-tkr out.mgz 
  -3.43750    0.00000    0.00000  110.00000 
   0.00000    0.00000    4.00000  -70.00000 
   0.00000   -3.43750    0.00000  110.00000 
   0.00000    0.00000    0.00000    1.00000 

and then with nibabel:

In [1]: import nibabel as nib

In [2]: hdr = nib.load('out.mgz').get_header()

In [5]: hdr.get_zooms()
Out[5]: (3.4375, 3.4374998, 4.0)

In [6]: print hdr.get_affine()
[[ -3.42805386e+00   2.18914479e-01   1.51402012e-01   1.01833700e+02]
 [ -2.47560844e-01  -3.27652526e+00  -1.17492497e+00   1.34840126e+02]
 [  5.97159155e-02  -1.01629841e+00   3.82055354e+00  -3.09163036e+01]
 [  0.00000000e+00   0.00000000e+00   0.00000000e+00   1.00000000e+00]]

In [7]: size = nib.load('out.mgz').shape

In [10]: import numpy as np

In [11]: (np.array(size) * np.array(hdr.get_zooms()))/2
Out[11]: array([ 110.        ,  109.99999237,   70.        ])

so we have all the components needed:

Out[5]: (3.4375, 3.4374998, 4.0)
Out[11]: array([ 110.        ,  109.99999237,   70.        ])

to create:

(devpype)[satra@warpspeed temp]$ mri_info --vox2ras-tkr out.mgz 
  -3.43750    0.00000    0.00000  110.00000 
   0.00000    0.00000    4.00000  -70.00000 
   0.00000   -3.43750    0.00000  110.00000 
   0.00000    0.00000    0.00000    1.00000 

@larsoner
Copy link
Contributor Author

Here is the same process for the test.mgz file (which is <100kB if you want to try it):

larsoner@bunk:~/custombuilds/nibabel/nibabel/tests/data$ mri_info test.mgz 
Volume information for test.mgz
          type: MGH
    dimensions: 3 x 4 x 5 x 2
   voxel sizes: 1.0000, 1.0000, 1.0000
          type: FLOAT (3)
           fov: 5.000
           dof: 0
        xstart: -1.5, xend: 1.5
        ystart: -2.0, yend: 2.0
        zstart: -2.5, zend: 2.5
            TR: 2.00 msec, TE: 0.00 msec, TI: 0.00 msec, flip angle: 0.00 degrees
       nframes: 2
       PhEncDir: UNKNOWN
ras xform present
    xform info: x_r =   1.0000, y_r =   2.0000, z_r =   3.0000, c_r =     0.0000
              : x_a =   2.0000, y_a =   3.0000, z_a =   1.0000, c_a =     0.0000
              : x_s =   3.0000, y_s =   1.0000, z_s =   2.0000, c_s =     0.0000

talairach xfm : 
Orientation   : SAR
Primary Slice Direction: sagittal

voxel to ras transform:
                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

voxel-to-ras determinant -18

ras to voxel transform:
               -0.2778   0.0556   0.3889     1.5000
                0.0556   0.3889  -0.2778     2.0000
                0.3889  -0.2778   0.0556     2.5000
                0.0000   0.0000   0.0000     1.0000

and then with nibabel:

>>> import nibabel as nib
>>> hdr = nib.load('test.mgz').get_header()
>>> hdr.get_zooms()
(3.7416575, 3.7416575, 3.7416575)
>>> print hdr.get_affine()
[[  1.           2.           3.00000024 -13.00000041]
 [  2.           3.00000024   1.         -11.50000034]
 [  3.00000024   1.           2.         -11.50000029]
 [  0.           0.           0.           1.        ]]
>>> size = nib.load('test.mgz').shape
>>> import numpy as np
>>> (np.array(size)[:3] * np.array(hdr.get_zooms()))/2
array([ 5.61248624,  7.48331499,  9.35414374])

These would be the corresponding components:

(3.7416575, 3.7416575, 3.7416575)
array([ 5.61248624,  7.48331499,  9.35414374])

But they don't create:

mri_info --vox2ras-tkr test.mgz 
  -1.00000    0.00000    0.00000    1.50000 
   0.00000    0.00000    1.00000   -2.50000 
   0.00000   -1.00000    0.00000    2.00000 
   0.00000    0.00000    0.00000    1.00000

Thus my confusion...

@larsoner
Copy link
Contributor Author

I can understand the last column mismatch, since I expected those values just to be np.array(size)[:3] / 2. This would give a correct value for my example, but an incorrect one for yours (!), even though I thought I was following the freesurfer coordinates .ppt. The upper-left 3x3 mismatch seems to be driven by mri_info saying that the voxel sizes are 1.0000, 1.0000, 1.0000, while nibabel says they are (3.7416575, 3.7416575, 3.7416575), which I can't reconcile.

@larsoner
Copy link
Contributor Author

Actually, both can thus be explained by the mismatch between hdr.get_zooms() and mri_info --csize/rsize/ssize. If those matched [(1, 1, 1) I guess is correct?], then they would match in both examples. So I guess nibabel is mis-reading the zooms, then...?

@larsoner
Copy link
Contributor Author

Either that or zooms a.k.a. deltas are not the same thing as csize/rsize/ssize, but I don't know enough about this stuff to say either way.

@agramfort
Copy link
Contributor

sorry I cannot really help ... but I see you moving forward :)

@larsoner
Copy link
Contributor Author

I consider this most recent revelation somewhat of a lateral move, actually :)

@agramfort
Copy link
Contributor

at least you're moving :)

@matthew-brett
Copy link
Member

I see from the nibabel code that it's getting the zooms from the the Spacing X, Y, Z fields:

http://surfer.nmr.mgh.harvard.edu/fswiki/FsTutorial/MghFormat

Is that incorrect?

This is what I get for the header info from nibabel:

import nibabel.freesurfer as nif
img = nif.load('nibabel/tests/data/test.mgz')
print(img.get_header())
Dims: (3, 4, 5, 2)
MRI Type: MRI_FLOAT
goodRASFlag: 1
delta: [ 3.7416575  3.7416575  3.7416575]
Mdc: 
[[ 0.26726124  0.53452247  0.80178374]
 [ 0.53452247  0.80178374  0.26726124]
 [ 0.80178374  0.26726124  0.53452247]]
Pxyz_c: [ 0.  0.  0.]
mrparms: [ 2.  0.  0.  0.]

@larsoner
Copy link
Contributor Author

@matthew-brett the issue seems to be that those delta values don't match the info obtained using mri_info --res or mri_info --csize/rsize/ssize, and I don't know enough about FreeSurfer and its file I/O formatting to know what the underlying issue is.

@larsoner
Copy link
Contributor Author

@matthew-brett I think I figured it out. The issue is that in update_header, hdr['delta'] is updated, and this sets the zooms to be 3.74 for each dimension. I'm not sure if we need a separate entry to hold the original value, or what...

@larsoner
Copy link
Contributor Author

@matthew-brett, this is what I now get:

>>> print(img.get_header())
Dims: (3, 4, 5, 2)
MRI Type: MRI_FLOAT
goodRASFlag: 1
delta: [ 1.  1.  1.]
Mdc: 
[[ 1.  2.  3.]
 [ 2.  3.  1.]
 [ 3.  1.  2.]]
Pxyz_c: [ 0.  0.  0.]
mrparms: [ 2.  0.  0.  0.]

Which seems to agree better with mri_info:

Volume information for nibabel/tests/data/test.mgz
          type: MGH
    dimensions: 3 x 4 x 5 x 2
   voxel sizes: 1.0000, 1.0000, 1.0000
          type: FLOAT (3)
           fov: 5.000
           dof: 0
        xstart: -1.5, xend: 1.5
        ystart: -2.0, yend: 2.0
        zstart: -2.5, zend: 2.5
            TR: 2.00 msec, TE: 0.00 msec, TI: 0.00 msec, flip angle: 0.00 degrees
       nframes: 2
       PhEncDir: UNKNOWN
ras xform present
    xform info: x_r =   1.0000, y_r =   2.0000, z_r =   3.0000, c_r =     0.0000
              : x_a =   2.0000, y_a =   3.0000, z_a =   1.0000, c_a =     0.0000
              : x_s =   3.0000, y_s =   1.0000, z_s =   2.0000, c_s =     0.0000

talairach xfm : 
Orientation   : SAR
Primary Slice Direction: sagittal

voxel to ras transform:
                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

voxel-to-ras determinant -18

ras to voxel transform:
               -0.2778   0.0556   0.3889     1.5000
                0.0556   0.3889  -0.2778     2.0000
                0.3889  -0.2778   0.0556     2.5000
                0.0000   0.0000   0.0000     1.0000

WDYT? To be perfectly honest, I don't understand why the deltas were being re-calculated, but my recent commit makes it so they aren't and it seems to match better...

@@ -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...

@matthew-brett
Copy link
Member

The reason for the update_header command is to deal with the situation where the user updates the affine.

The voxel sizes as usually defined are given by np.sqrt(np.sum(affine ** 2, axis=0))

http://mail.scipy.org/pipermail/nipy-devel/2011-March/005696.html

I'm trying to understand how this relates to the sizes in the MGH header.

@larsoner
Copy link
Contributor Author

@matthew-brett assuming the comments I've made in the test_mghformat.py are clear enough, this should be good for merge now.

@matthew-brett
Copy link
Member

Let's try and stick as close as possible to what mri_info does... I think the best thing is to do what the analyze converters do, and only reset the affine information in 'update_header', if the affine required is significantly different from the affine stored in the header... Sound reasonable?

@larsoner
Copy link
Contributor Author

That sounds reasonable, but I'm not sure how to define “significantly
different" in this context. Do you want to branch my PR and work from it
yourself? Might be faster than trying to explain in detail what you mean
only to have me probably bungle the implementation :)
On Mar 27, 2013 2:44 PM, "Matthew Brett" [email protected] wrote:

Let's try and stick as close as possible to what mri_info does... I think
the best thing is to do what the analyze converters do, and only reset the
affine information in 'update_header', if the affine required is
significantly different from the affine stored in the header... Sound
reasonable?


Reply to this email directly or view it on GitHubhttps://github.com//pull/164#issuecomment-15555081
.

@larsoner
Copy link
Contributor Author

By the way, the output now matches that of MRIRead.m in the Freesurfer MATLAB toolbox, so maybe we're okay...

@larsoner
Copy link
Contributor Author

But there is still something funny going on, since the order shouldn't be 1, 2, 0, it should be 0, 2, 1 in the vox2ras-tkr matrix. This is frustrating :(

larsoner and others added 3 commits March 27, 2013 15:52
We were having trouble with the header updating on load of an image, and
changing the underlying header information.

nipy#164

Use the Analyze trick to only update the header if the affine has
changed.
No we are not resetting the deltas on loading the image, we should get
zooms of 1 from the test image.
@matthew-brett
Copy link
Member

larsoner#1

But I think I've messed up your test :)

@matthew-brett
Copy link
Member

It's annoying that the test.mgz file has a symmetric 3x3 component. I wonder if we're transposing it or something.

@larsoner
Copy link
Contributor Author

Doug said that the test.mgz file is basically invalid/corrupted. Perhaps we should replace the direction cosines with valid ones in this PR, too?

@larsoner
Copy link
Contributor Author

By the way, it "breaking" the test actually probably means it's fixed. I had changed it to match the output of MRILoad.m, which wasn't the same as mri_info. What you did should make it match mri_info --vox2ras-tkr, I'll update the test to reflect that.

@larsoner
Copy link
Contributor Author

Well, tests pass over here with a match to mri_info outputs, so I'm happy. WDYT?

@larsoner
Copy link
Contributor Author

I like the comment quote, by the way: Luckily the test.mgz dataset had a bad set of cosine vectors.... I see you're an optimist :)

We weren't previously testing the transpose / not-transpose of the
cosine vectors because the test dataset had a symmetical 3x3 in the
affine.
@matthew-brett
Copy link
Member

:)

Added a test for cosine vector order... Are we happy now? Does it all make sense?

@larsoner
Copy link
Contributor Author

Yeah, makes sense to me. I'm satisfied with it.

@larsoner
Copy link
Contributor Author

Tests pass over here, by the way.

@matthew-brett
Copy link
Member

Thanks for taking the time to track all this stuff down...

@matthew-brett
Copy link
Member

Just to check - do you now understand this thing you asked about before? Sorry, I don't...

Based on the powerpoint linked from this page (http://surfer.nmr.mgh.harvard.edu/fswiki/CoordinateSystems)
I thought that the [-1, -1, 1] triplet that is used to fill in the upper left corner of the matrix would need to
be [-delta[0], -delta[2], delta[1]], but that wasn't what the output of mri_info --vox2ras-tkr gave me. 

@larsoner
Copy link
Contributor Author

Yeah, turns out that it was [-delta[0], -delta[2], delta[1]], consistent with the Freesurfer powerpoint, but we had been getting different values (which is fixed by this PR).

@matthew-brett
Copy link
Member

OK - makes sense.

matthew-brett added a commit that referenced this pull request Mar 28, 2013
ENH: Adding vox2ras_tkr

This adds the support requested in issue #163 for access to the vox2ras_tkr output from the MGH format.
@matthew-brett matthew-brett merged commit 4ef5573 into nipy:master Mar 28, 2013
@larsoner larsoner deleted the vox2ras_tkr branch October 24, 2014 19:02
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants