[RTM] ENH: Register FreeSurfer template to FMRIPREP template

.circle/data/ds005_outputs.txt

@@ -23,7 +23,8 @@ ds005/out/fmriprep/sub-01/anat/sub-01_T1w_space-MNI152NLin2009cAsym_class-GM_pro
 ds005/out/fmriprep/sub-01/anat/sub-01_T1w_space-MNI152NLin2009cAsym_class-WM_probtissue.nii.gz
 ds005/out/fmriprep/sub-01/anat/sub-01_T1w_space-MNI152NLin2009cAsym_dtissue.nii.gz
 ds005/out/fmriprep/sub-01/anat/sub-01_T1w_space-MNI152NLin2009cAsym_preproc.nii.gz
-ds005/out/fmriprep/sub-01/anat/sub-01_T1w_space-MNI152NLin2009cAsym_warp.h5
+ds005/out/fmriprep/sub-01/anat/sub-01_T1w_target-fsnative_affine.txt
+ds005/out/fmriprep/sub-01/anat/sub-01_T1w_target-MNI152NLin2009cAsym_warp.h5
 ds005/out/fmriprep/sub-01/func
 ds005/out/fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-01_bold_AROMAnoiseICs.csv
 ds005/out/fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-01_bold_confounds.tsv

.circle/data/ds054_outputs.txt

@@ -15,7 +15,7 @@ ds054/out/fmriprep/sub-100185/anat/sub-100185_T1w_space-MNI152NLin2009cAsym_clas
 ds054/out/fmriprep/sub-100185/anat/sub-100185_T1w_space-MNI152NLin2009cAsym_class-WM_probtissue.nii.gz
 ds054/out/fmriprep/sub-100185/anat/sub-100185_T1w_space-MNI152NLin2009cAsym_dtissue.nii.gz
 ds054/out/fmriprep/sub-100185/anat/sub-100185_T1w_space-MNI152NLin2009cAsym_preproc.nii.gz
-ds054/out/fmriprep/sub-100185/anat/sub-100185_T1w_space-MNI152NLin2009cAsym_warp.h5
+ds054/out/fmriprep/sub-100185/anat/sub-100185_T1w_target-MNI152NLin2009cAsym_warp.h5
 ds054/out/fmriprep/sub-100185/func
 ds054/out/fmriprep/sub-100185/func/sub-100185_task-machinegame_run-01_bold_confounds.tsv
 ds054/out/fmriprep/sub-100185/func/sub-100185_task-machinegame_run-01_bold_space-MNI152NLin2009cAsym_brainmask.nii.gz

docs/workflows.rst

@@ -429,12 +429,12 @@ Derivatives
 There are additional files, called "Derivatives", outputted to ``<output dir>/fmriprep/sub-<subject_label>/``.
 See the `BIDS Derivatives`_ spec for more information.
 
-Derivatives related to t1w files are in the ``anat`` subfolder:
+Derivatives related to T1w files are in the ``anat`` subfolder:
 
 - ``*T1w_brainmask.nii.gz`` Brain mask derived using ANTS or AFNI, depending on the command flag ``--skull-strip-ants``
 - ``*T1w_space-MNI152NLin2009cAsym_brainmask.nii.gz`` Same as above, but in MNI space.
 - ``*T1w_dtissue.nii.gz`` Tissue class map derived using FAST.
-- ``*T1w_preproc.nii.gz`` Bias field corrected t1w file, using ANTS' N4BiasFieldCorrection
+- ``*T1w_preproc.nii.gz`` Bias field corrected T1w file, using ANTS' N4BiasFieldCorrection
 - ``*T1w_smoothwm.[LR].surf.gii`` Smoothed GrayWhite surfaces
 - ``*T1w_pial.[LR].surf.gii`` Pial surfaces
 - ``*T1w_midthickness.[LR].surf.gii`` MidThickness surfaces
@@ -443,7 +443,8 @@ Derivatives related to t1w files are in the ``anat`` subfolder:
 - ``*T1w_space-MNI152NLin2009cAsym_class-CSF_probtissue.nii.gz``
 - ``*T1w_space-MNI152NLin2009cAsym_class-GM_probtissue.nii.gz``
 - ``*T1w_space-MNI152NLin2009cAsym_class-WM_probtissue.nii.gz`` Probability tissue maps, transformed into MNI space
-- ``*T1w_target-MNI152NLin2009cAsym_warp.h5`` Composite (warp and affine) transform to transform t1w into MNI space
+- ``*T1w_target-MNI152NLin2009cAsym_warp.h5`` Composite (warp and affine) transform to transform T1w into MNI space
+- ``*T1w_target-fsnative_affine.txt`` Affine transform to transform T1w into ``fsnative`` space
 
 Derivatives related to EPI files are in the ``func`` subfolder.
 

fmriprep/interfaces/surf.py

@@ -13,13 +13,15 @@
 import numpy as np
 import nibabel as nb
 
-from niworkflows.nipype.interfaces.base import BaseInterfaceInputSpec, TraitedSpec, File, traits
+from niworkflows.nipype.interfaces.base import (
+    BaseInterfaceInputSpec, TraitedSpec, File, traits, isdefined)
 
 from niworkflows.interfaces.base import SimpleInterface
 
 
 class NormalizeSurfInputSpec(BaseInterfaceInputSpec):
     in_file = File(mandatory=True, exists=True, desc='Freesurfer-generated GIFTI file')
+    transform_file = File(exists=True, desc='FSL or LTA affine transform file')
 
 
 class NormalizeSurfOutputSpec(TraitedSpec):
@@ -60,7 +62,10 @@ class NormalizeSurf(SimpleInterface):
     output_spec = NormalizeSurfOutputSpec
 
     def _run_interface(self, runtime):
-        self._results['out_file'] = normalize_surfs(self.inputs.in_file)
+        transform_file = self.inputs.transform_file
+        if not isdefined(transform_file):
+            transform_file = None
+        self._results['out_file'] = normalize_surfs(self.inputs.in_file, transform_file)
         return runtime
 
 
@@ -176,7 +181,7 @@ def _run_interface(self, runtime):
         return runtime
 
 
-def normalize_surfs(in_file):
+def normalize_surfs(in_file, transform_file):
     """ Re-center GIFTI coordinates to fit align to native T1 space
 
     For midthickness surfaces, add MidThickness metadata
@@ -188,16 +193,17 @@ def normalize_surfs(in_file):
     """
 
     img = nb.load(in_file)
+    transform = load_transform(transform_file)
     pointset = img.get_arrays_from_intent('NIFTI_INTENT_POINTSET')[0]
-    coords = pointset.data
+    coords = pointset.data.T
     c_ras_keys = ('VolGeomC_R', 'VolGeomC_A', 'VolGeomC_S')
-    ras = np.array([float(pointset.metadata[key])
+    ras = np.array([[float(pointset.metadata[key])]
                     for key in c_ras_keys])
-    # Apply C_RAS translation to coordinates
-    pointset.data = (coords + ras).astype(coords.dtype)
+    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')
+    secondary = nb.gifti.GiftiNVPairs('AnatomicalStructureSecondary', 'MidThickness')
     geom_type = nb.gifti.GiftiNVPairs('GeometricType', 'Anatomical')
     has_ass = has_geo = False
     for nvpair in pointset.meta.data:
@@ -218,3 +224,32 @@ def normalize_surfs(in_file):
             pointset.meta.data.insert(2, geom_type)
     img.to_filename(fname)
     return os.path.abspath(fname)
+
+
+def load_transform(fname):
+    """Load affine transform from file
+
+    Parameters
+    ----------
+    fname : str or None
+        Filename of an LTA or FSL-style MAT transform file.
+        If ``None``, return an identity transform
+
+    Returns
+    -------
+    affine : (4, 4) numpy.ndarray
+    """
+    if fname is None:
+        return np.eye(4)
+
+    if fname.endswith('.mat'):
+        return np.loadtxt(fname)
+    elif fname.endswith('.lta'):
+        with open(fname, 'rb') as fobj:
+            for line in fobj:
+                if line.startswith(b'1 4 4'):
+                    break
+            lines = fobj.readlines()[:4]
+        return np.genfromtxt(lines)
+
+    raise ValueError("Unknown transform type; pass FSL (.mat) or LTA (.lta)")

fmriprep/workflows/anatomical.py

@@ -157,7 +157,7 @@ def init_anat_preproc_wf(skull_strip_ants, skull_strip_template, output_spaces,
         subject_id
             FreeSurfer subject ID
         t1_2_fsnative_reverse_transform
-            Affine transform from FreeSurfer subject space to T1w space
+            LTA-style affine matrix translating from FreeSurfer-conformed subject space to T1w
         surfaces
             GIFTI surfaces (gray/white boundary, midthickness, pial, inflated)
 
@@ -367,6 +367,7 @@ def len_above_thresh(in_list, threshold, longitudinal):
             ('mni_mask', 'inputnode.mni_mask'),
             ('mni_seg', 'inputnode.mni_seg'),
             ('mni_tpms', 'inputnode.mni_tpms'),
+            ('t1_2_fsnative_reverse_transform', 'inputnode.t1_2_fsnative_reverse_transform'),
             ('surfaces', 'inputnode.surfaces'),
         ]),
     ])
@@ -515,7 +516,7 @@ def init_surface_recon_wf(omp_nthreads, hires, name='surface_recon_wf'):
         subject_id
             FreeSurfer subject ID
         t1_2_fsnative_reverse_transform
-            FSL-style affine matrix translating from FreeSurfer T1.mgz to T1w
+            LTA-style affine matrix translating from FreeSurfer-conformed subject space to T1w
         surfaces
             GIFTI surfaces for gray/white matter boundary, pial surface,
             midthickness (or graymid) surface, and inflated surfaces
@@ -555,10 +556,8 @@ def init_surface_recon_wf(omp_nthreads, hires, name='surface_recon_wf'):
     skull_strip_extern = pe.Node(FSInjectBrainExtracted(), name='skull_strip_extern',
                                  run_without_submitting=True)
 
-    fs_transform = pe.Node(
-        fs.Tkregister2(fsl_out='freesurfer2subT1.mat', reg_header=True),
-        name='fs_transform')
-    fsl2lta = pe.Node(fs.utils.LTAConvert(out_lta=True), name='fsl2lta')
+    fsnative_2_t1_xfm = pe.Node(fs.RobustRegister(auto_sens=True, est_int_scale=True),
+                                name='fsnative_2_t1_xfm')
 
     autorecon_resume_wf = init_autorecon_resume_wf(omp_nthreads=omp_nthreads)
     gifti_surface_wf = init_gifti_surface_wf()
@@ -587,17 +586,16 @@ def init_surface_recon_wf(omp_nthreads, hires, name='surface_recon_wf'):
         (recon_config, autorecon_resume_wf, [('use_t2w', 'inputnode.use_T2')]),
         # Construct transform from FreeSurfer conformed image to FMRIPREP
         # reoriented image
-        (inputnode, fs_transform, [('t1w', 'target_image')]),
-        (autorecon1, fs_transform, [('T1', 'moving_image')]),
-        (inputnode, fsl2lta, [('t1w', 'target_file')]),
-        (autorecon1, fsl2lta, [('T1', 'source_file')]),
-        (fs_transform, fsl2lta, [('fsl_file', 'in_fsl')]),
+        (inputnode, fsnative_2_t1_xfm, [('t1w', 'target_file')]),
+        (autorecon1, fsnative_2_t1_xfm, [('T1', 'source_file')]),
+        (fsnative_2_t1_xfm, gifti_surface_wf, [
+            ('out_reg_file', 'inputnode.t1_2_fsnative_reverse_transform')]),
         # Output
         (autorecon_resume_wf, outputnode, [('outputnode.subjects_dir', 'subjects_dir'),
                                            ('outputnode.subject_id', 'subject_id'),
                                            ('outputnode.out_report', 'out_report')]),
         (gifti_surface_wf, outputnode, [('outputnode.surfaces', 'surfaces')]),
-        (fsl2lta, outputnode, [('out_lta', 't1_2_fsnative_reverse_transform')]),
+        (fsnative_2_t1_xfm, outputnode, [('out_reg_file', 't1_2_fsnative_reverse_transform')]),
     ])
 
     return workflow
@@ -750,6 +748,8 @@ def init_gifti_surface_wf(name='gifti_surface_wf'):
             FreeSurfer SUBJECTS_DIR
         subject_id
             FreeSurfer subject ID
+        t1_2_fsnative_reverse_transform
+            LTA formatted affine transform file (inverse)
 
     **Outputs**
 
@@ -760,7 +760,9 @@ def init_gifti_surface_wf(name='gifti_surface_wf'):
     """
     workflow = pe.Workflow(name=name)
 
-    inputnode = pe.Node(niu.IdentityInterface(['subjects_dir', 'subject_id']), name='inputnode')
+    inputnode = pe.Node(niu.IdentityInterface(['subjects_dir', 'subject_id',
+                                               't1_2_fsnative_reverse_transform']),
+                        name='inputnode')
     outputnode = pe.Node(niu.IdentityInterface(['surfaces']), name='outputnode')
 
     get_surfaces = pe.Node(nio.FreeSurferSource(), name='get_surfaces')
@@ -795,6 +797,7 @@ def init_gifti_surface_wf(name='gifti_surface_wf'):
         (save_midthickness, surface_list, [('out_file', 'in4')]),
         (surface_list, fs_2_gii, [('out', 'in_file')]),
         (fs_2_gii, fix_surfs, [('converted', 'in_file')]),
+        (inputnode, fix_surfs, [('t1_2_fsnative_reverse_transform', 'transform_file')]),
         (fix_surfs, outputnode, [('out_file', 'surfaces')]),
     ])
     return workflow
@@ -870,7 +873,8 @@ def init_anat_derivatives_wf(output_dir, output_spaces, template, freesurfer,
         niu.IdentityInterface(
             fields=['source_file', 't1_preproc', 't1_mask', 't1_seg', 't1_tpms',
                     't1_2_mni_forward_transform', 't1_2_mni', 'mni_mask',
-                    'mni_seg', 'mni_tpms', 'surfaces']),
+                    'mni_seg', 'mni_tpms',
+                    't1_2_fsnative_reverse_transform', 'surfaces']),
         name='inputnode')
 
     ds_t1_preproc = pe.Node(
@@ -913,10 +917,18 @@ def init_anat_derivatives_wf(output_dir, output_spaces, template, freesurfer,
         name='ds_mni_tpms', run_without_submitting=True)
     ds_mni_tpms.inputs.extra_values = ['CSF', 'GM', 'WM']
 
+    # Transforms
+    suffix_fmt = 'target-{}_{}'.format
     ds_t1_mni_warp = pe.Node(
         DerivativesDataSink(base_directory=output_dir, suffix=suffix_fmt(template, 'warp')),
         name='ds_t1_mni_warp', run_without_submitting=True)
 
+    lta_2_itk = pe.Node(fs.utils.LTAConvert(out_itk=True, invert=True), name='lta_2_itk')
+
+    ds_t1_fsnative = pe.Node(
+        DerivativesDataSink(base_directory=output_dir, suffix=suffix_fmt('fsnative', 'affine')),
+        name='ds_t1_fsnative', run_without_submitting=True)
+
     name_surfs = pe.MapNode(GiftiNameSource(pattern=r'(?P<LR>[lr])h.(?P<surf>.+)_converted.gii',
                                             template='{surf}.{LR}.surf'),
                             iterfield='in_file',
@@ -940,6 +952,9 @@ def init_anat_derivatives_wf(output_dir, output_spaces, template, freesurfer,
 
     if freesurfer:
         workflow.connect([
+            (inputnode, lta_2_itk, [('t1_2_fsnative_reverse_transform', 'in_lta')]),
+            (inputnode, ds_t1_fsnative, [('source_file', 'source_file')]),
+            (lta_2_itk, ds_t1_fsnative, [('out_itk', 'in_file')]),
             (inputnode, name_surfs, [('surfaces', 'in_file')]),
             (inputnode, ds_surfs, [('source_file', 'source_file'),
                                    ('surfaces', 'in_file')]),