Merge pull request #1769 from FCP-INDI/feature/total_readout_time

Ingress and use TotalReadoutTime from epi fmap meta-data in FSL topup

Author: sgiavasis

CHANGELOG.md

@@ -18,6 +18,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 - Added the ability to downsample to 10K or 2K resolution for freesurfer runs
+- Added the ability to ingress TotalReadoutTime from epi field map meta-data from the JSON sidecars.
+- Added the ability to use TotalReadoutTime of epi field maps in the calculation of FSL topup distortion correction.
 - Difference method (``-``) for ``CPAC.utils.configuration.Configuration`` instances
 
 ### Changed
@@ -35,6 +37,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Fixed
 - Fixed [bug](https://github.com/FCP-INDI/C-PAC/issues/1795) that was causing `cpac run` to fail when passing a manual random seed via `--random_seed`.
 - Replaces ``DwellTime`` with ``EffectiveEchoSpacing`` for FSL usage of the term
+- Fixed an issue that was causing some epi field maps to not be ingressed if the BIDS tags were not in the correct order.
 
 ## [v1.8.4] - 2022-06-27
 

CPAC/distortion_correction/distortion_correction.py

@@ -501,10 +501,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-mean_bold",
                  "space-bold_desc-brain_mask",
                  "blip-warp"]}
@@ -605,7 +607,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,
@@ -621,12 +625,24 @@ 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')
+    
+    if strat_pool.check_rpool('epi-1-dwell') and \
+            strat_pool.check_rpool('epi-2-dwell'):
 
-    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-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')
+        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') and \
+            strat_pool.check_rpool('epi-2-total-readout'):
+
+        node, out = strat_pool.get_data('epi-1-total-readout')
+        wf.connect(node, out, phase_encoding, 'ro_time_one')
+    
+        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
     _______
@@ -170,29 +170,44 @@ def phase_encode(unwarp_dir, phase_one, phase_two, dwell_time_one,
 
     """
 
+    meta_data = [dwell_time_one, dwell_time_two,
+                 ro_time_one, ro_time_two]
+    if not any(meta_data):
+        raise Exception("\n[!] Blip-FSL-TOPUP workflow: neither "
+                        "TotalReadoutTime nor DwellTime is present in the "
+                        "epi field map meta-data.")
+
     # create text file
     acq_params = os.path.join(os.getcwd(), "acqparams.txt")
 
     if isinstance(unwarp_dir, bytes):
         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)
-        ro_times = [f"-1 0 0 {ro_time_one}", f"1 0 0 {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)
+        elif ro_time_one and ro_time_two:
+            ro_times = [f"-1 0 0 {ro_time_one}", f"1 0 0 {ro_time_two}"]
+        else:
+            raise Exception("[!] No dwell time or total readout time "
+                            "present for the acq-fMRI EPI field maps.")
     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)
-        ro_times = [f"0 -1 0 {ro_time_one}", f"0 1 0 {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)
+        elif ro_time_one and ro_time_two:
+            ro_times = [f"0 -1 0 {ro_time_one}", f"0 1 0 {ro_time_two}"]
+        else:
+            raise Exception("[!] No dwell time or total readout time "
+                            "present for the acq-fMRI EPI field maps.")
     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/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

@@ -350,8 +350,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_delta_te_and_asym_ratio(effective_echo_spacing: float,
@@ -484,7 +485,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{name_suffix}')
@@ -498,6 +500,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")