Merge remote-tracking branch 'origin/192-introduce-func2func-workflow…

…-for-functional-to-sbref-registration'

# Conflicts:
#	pipelines/rcpl.py

PUMI/pipelines/func/func2func.py

@@ -0,0 +1,101 @@
+from nipype.interfaces import fsl
+from PUMI.engine import AnatPipeline, NestedNode as Node, QcPipeline
+from PUMI.pipelines.multimodal.image_manipulation import pick_volume, vol2png
+
+
+@QcPipeline(inputspec_fields=['registered_func', 'func_2'],
+            outputspec_fields=['out_file'])
+def qc(wf, registered_func_volume='middle', func_2_volume='middle', **kwargs):
+    """
+    Creates quality check images for func2func workflow.
+
+    Parameters:
+        registered_func_volume: Select which volume from registered image to use as overlay
+                            (i.e., 'first', 'middle', 'last', 'mean' or an arbitrary volume).
+        func_2_volume: Select which volume from func_2 to use as background.
+                            (i.e., 'first', 'middle', 'last', 'mean' or an arbitrary volume).
+
+    Inputs:
+        registered_func (str): Path to registered functional image.
+        func_2 (str): Path to 4D reference image (e.g., sbref).
+    """
+
+    one_registered_func_volume = pick_volume('one_registered_func_volume', volume=registered_func_volume, sink=False)
+    wf.connect('inputspec', 'registered_func', one_registered_func_volume, 'in_file')
+
+    one_func_2_volume = pick_volume('one_func_2_volume', volume=func_2_volume, sink=False)
+    wf.connect('inputspec', 'func_2', one_func_2_volume, 'in_file')
+
+    qc_vol2png = vol2png(wf.name + '_vol2png')
+    wf.connect(one_registered_func_volume, 'out_file', qc_vol2png, 'overlay_image')
+    wf.connect(one_func_2_volume, 'out_file', qc_vol2png, 'bg_image')
+
+    # Output
+    wf.connect(qc_vol2png, 'out_file', 'outputspec', 'out_file')
+
+    # Sinking
+    wf.connect(qc_vol2png, 'out_file', 'sinker', wf.name)
+
+
+@AnatPipeline(inputspec_fields=['func_1', 'func_2'],
+              outputspec_fields=['out_file', 'func_to_func_2_matrix', 'inv_func_to_func_2_matrix'])
+def func2func(wf, func_1_volume='mean', func_2_volume='mean', **kwargs):
+    """
+    Register one functional image to another functional image.
+
+    Parameters:
+        func_1_volume (str): Select which volume from func_1 to use for calculation of transformation matrix
+                            (i.e., 'first', 'middle', 'last', 'mean' or an arbitrary volume).
+        func_2_volume (str): Select which volume from func_2 to use as reference for calculation of transformation matrix.
+                            (i.e., 'first', 'middle', 'last', 'mean' or an arbitrary volume).
+
+    Inputs:
+        func_1 (str): Path to 4D functional image (e.g., rsfMRI scan).
+        func_2 (str): Path to 4D reference image (e.g., sbref).
+
+    Returns:
+        out_file (str): Registered functional image.
+        func_to_func_2_matrix: Transformation matrix (func_1 to func_2).
+        inv_func_to_func_2_matrix: Inverse transformation matrix (func_2 to func_1).
+    """
+
+    # Select volumes from func_1
+    func_1_volume_selected = pick_volume('func_1_volume_selected', volume=func_1_volume)
+    wf.connect('inputspec', 'func_1', func_1_volume_selected, 'in_file')
+
+    # Select volumes from func_2
+    func_2_volume_selected = pick_volume('func_2_volume_selected', volume=func_2_volume)
+    wf.connect('inputspec', 'func_2', func_2_volume_selected, 'in_file')
+
+    # Linear registration (FLIRT) to register selected volumes from func_1 and func_2
+    linear_func_to_func_2 = Node(interface=fsl.FLIRT(), name='linear_func_to_func_2')
+    linear_func_to_func_2.inputs.cost = 'corratio'
+    linear_func_to_func_2.inputs.dof = 6
+    linear_func_to_func_2.inputs.out_matrix_file = "func_to_func_2.mat"
+    wf.connect(func_1_volume_selected, 'out_file', linear_func_to_func_2, 'in_file')
+    wf.connect(func_2_volume_selected, 'out_file', linear_func_to_func_2, 'reference')
+
+    # Calculate the inverse of the transformation matrix
+    invert_matrix = Node(interface=fsl.utils.ConvertXFM(), name='invert_matrix')
+    invert_matrix.inputs.invert_xfm = True
+    wf.connect(linear_func_to_func_2, 'out_matrix_file', invert_matrix, 'in_file')
+
+    # Apply transformation matrix to the entire 4D func_1 image
+    apply_xfm_4d = Node(interface=fsl.ApplyXFM(), name='apply_xfm_4d')
+    apply_xfm_4d.inputs.interp = 'trilinear'
+    wf.connect('inputspec', 'func_1', apply_xfm_4d, 'in_file')  # entire 4D func_1 image
+    wf.connect(linear_func_to_func_2, 'out_matrix_file', apply_xfm_4d, 'in_matrix_file')
+    wf.connect('inputspec', 'func_2', apply_xfm_4d, 'reference')  # entire 4D func_2 image
+
+    # QC
+    func2func_qc = qc('func2func_qc')
+    wf.connect(apply_xfm_4d, 'out_file', func2func_qc, 'registered_func')
+    wf.connect('inputspec', 'func_2', func2func_qc, 'func_2')
+
+    # Sinking
+    wf.connect(apply_xfm_4d, 'out_file', 'sinker', 'registered_func')
+
+    # Outputs
+    wf.connect(apply_xfm_4d, 'out_file', 'outputspec', 'out_file')
+    wf.connect(linear_func_to_func_2, 'out_matrix_file', 'outputspec', 'func_to_func_2_matrix')
+    wf.connect(invert_matrix, 'out_file', 'outputspec', 'inv_func_to_func_2_matrix')

pipelines/rcpl.py

@@ -8,7 +8,7 @@
 from PUMI.pipelines.func.compcor import anat_noise_roi, compcor
 from PUMI.pipelines.anat.func_to_anat import func2anat
 from nipype.interfaces import utility
-
+from PUMI.pipelines.func.func2func import func2func
 from PUMI.pipelines.func.deconfound import fieldmap_correction_topup
 from PUMI.pipelines.func.func_proc import func_proc_despike_afni
 from PUMI.pipelines.func.timeseries_extractor import pick_atlas, extract_timeseries_nativespace
@@ -392,6 +392,11 @@ def merge_predictions(rpn_out_file, rcpl_out_file):
         acquisition='bold',
         suffix='epi',
         extension='.json'
+    ),
+    'sbref': dict(
+        datatype='func',
+        suffix="sbref",
+        extension=['nii', 'nii.gz']
     )
 })
 def rcpl(wf, bbr=True, **kwargs):
@@ -407,17 +412,24 @@ def rcpl(wf, bbr=True, **kwargs):
     reorient_fmap_wf = Node(Reorient2Std(output_type='NIFTI_GZ'), name="reorient_fmap_wf")
     wf.connect('inputspec', 'fmap', reorient_fmap_wf, 'in_file')
 
+    reorient_sbref_wf = Node(Reorient2Std(output_type='NIFTI_GZ'), name="reorient_sbref_wf")
+    wf.connect('inputspec', 'sbref', reorient_sbref_wf, 'in_file')
+
+    anatomical_preprocessing_wf = anat_proc(name='anatomical_preprocessing_wf', bet_tool='deepbet')
+    wf.connect(reorient_struct_wf, 'out_file', anatomical_preprocessing_wf, 'in_file')
+
     fieldmap_corr = fieldmap_correction_topup('fieldmap_corr')
     wf.connect(reorient_func_wf, 'out_file', fieldmap_corr, 'main')
     wf.connect('inputspec', 'bold_json', fieldmap_corr, 'main_json')
     wf.connect(reorient_fmap_wf, 'out_file', fieldmap_corr, 'fmap')
     wf.connect('inputspec', 'fmap_json', fieldmap_corr, 'fmap_json')
 
-    anatomical_preprocessing_wf = anat_proc(name='anatomical_preprocessing_wf', bet_tool='deepbet')
-    wf.connect(reorient_struct_wf, 'out_file', anatomical_preprocessing_wf, 'in_file')
+    func2sbref = func2func('func2sbref')
+    wf.connect(fieldmap_corr, 'out_file', func2sbref, 'func_1')
+    wf.connect(reorient_sbref_wf, 'out_file', func2sbref, 'func_2')
 
     func2anat_wf = func2anat(name='func2anat_wf', bbr=bbr)
-    wf.connect(fieldmap_corr, 'out_file', func2anat_wf, 'func')
+    wf.connect(func2sbref, 'out_file', func2anat_wf, 'func')
     wf.connect(anatomical_preprocessing_wf, 'brain', func2anat_wf, 'head')
     wf.connect(anatomical_preprocessing_wf, 'probmap_wm', func2anat_wf, 'anat_wm_segmentation')
     wf.connect(anatomical_preprocessing_wf, 'probmap_csf', func2anat_wf, 'anat_csf_segmentation')
@@ -429,7 +441,7 @@ def rcpl(wf, bbr=True, **kwargs):
     wf.connect(func2anat_wf, 'ventricle_mask_in_funcspace', compcor_roi_wf, 'ventricle_mask')
 
     func_proc_wf = func_proc_despike_afni('func_proc_wf', bet_tool='deepbet', deepbet_n_dilate=2)
-    wf.connect(fieldmap_corr, 'out_file', func_proc_wf, 'func')
+    wf.connect(func2sbref, 'out_file', func_proc_wf, 'func')
     wf.connect(compcor_roi_wf, 'out_file', func_proc_wf, 'cc_noise_roi')
 
     pick_atlas_wf = mist_atlas('pick_atlas_wf')
@@ -463,11 +475,11 @@ def rcpl(wf, bbr=True, **kwargs):
 
     predict_pain_sensitivity_rpn_wf = predict_pain_sensitivity_rpn('predict_pain_sensitivity_rpn_wf')
     wf.connect(calculate_connectivity_wf, 'features', predict_pain_sensitivity_rpn_wf, 'X')
-    wf.connect(fieldmap_corr, 'out_file', predict_pain_sensitivity_rpn_wf, 'in_file')
+    wf.connect(func2sbref, 'out_file', predict_pain_sensitivity_rpn_wf, 'in_file')
 
     predict_pain_sensitivity_rcpl_wf = predict_pain_sensitivity_rcpl('predict_pain_sensitivity_rcpl_wf')
     wf.connect(calculate_connectivity_wf, 'features', predict_pain_sensitivity_rcpl_wf, 'X')
-    wf.connect(fieldmap_corr, 'out_file', predict_pain_sensitivity_rcpl_wf, 'in_file')
+    wf.connect(func2sbref, 'out_file', predict_pain_sensitivity_rcpl_wf, 'in_file')
 
     collect_pain_predictions_wf = collect_pain_predictions('collect_pain_predictions_wf')
     wf.connect(predict_pain_sensitivity_rpn_wf, 'out_file', collect_pain_predictions_wf, 'rpn_out_file')