Read in TotalReadoutTime from field map metadata, if it exists. Pass TotalReadoutTime to FSL TOPUP if present.

Author: sgiavasis

CPAC/distortion_correction/distortion_correction.py

@@ -486,10 +486,12 @@ def distcor_blip_fsl_topup(wf, cfg, strat_pool, pipe_num, opt=None):
                 "epi-1-pedir",
                 "epi-1-TE",
                 "epi-1-dwell",
+                "epi-1-total-readout",
                 "epi-2",
                 "epi-2-pedir",
                 "epi-2-TE",
-                "epi-2-dwell"],
+                "epi-2-dwell",
+                "epi-2-total-readout"],
      "outputs": ["desc-reginput_bold",
                  "space-bold_desc-brain_mask",
                  "blip-warp"]}
@@ -588,7 +590,9 @@ def distcor_blip_fsl_topup(wf, cfg, strat_pool, pipe_num, opt=None):
                 "phase_one",
                 "phase_two",
                 "dwell_time_one",
-                "dwell_time_two"
+                "dwell_time_two",
+                "ro_time_one",
+                "ro_time_two"
             ],
             output_names=["acq_params"],
             function=phase_encode,
@@ -604,12 +608,22 @@ def distcor_blip_fsl_topup(wf, cfg, strat_pool, pipe_num, opt=None):
 
     node, out = strat_pool.get_data('pe-direction')
     wf.connect(node, out, phase_encoding, 'unwarp_dir')
-        
-    node, out = strat_pool.get_data('epi-1-dwell')
-    wf.connect(node, out, phase_encoding, 'dwell_time_one')
 
-    node, out = strat_pool.get_data('epi-2-dwell')
-    wf.connect(node, out, phase_encoding, 'dwell_time_two')
+    if strat_pool.check_rpool('epi-1-dwell'):
+        node, out = strat_pool.get_data('epi-1-dwell')
+        wf.connect(node, out, phase_encoding, 'dwell_time_one')
+    
+    if strat_pool.check_rpool('epi-2-dwell'):
+        node, out = strat_pool.get_data('epi-2-dwell')
+        wf.connect(node, out, phase_encoding, 'dwell_time_two')
+
+    if strat_pool.check_rpool('epi-1-total-readout'):
+        node, out = strat_pool.get_data('epi-1-total-readout')
+        wf.connect(node, out, phase_encoding, 'ro_time_one')
+    
+    if strat_pool.check_rpool('epi-2-total-readout'):
+        node, out = strat_pool.get_data('epi-2-total-readout')
+        wf.connect(node, out, phase_encoding, 'ro_time_two')
 
     topup_imports = ["import os",
                      "import subprocess"]

CPAC/distortion_correction/utils.py

@@ -140,11 +140,9 @@ def gradient_distortion_correction(wf, inp_image, name):
     return (wf, out_warpmask, out_applywarp)
 
 
-def phase_encode(unwarp_dir, phase_one, phase_two, dwell_time_one, 
-                 dwell_time_two):
-    """
-
-    Calculate readout time and populate parameter file
+def phase_encode(unwarp_dir, phase_one, phase_two, dwell_time_one=None, 
+                 dwell_time_two=None, ro_time_one=None, ro_time_two=None):
+    """Calculate readout time and populate parameter file
 
     Parameters
     __________
@@ -159,8 +157,10 @@ def phase_encode(unwarp_dir, phase_one, phase_two, dwell_time_one,
         echo spacing of phase one
     dwell_time_two
         echo spacing of phase two
-    fsl_dir
-        FSL directory
+    ro_time_one
+        total readout time of phase one
+    ro_time_two
+        total readout time of phase two
 
     Returns
     _______
@@ -177,22 +177,30 @@ def phase_encode(unwarp_dir, phase_one, phase_two, dwell_time_one,
         unwarp_dir = unwarp_dir.decode()
 
     if unwarp_dir in ["x", "x-", "-x","i","-i","i-"]:
-        dim = nibabel.load(phase_one).shape[0]
-        n_PE_steps = dim - 1
-        ro_time_one = np.round(dwell_time_one * n_PE_steps, 6)
-        ro_time_two = np.round(dwell_time_two * n_PE_steps, 6)
+        if not ro_time_one or not ro_time_two:  
+            if dwell_time_one and dwell_time_two:
+                dim = nibabel.load(phase_one).shape[0]
+                n_PE_steps = dim - 1
+                ro_time_one = np.round(dwell_time_one * n_PE_steps, 6)
+                ro_time_two = np.round(dwell_time_two * n_PE_steps, 6)
+            else:
+                raise Exception("[!] No dwell time or total readout time "
+                                "present for the acq-fMRI EPI field maps.")
         ro_times = [f"-1 0 0 {ro_time_one}", f"1 0 0 {ro_time_two}"]
     elif unwarp_dir in ["y", "y-", "-y","j","-j","j-"]:
-        dim = nibabel.load(phase_one).shape[1]
-        n_PE_steps = dim - 1
-        ro_time_one = np.round(dwell_time_one * n_PE_steps, 6)
-        ro_time_two = np.round(dwell_time_two * n_PE_steps, 6)
+        if not ro_time_one or not ro_time_two:
+            if dwell_time_one and dwell_time_two:
+                dim = nibabel.load(phase_one).shape[1]
+                n_PE_steps = dim - 1
+                ro_time_one = np.round(dwell_time_one * n_PE_steps, 6)
+                ro_time_two = np.round(dwell_time_two * n_PE_steps, 6)
+            else:
+                raise Exception("[!] No dwell time or total readout time "
+                                "present for the acq-fMRI EPI field maps.")
         ro_times = [f"0 -1 0 {ro_time_one}", f"0 1 0 {ro_time_two}"]
     else:
         raise Exception(f"unwarp_dir={unwarp_dir} is unsupported.")
 
-
-
     # get number of volumes
     dims = [
         int(subprocess.check_output([f"fslval", phase_one, "dim4"]).decode(sys.stdout.encoding)),

CPAC/utils/bids_utils.py

@@ -702,7 +702,7 @@ def collect_bids_files_configs(bids_dir, aws_input_creds=''):
     file_paths = []
     config_dict = {}
 
-    suffixes = ['T1w', 'T2w', 'bold', 'acq-fMRI_epi', 'phasediff', 'magnitude',
+    suffixes = ['T1w', 'T2w', 'bold', 'epi', 'phasediff', 'magnitude',
                 'magnitude1', 'magnitude2']
 
     if bids_dir.lower().startswith("s3://"):

CPAC/utils/build_data_config.py

@@ -564,11 +564,11 @@ def get_BIDS_data_dct(bids_base_dir, file_list=None, anat_scan=None,
     fmap_pedir_sess = os.path.join(bids_base_dir,
                                    "sub-{participant}/ses-{session}/fmap/"
                                    "sub-{participant}_ses-{session}/"
-                                   "dir-*_acq-fMRI_epi.nii.gz")
+                                   "*acq-fMR*_epi.nii.gz")
 
     fmap_pedir = os.path.join(bids_base_dir,
                               "sub-{participant}/fmap/sub-{participant}"
-                              "_dir-*_acq-fMRI_epi.nii.gz")
+                              "*acq-fMR*_epi.nii.gz")
 
     sess_glob = os.path.join(bids_base_dir, "sub-*/ses-*/*")
 
@@ -582,7 +582,7 @@ def get_BIDS_data_dct(bids_base_dir, file_list=None, anat_scan=None,
 
     fmap_pedir_scan_glob = os.path.join(bids_base_dir,
                                         "sub-*fmap/"
-                                        "sub-*_dir-*_acq-fMRI_epi.nii.gz")
+                                        "sub-*_*acq-fMR*_epi.nii.gz")
 
     part_tsv_glob = os.path.join(bids_base_dir, "*participants.tsv")
 
@@ -648,7 +648,7 @@ def get_BIDS_data_dct(bids_base_dir, file_list=None, anat_scan=None,
                 fmap_mag = os.path.join(bids_base_dir,
                                         "sub-{participant}/fmap/sub-{participant}"
                                         "*magnitud*.nii.gz")
-
+            
             '''
             if fnmatch.fnmatch(filepath, fmap_pedir_scan_glob):
                 # check if there is a scan level for the fmap magnitude files

CPAC/utils/datasource.py

@@ -335,8 +335,9 @@ def get_fmap_phasediff_metadata(data_config_scan_params):
     echo_time = data_config_scan_params.get("EchoTime")
     dwell_time = data_config_scan_params.get("DwellTime")
     pe_direction = data_config_scan_params.get("PhaseEncodingDirection")
+    total_readout = data_config_scan_params.get("TotalReadoutTime")
 
-    return (echo_time, dwell_time, pe_direction)
+    return (echo_time, dwell_time, pe_direction, total_readout)
 
 
 def calc_deltaTE_and_asym_ratio(dwell_time, echo_time_one, echo_time_two,
@@ -445,7 +446,8 @@ def ingress_func_metadata(wf, cfg, rpool, sub_dict, subject_id,
                 input_names=['data_config_scan_params'],
                 output_names=['echo_time',
                               'dwell_time',
-                              'pe_direction'],
+                              'pe_direction',
+                              'total_readout'],
                 function=get_fmap_phasediff_metadata,
                 imports=get_fmap_metadata_imports),
                 name=f'{key}_get_metadata')
@@ -459,6 +461,8 @@ def ingress_func_metadata(wf, cfg, rpool, sub_dict, subject_id,
                            'dwell_time', {}, "", "fmap_dwell_ingress")
             rpool.set_data(f'{key}-pedir', get_fmap_metadata,
                            'pe_direction', {}, "", "fmap_pedir_ingress")
+            rpool.set_data(f'{key}-total-readout', get_fmap_metadata,
+                           'total_readout', {}, "", "fmap_readout_ingress")
 
             fmap_TE_list.append(f"{key}-TE")