Merge pull request #1925 from FCP-INDI/develop

Update main with recent changes to develop (FS ingress fix)

Author: sgiavasis

.github/workflows/smoke_test_participant.yml

@@ -122,14 +122,6 @@ jobs:
           - ''
         participant:
           - 032102 032106 032164 032167 032130 032128 2215 2312 032191 032195
-        include:
-          # requires T2w
-          - preconfig: monkey-ABCD
-            variant: ''
-            participant: 032102 032106 032130 032128 032191 032195
-          - preconfig: monkey-ABCD
-            variant: ABCD-HCP
-            participant: 032102 032106 032130 032128 032191 032195
     steps:
       - name: Get C-PAC
         run: |

CHANGELOG.md

@@ -30,6 +30,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Added NodeBlock information to `pypeline.log` when verbose debugging is on
 - Added the ability to ingress FreeSurfer data into CPAC
 - Added the ability to toggle FreeSurfer derived masks for brain extraction
+- Added an optional volume center to FD-J calculation
+- Added new preconfig `abcd-prep`, which performs minimal preprocessing on the T1w data in preparation for Freesurfer Recon-All
 
 ### Changed
 - Freesurfer output directory ingress moved to the data configuration YAML

CPAC/anat_preproc/anat_preproc.py

@@ -16,7 +16,8 @@
     mri_convert, \
     wb_command, \
     fslmaths_command, \
-    VolumeRemoveIslands
+    VolumeRemoveIslands, \
+    normalize_wmparc
 from CPAC.pipeline.engine import flatten_list
 from CPAC.utils.docs import docstring_parameter
 from CPAC.utils.interfaces.fsl import Merge as fslMerge
@@ -914,16 +915,44 @@ def freesurfer_abcd_brain_connector(wf, cfg, strat_pool, pipe_num, opt):
     ### ABCD harmonization - anatomical brain mask generation ###
     # Ref: https://github.com/DCAN-Labs/DCAN-HCP/blob/master/PostFreeSurfer/PostFreeSurferPipeline.sh#L151-L156
 
+    
     wmparc_to_nifti = pe.Node(util.Function(input_names=['in_file',
                                                          'reslice_like',
                                                          'args'],
                                             output_names=['out_file'],
                                             function=mri_convert),
                               name=f'wmparc_to_nifti_{pipe_num}')
-    wmparc_to_nifti.inputs.args = '-rt nearest'
+    
+    # Register wmparc file if ingressing FreeSurfer data
+    if strat_pool.check_rpool('pipeline-fs_xfm'):
+
+        wmparc_to_native = pe.Node(util.Function(input_names=['source_file',
+                                                            'target_file',
+                                                            'xfm',
+                                                            'out_file'],
+                                                output_names=['transformed_file'],
+                                                function=normalize_wmparc),
+                                        name=f'wmparc_to_native_{pipe_num}')
+        
+        wmparc_to_native.inputs.out_file = 'wmparc_warped.mgz'
+
+        node, out = strat_pool.get_data('pipeline-fs_wmparc')
+        wf.connect(node, out, wmparc_to_native, 'source_file')
+
+        node, out = strat_pool.get_data('pipeline-fs_raw-average')
+        wf.connect(node, out, wmparc_to_native, 'target_file')
 
-    node, out = strat_pool.get_data('pipeline-fs_wmparc')
-    wf.connect(node, out, wmparc_to_nifti, 'in_file')
+        node, out = strat_pool.get_data('pipeline-fs_xfm')
+        wf.connect(node, out, wmparc_to_native, 'xfm')
+
+        wf.connect(wmparc_to_native, 'transformed_file', wmparc_to_nifti, 'in_file')
+    
+    else:
+        
+        node, out = strat_pool.get_data('pipeline-fs_wmparc')
+        wf.connect(node, out, wmparc_to_nifti, 'in_file')
+
+    wmparc_to_nifti.inputs.args = '-rt nearest'
 
     node, out = strat_pool.get_data('desc-preproc_T1w')
     wf.connect(node, out, wmparc_to_nifti, 'reslice_like')
@@ -1244,21 +1273,24 @@ def acpc_align_head(wf, cfg, strat_pool, pipe_num, opt=None):
                 ["anatomical_preproc", "run"]],
      "option_key": "None",
      "option_val": "None",
-     "inputs": ["desc-head_T1w",
-                "desc-preproc_T1w",
+     "inputs": ["desc-preproc_T1w",
+                "desc-head_T1w",
                 "T1w-ACPC-template"],
      "outputs": ["desc-head_T1w",
                  "desc-preproc_T1w",
                  "from-T1w_to-ACPC_mode-image_desc-aff2rig_xfm"]}
     '''
-
+    
+    # NOTE: desc-preproc_T1w and desc-head_T1w are both needed in this 
+    # nodeblock despite the list called in get_data() below
+    
     acpc_align = acpc_alignment(config=cfg,
                                 acpc_target=cfg.anatomical_preproc[
                                     'acpc_alignment']['acpc_target'],
                                 mask=False,
                                 wf_name=f'acpc_align_{pipe_num}')
 
-    node, out = strat_pool.get_data(['desc-head_T1w', 'desc-preproc_T1w'])
+    node, out = strat_pool.get_data(['desc-preproc_T1w','desc-head_T1w'])
     wf.connect(node, out, acpc_align, 'inputspec.anat_leaf')
 
     node, out = strat_pool.get_data('T1w-ACPC-template')
@@ -1772,6 +1804,8 @@ def brain_mask_freesurfer_abcd(wf, cfg, strat_pool, pipe_num, opt=None):
      "option_val": "FreeSurfer-ABCD",
      "inputs": ["desc-preproc_T1w",
                 "pipeline-fs_wmparc",
+                "pipeline-fs_raw-average",
+                "pipeline-fs_xfm",
                 "freesurfer-subject-dir"],
      "outputs": ["space-T1w_desc-brain_mask"]}
     '''
@@ -1813,6 +1847,8 @@ def brain_mask_acpc_freesurfer_abcd(wf, cfg, strat_pool, pipe_num, opt=None):
      "option_val": "FreeSurfer-ABCD",
      "inputs": ["desc-preproc_T1w",
                 "pipeline-fs_wmparc",
+                "pipeline-fs_raw-average",
+                "pipeline-fs_xfm",
                 "freesurfer-subject-dir"],
      "outputs": ["space-T1w_desc-acpcbrain_mask"]}
     '''
@@ -2533,10 +2569,10 @@ def brain_extraction_temp_T2(wf, cfg, strat_pool, pipe_num, opt=None):
 def freesurfer_abcd_preproc(wf, cfg, strat_pool, pipe_num, opt=None):
     '''
     {"name": "freesurfer_abcd_preproc",
-     "config": ["anatomical_preproc", "brain_extraction"],
-     "switch": "None",
-     "option_key": "using",
-     "option_val": "FreeSurfer-ABCD",
+     "config": ["surface_analysis", "abcd_prefreesurfer_prep"],
+     "switch": ["run"],
+     "option_key": "None",
+     "option_val": "None",
      "inputs": ["desc-preproc_T1w",
                 "T1w-template",
                 "T1w-brain-template-mask",
@@ -2545,16 +2581,17 @@ def freesurfer_abcd_preproc(wf, cfg, strat_pool, pipe_num, opt=None):
                 "freesurfer-subject-dir"],
      "outputs": ["desc-restore_T1w",
                  "desc-restore-brain_T1w",
+                 "desc-ABCDpreproc_T1w",
                  "pipeline-fs_desc-fast_biasfield",
                  "pipeline-fs_hemi-L_desc-surface_curv",
-                 "pipeline-fs_hemi-R_desc-surface_curv", 
+                 "pipeline-fs_hemi-R_desc-surface_curv",
                  "pipeline-fs_hemi-L_desc-surfaceMesh_pial",
                  "pipeline-fs_hemi-R_desc-surfaceMesh_pial",
                  "pipeline-fs_hemi-L_desc-surfaceMesh_smoothwm",
                  "pipeline-fs_hemi-R_desc-surfaceMesh_smoothwm",
                  "pipeline-fs_hemi-L_desc-surfaceMesh_sphere",
                  "pipeline-fs_hemi-R_desc-surfaceMesh_sphere",
-                 "pipeline-fs_hemi-L_desc-surfaceMap_sulc", 
+                 "pipeline-fs_hemi-L_desc-surfaceMap_sulc",
                  "pipeline-fs_hemi-R_desc-surfaceMap_sulc",
                  "pipeline-fs_hemi-L_desc-surfaceMap_thickness",
                  "pipeline-fs_hemi-R_desc-surfaceMap_thickness",
@@ -2565,7 +2602,6 @@ def freesurfer_abcd_preproc(wf, cfg, strat_pool, pipe_num, opt=None):
                  "pipeline-fs_wmparc",
                  "freesurfer-subject-dir"]}
     '''
-
     # fnirt-based brain extraction
     brain_extraction = fnirt_based_brain_extraction(config=cfg,
                                                     wf_name=f'fnirt_based_brain_extraction_{pipe_num}')
@@ -2657,52 +2693,15 @@ def freesurfer_abcd_preproc(wf, cfg, strat_pool, pipe_num, opt=None):
     wf.connect(average_brain, 'out_stat',
                normalize_head, 'number')
 
-    if strat_pool.check_rpool('freesurfer-subject-dir'):
-        outputs = {
+    outputs = {
             'desc-restore_T1w': (fast_correction, 'outputspec.anat_restore'),
             'desc-restore-brain_T1w': (fast_correction,
                                        'outputspec.anat_brain_restore'),
-            'pipeline-fs_desc-fast_biasfield': (fast_correction, 'outputspec.bias_field')}
-        return (wf, outputs)
-
-    else:
-        ### recon-all -all step ###
-        reconall = pe.Node(interface=freesurfer.ReconAll(),
-                           name=f'anat_freesurfer_{pipe_num}',
-                           mem_gb=2.7)
-
-        sub_dir = cfg.pipeline_setup['working_directory']['path']
-        freesurfer_subject_dir = os.path.join(sub_dir,
-                                              'cpac_'+cfg['subject_id'],
-                                              f'anat_preproc_freesurfer_{pipe_num}',
-                                              'anat_freesurfer')
-
-        # create the directory for FreeSurfer node
-        if not os.path.exists(freesurfer_subject_dir):
-            os.makedirs(freesurfer_subject_dir)
-
-        reconall.inputs.directive = 'all'
-        reconall.inputs.subjects_dir = freesurfer_subject_dir
-        reconall.inputs.openmp = cfg.pipeline_setup['system_config']['num_OMP_threads']
-
-        wf.connect(normalize_head, 'out_file',
-                   reconall, 'T1_files')
-
-        wf, hemisphere_outputs = freesurfer_hemispheres(wf, reconall, pipe_num)
-
-        outputs = {
-        'desc-restore_T1w': (fast_correction, 'outputspec.anat_restore'),
-        'desc-restore-brain_T1w': (fast_correction,
-                                   'outputspec.anat_brain_restore'),
-        'pipeline-fs_desc-fast_biasfield': (fast_correction, 'outputspec.bias_field'),
-        'pipeline-fs_wmparc': (reconall, 'wmparc'),
-        'freesurfer-subject-dir': (reconall, 'subjects_dir'),
-        **hemisphere_outputs
-    }
-
+            'pipeline-fs_desc-fast_biasfield': (fast_correction, 'outputspec.bias_field'),
+            'desc-ABCDpreproc_T1w': (normalize_head, 'out_file')
+            }
     return (wf, outputs)
 
-
 # we're grabbing the postproc outputs and appending them to
 # the reconall outputs
 @docstring_parameter(postproc_outputs=str(flatten_list(freesurfer_abcd_preproc, 'outputs')
@@ -2715,7 +2714,8 @@ def freesurfer_reconall(wf, cfg, strat_pool, pipe_num, opt=None):
       "switch": ["run_reconall"],
       "option_key": "None",
       "option_val": "None",
-      "inputs": ["desc-preproc_T1w"],
+      "inputs": [["desc-ABCDpreproc_T1w",
+                 "desc-preproc_T1w"]],
       "outputs": ["freesurfer-subject-dir",
                   "pipeline-fs_raw-average",
                   "pipeline-fs_subcortical-seg",
@@ -2748,7 +2748,7 @@ def freesurfer_reconall(wf, cfg, strat_pool, pipe_num, opt=None):
         reconall.inputs.args = cfg.surface_analysis['freesurfer'][
             'reconall_args']
 
-    node, out = strat_pool.get_data("desc-preproc_T1w")
+    node, out = strat_pool.get_data(["desc-ABCDpreproc_T1w","desc-preproc_T1w"])
     wf.connect(node, out, reconall, 'T1_files')
 
     wf, hemisphere_outputs = freesurfer_hemispheres(wf, reconall, pipe_num)
@@ -2983,9 +2983,9 @@ def fast_bias_field_correction(config=None, wf_name='fast_bias_field_correction'
 def correct_restore_brain_intensity_abcd(wf, cfg, strat_pool, pipe_num, opt=None):
     '''
     {"name": "correct_restore_brain_intensity_abcd",
-     "config": ["anatomical_preproc", "brain_extraction"],
-     "switch": "None",
-     "option_key": "using",
+     "config": "None",
+     "switch": ["anatomical_preproc", "brain_extraction", "run"],
+     "option_key": ["anatomical_preproc", "brain_extraction", "using"],
      "option_val": "FreeSurfer-ABCD",
      "inputs": [("desc-preproc_T1w",
                  "desc-n4_T1w",

CPAC/anat_preproc/utils.py

@@ -420,6 +420,15 @@ def fslmaths_command(in_file, number, out_file_suffix):
 
     return out_file
 
+def normalize_wmparc(source_file, target_file, xfm, out_file):
+    from CPAC.utils.monitoring.custom_logging import log_subprocess
+    import os
+
+    cmd = ['mri_vol2vol', '--mov', source_file, \
+                '--targ', target_file, '--o', out_file, '--lta', xfm]
+    log_subprocess(cmd)
+    output = os.path.join(os.getcwd(), out_file)
+    return output
 
 """This module provides interfaces for workbench -volume-remove-islands commands"""
 from nipype.interfaces.base import TraitedSpec, File, traits, CommandLineInputSpec

CPAC/generate_motion_statistics/generate_motion_statistics.py

@@ -324,21 +324,31 @@ def calculate_FD_P(in_file):
     return out_file
 
 
-def calculate_FD_J(in_file, motion_correct_tool='3dvolreg'):
+def calculate_FD_J(in_file, motion_correct_tool='3dvolreg', center=None):
     """
     Method to calculate framewise displacement as per Jenkinson et al. 2002
 
     Parameters
     ----------
     in_file : string
-        matrix transformations from volume alignment file path
+        matrix transformations from volume alignment file path if
+        motion_correct_tool is '3dvolreg', or FDRMS (*_rel.rms) output if
+        motion_correct_tool is 'mcflirt'.
+    motion_correct_tool : string
+        motion correction tool used, '3dvolreg' or 'mcflirt'.
+    center : ndarray
+        optional volume center for the calculation.
 
     Returns
     -------
     out_file : string
         Frame-wise displacement file path
 
     """
+    if center is None:
+        center = np.zeros((3, 1))
+    else:
+        center = np.asarray(center).reshape((3, 1))
 
     if motion_correct_tool == '3dvolreg':
         pm_ = np.genfromtxt(in_file)
@@ -359,7 +369,7 @@ def calculate_FD_J(in_file, motion_correct_tool='3dvolreg'):
 
             M = np.dot(T_rb, np.linalg.inv(T_rb_prev)) - np.eye(4)
             A = M[0:3, 0:3]
-            b = M[0:3, 3]
+            b = M[0:3, 3:4] + A @ center
 
             fd[i] = np.sqrt(
                 (rmax * rmax / 5) * np.trace(np.dot(A.T, A)) + np.dot(b.T, b)
@@ -370,13 +380,39 @@ def calculate_FD_J(in_file, motion_correct_tool='3dvolreg'):
     elif motion_correct_tool == 'mcflirt':
         rel_rms = np.loadtxt(in_file)
         fd = np.append(0, rel_rms)
+    
+    else:
+        raise ValueError(f"motion_correct_tool {motion_correct_tool} not supported")
 
     out_file = os.path.join(os.getcwd(), 'FD_J.1D')
     np.savetxt(out_file, fd, fmt='%.8f')
 
     return out_file
 
 
+def find_volume_center(img_file):
+    """
+    Find the center of mass of a Nifti image volume
+
+    Parameters
+    ----------
+    img_file : string (nifti file)
+        path to nifti volume image
+
+    Returns
+    -------
+    center : ndarray 
+        volume center of mass vector
+    """
+    img = nb.load(img_file)
+    dim = np.array(img.header["dim"][1:4])
+    pixdim = np.array(img.header["pixdim"][1:4])
+    # Calculation follows MCFLIRT
+    # https://github.com/fithisux/FSL/blob/7aa2932949129f5c61af912ea677d4dbda843895/src/mcflirt/mcflirt.cc#L479
+    center = 0.5 * (dim - 1) * pixdim
+    return center
+
+
 def gen_motion_parameters(subject_id, scan_id, movement_parameters,
                           max_displacement, motion_correct_tool):
     """

CPAC/pipeline/__init__.py

@@ -25,7 +25,8 @@
                         x in ALL_PIPELINE_CONFIGS if 'pipeline_config' in x]
 ALL_PIPELINE_CONFIGS.sort()
 AVAILABLE_PIPELINE_CONFIGS = [preconfig for preconfig in ALL_PIPELINE_CONFIGS
-                              if preconfig != 'benchmark-ANTS' and
+                              if preconfig not in
+                              ['benchmark-ANTS', 'monkey-ABCD'] and
                               not preconfig.startswith('regtest-')]