Fsnative to fmriprep fsnative space

[ntardiff]

I am trying to understand how fmriprep brings the freesurfer fsnative surfaces into T1w space giftis so that I could run my own transformations. I've come up w/ something that seems to work in practice but conceptually seems wrong, and I'm wondering if anyone could provide insight into how this is working. As an example using lh.pial, and assuming fMRIprep has been run:

1. mris_convert lh.pial lh.pial.surf.gii
2. run code adapted from: https://fmriprep.org/en/1.0.3/_modules/fmriprep/interfaces/surf.html

import numpy as np
import nitransforms as nt
import nibabel as nb

in_surf = os.path.join(base_path,surf_path,'lh.pial.surf.gii')
tx_file = os.path.join(base_path,anat_path,'sub-XXX_ses-XX_acq-highresAP_from-fsnative_to-T1w_mode-image_xfm.txt')

transform = nt.io.itk.ITKLinearTransform.from_filename(tx_file).to_ras()
transform = np.linalg.inv(transform) 

img = nb.load(in_surf)

pointset = img.get_arrays_from_intent('NIFTI_INTENT_POINTSET')[0]
coords = pointset.data.T
c_ras_keys = ('VolGeomC_R', 'VolGeomC_A', 'VolGeomC_S')
ras = np.array([[float(pointset.metadata[key])]
                for key in c_ras_keys])
ones = np.ones((1, coords.shape[1]), dtype=coords.dtype)
# Apply C_RAS translation to coordinates, then transform
pointset.data = transform.dot(np.vstack((coords + ras, ones)))[:3].T.astype(coords.dtype)

secondary = nb.gifti.GiftiNVPairs('AnatomicalStructureSecondary', 'MidThickness')
geom_type = nb.gifti.GiftiNVPairs('GeometricType', 'Anatomical')
has_ass = has_geo = False
for nvpair in pointset.meta.data:
    # Remove C_RAS translation from metadata to avoid double-dipping in FreeSurfer
    if nvpair.name in c_ras_keys:
        nvpair.value = '0.000000'
    # Check for missing metadata
    elif nvpair.name == secondary.name:
        has_ass = True
    elif nvpair.name == geom_type.name:
        has_geo = True

img.to_filename(out_file)

After running this, my output surface has vertices matching sub-XX_ses-XX_acq-highresAP_hemi-L_pial.surf.gii.

BUT:

It is not clear to me why one should need to invert the stored transform, since shouldn’t it already be from fsnative to T1w coords? But the code I based this off of (fmriprep.interfaces.surf.normalize_surfs()) uses the LTA transform, which I believe should be from T1w to fsnative coords (see: nitransforms.io.lta — nitransforms 24.1.2 documentation).

Am I confusing something here?

Thanks,
Nathan

[ellisdg]

I'm not sure why it is labeled native to T1 when the inverse of the transform is being applied. Looking at the code, the inverse transform seems to be applied by smriprep using the normalize_surfs function. Retracing how fmriprep produces the surfaces:

1. fmriprep.workflows.base.init_single_subject_wf calls the init_anat_fit_wf from smriprep which runs the processing of the anatomical data including surface reconstruction.
2. smriprep.workflows.anatomical.init_anat_fit_wf calls the init_gifti_surfaces_wf and passes the surfaces from the FreeSurfer reconstruction for the purpose of converting them to GIFTI files in T1w space.
3. smriprep.workflows.surfaces.init_gifti_surfaces_wf calls mris_convert with to_scanner=True (or '--to-scanner' on the commandline) which converts the surfaces to GIFTI format in native space.
4. smriprep.workflows.surfaces.init_gifti_surfaces_wf calls the NormalizeSurf interface.
5. smriprep.interfaces.surf.NormalizeSurf calls the normalize_surfs function which applies the inverse transform and removes the VolGeom metadata.

I know this doesn't answer your question, but hopefully it helps with your applying your own transformations. If you are ok with using mris_convert, you can simply call mris_convert --to-scanner h.pial lh.pial.surf.gii instead of apply the C_RAS translation yourself, as that is what fmriprep does.

For reference, here are some links I found helpful:

• Chris Markiewicz wrote some python code on how to convert surface coordinates to native space.
• Relevant discussion on neurostars: https://neurostars.org/t/transform-freesurfer-mesh-coordinates-to-volume-coordinates/23693
• niworkflows seems to have the code which you linked as previously being in fmriprep (see here)