Support for reading freesurfer's realras property

[jduemose]

Hi

Does anyone know why nibabel.freesurfer.io.read_geometry only allows either [20] or [2, 0, 20] as the first entry in the metadata section of a FreeSurfer surface file (referred to as head in the metadata dict)? The reason I am asking is that I am interested in being able to tell whether a surface is stored in surface RAS (tkregister space) or scanner RAS. Currently, it seems that read_geometry only supports reading the metadata when a surface is in surface RAS.

To test this, I tried to convert a surface file (created during a recon-all run) to scanner RAS like this

mris_convert --to-scanner lh.white lh.white.scanner

nibabel reads the metadata of lh.white just fine, however, for lh.white.scanner I get

.../nibabel/freesurfer/io.py:62: UserWarning: Unknown extension code. warnings.warn('Unknown extension code.') .../nibabel/freesurfer/io.py:188: UserWarning: No volume information contained in the file warnings.warn('No volume information contained in the file')

It turns out that it reads the head property just fine in both cases. For lh.white it is [2, 0, 20] but for lh.white.scanner it is [2, 1, 20] which, from looking into the freesurfer code (e.g., here https://github.com/freesurfer/freesurfer/blob/920f33cade45b901f702192ace64b37ef2c4b3e1/utils/mrisurf_io.cpp#L5348), is what I would expect as it seems that these three numbers signify [TAG_OLD_USEREALRAS (= 2), value associated with this tag (= 0 or 1), TAG_OLD_SURF_GEOM (= 20)]. From the FS code, it is difficult to tell whether this is always the pattern or not but I was thinking that it might make sense to add [2, 1, 20] as an allowed value of head metadata - and perhaps even add it as an entry in the returned metadata dict (e.g., realras = True/False). Any thoughts? I can do a PR if needed.

[effigies]

The reason is that those were the only types of data anybody had available (or found to test with) when the feature was added. As the warning says, [2, 1, 20] is unknown to nibabel.

You're definitely welcome to propose a patch. I think just permitting [2, 1, 20] makes sense, as it encodes what needs encoding. The freesurfer.io functions are low-level data access and not intended to provide interpretation, such as whether the coordinates are in any particular space. If we had a surface data structure, similar to FreeSurfer's internal structure, then it would make sense to add the field as a boolean value.

[jduemose]

Alright, I'll make a PR with minimal changes, i.e., allowing [2,1,20] to be read/written. I'll just interpret the head array myself for now.

However, since reading the footer of an FS surface file is reading (tag, value) pairs, I think it would make sense to return a mapping of tag and values. Right now you get something like this when reading the metadata

OrderedDict([
              ('head', array([ 2,  0, 20], dtype=int32)),
              ('valid', '1  # volume info valid'),
              ('filename', 'wm.mgz'),
              ('volume', array([256, 256, 256])),
              ('voxelsize', array([1., 1., 1.])),
              ('xras', array([-1.00000012,  0.        ,  0.        ])),
              ('yras', array([ 0.        ,  0.        , -1.00000012])),
              ('zras', array([0.        , 1.00000012, 0.        ])),
              ('cras', array([  2.58947754,  27.69668579, -28.27163696]))]))

where head is (tag, value, tag) and then the values of latter tag (volume geometry) are expanded to the dict. head=[2,0,20] is hardly interpretable. It would be straight-forward to return something like this instead (or perhaps a typeddict or named tuple or whatever)

OrderedDict([
              ('real_ras', False),
              ('vol_geom', OrderedDict([
                            ('valid', '1  # volume info valid'),
                            ('filename', 'wm.mgz'),
                            ('volume', array([256, 256, 256])),
                            ('voxelsize', array([1., 1., 1.])),
                            ('xras', array([-1.00000012,  0.        ,  0.        ])),
                            ('yras', array([ 0.        ,  0.        , -1.00000012])),
                            ('zras', array([0.        , 1.00000012, 0.        ])),
                            ('cras', array([  2.58947754,  27.69668579, -28.27163696]))]))]))

For this, the API of the functions should probably be modified so it would be breaking. Perhaps for next major release. If you prefer to keep it low level, I understand. I just think that the current structure of the metadata is puzzling.

[jduemose]

PR is open #1420