BUG: write_geometry and read_geometry don't seem to be in the same coordinate system

[kingjr]

As identified in mne-tools/mne-python#3281,

nibabel.freesurfer.io.write_geometry doesn't seem to write in the correct coordinate system.

Weirdly, this is visible in freeviewer, but not in the python viz functions we use in mne-python.

To reproduce the error (you have to copy/paste subjects/sample into subjects/sample_test since it modifies the files):

import os.path as op
from nibabel.freesurfer.io import read_geometry, write_geometry

 # I'm not sure where nibabel example is stored, so I'm using MNE.
from mne.datasets import sample
data_path = sample.data_path()

# The paths to freesurfer reconstructions
subjects_dir = data_path + '/subjects'
subject = 'sample_test'

fname = op.join(subjects_dir, subject, 'bem', 'inner_skull.surf')
coord, faces = read_geometry(fname)
write_geometry(fname, coord, faces)

in freeview, we now see a misalignment:

Weirdly, it doesn't appear when plotting in in python:

mne.viz.plot_bem plots as expected:

[effigies]

If I remember correctly, there is extra data at the end of a geometry file that doesn't appear in the spec, so we couldn't reproduce the original file completely. We do verify that, to the extent that load-save-load loop preserves the same numpy data structure. I suppose a possible solution is to make read_geometry return any extra data for re-appending, and allow write_geometry to append arbitrary bytes.

Another possibility is that freeview interprets triangle and quad files slightly differently, while nibabel reads them the same, and only outputs in triangle files. Could you give the first three bytes of the file in question?

import numpy as np
with open(fname, 'rb') as fobj:
    print(np.fromfile(fobj, ">u1", 3))

[kingjr]

Both the original and the modified bem files indicate [255 255 254]

[effigies]

From a quick look, it looks like bem surfaces are no longer supported by FreeSurfer, so it's possible that they're using an old format. Can you give several examples of bem surfaces? Perhaps the extra data at the end will have a pattern we can use.

[kingjr]

Here are three non modified :
https://dl.dropboxusercontent.com/u/3518236/outer_skull.surf
https://dl.dropboxusercontent.com/u/3518236/inner_skull.surf
https://dl.dropboxusercontent.com/u/3518236/outer_skin.surf

and here is the file after read/writing:
https://dl.dropboxusercontent.com/u/3518236/inner_skull_mod.surf

[agramfort]

pb confirmed. Freeview says: "Surface does not contain valid volume geometry information. It may not align with volumes and other surfaces." cc @jaeilepp this is related to your current efforts.

[effigies]

Okay, this is the extra for all three original files:

"""\x00\x00\x00\x14valid = 1  # volume info valid
filename = /cluster/fusion/mluessi/MNE-sample-data-build/MNE-sample-data/subjects/sample/mri/flash/parameter_maps/flash5_reg.mgz
volume = 256 256 256
voxelsize = 1.0000 1.0000 1.0000
xras   = -1.0000 0.0000 0.0000
yras   = 0.0000 0.0000 -1.0000
zras   = -0.0000 1.0000 0.0000
cras   = -5.2736 9.0391 -27.2880
"""

I'm guessing that the first four bytes are an extension code, since it doesn't indicate the length of the remaining data. The rest seems straightforward. So the question is mostly what to do. Treat it as an opaque block that should be passed back to write_geometry, or parse it and make it available to nibabel users?

[agramfort]

i guess we could return a volume_info that can also be passed on write. looking at different datasets only cras changes. I think it would be better to parse eg with volume_info being a dict with keys (volume, voxelsize etc.) and values as ndarray of shape (3,)

[effigies]

That seems reasonable to me. I'd probably go with an ordered dict, and for the first 4 bytes, just check that they're b'\x00\x00\x00\x14' and raise a warning if not.

[agramfort]

Are you taking a stab at it?

On Sun, Jun 5, 2016 at 2:05 PM -0700, "Chris Markiewicz" notifications@github.com wrote:

That seems reasonable to me. I'd probably go with an ordered dict, and for the first 4 bytes, just check that they're b'\x00\x00\x00\x14' and raise a warning if not.

—
You are receiving this because you commented.
Reply to this email directly, view it on GitHub, or mute the thread.

[effigies]

I've poked around a little in 28ee44a, and folks are welcome to work off of that to put together a proper PR. I found some other post-blocks, when starting on testing, so a bit more thinking may be required:

fsaverage/surf/lh.orig

"""\x00\x00\x00\x02\x00\x00\x00\x00\x00\x00\x00\x14valid = 0  # volume info invalid
filename = \x7f
volume = 0 3 0
voxelsize = -1.9991 0.0000 0.0000
xras   = 0.0000 -0.0000 0.0000
yras   = 0.0000 0.0000 0.0000
zras   = 0.0000 0.0000 0.0000
cras   = -1.9991 0.0000 -1.9991
\x00\x00\x00 \x00\x00\x00\x00\x00\x00\x00\x04G\x92l\xbe"""

fsaverage/surf/lh.sphere

"""\x00\x00\x00\x02\x00\x00\x00\x00\x00\x00\x00\x14valid = 0  # volume info invalid
filename = 
volume = 0 0 0
voxelsize = 0.0000 0.0000 0.0000
xras   = 0.0000 0.0000 0.0000
yras   = 0.0000 0.0000 0.0000
zras   = 0.0000 0.0000 0.0000
cras   = 0.0000 0.0000 0.0000
"""

fsaverage/surf/lh.inflated

"""\x00\x00\x00\x02\x00\x00\x00\x00\x00\x00\x00\x14valid = 0  # volume info invalid
filename = \x7f
volume = 0 3 0
voxelsize = -1.9991 0.0000 0.0000
xras   = 0.0000 -0.0000 0.0000
yras   = 0.0000 0.0000 0.0000
zras   = 0.0000 0.0000 0.0000
cras   = -1.9991 0.0000 -1.9991
\x00\x00\x00 \x00\x00\x00\x00\x00\x00\x00\x04G\xa0\x95\xfb\x00\x00\x00\x03\x00\x00\x00\x00\x00\x00\x01]mris_smooth -n 5 ./lh.white_avg ./lh.smoothwm ProgramVersion: $Name: stable4 $  TimeStamp: 2007/08/14-17:58:44-GMT  BuildTimeStamp: Aug 14 2007 05:37:36  CVS: $Id: mris_smooth.c,v 1.21 2007/01/01 18:16:19 fischl Exp $  User: nicks  Machine: node0354  Platform: Linux  PlatformVersion: 2.6.9-42.0.3.ELsmp  CompilerName: GCC  CompilerVersion: 30400  \x00\x00\x00\x00\x03\x00\x00\x00\x00\x00\x00\x01[mris_inflate -n 15 -f 1e-8 lh.smoothwm lh.inflated ProgramVersion: $Name:  $  TimeStamp: 2007/08/16-01:17:11-GMT  BuildTimeStamp: Aug  9 2007 18:37:20  CVS: $Id: mris_inflate.c,v 1.35 2007/08/16 00:49:14 fischl Exp $  User: nicks  Machine: minerva  Platform: Linux  PlatformVersion: 2.6.9-42.0.8.ELsmp  CompilerName: GCC  CompilerVersion: 30400  \x00"""

[agramfort]

cool. I commented on your commit. Can you open a PR?

[effigies]

I can, but I don't have the availability to quickly respond to feedback and get it in mergeable shape. It might make more sense for someone else to fork that branch themselves. (My impression was this is something that needs doing quickly.)

Let me know what works best for people.

[agramfort]

see https://github.com/nipy/nibabel/pull/458/files