I don't use mgz format - can you attach some sample data in both mgz and nii format that illustrates the issue? (or a link to a location where the data can be downloaded).

Thanks

Thanks for looking into this Steve.

take a look in
https://www.dropbox.com/sh/j1gtjueo9cemw75/15MXkJYPlw

there are the T1.mgz and T1.mgh files from the BERT subject that comes with freesurfer as well as the T1.nii.gz file that I made from the T1.mgz file.

There is also a simple nibabel script that seems to not have the same issues.

I imagine that the culprit is lines 309-320 of itkMGHImageIO.cxx but I haven't really dug into this.

The problem comes from not converting the center read from the MGH header from RAS to LAS...

I added the following before the loop in line 309 of itkMGHImageIO.cxx, which I think solved the problem

//convert C to from RAS to LPS
c[0] = -1;
c[1] = -1;

we should also remove the comment that this only works in some orientations as I think it should work all the time.

Issue targeted for 4.2.0, re-target if it can't be addressed by sept 1st.

I haven't looked at freesurfer dataformats in quite a while so I would not be surprised if the formats have been updated but the changes are not reflected in the reader used in slicer. The reader was provided by the MGH freesurfer team, so hopefully it accounted for the coordinate transforms as they existed at the time. Also I know that several people had used the combination of slicer and freesurfer, but I never heard of a systematic issue like this.

I was able to dig up some old data that I have and was able to load it into slicer to make the attached screenshot from today. As you can see, these line up perfectly:

https://dl.dropbox.com/u/1686930/2422-freesurfer/lh.pial
https://dl.dropbox.com/u/1686930/2422-freesurfer/aseg_edited.mgh

So I'm hesitant to change the code unless we know what is different in the datasets we are working with.

So there are two reasons I think that this is a bug and not related to changes in the file format specs.

1) the aseg_edited.mgh file you uploaded overlays and is unaffected by this bug because the center is (0,0,0).

In [8]: im = nibabel.load('aseg_edited.mgh')
In [9]: print im.get_header()
<class 'nibabel.freesurfer.mghformat.MGHHeader'>
Dims: (256, 256, 256)
MRI Type: MRI_UCHAR
goodRASFlag: 1
delta: [ 1. 1. 1.]
Mdc:
[[-1. 0. 0.]
[ 0. 0. -1.]
[ 0. 1. 0.]]
Pxyz_c: [ 0. 0. 0.] <--- c in itkMGHImageIO
mrparms: [ 0. 0. 0. 0.]

if you put a print in itkMGHImageIO.cxx you'll see that c is 0,0,0.
this translates to an RAS location of the (128,128,128) voxel of 0,0,0, which is exactly the world coordinate frame expected by the lh.pial file.

When c is non zero, this bug shows up which causes problems, however even if we fix it there will still be issues with overlaying surfaces.

2) If we look at itkMGHImageIO.cxx it seems pretty clear that c needs to be expressed in LPS as opposed to RAS. Lines 263 and 264 convert matrix into LPS. matrix is then used with 'c' to compute the origin in lines 316-318, but c is still in RAS. This will cause problems whenever c is non-zero. Its also pretty clear from the comment in lines 305-309 that the authors knew that the proposed code didn't cover all the cases.

We can see this problem in the bert subject that comes with freesurfer so I don't think that it is related to some recent change in file format.

I added Luke's suggested fix:

https://github.com/pieper/Slicer/commit/0bc1b6da85156e3e41bd627ad1df6fe19d196fd3

It's available in this topic branch:

https://github.com/pieper/Slicer/tree/2422-mgh-format-center

There is no regression with respect to the data that I have (results are still as shown on the 2012-08-22 screen shot) and the logic of the fix seems sound to me.

My only concern is that the reference data from the BERT subject still does not load correctly so I think there are more fixes needed.

topic branch merged as r20837