How to load FreeSurfer surface files in correct position using Python script

Development
1

Hello, am writing a module which uses FreeSurfer segmentations and surface files. Empirically and based on another discourse post (Nifti file and freesurfer pial file mismatch in slicer) I had previously developed a transform to put the pial surface files back in the coordinate system of the segmented labelmap, based on the idea of ‘uncentering’.

def createSurfTransform(self, referenceNode, transformName='FreeSurfer_to_Parcellation_Transform'):
  '''Create the transform which takes the FreeSurfer surfaces into the coordinate system of the 
  reference node.  The transform is a translation from the origin to the center of the reference 
   volume. 
   '''
 
   bounds = [0,0,0,0,0,0] # xmin,xmax,ymin,ymax,zmin,zmax, preallocated
 
   referenceNode.GetBounds(bounds) # this is always the bounds of the untransformed object
   tr = [(bounds[1]+bounds[0])/2,(bounds[3]+bounds[2])/2, (bounds[5]+bounds[4])/2]
   # Create new transformation matrix with this translation
   transformMatrix = vtk.vtkMatrix4x4()
   transformMatrix.SetElement(0,3,tr[0])
   transformMatrix.SetElement(1,3,tr[1])
   transformMatrix.SetElement(2,3,tr[2])
   # Create transform node with given name
  transformNode = slicer.mrmlScene.AddNewNodeByClass('vtkMRMLTransformNode',transformName) 
   transformNode.SetAndObserveMatrixTransformToParent(transformMatrix) 
   return transformNode

However, the transform this code generates seems to differ from the one generated by SlicerFreeSurfer by exactly 0.5 in each dimension. The SlicerFreeSurfer transform does appear to line up the surface better than my transform. I see there is a note in the code about a shift of 0.5 to move to pixel centers (https://github.com/PerkLab/SlicerFreeSurfer/blob/910ee2c8adbce1c2d472e658f4cbffa982618cd9/FreeSurferImporter/MRML/vtkMRMLFreeSurferModelStorageNode.cxx#L37-L39). However, I don’t see any actual code associated with this comment. Also, I don’t understand why the shift would be 0.5 rather than 0.5*voxel dimension. The reference volumes I am using are usually approximately but not exactly 1mm isotropic voxels. For example, the one I am looking at now has voxels which are .97x.97x1mm, but the offset between the transform I generate and the one SlicerFreeSurfer generates is exactly 0.5mm in all 3 dimensions.

I don’t currently understand the logic of why a half-voxel shift is needed, and it sounds wrong to me if the shift is half a mm regardless of voxel size. I don’t argue with the fact that the surfaces line up better with the image with an offset than without, but I’d like to understand why.

In a related issue, I don’t see how to load a FreeSurfer surface file from a scripted module and specify the reference volume. I’d like to have the same thing happen as dropping a (for example) “lh.pial” file on the Slicer window and choosing the appropriate reference volume node. I have tried

slicer.util.loadNodeFromFile('lh.pial', 'FreeSurferModelFile')

but while this loads the model, it does not generate the transform. I am guessing that I could specify the reference volume somehow through the properties argument, but I haven’t been able to guess the right property name. I’ve tried properties={‘ReferenceVolume’: ‘myReferenceNodeName’}, and similar, trying ‘Referencevolume’, ‘Reference volume’, and I think a few other variations, none of which seemed to have any effect.

Bottom line: I’d love to understand why the SlicerFreeSurfer-generated transform is correct, but if I can’t manage that, I’d at least like to be able to be able to generate the same transform from a scripted module that is generated by dragging and dropping a specified surface file onto the slicer main window and choosing the appropriate reference volume.

2

Offset of 0.5 is appropriate if the value is in voxel coordinate system.

From this data set you won’t learn much. Could you create a data set with FreeSurfer that has much larger voxels than 1mm so that we can see subpixel size misalignments more clearly and make sure we don’t mistake offsets in voxel/physical coordinate systems?

To make FreeSurfer model loader take into account the reference volume node, specify referenceVolumeID in additional property:

modelNode = slicer.util.loadNodeFromFile('lh.pial', 'FreeSurferModelFile', {'referenceVolumeID': volumeNode.GetID()})

Alternatively, whenever you want to transform a model to RAS coordinate system you can call slicer.modules.freesurferimporter.logic().TransformFreeSurferModelToWorld().

2 Likes
3

Thanks, I’ll try processing a volume with larger and asymmetric voxels, and we can see how it comes out. I should have something to post back tomorrow.

4

OK, I tried resampling a volume to 1.5mm voxels and re-running FreeSurfer’s recon-all. The automated transform generated by FreeSurferImporter differs in it’s offset from my code by approximately 0.75 mm rather than approximately 0.5 mm, so it seems clear that it is using half a voxel rather than always half a millimeter. It is again true, that the FreeSurferImporter-generated transform does a better job of aligning the model with the labelmap volume; my transform results in outlines which are a little too right, anterior, and superior. I have also recalled that a volume’s IJK origin is generally the center of the 0,0,0 voxel, rather than the corner of the imaged volume, and that this is why a half-voxel offset is conceptually needed. So, I am satisfied that the scaling behavior is appropriate and I’ll just use the transform generated by FreeSurferImporter from now on. Thanks!

1 Like
5

In Slicer >=5.6, the reference volume is no longer needed, and will indeed cause a failure to load. The correct version of this code today is simpler:

modelNode = slicer.util.loadNodeFromFile('lh.pial')

The correct coordinates are applied automatically during loading (which now happens via NiBabelIO, see BUG: Fix FreeSurfer surface loading transformation by switching to Ni… · PerkLab/SlicerFreeSurfer@12c8ffd (github.com)). I’m just adding this note here because I just finished debugging why this code no longer worked in Slicer 5.6, and perhaps it will be helpful for others.

1 Like