From a783effa5e2196b4c4bbc6eef054b486aef1c5bc Mon Sep 17 00:00:00 2001 From: Matt Cieslak Date: Mon, 4 Mar 2024 15:40:11 -0500 Subject: [PATCH 1/5] make a json for the merged dwis --- qsiprep/interfaces/dwi_merge.py | 24 +++++++++++++++++++++++- qsiprep/interfaces/eddy.py | 24 ++++++++++++++++++++++++ qsiprep/workflows/dwi/base.py | 1 + qsiprep/workflows/dwi/fsl.py | 7 ++++++- qsiprep/workflows/dwi/hmc_sdc.py | 1 + qsiprep/workflows/dwi/merge.py | 6 ++++++ qsiprep/workflows/dwi/pre_hmc.py | 7 +++++++ 7 files changed, 68 insertions(+), 2 deletions(-) diff --git a/qsiprep/interfaces/dwi_merge.py b/qsiprep/interfaces/dwi_merge.py index 2f89d5d6..45b6a39c 100644 --- a/qsiprep/interfaces/dwi_merge.py +++ b/qsiprep/interfaces/dwi_merge.py @@ -48,6 +48,7 @@ class MergeDWIsInputSpec(BaseInterfaceInputSpec): carpetplot_data = InputMultiObject( File(exists=True), mandatory=False, desc="list of carpetplot_data files" ) + scan_metadata = traits.Dict(desc="Dict of metadata for the to-be-combined scans") class MergeDWIsOutputSpec(TraitedSpec): @@ -55,7 +56,7 @@ class MergeDWIsOutputSpec(TraitedSpec): out_bval = File(desc="the merged bval file") out_bvec = File(desc="the merged bvec file") original_images = traits.List() - merged_metadata = traits.Dict() + merged_metadata = File(exists=True) merged_denoising_confounds = File(exists=True) merged_b0_ref = File(exists=True) merged_raw_dwi = File(exists=True, mandatory=False) @@ -79,6 +80,17 @@ def _run_interface(self, runtime): self.inputs.harmonize_b0_intensities, ) + # Create a merged metadata json file for + if isdefined(self.inputs.scan_metadata): + combined_metadata = combine_metadata( + self.inputs.dwi_files, + self.inputs.scan_metadata, + ) + merged_metadata_file = op.join(runtime.cwd, "merged_metadata.json") + with open(merged_metadata_file, "w") as merged_meta_f: + json.dump(combined_metadata, merged_meta_f, indent=4) + self._results["merged_metadata"] = merged_metadata_file + # Get basic qc / provenance per volume provenance_df = create_provenance_dataframe( self.inputs.bids_dwi_files, to_concat, b0_means, corrections @@ -151,6 +163,16 @@ def _run_interface(self, runtime): return runtime +def combine_metadata(scan_list, metadata_dict): + """Create a merged metadata dictionary. + + Most importantly, combine the slice timings in some way. Here the slice + metadata is taken from the first scan. + + """ + return metadata_dict[scan_list[0]] + + class AveragePEPairsInputSpec(MergeDWIsInputSpec): original_bvec_files = InputMultiObject( File(exists=True), mandatory=True, desc="list of original bvec files" diff --git a/qsiprep/interfaces/eddy.py b/qsiprep/interfaces/eddy.py index 644c5434..b5b14364 100644 --- a/qsiprep/interfaces/eddy.py +++ b/qsiprep/interfaces/eddy.py @@ -9,6 +9,7 @@ """ import os import os.path as op +import json import nibabel as nb import numpy as np @@ -45,6 +46,8 @@ class GatherEddyInputsInputSpec(BaseInterfaceInputSpec): topup_max_b0s_per_spec = traits.CInt(1, usedefault=True) topup_requested = traits.Bool(False, usedefault=True) raw_image_sdc = traits.Bool(True, usedefault=True) + eddy_config = File(exists=True, mandatory=True) + json_file = File(exists=True) class GatherEddyInputsOutputSpec(TraitedSpec): @@ -60,6 +63,8 @@ class GatherEddyInputsOutputSpec(TraitedSpec): forward_transforms = traits.List() forward_warps = traits.List() topup_report = traits.Str(desc="description of where data came from") + json_file = File(exists=True) + multiband_factor = traits.Int() class GatherEddyInputs(SimpleInterface): @@ -125,6 +130,21 @@ def _run_interface(self, runtime): # these have already had HMC, SDC applied self._results["forward_transforms"] = [] self._results["forward_warps"] = [] + + # Based on the eddy config, determine whether to send a json argument + with open(self.inputs.eddy_config, "r") as eddy_cfg_f: + eddy_config = json.load(eddy_cfg_f) + # json file is only allowed if mporder is defined + if "mporder" in eddy_config: + self._results["json_file"] = self.inputs.json_file + + # Get the multiband factor from the json file and convert it to an arg + if isdefined(self.inputs.json_file): + with open(self.inputs.json_file) as json_f: + input_sidecar = json.load(json_f) + if "MultibandAccelerationFactor" in input_sidecar: + self._results["multiband_factor"] = input_sidecar["MultibandAccelerationFactor"] + return runtime @@ -241,6 +261,10 @@ def _format_arg(self, name, spec, value): if name == "field": pth, fname, _ = split_filename(value) return spec.argstr % op.join(pth, fname) + if name == "json": + if isdefined(self.inputs.mporder): + return spec.argstr % value + return "" return super(ExtendedEddy, self)._format_arg(name, spec, value) diff --git a/qsiprep/workflows/dwi/base.py b/qsiprep/workflows/dwi/base.py index e08461b4..bff65cdf 100644 --- a/qsiprep/workflows/dwi/base.py +++ b/qsiprep/workflows/dwi/base.py @@ -426,6 +426,7 @@ def init_dwi_preproc_wf( ('outputnode.dwi_file', 'inputnode.dwi_file'), ('outputnode.bval_file', 'inputnode.bval_file'), ('outputnode.bvec_file', 'inputnode.bvec_file'), + ('outputnode.json_file', 'inputnode.json_file'), ('outputnode.original_files', 'inputnode.original_files')]), (inputnode, hmc_wf, [ ('t1_brain', 'inputnode.t1_brain'), diff --git a/qsiprep/workflows/dwi/fsl.py b/qsiprep/workflows/dwi/fsl.py index adbf9e82..ec7752fd 100644 --- a/qsiprep/workflows/dwi/fsl.py +++ b/qsiprep/workflows/dwi/fsl.py @@ -91,6 +91,8 @@ def init_fsl_hmc_wf( bvec file bval_file: str bval file + json_file: str + path to sidecar json file for dwi_file b0_indices: list Indexes into ``dwi_files`` that correspond to b=0 volumes b0_images: list @@ -116,6 +118,7 @@ def init_fsl_hmc_wf( "dwi_file", "bvec_file", "bval_file", + "json_file", "b0_indices", "b0_images", "original_files", @@ -224,7 +227,9 @@ def init_fsl_hmc_wf( ('outputnode.ref_image_brain', 'dwi_file')]), (gather_inputs, eddy, [ ('eddy_indices', 'in_index'), - ('eddy_acqp', 'in_acqp')]), + ('eddy_acqp', 'in_acqp'), + ('json_file', 'json'), + ('multiband_factor', 'multiband_factor')]), (inputnode, eddy, [ ('dwi_file', 'in_file'), ('bval_file', 'in_bval'), diff --git a/qsiprep/workflows/dwi/hmc_sdc.py b/qsiprep/workflows/dwi/hmc_sdc.py index c97e9ccc..90795da0 100644 --- a/qsiprep/workflows/dwi/hmc_sdc.py +++ b/qsiprep/workflows/dwi/hmc_sdc.py @@ -82,6 +82,7 @@ def init_qsiprep_hmcsdc_wf( "dwi_file", "bvec_file", "bval_file", + "json_file", "rpe_b0", "t2w_unfatsat", "original_files", diff --git a/qsiprep/workflows/dwi/merge.py b/qsiprep/workflows/dwi/merge.py index 4285e607..a19867ff 100644 --- a/qsiprep/workflows/dwi/merge.py +++ b/qsiprep/workflows/dwi/merge.py @@ -18,6 +18,7 @@ from ...engine import Workflow from ...interfaces import ConformDwi, DerivativesDataSink +from ...interfaces.bids import get_metadata_for_nifti from ...interfaces.dipy import Patch2Self from ...interfaces.dwi_merge import MergeDWIs, PhaseToRad, StackConfounds from ...interfaces.gradients import ExtractB0s @@ -117,6 +118,8 @@ def init_merge_and_denoise_wf( bvals from merged images merged_bvec bvecs from merged images + merged_json + JSON file containing slice timings for slice2vol noise_image image(s) created by ``dwidenoise`` original_files @@ -133,6 +136,7 @@ def init_merge_and_denoise_wf( "merged_raw_image", "merged_bval", "merged_bvec", + "merged_json", "noise_images", "bias_images", "denoising_confounds", @@ -151,6 +155,7 @@ def init_merge_and_denoise_wf( bids_dwi_files=raw_dwi_files, b0_threshold=b0_threshold, harmonize_b0_intensities=not no_b0_harmonization, + scan_metadata={scan: get_metadata_for_nifti(scan) for scan in raw_dwi_files}, ), name="merge_dwis", n_procs=omp_nthreads, @@ -285,6 +290,7 @@ def init_merge_and_denoise_wf( ('original_images', 'original_files'), ('out_bval', 'merged_bval'), ('out_bvec', 'merged_bvec'), + ('merged_metadata', 'merged_json') ]), ]) # fmt:skip diff --git a/qsiprep/workflows/dwi/pre_hmc.py b/qsiprep/workflows/dwi/pre_hmc.py index 074fa439..14a1b83f 100644 --- a/qsiprep/workflows/dwi/pre_hmc.py +++ b/qsiprep/workflows/dwi/pre_hmc.py @@ -103,6 +103,8 @@ def init_dwi_pre_hmc_wf( a bvec file bval_file a bval files + sidecar_file + a json sidecar file for the scan data b0_indices list of the positions of the b0 images in the dwi series b0_images @@ -121,6 +123,7 @@ def init_dwi_pre_hmc_wf( "dwi_file", "bval_file", "bvec_file", + "json_file", "original_files", "denoising_confounds", "noise_images", @@ -296,6 +299,9 @@ def init_dwi_pre_hmc_wf( (pm_raw_images, raw_rpe_concat, [('out', 'in_files')]), (raw_rpe_concat, outputnode, [('out_file', 'raw_concatenated')]), + # Send the slice timings from "plus" to the next steps + (merge_plus, outputnode, [('outputnode.merged_json', 'json_file')]), + # Connect to the QC calculator (raw_rpe_concat, qc_wf, [('out_file', 'inputnode.dwi_file')]), (rpe_concat, qc_wf, [ @@ -333,6 +339,7 @@ def init_dwi_pre_hmc_wf( ('outputnode.merged_image', 'dwi_file'), ('outputnode.merged_bval', 'bval_file'), ('outputnode.merged_bvec', 'bvec_file'), + ('outputnode.merged_json', 'json_file'), ('outputnode.bias_images', 'bias_images'), ('outputnode.noise_images', 'noise_images'), ('outputnode.validation_reports', 'validation_reports'), From 42f6409075832b3b235218a0fdc7087308fc5e76 Mon Sep 17 00:00:00 2001 From: Matt Cieslak Date: Mon, 4 Mar 2024 20:34:28 -0500 Subject: [PATCH 2/5] now a unicode error? --- qsiprep/interfaces/dwi_merge.py | 2 +- qsiprep/workflows/dwi/fsl.py | 11 +++++++---- 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/qsiprep/interfaces/dwi_merge.py b/qsiprep/interfaces/dwi_merge.py index 45b6a39c..1aed654a 100644 --- a/qsiprep/interfaces/dwi_merge.py +++ b/qsiprep/interfaces/dwi_merge.py @@ -83,7 +83,7 @@ def _run_interface(self, runtime): # Create a merged metadata json file for if isdefined(self.inputs.scan_metadata): combined_metadata = combine_metadata( - self.inputs.dwi_files, + self.inputs.bids_dwi_files, self.inputs.scan_metadata, ) merged_metadata_file = op.join(runtime.cwd, "merged_metadata.json") diff --git a/qsiprep/workflows/dwi/fsl.py b/qsiprep/workflows/dwi/fsl.py index ec7752fd..4382f0fb 100644 --- a/qsiprep/workflows/dwi/fsl.py +++ b/qsiprep/workflows/dwi/fsl.py @@ -168,10 +168,6 @@ def init_fsl_hmc_wf( ) workflow = Workflow(name=name) - gather_inputs = pe.Node( - GatherEddyInputs(b0_threshold=b0_threshold, raw_image_sdc=raw_image_sdc), - name="gather_inputs", - ) if eddy_config is None: # load from the defaults eddy_cfg_file = pkgr_fn("qsiprep.data", "eddy_params.json") @@ -180,6 +176,13 @@ def init_fsl_hmc_wf( with open(eddy_cfg_file, "r") as f: eddy_args = json.load(f) + + gather_inputs = pe.Node( + GatherEddyInputs( + b0_threshold=b0_threshold, raw_image_sdc=raw_image_sdc, eddy_config=eddy_cfg_file + ), + name="gather_inputs", + ) enhance_pre_sdc = pe.Node(EnhanceB0(), name="enhance_pre_sdc") # Run in parallel if possible From 5b98c727d68b1981ca6e60062c72523fac1ce25c Mon Sep 17 00:00:00 2001 From: Matt Cieslak Date: Tue, 5 Mar 2024 10:24:33 -0500 Subject: [PATCH 3/5] actually connect the json --- qsiprep/workflows/dwi/fsl.py | 1 + 1 file changed, 1 insertion(+) diff --git a/qsiprep/workflows/dwi/fsl.py b/qsiprep/workflows/dwi/fsl.py index 4382f0fb..d5cb49e1 100644 --- a/qsiprep/workflows/dwi/fsl.py +++ b/qsiprep/workflows/dwi/fsl.py @@ -216,6 +216,7 @@ def init_fsl_hmc_wf( ('dwi_file', 'dwi_file'), ('bval_file', 'bval_file'), ('bvec_file', 'bvec_file'), + ('json_file', 'json_file'), ('original_files', 'original_files')]), (inputnode, pre_eddy_b0_ref_wf, [ ('t1_brain', 'inputnode.t1_brain'), From ca992724218cd50dcc636ea064d187a4a643068f Mon Sep 17 00:00:00 2001 From: Matt Cieslak Date: Tue, 5 Mar 2024 14:25:10 -0500 Subject: [PATCH 4/5] [FIX] Prevent undecodable terminal output in calc_connectivity (#711) closes #664 --- qsiprep/interfaces/dsi_studio.py | 23 ++++++++++++++--------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/qsiprep/interfaces/dsi_studio.py b/qsiprep/interfaces/dsi_studio.py index 1693ebdc..c3484478 100644 --- a/qsiprep/interfaces/dsi_studio.py +++ b/qsiprep/interfaces/dsi_studio.py @@ -361,6 +361,7 @@ class DSIStudioConnectivityMatrix(CommandLine): input_spec = DSIStudioConnectivityMatrixInputSpec output_spec = DSIStudioConnectivityMatrixOutputSpec _cmd = "dsi_studio --action=ana " + _terminal_output = "file" def _post_run_hook(self, runtime): atlas_config = self.inputs.atlas_config @@ -448,18 +449,22 @@ def _run_interface(self, runtime): workflow.config["execution"]["stop_on_first_crash"] = "true" workflow.config["execution"]["remove_unnecessary_outputs"] = "false" workflow.base_dir = runtime.cwd + plugin_settings = {} if num_threads > 1: - plugin_settings = { - "plugin": "MultiProc", - "plugin_args": { - "raise_insufficient": False, - "maxtasksperchild": 1, - "n_procs": num_threads, - }, + plugin_settings["plugin"] = "MultiProc" + plugin_settings["plugin_args"] = { + "raise_insufficient": False, + "maxtasksperchild": 1, + "n_procs": num_threads, } - wf_result = workflow.run(**plugin_settings) else: - wf_result = workflow.run() + plugin_settings["plugin"] = "Linear" + + workflow.config["execution"] = { + "stop_on_first_crash": "True", + "remove_unnecessary_outputs": "False", + } + wf_result = workflow.run(**plugin_settings) (merge_node,) = [ node for node in list(wf_result.nodes) if node.name.endswith("merge_mats") ] From 2d14337f5ef29e6ad2fd2d497352157d7af1e924 Mon Sep 17 00:00:00 2001 From: Matt Cieslak Date: Tue, 5 Mar 2024 14:29:43 -0500 Subject: [PATCH 5/5] address review --- qsiprep/interfaces/dwi_merge.py | 27 ++++++++++++++++++++++----- qsiprep/interfaces/eddy.py | 9 +-------- 2 files changed, 23 insertions(+), 13 deletions(-) diff --git a/qsiprep/interfaces/dwi_merge.py b/qsiprep/interfaces/dwi_merge.py index 1aed654a..fe1b7e5f 100644 --- a/qsiprep/interfaces/dwi_merge.py +++ b/qsiprep/interfaces/dwi_merge.py @@ -88,7 +88,7 @@ def _run_interface(self, runtime): ) merged_metadata_file = op.join(runtime.cwd, "merged_metadata.json") with open(merged_metadata_file, "w") as merged_meta_f: - json.dump(combined_metadata, merged_meta_f, indent=4) + json.dump(combined_metadata, merged_meta_f, sort_keys=True, indent=4) self._results["merged_metadata"] = merged_metadata_file # Get basic qc / provenance per volume @@ -163,14 +163,31 @@ def _run_interface(self, runtime): return runtime -def combine_metadata(scan_list, metadata_dict): +def combine_metadata(scan_list, metadata_dict, merge_method="first"): """Create a merged metadata dictionary. - Most importantly, combine the slice timings in some way. Here the slice - metadata is taken from the first scan. + Most importantly, combine the slice timings in some way. + + Parameters + ---------- + scan_list: list + List of BIDS inputs in the order in which they'll be concatenated + medadata_dict: dict + Mapping keys (values in ``scan_list``) to BIDS metadata dictionaries + merge_method: str + How to combine the metadata when multiple scans are being concatenated. + If "first" the metadata from the first scan is selected. Someday other + methods like "average" may be added. + + Returns + ------- + metadata: dict + A BIDS metadata dictionary """ - return metadata_dict[scan_list[0]] + if merge_method == "first": + return metadata_dict[scan_list[0]] + raise NotImplementedError(f"merge_method '{merge_method}' is not implemented") class AveragePEPairsInputSpec(MergeDWIsInputSpec): diff --git a/qsiprep/interfaces/eddy.py b/qsiprep/interfaces/eddy.py index b5b14364..dd669a44 100644 --- a/qsiprep/interfaces/eddy.py +++ b/qsiprep/interfaces/eddy.py @@ -7,9 +7,9 @@ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ """ +import json import os import os.path as op -import json import nibabel as nb import numpy as np @@ -138,13 +138,6 @@ def _run_interface(self, runtime): if "mporder" in eddy_config: self._results["json_file"] = self.inputs.json_file - # Get the multiband factor from the json file and convert it to an arg - if isdefined(self.inputs.json_file): - with open(self.inputs.json_file) as json_f: - input_sidecar = json.load(json_f) - if "MultibandAccelerationFactor" in input_sidecar: - self._results["multiband_factor"] = input_sidecar["MultibandAccelerationFactor"] - return runtime