[RTM] Replace FLIRT-BBR with FreeSurfer bbregister when FS preprocessing is enabled

docs/links.rst

@@ -3,3 +3,5 @@
 .. _Usage: usage.html
 .. _Installation: installation.html
 .. _workflows: workflows.html
+.. _FreeSurfer: https://surfer.nmr.mgh.harvard.edu/
+.. _`submillimeter reconstruction`: https://surfer.nmr.mgh.harvard.edu/fswiki/SubmillimeterRecon

docs/workflows.rst

@@ -12,6 +12,8 @@ and at least one functional (EPI) file, but no SBRefs or fieldmaps.
 To force using this pipeline on datasets that do include fieldmaps and SBRefs
 use the ``--ignore fieldmaps`` flag.
 
+Several steps are added or modified if `Surface preprocessing`_ is enabled.
+
 What It Does
 ------------
 High-level view of the basic pipeline:
@@ -49,14 +51,9 @@ warp to the MNI space.
 
     Animation showing T1 to MNI normalization (ANTs)
 
-Additionally, FreeSurfer surfaces are reconstructed from T1-weighted structural
+If enabled, FreeSurfer surfaces are reconstructed from T1-weighted structural
 image(s), using the ANTs-extracted brain mask.
-This feature may be disabled with the ``--no-freesurfer`` flag.
-
-.. figure:: _static/reconall.svg
-    :scale: 100%
-
-    Surface reconstruction (FreeSurfer)
+See Reconstruction_ for details.
 
 EPI_HMC
 ~~~~~~~
@@ -75,7 +72,7 @@ is used to perform skullstripping of the mean EPI image.
     Brain extraction (nilearn).
 
 ref_epi_t1_registration
-~~~~~~~~~~~~~~~~~~~~
+~~~~~~~~~~~~~~~~~~~~~~~
 
 .. image:: ref_epi_t1_registration.dot.png
     :scale: 100%
@@ -88,14 +85,17 @@ function to find the transform that maps the EPI space into the T1-space.
 
     Animation showing EPI to T1 registration (FSL FLIRT with BBR)
 
+If surface processing is enabled, ``bbregister`` is used instead.
+See `Boundary-based Registration (BBR)`_ for details.
+
 EPIMNITransformation
 ~~~~~~~~~~~~~~~~~~~~
 
 .. image:: EPIMNITransformation.dot.png
     :scale: 100%
 
 The EPIMNITransformation sub-workflow uses the transform from
-`EPIMeanNormalization`_ and a t1-to-MNI transform from `t1w_preprocessing`_ to
+`EPIMeanNormalization`_ and a T1-to-MNI transform from `t1w_preprocessing`_ to
 map the EPI image to standardized MNI space.
 It also maps the t1w-based mask to MNI space.
 
@@ -151,6 +151,82 @@ Derivatives related to EPI files are in the ``func`` subfolder:
 - ``*bold_space-T1w_preproc.nii.gz`` Motion-corrected (using MCFLIRT for estimation and ANTs for interpolation) EPI file in T1w space
 - ``*bold_space-MNI152NLin2009cAsym_preproc.nii.gz`` Same as above, but in MNI space
 
+
+Surface preprocessing
+=====================
+
+``fmriprep`` uses FreeSurfer_ to reconstruct surfaces from T1/T2-weighted
+structural images.
+If enabled, several steps in the ``fmriprep`` pipeline are added or replaced.
+All surface preprocessing may be disabled with the ``--no-freesurfer`` flag.
+
+Reconstruction
+--------------
 If FreeSurfer reconstruction is performed, the reconstructed subject is placed in
-``<output dir>/freesurfer/sub-<subject_label>/``.
+``<output dir>/freesurfer/sub-<subject_label>/`` (see `FreeSurfer Derivatives`_).
+
+Surface reconstruction is performed in three phases.
+The first phase initializes the subject with T1- and T2-weighted (if available)
+structural images and performs basic reconstruction (``autorecon1``) with the
+exception of skull-stripping.
+For example, a subject with only one session with T1 and T2-weighted images
+would be processed by the following command::
+
+    $ recon-all -sd <output dir>/freesurfer -subjid sub-<subject_label> \
+        -i <bids-root>/sub-<subject_label>/anat/sub-<subject_label>_T1w.nii.gz \
+        -T2 <bids-root>/sub-<subject_label>/anat/sub-<subject_label>_T2w.nii.gz \
+        -autorecon1 \
+        -noskullstrip
+
+The second phase imports the brainmask calculated in the t1w_preprocessing_
+sub-workflow.
+The final phase resumes reconstruction, using the T2-weighted image to assist
+in finding the pial surface, if available::
+
+    $ recon-all -sd <output dir>/freesurfer -subjid sub-<subject_label> \
+        -all -T2pial
+
+Reconstructed white and pial surfaces are included in the report.
+
+.. figure:: _static/reconall.svg
+    :scale: 100%
+
+    Surface reconstruction (FreeSurfer)
+
+If T1-weighted voxel sizes are less 1mm in all dimensions (rounding to nearest
+.1mm), `submillimeter reconstruction`_ is used.
+
+In order to bypass reconstruction in ``fmriprep``, place existing reconstructed
+subjects in ``<output dir>/freesurfer`` prior to the run.
+``fmriprep`` will perform any missing ``recon-all`` steps, but will not perform
+any steps whose outputs already exist.
+
+Boundary-based Registration (BBR)
+---------------------------------
+The mean EPI image of each run is aligned to the reconstructed subject using
+the gray/white matter boundary (FreeSurfer's ``?h.white`` surfaces).
+
+If FreeSurfer processing is disabled, FLIRT is performed with the BBR cost
+function, using the FAST segmentation to establish the gray/white matter
+boundary.
+
+FreeSurfer Derivatives
+----------------------
+
+A FreeSurfer subjects directory is created in ``<output dir>/freesurfer``.
+
+::
+
+    freesurfer/
+        fsaverage/
+            mri/
+            surf/
+            ...
+        sub-<subject_label>/
+            mri/
+            surf/
+            ...
+        ...
 
+A copy of the ``fsaverage`` subject distributed with the running version of
+FreeSurfer is copied into this subjects directory.

fmriprep/viz/config.json

@@ -52,7 +52,13 @@
                 "name": "epi_mean_t1_registration/flt_bbr",
                 "file_pattern": "func/.*_bold_flt_bbr",
                 "title": "EPI",
-                "description": "No SBRef was found so the mean epi was used to generate transformations from EPI space to T1 Space"
+                "description": "No SBRef was found so the mean EPI was used to generate transformations from EPI space to T1 Space - FAST segmentation used for BBR"
+            },
+            {
+                "name": "epi_mean_t1_registration/bbr",
+                "file_pattern": "func/.*_bold_bbr",
+                "title": "EPI",
+                "description": "No SBRef was found so the mean EPI was used to generate transformations from EPI space to T1 Space"
             },
             {
                 "name": "epi/acompcor",
@@ -76,6 +82,12 @@
                 "name": "sbref_t1_registration/flt_bbr",
                 "file_pattern": "func/.*sbref.*flt_bbr",
                 "title": "SBRef to T1 registration",
+                "description": "FreeSurfer surfaces not found - FAST segmentation used for BBR"
+            },
+            {
+                "name": "sbref_t1_registration/bbr",
+                "file_pattern": "func/.*sbref.*bbr",
+                "title": "SBRef to T1 registration",
                 "description": " "
             },
             {

fmriprep/workflows/anatomical.py

@@ -39,7 +39,8 @@ def t1w_preprocessing(name='t1w_preprocessing', settings=None):
     outputnode = pe.Node(niu.IdentityInterface(
         fields=['t1_seg', 't1_tpms', 'bias_corrected_t1', 't1_brain', 't1_mask',
                 't1_2_mni', 't1_2_mni_forward_transform',
-                't1_2_mni_reverse_transform']), name='outputnode')
+                't1_2_mni_reverse_transform', 'subject_id',
+                'fs_2_t1_transform']), name='outputnode')
 
     # 0. Align and merge if several T1w images are provided
     t1wmrg = pe.Node(IntraModalMerge(), name='MergeT1s')
@@ -180,6 +181,11 @@ def inject_skullstripped(subjects_dir, subject_id, skullstripped):
             name='Reconstruction2')
         reconall.interface.num_threads = nthreads
 
+        fs_transform = pe.Node(
+            freesurfer.utils.Tkregister2(fsl_out='freesurfer2subT1.mat',
+                                         reg_header=True),
+            name='FreeSurferTransform')
+
         recon_report = pe.Node(
             DerivativesDataSink(base_directory=settings['reportlets_dir'],
                                 suffix='reconall'),
@@ -275,9 +281,13 @@ def inject_skullstripped(subjects_dir, subject_id, skullstripped):
             (injector, reconall, [('subjects_dir', 'subjects_dir'),
                                   ('subject_id', 'subject_id')]),
             (recon_config, reconall, [('use_T2', 'use_T2')]),
-            (inputnode, recon_report, [
+            (t1wmrg, fs_transform, [('out_avg', 'target_image')]),
+            (autorecon1, fs_transform, [('T1', 'moving_image')]),
+            (recon_config, recon_report, [
                 (('t1w', fix_multi_T1w_source_name), 'source_file')]),
             (reconall, recon_report, [('out_report', 'in_file')]),
+            (reconall, outputnode, [('subject_id', 'subject_id')]),
+            (fs_transform, outputnode, [('fsl_file', 'fs_2_t1_transform')]),
             ])
 
     # Write corrected file in the designated output dir

fmriprep/workflows/base.py

@@ -108,7 +108,7 @@ def basic_fmap_sbref_wf(subject_data, settings, name='fMRI_prep'):
     sbref_pre = sbref_preprocess(settings=settings)
 
     # Register SBRef to T1
-    sbref_t1 = ref_epi_t1_registration(reportlet_suffix='sbref_t1_flt_bbr',
+    sbref_t1 = ref_epi_t1_registration(reportlet_suffix='sbref_t1_bbr',
                                        inv_ds_suffix='target-sbref_affine',
                                        settings=settings)
 
@@ -175,6 +175,9 @@ def basic_fmap_sbref_wf(subject_data, settings, name='fMRI_prep'):
     if settings['freesurfer']:
         workflow.connect([
             (inputnode, t1w_pre, [('subjects_dir', 'inputnode.subjects_dir')]),
+            (inputnode, sbref_t1, [('subjects_dir', 'inputnode.subjects_dir')]),
+            (t1w_pre, sbref_t1, [('outputnode.subject_id', 'inputnode.subject_id'),
+                                 ('outputnode.fs_2_t1_transform', 'inputnode.fs_2_t1_transform')]),
             ])
 
     return workflow
@@ -214,7 +217,7 @@ def basic_wf(subject_data, settings, name='fMRI_prep'):
     hmcwf.get_node('inputnode').iterables = ('epi', subject_data['func'])
 
     # mean EPI registration to T1w
-    epi_2_t1 = ref_epi_t1_registration(reportlet_suffix='flt_bbr',
+    epi_2_t1 = ref_epi_t1_registration(reportlet_suffix='bbr',
                                        inv_ds_suffix='target-meanBOLD_affine',
                                        settings=settings)
 
@@ -260,6 +263,9 @@ def basic_wf(subject_data, settings, name='fMRI_prep'):
     if settings['freesurfer']:
         workflow.connect([
             (inputnode, t1w_pre, [('subjects_dir', 'inputnode.subjects_dir')]),
+            (inputnode, epi_2_t1, [('subjects_dir', 'inputnode.subjects_dir')]),
+            (t1w_pre, epi_2_t1, [('outputnode.subject_id', 'inputnode.subject_id'),
+                                 ('outputnode.fs_2_t1_transform', 'inputnode.fs_2_t1_transform')]),
             ])
 
     return workflow

fmriprep/workflows/epi.py

@@ -15,9 +15,10 @@
 from nipype.interfaces import ants
 from nipype.interfaces import c3
 from nipype.interfaces import fsl
+from nipype.interfaces import freesurfer as fs
 from nipype.interfaces import utility as niu
 from niworkflows.interfaces.masks import ComputeEPIMask, BETRPT
-from niworkflows.interfaces.registration import FLIRTRPT
+from niworkflows.interfaces.registration import FLIRTRPT, BBRegisterRPT
 from niworkflows.data import get_mni_icbm152_nlin_asym_09c
 
 from fmriprep.interfaces import DerivativesDataSink, FormatHMCParam
@@ -93,7 +94,8 @@ def ref_epi_t1_registration(reportlet_suffix, inv_ds_suffix, name='ref_epi_t1_re
     inputnode = pe.Node(
         niu.IdentityInterface(fields=['name_source', 'ref_epi', 'ref_epi_mask',
                                       'bias_corrected_t1', 't1_brain', 't1_mask',
-                                      't1_seg', 't1w', 'epi_split', 'hmc_xforms']),
+                                      't1_seg', 't1w', 'epi_split', 'hmc_xforms',
+                                      'subjects_dir', 'subject_id', 'fs_2_t1_transform']),
         name='inputnode'
     )
     outputnode = pe.Node(
@@ -112,16 +114,47 @@ def ref_epi_t1_registration(reportlet_suffix, inv_ds_suffix, name='ref_epi_t1_re
 
     explicit_mask_epi = pe.Node(fsl.ApplyMask(), name="explicit_mask_epi")
 
-    flt_bbr_init = pe.Node(
-        FLIRTRPT(generate_report=True, dof=6),
-        name='flt_bbr_init'
-    )
-    flt_bbr = pe.Node(
-        FLIRTRPT(generate_report=True, dof=6, cost_func='bbr'),
-        name='flt_bbr'
-    )
-    flt_bbr.inputs.schedule = op.join(os.getenv('FSLDIR'),
-                                      'etc/flirtsch/bbr.sch')
+    if settings['freesurfer']:
+        bbregister = pe.Node(
+            BBRegisterRPT(
+                contrast_type='t2',
+                init='fsl',
+                registered_file=True,
+                out_fsl_file=True,
+                generate_report=True),
+            name='bbregister'
+            )
+
+        def apply_fs_transform(fs_2_t1_transform, bbreg_transform):
+            import os
+            import numpy as np
+            out_file = os.path.abspath('transform.mat')
+            fs_xfm = np.loadtxt(fs_2_t1_transform)
+            bbrxfm = np.loadtxt(bbreg_transform)
+            out_xfm = fs_xfm.dot(bbrxfm)
+            assert np.allclose(out_xfm[3], [0, 0, 0, 1])
+            out_xfm[3] = [0, 0, 0, 1]
+            np.savetxt(out_file, out_xfm, fmt='%.12g')
+            return out_file
+
+        transformer = pe.Node(
+            niu.Function(
+                function=apply_fs_transform,
+                input_names=['fs_2_t1_transform', 'bbreg_transform'],
+                output_names=['out_file']),
+            name='BBRegTransform')
+    else:
+        flt_bbr_init = pe.Node(
+            FLIRTRPT(generate_report=True, dof=6),
+            name='flt_bbr_init'
+        )
+        flt_bbr = pe.Node(
+            FLIRTRPT(generate_report=True, dof=6, cost_func='bbr'),
+            name='flt_bbr'
+        )
+        flt_bbr.inputs.schedule = op.join(os.getenv('FSLDIR'),
+                                          'etc/flirtsch/bbr.sch')
+        reportlet_suffix = reportlet_suffix.replace('bbr', 'flt_bbr')
 
     # make equivalent warp fields
     invt_bbr = pe.Node(fsl.ConvertXFM(invert_xfm=True), name='Flirt_BBR_Inv')
@@ -144,24 +177,14 @@ def ref_epi_t1_registration(reportlet_suffix, inv_ds_suffix, name='ref_epi_t1_re
         (inputnode, explicit_mask_epi, [('ref_epi', 'in_file'),
                                         ('ref_epi_mask', 'mask_file')
                                         ]),
-        (explicit_mask_epi, flt_bbr_init, [('out_file', 'in_file')]),
-        (inputnode, flt_bbr_init, [('t1_brain', 'reference')]),
         (inputnode, fsl2itk_fwd, [('bias_corrected_t1', 'reference_file'),
                                   ('ref_epi', 'source_file')]),
         (inputnode, fsl2itk_inv, [('ref_epi', 'reference_file'),
                                   ('bias_corrected_t1', 'source_file')]),
-        (flt_bbr_init, flt_bbr, [('out_matrix_file', 'in_matrix_file')]),
-        (inputnode, flt_bbr, [('t1_brain', 'reference')]),
-        (explicit_mask_epi, flt_bbr, [('out_file', 'in_file')]),
-        (wm_mask, flt_bbr, [('out_file', 'wm_seg')]),
-        (flt_bbr, invt_bbr, [('out_matrix_file', 'in_file')]),
         (invt_bbr, outputnode, [('out_file', 'mat_t1_to_epi')]),
-        (flt_bbr, outputnode, [('out_matrix_file', 'mat_epi_to_t1')]),
-        (flt_bbr, fsl2itk_fwd, [('out_matrix_file', 'transform_file')]),
         (invt_bbr, fsl2itk_inv, [('out_file', 'transform_file')]),
         (fsl2itk_fwd, outputnode, [('itk_transform', 'itk_epi_to_t1')]),
         (fsl2itk_inv, outputnode, [('itk_transform', 'itk_t1_to_epi')]),
-        (flt_bbr, ds_report, [('out_report', 'in_file')]),
         (inputnode, ds_report, [(('name_source', _first), 'source_file')])
     ])
 
@@ -227,6 +250,32 @@ def ref_epi_t1_registration(reportlet_suffix, inv_ds_suffix, name='ref_epi_t1_re
              [(('name_source', _first), 'source_file')]),
             (merge, ds_t1w, [('merged_file', 'in_file')]),
             (mask_t1w_tfm, ds_t1w_mask, [('output_image', 'in_file')]),
+            ])
+
+    if settings['freesurfer']:
+        workflow.connect([
+            (inputnode, bbregister, [('subjects_dir', 'subjects_dir'),
+                                     ('subject_id', 'subject_id')]),
+            (explicit_mask_epi, bbregister, [('out_file', 'source_file')]),
+            (inputnode, transformer, [('fs_2_t1_transform', 'fs_2_t1_transform')]),
+            (bbregister, transformer, [('out_fsl_file', 'bbreg_transform')]),
+            (transformer, invt_bbr, [('out_file', 'in_file')]),
+            (transformer, outputnode, [('out_file', 'mat_epi_to_t1')]),
+            (transformer, fsl2itk_fwd, [('out_file', 'transform_file')]),
+            (bbregister, ds_report, [('out_report', 'in_file')]),
+        ])
+    else:
+        workflow.connect([
+            (explicit_mask_epi, flt_bbr_init, [('out_file', 'in_file')]),
+            (inputnode, flt_bbr_init, [('t1_brain', 'reference')]),
+            (flt_bbr_init, flt_bbr, [('out_matrix_file', 'in_matrix_file')]),
+            (inputnode, flt_bbr, [('t1_brain', 'reference')]),
+            (explicit_mask_epi, flt_bbr, [('out_file', 'in_file')]),
+            (wm_mask, flt_bbr, [('out_file', 'wm_seg')]),
+            (flt_bbr, invt_bbr, [('out_matrix_file', 'in_file')]),
+            (flt_bbr, outputnode, [('out_matrix_file', 'mat_epi_to_t1')]),
+            (flt_bbr, fsl2itk_fwd, [('out_matrix_file', 'transform_file')]),
+            (flt_bbr, ds_report, [('out_report', 'in_file')]),
         ])
 
     return workflow