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/data/pipelines/mrtrix_singleshell_ss3t_ACT-fast.yaml b/qsirecon/data/pipelines/mrtrix_singleshell_ss3t_ACT-fast.yaml index 2d6ffe4d..aba3306d 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: 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_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: 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 2be923b2..e1dbdd17 100644 --- a/qsirecon/interfaces/gradients.py +++ b/qsirecon/interfaces/gradients.py @@ -2,12 +2,12 @@ import logging -import nibabel as nb import numpy as np from nilearn import image as nim from nipype.interfaces.base import ( BaseInterfaceInputSpec, File, + InputMultiObject, SimpleInterface, TraitedSpec, isdefined, @@ -18,33 +18,214 @@ 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) + 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( + 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) + + +class GradientSelect(SimpleInterface): + input_spec = _GradientSelectInputSpec + output_spec = _GradientSelectOutputSpec + + def _run_interface(self, runtime): + """Find shells in the input data and select""" + bvals = np.loadtxt(self.inputs.bval_file) + + 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, max_distance) + + # Make sure the shells are unique (no overlap when accounting for allowed distances) + 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 " + "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"Unable to unambiguously select b={requested_bval}. " + f"Instead select from {bval_df.assigned_shells.unique().tolist()}" + ) + + numeric_bvals.append(ok_shells.iloc[0]) + + # 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 np.array(sorted(numeric_bvals)) + + +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(X) + shells = agg_cluster.labels_ + + # Introduce a new check: + score = silhouette_score(X, shells) + if score < 0.8: + 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)) + 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( + shell_df["assignment"].tolist(), shell_df["bvalue"].tolist() + ) + bval_df["shell_num"] = bval_df["assigned_shell"].rank(method="dense") + bval_df.drop(columns=["assignment"], inplace=True) + + 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() -class RemoveDuplicatesOutputSpec(TraitedSpec): +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) 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 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 @@ -77,32 +258,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 @@ -176,3 +350,56 @@ 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 + 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) + 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/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_interfaces.py b/qsirecon/tests/test_interfaces.py new file mode 100644 index 00000000..81ee2dbb --- /dev/null +++ b/qsirecon/tests/test_interfaces.py @@ -0,0 +1,80 @@ +"""Command-line interface tests.""" + +from pathlib import Path + +import nibabel as nb +import numpy as np +import pytest + +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.interfaces +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" + 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 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", + 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() + + +@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")) + 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, 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"] + 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, 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) 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/utils.py b/qsirecon/workflows/recon/utils.py index 1de9feb4..1196f762 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 GradientSelect, RemoveDuplicates from ...interfaces.interchange import recon_workflow_input_fields from ...interfaces.mrtrix import MRTrixGradientTable from ...interfaces.recon_scalars import OrganizeScalarData @@ -72,10 +72,49 @@ 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 + + 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'), + ('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