RF: Simplify MGHImage and add footer fields #569

effigies · 2017-10-25T14:36:46Z

This PR includes a number of simplifications to the MGHHeader and MGHImage data structures.

The most significant is subclassing MGHHeader from LabeledWrapStruct. It appeared to have been a reimplementation of LabeledWrapStruct in many places without the endianness option. I have simply raised an exception on attempts to instantiate as little-endian.

The next most significant change is moving from Mdc (matrix of direction cosines) and Pxyz_c (center point in world coordinates) to the more explicit x/y/z/c_ras column vectors. These are consistent with most FreeSurfer literature and program outputs, and correctly encode the shape of the direction cosines as column-major rather than row-major (needing to be transposed when reading/writing).

Somewhat confusingly, there are two ways of referring to the voxel and world coordinates, and XYZ means world in one, and voxel in the other, and in both cases the alternative includes "R" and "S".

Voxel	World
Column Row Slice	XYZ
XYZ	Right Anterior Superior

As CRS/XYZ appears to be FreeSurfer-unique, I'd prefer to settle on XYZ/RAS, and perhaps consider non-XYZ variable names when they can reduce confusion.

TODO:

Update tests with new header fields
Annotate header dtype with descriptions of fields (possibly further updating field names, shapes)
Improve footer handling (not targeting full features, but access to TR would be good)
Fix get_data_shape handling (see ENH: Add image slicing #550)

coveralls · 2017-10-28T15:50:11Z

Coverage increased (+0.02%) to 96.27% when pulling 8870068 on effigies:rf/mgh into b1bf5e7 on nipy:master.

codecov-io · 2017-10-28T15:50:17Z

Codecov Report

Merging #569 into master will increase coverage by 0.09%.
The diff coverage is 97.83%.

@@            Coverage Diff             @@
##           master     #569      +/-   ##
==========================================
+ Coverage   94.32%   94.41%   +0.09%     
==========================================
  Files         177      177              
  Lines       24697    24884     +187     
  Branches     2640     2656      +16     
==========================================
+ Hits        23295    23495     +200     
+ Misses        925      913      -12     
+ Partials      477      476       -1

Impacted Files	Coverage Δ
nibabel/spatialimages.py	`96.05% <100%> (ø)`	⬆️
nibabel/tests/test_spatialimages.py	`99.7% <100%> (ø)`	⬆️
nibabel/freesurfer/tests/test_mghformat.py	`98.44% <97.31%> (-1.56%)`	⬇️
nibabel/freesurfer/mghformat.py	`95.53% <98.42%> (+7.03%)`	⬆️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 2817d1b...8d720fe. Read the comment docs.

coveralls · 2017-10-30T14:30:45Z

Coverage increased (+0.06%) to 96.311% when pulling 21cbbcf on effigies:rf/mgh into b1bf5e7 on nipy:master.

effigies · 2017-10-30T14:39:14Z

This is ready for review, if anybody's got a few minutes.

matthew-brett

Main question - is it really worth renaming Mdc and the other fields?

matthew-brett · 2017-10-30T14:46:28Z

nibabel/freesurfer/mghformat.py

@@ -84,6 +92,7 @@ class MGHHeader(object):

    def __init__(self,
                 binaryblock=None,
+                 endianness='>',


Does it make sense to allow passing this flag? It also changes the API.

I did this because of places where LabeledWrapStruct calls __init__ assuming the type signature of the base class. I can instead make sure to override any methods that make such calls, just figured this would be a little shorter.

I played with the idea of making WrapStruct have a parent or similar that had forced endianness, but I seem to remember it got a bit messy.

matthew-brett · 2017-10-30T14:47:40Z

nibabel/freesurfer/mghformat.py

-
-    def check_fix(self):
-        ''' Pass. maybe for now'''
+        return klass(hdr_str + ftr_str, check=check)

    def get_affine(self):
        ''' Get the affine transform from the header information.


Blank line after first line (sorry, I know that wasn't present before either).

matthew-brett · 2017-10-30T14:48:50Z

nibabel/freesurfer/mghformat.py

+        MdcD = np.hstack((hdr['x_ras'], hdr['y_ras'], hdr['z_ras'])) * hdr['voxelsize']
+        vol_center = MdcD.dot(hdr['dims'][:3].reshape(-1, 1)) / 2
+        affine[:3, :3] = MdcD
+        affine[:3, [3]] = hdr['c_ras'] - vol_center


3 rather than [3] ?

The RHS is a column, so the [3] keeps the LHS expecting a column. Is there a trick I don't know here?

matthew-brett · 2017-10-30T14:50:19Z

nibabel/freesurfer/mghformat.py

-        hdr = self._header_data
-        zooms = hdr['delta']
-        return tuple(zooms[:])
+        # Do not return time zoom (TR) if 3D image


Maybe worth a note in the docstring about the time zoom and 3D?

matthew-brett · 2017-10-30T14:52:04Z

nibabel/freesurfer/mghformat.py

+    ('dims', '>i4', (4,)),          # 4; width, height, depth, nframes
+    ('type', '>i4'),                # 20; data type
+    ('dof', '>i4'),                 # 24; degrees of freedom
+    ('ras_good', '>i2'),            # 28; *_ras fields valid


I guess these renamings are API changes - our our Freesurfer colleagues happy with these name changes? Any way to offer the old names as alternatives for a while?

Yeah, I can re-add a _header_data structure (which I've replaced with LabeledWrapStruct._structarr) that will expose the old fields and raise a DeprecationWarning.

matthew-brett · 2017-10-30T14:55:19Z

nibabel/freesurfer/mghformat.py

+    ('dof', '>i4'),                 # 24; degrees of freedom
+    ('ras_good', '>i2'),            # 28; *_ras fields valid
+    ('voxelsize', '>f4', (3,)),     # 30; zooms (X, Y, Z)
+    ('x_ras', '>f4', (3, 1)),       # 42; X direction cosine column


(3, 1) rather than (3,)?

Is this split into three columns really useful? Compared to the native 3, 3 size?

By including all three columns as a single, row-ordered Mdc matrix, we're exposing a transposed matrix. It was being transposed correctly in and out, but if someone does decide to use it directly, they need to know this. I prefer explicit to silently misleading...

I could use (3,) if that's important. But it will mean users knowing that they're columns and not row vectors (can add to the docs).

For context, this is the comparable portion of the C data structure, which is just 16 individual float fields: https://github.com/freesurfer/freesurfer/blob/5367eec84621b331271e97f4e1f1502c3f6d0d87/include/mri.h#L184-L191

OK - fine to leave as is (3,1) etc.

matthew-brett · 2017-10-30T15:00:16Z

nibabel/freesurfer/mghformat.py

+        hdr = self._structarr
+        MdcD = np.hstack((hdr['x_ras'], hdr['y_ras'], hdr['z_ras'])) * hdr['voxelsize']
+        vol_center = MdcD.dot(hdr['dims'][:3].reshape(-1, 1)) / 2
+        affine[:3, :3] = MdcD


return nib.affines.from_matvec(MdcD, hdr['c_ras'] - vol_center)

?

matthew-brett · 2017-10-30T15:10:29Z

nibabel/freesurfer/mghformat.py

+
+        Ignores byte order; always big endian
+        '''
+        structarr = super(MGHHeader, klass).default_structarr()


Here might be worth checking and error if endian not big.

matthew-brett · 2017-10-30T15:12:10Z

nibabel/freesurfer/mghformat.py

+
+        # Assign after we've had a chance to raise exceptions
+        hdr['voxelsize'] = voxelsize
+        hdr['x_ras'] = Mdc[:, [0]]


hdr['x_ras'], hdr['y_ras'], hdr['z_ras'] = Mdc.T

?

effigies · 2017-10-30T16:17:39Z

Re: Renaming, the biggest motivation is that a lot of the field names seem only to align with one slideshow that documents the spaces/transforms, and not at all with either the FreeSurfer source code or the outputs of mri_info or other FreeSurfer tools, let alone how anybody outside FreeSurfer refers to these spaces.

So my strong preference is to move to a more consistent naming scheme (I'm open to alternative suggestions), but I'm also perfectly happy to provide a backwards-compatible interface to access the data using the old field names to get back the old shapes.

matthew-brett

I guess it would be good to get more feedback from Freesurferers out there.

matthew-brett · 2017-10-31T15:15:44Z

nibabel/freesurfer/mghformat.py

+    ('dof', '>i4'),                 # 24; degrees of freedom
+    ('ras_good', '>i2'),            # 28; *_ras fields valid
+    ('voxelsize', '>f4', (3,)),     # 30; zooms (X, Y, Z)
+    ('x_ras', '>f4', (3, 1)),       # 42; X direction cosine column


OK - fine to leave as is (3,1) etc.

matthew-brett · 2017-10-31T15:16:27Z

nibabel/freesurfer/mghformat.py

@@ -84,6 +92,7 @@ class MGHHeader(object):

    def __init__(self,
                 binaryblock=None,
+                 endianness='>',


I played with the idea of making WrapStruct have a parent or similar that had forced endianness, but I seem to remember it got a bit messy.

effigies · 2017-11-01T01:16:03Z

Okay. I've pinged the Python neuroimaging list and the FreeSurfer support list. Hopefully we can find somebody willing to share their opinions on the Internet. :-)

mwaskom · 2017-11-01T13:52:38Z

I don't have a strong opinion here, but I would guess it's more likely that a user has encountered the fscoordinates slide deck than read through the freesurfer source code.

satra · 2017-11-01T15:14:28Z

@effigies - perhaps chatting with doug would be a good idea as well. also folks at corticometrics may have some suggestions as well.

however, i do not have a strong opinion either. since i rarely do much with those files than reading into an array, doing some processing, and writing. or just visualization.

effigies · 2017-11-01T15:35:46Z

@satra Can I expect them to respond to the post on the FreeSurfer list, or is there a preferred way of contacting them more directly?

satra · 2017-11-01T22:46:34Z

@effigies - i hadn't seen that the mail was on the freesurfer list. hopefully they respond - otherwise it's your call :)

effigies · 2017-11-07T15:47:16Z

Okay, it sounds like there's a slight preference from @mwaskom to follow the FS documentary slideshow and a preference (of unknown strength) from @matthew-brett to avoid changing names without a Very Good Reason. I'll defer to this for now and revert Mdc and Pxyz_c, because I don't think it's critical to my more pressing goals here of getting zooms for 4D images to include TR (to match the behavior of Analyze-like images) and enforcing the 3/4D constraint.

However, I do think I would like to move to more "standard" FreeSurfer naming, because this construct ("volume geometry", or dimensions + x/y/z/c_ras) shows up in several places, including m3z warp files, LTA source and destination geometries (see #565), and in a text appendage to surface files:

nibabel/nibabel/freesurfer/io.py

Lines 50 to 73 in 2139ce0

    
           def _read_volume_info(fobj): 
        
               """Helper for reading the footer from a surface file.""" 
        
               volume_info = OrderedDict() 
        
               head = np.fromfile(fobj, '>i4', 1) 
        
               if not np.array_equal(head, [20]):  # Read two bytes more 
        
                   head = np.concatenate([head, np.fromfile(fobj, '>i4', 2)]) 
        
                   if not np.array_equal(head, [2, 0, 20]): 
        
                       warnings.warn("Unknown extension code.") 
        
                       return volume_info 
        
               volume_info['head'] = head 
        
               for key in ['valid', 'filename', 'volume', 'voxelsize', 'xras', 'yras', 
        
                           'zras', 'cras']: 
        
                   pair = fobj.readline().decode('utf-8').split('=') 
        
                   if pair[0].strip() != key or len(pair) != 2: 
        
                       raise IOError('Error parsing volume info.') 
        
                   if key in ('valid', 'filename'): 
        
                       volume_info[key] = pair[1].strip() 
        
                   elif key == 'volume': 
        
                       volume_info[key] = np.array(pair[1].split()).astype(int) 
        
                   else: 
        
                       volume_info[key] = np.array(pair[1].split()).astype(float) 
        
               # Ignore the rest 
        
               return volume_info

nibabel/nibabel/freesurfer/io.py

Lines 462 to 488 in 2139ce0

    
           def _serialize_volume_info(volume_info): 
        
               """Helper for serializing the volume info.""" 
        
               keys = ['head', 'valid', 'filename', 'volume', 'voxelsize', 'xras', 'yras', 
        
                       'zras', 'cras'] 
        
               diff = set(volume_info.keys()).difference(keys) 
        
               if len(diff) > 0: 
        
                   raise ValueError('Invalid volume info: %s.' % diff.pop()) 
        
               strings = list() 
        
               for key in keys: 
        
                   if key == 'head': 
        
                       if not (np.array_equal(volume_info[key], [20]) or np.array_equal( 
        
                               volume_info[key], [2, 0, 20])): 
        
                           warnings.warn("Unknown extension code.") 
        
                       strings.append(np.array(volume_info[key], dtype='>i4').tostring()) 
        
                   elif key in ('valid', 'filename'): 
        
                       val = volume_info[key] 
        
                       strings.append('{0} = {1}\n'.format(key, val).encode('utf-8')) 
        
                   elif key == 'volume': 
        
                       val = volume_info[key] 
        
                       strings.append('{0} = {1} {2} {3}\n'.format( 
        
                           key, val[0], val[1], val[2]).encode('utf-8')) 
        
                   else: 
        
                       val = volume_info[key] 
        
                       strings.append('{0} = {1:0.10g} {2:0.10g} {3:0.10g}\n'.format( 
        
                           key.ljust(6), val[0], val[1], val[2]).encode('utf-8')) 
        
               return b''.join(strings)

But I think it does make sense to postpone the discussion of field naming to an actual attempt to provide a common interface to volume geometries (which is actually re-implemented in several FreeSurfer data structures).

coveralls · 2017-11-09T16:27:27Z

Coverage increased (+0.004%) to 96.257% when pulling 1bad727 on effigies:rf/mgh into b1bf5e7 on nipy:master.

matthew-brett · 2017-11-10T14:36:00Z

Chris - sorry - I should have said - I don't care too much about the name changes, especially if there's a deprecation route. I should have said before, but if you think these names are more canonical, that's OK - we can always explain in the docstring and comments what the relationship is to the slide deck.

effigies · 2017-11-10T17:11:49Z

Okay. I think I'll still push it off for now. If/when we add support for LTA/M3Z etc, it will be time enough to consider a consistent naming scheme, rather than deprecating the old one now and finding later that we would have preferred something slightly different.

I will add docs though noting that Mdc is a transposed matrix. (It's a little annoying that I can't specify row/column order in the dtype.)

coveralls · 2017-11-23T01:53:45Z

Coverage increased (+0.003%) to 96.258% when pulling 79a5876 on effigies:rf/mgh into d420dcf on nipy:master.

coveralls · 2017-11-23T15:52:26Z

Coverage increased (+0.003%) to 96.258% when pulling 821b6b6 on effigies:rf/mgh into 2817d1b on nipy:master.

coveralls · 2017-12-09T14:56:58Z

Coverage increased (+0.004%) to 96.258% when pulling fede660 on effigies:rf/mgh into 2817d1b on nipy:master.

…block

coveralls · 2017-12-09T15:39:07Z

Coverage decreased (-0.02%) to 96.236% when pulling 0d4cdd0 on effigies:rf/mgh into 2817d1b on nipy:master.

coveralls · 2017-12-09T16:55:18Z

Coverage increased (+0.04%) to 96.293% when pulling 664f11c on effigies:rf/mgh into 2817d1b on nipy:master.

effigies · 2017-12-09T18:29:29Z

@matthew-brett Got back to this. I think I've addressed all of your comments, and I've fleshed out some more tests. Ready for another review.

coveralls · 2017-12-09T18:58:35Z

Coverage increased (+0.08%) to 96.331% when pulling fa73ad4 on effigies:rf/mgh into 2817d1b on nipy:master.

effigies · 2017-12-18T15:07:18Z

Hi @matthew-brett, just bugging you about this. No worries if you're busy.

matthew-brett

Just a few small comments.

matthew-brett · 2017-12-18T16:25:54Z

nibabel/freesurfer/mghformat.py


    def set_zooms(self, zooms):
        ''' Set zooms into header fields

-        See docstring for ``get_zooms`` for examples
+        Sets the spaing of voxels in the x, y, and z dimensions.


Typo 'spacing'

matthew-brett · 2017-12-18T16:26:36Z

nibabel/freesurfer/mghformat.py


    def set_zooms(self, zooms):
        ''' Set zooms into header fields

-        See docstring for ``get_zooms`` for examples
+        Sets the spaing of voxels in the x, y, and z dimensions.
+        For four-dimensional files, a temporal zoom (repetition time, or TR, in


Ouch - it's really in milliseconds? That's unfortunate (compared to nifti).

Annoyingly. Fortunately #567 proposes a way to get normalized zooms. (It was working on that that led to this PR in the first place.)

Reference: https://surfer.nmr.mgh.harvard.edu/fswiki/FsTutorial/MghFormat#line-82

Add reference to docstring?

coveralls · 2017-12-18T17:02:27Z

Coverage increased (+0.08%) to 96.331% when pulling c066cea on effigies:rf/mgh into 2817d1b on nipy:master.

matthew-brett · 2017-12-18T17:13:33Z

Thanks for your patience, in it goes.

coveralls · 2017-12-18T17:24:15Z

Coverage increased (+0.08%) to 96.331% when pulling 8d720fe on effigies:rf/mgh into 2817d1b on nipy:master.

matthew-brett · 2017-12-18T17:51:20Z

Oops - a couple of PEP8 failures : https://travis-ci.org/nipy/nibabel/builds/318192252

effigies · 2017-12-18T17:53:41Z

Fixed in 7afc6d2. Feel free to cherry-pick into master if #550 isn't ready to merge.

effigies force-pushed the rf/mgh branch 3 times, most recently from 5241b1b to 832f584 Compare October 28, 2017 14:43

effigies changed the title ~~[WIP] RF: Simplify MGHImage~~ RF: Simplify MGHImage and add footer fields Oct 28, 2017

matthew-brett reviewed Oct 30, 2017

View reviewed changes

matthew-brett reviewed Oct 31, 2017

View reviewed changes

effigies force-pushed the rf/mgh branch from 1bad727 to 79a5876 Compare November 23, 2017 01:17

effigies added 8 commits November 23, 2017 09:09

FIX/ENH: get/set_data_shape more robustly

bb56000

RF: Subclass MGHHeader from LabeledWrapStruct

d7b472a

RF: Remove one-line private method

415ee10

FIX: goodRASFlag check, update

4e520ac

RF: Remove unused _load_cache

a83f80b

FIX: Label x/y/z_ras correctly

ae5c085

RF: Make direction cosines explicit, not transposed

6b6e512

RF: Pxyz_c to c_ras

971c98a

effigies added 4 commits November 23, 2017 09:09

ENH: Check for expected byte order

16c8bb6

ENH: Re-add _header_data with deprecation warning

743db3b

FIX: Return on get, not on set

495bbce

REVERT renamed/reshaped header fields

821b6b6

effigies force-pushed the rf/mgh branch from 79a5876 to 821b6b6 Compare November 23, 2017 14:11

RF/TEST: Check endianness rather than ignoring

fede660

RF: Remove endianness, override copy, as_byteswapped, diagnose_binary…

0d4cdd0

…block

effigies added 2 commits December 9, 2017 11:21

FIX: Do not double-set attributes

93f3c70

TEST: TestWrapStruct features, set_zooms, deprecated fields

664f11c

effigies added 2 commits December 9, 2017 13:00

TEST: Add and test version check

3f62445

TEST: Test MGHHeader._HeaderData fully

fa73ad4

effigies force-pushed the rf/mgh branch from 30fc9ab to fa73ad4 Compare December 9, 2017 18:28

matthew-brett reviewed Dec 18, 2017

View reviewed changes

STY: Typo

c066cea

DOC: Add reference to FreeSurfer wiki

8d720fe

matthew-brett merged commit b646d5f into nipy:master Dec 18, 2017

effigies deleted the rf/mgh branch December 18, 2017 17:19

effigies mentioned this pull request Dec 18, 2017

ENH: Add image slicing #550

Merged

RF: Simplify MGHImage and add footer fields #569

RF: Simplify MGHImage and add footer fields #569

Conversation

effigies commented Oct 25, 2017 • edited Loading

coveralls commented Oct 28, 2017

codecov-io commented Oct 28, 2017 • edited Loading

Codecov Report

coveralls commented Oct 30, 2017

effigies commented Oct 30, 2017

matthew-brett left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

effigies commented Oct 30, 2017

matthew-brett left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

effigies commented Nov 1, 2017

mwaskom commented Nov 1, 2017 • edited Loading

satra commented Nov 1, 2017

effigies commented Nov 1, 2017 • edited Loading

satra commented Nov 1, 2017

effigies commented Nov 7, 2017

coveralls commented Nov 9, 2017

matthew-brett commented Nov 10, 2017

effigies commented Nov 10, 2017

coveralls commented Nov 23, 2017

coveralls commented Nov 23, 2017

coveralls commented Dec 9, 2017

coveralls commented Dec 9, 2017

coveralls commented Dec 9, 2017

effigies commented Dec 9, 2017

coveralls commented Dec 9, 2017 • edited Loading

effigies commented Dec 18, 2017

matthew-brett left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

coveralls commented Dec 18, 2017

matthew-brett commented Dec 18, 2017

coveralls commented Dec 18, 2017

matthew-brett commented Dec 18, 2017

effigies commented Dec 18, 2017

effigies commented Oct 25, 2017 •

edited

Loading

codecov-io commented Oct 28, 2017 •

edited

Loading

mwaskom commented Nov 1, 2017 •

edited

Loading

effigies commented Nov 1, 2017 •

edited

Loading

coveralls commented Dec 9, 2017 •

edited

Loading