From 32172f8d51b99d255ffa675cd9f2cccab9510258 Mon Sep 17 00:00:00 2001 From: Matt Cieslak Date: Wed, 4 Dec 2024 12:00:44 -0500 Subject: [PATCH 01/13] initial commit [skip ci] --- qsirecon/interfaces/gradients.py | 83 +++++++++++- qsirecon/workflows/recon/fod_autotrack.py | 156 ++++++++++++++++++++++ qsirecon/workflows/recon/utils.py | 33 ++++- 3 files changed, 267 insertions(+), 5 deletions(-) create mode 100644 qsirecon/workflows/recon/fod_autotrack.py diff --git a/qsirecon/interfaces/gradients.py b/qsirecon/interfaces/gradients.py index 2be923b2..bd4a34cd 100644 --- a/qsirecon/interfaces/gradients.py +++ b/qsirecon/interfaces/gradients.py @@ -7,6 +7,7 @@ from nilearn import image as nim from nipype.interfaces.base import ( BaseInterfaceInputSpec, + InputMultiObject, File, SimpleInterface, TraitedSpec, @@ -18,7 +19,81 @@ LOGGER = logging.getLogger("nipype.interface") -class RemoveDuplicatesInputSpec(BaseInterfaceInputSpec): +class _GradientSelectInputSpec(BaseInterfaceInputSpec): + dwi_file = File(exists=True, mandatory=True) + bval_file = File(exists=True, mandatory=True) + bvec_file = File(exists=True, mandatory=True) + bval_distance_cutoff = traits.Float(5.0, usedefault=True) + expected_n_input_shells = traits.Int() + expected_n_output_shells = traits.Int() + requested_shell_bvals = InputMultiObject(traits.CInt(), mandatory=True) + + +class _GradientSelectOutputSpec(TraitedSpec): + dwi_file = File(exists=True) + bval_file = File(exists=True) + bvec_file = File(exists=True) + local_bvec_file = File(exists=True) + + +class GradientSelect(SimpleInterface): + input_spec = _GradientSelectInputSpec + output_spec = _GradientSelectOutputSpec + + def _run_interface(self, runtime): + """Find shells in the input data and select""" + + from sklearn.cluster import AgglomerativeClustering + from sklearn.d + max_distance = self.inputs.bval_distance_cutoff + bvals = np.loadtxt(self.inputs.bval_file) + agg_cluster = AgglomerativeClustering( + n_clusters=None, + distance_threshold=2 * max_distance, + linkage="complete", + ).fit(bvals.reshape(-1, 1)) + + # Check that the correct number of input shells are detected + if isdefined(self.inputs.expected_n_input_shells): + if not self.inputs.expected_n_input_shells == agg_cluster.n_clusters_: + howmuch, newthresh = ( + ("too many", "higher") + if agg_cluster.n_clusters_ > self.inputs.expected_n_input_shells + else ("too few", "lower") + ) + raise Exception( + f"Expected to find {self.inputs.expected_n_input_shells} shells in " + f"the input data. Instead we found {agg_cluster.n_clusters_}. Having " + f"{howmuch} shells detected means you may need to adjust the" + f"bval_distance_cutoff parameter to be {newthresh}." + ) + + # Ensure that at lease 1 b>0 and one b=0 shell are selected + select_shell_bs = self.inputs.requested_shell_bvals + if len(select_shell_bs) < 2: + if select_shell_bs[0] < max_distance: + LOGGER.critical("No effectively b>0 shells are selected.") + else: + LOGGER.warning("Adding a b=0 shell to the selection.") + select_shell_bs.append(0) + + # Make sure the shells are unique + select_shell_bs = np.array(sorted(select_shell_bs)) + shell_distances = + + + + return runtime + + +def _select_shell_by_bval(desired_bval, distance_max, bvals, bval_assignments=None): + """Find indices where bvals are within distance_max of desired_bval. + + If bval_assignments is defined, + """ + + +class _RemoveDuplicatesInputSpec(BaseInterfaceInputSpec): dwi_file = File(exists=True, mandatory=True) bval_file = File(exists=True, mandatory=True) bvec_file = File(exists=True, mandatory=True) @@ -27,7 +102,7 @@ class RemoveDuplicatesInputSpec(BaseInterfaceInputSpec): expected_directions = traits.Int() -class RemoveDuplicatesOutputSpec(TraitedSpec): +class _RemoveDuplicatesOutputSpec(TraitedSpec): dwi_file = File(exists=True) bval_file = File(exists=True) bvec_file = File(exists=True) @@ -35,8 +110,8 @@ class RemoveDuplicatesOutputSpec(TraitedSpec): class RemoveDuplicates(SimpleInterface): - input_spec = RemoveDuplicatesInputSpec - output_spec = RemoveDuplicatesOutputSpec + input_spec = _RemoveDuplicatesInputSpec + output_spec = _RemoveDuplicatesOutputSpec def _run_interface(self, runtime): bvecs = np.loadtxt(self.inputs.bvec_file).T diff --git a/qsirecon/workflows/recon/fod_autotrack.py b/qsirecon/workflows/recon/fod_autotrack.py new file mode 100644 index 00000000..fea2b7ee --- /dev/null +++ b/qsirecon/workflows/recon/fod_autotrack.py @@ -0,0 +1,156 @@ +import logging +import nipype.pipeline.engine as pe +from nipype.interfaces import utility as niu + +from ... import config +from ...engine import Workflow +from ...interfaces.interchange import recon_workflow_input_fields + +from .dsi_studio import ( + init_dsi_studio_autotrack_wf, + init_dsi_studio_connectivity_wf, + init_dsi_studio_export_wf, + init_dsi_studio_recon_wf, + init_dsi_studio_tractography_wf, +) +from .mrtrix import ( + init_global_tractography_wf, + init_mrtrix_connectivity_wf, + init_mrtrix_csd_recon_wf, + init_mrtrix_tractography_wf, +) + +LOGGER = logging.getLogger("nipype.interface") + + +def init_singleshell_benchmarking_wf( + available_anatomical_data, name="_recon", qsirecon_suffix="SingleShellBenchmark", params={} +): + inputnode = pe.Node( + niu.IdentityInterface(fields=recon_workflow_input_fields + ["odf_rois"]), name="inputnode" + ) + outputnode = pe.Node( + niu.IdentityInterface(fields=["tck_files", "bundle_names", "recon_scalars"]), + name="outputnode", + ) + outputnode.inputs.recon_scalars = [] + workflow = Workflow(name=name) + omp_nthreads = config.nipype.omp_nthreads + + autotrack_params = params.get("autotrack", {}) + bundle_names = _get_dsi_studio_bundles(autotrack_params.get("track_id", "")) + bundle_desc = ( + "AutoTrack attempted to reconstruct the following bundles:\n * " + + "\n * ".join(bundle_names) + + "\n\n" + ) + LOGGER.info(bundle_desc) + + # First do a standard GQI reconstruction + gqi_params = params.get("gqi_recon", {}) + initial_gqi_wf = init_dsi_studio_recon_wf( + available_anatomical_data=available_anatomical_data, + name="initial_gqi", + qsirecon_suffix=f"{qsirecon_suffix}_part-GQI", + params=gqi_params, + ) + + # Do an SS3T recon to feed to autotrack + ss3t_params = params.get("ss3t_recon") + ss3t_wf = init_mrtrix_csd_recon_wf( + available_anatomical_data=available_anatomical_data, + name="ss3t_recon", + qsirecon_suffix=f"{qsirecon_suffix}_part-SS3T", + params=ss3t_params, + ) + ss3t_to_fib = pe.Node(FODtoFIBGZ(), name="ss3t_to_fib") + + # For comparison, also do a regular CSD + csd_params = params.get("csd_recon") + csd_wf = init_mrtrix_csd_recon_wf( + available_anatomical_data=available_anatomical_data, + name="csd_recon", + qsirecon_suffix=f"{qsirecon_suffix}_part-CSD", + params=csd_params, + ) + csd_to_fib = pe.Node(FODtoFIBGZ(), name="csd_to_fib") + + # Run autotrack! + gqi_autotrack = pe.Node( + AutoTrack(num_threads=omp_nthreads, **params), name="gqi_autotrack", n_procs=omp_nthreads + ) + ss3t_autotrack = pe.Node( + AutoTrack(num_threads=omp_nthreads, **params), name="ss3t_autotrack", n_procs=omp_nthreads + ) + csd_autotrack = pe.Node( + AutoTrack(num_threads=omp_nthreads, **params), name="csd_autotrack", n_procs=omp_nthreads + ) + + # Create a single output + aggregate_gqi_atk_results = pe.Node( + AggregateAutoTrackResults(expected_bundles=bundle_names), name="aggregate_gqi_atk_results" + ) + aggregate_ss3t_atk_results = pe.Node( + AggregateAutoTrackResults(expected_bundles=bundle_names), name="aggregate_ss3t_atk_results" + ) + aggregate_csd_atk_results = pe.Node( + AggregateAutoTrackResults(expected_bundles=bundle_names), name="aggregate_csd_atk_results" + ) + + convert_gqi_to_tck = pe.MapNode(DSIStudioTrkToTck(), name="convert_gqi_to_tck", iterfield="trk_file") + convert_ss3t_to_tck = pe.MapNode(DSIStudioTrkToTck(), name="convert_ss3t_to_tck", iterfield="trk_file") + convert_csd_to_tck = pe.MapNode(DSIStudioTrkToTck(), name="convert_csd_to_tck", iterfield="trk_file") + + # Save the bundle csv + ds_gqi_bundle_csv = pe.Node( + ReconDerivativesDataSink(suffix="bundlestats", qsirecon_suffix=f"{qsirecon_suffix}_part-GQI"), + name="ds_gqi_bundle_csv", + run_without_submitting=True, + ) + ds_ss3t_bundle_csv = pe.Node( + ReconDerivativesDataSink(suffix="bundlestats", qsirecon_suffix=f"{qsirecon_suffix}_part-SS3T"), + name="ds_ss3t_bundle_csv", + run_without_submitting=True, + ) + ds_csd_bundle_csv = pe.Node( + ReconDerivativesDataSink(suffix="bundlestats", qsirecon_suffix=f"{qsirecon_suffix}_part-CSD"), + name="ds_csd_bundle_csv", + run_without_submitting=True, + ) + + # Save the mapping file. We're only using the mapping from GQI + ds_mapping = pe.Node( + ReconDerivativesDataSink(suffix="mapping", qsirecon_suffix=qsirecon_suffix), + name="ds_mapping", + run_without_submitting=True, + ) + + workflow.connect([ + # Connect the qsiprep inputs to the recon workflows we're creating here + # (normally this is done in build_workflow()) + (inputnode, initial_gqi_wf, + _as_connections(recon_workflow_input_fields, dest_prefix='inputnode.')), + (inputnode, ss3t_wf, + _as_connections(recon_workflow_input_fields, dest_prefix='inputnode.')), + (inputnode, csd_wf, + _as_connections(recon_workflow_input_fields, dest_prefix='inputnode.')), + + # Convert the sh mifs from csd and ss3t into fib files + (csd_wf, ss3t_to_fib, [("outputnode.mif_file", "mif_file")]), + (initial_gqi_wf, csd_to_fib, [("outputnode.fibgz", "fib_file")]), + (ss3t_wf, ss3t_to_fib, [("outputnode.mif_file", "mif_file")]), + (initial_gqi_wf, ss3t_to_fib, [("outputnode.fibgz", "fib_file")]), + + # Send the fib files to autotrack. Use the map file from gqi in ss3t and csd + (initial_gqi_wf, gqi_autotrack, [("outputnode.fibgz", "fib_file")]), + (ss3t_to_fib, ss3t_autotrack, [("fib_file", "fib_file")]), + (gqi_autotrack, ss3t_autotrack, [("map_file", "map_file")]), + (csd_to_fib, csd_autotrack, [("fib_file", "fib_file")]), + (gqi_autotrack, csd_autotrack, [("map_file", "map_file")]), + ]) # fmt:skip + + return workflow + + +def _as_connections(attr_list, src_prefix="", dest_prefix=""): + return [(src_prefix + item, dest_prefix + item) for item in attr_list] \ No newline at end of file diff --git a/qsirecon/workflows/recon/utils.py b/qsirecon/workflows/recon/utils.py index 1de9feb4..aa1c3863 100644 --- a/qsirecon/workflows/recon/utils.py +++ b/qsirecon/workflows/recon/utils.py @@ -14,7 +14,7 @@ from ...interfaces import ConformDwi from ...interfaces.bids import DerivativesDataSink -from ...interfaces.gradients import RemoveDuplicates +from ...interfaces.gradients import RemoveDuplicates, GradientSelect from ...interfaces.interchange import recon_workflow_input_fields from ...interfaces.mrtrix import MRTrixGradientTable from ...interfaces.recon_scalars import OrganizeScalarData @@ -82,6 +82,37 @@ def init_discard_repeated_samples_wf( return workflow +def init_gradient_select_wf( + inputs_dict, + name="gradient_select_wf", + qsirecon_suffix="", + params={}, +): + """Remove a sample if a similar direction/gradient has already been sampled.""" + inputnode = pe.Node( + niu.IdentityInterface(fields=recon_workflow_input_fields), name="inputnode" + ) + outputnode = pe.Node( + niu.IdentityInterface(fields=["dwi_file", "bval_file", "bvec_file", "local_bvec_file"]), + name="outputnode", + ) + workflow = Workflow(name=name) + + gradient_select = pe.Node(GradientSelect(**params), name="gradient_select") + workflow.connect([ + (inputnode, gradient_select, [ + ('dwi_file', 'dwi_file'), + ('bval_file', 'bval_file'), + ('bvec_file', 'bvec_file')]), + (gradient_select, outputnode, [ + ('dwi_file', 'dwi_file'), + ('bval_file', 'bval_file'), + ('bvec_file', 'bvec_file')]) + ]) # fmt:skip + + return workflow + + def init_scalar_output_wf( name="scalar_output_wf", ): From 67a3b0adebf72607d4857658e23845b7e777441c Mon Sep 17 00:00:00 2001 From: Matt Cieslak Date: Wed, 4 Dec 2024 23:25:23 -0500 Subject: [PATCH 02/13] continued [skip ci] --- .../mrtrix_singleshell_ss3t_ACT-fast.yaml | 12 +- qsirecon/interfaces/gradients.py | 271 +++++++++++++----- qsirecon/workflows/recon/build_workflow.py | 9 +- qsirecon/workflows/recon/fod_autotrack.py | 28 +- qsirecon/workflows/recon/utils.py | 2 +- 5 files changed, 244 insertions(+), 78 deletions(-) diff --git a/qsirecon/data/pipelines/mrtrix_singleshell_ss3t_ACT-fast.yaml b/qsirecon/data/pipelines/mrtrix_singleshell_ss3t_ACT-fast.yaml index 2d6ffe4d..c8348905 100644 --- a/qsirecon/data/pipelines/mrtrix_singleshell_ss3t_ACT-fast.yaml +++ b/qsirecon/data/pipelines/mrtrix_singleshell_ss3t_ACT-fast.yaml @@ -2,8 +2,18 @@ anatomical: - mrtrix_5tt_fast name: mrtrix_singleshell_ss3t_fast nodes: -- action: csd + +- action: input: qsirecon + name: select_single_shell + parameters: + requested_shells: + - 0 + - highest + bval_distance_cutoff: 100 + +- action: csd + input: select_single_shell name: ss3t_csd parameters: fod: diff --git a/qsirecon/interfaces/gradients.py b/qsirecon/interfaces/gradients.py index bd4a34cd..2e9a6d3c 100644 --- a/qsirecon/interfaces/gradients.py +++ b/qsirecon/interfaces/gradients.py @@ -2,13 +2,12 @@ import logging -import nibabel as nb import numpy as np from nilearn import image as nim from nipype.interfaces.base import ( BaseInterfaceInputSpec, - InputMultiObject, File, + InputMultiObject, SimpleInterface, TraitedSpec, isdefined, @@ -23,16 +22,21 @@ class _GradientSelectInputSpec(BaseInterfaceInputSpec): dwi_file = File(exists=True, mandatory=True) bval_file = File(exists=True, mandatory=True) bvec_file = File(exists=True, mandatory=True) + b_file = File(exists=True, mandatory=True) + btable_file = File(exists=True, mandatory=True) bval_distance_cutoff = traits.Float(5.0, usedefault=True) expected_n_input_shells = traits.Int() - expected_n_output_shells = traits.Int() - requested_shell_bvals = InputMultiObject(traits.CInt(), mandatory=True) + requested_shells = InputMultiObject( + traits.Either(traits.CInt(), traits.Enum("highest", "lowest")), mandatory=True + ) class _GradientSelectOutputSpec(TraitedSpec): dwi_file = File(exists=True) bval_file = File(exists=True) bvec_file = File(exists=True) + b_file = File(exists=True) + btable_file = File(exists=True) local_bvec_file = File(exists=True) @@ -42,61 +46,150 @@ class GradientSelect(SimpleInterface): def _run_interface(self, runtime): """Find shells in the input data and select""" - - from sklearn.cluster import AgglomerativeClustering - from sklearn.d - max_distance = self.inputs.bval_distance_cutoff bvals = np.loadtxt(self.inputs.bval_file) - agg_cluster = AgglomerativeClustering( - n_clusters=None, - distance_threshold=2 * max_distance, - linkage="complete", - ).fit(bvals.reshape(-1, 1)) - - # Check that the correct number of input shells are detected - if isdefined(self.inputs.expected_n_input_shells): - if not self.inputs.expected_n_input_shells == agg_cluster.n_clusters_: - howmuch, newthresh = ( - ("too many", "higher") - if agg_cluster.n_clusters_ > self.inputs.expected_n_input_shells - else ("too few", "lower") + + selected_indices = _select_gradients( + requested_gradients=self.inputs.requested_shells, + max_distance=self.inputs.bval_distance_cutoff, + original_bvals=bvals, + expected_n_input_shells=self.inputs.expected_n_input_shells or None, + ) + + # Make a new dwi and associated gradient files + new_bval, new_bvec, new_b, new_btable, new_nifti = subset_dwi( + original_bval=self.inputs.bval_file, + original_bvec=self.inputs.bvec_file, + original_b=self.inputs.b_file or None, + original_btable=self.inputs.btable_file or None, + original_nifti=self.inputs.dwi_file, + indices=selected_indices, + newdir=runtime.cwd, + suffix="_selected", + ) + self._results["bval_file"] = new_bval + self._results["bvec_file"] = new_bvec + self._results["dwi_file"] = new_nifti + if new_b: + self._results["b_file"] = new_b + if new_btable: + self._results["btable_file"] = new_btable + + return runtime + + +def _select_gradients( + requested_gradients, + max_distance, + original_bvals, + expected_n_input_shells=None, +): + """Find indices where bvals are within distance_max of desired_bval.""" + from sklearn.metrics import pairwise_distances + + # Get the bvals clustered into assigned shells. Also in a dataframe + bval_df = _find_shells(original_bvals, max_distance) + + # Check that the correct number of input shells are detected + if expected_n_input_shells is not None: + n_found_shells = len(bval_df["assigned_shell"].unique()) + if not expected_n_input_shells == n_found_shells: + howmuch, newthresh = ( + ("too many", "higher") + if n_found_shells > expected_n_input_shells + else ("too few", "lower") + ) + raise Exception( + f"Expected to find {expected_n_input_shells} shells in " + f"the input data. Instead we found {n_found_shells}. Having " + f"{howmuch} shells detected means you may need to adjust the" + f"bval_distance_cutoff parameter to be {newthresh}." + ) + + # Ensure that at lease 1 b>0 and one b=0 shell are selected + select_shell_bs = _parse_shell_selection(requested_gradients, bval_df) + + # Make sure the shells are unique (no overlap when accounting for allowed distances) + shell_distances = pairwise_distances(select_shell_bs) + if np.any(shell_distances[np.triu_indices_from(shell_distances, k=1)] < max_distance): + raise Exception( + "Shells and bval_distance_cutoff have overlap. Choose a lower " + "bval_distance_cutoff or more separated shells." + ) + + selected_indices = np.flatnonzero(bval_df["assigned_shell"].isin(select_shell_bs)) + + return selected_indices + + +def _parse_shell_selection(requested_bvals, bval_df, max_distance): + """Turn bval requests into numbers. Add a 0 if none were originally included.""" + numeric_bvals = [] + for requested_bval in requested_bvals: + if requested_bval == "highest": + highest_shell = bval_df["assigned_shell"].max() + LOGGER.info(f"Selecting b={highest_shell} as the highest shell") + numeric_bvals.append(highest_shell) + elif requested_bval == "lowest": + lowest_shell = bval_df["assigned_shell"][ + bval_df["assigned_shell"] > max_distance + ].min() + LOGGER.info(f"Selecting b={lowest_shell} as the lowest b>0 shell") + numeric_bvals.append(lowest_shell) + else: + # Find the closest detected shell + ok_shells = bval_df["assigned_shell"][ + np.logical_and( + bval_df["assigned_shell"] >= requested_bval - max_distance, + bval_df["assigned_shell"] <= requested_bval + max_distance, ) + ] + if len(ok_shells.unique()) > 1: raise Exception( - f"Expected to find {self.inputs.expected_n_input_shells} shells in " - f"the input data. Instead we found {agg_cluster.n_clusters_}. Having " - f"{howmuch} shells detected means you may need to adjust the" - f"bval_distance_cutoff parameter to be {newthresh}." + f"Unable to unambiguously select b={requested_bval}." + f"instead select from {bval_df.assigned_shells.unique().tolist()}" ) - # Ensure that at lease 1 b>0 and one b=0 shell are selected - select_shell_bs = self.inputs.requested_shell_bvals - if len(select_shell_bs) < 2: - if select_shell_bs[0] < max_distance: - LOGGER.critical("No effectively b>0 shells are selected.") - else: - LOGGER.warning("Adding a b=0 shell to the selection.") - select_shell_bs.append(0) + numeric_bvals.append(ok_shells.iloc[0]) - # Make sure the shells are unique - select_shell_bs = np.array(sorted(select_shell_bs)) - shell_distances = + # Make sure there is a 0 in the list + if len(numeric_bvals) < 2: + if numeric_bvals[0] < max_distance: + LOGGER.critical("No effectively b>0 shells are selected.") + else: + LOGGER.warning("Adding a b=0 shell to the selection.") + numeric_bvals.append(0) + return sorted(numeric_bvals) - return runtime +def _find_shells(bvals, max_distance): + + import pandas as pd + from sklearn.cluster import AgglomerativeClustering + agg_cluster = AgglomerativeClustering( + n_clusters=None, + distance_threshold=2 * max_distance, + linkage="complete", + ).fit(bvals.reshape(-1, 1)) -def _select_shell_by_bval(desired_bval, distance_max, bvals, bval_assignments=None): - """Find indices where bvals are within distance_max of desired_bval. + bval_df = pd.DataFrame({"bvalue": bvals, "assignment": agg_cluster.labels_}) + shell_df = bval_df.groupby("assignment", as_index=False).agg({"bvalue": "median"}) + bval_df["assigned_shell"] = bval_df["assignment"].replace( + shell_df["assignment"].tolist(), shell_df["bvalue"].tolist() + ) + bval_df["shell_num"] = bval_df["assigned_shell"].rank(method="dense") + bval_df.drop("assigned_shell", inplace=True) - If bval_assignments is defined, - """ + return bval_df class _RemoveDuplicatesInputSpec(BaseInterfaceInputSpec): dwi_file = File(exists=True, mandatory=True) bval_file = File(exists=True, mandatory=True) bvec_file = File(exists=True, mandatory=True) + b_file = File(exists=True) + btable_file = File(exists=True) local_bvec_file = File(exists=True) distance_cutoff = traits.Float(5.0, usedefault=True) expected_directions = traits.Int() @@ -106,6 +199,8 @@ class _RemoveDuplicatesOutputSpec(TraitedSpec): dwi_file = File(exists=True) bval_file = File(exists=True) bvec_file = File(exists=True) + b_file = File(exists=True) + btable_file = File(exists=True) local_bvec_file = File(exists=True) @@ -116,10 +211,8 @@ class RemoveDuplicates(SimpleInterface): def _run_interface(self, runtime): bvecs = np.loadtxt(self.inputs.bvec_file).T bvals = np.loadtxt(self.inputs.bval_file).squeeze() - orig_bvals = bvals.copy() bvals = np.sqrt(bvals - bvals.min()) bvals = bvals / bvals.max() * 100 - original_image = nb.load(self.inputs.dwi_file) cutoff = self.inputs.distance_cutoff scaled_bvecs = bvals[:, np.newaxis] * bvecs @@ -152,32 +245,25 @@ def is_unique_sample(vec): "Expected %d unique samples but found %d", expected, len(seen_vecs) ) - # Are all the directions unique? - if len(ok_vecs) == len(bvals): - self._results["dwi_file"] = self.inputs.dwi_file - self._results["bval_file"] = self.inputs.bval_file - self._results["bvec_file"] = self.inputs.bvec_file - self._results["local_bvec_file"] = self.inputs.local_bvec_file - - return runtime - # Extract the unique samples - output_bval = fname_presuffix(self.inputs.bval_file, newpath=runtime.cwd, suffix="_unique") - output_bvec = fname_presuffix(self.inputs.bvec_file, newpath=runtime.cwd, suffix="_unique") - output_nii = fname_presuffix(self.inputs.dwi_file, newpath=runtime.cwd, suffix="_unique") - unique_indices = np.array(ok_vecs) - unique_bvals = orig_bvals[unique_indices] - np.savetxt(output_bval, unique_bvals, fmt="%d", newline=" ") - unique_bvecs = bvecs[unique_indices] - np.savetxt(output_bvec, unique_bvecs.T, fmt="%.8f") - unique_data = original_image.get_fdata()[..., unique_indices] - nb.Nifti1Image(unique_data, original_image.affine, original_image.header).to_filename( - output_nii + new_bval, new_bvec, new_b, new_btable, new_nifti = subset_dwi( + original_bval=self.inputs.bval_file, + original_bvec=self.inputs.bvec_file, + original_b=self.inputs.b_file or None, + original_btable=self.inputs.btable_file or None, + original_nifti=self.inputs.dwi_file, + indices=np.array(ok_vecs), + newdir=runtime.cwd, + suffix="_unique", ) - self._results["bval_file"] = output_bval - self._results["bvec_file"] = output_bvec - self._results["dwi_file"] = output_nii - # TODO: support local bvecs + self._results["bval_file"] = new_bval + self._results["bvec_file"] = new_bvec + self._results["dwi_file"] = new_nifti + if new_b: + self._results["b_file"] = new_b + if new_btable: + self._results["btable_file"] = new_btable + return runtime @@ -251,3 +337,54 @@ def concatenate_bvecs(input_files): if not stacked.shape[1] == 3: stacked = stacked.T return stacked + + +def subset_dwi( + original_bval, + original_bvec, + original_b, + original_btable, + original_nifti, + indices, + newdir, + suffix="_bsel", +): + """Create a subset of a dwi based on a set of indices.""" + bvals = np.loadtxt(original_bval) + if np.all(indices == np.arange(len(bvals))): + return original_bval, original_bvec, original_b, original_btable, original_nifti + + # Subset and write the bval + new_bval = fname_presuffix(original_bval, newpath=newdir, suffix=suffix) + subsetted_bval_data = np.loadtxt(original_bval).squeeze()[indices] + np.savetxt(new_bval, subsetted_bval_data, fmt="%d", newline=" ") + + # Subset and write the bvec + new_bvec = fname_presuffix(original_bvec, newpath=newdir, suffix=suffix) + selected_bvecs = np.loadtxt(original_bvec)[:, indices] + np.savetxt(new_bvec, selected_bvecs, fmt="%.8f") + + # Subset and write the dwi nifti + new_nifti = fname_presuffix(original_nifti, newpath=newdir, suffix=suffix) + nim.index_img(original_nifti, indices).to_filename(new_nifti) + + # If there is a b file, subset and write it + if original_b is not None: + new_b = fname_presuffix(original_b, newpath=newdir, suffix=suffix) + _select_lines(original_b, new_b, indices) + + # If there is a dsi studio btable file, subset and write it + if original_btable is not None: + new_btable = fname_presuffix(original_btable, newpath=newdir, suffix=suffix) + _select_lines(original_btable, new_btable) + return new_bval, new_bvec, new_b, new_btable, new_nifti + + +def _select_lines(in_file, out_file, indices): + + with open(in_file, "r") as in_f: + in_lines = in_f.readlines() + new_lines = [in_lines[lineno] for lineno in indices] + + with open(out_file, "w") as out_f: + out_f.writelines(new_lines) diff --git a/qsirecon/workflows/recon/build_workflow.py b/qsirecon/workflows/recon/build_workflow.py index cc891402..844e7ddb 100644 --- a/qsirecon/workflows/recon/build_workflow.py +++ b/qsirecon/workflows/recon/build_workflow.py @@ -27,7 +27,12 @@ from .scalar_mapping import init_scalar_to_bundle_wf, init_scalar_to_template_wf from .steinhardt import init_steinhardt_order_param_wf from .tortoise import init_tortoise_estimator_wf -from .utils import init_conform_dwi_wf, init_discard_repeated_samples_wf, init_test_wf +from .utils import ( + init_conform_dwi_wf, + init_discard_repeated_samples_wf, + init_gradient_select_wf, + init_test_wf, +) def _check_repeats(nodelist): @@ -272,6 +277,8 @@ def workflow_from_spec(inputs_dict, node_spec): else: if node_spec["action"] == "discard_repeated_samples": return init_discard_repeated_samples_wf(**kwargs) + if node_spec["action"] == "select_gradients": + return init_gradient_select_wf(**kwargs) if node_spec["action"] == "conform": return init_conform_dwi_wf(**kwargs) if node_spec["action"] == "mif_to_fib": diff --git a/qsirecon/workflows/recon/fod_autotrack.py b/qsirecon/workflows/recon/fod_autotrack.py index fea2b7ee..757f1ed7 100644 --- a/qsirecon/workflows/recon/fod_autotrack.py +++ b/qsirecon/workflows/recon/fod_autotrack.py @@ -1,11 +1,11 @@ import logging + import nipype.pipeline.engine as pe from nipype.interfaces import utility as niu from ... import config from ...engine import Workflow from ...interfaces.interchange import recon_workflow_input_fields - from .dsi_studio import ( init_dsi_studio_autotrack_wf, init_dsi_studio_connectivity_wf, @@ -97,23 +97,35 @@ def init_singleshell_benchmarking_wf( AggregateAutoTrackResults(expected_bundles=bundle_names), name="aggregate_csd_atk_results" ) - convert_gqi_to_tck = pe.MapNode(DSIStudioTrkToTck(), name="convert_gqi_to_tck", iterfield="trk_file") - convert_ss3t_to_tck = pe.MapNode(DSIStudioTrkToTck(), name="convert_ss3t_to_tck", iterfield="trk_file") - convert_csd_to_tck = pe.MapNode(DSIStudioTrkToTck(), name="convert_csd_to_tck", iterfield="trk_file") + convert_gqi_to_tck = pe.MapNode( + DSIStudioTrkToTck(), name="convert_gqi_to_tck", iterfield="trk_file" + ) + convert_ss3t_to_tck = pe.MapNode( + DSIStudioTrkToTck(), name="convert_ss3t_to_tck", iterfield="trk_file" + ) + convert_csd_to_tck = pe.MapNode( + DSIStudioTrkToTck(), name="convert_csd_to_tck", iterfield="trk_file" + ) # Save the bundle csv ds_gqi_bundle_csv = pe.Node( - ReconDerivativesDataSink(suffix="bundlestats", qsirecon_suffix=f"{qsirecon_suffix}_part-GQI"), + ReconDerivativesDataSink( + suffix="bundlestats", qsirecon_suffix=f"{qsirecon_suffix}_part-GQI" + ), name="ds_gqi_bundle_csv", run_without_submitting=True, ) ds_ss3t_bundle_csv = pe.Node( - ReconDerivativesDataSink(suffix="bundlestats", qsirecon_suffix=f"{qsirecon_suffix}_part-SS3T"), + ReconDerivativesDataSink( + suffix="bundlestats", qsirecon_suffix=f"{qsirecon_suffix}_part-SS3T" + ), name="ds_ss3t_bundle_csv", run_without_submitting=True, ) ds_csd_bundle_csv = pe.Node( - ReconDerivativesDataSink(suffix="bundlestats", qsirecon_suffix=f"{qsirecon_suffix}_part-CSD"), + ReconDerivativesDataSink( + suffix="bundlestats", qsirecon_suffix=f"{qsirecon_suffix}_part-CSD" + ), name="ds_csd_bundle_csv", run_without_submitting=True, ) @@ -153,4 +165,4 @@ def init_singleshell_benchmarking_wf( def _as_connections(attr_list, src_prefix="", dest_prefix=""): - return [(src_prefix + item, dest_prefix + item) for item in attr_list] \ No newline at end of file + return [(src_prefix + item, dest_prefix + item) for item in attr_list] diff --git a/qsirecon/workflows/recon/utils.py b/qsirecon/workflows/recon/utils.py index aa1c3863..eb7efa18 100644 --- a/qsirecon/workflows/recon/utils.py +++ b/qsirecon/workflows/recon/utils.py @@ -14,7 +14,7 @@ from ...interfaces import ConformDwi from ...interfaces.bids import DerivativesDataSink -from ...interfaces.gradients import RemoveDuplicates, GradientSelect +from ...interfaces.gradients import GradientSelect, RemoveDuplicates from ...interfaces.interchange import recon_workflow_input_fields from ...interfaces.mrtrix import MRTrixGradientTable from ...interfaces.recon_scalars import OrganizeScalarData From 94a73fc026df0fd3915398ceea1f79a17d4c50a1 Mon Sep 17 00:00:00 2001 From: Matt Cieslak Date: Thu, 5 Dec 2024 09:15:33 -0500 Subject: [PATCH 03/13] use shell select for ss3t --- .../pipelines/mrtrix_singleshell_ss3t_ACT-fast.yaml | 2 +- .../pipelines/mrtrix_singleshell_ss3t_ACT-hsvs.yaml | 10 +++++++++- .../pipelines/mrtrix_singleshell_ss3t_noACT.yaml | 12 +++++++++++- 3 files changed, 21 insertions(+), 3 deletions(-) diff --git a/qsirecon/data/pipelines/mrtrix_singleshell_ss3t_ACT-fast.yaml b/qsirecon/data/pipelines/mrtrix_singleshell_ss3t_ACT-fast.yaml index c8348905..aba3306d 100644 --- a/qsirecon/data/pipelines/mrtrix_singleshell_ss3t_ACT-fast.yaml +++ b/qsirecon/data/pipelines/mrtrix_singleshell_ss3t_ACT-fast.yaml @@ -3,7 +3,7 @@ anatomical: name: mrtrix_singleshell_ss3t_fast nodes: -- action: +- action: select_gradients input: qsirecon name: select_single_shell parameters: diff --git a/qsirecon/data/pipelines/mrtrix_singleshell_ss3t_ACT-hsvs.yaml b/qsirecon/data/pipelines/mrtrix_singleshell_ss3t_ACT-hsvs.yaml index f882cd96..2eb4cd63 100644 --- a/qsirecon/data/pipelines/mrtrix_singleshell_ss3t_ACT-hsvs.yaml +++ b/qsirecon/data/pipelines/mrtrix_singleshell_ss3t_ACT-hsvs.yaml @@ -2,8 +2,16 @@ anatomical: - mrtrix_5tt_hsvs name: mrtrix_singleshell_ss3_hsvst nodes: -- action: csd +- action: select_gradients input: qsirecon + name: select_single_shell + parameters: + requested_shells: + - 0 + - highest + bval_distance_cutoff: 100 +- action: csd + input: select_single_shell name: ss3t_csd parameters: fod: diff --git a/qsirecon/data/pipelines/mrtrix_singleshell_ss3t_noACT.yaml b/qsirecon/data/pipelines/mrtrix_singleshell_ss3t_noACT.yaml index 543af02d..75b2ce80 100644 --- a/qsirecon/data/pipelines/mrtrix_singleshell_ss3t_noACT.yaml +++ b/qsirecon/data/pipelines/mrtrix_singleshell_ss3t_noACT.yaml @@ -1,8 +1,18 @@ anatomical: [] name: mrtrix_singleshell_ss3t_noACT nodes: -- action: csd + +- action: select_gradients input: qsirecon + name: select_single_shell + parameters: + requested_shells: + - 0 + - highest + bval_distance_cutoff: 100 + +- action: csd + input: select_single_shell name: ss3t_csd parameters: fod: From 28ea2c35ffe69eb25da2ba2d56716249d940eb3b Mon Sep 17 00:00:00 2001 From: Matt Cieslak Date: Thu, 5 Dec 2024 09:17:53 -0500 Subject: [PATCH 04/13] save fod_autotrack for another PR --- qsirecon/workflows/recon/fod_autotrack.py | 168 ---------------------- 1 file changed, 168 deletions(-) delete mode 100644 qsirecon/workflows/recon/fod_autotrack.py diff --git a/qsirecon/workflows/recon/fod_autotrack.py b/qsirecon/workflows/recon/fod_autotrack.py deleted file mode 100644 index 757f1ed7..00000000 --- a/qsirecon/workflows/recon/fod_autotrack.py +++ /dev/null @@ -1,168 +0,0 @@ -import logging - -import nipype.pipeline.engine as pe -from nipype.interfaces import utility as niu - -from ... import config -from ...engine import Workflow -from ...interfaces.interchange import recon_workflow_input_fields -from .dsi_studio import ( - init_dsi_studio_autotrack_wf, - init_dsi_studio_connectivity_wf, - init_dsi_studio_export_wf, - init_dsi_studio_recon_wf, - init_dsi_studio_tractography_wf, -) -from .mrtrix import ( - init_global_tractography_wf, - init_mrtrix_connectivity_wf, - init_mrtrix_csd_recon_wf, - init_mrtrix_tractography_wf, -) - -LOGGER = logging.getLogger("nipype.interface") - - -def init_singleshell_benchmarking_wf( - available_anatomical_data, name="_recon", qsirecon_suffix="SingleShellBenchmark", params={} -): - inputnode = pe.Node( - niu.IdentityInterface(fields=recon_workflow_input_fields + ["odf_rois"]), name="inputnode" - ) - outputnode = pe.Node( - niu.IdentityInterface(fields=["tck_files", "bundle_names", "recon_scalars"]), - name="outputnode", - ) - outputnode.inputs.recon_scalars = [] - workflow = Workflow(name=name) - omp_nthreads = config.nipype.omp_nthreads - - autotrack_params = params.get("autotrack", {}) - bundle_names = _get_dsi_studio_bundles(autotrack_params.get("track_id", "")) - bundle_desc = ( - "AutoTrack attempted to reconstruct the following bundles:\n * " - + "\n * ".join(bundle_names) - + "\n\n" - ) - LOGGER.info(bundle_desc) - - # First do a standard GQI reconstruction - gqi_params = params.get("gqi_recon", {}) - initial_gqi_wf = init_dsi_studio_recon_wf( - available_anatomical_data=available_anatomical_data, - name="initial_gqi", - qsirecon_suffix=f"{qsirecon_suffix}_part-GQI", - params=gqi_params, - ) - - # Do an SS3T recon to feed to autotrack - ss3t_params = params.get("ss3t_recon") - ss3t_wf = init_mrtrix_csd_recon_wf( - available_anatomical_data=available_anatomical_data, - name="ss3t_recon", - qsirecon_suffix=f"{qsirecon_suffix}_part-SS3T", - params=ss3t_params, - ) - ss3t_to_fib = pe.Node(FODtoFIBGZ(), name="ss3t_to_fib") - - # For comparison, also do a regular CSD - csd_params = params.get("csd_recon") - csd_wf = init_mrtrix_csd_recon_wf( - available_anatomical_data=available_anatomical_data, - name="csd_recon", - qsirecon_suffix=f"{qsirecon_suffix}_part-CSD", - params=csd_params, - ) - csd_to_fib = pe.Node(FODtoFIBGZ(), name="csd_to_fib") - - # Run autotrack! - gqi_autotrack = pe.Node( - AutoTrack(num_threads=omp_nthreads, **params), name="gqi_autotrack", n_procs=omp_nthreads - ) - ss3t_autotrack = pe.Node( - AutoTrack(num_threads=omp_nthreads, **params), name="ss3t_autotrack", n_procs=omp_nthreads - ) - csd_autotrack = pe.Node( - AutoTrack(num_threads=omp_nthreads, **params), name="csd_autotrack", n_procs=omp_nthreads - ) - - # Create a single output - aggregate_gqi_atk_results = pe.Node( - AggregateAutoTrackResults(expected_bundles=bundle_names), name="aggregate_gqi_atk_results" - ) - aggregate_ss3t_atk_results = pe.Node( - AggregateAutoTrackResults(expected_bundles=bundle_names), name="aggregate_ss3t_atk_results" - ) - aggregate_csd_atk_results = pe.Node( - AggregateAutoTrackResults(expected_bundles=bundle_names), name="aggregate_csd_atk_results" - ) - - convert_gqi_to_tck = pe.MapNode( - DSIStudioTrkToTck(), name="convert_gqi_to_tck", iterfield="trk_file" - ) - convert_ss3t_to_tck = pe.MapNode( - DSIStudioTrkToTck(), name="convert_ss3t_to_tck", iterfield="trk_file" - ) - convert_csd_to_tck = pe.MapNode( - DSIStudioTrkToTck(), name="convert_csd_to_tck", iterfield="trk_file" - ) - - # Save the bundle csv - ds_gqi_bundle_csv = pe.Node( - ReconDerivativesDataSink( - suffix="bundlestats", qsirecon_suffix=f"{qsirecon_suffix}_part-GQI" - ), - name="ds_gqi_bundle_csv", - run_without_submitting=True, - ) - ds_ss3t_bundle_csv = pe.Node( - ReconDerivativesDataSink( - suffix="bundlestats", qsirecon_suffix=f"{qsirecon_suffix}_part-SS3T" - ), - name="ds_ss3t_bundle_csv", - run_without_submitting=True, - ) - ds_csd_bundle_csv = pe.Node( - ReconDerivativesDataSink( - suffix="bundlestats", qsirecon_suffix=f"{qsirecon_suffix}_part-CSD" - ), - name="ds_csd_bundle_csv", - run_without_submitting=True, - ) - - # Save the mapping file. We're only using the mapping from GQI - ds_mapping = pe.Node( - ReconDerivativesDataSink(suffix="mapping", qsirecon_suffix=qsirecon_suffix), - name="ds_mapping", - run_without_submitting=True, - ) - - workflow.connect([ - # Connect the qsiprep inputs to the recon workflows we're creating here - # (normally this is done in build_workflow()) - (inputnode, initial_gqi_wf, - _as_connections(recon_workflow_input_fields, dest_prefix='inputnode.')), - (inputnode, ss3t_wf, - _as_connections(recon_workflow_input_fields, dest_prefix='inputnode.')), - (inputnode, csd_wf, - _as_connections(recon_workflow_input_fields, dest_prefix='inputnode.')), - - # Convert the sh mifs from csd and ss3t into fib files - (csd_wf, ss3t_to_fib, [("outputnode.mif_file", "mif_file")]), - (initial_gqi_wf, csd_to_fib, [("outputnode.fibgz", "fib_file")]), - (ss3t_wf, ss3t_to_fib, [("outputnode.mif_file", "mif_file")]), - (initial_gqi_wf, ss3t_to_fib, [("outputnode.fibgz", "fib_file")]), - - # Send the fib files to autotrack. Use the map file from gqi in ss3t and csd - (initial_gqi_wf, gqi_autotrack, [("outputnode.fibgz", "fib_file")]), - (ss3t_to_fib, ss3t_autotrack, [("fib_file", "fib_file")]), - (gqi_autotrack, ss3t_autotrack, [("map_file", "map_file")]), - (csd_to_fib, csd_autotrack, [("fib_file", "fib_file")]), - (gqi_autotrack, csd_autotrack, [("map_file", "map_file")]), - ]) # fmt:skip - - return workflow - - -def _as_connections(attr_list, src_prefix="", dest_prefix=""): - return [(src_prefix + item, dest_prefix + item) for item in attr_list] From 1488112083328a5b19068bc9c5fa258a88d49dd4 Mon Sep 17 00:00:00 2001 From: Matt Cieslak Date: Thu, 5 Dec 2024 10:40:51 -0500 Subject: [PATCH 05/13] set up tests --- .circleci/config.yml | 28 ++++++++++++++++++++++++++++ qsirecon/tests/test_interfaces.py | 23 +++++++++++++++++++++++ 2 files changed, 51 insertions(+) create mode 100644 qsirecon/tests/test_interfaces.py diff --git a/.circleci/config.yml b/.circleci/config.yml index dec58598..937bde1a 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -249,6 +249,25 @@ jobs: - store_artifacts: path: /src/qsirecon/.circleci/out/multises_post1_qsiprep_reportsession/ + Recon_Interfaces: + <<: *dockersetup + steps: + - checkout + - restore_cache: + key: multishell_output-01 + - run: *runinstall + - run: + name: Pytest the recon interfaces + no_output_timeout: 1h + command: | + pytest -rP -o log_cli=true -m "interfaces" --cov-config=/src/qsirecon/pyproject.toml --cov-append --cov-report term-missing --cov=qsirecon --data_dir=/src/qsirecon/.circleci/data --output_dir=/src/qsirecon/.circleci/out --working_dir=/src/qsirecon/.circleci/work qsirecon + mkdir /src/coverage + mv /src/qsirecon/.coverage /src/coverage/.coverage.recon_interfaces + - persist_to_workspace: + root: /src/coverage/ + paths: + - .coverage.recon_interfaces + Recon_MRtrix3: <<: *dockersetup steps: @@ -632,6 +651,13 @@ workflows: tags: only: /.*/ + - Recon_Interfaces: + requires: + - download_multishell_output + filters: + tags: + only: /.*/ + - Recon_AutoTrack: requires: - download_multishell_output @@ -693,6 +719,7 @@ workflows: - Recon_3Tissue_Singleshell_ACT - Recon_3Tissue_Singleshell_NoACT - Recon_MRtrix3 + - Recon_Interfaces - Recon_AutoTrack - Recon_Tortoise - Recon_DIPY_MAPMRI @@ -717,6 +744,7 @@ workflows: - Recon_3Tissue_Singleshell_ACT - Recon_3Tissue_Singleshell_NoACT - Recon_MRtrix3 + - Recon_Interfaces - Recon_msmt_Multishell_HSVS - Recon_AutoTrack - Recon_Tortoise diff --git a/qsirecon/tests/test_interfaces.py b/qsirecon/tests/test_interfaces.py new file mode 100644 index 00000000..99faf429 --- /dev/null +++ b/qsirecon/tests/test_interfaces.py @@ -0,0 +1,23 @@ +"""Command-line interface tests.""" + +import os + +import pytest + +from qsirecon.tests.utils import download_test_data + + +@pytest.mark.integration +@pytest.mark.interfaces +def test_shell_selection(data_dir, working_dir): + """Run reconstruction workflow tests. + + + """ + TEST_NAME = "mrtrix_singleshell_ss3t_act" + + dataset_dir = download_test_data("multishell_output", data_dir) + dataset_dir = os.path.join(dataset_dir, "qsiprep") + + + From 7df58c3e91d93d2fde28a7fafcc8fad6c2245e09 Mon Sep 17 00:00:00 2001 From: Matt Cieslak Date: Thu, 5 Dec 2024 14:25:57 -0500 Subject: [PATCH 06/13] Update qsirecon/interfaces/gradients.py Co-authored-by: Taylor Salo --- qsirecon/interfaces/gradients.py | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/qsirecon/interfaces/gradients.py b/qsirecon/interfaces/gradients.py index 2e9a6d3c..1df7ba33 100644 --- a/qsirecon/interfaces/gradients.py +++ b/qsirecon/interfaces/gradients.py @@ -166,14 +166,21 @@ def _find_shells(bvals, max_distance): import pandas as pd from sklearn.cluster import AgglomerativeClustering + from sklearn.metrics import silhouette_score + X = bvals.reshape(-1, 1) agg_cluster = AgglomerativeClustering( n_clusters=None, distance_threshold=2 * max_distance, linkage="complete", - ).fit(bvals.reshape(-1, 1)) - - bval_df = pd.DataFrame({"bvalue": bvals, "assignment": agg_cluster.labels_}) + ).fit(X) + shells = agg_cluster.labels_ + + score = silhouette_score(X, shells) + if score < 0.8: + print('Silhouette score is low. Is this is a DSI scheme?') + + bval_df = pd.DataFrame({"bvalue": bvals, "assignment": shells}) shell_df = bval_df.groupby("assignment", as_index=False).agg({"bvalue": "median"}) bval_df["assigned_shell"] = bval_df["assignment"].replace( shell_df["assignment"].tolist(), shell_df["bvalue"].tolist() From cc82bc62eff48a1882639a6f07215d11fbc3fa68 Mon Sep 17 00:00:00 2001 From: Matt Cieslak Date: Thu, 5 Dec 2024 14:55:06 -0500 Subject: [PATCH 07/13] add tests --- qsirecon/interfaces/gradients.py | 14 ++++----- qsirecon/tests/test_interfaces.py | 48 +++++++++++++++++++++++-------- qsirecon/workflows/recon/utils.py | 8 ++++++ 3 files changed, 51 insertions(+), 19 deletions(-) diff --git a/qsirecon/interfaces/gradients.py b/qsirecon/interfaces/gradients.py index 1df7ba33..671d83e0 100644 --- a/qsirecon/interfaces/gradients.py +++ b/qsirecon/interfaces/gradients.py @@ -22,8 +22,8 @@ class _GradientSelectInputSpec(BaseInterfaceInputSpec): dwi_file = File(exists=True, mandatory=True) bval_file = File(exists=True, mandatory=True) bvec_file = File(exists=True, mandatory=True) - b_file = File(exists=True, mandatory=True) - btable_file = File(exists=True, mandatory=True) + b_file = File(exists=True) + btable_file = File(exists=True) bval_distance_cutoff = traits.Float(5.0, usedefault=True) expected_n_input_shells = traits.Int() requested_shells = InputMultiObject( @@ -106,10 +106,10 @@ def _select_gradients( ) # Ensure that at lease 1 b>0 and one b=0 shell are selected - select_shell_bs = _parse_shell_selection(requested_gradients, bval_df) + select_shell_bs = _parse_shell_selection(requested_gradients, bval_df, max_distance) # Make sure the shells are unique (no overlap when accounting for allowed distances) - shell_distances = pairwise_distances(select_shell_bs) + shell_distances = pairwise_distances(select_shell_bs.reshape(-1, 1)) if np.any(shell_distances[np.triu_indices_from(shell_distances, k=1)] < max_distance): raise Exception( "Shells and bval_distance_cutoff have overlap. Choose a lower " @@ -159,7 +159,7 @@ def _parse_shell_selection(requested_bvals, bval_df, max_distance): LOGGER.warning("Adding a b=0 shell to the selection.") numeric_bvals.append(0) - return sorted(numeric_bvals) + return np.array(sorted(numeric_bvals)) def _find_shells(bvals, max_distance): @@ -186,7 +186,7 @@ def _find_shells(bvals, max_distance): shell_df["assignment"].tolist(), shell_df["bvalue"].tolist() ) bval_df["shell_num"] = bval_df["assigned_shell"].rank(method="dense") - bval_df.drop("assigned_shell", inplace=True) + bval_df.drop(columns=["assignment"], inplace=True) return bval_df @@ -383,7 +383,7 @@ def subset_dwi( # If there is a dsi studio btable file, subset and write it if original_btable is not None: new_btable = fname_presuffix(original_btable, newpath=newdir, suffix=suffix) - _select_lines(original_btable, new_btable) + _select_lines(original_btable, new_btable, indices) return new_bval, new_bvec, new_b, new_btable, new_nifti diff --git a/qsirecon/tests/test_interfaces.py b/qsirecon/tests/test_interfaces.py index 99faf429..2d233f85 100644 --- a/qsirecon/tests/test_interfaces.py +++ b/qsirecon/tests/test_interfaces.py @@ -1,23 +1,47 @@ """Command-line interface tests.""" -import os +from pathlib import Path +import nibabel as nb +import numpy as np import pytest +from qsirecon.interfaces.gradients import GradientSelect from qsirecon.tests.utils import download_test_data @pytest.mark.integration @pytest.mark.interfaces def test_shell_selection(data_dir, working_dir): - """Run reconstruction workflow tests. - - - """ - TEST_NAME = "mrtrix_singleshell_ss3t_act" - - dataset_dir = download_test_data("multishell_output", data_dir) - dataset_dir = os.path.join(dataset_dir, "qsiprep") - - - + """Run reconstruction workflow tests.""" + dwi_prefix = "sub-ABCD_acq-10per000_space-T1w_desc-preproc_dwi" + data_dir = "/home/matt/projects/qsirecon/.circleci/data" + dataset_dir = Path(download_test_data("multishell_output", data_dir)) + dataset_dir = Path(dataset_dir) / "multishell_output" / "qsiprep" + data_stem = str(dataset_dir / "sub-ABCD" / "dwi" / dwi_prefix) + + # Test the actual data (including a nifti) + grad_select = GradientSelect( + dwi_file=data_stem + ".nii.gz", + bval_file=data_stem + ".bval", + bvec_file=data_stem + ".bvec", + b_file=data_stem + ".b", + btable_file=data_stem + ".b", + bval_distance_cutoff=100, + expected_n_input_shells=5, + requested_shells=[0, "lowest", "highest"], + ) + grad_select.run() + + correct_n = 73 + sel_nii = nb.load(dwi_prefix + "_selected.nii.gz") + assert sel_nii.shape[3] == correct_n + sel_bval = np.loadtxt(dwi_prefix + "_selected.bval") + assert sel_bval.shape == (correct_n,) + sel_bvec = np.loadtxt(dwi_prefix + "_selected.bvec") + assert sel_bvec.shape == (3, correct_n) + sel_b = np.loadtxt(dwi_prefix + "_selected.b") + assert sel_b.shape[0] == correct_n + # There is no btable from this dataset because it was created + # before those were written in the outputs. + assert not Path(dwi_prefix + "_selected.txt").exists() diff --git a/qsirecon/workflows/recon/utils.py b/qsirecon/workflows/recon/utils.py index eb7efa18..1196f762 100644 --- a/qsirecon/workflows/recon/utils.py +++ b/qsirecon/workflows/recon/utils.py @@ -72,10 +72,14 @@ def init_discard_repeated_samples_wf( (inputnode, discard_repeats, [ ('dwi_file', 'dwi_file'), ('bval_file', 'bval_file'), + ('b_file', 'b_file'), + ('btable_file', 'btable_file'), ('bvec_file', 'bvec_file')]), (discard_repeats, outputnode, [ ('dwi_file', 'dwi_file'), ('bval_file', 'bval_file'), + ('b_file', 'b_file'), + ('btable_file', 'btable_file'), ('bvec_file', 'bvec_file')]) ]) # fmt:skip @@ -103,10 +107,14 @@ def init_gradient_select_wf( (inputnode, gradient_select, [ ('dwi_file', 'dwi_file'), ('bval_file', 'bval_file'), + ('b_file', 'b_file'), + ('btable_file', 'btable_file'), ('bvec_file', 'bvec_file')]), (gradient_select, outputnode, [ ('dwi_file', 'dwi_file'), ('bval_file', 'bval_file'), + ('b_file', 'b_file'), + ('btable_file', 'btable_file'), ('bvec_file', 'bvec_file')]) ]) # fmt:skip From 5a8c85cfbc90728012f1d791c24f20d9ea6a877d Mon Sep 17 00:00:00 2001 From: Matt Cieslak Date: Thu, 5 Dec 2024 14:56:04 -0500 Subject: [PATCH 08/13] lint --- qsirecon/interfaces/gradients.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/qsirecon/interfaces/gradients.py b/qsirecon/interfaces/gradients.py index 671d83e0..74f7591f 100644 --- a/qsirecon/interfaces/gradients.py +++ b/qsirecon/interfaces/gradients.py @@ -175,10 +175,10 @@ def _find_shells(bvals, max_distance): linkage="complete", ).fit(X) shells = agg_cluster.labels_ - + score = silhouette_score(X, shells) if score < 0.8: - print('Silhouette score is low. Is this is a DSI scheme?') + print("Silhouette score is low. Is this is a DSI scheme?") bval_df = pd.DataFrame({"bvalue": bvals, "assignment": shells}) shell_df = bval_df.groupby("assignment", as_index=False).agg({"bvalue": "median"}) From 9f57f92d45dcc98a71cc6c90dec69edd348c66e0 Mon Sep 17 00:00:00 2001 From: Matt Cieslak Date: Thu, 5 Dec 2024 17:25:04 -0500 Subject: [PATCH 09/13] check a whole bunch of schemes --- .circleci/config.yml | 4 ++-- qsirecon/data/schemes/DSIQ5.bval | 1 + qsirecon/data/schemes/DSIQ5.bvec | 3 +++ qsirecon/interfaces/gradients.py | 5 +++++ qsirecon/tests/data/HASC55-1.bval | 1 + qsirecon/tests/data/HASC55-2.bval | 1 + qsirecon/tests/data/HASC92.bval | 1 + qsirecon/tests/data/Q7.bval | 1 + qsirecon/tests/data/RAND57.bval | 1 + qsirecon/tests/test_cli.py | 2 +- qsirecon/tests/test_interfaces.py | 35 +++++++++++++++++++++++++------ 11 files changed, 46 insertions(+), 9 deletions(-) create mode 100644 qsirecon/data/schemes/DSIQ5.bval create mode 100644 qsirecon/data/schemes/DSIQ5.bvec create mode 100644 qsirecon/tests/data/HASC55-1.bval create mode 100644 qsirecon/tests/data/HASC55-2.bval create mode 100644 qsirecon/tests/data/HASC92.bval create mode 100644 qsirecon/tests/data/Q7.bval create mode 100644 qsirecon/tests/data/RAND57.bval diff --git a/.circleci/config.yml b/.circleci/config.yml index 937bde1a..dc9b10b4 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -468,7 +468,7 @@ jobs: steps: - checkout - restore_cache: - key: singleshell_output-01 + key: multishell_output-01 - run: *runinstall - run: name: Test the DIPY recon workflows @@ -688,7 +688,7 @@ workflows: - Recon_AMICO: requires: - - download_singleshell_output + - download_multishell_output filters: tags: only: /.*/ diff --git a/qsirecon/data/schemes/DSIQ5.bval b/qsirecon/data/schemes/DSIQ5.bval new file mode 100644 index 00000000..350886c7 --- /dev/null +++ b/qsirecon/data/schemes/DSIQ5.bval @@ -0,0 +1 @@ +0 0 195 200 195 395 390 390 395 395 395 590 590 590 0 590 785 790 800 995 995 980 985 980 985 980 995 995 0 985 985 980 1180 1185 1180 1180 1180 1180 1185 1185 1185 1180 0 1180 1180 1575 1585 1575 1585 1585 1585 1780 1780 1780 1780 1770 0 1780 1770 1770 1775 1780 1780 1795 1780 1770 1780 1970 1970 1970 0 1980 1990 1975 1970 1990 1990 1975 1980 1990 2170 2170 2185 2170 0 2170 2185 2170 2170 2170 2185 2170 2185 2370 2370 2370 2370 2575 0 2575 2580 2565 2565 2565 2580 2575 2575 2580 2580 2565 2770 2770 0 2775 2775 2765 2760 2760 2760 2770 2775 2770 2765 2765 2775 2775 0 2760 2775 2775 2765 2770 2770 2775 2770 2770 3165 3170 3190 3365 0 3365 3365 3365 3360 3385 3365 3385 3360 3360 3360 3360 3360 3360 0 3360 3365 3385 3360 3365 3365 3360 3360 3385 3360 3570 3560 3560 0 3560 3570 3570 3555 3560 3560 3580 3555 3560 3580 3560 3580 3570 0 3560 3580 3755 3765 3765 3765 3765 3765 3755 3765 3765 3755 3765 0 3755 3955 3975 3965 3965 3975 3950 3955 3965 3950 3975 3965 3975 0 4155 4155 4150 4170 4155 4160 4155 4150 4150 4170 4150 4170 4170 0 4170 4150 4170 4150 4160 4160 4150 4160 4170 4150 4170 4355 4350 0 4355 4350 4355 4350 4355 4355 4355 4355 4350 4355 4760 4760 4750 0 4760 4760 4750 4750 4750 4750 4750 4750 4750 4960 4955 4965 0 4960 4965 4985 4940 4965 4960 4940 4940 4960 4950 4940 4965 0 diff --git a/qsirecon/data/schemes/DSIQ5.bvec b/qsirecon/data/schemes/DSIQ5.bvec new file mode 100644 index 00000000..c278b6ec --- /dev/null +++ b/qsirecon/data/schemes/DSIQ5.bvec @@ -0,0 +1,3 @@ +0 0 0.00512677 0.999898 0.0051803 -0.712451 0.00365419 -0.00365419 0.712451 0.712451 -0.712451 -0.583226 -0.583226 0.583226 0 0.583226 0.00126569 0.00127221 0.999961 0.896745 0.897016 -0.00136882 0.450107 0.00136834 0.449748 -0.00136834 -0.896745 -0.897016 0 -0.449748 -0.450107 0.00136882 -0.41136 0.820265 0.41136 -0.411358 0.411358 -0.41136 -0.820265 0.820265 -0.820265 0.411358 0 0.41136 -0.411358 0.00089988 -0.709773 -0.00089988 -0.709553 0.709553 0.709773 0.669846 0.669866 0.00056305 0.669846 0.33559 0 -0.669866 -0.33559 -0.33559 0.00056114 -0.669846 -0.669866 0.999982 0.669866 0.33559 -0.669846 -0.00064193 0.00064193 -0.00064131 0 -0.318 -0.949989 0.318288 0.00064131 0.950113 0.949989 -0.318288 0.318 -0.950113 -0.303708 -0.303555 0.906895 -0.303555 0 -0.303708 0.906895 0.303555 0.303708 0.303555 -0.906895 0.303708 -0.906895 -0.580189 0.580189 -0.580189 0.580189 -0.556709 0 -0.556883 0.834008 -0.0005413 -0.00054141 0.0005413 -0.833828 0.556883 0.556709 0.833828 -0.834008 0.00054141 0.537016 0.537016 0 0.804246 0.804275 0.268996 -0.26901 0.26901 0.26901 0.536932 0.804275 0.536932 0.268996 -0.268996 -0.804275 -0.804246 0 -0.26901 -0.804275 -0.804246 -0.268996 -0.536932 -0.537016 0.804246 -0.537016 -0.536932 0.00031515 0.00031595 0.99999 -0.243569 0 -0.730103 -0.243705 0.243705 -0.00036051 -0.970951 -0.730103 -0.971013 -0.00036075 0.487374 -0.487278 0.487374 0.00036075 0.487278 0 -0.487374 0.730103 0.970951 -0.487278 0.243569 0.730103 0.00036051 0.487278 0.971013 -0.487374 -0.709081 -0.2369 0.23699 0 0.2369 0.708883 0.709081 -0.0003982 -0.23699 0.23699 -0.944344 0.0003982 0.2369 0.944344 -0.23699 0.944344 -0.708883 0 -0.2369 -0.944344 0.230883 -0.690449 0.690449 0.690559 -0.690559 0.690449 0.230883 -0.690559 -0.690449 -0.230883 0.690559 0 -0.230883 0.00033944 0.895814 0.448611 -0.448791 -0.895814 -0.00033959 -0.00033944 -0.448611 0.00033959 -0.895704 0.448791 0.895704 0 -0.438229 -0.438229 -0.219309 0.874687 0.438229 -0.438076 0.438229 0.21938 -0.21938 0.874715 0.219309 0.874687 -0.874687 0 -0.874715 -0.21938 0.874715 0.21938 0.438076 -0.438076 0.219309 0.438076 -0.874687 -0.219309 -0.874715 0.641837 0.428312 0 -0.641842 0.428312 0.641842 -0.428312 -0.641837 -0.641842 -0.641837 0.641837 -0.428312 0.641842 0.818493 0.818493 0.409953 0 -0.818493 -0.818493 0.409852 -0.409953 -0.409953 -0.409852 -0.409852 0.409953 0.409852 -0.60152 0.0002019 0.801575 0 0.60152 0.801392 0.999993 0.00028335 -0.801575 -0.60163 -0.00028335 -0.00028339 0.60163 0.00020149 0.00028339 -0.801392 0 +0 0 -0.010334 -0.0100979 -0.999933 0.00719586 -0.707102 0.707102 -0.00719586 -0.701685 -0.701685 -0.57439 0.57439 0.57439 0 -0.57439 -0.00636733 -0.999986 -0.00628066 -0.442514 -0.00541963 0.445206 -0.00544864 -0.895882 -0.893144 0.895882 -0.442514 0.00541963 0 -0.893144 0.00544864 -0.445206 0.405824 -0.404799 0.405824 -0.816489 0.816489 -0.405824 -0.404799 0.404799 0.404799 -0.816489 0 -0.405824 0.816489 -0.707331 0.0044666 0.707331 -0.704641 -0.704641 -0.0044666 -0.665159 0.330623 -0.999992 0.665159 -0.666288 0 0.330623 0.666288 -0.666288 -0.00450913 0.665159 -0.330623 -0.00446125 -0.330623 0.666288 -0.665159 0.313963 -0.313963 0.949658 0 -0.948084 -0.312263 -0.00416406 -0.949658 -0.00413844 -0.312263 0.00416406 -0.948084 0.00413844 0.299122 -0.904832 -0.298153 0.904832 0 -0.299122 0.298153 0.904832 0.299122 -0.904832 0.298153 -0.299122 -0.298153 -0.576047 0.576047 0.576047 -0.576047 -0.830701 0 0.00366618 -0.00365633 0.832732 0.554251 -0.832732 -0.552016 -0.00366618 -0.830701 -0.552016 0.00365633 -0.554251 -0.264935 0.264935 0 -0.532439 0.264526 0.802052 0.533761 0.533761 -0.533761 0.801151 -0.264526 -0.801151 -0.802052 0.802052 -0.264526 -0.532439 0 -0.533761 0.264526 0.532439 -0.802052 -0.801151 -0.264935 0.532439 0.264935 0.801151 -0.00316136 -0.999996 -0.00345014 -0.96988 0 -0.483373 0.00302909 -0.00302909 0.97081 -0.239261 0.483373 0.00329932 0.240241 0.727129 -0.483957 -0.727129 -0.240241 -0.483957 0 -0.727129 -0.483373 -0.239261 0.483957 -0.96988 0.483373 -0.97081 0.483957 -0.00329932 0.727129 0.00316998 0.943166 -0.233412 0 -0.943166 -0.70532 -0.00316998 0.707306 -0.233412 0.233412 -0.232741 -0.707306 0.943166 0.232741 0.233412 -0.232741 -0.70532 0 -0.943166 0.232741 0.688185 -0.686976 -0.686976 -0.227037 -0.227037 0.686976 -0.688185 0.227037 0.686976 0.688185 0.227037 0 -0.688185 -0.895083 -0.00303881 -0.893723 0.00282259 0.00303881 0.446171 0.895083 -0.893723 -0.446171 -0.444645 -0.00282259 -0.444645 0 0.216017 -0.216017 0.873293 0.434228 0.216017 0.872688 -0.216017 -0.435288 0.435288 0.215668 -0.873293 -0.434228 0.434228 0 0.215668 -0.435288 -0.215668 0.435288 0.872688 -0.872688 0.873293 -0.872688 -0.434228 -0.873293 -0.215668 -0.638607 0.63916 0 -0.424958 -0.63916 0.424958 -0.63916 0.638607 0.424958 -0.638607 0.638607 0.63916 -0.424958 -0.406373 0.406373 0.406943 0 -0.406373 0.406373 0.816437 -0.406943 0.406943 0.816437 -0.816437 -0.406943 -0.816437 -0.798855 -0.999997 -0.00274017 0 -0.798855 -0.598135 -0.00281054 -0.800461 0.00274017 0.00258182 0.800461 0.599696 -0.00258182 -0.00262611 -0.599696 -0.598135 0 +0 0 -0.999933 -0.0100979 -0.0103606 -0.701685 -0.707102 -0.707102 -0.701685 -0.00719583 -0.00719583 -0.57439 -0.57439 -0.57439 0 -0.57439 -0.999979 -0.00508885 -0.00628066 -0.00541502 -0.441964 -0.895427 -0.892958 -0.44429 -0.0045382 -0.44429 -0.00541502 -0.441964 0 -0.0045382 -0.892958 -0.895427 -0.816143 -0.404107 -0.816143 -0.405131 -0.405131 -0.816143 -0.404107 -0.404107 -0.404107 -0.405131 0 -0.816143 -0.405131 -0.706882 -0.704416 -0.706882 -0.0040203 -0.0040203 -0.704416 -0.329954 -0.664806 -0.00394135 -0.329954 -0.665912 0 -0.664806 -0.665912 -0.665912 -0.99999 -0.329954 -0.664806 -0.00390418 -0.664806 -0.665912 -0.329954 -0.949435 -0.949435 -0.313288 0 -0.00368058 -0.00365925 -0.947985 -0.313288 -0.311877 -0.00365925 -0.947985 -0.00368058 -0.311877 -0.904593 -0.298554 -0.297737 -0.298554 0 -0.904593 -0.297737 -0.298554 -0.904593 -0.298554 -0.297737 -0.904593 -0.297737 -0.575804 -0.575804 -0.575804 -0.575804 -0.00334024 0 -0.830583 -0.55174 -0.553675 -0.83235 -0.553675 -0.00311754 -0.830583 -0.00334024 -0.00311754 -0.55174 -0.83235 -0.800889 -0.800889 0 -0.264003 -0.532135 -0.533249 -0.801706 -0.801706 -0.801706 -0.264314 -0.532135 -0.264314 -0.533249 -0.533249 -0.532135 -0.264003 0 -0.801706 -0.532135 -0.264003 -0.533249 -0.264314 -0.800889 -0.264003 -0.800889 -0.264314 -0.999995 -0.00284357 -0.0028232 -0.00273906 0 -0.483012 -0.969845 -0.969845 -0.23985 -0.00272476 -0.483012 -0.239006 -0.970713 -0.483477 -0.726874 -0.483477 -0.970713 -0.726874 0 -0.483477 -0.483012 -0.00272476 -0.726874 -0.00273906 -0.483012 -0.23985 -0.726874 -0.239006 -0.483477 -0.70512 -0.233059 -0.943056 0 -0.233059 -0.00277213 -0.70512 -0.706908 -0.943056 -0.943056 -0.232477 -0.706908 -0.233059 -0.232477 -0.943056 -0.232477 -0.00277213 0 -0.233059 -0.232477 -0.687819 -0.226595 -0.226595 -0.686719 -0.686719 -0.226595 -0.687819 -0.686719 -0.226595 -0.687819 -0.686719 0 -0.687819 -0.4459 -0.444418 -0.00259589 -0.893632 -0.444418 -0.894947 -0.4459 -0.00259589 -0.894947 -0.00247554 -0.893632 -0.00247554 0 -0.872521 -0.872521 -0.435043 -0.215334 -0.872521 -0.215652 -0.872521 -0.873154 -0.873154 -0.434006 -0.435043 -0.215334 -0.215334 0 -0.434006 -0.873154 -0.434006 -0.873154 -0.215652 -0.215652 -0.435043 -0.215652 -0.215334 -0.435043 -0.434006 -0.424531 -0.638767 0 -0.638318 -0.638767 -0.638318 -0.638767 -0.424531 -0.638318 -0.424531 -0.424531 -0.638767 -0.638318 -0.406116 -0.406116 -0.816294 0 -0.406116 -0.406116 -0.406757 -0.816294 -0.816294 -0.406757 -0.406757 -0.816294 -0.406757 -0.0022984 -0.00242275 -0.597888 0 -0.0022984 -0.00229544 -0.00240919 -0.599385 -0.597888 -0.798771 -0.599385 -0.800228 -0.798771 -0.999997 -0.800228 -0.00229544 0 diff --git a/qsirecon/interfaces/gradients.py b/qsirecon/interfaces/gradients.py index 74f7591f..88ee54cd 100644 --- a/qsirecon/interfaces/gradients.py +++ b/qsirecon/interfaces/gradients.py @@ -180,6 +180,11 @@ def _find_shells(bvals, max_distance): if score < 0.8: print("Silhouette score is low. Is this is a DSI scheme?") + # Do the same check as mrtrix + max_shells = np.sqrt(np.sum(bvals > max_distance)) + if agg_cluster.n_clusters_ > max_shells: + raise Exception("Too many possible shells detected.") + bval_df = pd.DataFrame({"bvalue": bvals, "assignment": shells}) shell_df = bval_df.groupby("assignment", as_index=False).agg({"bvalue": "median"}) bval_df["assigned_shell"] = bval_df["assignment"].replace( diff --git a/qsirecon/tests/data/HASC55-1.bval b/qsirecon/tests/data/HASC55-1.bval new file mode 100644 index 00000000..b9abf862 --- /dev/null +++ b/qsirecon/tests/data/HASC55-1.bval @@ -0,0 +1 @@ +5 5 3395 3400 2595 4395 3795 2795 1995 4190 3600 3395 2795 1595 5 3790 4390 800 3400 3990 1195 3590 2195 4190 4000 2790 5000 5 1795 1795 4195 3395 1195 2795 595 3590 3395 1990 2795 4195 5 3390 3600 4395 4985 4195 3390 3990 3400 2590 3590 995 2790 5000 2395 2000 1795 2190 1195 1195 2595 3790 5 diff --git a/qsirecon/tests/data/HASC55-2.bval b/qsirecon/tests/data/HASC55-2.bval new file mode 100644 index 00000000..9d3ce231 --- /dev/null +++ b/qsirecon/tests/data/HASC55-2.bval @@ -0,0 +1 @@ +5 5 2000 3790 3795 4390 2195 2195 4195 3590 995 2795 3790 3395 5 3595 3995 2590 3595 4195 3590 4190 3390 3400 4195 2395 5000 5 1795 1795 3790 1600 3995 4190 1195 2790 3200 200 1990 3995 5 3395 4390 4985 395 3990 3395 3990 3395 200 3390 4390 5000 995 1795 1000 1195 2600 2195 4390 5 diff --git a/qsirecon/tests/data/HASC92.bval b/qsirecon/tests/data/HASC92.bval new file mode 100644 index 00000000..2272b355 --- /dev/null +++ b/qsirecon/tests/data/HASC92.bval @@ -0,0 +1 @@ +5 5 4195 3590 4790 4990 3790 2590 1800 4195 995 4990 2595 2795 5 5000 3395 1795 3595 2790 2190 4390 5000 1195 1195 2795 4990 5 1795 995 2195 4795 4985 5000 2795 1790 1595 4795 1995 800 5 2790 2190 3790 5000 1595 4990 2795 1195 2195 1800 2000 3595 5 4395 2000 2600 4990 2000 2195 4390 800 2795 1595 2790 4790 5 3395 1195 1795 2195 5000 4990 2395 4795 4795 1000 1195 995 5 2190 1990 4990 2795 3395 995 5000 2595 3400 1990 795 995 1795 1990 1195 1795 1795 2600 2395 4990 1000 1600 1995 5 diff --git a/qsirecon/tests/data/Q7.bval b/qsirecon/tests/data/Q7.bval new file mode 100644 index 00000000..6b192343 --- /dev/null +++ b/qsirecon/tests/data/Q7.bval @@ -0,0 +1 @@ +0 0 105 100 100 205 205 205 200 205 200 300 300 300 300 400 410 400 505 505 510 510 505 500 510 500 510 500 500 505 605 600 605 605 605 600 600 0 605 605 600 605 605 800 810 810 810 810 800 910 905 910 910 905 905 910 910 920 910 905 910 910 905 905 1015 1005 1005 1010 1005 1010 1010 1010 1020 1020 0 1015 1005 1115 1110 1110 1105 1105 1110 1105 1115 1115 1105 1110 1115 1210 1210 1210 1210 1305 1320 1305 1320 1320 1320 1315 1305 1315 1315 1315 1305 1415 1410 1415 1410 1410 1410 0 1410 1415 1415 1415 1415 1410 1415 1415 1410 1415 1410 1410 1415 1410 1410 1410 1415 1415 1615 1615 1630 1715 1715 1715 1715 1715 1715 1715 1730 1715 1715 1715 1730 1730 1715 1720 0 1720 1710 1715 1720 1715 1730 1720 1715 1710 1825 1825 1815 1825 1825 1830 1830 1810 1815 1810 1815 1815 1830 1815 1815 1815 1815 1830 1920 1915 1920 1920 1915 1920 1920 1920 1920 0 1915 1920 1915 2020 2030 2015 2015 2025 2030 2030 2015 2030 2015 2020 2025 2115 2115 2130 2120 2120 2130 2120 2120 2115 2120 2115 2120 2115 2115 2130 2120 2120 2115 2115 2130 2130 0 2130 2130 2130 2220 2225 2225 2220 2225 2225 2220 2225 2220 2225 2225 2225 2420 2420 2430 2420 2420 2420 2420 2430 2420 2430 2430 2420 2520 2525 2520 2520 2520 2535 2530 2530 2535 0 2535 2530 2545 2535 2525 2530 2620 2630 2635 2620 2620 2620 2630 2635 2625 2630 2635 2620 2645 2625 2620 2630 2630 2635 2630 2630 2645 2630 2630 2635 2620 2625 2635 2645 2620 2645 0 2635 2625 2630 2635 2630 2630 2745 2725 2745 2725 2725 2725 2725 2725 2745 2725 2725 2725 2725 2725 2745 2725 2945 2935 2930 2925 2925 2930 2935 2935 2935 2925 2930 2945 2925 2935 0 2935 2925 2935 2935 2935 2945 2925 2925 2930 2925 2930 2930 2935 2925 2925 2925 2930 2935 2925 2935 2930 2945 3030 3045 3030 3030 3045 3025 3045 3035 3045 3045 3045 3030 3025 3030 0 3025 3030 3035 3045 3030 3025 3035 3045 3035 3030 3225 3225 3245 3245 3245 3245 3345 3330 3330 3340 3330 3335 3335 3330 3345 3345 3330 3330 3335 3340 3330 3345 3340 3340 3335 0 3340 3340 3340 3340 3330 3430 3450 3435 3445 3440 3440 3430 3435 3435 3440 3440 3445 3435 3440 3450 3435 3435 3435 3450 3440 3435 3450 3430 3430 3530 3540 3540 3540 3550 3530 0 3535 3550 3550 3530 3550 3550 3535 3550 3530 3540 3550 3540 3535 3540 3540 3550 3540 3535 3635 3635 3640 3640 3645 3645 3635 3665 3640 3645 3640 3640 3640 3645 3635 3765 3765 0 3765 3745 3745 3765 3740 3740 3740 3740 3740 3740 3835 3840 3840 3850 3865 3840 3845 3840 3840 3840 3835 3840 3840 3840 3840 3850 3865 3840 3845 3840 3840 3850 3850 3845 3865 0 3865 3835 3845 3840 3850 3850 3850 3840 3840 3835 3850 4040 4040 4050 4050 4050 4065 4040 4065 4065 4040 4065 4050 4145 4140 4150 4155 4150 4150 4145 4145 4145 4145 4165 4150 0 4160 4145 4155 4155 4150 4150 4140 4140 4145 4160 4145 4140 4165 4140 4165 4160 4150 4155 4140 4145 4145 4165 4145 4140 4145 4165 4145 4165 4165 4160 4145 4140 4150 4165 4145 0 4145 4255 4260 4255 4260 4255 4255 4240 4250 4240 4240 4255 4250 4250 4255 4240 4250 4260 4240 4255 4240 4240 4240 4260 4255 4345 4345 4345 4345 4355 4355 4345 4345 4345 4345 0 4355 4355 4465 4445 4445 4450 4445 4465 4465 4465 4450 4450 4450 4445 4560 4555 4560 4560 4555 4560 4555 4570 4570 4555 4545 4545 4545 4545 4570 4545 4545 4555 4545 4545 4555 0 4570 4545 4545 4560 4560 4560 4560 4555 4545 4555 4560 4555 4555 4560 4545 4655 4670 4655 4655 4670 4655 4670 4650 4655 4655 4655 4670 4655 4670 4670 4645 4645 4650 4645 4650 0 4645 4670 4650 4670 4855 4855 4855 4855 4990 4960 4970 4970 4970 4970 4970 4955 4970 4960 4955 4955 4960 4955 4950 4950 4970 4955 4955 4950 4950 4960 4955 4955 4960 4960 4970 0 diff --git a/qsirecon/tests/data/RAND57.bval b/qsirecon/tests/data/RAND57.bval new file mode 100644 index 00000000..25e7f3ee --- /dev/null +++ b/qsirecon/tests/data/RAND57.bval @@ -0,0 +1 @@ +5 5 4000 4795 4790 4190 3590 4795 4195 4195 3595 3390 2790 5000 5 3995 595 400 1000 3790 395 3795 3795 3395 3200 4795 4985 5 4000 3595 4390 195 3390 200 3390 2790 5000 595 4195 3395 5 2795 3595 4195 2795 2795 3990 3595 1995 4000 4795 4195 4190 5 4790 4190 2600 4195 395 200 4195 3390 4395 4195 395 3190 4190 395 2795 595 2590 400 2795 3590 5 diff --git a/qsirecon/tests/test_cli.py b/qsirecon/tests/test_cli.py index 0733259c..a4cb4612 100644 --- a/qsirecon/tests/test_cli.py +++ b/qsirecon/tests/test_cli.py @@ -312,7 +312,7 @@ def test_amico_noddi(data_dir, output_dir, working_dir): """ TEST_NAME = "amico_noddi" - dataset_dir = download_test_data("singleshell_output", data_dir) + dataset_dir = download_test_data("multishell_output", data_dir) # XXX: Having to modify dataset_dirs is suboptimal. dataset_dir = os.path.join(dataset_dir, "qsiprep") out_dir = os.path.join(output_dir, TEST_NAME) diff --git a/qsirecon/tests/test_interfaces.py b/qsirecon/tests/test_interfaces.py index 2d233f85..17d80d8e 100644 --- a/qsirecon/tests/test_interfaces.py +++ b/qsirecon/tests/test_interfaces.py @@ -6,11 +6,11 @@ import numpy as np import pytest -from qsirecon.interfaces.gradients import GradientSelect -from qsirecon.tests.utils import download_test_data +from qsirecon.data import load as load_data +from qsirecon.interfaces.gradients import GradientSelect, _find_shells +from qsirecon.tests.utils import download_test_data, get_test_data_path -@pytest.mark.integration @pytest.mark.interfaces def test_shell_selection(data_dir, working_dir): """Run reconstruction workflow tests.""" @@ -20,19 +20,17 @@ def test_shell_selection(data_dir, working_dir): dataset_dir = Path(dataset_dir) / "multishell_output" / "qsiprep" data_stem = str(dataset_dir / "sub-ABCD" / "dwi" / dwi_prefix) - # Test the actual data (including a nifti) + # Test actual data (including a nifti) grad_select = GradientSelect( dwi_file=data_stem + ".nii.gz", bval_file=data_stem + ".bval", bvec_file=data_stem + ".bvec", b_file=data_stem + ".b", - btable_file=data_stem + ".b", bval_distance_cutoff=100, expected_n_input_shells=5, requested_shells=[0, "lowest", "highest"], ) grad_select.run() - correct_n = 73 sel_nii = nb.load(dwi_prefix + "_selected.nii.gz") assert sel_nii.shape[3] == correct_n @@ -45,3 +43,28 @@ def test_shell_selection(data_dir, working_dir): # There is no btable from this dataset because it was created # before those were written in the outputs. assert not Path(dwi_prefix + "_selected.txt").exists() + + # Check some other sequences we might run into + # ABCD has 4 shells + b=0s + abcd_bval = np.loadtxt(load_data("schemes/ABCD.bval")) + abcd_df = _find_shells(abcd_bval, 100) + assert abcd_df["shell_num"].max() == 5.0 + + # HCP has 3 shells + b=0s + hcp_bval = np.loadtxt(load_data("schemes/HCP.bval")) + hcp_df = _find_shells(hcp_bval, 100) + assert hcp_df["shell_num"].max() == 4.0 + + # DSIQ5 should raise an exception + dsi_bval = np.loadtxt(load_data("schemes/DSIQ5.bval")) + with pytest.raises(Exception): + _find_shells(dsi_bval, 100) + + # Some other assorted test schemes that should fail + nonshelled_schemes = ["HASC55-1", "HASC55-2", "HASC92", "RAND57"] + for scheme in nonshelled_schemes: + bval_file = Path(get_test_data_path()) / f"{scheme}.bval" + test_bvals = np.loadtxt(bval_file) + + with pytest.raises(Exception): + _find_shells(test_bvals, 100) From 1f2587cfa87014ada48d9913d7903fe38a9bea92 Mon Sep 17 00:00:00 2001 From: Matt Cieslak Date: Thu, 5 Dec 2024 17:30:57 -0500 Subject: [PATCH 10/13] add Q7 to tests. --- qsirecon/interfaces/gradients.py | 3 ++- qsirecon/tests/test_interfaces.py | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/qsirecon/interfaces/gradients.py b/qsirecon/interfaces/gradients.py index 88ee54cd..9b2e45c2 100644 --- a/qsirecon/interfaces/gradients.py +++ b/qsirecon/interfaces/gradients.py @@ -176,9 +176,10 @@ def _find_shells(bvals, max_distance): ).fit(X) shells = agg_cluster.labels_ + # Introduce a new check: score = silhouette_score(X, shells) if score < 0.8: - print("Silhouette score is low. Is this is a DSI scheme?") + raise Exception("Silhouette score is low. Is this is a DSI scheme?") # Do the same check as mrtrix max_shells = np.sqrt(np.sum(bvals > max_distance)) diff --git a/qsirecon/tests/test_interfaces.py b/qsirecon/tests/test_interfaces.py index 17d80d8e..c09466f1 100644 --- a/qsirecon/tests/test_interfaces.py +++ b/qsirecon/tests/test_interfaces.py @@ -61,7 +61,7 @@ def test_shell_selection(data_dir, working_dir): _find_shells(dsi_bval, 100) # Some other assorted test schemes that should fail - nonshelled_schemes = ["HASC55-1", "HASC55-2", "HASC92", "RAND57"] + nonshelled_schemes = ["HASC55-1", "HASC55-2", "HASC92", "RAND57", "Q7"] for scheme in nonshelled_schemes: bval_file = Path(get_test_data_path()) / f"{scheme}.bval" test_bvals = np.loadtxt(bval_file) From ede492afa8a53f309a537d4bc704509144c5b789 Mon Sep 17 00:00:00 2001 From: Matt Cieslak Date: Thu, 5 Dec 2024 17:48:41 -0500 Subject: [PATCH 11/13] set default none for btable and b --- .circleci/config.yml | 4 ++-- qsirecon/interfaces/gradients.py | 2 ++ qsirecon/tests/test_cli.py | 2 +- 3 files changed, 5 insertions(+), 3 deletions(-) diff --git a/.circleci/config.yml b/.circleci/config.yml index dc9b10b4..937bde1a 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -468,7 +468,7 @@ jobs: steps: - checkout - restore_cache: - key: multishell_output-01 + key: singleshell_output-01 - run: *runinstall - run: name: Test the DIPY recon workflows @@ -688,7 +688,7 @@ workflows: - Recon_AMICO: requires: - - download_multishell_output + - download_singleshell_output filters: tags: only: /.*/ diff --git a/qsirecon/interfaces/gradients.py b/qsirecon/interfaces/gradients.py index 9b2e45c2..ea89b439 100644 --- a/qsirecon/interfaces/gradients.py +++ b/qsirecon/interfaces/gradients.py @@ -382,11 +382,13 @@ def subset_dwi( nim.index_img(original_nifti, indices).to_filename(new_nifti) # If there is a b file, subset and write it + new_b = None if original_b is not None: new_b = fname_presuffix(original_b, newpath=newdir, suffix=suffix) _select_lines(original_b, new_b, indices) # If there is a dsi studio btable file, subset and write it + new_btable = None if original_btable is not None: new_btable = fname_presuffix(original_btable, newpath=newdir, suffix=suffix) _select_lines(original_btable, new_btable, indices) diff --git a/qsirecon/tests/test_cli.py b/qsirecon/tests/test_cli.py index a4cb4612..0733259c 100644 --- a/qsirecon/tests/test_cli.py +++ b/qsirecon/tests/test_cli.py @@ -312,7 +312,7 @@ def test_amico_noddi(data_dir, output_dir, working_dir): """ TEST_NAME = "amico_noddi" - dataset_dir = download_test_data("multishell_output", data_dir) + dataset_dir = download_test_data("singleshell_output", data_dir) # XXX: Having to modify dataset_dirs is suboptimal. dataset_dir = os.path.join(dataset_dir, "qsiprep") out_dir = os.path.join(output_dir, TEST_NAME) From ce3fb25b1df1a1166596f71cfec9e0ef816b7c90 Mon Sep 17 00:00:00 2001 From: Matt Cieslak Date: Fri, 6 Dec 2024 13:41:24 -0500 Subject: [PATCH 12/13] improve tests --- qsirecon/tests/test_interfaces.py | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/qsirecon/tests/test_interfaces.py b/qsirecon/tests/test_interfaces.py index c09466f1..81ee2dbb 100644 --- a/qsirecon/tests/test_interfaces.py +++ b/qsirecon/tests/test_interfaces.py @@ -12,7 +12,7 @@ @pytest.mark.interfaces -def test_shell_selection(data_dir, working_dir): +def test_shell_selection(data_dir): """Run reconstruction workflow tests.""" dwi_prefix = "sub-ABCD_acq-10per000_space-T1w_desc-preproc_dwi" data_dir = "/home/matt/projects/qsirecon/.circleci/data" @@ -44,6 +44,9 @@ def test_shell_selection(data_dir, working_dir): # before those were written in the outputs. assert not Path(dwi_prefix + "_selected.txt").exists() + +@pytest.mark.interfaces +def test_real_shells(): # Check some other sequences we might run into # ABCD has 4 shells + b=0s abcd_bval = np.loadtxt(load_data("schemes/ABCD.bval")) @@ -57,14 +60,21 @@ def test_shell_selection(data_dir, working_dir): # DSIQ5 should raise an exception dsi_bval = np.loadtxt(load_data("schemes/DSIQ5.bval")) - with pytest.raises(Exception): + with pytest.raises(Exception, match="Too many possible shells detected."): _find_shells(dsi_bval, 100) # Some other assorted test schemes that should fail - nonshelled_schemes = ["HASC55-1", "HASC55-2", "HASC92", "RAND57", "Q7"] + nonshelled_schemes = ["HASC55-1", "HASC55-2", "HASC92", "RAND57"] for scheme in nonshelled_schemes: bval_file = Path(get_test_data_path()) / f"{scheme}.bval" test_bvals = np.loadtxt(bval_file) - with pytest.raises(Exception): + with pytest.raises(Exception, match="Too many possible shells detected."): _find_shells(test_bvals, 100) + + # DSIQ7 actually _passes_ the mrtrix test, but failes silhouette + bval_file = Path(get_test_data_path()) / "Q7.bval" + test_bvals = np.loadtxt(bval_file) + + with pytest.raises(Exception, match="Silhouette score is low. Is this is a DSI scheme?"): + _find_shells(test_bvals, 100) From b43b9c260443ca328621733309553f316bc70bc6 Mon Sep 17 00:00:00 2001 From: Matt Cieslak Date: Fri, 6 Dec 2024 13:49:24 -0500 Subject: [PATCH 13/13] Update qsirecon/interfaces/gradients.py Co-authored-by: Taylor Salo --- qsirecon/interfaces/gradients.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/qsirecon/interfaces/gradients.py b/qsirecon/interfaces/gradients.py index ea89b439..e1dbdd17 100644 --- a/qsirecon/interfaces/gradients.py +++ b/qsirecon/interfaces/gradients.py @@ -145,8 +145,8 @@ def _parse_shell_selection(requested_bvals, bval_df, max_distance): ] if len(ok_shells.unique()) > 1: raise Exception( - f"Unable to unambiguously select b={requested_bval}." - f"instead select from {bval_df.assigned_shells.unique().tolist()}" + f"Unable to unambiguously select b={requested_bval}. " + f"Instead select from {bval_df.assigned_shells.unique().tolist()}" ) numeric_bvals.append(ok_shells.iloc[0])