From 97cc2b6ecbbd724461d63fec13c61c231376ce5b Mon Sep 17 00:00:00 2001 From: Arshiyas Date: Wed, 28 Dec 2022 14:21:23 +1100 Subject: [PATCH 01/24] dwi2mask: added algorithm that iteratively performs masking, RF estimation, CSD and bias field removal --- docs/reference/commands/dwi2mask.rst | 98 ++++++++++++- lib/mrtrix3/dwi2mask/iterative.py | 211 +++++++++++++++++++++++++++ 2 files changed, 308 insertions(+), 1 deletion(-) create mode 100644 lib/mrtrix3/dwi2mask/iterative.py diff --git a/docs/reference/commands/dwi2mask.rst b/docs/reference/commands/dwi2mask.rst index 64a60570ab..34d0558eb3 100644 --- a/docs/reference/commands/dwi2mask.rst +++ b/docs/reference/commands/dwi2mask.rst @@ -15,7 +15,7 @@ Usage dwi2mask algorithm [ options ] ... -- *algorithm*: Select the algorithm to be used to complete the script operation; additional details and options become available once an algorithm is nominated. Options are: 3dautomask, ants, b02template, consensus, fslbet, hdbet, legacy, mean, synthstrip, trace +- *algorithm*: Select the algorithm to be used to complete the script operation; additional details and options become available once an algorithm is nominated. Options are: 3dautomask, ants, b02template, consensus, fslbet, hdbet, iterative, legacy, mean, synthstrip, trace Description ----------- @@ -685,6 +685,102 @@ See the Mozilla Public License v. 2.0 for more details. For more details, see http://www.mrtrix.org/. +.. _dwi2mask_iterative: + +dwi2mask iterative +================== + +Synopsis +-------- + +Iteratively perform masking, RF estimation, CSD, bias field removal, and mask revision to derive a DWI brain mask + +Usage +----- + +:: + + dwi2mask iterative input output [ options ] + +- *input*: The input DWI series +- *output*: The output mask image + +Description +----------- + +DWI brain masking may be inaccurate due to residual bias field. This script first derives an initial brain mask using the legacy MRtrix3 dwi2mask heuristic (based on thresholded trace images), and then performs response function estimation, multi-tissue CSD (with a lower lmax than the dwi2fod default, for higher speed), and mtnormalise to remove any bias field, before the DWI brain mask is recalculated. These steps are performed iteratively, since each step may influence the others. + +Options +------- + +Options specific to the "iterative" algorithm +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-max_iters** the maximum number of iterations. The default is 2iterations since for some problematic data more iterations may lead to the masks diverging, but more iterations can lead to a better mask. + +- **-lmax** the maximum spherical harmonic order for the output FODs. The value is passed to the dwi2fod command and should be provided in the format which it expects (Default is "4,0,0" for multi-shell and "4,0" for single-shell data) + +Options for importing the diffusion gradient table +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-grad** Provide the diffusion gradient table in MRtrix format + +- **-fslgrad bvecs bvals** Provide the diffusion gradient table in FSL bvecs/bvals format + +Additional standard options for Python scripts +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-nocleanup** do not delete intermediate files during script execution, and do not delete scratch directory at script completion. + +- **-scratch /path/to/scratch/** manually specify the path in which to generate the scratch directory. + +- **-continue ** continue the script from a previous execution; must provide the scratch directory path, and the name of the last successfully-generated file. + +Standard options +^^^^^^^^^^^^^^^^ + +- **-info** display information messages. + +- **-quiet** do not display information messages or progress status. Alternatively, this can be achieved by setting the MRTRIX_QUIET environment variable to a non-empty string. + +- **-debug** display debugging messages. + +- **-force** force overwrite of output files. + +- **-nthreads number** use this number of threads in multi-threaded applications (set to 0 to disable multi-threading). + +- **-config key value** *(multiple uses permitted)* temporarily set the value of an MRtrix config file entry. + +- **-help** display this information page and exit. + +- **-version** display version information and exit. + +References +^^^^^^^^^^ + +Tournier, J.-D.; Smith, R. E.; Raffelt, D.; Tabbara, R.; Dhollander, T.; Pietsch, M.; Christiaens, D.; Jeurissen, B.; Yeh, C.-H. & Connelly, A. MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation. NeuroImage, 2019, 202, 116137 + +-------------- + + + +**Author:** Robert E. Smith (robert.smith@florey.edu.au) and Arshiya Sangchooli (asangchooli@student.unimelb.edu.au) + +**Copyright:** Copyright (c) 2008-2022 the MRtrix3 contributors. + +This Source Code Form is subject to the terms of the Mozilla Public +License, v. 2.0. If a copy of the MPL was not distributed with this +file, You can obtain one at http://mozilla.org/MPL/2.0/. + +Covered Software is provided under this License on an "as is" +basis, without warranty of any kind, either expressed, implied, or +statutory, including, without limitation, warranties that the +Covered Software is free of defects, merchantable, fit for a +particular purpose or non-infringing. +See the Mozilla Public License v. 2.0 for more details. + +For more details, see http://www.mrtrix.org/. + .. _dwi2mask_legacy: dwi2mask legacy diff --git a/lib/mrtrix3/dwi2mask/iterative.py b/lib/mrtrix3/dwi2mask/iterative.py new file mode 100644 index 0000000000..9fe5c49fa8 --- /dev/null +++ b/lib/mrtrix3/dwi2mask/iterative.py @@ -0,0 +1,211 @@ +# Copyright (c) 2008-2022 the MRtrix3 contributors. +# +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at http://mozilla.org/MPL/2.0/. +# +# Covered Software is provided under this License on an "as is" +# basis, without warranty of any kind, either expressed, implied, or +# statutory, including, without limitation, warranties that the +# Covered Software is free of defects, merchantable, fit for a +# particular purpose or non-infringing. +# See the Mozilla Public License v. 2.0 for more details. +# +# For more details, see http://www.mrtrix.org/. + +import math +from mrtrix3 import app, run, image, matrix + +DWIBIASCORRECT_MAX_ITERS = 2 + +# We use a threshold on the balanced tissue sum image as a replacement of dwi2mask +# within the iterative bias field correction / brain masking loop +TISSUESUM_THRESHOLD = 0.5 / math.sqrt(4.0 * math.pi) + +LMAXES_MULTI = '4,0,0' +LMAXES_SINGLE = '4,0' + +def usage(base_parser, subparsers): #pylint: disable=unused-variable + parser = subparsers.add_parser('iterative', parents=[base_parser]) + parser.set_author('Robert E. Smith (robert.smith@florey.edu.au) and Arshiya Sangchooli (asangchooli@student.unimelb.edu.au)') + parser.set_synopsis('Iteratively perform masking, RF estimation, CSD, bias field removal, and mask revision to derive a DWI brain mask') + parser.add_description('DWI brain masking may be inaccurate due to residual bias field. This script first derives an initial ' + 'brain mask using the legacy MRtrix3 dwi2mask heuristic (based on thresholded trace images), and then ' + 'performs response function estimation, multi-tissue CSD (with a lower lmax than the dwi2fod default, for ' + 'higher speed), and mtnormalise to remove any bias field, before the DWI brain mask is recalculated. ' + 'These steps are performed iteratively, since each step may influence the others.') + parser.add_argument('input', help='The input DWI series') + parser.add_argument('output', help='The output mask image') + options = parser.add_argument_group('Options specific to the "iterative" algorithm') + options.add_argument('-max_iters', + type=int, + default=DWIBIASCORRECT_MAX_ITERS, + help='the maximum number of iterations. The default is ' + str(DWIBIASCORRECT_MAX_ITERS) + + 'iterations since for some problematic data more iterations may lead to the masks ' + 'diverging, but more iterations can lead to a better mask.') + options.add_argument('-lmax', + type=str, + help='the maximum spherical harmonic order for the output FODs. The value is passed to ' + 'the dwi2fod command and should be provided in the format which it expects ' + '(Default is "' + str(LMAXES_MULTI) + '" for multi-shell and "' + str(LMAXES_SINGLE) + '" for single-shell data)') + + +def get_inputs(): #pylint: disable=unused-variable + pass + + +def needs_mean_bzero(): #pylint: disable=unused-variable + return False + + + + + +def execute(): #pylint: disable=unused-variable + # Determine whether we are working with single-shell or multi-shell data + bvalues = [ + int(round(float(value))) + for value in image.mrinfo('input.mif', 'shell_bvalues') \ + .strip().split()] + multishell = (len(bvalues) > 2) + + # Step 1: Initial DWI brain mask + dwi_mask_image = 'dwi_mask_init.mif' + app.console('Performing intial DWI brain masking with the legacy MRtrix3 algorithm') + run.command('dwi2mask legacy input.mif ' + dwi_mask_image) + + # Step 2: Combined RF estimation / CSD / mtnormalise / mask revision + class Tissue(object): #pylint: disable=useless-object-inheritance + def __init__(self, name, index): + self.name = name + iter_string = '_iter' + str(index) + self.tissue_rf = 'response_' + name + iter_string + '.txt' + self.fod_init = 'FODinit_' + name + iter_string + '.mif' + self.fod_norm = 'FODnorm_' + name + iter_string + '.mif' + + app.console('Commencing iterative DWI bias field correction and brain masking ' + 'with a maximum of ' + str(app.ARGS.max_iters) + ' iterations') + + dwi_image = 'input.mif' + + for iteration in range(0, app.ARGS.max_iters): + iter_string = '_iter' + str(iteration+1) + tissues = [Tissue('WM', iteration), + Tissue('GM', iteration), + Tissue('CSF', iteration)] + app.console('Iteration ' + str(iteration) + ', ' + + 'estimating response function using brain mask...') + run.command('dwi2response dhollander ' + + dwi_image + + ' -mask ' + + dwi_mask_image + + ' ' + + ' '.join(tissue.tissue_rf for tissue in tissues)) + + # Remove GM if we can't deal with it + if app.ARGS.lmax: + lmaxes = app.ARGS.lmax + else: + lmaxes = LMAXES_MULTI + if not multishell: + app.cleanup(tissues[1].tissue_rf) + tissues = tissues[::2] + lmaxes = LMAXES_SINGLE + + app.console('Iteration ' + str(iteration) + ', ' + + 'FOD estimation with lmaxes ' + lmaxes + '...') + run.command('dwi2fod msmt_csd ' + + dwi_image + + ' -lmax ' + lmaxes + + ' ' + + ' '.join(tissue.tissue_rf + ' ' + tissue.fod_init + for tissue in tissues)) + + app.console('Iteration ' + str(iteration) + ', ' + + 'running mtnormalise...') + field_path = 'field' + iter_string + '.mif' + factors_path = 'factors' + iter_string + '.txt' + run.command('maskfilter ' + dwi_mask_image + ' erode - |' + + ' mtnormalise -mask - -balanced' + + ' -check_norm ' + field_path + + ' -check_factors ' + factors_path + + ' ' + + ' '.join(tissue.fod_init + ' ' + tissue.fod_norm + for tissue in tissues)) + app.cleanup([tissue.fod_init for tissue in tissues]) + + app.console('Iteration ' + str(iteration) + ', ' + + 'applying estimated bias field and appropiate scaling factor...') + csf_rf = matrix.load_matrix(tissues[-1].tissue_rf) + csf_rf_bzero_lzero = csf_rf[0][0] + app.cleanup([tissue.tissue_rf for tissue in tissues]) + balance_factors = matrix.load_vector(factors_path) + csf_balance_factor = balance_factors[-1] + app.cleanup(factors_path) + scale_multiplier = (1000.0 * math.sqrt(4.0*math.pi)) / \ + (csf_rf_bzero_lzero / csf_balance_factor) + new_dwi_image = 'dwi' + iter_string + '.mif' + run.command('mrcalc ' + dwi_image + ' ' + + field_path + ' -div ' + + str(scale_multiplier) + ' -mult ' + + new_dwi_image) + app.cleanup(field_path) + app.cleanup(dwi_image) + dwi_image = new_dwi_image + + app.console('Iteration ' + str(iteration) + ', ' + + 'revising brain mask...') + new_dwi_mask_image = 'dwi_mask' + iter_string + '.mif' + run.command('mrconvert ' + + tissues[0].fod_norm + + ' -coord 3 0 - |' + + ' mrmath - ' + + ' '.join(tissue.fod_norm for tissue in tissues[1:]) + + ' sum - |' + + ' mrthreshold - -abs ' + + str(TISSUESUM_THRESHOLD) + + ' - |' + + ' maskfilter - connect -largest - |' + + ' mrcalc 1 - -sub - -datatype bit |' + + ' maskfilter - connect -largest - |' + + ' mrcalc 1 - -sub - -datatype bit |' + + ' maskfilter - clean - |' + + ' mrcalc - input_pos_mask.mif -mult ' + + new_dwi_mask_image + + ' -datatype bit') + app.cleanup([tissue.fod_norm for tissue in tissues]) + + app.console('Iteration ' + str(iteration) + ', ' + + 'comparing the new mask with the previous iteration...') + dwi_old_mask_count = image.statistics(dwi_mask_image, + mask=dwi_mask_image).count + dwi_new_mask_count = image.statistics(new_dwi_mask_image, + mask=new_dwi_mask_image).count + app.debug('Old mask: ' + str(dwi_old_mask_count)) + app.debug('New mask: ' + str(dwi_new_mask_count)) + dwi_mask_overlap_image = 'dwi_mask_overlap' + iter_string + '.mif' + run.command('mrcalc ' + + dwi_mask_image + + ' ' + + new_dwi_mask_image + + ' -mult ' + + dwi_mask_overlap_image) + app.cleanup(dwi_mask_image) + dwi_mask_image = new_dwi_mask_image + mask_overlap_count = image.statistics(dwi_mask_overlap_image, + mask=dwi_mask_overlap_image).count + app.debug('Mask overlap: ' + str(mask_overlap_count)) + dice_coefficient = 2.0 * mask_overlap_count / \ + (dwi_old_mask_count + dwi_new_mask_count) + app.console('Iteration ' + str(iteration) + + ' dice coefficient: ' + str(dice_coefficient) + '...') + if dice_coefficient > (1.0 - 1e-3): + app.console('Exiting iterative loop due to mask convergence') + break + + # Processing completed; export + app.console('Processing completed after ' + + str(iteration + 1) + + ' iterations, writing results to output directory') + + return dwi_mask_image From 2282ebee694f3a478f600ce5929d844eb3d369ee Mon Sep 17 00:00:00 2001 From: Arshiyas Date: Mon, 20 Feb 2023 14:37:19 +1100 Subject: [PATCH 02/24] Added additional algorithms --- lib/mrtrix3/dwi2mask/mtnorm.py | 125 ++++++++++++++++++++++ lib/mrtrix3/dwibiascorrect/mtnorm.py | 111 ++++++++++++++++++++ lib/mrtrix3/dwinormalise/mtnorm.py | 151 +++++++++++++++++++++++++++ 3 files changed, 387 insertions(+) create mode 100644 lib/mrtrix3/dwi2mask/mtnorm.py create mode 100644 lib/mrtrix3/dwibiascorrect/mtnorm.py create mode 100644 lib/mrtrix3/dwinormalise/mtnorm.py diff --git a/lib/mrtrix3/dwi2mask/mtnorm.py b/lib/mrtrix3/dwi2mask/mtnorm.py new file mode 100644 index 0000000000..8c52459724 --- /dev/null +++ b/lib/mrtrix3/dwi2mask/mtnorm.py @@ -0,0 +1,125 @@ +# Copyright (c) 2008-2022 the MRtrix3 contributors. +# +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at http://mozilla.org/MPL/2.0/. +# +# Covered Software is provided under this License on an "as is" +# basis, without warranty of any kind, either expressed, implied, or +# statutory, including, without limitation, warranties that the +# Covered Software is free of defects, merchantable, fit for a +# particular purpose or non-infringing. +# See the Mozilla Public License v. 2.0 for more details. +# +# For more details, see http://www.mrtrix.org/. + +import math +from mrtrix3 import app, run, image + +LMAXES_MULTI = '4,0,0' +LMAXES_SINGLE = '4,0' + +def usage(base_parser, subparsers): #pylint: disable=unused-variable + parser = subparsers.add_parser('mtnorm', parents=[base_parser]) + parser.set_author('Robert E. Smith (robert.smith@florey.edu.au) and Arshiya Sangchooli (asangchooli@student.unimelb.edu.au)') + parser.set_synopsis('Derives a DWI brain mask by calculating and then thresholding a tissue density image') + parser.add_description('This script first derives an initial brain mask using the legacy MRtrix3 dwi2mask heuristic ' + '(based on thresholded trace images), and then performs response function estimation and multi-tissue CSD ' + '(with a lower lmax than the dwi2fod default, for higher speed), before summing the derived tissue ODFs and ' + 'thresholding the resulting tissue density image to derive a DWI brain mask.') + parser.add_argument('input', help='The input DWI series') + parser.add_argument('output', help='The output mask image') + options = parser.add_argument_group('Options specific to the "mtnorm" algorithm') + options.add_argument('-threshold', + type=float, + default=0.5, + help='the threshold on the total tissue density sum image used to derive the brain mask. the default is 0.5') + options.add_argument('-lmax', + type=str, + help='the maximum spherical harmonic order for the output FODs. The value is passed to ' + 'the dwi2fod command and should be provided in the format which it expects ' + '(Default is "' + str(LMAXES_MULTI) + '" for multi-shell and "' + str(LMAXES_SINGLE) + '" for single-shell data)') + +def get_inputs(): #pylint: disable=unused-variable + pass + + +def needs_mean_bzero(): #pylint: disable=unused-variable + return False + + +def execute(): #pylint: disable=unused-variable + # Determine whether we are working with single-shell or multi-shell data + bvalues = [ + int(round(float(value))) + for value in image.mrinfo('input.mif', 'shell_bvalues') \ + .strip().split()] + multishell = (len(bvalues) > 2) + + # Step 1: Initial DWI brain mask + dwi_mask_image = 'dwi_mask_init.mif' + app.debug('Performing intial DWI brain masking with the legacy MRtrix3 algorithm') + run.command('dwi2mask legacy input.mif ' + dwi_mask_image) + + # Step 2: RF estimation / CSD / mtnormalise / mask revision + class Tissue(object): #pylint: disable=useless-object-inheritance + def __init__(self, name): + self.name = name + self.tissue_rf = 'response_' + name + '.txt' + self.fod = 'FOD_' + name + '.mif' + + dwi_image = 'input.mif' + tissues = [Tissue('WM'), Tissue('GM'), Tissue('CSF')] + + app.debug('Estimating response function using initial brain mask...') + run.command('dwi2response dhollander ' + + dwi_image + + ' -mask ' + + dwi_mask_image + + ' ' + + ' '.join(tissue.tissue_rf for tissue in tissues)) + + # Remove GM if we can't deal with it + if app.ARGS.lmax: + lmaxes = app.ARGS.lmax + else: + lmaxes = LMAXES_MULTI + if not multishell: + app.cleanup(tissues[1].tissue_rf) + tissues = tissues[::2] + lmaxes = LMAXES_SINGLE + + app.debug('FOD estimation with lmaxes ' + lmaxes + '...') + run.command('dwi2fod msmt_csd ' + + dwi_image + + ' -lmax ' + lmaxes + + ' ' + + ' '.join(tissue.tissue_rf + ' ' + tissue.fod + for tissue in tissues)) + + app.debug('Deriving final brain mask by thresholding tissue sum image...') + new_dwi_mask_image = 'dwi_mask' + '.mif' + + TISSUESUM_THRESHOLD = app.ARGS.threshold / math.sqrt(4.0 * math.pi) + + run.command('mrconvert ' + + tissues[0].fod + + ' -coord 3 0 - |' + + ' mrmath - ' + + ' '.join(tissue.fod for tissue in tissues[1:]) + + ' sum - |' + + ' mrthreshold - -abs ' + + str(TISSUESUM_THRESHOLD) + + ' - |' + + ' maskfilter - connect -largest - |' + + ' mrcalc 1 - -sub - -datatype bit |' + + ' maskfilter - connect -largest - |' + + ' mrcalc 1 - -sub - -datatype bit |' + + ' maskfilter - clean - |' + + ' mrcalc - input_pos_mask.mif -mult ' + + new_dwi_mask_image + + ' -datatype bit') + app.cleanup([tissue.fod for tissue in tissues]) + + return dwi_mask_image + diff --git a/lib/mrtrix3/dwibiascorrect/mtnorm.py b/lib/mrtrix3/dwibiascorrect/mtnorm.py new file mode 100644 index 0000000000..35874521dd --- /dev/null +++ b/lib/mrtrix3/dwibiascorrect/mtnorm.py @@ -0,0 +1,111 @@ +# Copyright (c) 2008-2022 the MRtrix3 contributors. +# +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at http://mozilla.org/MPL/2.0/. +# +# Covered Software is provided under this License on an "as is" +# basis, without warranty of any kind, either expressed, implied, or +# statutory, including, without limitation, warranties that the +# Covered Software is free of defects, merchantable, fit for a +# particular purpose or non-infringing. +# See the Mozilla Public License v. 2.0 for more details. +# +# For more details, see http://www.mrtrix.org/. + +import math +from mrtrix3 import app, run, image, path + + +LMAXES_MULTI = '4,0,0' +LMAXES_SINGLE = '4,0' + +def usage(base_parser, subparsers): #pylint: disable=unused-variable + parser = subparsers.add_parser('mtnorm', parents=[base_parser]) + parser.set_author('Robert E. Smith (robert.smith@florey.edu.au) and Arshiya Sangchooli (asangchooli@student.unimelb.edu.au)') + parser.set_synopsis('Performs DWI bias field correction') + parser.add_description('This script first performs response function estimation and multi-tissue CSD within a provided or derived ' + 'brain mask (with a lower lmax than the dwi2fod default, for higher speed) before calling mtnormalise and ' + 'utilizing the final estimated spatially varying intensity level used for normalisation to correct bias fields') + parser.add_argument('input', help='The input image series to be corrected') + parser.add_argument('output', help='The output corrected image series') + options = parser.add_argument_group('Options specific to the "mtnorm" algorithm') + options.add_argument('-lmax', + type=str, + help='the maximum spherical harmonic order for the output FODs. The value is passed to ' + 'the dwi2fod command and should be provided in the format which it expects ' + '(Default is "' + str(LMAXES_MULTI) + '" for multi-shell and "' + str(LMAXES_SINGLE) + '" for single-shell data)') + + +def check_output_paths(): #pylint: disable=unused-variable + pass + +def get_inputs(): #pylint: disable=unused-variable + pass + + +def execute(): #pylint: disable=unused-variable + # Determine whether we are working with single-shell or multi-shell data + bvalues = [ + int(round(float(value))) + for value in image.mrinfo('in.mif', 'shell_bvalues') \ + .strip().split()] + multishell = (len(bvalues) > 2) + + # RF estimation and multi-tissue CSD + class Tissue(object): #pylint: disable=useless-object-inheritance + def __init__(self, name): + self.name = name + self.tissue_rf = 'response_' + name + '.txt' + self.fod = 'FOD_' + name + '.mif' + self.fod_norm = 'FODnorm_' + name + '.mif' + + tissues = [Tissue('WM'), Tissue('GM'), Tissue('CSF')] + + app.debug('Estimating response function using brain mask...') + run.command('dwi2response dhollander ' + + 'in.mif ' + + '-mask mask.mif ' + + ' '.join(tissue.tissue_rf for tissue in tissues)) + + # Remove GM if we can't deal with it + if app.ARGS.lmax: + lmaxes = app.ARGS.lmax + else: + lmaxes = LMAXES_MULTI + if not multishell: + app.cleanup(tissues[1].tissue_rf) + tissues = tissues[::2] + lmaxes = LMAXES_SINGLE + + app.debug('FOD estimation with lmaxes ' + lmaxes + '...') + run.command('dwi2fod msmt_csd ' + + 'in.mif' + + ' -lmax ' + lmaxes + + ' ' + + ' '.join(tissue.tissue_rf + ' ' + tissue.fod + for tissue in tissues)) + + app.debug('correcting bias field...') + run.command('maskfilter' + + ' mask.mif' + + ' erode - |' + + ' mtnormalise -mask - -balanced' + + ' -check_norm field.mif ' + + ' '.join(tissue.fod + ' ' + tissue.fod_norm + for tissue in tissues)) + app.cleanup([tissue.fod for tissue in tissues]) + app.cleanup([tissue.fod_norm for tissue in tissues]) + app.cleanup([tissue.tissue_rf for tissue in tissues]) + + run.command('mrcalc ' + + 'in.mif ' + + 'field.mif ' + + '-div ' + + path.from_user(app.ARGS.output), + force=app.FORCE_OVERWRITE) + + if app.ARGS.bias: + run.command('mrconvert field.mif ' + + path.from_user(app.ARGS.bias), + mrconvert_keyval=path.from_user(app.ARGS.input, False), force=app.FORCE_OVERWRITE) diff --git a/lib/mrtrix3/dwinormalise/mtnorm.py b/lib/mrtrix3/dwinormalise/mtnorm.py new file mode 100644 index 0000000000..30b9bb3948 --- /dev/null +++ b/lib/mrtrix3/dwinormalise/mtnorm.py @@ -0,0 +1,151 @@ +# Copyright (c) 2008-2022 the MRtrix3 contributors. +# +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at http://mozilla.org/MPL/2.0/. +# +# Covered Software is provided under this License on an "as is" +# basis, without warranty of any kind, either expressed, implied, or +# statutory, including, without limitation, warranties that the +# Covered Software is free of defects, merchantable, fit for a +# particular purpose or non-infringing. +# See the Mozilla Public License v. 2.0 for more details. +# +# For more details, see http://www.mrtrix.org/. + +import math +from mrtrix3 import CONFIG, MRtrixError +from mrtrix3 import app, run, image, path, matrix, _version + +DEFAULT_TARGET_INTENSITY=1000 + +LMAXES_MULTI = '4,0,0' +LMAXES_SINGLE = '4,0' + + +def usage(base_parser, subparsers): #pylint: disable=unused-variable + parser = subparsers.add_parser('mtnorm', parents=[base_parser]) + parser.set_author('Robert E. Smith (robert.smith@florey.edu.au) and Arshiya Sangchooli (asangchooli@student.unimelb.edu.au)') + parser.set_synopsis('Performs a global DWI intensity normalisation, normalising mean b=0 CSF signal intensity') + parser.add_description('This script first performs response function estimation and multi-tissue CSD within a provided brain mask' + '(with a lower lmax than the dwi2fod default, for higher speed) before using the global scaling factor derived ' + 'by mtnormalise to scale the intensity of the DWI series such that the mean b=0 CSF signal intensity is ' + 'approximately close to some reference value (by default ' + str(DEFAULT_TARGET_INTENSITY) + ').' + 'Note that if the -mask command-line option is not specified, the MRtrix3 command dwi2mask will automatically be ' + 'called to derive a mask that will be passed to the relevant bias field estimation command. ' + 'More information on mask derivation from DWI data can be found at the following link: \n' + 'https://mrtrix.readthedocs.io/en/' + _version.__tag__ + '/dwi_preprocessing/masking.html') + parser.add_argument('input_dwi', help='The input DWI series') + parser.add_argument('output', help='The normalised DWI series') + options = parser.add_argument_group('Options specific to the "mtnorm" algorithm') + options.add_argument('-target', type=float, default=DEFAULT_TARGET_INTENSITY, + help='the threshold on the balanced tissue sum image used to derive the brain mask. the default is 0.5') + options.add_argument('-lmax', type=str, + help='the maximum spherical harmonic order for the output FODs. The value is passed to ' + 'the dwi2fod command and should be provided in the format which it expects ' + '(Default is "' + str(LMAXES_MULTI) + '" for multi-shell and "' + str(LMAXES_SINGLE) + '" for single-shell data)') + options.add_argument('-mask', metavar='image', help='Manually provide a mask image for normalisation') + app.add_dwgrad_import_options(parser) + +def check_output_paths(): #pylint: disable=unused-variable + app.check_output_path(app.ARGS.output) + + +def execute(): #pylint: disable=unused-variable + # Get input data into the scratch directory + grad_option = '' + if app.ARGS.grad: + grad_option = ' -grad ' + path.from_user(app.ARGS.grad) + elif app.ARGS.fslgrad: + grad_option = ' -fslgrad ' + path.from_user(app.ARGS.fslgrad[0]) + ' ' + path.from_user(app.ARGS.fslgrad[1]) + + app.make_scratch_dir() + app.goto_scratch_dir() + + run.command('mrconvert ' + + path.from_user(app.ARGS.input_dwi) + + ' ' + + path.to_scratch('input.mif') + + grad_option) + + # Get mask into scratch directory + if app.ARGS.mask: + run.command('mrconvert ' + path.from_user(app.ARGS.mask) + ' ' + path.to_scratch('mask.mif')) + if not image.match('input.mif', 'mask.mif', up_to_dim=3): + raise MRtrixError('Provided mask image does not match input DWI') + else: + app.debug('Deriving brain mask...') + run.command('dwi2mask ' + CONFIG.get('Dwi2maskAlgorithm', 'legacy') + ' input.mif mask.mif') + + # Determine whether we are working with single-shell or multi-shell data + bvalues = [ + int(round(float(value))) + for value in image.mrinfo('input.mif', 'shell_bvalues') \ + .strip().split()] + multishell = (len(bvalues) > 2) + + # RF estimation and multi-tissue CSD + class Tissue(object): #pylint: disable=useless-object-inheritance + def __init__(self, name): + self.name = name + self.tissue_rf = 'response_' + name + '.txt' + self.fod = 'FOD_' + name + '.mif' + self.fod_norm = 'FODnorm_' + name + '.mif' + + tissues = [Tissue('WM'), Tissue('GM'), Tissue('CSF')] + + app.debug('Estimating response function using initial brain mask...') + run.command('dwi2response dhollander ' + + 'input.mif ' + + '-mask mask.mif ' + + ' '.join(tissue.tissue_rf for tissue in tissues)) + + # Remove GM if we can't deal with it + if app.ARGS.lmax: + lmaxes = app.ARGS.lmax + else: + lmaxes = LMAXES_MULTI + if not multishell: + app.cleanup(tissues[1].tissue_rf) + tissues = tissues[::2] + lmaxes = LMAXES_SINGLE + + app.debug('FOD estimation with lmaxes ' + lmaxes + '...') + run.command('dwi2fod msmt_csd ' + + 'input.mif' + + ' -lmax ' + lmaxes + + ' ' + + ' '.join(tissue.tissue_rf + ' ' + tissue.fod + for tissue in tissues)) + + # Normalisation in brain mask + run.command('maskfilter' + + ' mask.mif' + + ' erode - |' + + ' mtnormalise -mask - -balanced' + + ' -check_factors factors.txt ' + + ' '.join(tissue.fod + ' ' + tissue.fod_norm + for tissue in tissues)) + app.cleanup([tissue.fod for tissue in tissues]) + app.cleanup([tissue.fod_norm for tissue in tissues]) + + app.debug('applying appropiate scaling factor to DWI series...') + csf_rf = matrix.load_matrix(tissues[-1].tissue_rf) + csf_rf_bzero_lzero = csf_rf[0][0] + app.cleanup([tissue.tissue_rf for tissue in tissues]) + balance_factors = matrix.load_vector('factors.txt') + csf_balance_factor = balance_factors[-1] + app.cleanup('factors.txt') + scale_multiplier = (1000.0 * math.sqrt(4.0*math.pi)) / (csf_rf_bzero_lzero / csf_balance_factor) + run.command('mrcalc ' + + 'input.mif ' + + str(scale_multiplier) + + ' -mult ' + + path.from_user(app.ARGS.output), + force=app.FORCE_OVERWRITE) + + + + + + From f592d89704deab1aa0201ef5399b7a27209ddd49 Mon Sep 17 00:00:00 2001 From: Arshiyas Date: Thu, 16 Mar 2023 10:18:56 +1100 Subject: [PATCH 03/24] New command dwibiasnormmask --- bin/dwibiasnormmask | 339 ++++++++++++++++++++++++++++++ lib/mrtrix3/dwi2mask/iterative.py | 211 ------------------- 2 files changed, 339 insertions(+), 211 deletions(-) create mode 100644 bin/dwibiasnormmask delete mode 100644 lib/mrtrix3/dwi2mask/iterative.py diff --git a/bin/dwibiasnormmask b/bin/dwibiasnormmask new file mode 100644 index 0000000000..a0f6cdcca4 --- /dev/null +++ b/bin/dwibiasnormmask @@ -0,0 +1,339 @@ +#!/usr/bin/python3 + +# Copyright (c) 2008-2022 the MRtrix3 contributors. +# +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at http://mozilla.org/MPL/2.0/. +# +# Covered Software is provided under this License on an "as is" +# basis, without warranty of any kind, either expressed, implied, or +# statutory, including, without limitation, warranties that the +# Covered Software is free of defects, merchantable, fit for a +# particular purpose or non-infringing. +# See the Mozilla Public License v. 2.0 for more details. +# +# For more details, see http://www.mrtrix.org/. + +import math + +DWIBIASCORRECT_MAX_ITERS = 2 +LMAXES_MULTI = '4,0,0' +LMAXES_SINGLE = '4,0' + + +def usage(cmdline): #pylint: disable=unused-variable + from mrtrix3 import app #pylint: disable=no-name-in-module, import-outside-toplevel + cmdline.set_author('Robert E. Smith (robert.smith@florey.edu.au) and Arshiya Sangchooli (asangchooli@student.unimelb.edu.au)') + cmdline.set_synopsis('Iteratively perform masking, RF estimation, CSD, bias field removal, and mask revision on a DWI series') + cmdline.add_description('DWI brain masking, response function estimation and bias field correction are inter-related steps and errors ' + 'in each may influence other steps. This script first derives a brain mask (either with a provided algorithm ' + ' or by thresholding a balanced tissue sum image), and then performs response function estimation, multi-tissue ' + 'CSD (with a lower lmax than the dwi2fod default, for speed), and mtnormalise to remove bias field, before the ' + 'DWI brain mask is recalculated. These steps are performed iteratively until either a maximum number of iterations ' + 'or until the brain masks converge.') + cmdline.add_argument('input', help='The input DWI series to be corrected') + cmdline.add_argument('output', help='The output corrected image series') + cmdline.add_argument('-output_mask', help='Output the final brain mask to the specified file') + cmdline.add_argument('-max_iters', type=int, default=DWIBIASCORRECT_MAX_ITERS, + help='The maximum number of iterations. The default is ' + str(DWIBIASCORRECT_MAX_ITERS) + + 'iterations. More iterations may lead to a better mask, but for some problematic data ' + 'this may lead to the masks diverging (in which case a warning is issued).') + cmdline.add_argument('-lmax', type=str, + help='The maximum spherical harmonic order for the output FODs. the value is passed to the dwi2fod command ' + 'and should be provided in the format which it expects (the default value is "' + str(LMAXES_MULTI) + + '" for multi-shell and "' + str(LMAXES_SINGLE) + '" for single-shell data)') + loop_mask_options = cmdline.add_argument_group('Options to specify masking algorithms or tissue sum image threshold in the iterative algorithm') + loop_mask_options.add_argument('-mask_algo', + help='Use a specific dwi2mask algorithm. The default is to use the legacy algorithm or the one specified in the config file.') + loop_mask_options.add_argument('-mask_algo_options', help='Options to pass to the selected dwi2mask masking algorithm') + loop_mask_options.add_argument('-tissue_sum', help='Use a threshold on the balanced tissue sum image instead of a dwi2mask algorithm') + loop_mask_options.add_argument('-threshold', type=float, default=0.5, + help='Threshold on the total tissue density image used to derive the brain mask. the default is 0.5') + init_mask_options = cmdline.add_argument_group('Options to specify initial masking algorithm') + init_mask_options.add_argument('-mask_init', metavar='image', + help='Provide an initial mask to the algorithm and skip the initial masking') + init_mask_options.add_argument('-mask_init_algo', + help='Use a specific dwi2mask algorithm for the initial mask. The default is to use the legacy algorithm or the one specified in the config file.') + init_mask_options.add_argument('-mask_init_options', help='Options to pass to the selected dwi2mask masking algorithm for initial masking') + cmdline.flag_mutually_exclusive_options( [ 'mask_algo', 'tissue_sum' ] ) + cmdline.flag_mutually_exclusive_options( [ 'mask_algo', 'threshold' ] ) + cmdline.flag_mutually_exclusive_options( [ 'mask_algo_options', 'tissue_sum' ] ) + cmdline.flag_mutually_exclusive_options( [ 'mask_algo_options', 'threshold' ] ) + cmdline.flag_mutually_exclusive_options( [ 'mask_init', 'mask_init_algo' ] ) + cmdline.flag_mutually_exclusive_options( [ 'mask_init', 'mask_init_options' ] ) + app.add_dwgrad_import_options(cmdline) + + +def execute(): #pylint: disable=unused-variable + from mrtrix3 import CONFIG, MRtrixError #pylint: disable=no-name-in-module, import-outside-toplevel + from mrtrix3 import app, run, image, matrix, path #pylint: disable=no-name-in-module, import-outside-toplevel + + app.check_output_path(app.ARGS.output) + app.make_scratch_dir() + + grad_import_option = app.read_dwgrad_import_options() + run.command('mrconvert ' + path.from_user(app.ARGS.input) + ' ' + + path.to_scratch('input.mif') + grad_import_option) + + if app.ARGS.mask_init: + run.command('mrconvert ' + path.from_user(app.ARGS.mask_init) + ' ' + + path.to_scratch('dwi_mask_init.mif') + ' -datatype bit') + + app.goto_scratch_dir() + + + # Function to get the appropriate masking command syntax (initial mask or in the loop) + def mask(image, mask, context): + if (context == 'tissue') and app.ARGS.tissue_sum: + TISSUESUM_THRESHOLD = app.ARGS.threshold / math.sqrt(4.0 * math.pi) + syntax = ('mrthreshold ' + + image + + ' -abs ' + + str(TISSUESUM_THRESHOLD) + + ' - |' + + ' maskfilter - connect -largest - |' + + ' mrcalc 1 - -sub - -datatype bit |' + + ' maskfilter - connect -largest - |' + + ' mrcalc 1 - -sub - -datatype bit |' + + ' maskfilter - clean - |' + + ' mrcalc - input_pos_mask.mif -mult ' + + mask + + ' -datatype bit') + elif (context == 'initial') and app.ARGS.mask_init_algo: + syntax = ('dwi2mask ' + + app.ARGS.mask_init_algo + + ((' ' + app.ARGS.mask_init_options) if app.ARGS.mask_init_options else '') + + ' ' + image + ' ' + mask) + else: + syntax = ('dwi2mask ' + + (app.ARGS.mask_algo if app.ARGS.mask_algo else CONFIG.get('Dwi2maskAlgorithm', 'legacy')) + + ((' ' + app.ARGS.mask_algo_options) if app.ARGS.mask_algo_options else '') + + ' ' + image + ' ' + mask) + return syntax + + # Check inputs + # Make sure it's actually a DWI that's been passed + dwi_header = image.Header('input.mif') + if len(dwi_header.size()) != 4: + raise MRtrixError('Input image must be a 4D image') + if 'dw_scheme' not in dwi_header.keyval(): + raise MRtrixError('No valid DW gradient scheme provided or present in image header') + if len(dwi_header.keyval()['dw_scheme']) != dwi_header.size()[3]: + raise MRtrixError('DW gradient scheme contains different number of entries (' + + str(len(dwi_header.keyval()['dw_scheme'])) + + ' to number of volumes in DWIs (' + dwi_header.size()[3] + ')') + + # Determine whether we are working with single-shell or multi-shell data + bvalues = [ + int(round(float(value))) + for value in image.mrinfo('input.mif', 'shell_bvalues') \ + .strip().split()] + multishell = (len(bvalues) > 2) + + + # Generate an initial brain mask if required, or check the initial mask if provided by the user + if app.ARGS.mask_init: + if not image.match('input.mif', 'dwi_mask_init.mif', up_to_dim=3): + raise MRtrixError('Provided mask image does not match input DWI') + else: + app.debug('Performing intial DWI brain masking') + mask_script = mask('input.mif', 'dwi_mask_init.mif', 'initial') + run.command(mask_script) + + dwi_mask_image = 'dwi_mask_init.mif' + + + # Step 2: Combined RF estimation / CSD / mtnormalise / mask revision + class Tissue(object): #pylint: disable=useless-object-inheritance + def __init__(self, name, index): + self.name = name + iter_string = '_iter' + str(index) + self.tissue_rf = 'response_' + name + iter_string + '.txt' + self.fod_init = 'FODinit_' + name + iter_string + '.mif' + self.fod_norm = 'FODnorm_' + name + iter_string + '.mif' + + + app.console('Commencing iterative DWI bias field correction and brain masking ' + 'with a maximum of ' + str(app.ARGS.max_iters) + ' iterations') + + dwi_image = 'input.mif' + iteration = 0 + step = 'initialisation' + dice_coefficient = 0.0 + + def msg(): + return 'Iteration {0}; {1} step; previous Dice coefficient {2}' \ + .format(iteration, step, dice_coefficient) + progress = app.ProgressBar(msg) + + for iteration in range(0, app.ARGS.max_iters): + iter_string = '_iter' + str(iteration+1) + tissues = [Tissue('WM', iteration), + Tissue('GM', iteration), + Tissue('CSF', iteration)] + + step = 'dwi2response' + progress.increment() + run.command('dwi2response dhollander ' + + dwi_image + + ' -mask ' + + dwi_mask_image + + ' ' + + ' '.join(tissue.tissue_rf for tissue in tissues)) + + + # Determine the lmax values to use for the CSD, also remove GM if we can't deal with it + if multishell: + lmaxes = LMAXES_MULTI + else: + lmaxes = LMAXES_SINGLE + app.cleanup(tissues[1].tissue_rf) + tissues = tissues[::2] + if app.ARGS.lmax: + lmaxes = app.ARGS.lmax + if (lmaxes.count(",") + 1) != len(tissues): + raise MRtrixError('Number of lmax values provided is not appropriate for number of shells') + + step = 'dwi2fod' + progress.increment() + app.debug('Performing CSD with lmax values: ' + lmaxes) + run.command('dwi2fod msmt_csd ' + + dwi_image + + ' -lmax ' + lmaxes + + ' ' + + ' '.join(tissue.tissue_rf + ' ' + tissue.fod_init + for tissue in tissues)) + + step = 'mtnormalise' + progress.increment() + field_path = 'field' + iter_string + '.mif' + factors_path = 'factors' + iter_string + '.txt' + + # It might not be nescessary to erode the mask if masking is already restrictive + if 'synthstrip' in mask_script: + b = float(mask_script.strip().split()[mask_script.strip().split().index('-b') + 1]) + if b > 0: + app.debug('Masking is likely restrictive enough, not eroding mask') + run.command('maskfilter ' + dwi_mask_image + ' erode temp_mask.mif') + temp_mask = 'temp_mask.mif' + else: + app.debug('Masking may not be restrictive enough, eroding mask') + temp_mask = dwi_mask_image + + run.command('mtnormalise -mask ' + temp_mask + ' -balanced' + + ' -check_norm ' + field_path + + ' -check_factors ' + factors_path + + ' ' + + ' '.join(tissue.fod_init + ' ' + tissue.fod_norm + for tissue in tissues)) + app.cleanup([tissue.fod_init for tissue in tissues]) + app.cleanup(temp_mask) + + app.console('Iteration ' + str(iteration+1) + ', ' + + 'applying estimated bias field and appropiate scaling factor...') + csf_rf = matrix.load_matrix(tissues[-1].tissue_rf) + csf_rf_bzero_lzero = csf_rf[0][0] + app.cleanup([tissue.tissue_rf for tissue in tissues]) + balance_factors = matrix.load_vector(factors_path) + csf_balance_factor = balance_factors[-1] + app.cleanup(factors_path) + scale_multiplier = (1000.0 * math.sqrt(4.0*math.pi)) / \ + (csf_rf_bzero_lzero / csf_balance_factor) + new_dwi_image = 'dwi' + iter_string + '.mif' + run.command('mrcalc ' + dwi_image + ' ' + + field_path + ' -div ' + + str(scale_multiplier) + ' -mult ' + + new_dwi_image) + app.cleanup(field_path) + + old_dwi_image = dwi_image + dwi_image = new_dwi_image + + tissue_sum_image = 'tissue_sum' + iter_string + '.mif' + + app.console('Iteration ' + str(iteration+1) + ', ' + + 'revising brain mask...') + new_dwi_mask_image = 'dwi_mask' + iter_string + '.mif' + run.command('mrconvert ' + + tissues[0].fod_norm + + ' -coord 3 0 - |' + + ' mrmath - ' + + ' '.join(tissue.fod_norm for tissue in tissues[1:]) + + ' sum ' + + tissue_sum_image) + + mask_script = mask(tissue_sum_image, new_dwi_mask_image, 'tissue') + run.command(mask_script) + app.cleanup([tissue.fod_norm for tissue in tissues]) + + step = 'mrcalc_mask' + progress.increment() + dwi_old_mask_count = image.statistics(dwi_mask_image, + mask=dwi_mask_image).count + dwi_new_mask_count = image.statistics(new_dwi_mask_image, + mask=new_dwi_mask_image).count + app.debug('Old mask: ' + str(dwi_old_mask_count)) + app.debug('New mask: ' + str(dwi_new_mask_count)) + dwi_mask_overlap_image = 'dwi_mask_overlap' + iter_string + '.mif' + run.command('mrcalc ' + + dwi_mask_image + + ' ' + + new_dwi_mask_image + + ' -mult ' + + dwi_mask_overlap_image) + + old_dwi_mask_image = dwi_mask_image + dwi_mask_image = new_dwi_mask_image + + mask_overlap_count = image.statistics(dwi_mask_overlap_image, + mask=dwi_mask_overlap_image).count + app.debug('Mask overlap: ' + str(mask_overlap_count)) + + # Break if either the mask has converged properly, or if it starts to diverge after some convergence + old_dice_coefficient = dice_coefficient + dice_coefficient = 2.0 * mask_overlap_count / \ + (dwi_old_mask_count + dwi_new_mask_count) + + if dice_coefficient > (1.0 - 1e-3): + app.console('Exiting iterative loop due to mask convergence. Iteration: ' + str(iteration + 1)) + app.cleanup(old_dwi_image) + app.cleanup(old_dwi_mask_image) + break + + if (dice_coefficient < old_dice_coefficient) and (old_dice_coefficient > (1.0 - 1e-2)): + app.console('[WARNING] Mask divergence at iteration: ' + str(iteration + 1) + + ' using previous mask with dice coefficient: ' + str(old_dice_coefficient)) + app.cleanup(dwi_mask_image) + app.cleanup(dwi_image) + dwi_image = old_dwi_image + dwi_mask_image = old_dwi_mask_image + break + + if iteration == (DWIBIASCORRECT_MAX_ITERS - 1): + app.console('[WARNING] Maximum iterations without convergence. Final dice coefficient: ' + + str(dice_coefficient)) + app.cleanup(old_dwi_image) + app.cleanup(old_dwi_mask_image) + + + # Processing completed; export + app.console('Processing completed after ' + + str(iteration + 1) + + ' iterations, writing results to output directory') + + run.command('mrconvert ' + + dwi_image + + path.from_user(app.ARGS.output), + force=app.FORCE_OVERWRITE) + + if app.ARGS.output_mask: + run.command('mrconvert ' + + dwi_mask_image + + path.from_user(app.ARGS.output_mask), + force=app.FORCE_OVERWRITE) + + +# Execute the script +import mrtrix3 #pylint: disable=wrong-import-position +mrtrix3.execute() #pylint: disable=no-member \ No newline at end of file diff --git a/lib/mrtrix3/dwi2mask/iterative.py b/lib/mrtrix3/dwi2mask/iterative.py deleted file mode 100644 index 9fe5c49fa8..0000000000 --- a/lib/mrtrix3/dwi2mask/iterative.py +++ /dev/null @@ -1,211 +0,0 @@ -# Copyright (c) 2008-2022 the MRtrix3 contributors. -# -# This Source Code Form is subject to the terms of the Mozilla Public -# License, v. 2.0. If a copy of the MPL was not distributed with this -# file, You can obtain one at http://mozilla.org/MPL/2.0/. -# -# Covered Software is provided under this License on an "as is" -# basis, without warranty of any kind, either expressed, implied, or -# statutory, including, without limitation, warranties that the -# Covered Software is free of defects, merchantable, fit for a -# particular purpose or non-infringing. -# See the Mozilla Public License v. 2.0 for more details. -# -# For more details, see http://www.mrtrix.org/. - -import math -from mrtrix3 import app, run, image, matrix - -DWIBIASCORRECT_MAX_ITERS = 2 - -# We use a threshold on the balanced tissue sum image as a replacement of dwi2mask -# within the iterative bias field correction / brain masking loop -TISSUESUM_THRESHOLD = 0.5 / math.sqrt(4.0 * math.pi) - -LMAXES_MULTI = '4,0,0' -LMAXES_SINGLE = '4,0' - -def usage(base_parser, subparsers): #pylint: disable=unused-variable - parser = subparsers.add_parser('iterative', parents=[base_parser]) - parser.set_author('Robert E. Smith (robert.smith@florey.edu.au) and Arshiya Sangchooli (asangchooli@student.unimelb.edu.au)') - parser.set_synopsis('Iteratively perform masking, RF estimation, CSD, bias field removal, and mask revision to derive a DWI brain mask') - parser.add_description('DWI brain masking may be inaccurate due to residual bias field. This script first derives an initial ' - 'brain mask using the legacy MRtrix3 dwi2mask heuristic (based on thresholded trace images), and then ' - 'performs response function estimation, multi-tissue CSD (with a lower lmax than the dwi2fod default, for ' - 'higher speed), and mtnormalise to remove any bias field, before the DWI brain mask is recalculated. ' - 'These steps are performed iteratively, since each step may influence the others.') - parser.add_argument('input', help='The input DWI series') - parser.add_argument('output', help='The output mask image') - options = parser.add_argument_group('Options specific to the "iterative" algorithm') - options.add_argument('-max_iters', - type=int, - default=DWIBIASCORRECT_MAX_ITERS, - help='the maximum number of iterations. The default is ' + str(DWIBIASCORRECT_MAX_ITERS) + - 'iterations since for some problematic data more iterations may lead to the masks ' - 'diverging, but more iterations can lead to a better mask.') - options.add_argument('-lmax', - type=str, - help='the maximum spherical harmonic order for the output FODs. The value is passed to ' - 'the dwi2fod command and should be provided in the format which it expects ' - '(Default is "' + str(LMAXES_MULTI) + '" for multi-shell and "' + str(LMAXES_SINGLE) + '" for single-shell data)') - - -def get_inputs(): #pylint: disable=unused-variable - pass - - -def needs_mean_bzero(): #pylint: disable=unused-variable - return False - - - - - -def execute(): #pylint: disable=unused-variable - # Determine whether we are working with single-shell or multi-shell data - bvalues = [ - int(round(float(value))) - for value in image.mrinfo('input.mif', 'shell_bvalues') \ - .strip().split()] - multishell = (len(bvalues) > 2) - - # Step 1: Initial DWI brain mask - dwi_mask_image = 'dwi_mask_init.mif' - app.console('Performing intial DWI brain masking with the legacy MRtrix3 algorithm') - run.command('dwi2mask legacy input.mif ' + dwi_mask_image) - - # Step 2: Combined RF estimation / CSD / mtnormalise / mask revision - class Tissue(object): #pylint: disable=useless-object-inheritance - def __init__(self, name, index): - self.name = name - iter_string = '_iter' + str(index) - self.tissue_rf = 'response_' + name + iter_string + '.txt' - self.fod_init = 'FODinit_' + name + iter_string + '.mif' - self.fod_norm = 'FODnorm_' + name + iter_string + '.mif' - - app.console('Commencing iterative DWI bias field correction and brain masking ' - 'with a maximum of ' + str(app.ARGS.max_iters) + ' iterations') - - dwi_image = 'input.mif' - - for iteration in range(0, app.ARGS.max_iters): - iter_string = '_iter' + str(iteration+1) - tissues = [Tissue('WM', iteration), - Tissue('GM', iteration), - Tissue('CSF', iteration)] - app.console('Iteration ' + str(iteration) + ', ' - + 'estimating response function using brain mask...') - run.command('dwi2response dhollander ' - + dwi_image - + ' -mask ' - + dwi_mask_image - + ' ' - + ' '.join(tissue.tissue_rf for tissue in tissues)) - - # Remove GM if we can't deal with it - if app.ARGS.lmax: - lmaxes = app.ARGS.lmax - else: - lmaxes = LMAXES_MULTI - if not multishell: - app.cleanup(tissues[1].tissue_rf) - tissues = tissues[::2] - lmaxes = LMAXES_SINGLE - - app.console('Iteration ' + str(iteration) + ', ' - + 'FOD estimation with lmaxes ' + lmaxes + '...') - run.command('dwi2fod msmt_csd ' - + dwi_image - + ' -lmax ' + lmaxes - + ' ' - + ' '.join(tissue.tissue_rf + ' ' + tissue.fod_init - for tissue in tissues)) - - app.console('Iteration ' + str(iteration) + ', ' - + 'running mtnormalise...') - field_path = 'field' + iter_string + '.mif' - factors_path = 'factors' + iter_string + '.txt' - run.command('maskfilter ' + dwi_mask_image + ' erode - |' - + ' mtnormalise -mask - -balanced' - + ' -check_norm ' + field_path - + ' -check_factors ' + factors_path - + ' ' - + ' '.join(tissue.fod_init + ' ' + tissue.fod_norm - for tissue in tissues)) - app.cleanup([tissue.fod_init for tissue in tissues]) - - app.console('Iteration ' + str(iteration) + ', ' - + 'applying estimated bias field and appropiate scaling factor...') - csf_rf = matrix.load_matrix(tissues[-1].tissue_rf) - csf_rf_bzero_lzero = csf_rf[0][0] - app.cleanup([tissue.tissue_rf for tissue in tissues]) - balance_factors = matrix.load_vector(factors_path) - csf_balance_factor = balance_factors[-1] - app.cleanup(factors_path) - scale_multiplier = (1000.0 * math.sqrt(4.0*math.pi)) / \ - (csf_rf_bzero_lzero / csf_balance_factor) - new_dwi_image = 'dwi' + iter_string + '.mif' - run.command('mrcalc ' + dwi_image + ' ' - + field_path + ' -div ' - + str(scale_multiplier) + ' -mult ' - + new_dwi_image) - app.cleanup(field_path) - app.cleanup(dwi_image) - dwi_image = new_dwi_image - - app.console('Iteration ' + str(iteration) + ', ' - + 'revising brain mask...') - new_dwi_mask_image = 'dwi_mask' + iter_string + '.mif' - run.command('mrconvert ' - + tissues[0].fod_norm - + ' -coord 3 0 - |' - + ' mrmath - ' - + ' '.join(tissue.fod_norm for tissue in tissues[1:]) - + ' sum - |' - + ' mrthreshold - -abs ' - + str(TISSUESUM_THRESHOLD) - + ' - |' - + ' maskfilter - connect -largest - |' - + ' mrcalc 1 - -sub - -datatype bit |' - + ' maskfilter - connect -largest - |' - + ' mrcalc 1 - -sub - -datatype bit |' - + ' maskfilter - clean - |' - + ' mrcalc - input_pos_mask.mif -mult ' - + new_dwi_mask_image - + ' -datatype bit') - app.cleanup([tissue.fod_norm for tissue in tissues]) - - app.console('Iteration ' + str(iteration) + ', ' - + 'comparing the new mask with the previous iteration...') - dwi_old_mask_count = image.statistics(dwi_mask_image, - mask=dwi_mask_image).count - dwi_new_mask_count = image.statistics(new_dwi_mask_image, - mask=new_dwi_mask_image).count - app.debug('Old mask: ' + str(dwi_old_mask_count)) - app.debug('New mask: ' + str(dwi_new_mask_count)) - dwi_mask_overlap_image = 'dwi_mask_overlap' + iter_string + '.mif' - run.command('mrcalc ' - + dwi_mask_image - + ' ' - + new_dwi_mask_image - + ' -mult ' - + dwi_mask_overlap_image) - app.cleanup(dwi_mask_image) - dwi_mask_image = new_dwi_mask_image - mask_overlap_count = image.statistics(dwi_mask_overlap_image, - mask=dwi_mask_overlap_image).count - app.debug('Mask overlap: ' + str(mask_overlap_count)) - dice_coefficient = 2.0 * mask_overlap_count / \ - (dwi_old_mask_count + dwi_new_mask_count) - app.console('Iteration ' + str(iteration) - + ' dice coefficient: ' + str(dice_coefficient) + '...') - if dice_coefficient > (1.0 - 1e-3): - app.console('Exiting iterative loop due to mask convergence') - break - - # Processing completed; export - app.console('Processing completed after ' - + str(iteration + 1) - + ' iterations, writing results to output directory') - - return dwi_mask_image From 78a3154fe55aeefda43e9cd4ba32be6d951b1a70 Mon Sep 17 00:00:00 2001 From: Arshiyas Date: Mon, 20 Feb 2023 18:23:45 +1100 Subject: [PATCH 04/24] dwibiasnormmask: Integrate SynthStrip --- bin/dwibiasnormmask | 140 ++++++++++++++++++++------------------------ 1 file changed, 65 insertions(+), 75 deletions(-) mode change 100644 => 100755 bin/dwibiasnormmask diff --git a/bin/dwibiasnormmask b/bin/dwibiasnormmask old mode 100644 new mode 100755 index a0f6cdcca4..ae4f436674 --- a/bin/dwibiasnormmask +++ b/bin/dwibiasnormmask @@ -15,7 +15,7 @@ # # For more details, see http://www.mrtrix.org/. -import math +import math, shutil DWIBIASCORRECT_MAX_ITERS = 2 LMAXES_MULTI = '4,0,0' @@ -31,7 +31,7 @@ def usage(cmdline): #pylint: disable=unused-variable ' or by thresholding a balanced tissue sum image), and then performs response function estimation, multi-tissue ' 'CSD (with a lower lmax than the dwi2fod default, for speed), and mtnormalise to remove bias field, before the ' 'DWI brain mask is recalculated. These steps are performed iteratively until either a maximum number of iterations ' - 'or until the brain masks converge.') + 'or until the brain masks converge. If SynthStrip is installed, it will be used to derive brain masks.') cmdline.add_argument('input', help='The input DWI series to be corrected') cmdline.add_argument('output', help='The output corrected image series') cmdline.add_argument('-output_mask', help='Output the final brain mask to the specified file') @@ -43,25 +43,11 @@ def usage(cmdline): #pylint: disable=unused-variable help='The maximum spherical harmonic order for the output FODs. the value is passed to the dwi2fod command ' 'and should be provided in the format which it expects (the default value is "' + str(LMAXES_MULTI) + '" for multi-shell and "' + str(LMAXES_SINGLE) + '" for single-shell data)') - loop_mask_options = cmdline.add_argument_group('Options to specify masking algorithms or tissue sum image threshold in the iterative algorithm') - loop_mask_options.add_argument('-mask_algo', - help='Use a specific dwi2mask algorithm. The default is to use the legacy algorithm or the one specified in the config file.') - loop_mask_options.add_argument('-mask_algo_options', help='Options to pass to the selected dwi2mask masking algorithm') - loop_mask_options.add_argument('-tissue_sum', help='Use a threshold on the balanced tissue sum image instead of a dwi2mask algorithm') - loop_mask_options.add_argument('-threshold', type=float, default=0.5, + mask_options = cmdline.add_argument_group('Options to specify an initial mask or tissue sum image threshold in the iterative algorithm') + mask_options.add_argument('-threshold', type=float, default=0.5, help='Threshold on the total tissue density image used to derive the brain mask. the default is 0.5') - init_mask_options = cmdline.add_argument_group('Options to specify initial masking algorithm') - init_mask_options.add_argument('-mask_init', metavar='image', - help='Provide an initial mask to the algorithm and skip the initial masking') - init_mask_options.add_argument('-mask_init_algo', - help='Use a specific dwi2mask algorithm for the initial mask. The default is to use the legacy algorithm or the one specified in the config file.') - init_mask_options.add_argument('-mask_init_options', help='Options to pass to the selected dwi2mask masking algorithm for initial masking') - cmdline.flag_mutually_exclusive_options( [ 'mask_algo', 'tissue_sum' ] ) - cmdline.flag_mutually_exclusive_options( [ 'mask_algo', 'threshold' ] ) - cmdline.flag_mutually_exclusive_options( [ 'mask_algo_options', 'tissue_sum' ] ) - cmdline.flag_mutually_exclusive_options( [ 'mask_algo_options', 'threshold' ] ) - cmdline.flag_mutually_exclusive_options( [ 'mask_init', 'mask_init_algo' ] ) - cmdline.flag_mutually_exclusive_options( [ 'mask_init', 'mask_init_options' ] ) + mask_options.add_argument('-mask_init', metavar='image', + help='Provide an initial mask to the algorithm and skip the initial masking') app.add_dwgrad_import_options(cmdline) @@ -83,35 +69,6 @@ def execute(): #pylint: disable=unused-variable app.goto_scratch_dir() - # Function to get the appropriate masking command syntax (initial mask or in the loop) - def mask(image, mask, context): - if (context == 'tissue') and app.ARGS.tissue_sum: - TISSUESUM_THRESHOLD = app.ARGS.threshold / math.sqrt(4.0 * math.pi) - syntax = ('mrthreshold ' - + image - + ' -abs ' - + str(TISSUESUM_THRESHOLD) - + ' - |' - + ' maskfilter - connect -largest - |' - + ' mrcalc 1 - -sub - -datatype bit |' - + ' maskfilter - connect -largest - |' - + ' mrcalc 1 - -sub - -datatype bit |' - + ' maskfilter - clean - |' - + ' mrcalc - input_pos_mask.mif -mult ' - + mask - + ' -datatype bit') - elif (context == 'initial') and app.ARGS.mask_init_algo: - syntax = ('dwi2mask ' - + app.ARGS.mask_init_algo - + ((' ' + app.ARGS.mask_init_options) if app.ARGS.mask_init_options else '') - + ' ' + image + ' ' + mask) - else: - syntax = ('dwi2mask ' - + (app.ARGS.mask_algo if app.ARGS.mask_algo else CONFIG.get('Dwi2maskAlgorithm', 'legacy')) - + ((' ' + app.ARGS.mask_algo_options) if app.ARGS.mask_algo_options else '') - + ' ' + image + ' ' + mask) - return syntax - # Check inputs # Make sure it's actually a DWI that's been passed dwi_header = image.Header('input.mif') @@ -131,6 +88,10 @@ def execute(): #pylint: disable=unused-variable .strip().split()] multishell = (len(bvalues) > 2) + # Get the SynthStrip path if it's installed + synthstrip_cmd = shutil.which('mri_synthstrip') + if not synthstrip_cmd: + synthstrip_cmd = shutil.which('synthstrip-singularity') # Generate an initial brain mask if required, or check the initial mask if provided by the user if app.ARGS.mask_init: @@ -138,13 +99,22 @@ def execute(): #pylint: disable=unused-variable raise MRtrixError('Provided mask image does not match input DWI') else: app.debug('Performing intial DWI brain masking') - mask_script = mask('input.mif', 'dwi_mask_init.mif', 'initial') - run.command(mask_script) + if synthstrip_cmd: + run.command('dwiextract input.mif -bzero - | ' + 'mrmath - mean - -axis 3 | ' + 'mrconvert - bzero.nii -strides +1,+2,+3') + run.command(synthstrip_cmd + ' -i bzero.nii -b 0 --mask dwi_mask_init.nii') + run.command('mrconvert dwi_mask_init.nii dwi_mask_init.mif -datatype bit') + app.cleanup('bzero.nii') + app.cleanup('dwi_mask_init.nii') + else: + run.command('dwi2mask ' + + CONFIG.get('Dwi2maskAlgorithm', 'legacy') + + ' input.mif dwi_mask_init.mif') dwi_mask_image = 'dwi_mask_init.mif' - - # Step 2: Combined RF estimation / CSD / mtnormalise / mask revision + # Combined RF estimation / CSD / mtnormalise / mask revision class Tissue(object): #pylint: disable=useless-object-inheritance def __init__(self, name, index): self.name = name @@ -154,7 +124,7 @@ def execute(): #pylint: disable=unused-variable self.fod_norm = 'FODnorm_' + name + iter_string + '.mif' - app.console('Commencing iterative DWI bias field correction and brain masking ' + app.debug('Commencing iterative DWI bias field correction and brain masking ' 'with a maximum of ' + str(app.ARGS.max_iters) + ' iterations') dwi_image = 'input.mif' @@ -211,26 +181,23 @@ def execute(): #pylint: disable=unused-variable factors_path = 'factors' + iter_string + '.txt' # It might not be nescessary to erode the mask if masking is already restrictive - if 'synthstrip' in mask_script: - b = float(mask_script.strip().split()[mask_script.strip().split().index('-b') + 1]) - if b > 0: - app.debug('Masking is likely restrictive enough, not eroding mask') - run.command('maskfilter ' + dwi_mask_image + ' erode temp_mask.mif') - temp_mask = 'temp_mask.mif' - else: - app.debug('Masking may not be restrictive enough, eroding mask') - temp_mask = dwi_mask_image - - run.command('mtnormalise -mask ' + temp_mask + ' -balanced' + if synthstrip_cmd: + app.debug('Masking is likely restrictive enough, not eroding mask') + run.command('cp ' + dwi_mask_image + ' temp_mask.mif') + else: + app.debug('Masking may not be restrictive enough, eroding mask') + run.command('maskfilter ' + dwi_mask_image + ' erode temp_mask.mif') + + run.command('mtnormalise -mask temp_mask.mif -balanced' + ' -check_norm ' + field_path + ' -check_factors ' + factors_path + ' ' + ' '.join(tissue.fod_init + ' ' + tissue.fod_norm for tissue in tissues)) app.cleanup([tissue.fod_init for tissue in tissues]) - app.cleanup(temp_mask) + app.cleanup('temp_mask.mif') - app.console('Iteration ' + str(iteration+1) + ', ' + app.debug('Iteration ' + str(iteration+1) + ', ' + 'applying estimated bias field and appropiate scaling factor...') csf_rf = matrix.load_matrix(tissues[-1].tissue_rf) csf_rf_bzero_lzero = csf_rf[0][0] @@ -252,7 +219,7 @@ def execute(): #pylint: disable=unused-variable tissue_sum_image = 'tissue_sum' + iter_string + '.mif' - app.console('Iteration ' + str(iteration+1) + ', ' + app.debug('Iteration ' + str(iteration+1) + ', ' + 'revising brain mask...') new_dwi_mask_image = 'dwi_mask' + iter_string + '.mif' run.command('mrconvert ' @@ -263,10 +230,33 @@ def execute(): #pylint: disable=unused-variable + ' sum ' + tissue_sum_image) - mask_script = mask(tissue_sum_image, new_dwi_mask_image, 'tissue') - run.command(mask_script) app.cleanup([tissue.fod_norm for tissue in tissues]) + if synthstrip_cmd: + run.command('mrconvert ' + tissue_sum_image + ' tissue.nii -force') + run.command(synthstrip_cmd + ' -i tissue.nii -b 0 --mask dwi_mask_temp.nii') + run.command('mrconvert dwi_mask_temp.nii ' + new_dwi_mask_image + ' -datatype bit') + app.cleanup('tissue.nii') + app.cleanup('dwi_mask_temp.nii') + else: + TISSUESUM_THRESHOLD = app.ARGS.threshold / math.sqrt(4.0 * math.pi) + run.command('mrthreshold ' + + tissue_sum_image + + ' -abs ' + + str(TISSUESUM_THRESHOLD) + + ' - |' + + ' maskfilter - connect -largest - |' + + ' mrcalc 1 - -sub - -datatype bit |' + + ' maskfilter - connect -largest - |' + + ' mrcalc 1 - -sub - -datatype bit |' + + ' maskfilter - clean - |' + + ' mrcalc - input_pos_mask.mif -mult ' + + new_dwi_mask_image + + ' -datatype bit') + + + app.cleanup(tissue_sum_image) + step = 'mrcalc_mask' progress.increment() dwi_old_mask_count = image.statistics(dwi_mask_image, @@ -296,13 +286,13 @@ def execute(): #pylint: disable=unused-variable (dwi_old_mask_count + dwi_new_mask_count) if dice_coefficient > (1.0 - 1e-3): - app.console('Exiting iterative loop due to mask convergence. Iteration: ' + str(iteration + 1)) + app.debug('Exiting iterative loop due to mask convergence. Iteration: ' + str(iteration + 1)) app.cleanup(old_dwi_image) app.cleanup(old_dwi_mask_image) break if (dice_coefficient < old_dice_coefficient) and (old_dice_coefficient > (1.0 - 1e-2)): - app.console('[WARNING] Mask divergence at iteration: ' + str(iteration + 1) + app.debug('[WARNING] Mask divergence at iteration: ' + str(iteration + 1) + ' using previous mask with dice coefficient: ' + str(old_dice_coefficient)) app.cleanup(dwi_mask_image) app.cleanup(dwi_image) @@ -311,7 +301,7 @@ def execute(): #pylint: disable=unused-variable break if iteration == (DWIBIASCORRECT_MAX_ITERS - 1): - app.console('[WARNING] Maximum iterations without convergence. Final dice coefficient: ' + app.debug('[WARNING] Maximum iterations without convergence. Final dice coefficient: ' + str(dice_coefficient)) app.cleanup(old_dwi_image) app.cleanup(old_dwi_mask_image) @@ -323,13 +313,13 @@ def execute(): #pylint: disable=unused-variable + ' iterations, writing results to output directory') run.command('mrconvert ' - + dwi_image + + dwi_image + ' ' + path.from_user(app.ARGS.output), force=app.FORCE_OVERWRITE) if app.ARGS.output_mask: run.command('mrconvert ' - + dwi_mask_image + + dwi_mask_image + ' ' + path.from_user(app.ARGS.output_mask), force=app.FORCE_OVERWRITE) From 3fcd8cbca96e67a309371094f3dd5f525d34b3cc Mon Sep 17 00:00:00 2001 From: Arshiyas Date: Thu, 16 Mar 2023 10:38:45 +1100 Subject: [PATCH 05/24] Docs: Generate for dwibiasnormmask and mtnorm algorithms --- bin/dwibiasnormmask | 24 ++--- docs/reference/commands/dwi2mask.rst | 108 ++++++++++---------- docs/reference/commands/dwibiascorrect.rst | 103 ++++++++++++++++++- docs/reference/commands/dwibiasnormmask.rst | 102 ++++++++++++++++++ docs/reference/commands/dwinormalise.rst | 101 +++++++++++++++++- docs/reference/commands_list.rst | 2 + lib/mrtrix3/dwi2mask/mtnorm.py | 7 +- lib/mrtrix3/dwibiascorrect/mtnorm.py | 6 +- lib/mrtrix3/dwinormalise/mtnorm.py | 12 +-- 9 files changed, 380 insertions(+), 85 deletions(-) create mode 100644 docs/reference/commands/dwibiasnormmask.rst diff --git a/bin/dwibiasnormmask b/bin/dwibiasnormmask index ae4f436674..4efbebdc7c 100755 --- a/bin/dwibiasnormmask +++ b/bin/dwibiasnormmask @@ -41,7 +41,7 @@ def usage(cmdline): #pylint: disable=unused-variable 'this may lead to the masks diverging (in which case a warning is issued).') cmdline.add_argument('-lmax', type=str, help='The maximum spherical harmonic order for the output FODs. the value is passed to the dwi2fod command ' - 'and should be provided in the format which it expects (the default value is "' + str(LMAXES_MULTI) + + 'and should be provided in the format which it expects (the default value is "' + str(LMAXES_MULTI) + '" for multi-shell and "' + str(LMAXES_SINGLE) + '" for single-shell data)') mask_options = cmdline.add_argument_group('Options to specify an initial mask or tissue sum image threshold in the iterative algorithm') mask_options.add_argument('-threshold', type=float, default=0.5, @@ -54,7 +54,7 @@ def usage(cmdline): #pylint: disable=unused-variable def execute(): #pylint: disable=unused-variable from mrtrix3 import CONFIG, MRtrixError #pylint: disable=no-name-in-module, import-outside-toplevel from mrtrix3 import app, run, image, matrix, path #pylint: disable=no-name-in-module, import-outside-toplevel - + app.check_output_path(app.ARGS.output) app.make_scratch_dir() @@ -133,8 +133,8 @@ def execute(): #pylint: disable=unused-variable dice_coefficient = 0.0 def msg(): - return 'Iteration {0}; {1} step; previous Dice coefficient {2}' \ - .format(iteration, step, dice_coefficient) + return 'Iteration {0}; {1} step; previous Dice coefficient {2}' \ + .format(iteration, step, dice_coefficient) progress = app.ProgressBar(msg) for iteration in range(0, app.ARGS.max_iters): @@ -142,7 +142,7 @@ def execute(): #pylint: disable=unused-variable tissues = [Tissue('WM', iteration), Tissue('GM', iteration), Tissue('CSF', iteration)] - + step = 'dwi2response' progress.increment() run.command('dwi2response dhollander ' @@ -174,12 +174,12 @@ def execute(): #pylint: disable=unused-variable + ' ' + ' '.join(tissue.tissue_rf + ' ' + tissue.fod_init for tissue in tissues)) - + step = 'mtnormalise' progress.increment() field_path = 'field' + iter_string + '.mif' factors_path = 'factors' + iter_string + '.txt' - + # It might not be nescessary to erode the mask if masking is already restrictive if synthstrip_cmd: app.debug('Masking is likely restrictive enough, not eroding mask') @@ -239,11 +239,11 @@ def execute(): #pylint: disable=unused-variable app.cleanup('tissue.nii') app.cleanup('dwi_mask_temp.nii') else: - TISSUESUM_THRESHOLD = app.ARGS.threshold / math.sqrt(4.0 * math.pi) - run.command('mrthreshold ' + tissuesum_threshold = app.ARGS.threshold / math.sqrt(4.0 * math.pi) + run.command('mrthreshold ' + tissue_sum_image + ' -abs ' - + str(TISSUESUM_THRESHOLD) + + str(tissuesum_threshold) + ' - |' + ' maskfilter - connect -largest - |' + ' mrcalc 1 - -sub - -datatype bit |' @@ -279,7 +279,7 @@ def execute(): #pylint: disable=unused-variable mask_overlap_count = image.statistics(dwi_mask_overlap_image, mask=dwi_mask_overlap_image).count app.debug('Mask overlap: ' + str(mask_overlap_count)) - + # Break if either the mask has converged properly, or if it starts to diverge after some convergence old_dice_coefficient = dice_coefficient dice_coefficient = 2.0 * mask_overlap_count / \ @@ -326,4 +326,4 @@ def execute(): #pylint: disable=unused-variable # Execute the script import mrtrix3 #pylint: disable=wrong-import-position -mrtrix3.execute() #pylint: disable=no-member \ No newline at end of file +mrtrix3.execute() #pylint: disable=no-member diff --git a/docs/reference/commands/dwi2mask.rst b/docs/reference/commands/dwi2mask.rst index 34d0558eb3..05f2c774d8 100644 --- a/docs/reference/commands/dwi2mask.rst +++ b/docs/reference/commands/dwi2mask.rst @@ -15,7 +15,7 @@ Usage dwi2mask algorithm [ options ] ... -- *algorithm*: Select the algorithm to be used to complete the script operation; additional details and options become available once an algorithm is nominated. Options are: 3dautomask, ants, b02template, consensus, fslbet, hdbet, iterative, legacy, mean, synthstrip, trace +- *algorithm*: Select the algorithm to be used to complete the script operation; additional details and options become available once an algorithm is nominated. Options are: 3dautomask, ants, b02template, consensus, fslbet, hdbet, legacy, mean, mtnorm, synthstrip, trace Description ----------- @@ -685,40 +685,30 @@ See the Mozilla Public License v. 2.0 for more details. For more details, see http://www.mrtrix.org/. -.. _dwi2mask_iterative: +.. _dwi2mask_legacy: -dwi2mask iterative -================== +dwi2mask legacy +=============== Synopsis -------- -Iteratively perform masking, RF estimation, CSD, bias field removal, and mask revision to derive a DWI brain mask +Use the legacy MRtrix3 dwi2mask heuristic (based on thresholded trace images) Usage ----- :: - dwi2mask iterative input output [ options ] + dwi2mask legacy input output [ options ] - *input*: The input DWI series - *output*: The output mask image -Description ------------ - -DWI brain masking may be inaccurate due to residual bias field. This script first derives an initial brain mask using the legacy MRtrix3 dwi2mask heuristic (based on thresholded trace images), and then performs response function estimation, multi-tissue CSD (with a lower lmax than the dwi2fod default, for higher speed), and mtnormalise to remove any bias field, before the DWI brain mask is recalculated. These steps are performed iteratively, since each step may influence the others. - Options ------- -Options specific to the "iterative" algorithm -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -- **-max_iters** the maximum number of iterations. The default is 2iterations since for some problematic data more iterations may lead to the masks diverging, but more iterations can lead to a better mask. - -- **-lmax** the maximum spherical harmonic order for the output FODs. The value is passed to the dwi2fod command and should be provided in the format which it expects (Default is "4,0,0" for multi-shell and "4,0" for single-shell data) +- **-clean_scale** the maximum scale used to cut bridges. A certain maximum scale cuts bridges up to a width (in voxels) of 2x the provided scale. Setting this to 0 disables the mask cleaning step. (Default: 2) Options for importing the diffusion gradient table ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -764,9 +754,9 @@ Tournier, J.-D.; Smith, R. E.; Raffelt, D.; Tabbara, R.; Dhollander, T.; Pietsch -**Author:** Robert E. Smith (robert.smith@florey.edu.au) and Arshiya Sangchooli (asangchooli@student.unimelb.edu.au) +**Author:** Robert E. Smith (robert.smith@florey.edu.au) -**Copyright:** Copyright (c) 2008-2022 the MRtrix3 contributors. +**Copyright:** Copyright (c) 2008-2023 the MRtrix3 contributors. This Source Code Form is subject to the terms of the Mozilla Public License, v. 2.0. If a copy of the MPL was not distributed with this @@ -781,22 +771,22 @@ See the Mozilla Public License v. 2.0 for more details. For more details, see http://www.mrtrix.org/. -.. _dwi2mask_legacy: +.. _dwi2mask_mean: -dwi2mask legacy -=============== +dwi2mask mean +============= Synopsis -------- -Use the legacy MRtrix3 dwi2mask heuristic (based on thresholded trace images) +Generate a mask based on simply averaging all volumes in the DWI series Usage ----- :: - dwi2mask legacy input output [ options ] + dwi2mask mean input output [ options ] - *input*: The input DWI series - *output*: The output mask image @@ -804,6 +794,11 @@ Usage Options ------- +Options specific to the 'mean' algorithm +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-shells** Comma separated list of shells to be included in the volume averaging + - **-clean_scale** the maximum scale used to cut bridges. A certain maximum scale cuts bridges up to a width (in voxels) of 2x the provided scale. Setting this to 0 disables the mask cleaning step. (Default: 2) Options for importing the diffusion gradient table @@ -850,7 +845,7 @@ Tournier, J.-D.; Smith, R. E.; Raffelt, D.; Tabbara, R.; Dhollander, T.; Pietsch -**Author:** Robert E. Smith (robert.smith@florey.edu.au) +**Author:** Warda Syeda (wtsyeda@unimelb.edu.au) **Copyright:** Copyright (c) 2008-2023 the MRtrix3 contributors. @@ -867,35 +862,46 @@ See the Mozilla Public License v. 2.0 for more details. For more details, see http://www.mrtrix.org/. -.. _dwi2mask_mean: +.. _dwi2mask_synthstrip: -dwi2mask mean -============= +dwi2mask synthstrip +=================== Synopsis -------- -Generate a mask based on simply averaging all volumes in the DWI series +Use the FreeSurfer Synthstrip method on the mean b=0 image Usage ----- :: - dwi2mask mean input output [ options ] + dwi2mask synthstrip input output [ options ] - *input*: The input DWI series - *output*: The output mask image +Description +----------- + +This algorithm requires that the SynthStrip method be installed and available via PATH; this could be via Freesufer 7.3.0 or later, or the dedicated Singularity container. + Options ------- -Options specific to the 'mean' algorithm -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +Options specific to the 'Synthstrip' algorithm +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -- **-shells** Comma separated list of shells to be included in the volume averaging +- **-stripped** The output stripped image -- **-clean_scale** the maximum scale used to cut bridges. A certain maximum scale cuts bridges up to a width (in voxels) of 2x the provided scale. Setting this to 0 disables the mask cleaning step. (Default: 2) +- **-gpu** Use the GPU + +- **-model file** Alternative model weights + +- **-nocsf** Compute the immediate boundary of brain matter excluding surrounding CSF + +- **-border** Control the boundary distance from the brain Options for importing the diffusion gradient table ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -935,13 +941,15 @@ Standard options References ^^^^^^^^^^ +* A. Hoopes, J. S. Mora, A. V. Dalca, B. Fischl, M. Hoffmann. SynthStrip: Skull-Stripping for Any Brain Image. NeuroImage, 2022, 260, 119474 + Tournier, J.-D.; Smith, R. E.; Raffelt, D.; Tabbara, R.; Dhollander, T.; Pietsch, M.; Christiaens, D.; Jeurissen, B.; Yeh, C.-H. & Connelly, A. MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation. NeuroImage, 2019, 202, 116137 -------------- -**Author:** Warda Syeda (wtsyeda@unimelb.edu.au) +**Author:** Ruobing Chen (chrc@student.unimelb.edu.au) **Copyright:** Copyright (c) 2008-2023 the MRtrix3 contributors. @@ -958,22 +966,22 @@ See the Mozilla Public License v. 2.0 for more details. For more details, see http://www.mrtrix.org/. -.. _dwi2mask_synthstrip: +.. _dwi2mask_mtnorm: -dwi2mask synthstrip -=================== +dwi2mask mtnorm +=============== Synopsis -------- -Use the FreeSurfer Synthstrip method on the mean b=0 image +Derives a DWI brain mask by calculating and then thresholding a tissue density image Usage ----- :: - dwi2mask synthstrip input output [ options ] + dwi2mask mtnorm input output [ options ] - *input*: The input DWI series - *output*: The output mask image @@ -981,23 +989,17 @@ Usage Description ----------- -This algorithm requires that the SynthStrip method be installed and available via PATH; this could be via Freesufer 7.3.0 or later, or the dedicated Singularity container. +This script first derives an initial brain mask using the legacy MRtrix3 dwi2mask heuristic (based on thresholded trace images), and then performs response function estimation and multi-tissue CSD (with a lower lmax than the dwi2fod default, for higher speed), before summing the derived tissue ODFs and thresholding the resulting tissue density image to derive a DWI brain mask. Options ------- -Options specific to the 'Synthstrip' algorithm -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -- **-stripped** The output stripped image - -- **-gpu** Use the GPU - -- **-model file** Alternative model weights +Options specific to the "mtnorm" algorithm +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -- **-nocsf** Compute the immediate boundary of brain matter excluding surrounding CSF +- **-threshold** the threshold on the total tissue density sum image used to derive the brain mask. the default is 0.5 -- **-border** Control the boundary distance from the brain +- **-lmax** the maximum spherical harmonic order for the output FODs. The value is passed to the dwi2fod command and should be provided in the format which it expects (Default is "4,0,0" for multi-shell and "4,0" for single-shell data) Options for importing the diffusion gradient table ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -1037,15 +1039,13 @@ Standard options References ^^^^^^^^^^ -* A. Hoopes, J. S. Mora, A. V. Dalca, B. Fischl, M. Hoffmann. SynthStrip: Skull-Stripping for Any Brain Image. NeuroImage, 2022, 260, 119474 - Tournier, J.-D.; Smith, R. E.; Raffelt, D.; Tabbara, R.; Dhollander, T.; Pietsch, M.; Christiaens, D.; Jeurissen, B.; Yeh, C.-H. & Connelly, A. MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation. NeuroImage, 2019, 202, 116137 -------------- -**Author:** Ruobing Chen (chrc@student.unimelb.edu.au) +**Author:** Robert E. Smith (robert.smith@florey.edu.au) and Arshiya Sangchooli (asangchooli@student.unimelb.edu.au) **Copyright:** Copyright (c) 2008-2023 the MRtrix3 contributors. diff --git a/docs/reference/commands/dwibiascorrect.rst b/docs/reference/commands/dwibiascorrect.rst index 1feee0e948..7d8444ed53 100644 --- a/docs/reference/commands/dwibiascorrect.rst +++ b/docs/reference/commands/dwibiascorrect.rst @@ -15,7 +15,7 @@ Usage dwibiascorrect algorithm [ options ] ... -- *algorithm*: Select the algorithm to be used to complete the script operation; additional details and options become available once an algorithm is nominated. Options are: ants, fsl +- *algorithm*: Select the algorithm to be used to complete the script operation; additional details and options become available once an algorithm is nominated. Options are: ants, fsl, mtnorm Description ----------- @@ -296,3 +296,104 @@ See the Mozilla Public License v. 2.0 for more details. For more details, see http://www.mrtrix.org/. +.. _dwibiascorrect_mtnorm: + +dwibiascorrect mtnorm +===================== + +Synopsis +-------- + +Performs DWI bias field correction + +Usage +----- + +:: + + dwibiascorrect mtnorm input output [ options ] + +- *input*: The input image series to be corrected +- *output*: The output corrected image series + +Description +----------- + +This script first performs response function estimation and multi-tissue CSD within a provided or derived brain mask (with a lower lmax than the dwi2fod default, for higher speed) before calling mtnormalise and utilizing the final estimated spatially varying intensity level used for normalisation to correct bias fields + +Options +------- + +Options specific to the "mtnorm" algorithm +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-lmax** the maximum spherical harmonic order for the output FODs. The value is passed to the dwi2fod command and should be provided in the format which it expects (Default is "4,0,0" for multi-shell and "4,0" for single-shell data) + +Options for importing the diffusion gradient table +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-grad** Provide the diffusion gradient table in MRtrix format + +- **-fslgrad bvecs bvals** Provide the diffusion gradient table in FSL bvecs/bvals format + +Options common to all dwibiascorrect algorithms +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-mask image** Manually provide a mask image for bias field estimation + +- **-bias image** Output the estimated bias field + +Additional standard options for Python scripts +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-nocleanup** do not delete intermediate files during script execution, and do not delete scratch directory at script completion. + +- **-scratch /path/to/scratch/** manually specify the path in which to generate the scratch directory. + +- **-continue ** continue the script from a previous execution; must provide the scratch directory path, and the name of the last successfully-generated file. + +Standard options +^^^^^^^^^^^^^^^^ + +- **-info** display information messages. + +- **-quiet** do not display information messages or progress status. Alternatively, this can be achieved by setting the MRTRIX_QUIET environment variable to a non-empty string. + +- **-debug** display debugging messages. + +- **-force** force overwrite of output files. + +- **-nthreads number** use this number of threads in multi-threaded applications (set to 0 to disable multi-threading). + +- **-config key value** *(multiple uses permitted)* temporarily set the value of an MRtrix config file entry. + +- **-help** display this information page and exit. + +- **-version** display version information and exit. + +References +^^^^^^^^^^ + +Tournier, J.-D.; Smith, R. E.; Raffelt, D.; Tabbara, R.; Dhollander, T.; Pietsch, M.; Christiaens, D.; Jeurissen, B.; Yeh, C.-H. & Connelly, A. MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation. NeuroImage, 2019, 202, 116137 + +-------------- + + + +**Author:** Robert E. Smith (robert.smith@florey.edu.au) and Arshiya Sangchooli (asangchooli@student.unimelb.edu.au) + +**Copyright:** Copyright (c) 2008-2023 the MRtrix3 contributors. + +This Source Code Form is subject to the terms of the Mozilla Public +License, v. 2.0. If a copy of the MPL was not distributed with this +file, You can obtain one at http://mozilla.org/MPL/2.0/. + +Covered Software is provided under this License on an "as is" +basis, without warranty of any kind, either expressed, implied, or +statutory, including, without limitation, warranties that the +Covered Software is free of defects, merchantable, fit for a +particular purpose or non-infringing. +See the Mozilla Public License v. 2.0 for more details. + +For more details, see http://www.mrtrix.org/. + diff --git a/docs/reference/commands/dwibiasnormmask.rst b/docs/reference/commands/dwibiasnormmask.rst new file mode 100644 index 0000000000..ba261fb0a0 --- /dev/null +++ b/docs/reference/commands/dwibiasnormmask.rst @@ -0,0 +1,102 @@ +.. _dwibiasnormmask: + +dwibiasnormmask +=============== + +Synopsis +-------- + +Iteratively perform masking, RF estimation, CSD, bias field removal, and mask revision on a DWI series + +Usage +----- + +:: + + dwibiasnormmask input output [ options ] + +- *input*: The input DWI series to be corrected +- *output*: The output corrected image series + +Description +----------- + +DWI brain masking, response function estimation and bias field correction are inter-related steps and errors in each may influence other steps. This script first derives a brain mask (either with a provided algorithm or by thresholding a balanced tissue sum image), and then performs response function estimation, multi-tissue CSD (with a lower lmax than the dwi2fod default, for speed), and mtnormalise to remove bias field, before the DWI brain mask is recalculated. These steps are performed iteratively until either a maximum number of iterations or until the brain masks converge. If SynthStrip is installed, it will be used to derive brain masks. + +Options +------- + +- **-output_mask** Output the final brain mask to the specified file + +- **-max_iters** The maximum number of iterations. The default is 2iterations. More iterations may lead to a better mask, but for some problematic data this may lead to the masks diverging (in which case a warning is issued). + +- **-lmax** The maximum spherical harmonic order for the output FODs. the value is passed to the dwi2fod command and should be provided in the format which it expects (the default value is "4,0,0" for multi-shell and "4,0" for single-shell data) + +Options for importing the diffusion gradient table +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-grad** Provide the diffusion gradient table in MRtrix format + +- **-fslgrad bvecs bvals** Provide the diffusion gradient table in FSL bvecs/bvals format + +Options to specify an initial mask or tissue sum image threshold in the iterative algorithm +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-threshold** Threshold on the total tissue density image used to derive the brain mask. the default is 0.5 + +- **-mask_init image** Provide an initial mask to the algorithm and skip the initial masking + +Additional standard options for Python scripts +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-nocleanup** do not delete intermediate files during script execution, and do not delete scratch directory at script completion. + +- **-scratch /path/to/scratch/** manually specify the path in which to generate the scratch directory. + +- **-continue ** continue the script from a previous execution; must provide the scratch directory path, and the name of the last successfully-generated file. + +Standard options +^^^^^^^^^^^^^^^^ + +- **-info** display information messages. + +- **-quiet** do not display information messages or progress status. Alternatively, this can be achieved by setting the MRTRIX_QUIET environment variable to a non-empty string. + +- **-debug** display debugging messages. + +- **-force** force overwrite of output files. + +- **-nthreads number** use this number of threads in multi-threaded applications (set to 0 to disable multi-threading). + +- **-config key value** *(multiple uses permitted)* temporarily set the value of an MRtrix config file entry. + +- **-help** display this information page and exit. + +- **-version** display version information and exit. + +References +^^^^^^^^^^ + +Tournier, J.-D.; Smith, R. E.; Raffelt, D.; Tabbara, R.; Dhollander, T.; Pietsch, M.; Christiaens, D.; Jeurissen, B.; Yeh, C.-H. & Connelly, A. MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation. NeuroImage, 2019, 202, 116137 + +-------------- + + + +**Author:** Robert E. Smith (robert.smith@florey.edu.au) and Arshiya Sangchooli (asangchooli@student.unimelb.edu.au) + +**Copyright:** Copyright (c) 2008-2023 the MRtrix3 contributors. + +This Source Code Form is subject to the terms of the Mozilla Public +License, v. 2.0. If a copy of the MPL was not distributed with this +file, You can obtain one at http://mozilla.org/MPL/2.0/. + +Covered Software is provided under this License on an "as is" +basis, without warranty of any kind, either expressed, implied, or +statutory, including, without limitation, warranties that the +Covered Software is free of defects, merchantable, fit for a +particular purpose or non-infringing. +See the Mozilla Public License v. 2.0 for more details. + +For more details, see http://www.mrtrix.org/. + diff --git a/docs/reference/commands/dwinormalise.rst b/docs/reference/commands/dwinormalise.rst index a228bb5df9..0039a4ef91 100644 --- a/docs/reference/commands/dwinormalise.rst +++ b/docs/reference/commands/dwinormalise.rst @@ -15,7 +15,7 @@ Usage dwinormalise algorithm [ options ] ... -- *algorithm*: Select the algorithm to be used to complete the script operation; additional details and options become available once an algorithm is nominated. Options are: group, individual +- *algorithm*: Select the algorithm to be used to complete the script operation; additional details and options become available once an algorithm is nominated. Options are: group, individual, mtnorm Description ----------- @@ -257,3 +257,102 @@ See the Mozilla Public License v. 2.0 for more details. For more details, see http://www.mrtrix.org/. +.. _dwinormalise_mtnorm: + +dwinormalise mtnorm +=================== + +Synopsis +-------- + +Performs a global DWI intensity normalisation, normalising mean b=0 CSF signal intensity + +Usage +----- + +:: + + dwinormalise mtnorm input_dwi output [ options ] + +- *input_dwi*: The input DWI series +- *output*: The normalised DWI series + +Description +----------- + +This script first performs response function estimation and multi-tissue CSD within a provided brain mask(with a lower lmax than the dwi2fod default, for higher speed) before using the global scaling factor derived by mtnormalise to scale the intensity of the DWI series such that the mean b=0 CSF signal intensity is approximately close to some reference value (by default 1000).Note that if the -mask command-line option is not specified, the MRtrix3 command dwi2mask will automatically be called to derive a mask that will be passed to the relevant bias field estimation command. More information on mask derivation from DWI data can be found at the following link: +https://mrtrix.readthedocs.io/en/3.0.4/dwi_preprocessing/masking.html + +Options +------- + +Options for importing the diffusion gradient table +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-grad** Provide the diffusion gradient table in MRtrix format + +- **-fslgrad bvecs bvals** Provide the diffusion gradient table in FSL bvecs/bvals format + +Options specific to the "mtnorm" algorithm +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-target** the threshold on the balanced tissue sum image used to derive the brain mask. the default is 0.5 + +- **-lmax** the maximum spherical harmonic order for the output FODs. The value is passed to the dwi2fod command and should be provided in the format which it expects (Default is "4,0,0" for multi-shell and "4,0" for single-shell data) + +- **-mask image** Manually provide a mask image for normalisation + +Additional standard options for Python scripts +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-nocleanup** do not delete intermediate files during script execution, and do not delete scratch directory at script completion. + +- **-scratch /path/to/scratch/** manually specify the path in which to generate the scratch directory. + +- **-continue ** continue the script from a previous execution; must provide the scratch directory path, and the name of the last successfully-generated file. + +Standard options +^^^^^^^^^^^^^^^^ + +- **-info** display information messages. + +- **-quiet** do not display information messages or progress status. Alternatively, this can be achieved by setting the MRTRIX_QUIET environment variable to a non-empty string. + +- **-debug** display debugging messages. + +- **-force** force overwrite of output files. + +- **-nthreads number** use this number of threads in multi-threaded applications (set to 0 to disable multi-threading). + +- **-config key value** *(multiple uses permitted)* temporarily set the value of an MRtrix config file entry. + +- **-help** display this information page and exit. + +- **-version** display version information and exit. + +References +^^^^^^^^^^ + +Tournier, J.-D.; Smith, R. E.; Raffelt, D.; Tabbara, R.; Dhollander, T.; Pietsch, M.; Christiaens, D.; Jeurissen, B.; Yeh, C.-H. & Connelly, A. MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation. NeuroImage, 2019, 202, 116137 + +-------------- + + + +**Author:** Robert E. Smith (robert.smith@florey.edu.au) and Arshiya Sangchooli (asangchooli@student.unimelb.edu.au) + +**Copyright:** Copyright (c) 2008-2023 the MRtrix3 contributors. + +This Source Code Form is subject to the terms of the Mozilla Public +License, v. 2.0. If a copy of the MPL was not distributed with this +file, You can obtain one at http://mozilla.org/MPL/2.0/. + +Covered Software is provided under this License on an "as is" +basis, without warranty of any kind, either expressed, implied, or +statutory, including, without limitation, warranties that the +Covered Software is free of defects, merchantable, fit for a +particular purpose or non-infringing. +See the Mozilla Public License v. 2.0 for more details. + +For more details, see http://www.mrtrix.org/. + diff --git a/docs/reference/commands_list.rst b/docs/reference/commands_list.rst index 9cd39f485d..545a0da124 100644 --- a/docs/reference/commands_list.rst +++ b/docs/reference/commands_list.rst @@ -36,6 +36,7 @@ List of MRtrix3 commands commands/dwi2response commands/dwi2tensor commands/dwibiascorrect + commands/dwibiasnormmask commands/dwicat commands/dwidenoise commands/dwiextract @@ -167,6 +168,7 @@ List of MRtrix3 commands |python.png|, :ref:`dwi2response`, "Estimate response function(s) for spherical deconvolution" |cpp.png|, :ref:`dwi2tensor`, "Diffusion (kurtosis) tensor estimation" |python.png|, :ref:`dwibiascorrect`, "Perform B1 field inhomogeneity correction for a DWI volume series" + |python.png|, :ref:`dwibiasnormmask`, "Iteratively perform masking, RF estimation, CSD, bias field removal, and mask revision on a DWI series" |python.png|, :ref:`dwicat`, "Concatenating multiple DWI series accounting for differential intensity scaling" |cpp.png|, :ref:`dwidenoise`, "dMRI noise level estimation and denoising using Marchenko-Pastur PCA" |cpp.png|, :ref:`dwiextract`, "Extract diffusion-weighted volumes, b=0 volumes, or certain shells from a DWI dataset" diff --git a/lib/mrtrix3/dwi2mask/mtnorm.py b/lib/mrtrix3/dwi2mask/mtnorm.py index 8c52459724..014b679ff9 100644 --- a/lib/mrtrix3/dwi2mask/mtnorm.py +++ b/lib/mrtrix3/dwi2mask/mtnorm.py @@ -70,7 +70,7 @@ def __init__(self, name): dwi_image = 'input.mif' tissues = [Tissue('WM'), Tissue('GM'), Tissue('CSF')] - + app.debug('Estimating response function using initial brain mask...') run.command('dwi2response dhollander ' + dwi_image @@ -100,7 +100,7 @@ def __init__(self, name): app.debug('Deriving final brain mask by thresholding tissue sum image...') new_dwi_mask_image = 'dwi_mask' + '.mif' - TISSUESUM_THRESHOLD = app.ARGS.threshold / math.sqrt(4.0 * math.pi) + tissuesum_threshold = app.ARGS.threshold / math.sqrt(4.0 * math.pi) run.command('mrconvert ' + tissues[0].fod @@ -109,7 +109,7 @@ def __init__(self, name): + ' '.join(tissue.fod for tissue in tissues[1:]) + ' sum - |' + ' mrthreshold - -abs ' - + str(TISSUESUM_THRESHOLD) + + str(tissuesum_threshold) + ' - |' + ' maskfilter - connect -largest - |' + ' mrcalc 1 - -sub - -datatype bit |' @@ -122,4 +122,3 @@ def __init__(self, name): app.cleanup([tissue.fod for tissue in tissues]) return dwi_mask_image - diff --git a/lib/mrtrix3/dwibiascorrect/mtnorm.py b/lib/mrtrix3/dwibiascorrect/mtnorm.py index 35874521dd..22d28b3af8 100644 --- a/lib/mrtrix3/dwibiascorrect/mtnorm.py +++ b/lib/mrtrix3/dwibiascorrect/mtnorm.py @@ -13,10 +13,8 @@ # # For more details, see http://www.mrtrix.org/. -import math from mrtrix3 import app, run, image, path - LMAXES_MULTI = '4,0,0' LMAXES_SINGLE = '4,0' @@ -61,7 +59,7 @@ def __init__(self, name): self.fod_norm = 'FODnorm_' + name + '.mif' tissues = [Tissue('WM'), Tissue('GM'), Tissue('CSF')] - + app.debug('Estimating response function using brain mask...') run.command('dwi2response dhollander ' + 'in.mif ' @@ -85,7 +83,7 @@ def __init__(self, name): + ' ' + ' '.join(tissue.tissue_rf + ' ' + tissue.fod for tissue in tissues)) - + app.debug('correcting bias field...') run.command('maskfilter' + ' mask.mif' diff --git a/lib/mrtrix3/dwinormalise/mtnorm.py b/lib/mrtrix3/dwinormalise/mtnorm.py index 30b9bb3948..4526e16888 100644 --- a/lib/mrtrix3/dwinormalise/mtnorm.py +++ b/lib/mrtrix3/dwinormalise/mtnorm.py @@ -58,10 +58,10 @@ def execute(): #pylint: disable=unused-variable grad_option = ' -grad ' + path.from_user(app.ARGS.grad) elif app.ARGS.fslgrad: grad_option = ' -fslgrad ' + path.from_user(app.ARGS.fslgrad[0]) + ' ' + path.from_user(app.ARGS.fslgrad[1]) - + app.make_scratch_dir() app.goto_scratch_dir() - + run.command('mrconvert ' + path.from_user(app.ARGS.input_dwi) + ' ' @@ -93,7 +93,7 @@ def __init__(self, name): self.fod_norm = 'FODnorm_' + name + '.mif' tissues = [Tissue('WM'), Tissue('GM'), Tissue('CSF')] - + app.debug('Estimating response function using initial brain mask...') run.command('dwi2response dhollander ' + 'input.mif ' @@ -143,9 +143,3 @@ def __init__(self, name): + ' -mult ' + path.from_user(app.ARGS.output), force=app.FORCE_OVERWRITE) - - - - - - From 0b82d20885658cbdee43f2f7360dd403477417ad Mon Sep 17 00:00:00 2001 From: Arshiyas Date: Tue, 21 Feb 2023 16:04:31 +1100 Subject: [PATCH 06/24] Some tests to check dwibiasnormmask, no new files --- testing/scripts/tests/dwibiasnormmask | 4 ++++ 1 file changed, 4 insertions(+) create mode 100644 testing/scripts/tests/dwibiasnormmask diff --git a/testing/scripts/tests/dwibiasnormmask b/testing/scripts/tests/dwibiasnormmask new file mode 100644 index 0000000000..5c5c2e957a --- /dev/null +++ b/testing/scripts/tests/dwibiasnormmask @@ -0,0 +1,4 @@ +mkdir -p ../tmp/dwibiasnormmask && mrconvert BIDS/sub-01/dwi/sub-01_dwi.nii.gz -fslgrad BIDS/sub-01/dwi/sub-01_dwi.bvec BIDS/sub-01/dwi/sub-01_dwi.bval tmp-sub-01_dwi.mif -force && dwibiasnormmask tmp-sub-01_dwi.mif ../tmp/dwibiasnormmask/default_out.mif -output_mask ../tmp/dwibiasnormmask/default_mask.mif -force +dwibiasnormmask tmp-sub-01_dwi.mif -max_iters 3 ../tmp/dwibiasnormmask/3iters_out.mif -output_mask ../tmp/dwibiasnormmask/3iters_mask.mif -force +dwibiasnormmask tmp-sub-01_dwi.mif -lmax 6,0,0 ../tmp/dwibiasnormmask/lmax600_out.mif -output_mask ../tmp/dwibiasnormmask/lmax600_mask.mif -force + From fd3ba556b0fea65677dd1ac22b9103c43dedf6af Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Wed, 22 Feb 2023 12:27:43 +1100 Subject: [PATCH 07/24] Large modifications to dwibiasnormmask --- bin/dwibiasnormmask | 342 ++++++++++++-------- docs/reference/commands/dwibiasnormmask.rst | 46 ++- docs/reference/commands_list.rst | 2 +- lib/mrtrix3/app.py | 2 +- 4 files changed, 248 insertions(+), 144 deletions(-) diff --git a/bin/dwibiasnormmask b/bin/dwibiasnormmask index 4efbebdc7c..bf0d755361 100755 --- a/bin/dwibiasnormmask +++ b/bin/dwibiasnormmask @@ -15,55 +15,124 @@ # # For more details, see http://www.mrtrix.org/. -import math, shutil +import math, os, shutil DWIBIASCORRECT_MAX_ITERS = 2 -LMAXES_MULTI = '4,0,0' -LMAXES_SINGLE = '4,0' - +LMAXES_MULTI = [4,0,0] +LMAXES_SINGLE = [4,0] +REFERENCE_INTENSITY = 1000.0 +# TODO Add config file entry to specify user preferred default +MASK_ALGOS = ['dwi2mask', 'fslbet', 'hdbet', 'mrthreshold', 'synthstrip', 'threshold'] +MASK_ALGO_DEFAULT = 'threshold' +DICE_COEFF_DEFAULT = 1.0 - 1e-3 def usage(cmdline): #pylint: disable=unused-variable from mrtrix3 import app #pylint: disable=no-name-in-module, import-outside-toplevel cmdline.set_author('Robert E. Smith (robert.smith@florey.edu.au) and Arshiya Sangchooli (asangchooli@student.unimelb.edu.au)') - cmdline.set_synopsis('Iteratively perform masking, RF estimation, CSD, bias field removal, and mask revision on a DWI series') - cmdline.add_description('DWI brain masking, response function estimation and bias field correction are inter-related steps and errors ' - 'in each may influence other steps. This script first derives a brain mask (either with a provided algorithm ' - ' or by thresholding a balanced tissue sum image), and then performs response function estimation, multi-tissue ' - 'CSD (with a lower lmax than the dwi2fod default, for speed), and mtnormalise to remove bias field, before the ' - 'DWI brain mask is recalculated. These steps are performed iteratively until either a maximum number of iterations ' - 'or until the brain masks converge. If SynthStrip is installed, it will be used to derive brain masks.') - cmdline.add_argument('input', help='The input DWI series to be corrected') - cmdline.add_argument('output', help='The output corrected image series') - cmdline.add_argument('-output_mask', help='Output the final brain mask to the specified file') - cmdline.add_argument('-max_iters', type=int, default=DWIBIASCORRECT_MAX_ITERS, - help='The maximum number of iterations. The default is ' + str(DWIBIASCORRECT_MAX_ITERS) + - 'iterations. More iterations may lead to a better mask, but for some problematic data ' - 'this may lead to the masks diverging (in which case a warning is issued).') - cmdline.add_argument('-lmax', type=str, - help='The maximum spherical harmonic order for the output FODs. the value is passed to the dwi2fod command ' - 'and should be provided in the format which it expects (the default value is "' + str(LMAXES_MULTI) + - '" for multi-shell and "' + str(LMAXES_SINGLE) + '" for single-shell data)') - mask_options = cmdline.add_argument_group('Options to specify an initial mask or tissue sum image threshold in the iterative algorithm') - mask_options.add_argument('-threshold', type=float, default=0.5, - help='Threshold on the total tissue density image used to derive the brain mask. the default is 0.5') - mask_options.add_argument('-mask_init', metavar='image', - help='Provide an initial mask to the algorithm and skip the initial masking') + cmdline.set_synopsis('Perform a combination of bias field correction, intensity normalisation, and mask derivation, for DWI data') + cmdline.add_description('DWI bias field correction, intensity normalisation and masking are inter-related steps, and errors ' + 'in each may influence other steps. This script is designed to perform all of these steps in an integrated ' + 'iterative fashion, with the intention of making all steps more robust.') + cmdline.add_description('The operation of the algorithm is as follows. An initial mask is defined, either using the default dwi2mask ' + 'algorithm or as provided by the user. Based on this mask, a sequence of response function estimation, ' + 'multi-shell multi-tissue CSD, bias field correction (using the mtnormalise command), and intensity ' + 'normalisation is performed. The default dwi2mask algorithm is then re-executed on the bias-field-corrected ' + 'DWI series. This sequence of steps is then repeated based on the revised mask, until either a convergence ' + 'criterion or some number of maximum iterations is reached.') + cmdline.add_description('The MRtrix3 mtnormalise command is used to estimate information relating to bias field and intensity ' + 'normalisation. However its usage in this context is different to its conventional usage. Firstly, ' + 'while the corrected ODF images are typically used directly following invocation of this command, ' + 'here the estimated bias field and scaling factors are instead used to apply the relevant corrections to ' + 'the originating DWI data. Secondly, the global intensity scaling that is calculated and applied is ' + 'typically based on achieving close to a unity sum of tissue signal fractions throughout the masked region. ' + 'Here, it is instead the b=0 signal in CSF that forms the reference for this global intensity scaling; ' + 'this is calculated based on the estimated CSF response function and the tissue-specific intensity ' + 'scaling (this is calculated internally by mtnormalise as part of its optimisation process, but typically ' + 'subsequently discarded in favour of a single scaling factor for all tissues)') + cmdline.add_description('The ODFs estimated within this optimisation procedure are by default of lower maximal spherical harmonic ' + 'degree than what would be advised for analysis. This is done for computational efficiency. This ' + 'behaviour can be modified through the -lmax command-line option.') + cmdline.add_description('By default, the optimisation procedure will terminate after only two iterations. This is done because ' + 'it has been observed for some data / configurations that additional iterations can lead to unstable ' + 'divergence and erroneous results for bias field estimation and masking. For other configurations, ' + 'it may be preferable to use a greater number of iterations, and allow the iterative algorithm to ' + 'converge to a stable solution. This can be controlled via the -max_iters command-line option.') + cmdline.add_description('Within the optimisation algorithm, derivation of the mask may potentially be performed differently to ' + 'a conventional mask derivation that is based on a DWI series (where, in many instances, it is actually ' + 'only the mean b=0 image that is used). Here, the image corresponding to the sum of tissue signal fractions ' + 'following spherical deconvolution / bias field correction / intensity normalisation is also available, ' + 'and this can potentially be used for mask derivation. Available options are as follows. ' + '"dwi2mask": Use the MRtrix3 command dwi2mask on the bias-field-corrected DWI series ' + '(ie. do not use the ODF tissue sum image for mask derivation); ' + 'the algorithm to be invoked can be controlled by the user via the MRtrix config file entry "Dwi2maskAlgorithm". ' + '"fslbet": Invoke the FSL command "bet". ' + '"hdbet": Invoke the HD-BET command. ' + '"mrthreshold": Invoke the MRtrix3 command "mrthreshold", where an appropriate threshold value will be determined automatically ' + ' (and some heuristic cleanup of the resulting mask will be performed). ' + '"synthstrip": Invoke the FreeSurfer SynthStrip method. ' + '"threshold": Apply a fixed partial volume threshold of 0.5 to the tissue ODF sum image ' + ' (and some heuristic cleanup of the resulting mask will be performed).') + cmdline.add_argument('input', help='The input DWI series to be corrected') + cmdline.add_argument('output_dwi', help='The output corrected DWI series') + cmdline.add_argument('output_mask', help='The output DWI mask') + output_options = cmdline.add_argument_group('Options that modulate the outputs of the script') + output_options.add_argument('-output_bias', metavar='image', help='Export the final estimated bias field to an image') + output_options.add_argument('-output_tissuesum', metavar='image', help='Export the tissue sum image that was used to generate the final mask') + output_options.add_argument('-reference', type=float, metavar='value', default=REFERENCE_INTENSITY, help='Set the target CSF b=0 intensity in the output DWI series (default: ' + str(REFERENCE_INTENSITY) + ')') + internal_options = cmdline.add_argument_group('Options relevant to the internal optimisation procedure') + internal_options.add_argument('-dice', type=float, default=DICE_COEFF_DEFAULT, metavar='value', + help='Set the Dice coefficient threshold for similarity of masks between sequential iterations that will ' + 'result in termination due to convergence; default = ' + str(DICE_COEFF_DEFAULT)) + internal_options.add_argument('-init_mask', metavar='image', + help='Provide an initial mask for the first iteration of the algorithm ' + '(if not provided, the default dwi2mask algorithm will be used)') + internal_options.add_argument('-max_iters', type=int, default=DWIBIASCORRECT_MAX_ITERS, metavar='count', + help='The maximum number of iterations (see Description); default is ' + str(DWIBIASCORRECT_MAX_ITERS) + '; ' + 'set to 0 to proceed until convergence') + internal_options.add_argument('-mask_algo', choices=MASK_ALGOS, default=MASK_ALGO_DEFAULT, metavar='algorithm', + help='The algorithm to use for mask estimation, potentially based on the ODF sum image (see Description); default: ' + MASK_ALGO_DEFAULT) + internal_options.add_argument('-lmax', metavar='value', + help='The maximum spherical harmonic degree for the estimated FODs (see Description); ' + 'defaults are "' + ','.join(str(item) for item in LMAXES_MULTI) + '" for multi-shell and "' + ','.join(str(item) for item in LMAXES_SINGLE) + '" for single-shell data)') app.add_dwgrad_import_options(cmdline) + def execute(): #pylint: disable=unused-variable from mrtrix3 import CONFIG, MRtrixError #pylint: disable=no-name-in-module, import-outside-toplevel - from mrtrix3 import app, run, image, matrix, path #pylint: disable=no-name-in-module, import-outside-toplevel - - app.check_output_path(app.ARGS.output) + from mrtrix3 import app, fsl, image, matrix, path, run #pylint: disable=no-name-in-module, import-outside-toplevel + + # Check user inputs + if app.ARGS.max_iters < 0: + raise MRtrixError('Maximum number of iterations must be a non-negative integer') + lmax = None + if app.ARGS.lmax: + try: + lmax = [int(i) for i in app.ARGS.lmax.split(',')] + except ValueError: + raise MRtrixError('Values provided to -lmax option must be a comma-separated list of integers') + if any(value < 0 or value % 2 for value in lmax): + raise MRtrixError('lmax values must be non-negative even integers') + if len(lmax) not in [2, 3]: + raise MRtrixError('Length of lmax vector expected to be either 2 or 3') + if app.ARGS.dice <= 0.0 or app.ARGS.dice > 1.0: + raise MRtrixError('Dice coefficient for convergence detection must lie in the range (0.0, 1.0]') + + # Check mask algorithm, including availability of external software if necessary + for mask_algo, software, command in [('fslbet', 'FSL', 'bet'), ('hdbet', 'HD-BET', 'hd-bet'), ('synthstrip', 'FreeSurfer', 'mri_synthstrip')]: + if app.ARGS.mask_algo == mask_algo and not shutil.which(command): + raise MRtrixError(software + ' command "' + command + '" not found; cannot use for internal mask calculations') + + app.check_output_path(app.ARGS.output_dwi) + app.check_output_path(app.ARGS.output_mask) app.make_scratch_dir() grad_import_option = app.read_dwgrad_import_options() run.command('mrconvert ' + path.from_user(app.ARGS.input) + ' ' + path.to_scratch('input.mif') + grad_import_option) - if app.ARGS.mask_init: - run.command('mrconvert ' + path.from_user(app.ARGS.mask_init) + ' ' + if app.ARGS.init_mask: + run.command('mrconvert ' + path.from_user(app.ARGS.init_mask) + ' ' + path.to_scratch('dwi_mask_init.mif') + ' -datatype bit') app.goto_scratch_dir() @@ -87,32 +156,26 @@ def execute(): #pylint: disable=unused-variable for value in image.mrinfo('input.mif', 'shell_bvalues') \ .strip().split()] multishell = (len(bvalues) > 2) + if lmax is None: + lmax = LMAXES_MULTI if multishell else LMAXES_SINGLE + elif len(lmax) == 3 and not multishell: + raise MRtrixError('User specified 3 lmax values for three-tissue decomposition, but input DWI is not multi-shell') - # Get the SynthStrip path if it's installed - synthstrip_cmd = shutil.which('mri_synthstrip') - if not synthstrip_cmd: - synthstrip_cmd = shutil.which('synthstrip-singularity') + # Create a mask of voxels where the input data contain positive values; + # we want to make sure that these never end up included in the output mask + run.command('mrmath input.mif max - -axis 3 | ' + 'mrthreshold - -abs 0 -comparison gt input_pos_mask.mif') # Generate an initial brain mask if required, or check the initial mask if provided by the user - if app.ARGS.mask_init: + if app.ARGS.init_mask: if not image.match('input.mif', 'dwi_mask_init.mif', up_to_dim=3): raise MRtrixError('Provided mask image does not match input DWI') else: app.debug('Performing intial DWI brain masking') - if synthstrip_cmd: - run.command('dwiextract input.mif -bzero - | ' - 'mrmath - mean - -axis 3 | ' - 'mrconvert - bzero.nii -strides +1,+2,+3') - run.command(synthstrip_cmd + ' -i bzero.nii -b 0 --mask dwi_mask_init.nii') - run.command('mrconvert dwi_mask_init.nii dwi_mask_init.mif -datatype bit') - app.cleanup('bzero.nii') - app.cleanup('dwi_mask_init.nii') - else: - run.command('dwi2mask ' - + CONFIG.get('Dwi2maskAlgorithm', 'legacy') - + ' input.mif dwi_mask_init.mif') + run.command('dwi2mask ' + + CONFIG.get('Dwi2maskAlgorithm', 'legacy') + + ' input.mif dwi_mask_init.mif') - dwi_mask_image = 'dwi_mask_init.mif' # Combined RF estimation / CSD / mtnormalise / mask revision class Tissue(object): #pylint: disable=useless-object-inheritance @@ -124,21 +187,24 @@ def execute(): #pylint: disable=unused-variable self.fod_norm = 'FODnorm_' + name + iter_string + '.mif' - app.debug('Commencing iterative DWI bias field correction and brain masking ' - 'with a maximum of ' + str(app.ARGS.max_iters) + ' iterations') + app.debug('Commencing iterative DWI bias field correction and brain masking with ' + + ('a maximum of ' + str(app.ARGS.max_iters) if app.ARGS.max_iters else 'no limit on number of ') + ' iterations') dwi_image = 'input.mif' + dwi_mask_image = 'dwi_mask_init.mif' + tissue_sum_image = None iteration = 0 step = 'initialisation' - dice_coefficient = 0.0 + prev_dice_coefficient = 0.0 def msg(): return 'Iteration {0}; {1} step; previous Dice coefficient {2}' \ - .format(iteration, step, dice_coefficient) + .format(iteration, step, prev_dice_coefficient) progress = app.ProgressBar(msg) - for iteration in range(0, app.ARGS.max_iters): - iter_string = '_iter' + str(iteration+1) + iteration = 1 + while True: + iter_string = '_iter' + str(iteration) tissues = [Tissue('WM', iteration), Tissue('GM', iteration), Tissue('CSF', iteration)] @@ -153,51 +219,42 @@ def execute(): #pylint: disable=unused-variable + ' '.join(tissue.tissue_rf for tissue in tissues)) - # Determine the lmax values to use for the CSD, also remove GM if we can't deal with it - if multishell: - lmaxes = LMAXES_MULTI - else: - lmaxes = LMAXES_SINGLE + # Immediately remove GM if we can't deal with it + if not multishell: app.cleanup(tissues[1].tissue_rf) tissues = tissues[::2] - if app.ARGS.lmax: - lmaxes = app.ARGS.lmax - if (lmaxes.count(",") + 1) != len(tissues): - raise MRtrixError('Number of lmax values provided is not appropriate for number of shells') step = 'dwi2fod' progress.increment() - app.debug('Performing CSD with lmax values: ' + lmaxes) + app.debug('Performing CSD with lmax values: ' + ','.join(str(item) for item in lmax)) run.command('dwi2fod msmt_csd ' + dwi_image - + ' -lmax ' + lmaxes + + ' -lmax ' + ','.join(str(item) for item in lmax) + ' ' + ' '.join(tissue.tissue_rf + ' ' + tissue.fod_init for tissue in tissues)) + step = 'maskfilter' + progress.increment() + eroded_mask = os.path.splitext(dwi_mask_image)[0] + '_eroded.mif' + run.command('maskfilter ' + dwi_mask_image + ' erode ' + eroded_mask) + step = 'mtnormalise' progress.increment() field_path = 'field' + iter_string + '.mif' factors_path = 'factors' + iter_string + '.txt' - # It might not be nescessary to erode the mask if masking is already restrictive - if synthstrip_cmd: - app.debug('Masking is likely restrictive enough, not eroding mask') - run.command('cp ' + dwi_mask_image + ' temp_mask.mif') - else: - app.debug('Masking may not be restrictive enough, eroding mask') - run.command('maskfilter ' + dwi_mask_image + ' erode temp_mask.mif') - - run.command('mtnormalise -mask temp_mask.mif -balanced' + run.command('mtnormalise -balanced' + + ' -mask ' + eroded_mask + ' -check_norm ' + field_path + ' -check_factors ' + factors_path + ' ' + ' '.join(tissue.fod_init + ' ' + tissue.fod_norm for tissue in tissues)) app.cleanup([tissue.fod_init for tissue in tissues]) - app.cleanup('temp_mask.mif') + app.cleanup(eroded_mask) - app.debug('Iteration ' + str(iteration+1) + ', ' + app.debug('Iteration ' + str(iteration) + ', ' + 'applying estimated bias field and appropiate scaling factor...') csf_rf = matrix.load_matrix(tissues[-1].tissue_rf) csf_rf_bzero_lzero = csf_rf[0][0] @@ -205,7 +262,7 @@ def execute(): #pylint: disable=unused-variable balance_factors = matrix.load_vector(factors_path) csf_balance_factor = balance_factors[-1] app.cleanup(factors_path) - scale_multiplier = (1000.0 * math.sqrt(4.0*math.pi)) / \ + scale_multiplier = (app.ARGS.reference * math.sqrt(4.0*math.pi)) / \ (csf_rf_bzero_lzero / csf_balance_factor) new_dwi_image = 'dwi' + iter_string + '.mif' run.command('mrcalc ' + dwi_image + ' ' @@ -217,33 +274,53 @@ def execute(): #pylint: disable=unused-variable old_dwi_image = dwi_image dwi_image = new_dwi_image + old_tissue_sum_image = tissue_sum_image tissue_sum_image = 'tissue_sum' + iter_string + '.mif' - app.debug('Iteration ' + str(iteration+1) + ', ' + app.debug('Iteration ' + str(iteration) + ', ' + 'revising brain mask...') - new_dwi_mask_image = 'dwi_mask' + iter_string + '.mif' + step = 'masking' + progress.increment() + run.command('mrconvert ' + tissues[0].fod_norm + ' -coord 3 0 - |' + ' mrmath - ' + ' '.join(tissue.fod_norm for tissue in tissues[1:]) - + ' sum ' + + ' sum - | ' + + 'mrcalc - ' + str(math.sqrt(4.0 * math.pi)) + ' -mult ' + tissue_sum_image) - app.cleanup([tissue.fod_norm for tissue in tissues]) - if synthstrip_cmd: - run.command('mrconvert ' + tissue_sum_image + ' tissue.nii -force') - run.command(synthstrip_cmd + ' -i tissue.nii -b 0 --mask dwi_mask_temp.nii') - run.command('mrconvert dwi_mask_temp.nii ' + new_dwi_mask_image + ' -datatype bit') - app.cleanup('tissue.nii') - app.cleanup('dwi_mask_temp.nii') - else: - tissuesum_threshold = app.ARGS.threshold / math.sqrt(4.0 * math.pi) + new_dwi_mask_image = 'dwi_mask' + iter_string + '.mif' + tissue_sum_image_nii = None + new_dwi_mask_image_nii = None + if app.ARGS.mask_algo in ['fslbet', 'hdbet', 'synthstrip']: + tissue_sum_image_nii = os.path.splitext(tissue_sum_image)[0] + '.nii' + run.command('mrconvert ' + tissue_sum_image + ' ' + tissue_sum_image_nii) + new_dwi_mask_image_nii = os.path.splitext(new_dwi_mask_image)[0] + '.nii' + if app.ARGS.mask_algo == 'dwi2mask': + run.command('dwi2mask ' + CONFIG.get('Dwi2maskAlgorithm', 'legacy') + ' ' + new_dwi_image + ' ' + new_dwi_mask_image) + elif app.ARGS.mask_algo == 'fslbet': + run.command('bet ' + tissue_sum_image_nii + ' ' + new_dwi_mask_image_nii + ' -R -m') + app.cleanup(fsl.find_image(os.path.splitext(new_dwi_mask_image_nii)[0])) + new_dwi_mask_image_nii = fsl.find_image(os.path.splitext(new_dwi_mask_image_nii)[0] + '_mask') + run.command('mrcalc ' + new_dwi_mask_image_nii + ' input_pos_mask.mif -mult ' + new_dwi_mask_image) + elif app.ARGS.mask_algo == 'hdbet': + try: + run.command('hd-bet -i ' + tissue_sum_image_nii) + except run.MRtrixCmdError as e_gpu: + try: + run.command('hd-bet -i ' + tissue_sum_image_nii + ' -device cpu -mode fast -tta 0') + except run.MRtrixCmdError as e_cpu: + raise run.MRtrixCmdError('hd-bet', 1, e_gpu.stdout + e_cpu.stdout, e_gpu.stderr + e_cpu.stderr) + new_dwi_mask_image_nii = os.path.splitext(tissue_sum_image)[0] + '_bet_mask.nii.gz' + run.command('mrcalc ' + new_dwi_mask_image_nii + ' input_pos_mask.mif -mult ' + new_dwi_mask_image) + elif app.ARGS.mask_algo in ['mrthreshold', 'threshold']: + mrthreshold_abs_option = ' -abs 0.5' if app.ARGS.mask_algo == 'threshold' else '' run.command('mrthreshold ' + tissue_sum_image - + ' -abs ' - + str(tissuesum_threshold) + + mrthreshold_abs_option + ' - |' + ' maskfilter - connect -largest - |' + ' mrcalc 1 - -sub - -datatype bit |' @@ -253,74 +330,81 @@ def execute(): #pylint: disable=unused-variable + ' mrcalc - input_pos_mask.mif -mult ' + new_dwi_mask_image + ' -datatype bit') - - + elif app.ARGS.mask_algo == 'synthstrip': + run.command('mri_synthstrip -i ' + tissue_sum_image_nii + ' --mask ' + new_dwi_mask_image_nii) + run.command('mrcalc ' + new_dwi_mask_image_nii + ' input_pos_mask.mif -mult ' + new_dwi_mask_image) + else: + assert False app.cleanup(tissue_sum_image) + if tissue_sum_image_nii: + app.cleanup(tissue_sum_image_nii) + if new_dwi_mask_image_nii: + app.cleanup(new_dwi_mask_image_nii) - step = 'mrcalc_mask' + step = 'mask_comparison' progress.increment() dwi_old_mask_count = image.statistics(dwi_mask_image, mask=dwi_mask_image).count dwi_new_mask_count = image.statistics(new_dwi_mask_image, mask=new_dwi_mask_image).count - app.debug('Old mask: ' + str(dwi_old_mask_count)) - app.debug('New mask: ' + str(dwi_new_mask_count)) + app.debug('Old mask voxel count: ' + str(dwi_old_mask_count)) + app.debug('New mask voxel count: ' + str(dwi_new_mask_count)) dwi_mask_overlap_image = 'dwi_mask_overlap' + iter_string + '.mif' - run.command('mrcalc ' - + dwi_mask_image - + ' ' - + new_dwi_mask_image - + ' -mult ' - + dwi_mask_overlap_image) + run.command(['mrcalc', dwi_mask_image, new_dwi_mask_image, '-mult', dwi_mask_overlap_image]) old_dwi_mask_image = dwi_mask_image dwi_mask_image = new_dwi_mask_image mask_overlap_count = image.statistics(dwi_mask_overlap_image, mask=dwi_mask_overlap_image).count - app.debug('Mask overlap: ' + str(mask_overlap_count)) + app.debug('Mask overlap voxel count: ' + str(mask_overlap_count)) - # Break if either the mask has converged properly, or if it starts to diverge after some convergence - old_dice_coefficient = dice_coefficient - dice_coefficient = 2.0 * mask_overlap_count / \ - (dwi_old_mask_count + dwi_new_mask_count) + new_dice_coefficient = 2.0 * mask_overlap_count / \ + (dwi_old_mask_count + dwi_new_mask_count) - if dice_coefficient > (1.0 - 1e-3): - app.debug('Exiting iterative loop due to mask convergence. Iteration: ' + str(iteration + 1)) + if new_dice_coefficient > app.ARGS.dice: + progress.done() + app.console('Exiting loop after %d iterations due to mask convergence (Dice coefficient = %f)' % iteration, new_dice_coefficient) app.cleanup(old_dwi_image) app.cleanup(old_dwi_mask_image) + app.cleanup(old_tissue_sum_image) break - if (dice_coefficient < old_dice_coefficient) and (old_dice_coefficient > (1.0 - 1e-2)): - app.debug('[WARNING] Mask divergence at iteration: ' + str(iteration + 1) - + ' using previous mask with dice coefficient: ' + str(old_dice_coefficient)) - app.cleanup(dwi_mask_image) + if new_dice_coefficient < prev_dice_coefficient: + progress.done() + app.warn('Mask divergence at iteration %d (Dice coefficient = %f); ' % (iteration, new_dice_coefficient) + + ' using mask from previous iteration') app.cleanup(dwi_image) + app.cleanup(dwi_mask_image) + app.cleanup(tissue_sum_image) dwi_image = old_dwi_image dwi_mask_image = old_dwi_mask_image + tissue_sum_image = old_tissue_sum_image break - if iteration == (DWIBIASCORRECT_MAX_ITERS - 1): - app.debug('[WARNING] Maximum iterations without convergence. Final dice coefficient: ' - + str(dice_coefficient)) + if iteration == app.ARGS.max_iters: + progress.done() + app.console('Terminating due to reaching maximum %d iterations; final Dice coefficient = %f' % (iteration, new_dice_coefficient)) app.cleanup(old_dwi_image) app.cleanup(old_dwi_mask_image) + app.cleanup(old_tissue_sum_image) + break + iteration += 1 + prev_dice_coefficient = new_dice_coefficient - # Processing completed; export - app.console('Processing completed after ' - + str(iteration + 1) - + ' iterations, writing results to output directory') - - run.command('mrconvert ' - + dwi_image + ' ' - + path.from_user(app.ARGS.output), + run.command(['mrconvert', dwi_image, path.from_user(app.ARGS.output_dwi, False)], + mrconvert_keyval=path.from_user(app.ARGS.input, False), force=app.FORCE_OVERWRITE) if app.ARGS.output_mask: - run.command('mrconvert ' - + dwi_mask_image + ' ' - + path.from_user(app.ARGS.output_mask), + run.command(['mrconvert', dwi_mask_image, path.from_user(app.ARGS.output_mask, False)], + mrconvert_keyval=path.from_user(app.ARGS.input, False), + force=app.FORCE_OVERWRITE) + + if app.ARGS.output_tissuesum: + run.command(['mrconvert', tissue_sum_image, path.from_user(app.ARGS.output_tissuesum, False)], + mrconvert_keyval=path.from_user(app.ARGS.input, False), force=app.FORCE_OVERWRITE) diff --git a/docs/reference/commands/dwibiasnormmask.rst b/docs/reference/commands/dwibiasnormmask.rst index ba261fb0a0..23da1312db 100644 --- a/docs/reference/commands/dwibiasnormmask.rst +++ b/docs/reference/commands/dwibiasnormmask.rst @@ -6,31 +6,36 @@ dwibiasnormmask Synopsis -------- -Iteratively perform masking, RF estimation, CSD, bias field removal, and mask revision on a DWI series +Perform a combination of bias field correction, intensity normalisation, and mask derivation, for DWI data Usage ----- :: - dwibiasnormmask input output [ options ] + dwibiasnormmask input output_dwi output_mask [ options ] - *input*: The input DWI series to be corrected -- *output*: The output corrected image series +- *output_dwi*: The output corrected DWI series +- *output_mask*: The output DWI mask Description ----------- -DWI brain masking, response function estimation and bias field correction are inter-related steps and errors in each may influence other steps. This script first derives a brain mask (either with a provided algorithm or by thresholding a balanced tissue sum image), and then performs response function estimation, multi-tissue CSD (with a lower lmax than the dwi2fod default, for speed), and mtnormalise to remove bias field, before the DWI brain mask is recalculated. These steps are performed iteratively until either a maximum number of iterations or until the brain masks converge. If SynthStrip is installed, it will be used to derive brain masks. +DWI bias field correction, intensity normalisation and masking are inter-related steps, and errors in each may influence other steps. This script is designed to perform all of these steps in an integrated iterative fashion, with the intention of making all steps more robust. -Options -------- +The operation of the algorithm is as follows. An initial mask is defined, either using the default dwi2mask algorithm or as provided by the user. Based on this mask, a sequence of response function estimation, multi-shell multi-tissue CSD, bias field correction (using the mtnormalise command), and intensity normalisation is performed. The default dwi2mask algorithm is then re-executed on the bias-field-corrected DWI series. This sequence of steps is then repeated based on the revised mask, until either a convergence criterion or some number of maximum iterations is reached. -- **-output_mask** Output the final brain mask to the specified file +The MRtrix3 mtnormalise command is used to estimate information relating to bias field and intensity normalisation. However its usage in this context is different to its conventional usage. Firstly, while the corrected ODF images are typically used directly following invocation of this command, here the estimated bias field and scaling factors are instead used to apply the relevant corrections to the originating DWI data. Secondly, the global intensity scaling that is calculated and applied is typically based on achieving close to a unity sum of tissue signal fractions throughout the masked region. Here, it is instead the b=0 signal in CSF that forms the reference for this global intensity scaling; this is calculated based on the estimated CSF response function and the tissue-specific intensity scaling (this is calculated internally by mtnormalise as part of its optimisation process, but typically subsequently discarded in favour of a single scaling factor for all tissues) -- **-max_iters** The maximum number of iterations. The default is 2iterations. More iterations may lead to a better mask, but for some problematic data this may lead to the masks diverging (in which case a warning is issued). +The ODFs estimated within this optimisation procedure are by default of lower maximal spherical harmonic degree than what would be advised for analysis. This is done for computational efficiency. This behaviour can be modified through the -lmax command-line option. -- **-lmax** The maximum spherical harmonic order for the output FODs. the value is passed to the dwi2fod command and should be provided in the format which it expects (the default value is "4,0,0" for multi-shell and "4,0" for single-shell data) +By default, the optimisation procedure will terminate after only two iterations. This is done because it has been observed for some data / configurations that additional iterations can lead to unstable divergence and erroneous results for bias field estimation and masking. For other configurations, it may be preferable to use a greater number of iterations, and allow the iterative algorithm to converge to a stable solution. This can be controlled via the -max_iters command-line option. + +Within the optimisation algorithm, derivation of the mask may potentially be performed differently to a conventional mask derivation that is based on a DWI series (where, in many instances, it is actually only the mean b=0 image that is used). Here, the image corresponding to the sum of tissue signal fractions following spherical deconvolution / bias field correction / intensity normalisation is also available, and this can potentially be used for mask derivation. Available options are as follows. "dwi2mask": Use the MRtrix3 command dwi2mask on the bias-field-corrected DWI series (ie. do not use the ODF tissue sum image for mask derivation); the algorithm to be invoked can be controlled by the user via the MRtrix config file entry "Dwi2maskAlgorithm". "fslbet": Invoke the FSL command "bet". "hdbet": Invoke the HD-BET command. "mrthreshold": Invoke the MRtrix3 command "mrthreshold", where an appropriate threshold value will be determined automatically (and some heuristic cleanup of the resulting mask will be performed). "synthstrip": Invoke the FreeSurfer SynthStrip method. "threshold": Apply a fixed partial volume threshold of 0.5 to the tissue ODF sum image (and some heuristic cleanup of the resulting mask will be performed). + +Options +------- Options for importing the diffusion gradient table ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -39,12 +44,27 @@ Options for importing the diffusion gradient table - **-fslgrad bvecs bvals** Provide the diffusion gradient table in FSL bvecs/bvals format -Options to specify an initial mask or tissue sum image threshold in the iterative algorithm -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +Options relevant to the internal optimisation procedure +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-dice value** Set the Dice coefficient threshold for similarity of masks between sequential iterations that will result in termination due to convergence; default = 0.999 + +- **-init_mask image** Provide an initial mask for the first iteration of the algorithm (if not provided, the default dwi2mask algorithm will be used) + +- **-max_iters count** The maximum number of iterations (see Description); default is 2; set to 0 to proceed until convergence + +- **-mask_algo algorithm** The algorithm to use for mask estimation, potentially based on the ODF sum image (see Description); default: threshold + +- **-lmax value** The maximum spherical harmonic degree for the estimated FODs (see Description); defaults are "4,0,0" for multi-shell and "4,0" for single-shell data) + +Options that modulate the outputs of the script +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-output_bias image** Export the final estimated bias field to an image -- **-threshold** Threshold on the total tissue density image used to derive the brain mask. the default is 0.5 +- **-output_tissuesum image** Export the tissue sum image that was used to generate the final mask -- **-mask_init image** Provide an initial mask to the algorithm and skip the initial masking +- **-reference value** Set the target CSF b=0 intensity in the output DWI series (default: 1000.0) Additional standard options for Python scripts ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/docs/reference/commands_list.rst b/docs/reference/commands_list.rst index 545a0da124..38b99f9ed0 100644 --- a/docs/reference/commands_list.rst +++ b/docs/reference/commands_list.rst @@ -168,7 +168,7 @@ List of MRtrix3 commands |python.png|, :ref:`dwi2response`, "Estimate response function(s) for spherical deconvolution" |cpp.png|, :ref:`dwi2tensor`, "Diffusion (kurtosis) tensor estimation" |python.png|, :ref:`dwibiascorrect`, "Perform B1 field inhomogeneity correction for a DWI volume series" - |python.png|, :ref:`dwibiasnormmask`, "Iteratively perform masking, RF estimation, CSD, bias field removal, and mask revision on a DWI series" + |python.png|, :ref:`dwibiasnormmask`, "Perform a combination of bias field correction, intensity normalisation, and mask derivation, for DWI data" |python.png|, :ref:`dwicat`, "Concatenating multiple DWI series accounting for differential intensity scaling" |cpp.png|, :ref:`dwidenoise`, "dMRI noise level estimation and denoising using Marchenko-Pastur PCA" |cpp.png|, :ref:`dwiextract`, "Extract diffusion-weighted volumes, b=0 volumes, or certain shells from a DWI dataset" diff --git a/lib/mrtrix3/app.py b/lib/mrtrix3/app.py index 3a2e843c4a..5792207121 100644 --- a/lib/mrtrix3/app.py +++ b/lib/mrtrix3/app.py @@ -330,7 +330,7 @@ def goto_scratch_dir(): #pylint: disable=unused-variable # all intermediates, the resource will be retained; if not, it will be deleted (in particular # to dynamically free up storage space used by the script). def cleanup(items): #pylint: disable=unused-variable - if not DO_CLEANUP: + if not DO_CLEANUP or not items: return if isinstance(items, list): if len(items) == 1: From a027a60e9054e75e482821c86689b9c5300785b0 Mon Sep 17 00:00:00 2001 From: Arshiyas Date: Thu, 16 Mar 2023 10:43:05 +1100 Subject: [PATCH 08/24] dwibiasnormmask: Fix string formatting error Merge of 4313a76 and 68c532a. --- bin/dwibiasnormmask | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bin/dwibiasnormmask b/bin/dwibiasnormmask index bf0d755361..24682d09ed 100755 --- a/bin/dwibiasnormmask +++ b/bin/dwibiasnormmask @@ -364,7 +364,7 @@ def execute(): #pylint: disable=unused-variable if new_dice_coefficient > app.ARGS.dice: progress.done() - app.console('Exiting loop after %d iterations due to mask convergence (Dice coefficient = %f)' % iteration, new_dice_coefficient) + app.console('Exiting loop after %d iterations due to mask convergence (Dice coefficient = %f)' % (iteration, new_dice_coefficient)) app.cleanup(old_dwi_image) app.cleanup(old_dwi_mask_image) app.cleanup(old_tissue_sum_image) From 4803f031e80d761f7d27cdeb08d9378884db68ab Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Sat, 25 Feb 2023 23:38:18 +1100 Subject: [PATCH 09/24] dwibiasnormmask: New option -output_scale Exports total global scaling multiplier that was applied for intensity normalisation. Also includes some pylint fixes. --- bin/dwibiasnormmask | 37 +++++++++++++++------ docs/reference/commands/dwibiasnormmask.rst | 4 ++- 2 files changed, 29 insertions(+), 12 deletions(-) diff --git a/bin/dwibiasnormmask b/bin/dwibiasnormmask index 24682d09ed..df4d4b0784 100755 --- a/bin/dwibiasnormmask +++ b/bin/dwibiasnormmask @@ -76,9 +76,14 @@ def usage(cmdline): #pylint: disable=unused-variable cmdline.add_argument('output_dwi', help='The output corrected DWI series') cmdline.add_argument('output_mask', help='The output DWI mask') output_options = cmdline.add_argument_group('Options that modulate the outputs of the script') - output_options.add_argument('-output_bias', metavar='image', help='Export the final estimated bias field to an image') - output_options.add_argument('-output_tissuesum', metavar='image', help='Export the tissue sum image that was used to generate the final mask') - output_options.add_argument('-reference', type=float, metavar='value', default=REFERENCE_INTENSITY, help='Set the target CSF b=0 intensity in the output DWI series (default: ' + str(REFERENCE_INTENSITY) + ')') + output_options.add_argument('-output_bias', metavar='image', + help='Export the final estimated bias field to an image') + output_options.add_argument('-output_scale', metavar='file', + help='Write the scaling factor applied to the DWI series to a text file') + output_options.add_argument('-output_tissuesum', metavar='image', + help='Export the tissue sum image that was used to generate the final mask') + output_options.add_argument('-reference', type=float, metavar='value', default=REFERENCE_INTENSITY, + help='Set the target CSF b=0 intensity in the output DWI series (default: ' + str(REFERENCE_INTENSITY) + ')') internal_options = cmdline.add_argument_group('Options relevant to the internal optimisation procedure') internal_options.add_argument('-dice', type=float, default=DICE_COEFF_DEFAULT, metavar='value', help='Set the Dice coefficient threshold for similarity of masks between sequential iterations that will ' @@ -87,13 +92,13 @@ def usage(cmdline): #pylint: disable=unused-variable help='Provide an initial mask for the first iteration of the algorithm ' '(if not provided, the default dwi2mask algorithm will be used)') internal_options.add_argument('-max_iters', type=int, default=DWIBIASCORRECT_MAX_ITERS, metavar='count', - help='The maximum number of iterations (see Description); default is ' + str(DWIBIASCORRECT_MAX_ITERS) + '; ' - 'set to 0 to proceed until convergence') + help='The maximum number of iterations (see Description); default is ' + str(DWIBIASCORRECT_MAX_ITERS) + '; ' + 'set to 0 to proceed until convergence') internal_options.add_argument('-mask_algo', choices=MASK_ALGOS, default=MASK_ALGO_DEFAULT, metavar='algorithm', - help='The algorithm to use for mask estimation, potentially based on the ODF sum image (see Description); default: ' + MASK_ALGO_DEFAULT) - internal_options.add_argument('-lmax', metavar='value', - help='The maximum spherical harmonic degree for the estimated FODs (see Description); ' - 'defaults are "' + ','.join(str(item) for item in LMAXES_MULTI) + '" for multi-shell and "' + ','.join(str(item) for item in LMAXES_SINGLE) + '" for single-shell data)') + help='The algorithm to use for mask estimation, potentially based on the ODF sum image (see Description); default: ' + MASK_ALGO_DEFAULT) + internal_options.add_argument('-lmax', metavar='values', + help='The maximum spherical harmonic degree for the estimated FODs (see Description); ' + 'defaults are "' + ','.join(str(item) for item in LMAXES_MULTI) + '" for multi-shell and "' + ','.join(str(item) for item in LMAXES_SINGLE) + '" for single-shell data)') app.add_dwgrad_import_options(cmdline) @@ -109,14 +114,16 @@ def execute(): #pylint: disable=unused-variable if app.ARGS.lmax: try: lmax = [int(i) for i in app.ARGS.lmax.split(',')] - except ValueError: - raise MRtrixError('Values provided to -lmax option must be a comma-separated list of integers') + except ValueError as exc: + raise MRtrixError('Values provided to -lmax option must be a comma-separated list of integers') from exc if any(value < 0 or value % 2 for value in lmax): raise MRtrixError('lmax values must be non-negative even integers') if len(lmax) not in [2, 3]: raise MRtrixError('Length of lmax vector expected to be either 2 or 3') if app.ARGS.dice <= 0.0 or app.ARGS.dice > 1.0: raise MRtrixError('Dice coefficient for convergence detection must lie in the range (0.0, 1.0]') + if app.ARGS.reference <= 0.0: + raise MRtrixError('Reference intensity must be positive') # Check mask algorithm, including availability of external software if necessary for mask_algo, software, command in [('fslbet', 'FSL', 'bet'), ('hdbet', 'HD-BET', 'hd-bet'), ('synthstrip', 'FreeSurfer', 'mri_synthstrip')]: @@ -196,6 +203,7 @@ def execute(): #pylint: disable=unused-variable iteration = 0 step = 'initialisation' prev_dice_coefficient = 0.0 + total_scaling_factor = 1.0 def msg(): return 'Iteration {0}; {1} step; previous Dice coefficient {2}' \ @@ -368,6 +376,7 @@ def execute(): #pylint: disable=unused-variable app.cleanup(old_dwi_image) app.cleanup(old_dwi_mask_image) app.cleanup(old_tissue_sum_image) + total_scaling_factor *= scale_multiplier break if new_dice_coefficient < prev_dice_coefficient: @@ -388,11 +397,13 @@ def execute(): #pylint: disable=unused-variable app.cleanup(old_dwi_image) app.cleanup(old_dwi_mask_image) app.cleanup(old_tissue_sum_image) + total_scaling_factor *= scale_multiplier break iteration += 1 prev_dice_coefficient = new_dice_coefficient + run.command(['mrconvert', dwi_image, path.from_user(app.ARGS.output_dwi, False)], mrconvert_keyval=path.from_user(app.ARGS.input, False), force=app.FORCE_OVERWRITE) @@ -402,12 +413,16 @@ def execute(): #pylint: disable=unused-variable mrconvert_keyval=path.from_user(app.ARGS.input, False), force=app.FORCE_OVERWRITE) + if app.ARGS.output_scale: + matrix.save_vector(path.from_user(app.ARGS.output_scale, False), [total_scaling_factor]) + if app.ARGS.output_tissuesum: run.command(['mrconvert', tissue_sum_image, path.from_user(app.ARGS.output_tissuesum, False)], mrconvert_keyval=path.from_user(app.ARGS.input, False), force=app.FORCE_OVERWRITE) + # Execute the script import mrtrix3 #pylint: disable=wrong-import-position mrtrix3.execute() #pylint: disable=no-member diff --git a/docs/reference/commands/dwibiasnormmask.rst b/docs/reference/commands/dwibiasnormmask.rst index 23da1312db..56d66e83b6 100644 --- a/docs/reference/commands/dwibiasnormmask.rst +++ b/docs/reference/commands/dwibiasnormmask.rst @@ -55,13 +55,15 @@ Options relevant to the internal optimisation procedure - **-mask_algo algorithm** The algorithm to use for mask estimation, potentially based on the ODF sum image (see Description); default: threshold -- **-lmax value** The maximum spherical harmonic degree for the estimated FODs (see Description); defaults are "4,0,0" for multi-shell and "4,0" for single-shell data) +- **-lmax values** The maximum spherical harmonic degree for the estimated FODs (see Description); defaults are "4,0,0" for multi-shell and "4,0" for single-shell data) Options that modulate the outputs of the script ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - **-output_bias image** Export the final estimated bias field to an image +- **-output_scale file** Write the scaling factor applied to the DWI series to a text file + - **-output_tissuesum image** Export the tissue sum image that was used to generate the final mask - **-reference value** Set the target CSF b=0 intensity in the output DWI series (default: 1000.0) From df6732f4e5db42ad265fe85e808138ddab6ced13 Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Sat, 25 Feb 2023 23:41:32 +1100 Subject: [PATCH 10/24] Initial revisions of "mtnorm" Python algorithms --- docs/reference/commands/dwi2mask.rst | 16 ++- docs/reference/commands/dwibiascorrect.rst | 10 +- docs/reference/commands/dwinormalise.rst | 21 +-- lib/mrtrix3/dwi2mask/mtnorm.py | 142 ++++++++++++-------- lib/mrtrix3/dwibiascorrect/mtnorm.py | 105 ++++++++------- lib/mrtrix3/dwinormalise/mtnorm.py | 146 ++++++++++++--------- 6 files changed, 261 insertions(+), 179 deletions(-) diff --git a/docs/reference/commands/dwi2mask.rst b/docs/reference/commands/dwi2mask.rst index 05f2c774d8..b95bd071ef 100644 --- a/docs/reference/commands/dwi2mask.rst +++ b/docs/reference/commands/dwi2mask.rst @@ -974,7 +974,7 @@ dwi2mask mtnorm Synopsis -------- -Derives a DWI brain mask by calculating and then thresholding a tissue density image +Derives a DWI brain mask by calculating and then thresholding a sum-of-tissue-densities image Usage ----- @@ -989,7 +989,11 @@ Usage Description ----------- -This script first derives an initial brain mask using the legacy MRtrix3 dwi2mask heuristic (based on thresholded trace images), and then performs response function estimation and multi-tissue CSD (with a lower lmax than the dwi2fod default, for higher speed), before summing the derived tissue ODFs and thresholding the resulting tissue density image to derive a DWI brain mask. +This script attempts to identify brain voxels based on the total density of macroscopic tissues as estimated through multi-tissue deconvolution. Following response function estimation and multi-tissue deconvolution, the sum of tissue densities is thresholded at a fixed value (default is 0.5), and subsequent mask image cleaning operations are performed. + +The operation of this script is a subset of that performed by the script "dwibiasnormmask". Many users may find that comprehensive solution preferable; this dwi2mask algorithm is nevertheless provided to demonstrate specifically the mask estimation portion of that command. + +The ODFs estimated within this optimisation procedure are by default of lower maximal spherical harmonic degree than what would be advised for analysis. This is done for computational efficiency. This behaviour can be modified through the -lmax command-line option. Options ------- @@ -997,9 +1001,13 @@ Options Options specific to the "mtnorm" algorithm ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -- **-threshold** the threshold on the total tissue density sum image used to derive the brain mask. the default is 0.5 +- **-init_mask image** Provide an initial brain mask, which will constrain the response function estimation (if omitted, the default dwi2mask algorithm will be used) + +- **-lmax values** The maximum spherical harmonic degree for the estimated FODs (see Description); defaults are "4,0,0" for multi-shell and "4,0" for single-shell data) + +- **-threshold value** the threshold on the total tissue density sum image used to derive the brain mask; default is 0.5 -- **-lmax** the maximum spherical harmonic order for the output FODs. The value is passed to the dwi2fod command and should be provided in the format which it expects (Default is "4,0,0" for multi-shell and "4,0" for single-shell data) +- **-tissuesum image** Export the tissue sum image that was used to generate the mask Options for importing the diffusion gradient table ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/docs/reference/commands/dwibiascorrect.rst b/docs/reference/commands/dwibiascorrect.rst index 7d8444ed53..667c675404 100644 --- a/docs/reference/commands/dwibiascorrect.rst +++ b/docs/reference/commands/dwibiascorrect.rst @@ -304,7 +304,7 @@ dwibiascorrect mtnorm Synopsis -------- -Performs DWI bias field correction +Perform DWI bias field correction using the "mtnormalise" command Usage ----- @@ -319,7 +319,11 @@ Usage Description ----------- -This script first performs response function estimation and multi-tissue CSD within a provided or derived brain mask (with a lower lmax than the dwi2fod default, for higher speed) before calling mtnormalise and utilizing the final estimated spatially varying intensity level used for normalisation to correct bias fields +This algorithm bases its operation almost entirely on the utilisation of multi-tissue decomposition information to estimate an underlying B1 receive field, as is implemented in the MRtrix3 command "mtnormalise". Its typical usage is however slightly different, in that the primary output of the command is not the bias-field-corrected FODs, but a bias-field-corrected version of the DWI series. + +The operation of this script is a subset of that performed by the script "dwibiasnormmask". Many users may find that comprehensive solution preferable; this dwibiascorrect algorithm is nevertheless provided to demonstrate specifically the bias field correction portion of that command. + +The ODFs estimated within this optimisation procedure are by default of lower maximal spherical harmonic degree than what would be advised for analysis. This is done for computational efficiency. This behaviour can be modified through the -lmax command-line option. Options ------- @@ -327,7 +331,7 @@ Options Options specific to the "mtnorm" algorithm ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -- **-lmax** the maximum spherical harmonic order for the output FODs. The value is passed to the dwi2fod command and should be provided in the format which it expects (Default is "4,0,0" for multi-shell and "4,0" for single-shell data) +- **-lmax values** The maximum spherical harmonic degree for the estimated FODs (see Description); defaults are "4,0,0" for multi-shell and "4,0" for single-shell data) Options for importing the diffusion gradient table ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/docs/reference/commands/dwinormalise.rst b/docs/reference/commands/dwinormalise.rst index 0039a4ef91..bf7c3ff2d3 100644 --- a/docs/reference/commands/dwinormalise.rst +++ b/docs/reference/commands/dwinormalise.rst @@ -265,23 +265,26 @@ dwinormalise mtnorm Synopsis -------- -Performs a global DWI intensity normalisation, normalising mean b=0 CSF signal intensity +Normalise a DWI series to the estimated b=0 CSF intensity Usage ----- :: - dwinormalise mtnorm input_dwi output [ options ] + dwinormalise mtnorm input output [ options ] -- *input_dwi*: The input DWI series +- *input*: The input DWI series - *output*: The normalised DWI series Description ----------- -This script first performs response function estimation and multi-tissue CSD within a provided brain mask(with a lower lmax than the dwi2fod default, for higher speed) before using the global scaling factor derived by mtnormalise to scale the intensity of the DWI series such that the mean b=0 CSF signal intensity is approximately close to some reference value (by default 1000).Note that if the -mask command-line option is not specified, the MRtrix3 command dwi2mask will automatically be called to derive a mask that will be passed to the relevant bias field estimation command. More information on mask derivation from DWI data can be found at the following link: -https://mrtrix.readthedocs.io/en/3.0.4/dwi_preprocessing/masking.html +This algorithm determines an appropriate global scaling factor to apply to a DWI series such that after the scaling is applied, the b=0 CSF intensity corresponds to some reference value (1000 by default). + +The operation of this script is a subset of that performed by the script "dwibiasnormmask". Many users may find that comprehensive solution preferable; this dwinormalise algorithm is nevertheless provided to demonstrate specifically the global intensituy normalisation portion of that command. + +The ODFs estimated within this optimisation procedure are by default of lower maximal spherical harmonic degree than what would be advised for analysis. This is done for computational efficiency. This behaviour can be modified through the -lmax command-line option. Options ------- @@ -296,11 +299,13 @@ Options for importing the diffusion gradient table Options specific to the "mtnorm" algorithm ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -- **-target** the threshold on the balanced tissue sum image used to derive the brain mask. the default is 0.5 +- **-lmax values** The maximum spherical harmonic degree for the estimated FODs (see Description); defaults are "4,0,0" for multi-shell and "4,0" for single-shell data) + +- **-mask image** Provide a mask image for relevant calculations (if not provided, the default dwi2mask algorithm will be used) -- **-lmax** the maximum spherical harmonic order for the output FODs. The value is passed to the dwi2fod command and should be provided in the format which it expects (Default is "4,0,0" for multi-shell and "4,0" for single-shell data) +- **-reference value** Set the target CSF b=0 intensity in the output DWI series (default: 1000) -- **-mask image** Manually provide a mask image for normalisation +- **-scale file** Write the scaling factor applied to the DWI series to a text file Additional standard options for Python scripts ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/lib/mrtrix3/dwi2mask/mtnorm.py b/lib/mrtrix3/dwi2mask/mtnorm.py index 014b679ff9..3c2452e3cb 100644 --- a/lib/mrtrix3/dwi2mask/mtnorm.py +++ b/lib/mrtrix3/dwi2mask/mtnorm.py @@ -14,34 +14,54 @@ # For more details, see http://www.mrtrix.org/. import math -from mrtrix3 import app, run, image +from mrtrix3 import MRtrixError +from mrtrix3 import app, image, path, run + + +LMAXES_MULTI = [4, 0, 0] +LMAXES_SINGLE = [4, 0] +THRESHOLD_DEFAULT = 0.5 + -LMAXES_MULTI = '4,0,0' -LMAXES_SINGLE = '4,0' def usage(base_parser, subparsers): #pylint: disable=unused-variable parser = subparsers.add_parser('mtnorm', parents=[base_parser]) parser.set_author('Robert E. Smith (robert.smith@florey.edu.au) and Arshiya Sangchooli (asangchooli@student.unimelb.edu.au)') - parser.set_synopsis('Derives a DWI brain mask by calculating and then thresholding a tissue density image') - parser.add_description('This script first derives an initial brain mask using the legacy MRtrix3 dwi2mask heuristic ' - '(based on thresholded trace images), and then performs response function estimation and multi-tissue CSD ' - '(with a lower lmax than the dwi2fod default, for higher speed), before summing the derived tissue ODFs and ' - 'thresholding the resulting tissue density image to derive a DWI brain mask.') - parser.add_argument('input', help='The input DWI series') + parser.set_synopsis('Derives a DWI brain mask by calculating and then thresholding a sum-of-tissue-densities image') + parser.add_description('This script attempts to identify brain voxels based on the total density of macroscopic ' + 'tissues as estimated through multi-tissue deconvolution. Following response function ' + 'estimation and multi-tissue deconvolution, the sum of tissue densities is thresholded at a ' + 'fixed value (default is ' + str(THRESHOLD_DEFAULT) + '), and subsequent mask image cleaning ' + 'operations are performed.') + parser.add_description('The operation of this script is a subset of that performed by the script "dwibiasnormmask". ' + 'Many users may find that comprehensive solution preferable; this dwi2mask algorithm is nevertheless ' + 'provided to demonstrate specifically the mask estimation portion of that command.') + parser.add_description('The ODFs estimated within this optimisation procedure are by default of lower maximal spherical harmonic ' + 'degree than what would be advised for analysis. This is done for computational efficiency. This ' + 'behaviour can be modified through the -lmax command-line option.') + parser.add_argument('input', help='The input DWI series') parser.add_argument('output', help='The output mask image') options = parser.add_argument_group('Options specific to the "mtnorm" algorithm') + options.add_argument('-init_mask', + metavar='image', + help='Provide an initial brain mask, which will constrain the response function estimation ' + '(if omitted, the default dwi2mask algorithm will be used)') + options.add_argument('-lmax', + metavar='values', + help='The maximum spherical harmonic degree for the estimated FODs (see Description); ' + 'defaults are "' + ','.join(str(item) for item in LMAXES_MULTI) + '" for multi-shell and "' + ','.join(str(item) for item in LMAXES_SINGLE) + '" for single-shell data)') options.add_argument('-threshold', type=float, - default=0.5, - help='the threshold on the total tissue density sum image used to derive the brain mask. the default is 0.5') - options.add_argument('-lmax', - type=str, - help='the maximum spherical harmonic order for the output FODs. The value is passed to ' - 'the dwi2fod command and should be provided in the format which it expects ' - '(Default is "' + str(LMAXES_MULTI) + '" for multi-shell and "' + str(LMAXES_SINGLE) + '" for single-shell data)') + metavar='value', + default=THRESHOLD_DEFAULT, + help='the threshold on the total tissue density sum image used to derive the brain mask; default is ' + str(THRESHOLD_DEFAULT)) + options.add_argument('-tissuesum', metavar='image', help='Export the tissue sum image that was used to generate the mask') + + def get_inputs(): #pylint: disable=unused-variable - pass + if app.ARGS.init_mask: + run.command(['mrconvert', path.from_user(app.ARGS.init_mask, False), path.to_scratch('init_mask.mif', False), '-datatype', 'bit']) def needs_mean_bzero(): #pylint: disable=unused-variable @@ -49,19 +69,32 @@ def needs_mean_bzero(): #pylint: disable=unused-variable def execute(): #pylint: disable=unused-variable + + # Verify user inputs + lmax = None + if app.ARGS.lmax: + try: + lmax = [int(i) for i in app.ARGS.lmax.split(',')] + except ValueError as exc: + raise MRtrixError('Values provided to -lmax option must be a comma-separated list of integers') from exc + if any(value < 0 or value % 2 for value in lmax): + raise MRtrixError('lmax values must be non-negative even integers') + if len(lmax) not in [2, 3]: + raise MRtrixError('Length of lmax vector expected to be either 2 or 3') + if app.ARGS.threshold <= 0.0 or app.ARGS.threshold >= 1.0: + raise MRtrixError('Tissue density sum threshold must lie within the range (0.0, 1.0)') + # Determine whether we are working with single-shell or multi-shell data bvalues = [ int(round(float(value))) for value in image.mrinfo('input.mif', 'shell_bvalues') \ .strip().split()] multishell = (len(bvalues) > 2) + if lmax is None: + lmax = LMAXES_MULTI if multishell else LMAXES_SINGLE + elif len(lmax) == 3 and not multishell: + raise MRtrixError('User specified 3 lmax values for three-tissue decomposition, but input DWI is not multi-shell') - # Step 1: Initial DWI brain mask - dwi_mask_image = 'dwi_mask_init.mif' - app.debug('Performing intial DWI brain masking with the legacy MRtrix3 algorithm') - run.command('dwi2mask legacy input.mif ' + dwi_mask_image) - - # Step 2: RF estimation / CSD / mtnormalise / mask revision class Tissue(object): #pylint: disable=useless-object-inheritance def __init__(self, name): self.name = name @@ -71,54 +104,51 @@ def __init__(self, name): dwi_image = 'input.mif' tissues = [Tissue('WM'), Tissue('GM'), Tissue('CSF')] - app.debug('Estimating response function using initial brain mask...') run.command('dwi2response dhollander ' + dwi_image - + ' -mask ' - + dwi_mask_image + + (' -mask init_mask.mif' if app.ARGS.init_mask else '') + ' ' + ' '.join(tissue.tissue_rf for tissue in tissues)) - # Remove GM if we can't deal with it - if app.ARGS.lmax: - lmaxes = app.ARGS.lmax - else: - lmaxes = LMAXES_MULTI - if not multishell: - app.cleanup(tissues[1].tissue_rf) - tissues = tissues[::2] - lmaxes = LMAXES_SINGLE - - app.debug('FOD estimation with lmaxes ' + lmaxes + '...') + # Immediately remove GM if we can't deal with it + if not multishell: + app.cleanup(tissues[1].tissue_rf) + tissues = tissues[::2] + run.command('dwi2fod msmt_csd ' + dwi_image - + ' -lmax ' + lmaxes + + ' -lmax ' + ','.join(str(item) for item in lmax) + ' ' + ' '.join(tissue.tissue_rf + ' ' + tissue.fod - for tissue in tissues)) - - app.debug('Deriving final brain mask by thresholding tissue sum image...') - new_dwi_mask_image = 'dwi_mask' + '.mif' - - tissuesum_threshold = app.ARGS.threshold / math.sqrt(4.0 * math.pi) + for tissue in tissues)) + tissue_sum_image = 'tissuesum.mif' run.command('mrconvert ' - + tissues[0].fod - + ' -coord 3 0 - |' - + ' mrmath - ' - + ' '.join(tissue.fod for tissue in tissues[1:]) - + ' sum - |' - + ' mrthreshold - -abs ' - + str(tissuesum_threshold) + + tissues[0].fod + + ' -coord 3 0 - |' + + ' mrmath - ' + + ' '.join(tissue.fod for tissue in tissues[1:]) + + ' sum - | ' + + 'mrcalc - ' + str(math.sqrt(4.0 * math.pi)) + ' -mult ' + + tissue_sum_image) + + mask_image = 'mask.mif' + run.command('mrthreshold ' + + tissue_sum_image + + ' -abs ' + + str(app.ARGS.threshold) + ' - |' + ' maskfilter - connect -largest - |' + ' mrcalc 1 - -sub - -datatype bit |' + ' maskfilter - connect -largest - |' + ' mrcalc 1 - -sub - -datatype bit |' - + ' maskfilter - clean - |' - + ' mrcalc - input_pos_mask.mif -mult ' - + new_dwi_mask_image - + ' -datatype bit') + + ' maskfilter - clean ' + + mask_image) app.cleanup([tissue.fod for tissue in tissues]) - return dwi_mask_image + if app.ARGS.tissuesum: + run.command(['mrconvert', tissue_sum_image, path.from_user(app.ARGS.tissuesum, False)], + mrconvert_keyval=path.from_user(app.ARGS.input, False), + force=app.FORCE_OVERWRITE) + + return mask_image diff --git a/lib/mrtrix3/dwibiascorrect/mtnorm.py b/lib/mrtrix3/dwibiascorrect/mtnorm.py index 22d28b3af8..0d54872ad2 100644 --- a/lib/mrtrix3/dwibiascorrect/mtnorm.py +++ b/lib/mrtrix3/dwibiascorrect/mtnorm.py @@ -13,42 +13,72 @@ # # For more details, see http://www.mrtrix.org/. -from mrtrix3 import app, run, image, path +from mrtrix3 import MRtrixError +from mrtrix3 import app, image, path, run -LMAXES_MULTI = '4,0,0' -LMAXES_SINGLE = '4,0' + +LMAXES_MULTI = [4, 0, 0] +LMAXES_SINGLE = [4, 0] def usage(base_parser, subparsers): #pylint: disable=unused-variable parser = subparsers.add_parser('mtnorm', parents=[base_parser]) parser.set_author('Robert E. Smith (robert.smith@florey.edu.au) and Arshiya Sangchooli (asangchooli@student.unimelb.edu.au)') - parser.set_synopsis('Performs DWI bias field correction') - parser.add_description('This script first performs response function estimation and multi-tissue CSD within a provided or derived ' - 'brain mask (with a lower lmax than the dwi2fod default, for higher speed) before calling mtnormalise and ' - 'utilizing the final estimated spatially varying intensity level used for normalisation to correct bias fields') - parser.add_argument('input', help='The input image series to be corrected') + parser.set_synopsis('Perform DWI bias field correction using the "mtnormalise" command') + parser.add_description('This algorithm bases its operation almost entirely on the utilisation of multi-tissue ' + 'decomposition information to estimate an underlying B1 receive field, as is implemented ' + 'in the MRtrix3 command "mtnormalise". Its typical usage is however slightly different, ' + 'in that the primary output of the command is not the bias-field-corrected FODs, but a ' + 'bias-field-corrected version of the DWI series.') + parser.add_description('The operation of this script is a subset of that performed by the script "dwibiasnormmask". ' + 'Many users may find that comprehensive solution preferable; this dwibiascorrect algorithm is ' + 'nevertheless provided to demonstrate specifically the bias field correction portion of that command.') + parser.add_description('The ODFs estimated within this optimisation procedure are by default of lower maximal spherical harmonic ' + 'degree than what would be advised for analysis. This is done for computational efficiency. This ' + 'behaviour can be modified through the -lmax command-line option.') + parser.add_argument('input', help='The input image series to be corrected') parser.add_argument('output', help='The output corrected image series') options = parser.add_argument_group('Options specific to the "mtnorm" algorithm') options.add_argument('-lmax', - type=str, - help='the maximum spherical harmonic order for the output FODs. The value is passed to ' - 'the dwi2fod command and should be provided in the format which it expects ' - '(Default is "' + str(LMAXES_MULTI) + '" for multi-shell and "' + str(LMAXES_SINGLE) + '" for single-shell data)') + metavar='values', + help='The maximum spherical harmonic degree for the estimated FODs (see Description); ' + 'defaults are "' + ','.join(str(item) for item in LMAXES_MULTI) + '" for multi-shell and "' + ','.join(str(item) for item in LMAXES_SINGLE) + '" for single-shell data)') + def check_output_paths(): #pylint: disable=unused-variable pass + + def get_inputs(): #pylint: disable=unused-variable pass + def execute(): #pylint: disable=unused-variable + + # Verify user inputs + lmax = None + if app.ARGS.lmax: + try: + lmax = [int(i) for i in app.ARGS.lmax.split(',')] + except ValueError as exc: + raise MRtrixError('Values provided to -lmax option must be a comma-separated list of integers') from exc + if any(value < 0 or value % 2 for value in lmax): + raise MRtrixError('lmax values must be non-negative even integers') + if len(lmax) not in [2, 3]: + raise MRtrixError('Length of lmax vector expected to be either 2 or 3') + # Determine whether we are working with single-shell or multi-shell data bvalues = [ int(round(float(value))) for value in image.mrinfo('in.mif', 'shell_bvalues') \ .strip().split()] multishell = (len(bvalues) > 2) + if lmax is None: + lmax = LMAXES_MULTI if multishell else LMAXES_SINGLE + elif len(lmax) == 3 and not multishell: + raise MRtrixError('User specified 3 lmax values for three-tissue decomposition, but input DWI is not multi-shell') # RF estimation and multi-tissue CSD class Tissue(object): #pylint: disable=useless-object-inheritance @@ -60,35 +90,24 @@ def __init__(self, name): tissues = [Tissue('WM'), Tissue('GM'), Tissue('CSF')] - app.debug('Estimating response function using brain mask...') - run.command('dwi2response dhollander ' - + 'in.mif ' - + '-mask mask.mif ' + run.command('dwi2response dhollander in.mif' + + (' -mask mask.mif' if app.ARGS.mask else '') + + ' ' + ' '.join(tissue.tissue_rf for tissue in tissues)) - # Remove GM if we can't deal with it - if app.ARGS.lmax: - lmaxes = app.ARGS.lmax - else: - lmaxes = LMAXES_MULTI - if not multishell: - app.cleanup(tissues[1].tissue_rf) - tissues = tissues[::2] - lmaxes = LMAXES_SINGLE - - app.debug('FOD estimation with lmaxes ' + lmaxes + '...') - run.command('dwi2fod msmt_csd ' - + 'in.mif' - + ' -lmax ' + lmaxes + # Immediately remove GM if we can't deal with it + if not multishell: + app.cleanup(tissues[1].tissue_rf) + tissues = tissues[::2] + + run.command('dwi2fod msmt_csd in.mif' + + ' -lmax ' + ','.join(str(item) for item in lmax) + ' ' + ' '.join(tissue.tissue_rf + ' ' + tissue.fod - for tissue in tissues)) + for tissue in tissues)) - app.debug('correcting bias field...') - run.command('maskfilter' - + ' mask.mif' - + ' erode - |' - + ' mtnormalise -mask - -balanced' + run.command('maskfilter mask.mif erode - | ' + + 'mtnormalise -mask - -balanced' + ' -check_norm field.mif ' + ' '.join(tissue.fod + ' ' + tissue.fod_norm for tissue in tissues)) @@ -96,14 +115,12 @@ def __init__(self, name): app.cleanup([tissue.fod_norm for tissue in tissues]) app.cleanup([tissue.tissue_rf for tissue in tissues]) - run.command('mrcalc ' - + 'in.mif ' - + 'field.mif ' - + '-div ' - + path.from_user(app.ARGS.output), + run.command('mrcalc in.mif field.mif -div - | ' + 'mrconvert - '+ path.from_user(app.ARGS.output), + mrconvert_keyval=path.from_user(app.ARGS.input, False), force=app.FORCE_OVERWRITE) if app.ARGS.bias: - run.command('mrconvert field.mif ' - + path.from_user(app.ARGS.bias), - mrconvert_keyval=path.from_user(app.ARGS.input, False), force=app.FORCE_OVERWRITE) + run.command('mrconvert field.mif ' + path.from_user(app.ARGS.bias), + mrconvert_keyval=path.from_user(app.ARGS.input, False), + force=app.FORCE_OVERWRITE) diff --git a/lib/mrtrix3/dwinormalise/mtnorm.py b/lib/mrtrix3/dwinormalise/mtnorm.py index 4526e16888..d472c62729 100644 --- a/lib/mrtrix3/dwinormalise/mtnorm.py +++ b/lib/mrtrix3/dwinormalise/mtnorm.py @@ -15,66 +15,91 @@ import math from mrtrix3 import CONFIG, MRtrixError -from mrtrix3 import app, run, image, path, matrix, _version +from mrtrix3 import app, image, matrix, path, run -DEFAULT_TARGET_INTENSITY=1000 -LMAXES_MULTI = '4,0,0' -LMAXES_SINGLE = '4,0' +REFERENCE_INTENSITY = 1000 + +LMAXES_MULTI = [4, 0, 0] +LMAXES_SINGLE = [4, 0] def usage(base_parser, subparsers): #pylint: disable=unused-variable parser = subparsers.add_parser('mtnorm', parents=[base_parser]) parser.set_author('Robert E. Smith (robert.smith@florey.edu.au) and Arshiya Sangchooli (asangchooli@student.unimelb.edu.au)') - parser.set_synopsis('Performs a global DWI intensity normalisation, normalising mean b=0 CSF signal intensity') - parser.add_description('This script first performs response function estimation and multi-tissue CSD within a provided brain mask' - '(with a lower lmax than the dwi2fod default, for higher speed) before using the global scaling factor derived ' - 'by mtnormalise to scale the intensity of the DWI series such that the mean b=0 CSF signal intensity is ' - 'approximately close to some reference value (by default ' + str(DEFAULT_TARGET_INTENSITY) + ').' - 'Note that if the -mask command-line option is not specified, the MRtrix3 command dwi2mask will automatically be ' - 'called to derive a mask that will be passed to the relevant bias field estimation command. ' - 'More information on mask derivation from DWI data can be found at the following link: \n' - 'https://mrtrix.readthedocs.io/en/' + _version.__tag__ + '/dwi_preprocessing/masking.html') - parser.add_argument('input_dwi', help='The input DWI series') + parser.set_synopsis('Normalise a DWI series to the estimated b=0 CSF intensity') + parser.add_description('This algorithm determines an appropriate global scaling factor to apply to a DWI series ' + 'such that after the scaling is applied, the b=0 CSF intensity corresponds to some ' + 'reference value (' + str(REFERENCE_INTENSITY) + ' by default).') + parser.add_description('The operation of this script is a subset of that performed by the script "dwibiasnormmask". ' + 'Many users may find that comprehensive solution preferable; this dwinormalise algorithm is ' + 'nevertheless provided to demonstrate specifically the global intensituy normalisation portion of that command.') + parser.add_description('The ODFs estimated within this optimisation procedure are by default of lower maximal spherical harmonic ' + 'degree than what would be advised for analysis. This is done for computational efficiency. This ' + 'behaviour can be modified through the -lmax command-line option.') + parser.add_argument('input', help='The input DWI series') parser.add_argument('output', help='The normalised DWI series') options = parser.add_argument_group('Options specific to the "mtnorm" algorithm') - options.add_argument('-target', type=float, default=DEFAULT_TARGET_INTENSITY, - help='the threshold on the balanced tissue sum image used to derive the brain mask. the default is 0.5') - options.add_argument('-lmax', type=str, - help='the maximum spherical harmonic order for the output FODs. The value is passed to ' - 'the dwi2fod command and should be provided in the format which it expects ' - '(Default is "' + str(LMAXES_MULTI) + '" for multi-shell and "' + str(LMAXES_SINGLE) + '" for single-shell data)') - options.add_argument('-mask', metavar='image', help='Manually provide a mask image for normalisation') + options.add_argument('-lmax', + metavar='values', + help='The maximum spherical harmonic degree for the estimated FODs (see Description); ' + 'defaults are "' + ','.join(str(item) for item in LMAXES_MULTI) + '" for multi-shell and "' + ','.join(str(item) for item in LMAXES_SINGLE) + '" for single-shell data)') + options.add_argument('-mask', + metavar='image', + help='Provide a mask image for relevant calculations ' + '(if not provided, the default dwi2mask algorithm will be used)') + options.add_argument('-reference', + type=float, + metavar='value', + default=REFERENCE_INTENSITY, + help='Set the target CSF b=0 intensity in the output DWI series ' + '(default: ' + str(REFERENCE_INTENSITY) + ')') + options.add_argument('-scale', + metavar='file', + help='Write the scaling factor applied to the DWI series to a text file') app.add_dwgrad_import_options(parser) + + def check_output_paths(): #pylint: disable=unused-variable app.check_output_path(app.ARGS.output) + def execute(): #pylint: disable=unused-variable - # Get input data into the scratch directory - grad_option = '' - if app.ARGS.grad: - grad_option = ' -grad ' + path.from_user(app.ARGS.grad) - elif app.ARGS.fslgrad: - grad_option = ' -fslgrad ' + path.from_user(app.ARGS.fslgrad[0]) + ' ' + path.from_user(app.ARGS.fslgrad[1]) - app.make_scratch_dir() - app.goto_scratch_dir() + # Verify user inputs + lmax = None + if app.ARGS.lmax: + try: + lmax = [int(i) for i in app.ARGS.lmax.split(',')] + except ValueError as exc: + raise MRtrixError('Values provided to -lmax option must be a comma-separated list of integers') from exc + if any(value < 0 or value % 2 for value in lmax): + raise MRtrixError('lmax values must be non-negative even integers') + if len(lmax) not in [2, 3]: + raise MRtrixError('Length of lmax vector expected to be either 2 or 3') + if app.ARGS.reference <= 0.0: + raise MRtrixError('Reference intensity must be positive') + + grad_option = app.read_dwgrad_import_options() + # Get input data into the scratch directory + app.make_scratch_dir() run.command('mrconvert ' - + path.from_user(app.ARGS.input_dwi) + + path.from_user(app.ARGS.input) + ' ' + path.to_scratch('input.mif') + grad_option) + if app.ARGS.mask: + run.command('mrconvert ' + path.from_user(app.ARGS.mask) + ' ' + path.to_scratch('mask.mif') + ' -datatype bit') + app.goto_scratch_dir() - # Get mask into scratch directory + # Make sure we have a valid mask available if app.ARGS.mask: - run.command('mrconvert ' + path.from_user(app.ARGS.mask) + ' ' + path.to_scratch('mask.mif')) if not image.match('input.mif', 'mask.mif', up_to_dim=3): raise MRtrixError('Provided mask image does not match input DWI') else: - app.debug('Deriving brain mask...') run.command('dwi2mask ' + CONFIG.get('Dwi2maskAlgorithm', 'legacy') + ' input.mif mask.mif') # Determine whether we are working with single-shell or multi-shell data @@ -83,6 +108,10 @@ def execute(): #pylint: disable=unused-variable for value in image.mrinfo('input.mif', 'shell_bvalues') \ .strip().split()] multishell = (len(bvalues) > 2) + if lmax is None: + lmax = LMAXES_MULTI if multishell else LMAXES_SINGLE + elif len(lmax) == 3 and not multishell: + raise MRtrixError('User specified 3 lmax values for three-tissue decomposition, but input DWI is not multi-shell') # RF estimation and multi-tissue CSD class Tissue(object): #pylint: disable=useless-object-inheritance @@ -94,52 +123,41 @@ def __init__(self, name): tissues = [Tissue('WM'), Tissue('GM'), Tissue('CSF')] - app.debug('Estimating response function using initial brain mask...') - run.command('dwi2response dhollander ' - + 'input.mif ' - + '-mask mask.mif ' + run.command('dwi2response dhollander input.mif -mask mask.mif ' + ' '.join(tissue.tissue_rf for tissue in tissues)) - # Remove GM if we can't deal with it - if app.ARGS.lmax: - lmaxes = app.ARGS.lmax - else: - lmaxes = LMAXES_MULTI - if not multishell: - app.cleanup(tissues[1].tissue_rf) - tissues = tissues[::2] - lmaxes = LMAXES_SINGLE - - app.debug('FOD estimation with lmaxes ' + lmaxes + '...') - run.command('dwi2fod msmt_csd ' - + 'input.mif' - + ' -lmax ' + lmaxes + # Immediately remove GM if we can't deal with it + if not multishell: + app.cleanup(tissues[1].tissue_rf) + tissues = tissues[::2] + + run.command('dwi2fod msmt_csd input.mif' + + ' -lmax ' + ','.join(str(item) for item in lmax) + ' ' + ' '.join(tissue.tissue_rf + ' ' + tissue.fod - for tissue in tissues)) + for tissue in tissues)) # Normalisation in brain mask - run.command('maskfilter' - + ' mask.mif' - + ' erode - |' + run.command('maskfilter mask.mif erode - |' + ' mtnormalise -mask - -balanced' + ' -check_factors factors.txt ' + ' '.join(tissue.fod + ' ' + tissue.fod_norm - for tissue in tissues)) + for tissue in tissues)) app.cleanup([tissue.fod for tissue in tissues]) app.cleanup([tissue.fod_norm for tissue in tissues]) - app.debug('applying appropiate scaling factor to DWI series...') csf_rf = matrix.load_matrix(tissues[-1].tissue_rf) - csf_rf_bzero_lzero = csf_rf[0][0] app.cleanup([tissue.tissue_rf for tissue in tissues]) + csf_rf_bzero_lzero = csf_rf[0][0] balance_factors = matrix.load_vector('factors.txt') - csf_balance_factor = balance_factors[-1] app.cleanup('factors.txt') + csf_balance_factor = balance_factors[-1] scale_multiplier = (1000.0 * math.sqrt(4.0*math.pi)) / (csf_rf_bzero_lzero / csf_balance_factor) - run.command('mrcalc ' - + 'input.mif ' - + str(scale_multiplier) - + ' -mult ' - + path.from_user(app.ARGS.output), + + run.command('mrcalc input.mif ' + str(scale_multiplier) + ' -mult - | ' + + 'mrconvert - ' + path.from_user(app.ARGS.output), + mrconvert_keyval=path.from_user(app.ARGS.input, False), force=app.FORCE_OVERWRITE) + + if app.ARGS.scale: + matrix.save_vector(path.from_user(app.ARGS.scale, False), [scale_multiplier]) From 4903e862286b46e9fdd8d1f433e6e9da547a05ea Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Sun, 26 Feb 2023 21:42:02 +1100 Subject: [PATCH 11/24] dwibiasnormmask & "mtnorm" algorithms: Add references --- bin/dwibiasnormmask | 12 ++++++++++++ docs/reference/commands/dwi2mask.rst | 4 ++++ docs/reference/commands/dwibiascorrect.rst | 6 ++++++ docs/reference/commands/dwibiasnormmask.rst | 8 ++++++++ docs/reference/commands/dwinormalise.rst | 8 ++++++++ lib/mrtrix3/dwi2mask/mtnorm.py | 6 ++++++ lib/mrtrix3/dwibiascorrect/mtnorm.py | 9 +++++++++ lib/mrtrix3/dwinormalise/mtnorm.py | 12 ++++++++++++ 8 files changed, 65 insertions(+) diff --git a/bin/dwibiasnormmask b/bin/dwibiasnormmask index df4d4b0784..eaca253619 100755 --- a/bin/dwibiasnormmask +++ b/bin/dwibiasnormmask @@ -72,6 +72,18 @@ def usage(cmdline): #pylint: disable=unused-variable '"synthstrip": Invoke the FreeSurfer SynthStrip method. ' '"threshold": Apply a fixed partial volume threshold of 0.5 to the tissue ODF sum image ' ' (and some heuristic cleanup of the resulting mask will be performed).') + cmdline.add_citation('Jeurissen, B; Tournier, J-D; Dhollander, T; Connelly, A & Sijbers, J. ' + 'Multi-tissue constrained spherical deconvolution for improved analysis of multi-shell diffusion MRI data. ' + 'NeuroImage, 2014, 103, 411-426') + cmdline.add_citation('Raffelt, D.; Dhollander, T.; Tournier, J.-D.; Tabbara, R.; Smith, R. E.; Pierre, E. & Connelly, A. ' + 'Bias Field Correction and Intensity Normalisation for Quantitative Analysis of Apparent Fibre Density. ' + 'In Proc. ISMRM, 2017, 26, 3541') + cmdline.add_citation('Dhollander, T.; Raffelt, D. & Connelly, A. ' + 'Unsupervised 3-tissue response function estimation from single-shell or multi-shell diffusion MR data without a co-registered T1 image. ' + 'ISMRM Workshop on Breaking the Barriers of Diffusion MRI, 2016, 5') + cmdline.add_citation('Dhollander, T.; Tabbara, R.; Rosnarho-Tornstrand, J.; Tournier, J.-D.; Raffelt, D. & Connelly, A. ' + 'Multi-tissue log-domain intensity and inhomogeneity normalisation for quantitative apparent fibre density. ' + 'In Proc. ISMRM, 2021, 29, 2472') cmdline.add_argument('input', help='The input DWI series to be corrected') cmdline.add_argument('output_dwi', help='The output corrected DWI series') cmdline.add_argument('output_mask', help='The output DWI mask') diff --git a/docs/reference/commands/dwi2mask.rst b/docs/reference/commands/dwi2mask.rst index b95bd071ef..9ca6e43f6f 100644 --- a/docs/reference/commands/dwi2mask.rst +++ b/docs/reference/commands/dwi2mask.rst @@ -1047,6 +1047,10 @@ Standard options References ^^^^^^^^^^ +* Jeurissen, B; Tournier, J-D; Dhollander, T; Connelly, A & Sijbers, J. Multi-tissue constrained spherical deconvolution for improved analysis of multi-shell diffusion MRI data. NeuroImage, 2014, 103, 411-426 + +* Dhollander, T.; Raffelt, D. & Connelly, A. Unsupervised 3-tissue response function estimation from single-shell or multi-shell diffusion MR data without a co-registered T1 image. ISMRM Workshop on Breaking the Barriers of Diffusion MRI, 2016, 5 + Tournier, J.-D.; Smith, R. E.; Raffelt, D.; Tabbara, R.; Dhollander, T.; Pietsch, M.; Christiaens, D.; Jeurissen, B.; Yeh, C.-H. & Connelly, A. MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation. NeuroImage, 2019, 202, 116137 -------------- diff --git a/docs/reference/commands/dwibiascorrect.rst b/docs/reference/commands/dwibiascorrect.rst index 667c675404..a4be844bb7 100644 --- a/docs/reference/commands/dwibiascorrect.rst +++ b/docs/reference/commands/dwibiascorrect.rst @@ -378,6 +378,12 @@ Standard options References ^^^^^^^^^^ +* Jeurissen, B; Tournier, J-D; Dhollander, T; Connelly, A & Sijbers, J. Multi-tissue constrained spherical deconvolution for improved analysis of multi-shell diffusion MRI data. NeuroImage, 2014, 103, 411-426 + +* Raffelt, D.; Dhollander, T.; Tournier, J.-D.; Tabbara, R.; Smith, R. E.; Pierre, E. & Connelly, A. Bias Field Correction and Intensity Normalisation for Quantitative Analysis of Apparent Fibre Density. In Proc. ISMRM, 2017, 26, 3541 + +* Dhollander, T.; Tabbara, R.; Rosnarho-Tornstrand, J.; Tournier, J.-D.; Raffelt, D. & Connelly, A. Multi-tissue log-domain intensity and inhomogeneity normalisation for quantitative apparent fibre density. In Proc. ISMRM, 2021, 29, 2472 + Tournier, J.-D.; Smith, R. E.; Raffelt, D.; Tabbara, R.; Dhollander, T.; Pietsch, M.; Christiaens, D.; Jeurissen, B.; Yeh, C.-H. & Connelly, A. MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation. NeuroImage, 2019, 202, 116137 -------------- diff --git a/docs/reference/commands/dwibiasnormmask.rst b/docs/reference/commands/dwibiasnormmask.rst index 56d66e83b6..3112be9b0f 100644 --- a/docs/reference/commands/dwibiasnormmask.rst +++ b/docs/reference/commands/dwibiasnormmask.rst @@ -99,6 +99,14 @@ Standard options References ^^^^^^^^^^ +* Jeurissen, B; Tournier, J-D; Dhollander, T; Connelly, A & Sijbers, J. Multi-tissue constrained spherical deconvolution for improved analysis of multi-shell diffusion MRI data. NeuroImage, 2014, 103, 411-426 + +* Raffelt, D.; Dhollander, T.; Tournier, J.-D.; Tabbara, R.; Smith, R. E.; Pierre, E. & Connelly, A. Bias Field Correction and Intensity Normalisation for Quantitative Analysis of Apparent Fibre Density. In Proc. ISMRM, 2017, 26, 3541 + +* Dhollander, T.; Raffelt, D. & Connelly, A. Unsupervised 3-tissue response function estimation from single-shell or multi-shell diffusion MR data without a co-registered T1 image. ISMRM Workshop on Breaking the Barriers of Diffusion MRI, 2016, 5 + +* Dhollander, T.; Tabbara, R.; Rosnarho-Tornstrand, J.; Tournier, J.-D.; Raffelt, D. & Connelly, A. Multi-tissue log-domain intensity and inhomogeneity normalisation for quantitative apparent fibre density. In Proc. ISMRM, 2021, 29, 2472 + Tournier, J.-D.; Smith, R. E.; Raffelt, D.; Tabbara, R.; Dhollander, T.; Pietsch, M.; Christiaens, D.; Jeurissen, B.; Yeh, C.-H. & Connelly, A. MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation. NeuroImage, 2019, 202, 116137 -------------- diff --git a/docs/reference/commands/dwinormalise.rst b/docs/reference/commands/dwinormalise.rst index bf7c3ff2d3..6ecfeba737 100644 --- a/docs/reference/commands/dwinormalise.rst +++ b/docs/reference/commands/dwinormalise.rst @@ -338,6 +338,14 @@ Standard options References ^^^^^^^^^^ +* Jeurissen, B; Tournier, J-D; Dhollander, T; Connelly, A & Sijbers, J. Multi-tissue constrained spherical deconvolution for improved analysis of multi-shell diffusion MRI data. NeuroImage, 2014, 103, 411-426 + +* Raffelt, D.; Dhollander, T.; Tournier, J.-D.; Tabbara, R.; Smith, R. E.; Pierre, E. & Connelly, A. Bias Field Correction and Intensity Normalisation for Quantitative Analysis of Apparent Fibre Density. In Proc. ISMRM, 2017, 26, 3541 + +* Dhollander, T.; Tabbara, R.; Rosnarho-Tornstrand, J.; Tournier, J.-D.; Raffelt, D. & Connelly, A. Multi-tissue log-domain intensity and inhomogeneity normalisation for quantitative apparent fibre density. In Proc. ISMRM, 2021, 29, 2472 + +* Dhollander, T.; Raffelt, D. & Connelly, A. Unsupervised 3-tissue response function estimation from single-shell or multi-shell diffusion MR data without a co-registered T1 image. ISMRM Workshop on Breaking the Barriers of Diffusion MRI, 2016, 5 + Tournier, J.-D.; Smith, R. E.; Raffelt, D.; Tabbara, R.; Dhollander, T.; Pietsch, M.; Christiaens, D.; Jeurissen, B.; Yeh, C.-H. & Connelly, A. MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation. NeuroImage, 2019, 202, 116137 -------------- diff --git a/lib/mrtrix3/dwi2mask/mtnorm.py b/lib/mrtrix3/dwi2mask/mtnorm.py index 3c2452e3cb..2c51515581 100644 --- a/lib/mrtrix3/dwi2mask/mtnorm.py +++ b/lib/mrtrix3/dwi2mask/mtnorm.py @@ -39,6 +39,12 @@ def usage(base_parser, subparsers): #pylint: disable=unused-variable parser.add_description('The ODFs estimated within this optimisation procedure are by default of lower maximal spherical harmonic ' 'degree than what would be advised for analysis. This is done for computational efficiency. This ' 'behaviour can be modified through the -lmax command-line option.') + parser.add_citation('Jeurissen, B; Tournier, J-D; Dhollander, T; Connelly, A & Sijbers, J. ' + 'Multi-tissue constrained spherical deconvolution for improved analysis of multi-shell diffusion MRI data. ' + 'NeuroImage, 2014, 103, 411-426') + parser.add_citation('Dhollander, T.; Raffelt, D. & Connelly, A. ' + 'Unsupervised 3-tissue response function estimation from single-shell or multi-shell diffusion MR data without a co-registered T1 image. ' + 'ISMRM Workshop on Breaking the Barriers of Diffusion MRI, 2016, 5') parser.add_argument('input', help='The input DWI series') parser.add_argument('output', help='The output mask image') options = parser.add_argument_group('Options specific to the "mtnorm" algorithm') diff --git a/lib/mrtrix3/dwibiascorrect/mtnorm.py b/lib/mrtrix3/dwibiascorrect/mtnorm.py index 0d54872ad2..8d5eb2911e 100644 --- a/lib/mrtrix3/dwibiascorrect/mtnorm.py +++ b/lib/mrtrix3/dwibiascorrect/mtnorm.py @@ -35,6 +35,15 @@ def usage(base_parser, subparsers): #pylint: disable=unused-variable parser.add_description('The ODFs estimated within this optimisation procedure are by default of lower maximal spherical harmonic ' 'degree than what would be advised for analysis. This is done for computational efficiency. This ' 'behaviour can be modified through the -lmax command-line option.') + parser.add_citation('Jeurissen, B; Tournier, J-D; Dhollander, T; Connelly, A & Sijbers, J. ' + 'Multi-tissue constrained spherical deconvolution for improved analysis of multi-shell diffusion MRI data. ' + 'NeuroImage, 2014, 103, 411-426') + parser.add_citation('Raffelt, D.; Dhollander, T.; Tournier, J.-D.; Tabbara, R.; Smith, R. E.; Pierre, E. & Connelly, A. ' + 'Bias Field Correction and Intensity Normalisation for Quantitative Analysis of Apparent Fibre Density. ' + 'In Proc. ISMRM, 2017, 26, 3541') + parser.add_citation('Dhollander, T.; Tabbara, R.; Rosnarho-Tornstrand, J.; Tournier, J.-D.; Raffelt, D. & Connelly, A. ' + 'Multi-tissue log-domain intensity and inhomogeneity normalisation for quantitative apparent fibre density. ' + 'In Proc. ISMRM, 2021, 29, 2472') parser.add_argument('input', help='The input image series to be corrected') parser.add_argument('output', help='The output corrected image series') options = parser.add_argument_group('Options specific to the "mtnorm" algorithm') diff --git a/lib/mrtrix3/dwinormalise/mtnorm.py b/lib/mrtrix3/dwinormalise/mtnorm.py index d472c62729..b255af84ce 100644 --- a/lib/mrtrix3/dwinormalise/mtnorm.py +++ b/lib/mrtrix3/dwinormalise/mtnorm.py @@ -37,6 +37,18 @@ def usage(base_parser, subparsers): #pylint: disable=unused-variable parser.add_description('The ODFs estimated within this optimisation procedure are by default of lower maximal spherical harmonic ' 'degree than what would be advised for analysis. This is done for computational efficiency. This ' 'behaviour can be modified through the -lmax command-line option.') + parser.add_citation('Jeurissen, B; Tournier, J-D; Dhollander, T; Connelly, A & Sijbers, J. ' + 'Multi-tissue constrained spherical deconvolution for improved analysis of multi-shell diffusion MRI data. ' + 'NeuroImage, 2014, 103, 411-426') + parser.add_citation('Raffelt, D.; Dhollander, T.; Tournier, J.-D.; Tabbara, R.; Smith, R. E.; Pierre, E. & Connelly, A. ' + 'Bias Field Correction and Intensity Normalisation for Quantitative Analysis of Apparent Fibre Density. ' + 'In Proc. ISMRM, 2017, 26, 3541') + parser.add_citation('Dhollander, T.; Tabbara, R.; Rosnarho-Tornstrand, J.; Tournier, J.-D.; Raffelt, D. & Connelly, A. ' + 'Multi-tissue log-domain intensity and inhomogeneity normalisation for quantitative apparent fibre density. ' + 'In Proc. ISMRM, 2021, 29, 2472') + parser.add_citation('Dhollander, T.; Raffelt, D. & Connelly, A. ' + 'Unsupervised 3-tissue response function estimation from single-shell or multi-shell diffusion MR data without a co-registered T1 image. ' + 'ISMRM Workshop on Breaking the Barriers of Diffusion MRI, 2016, 5') parser.add_argument('input', help='The input DWI series') parser.add_argument('output', help='The normalised DWI series') options = parser.add_argument_group('Options specific to the "mtnorm" algorithm') From 73704d652052f82b81c364d0bdec4e2f28841291 Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Thu, 16 Mar 2023 11:03:59 +1100 Subject: [PATCH 12/24] Docs: Resolve dwi2mask mtnorm and synthstrip Resolves between #2604 and #2583. --- docs/reference/commands/dwi2mask.rst | 72 ++++++++++++++-------------- 1 file changed, 36 insertions(+), 36 deletions(-) diff --git a/docs/reference/commands/dwi2mask.rst b/docs/reference/commands/dwi2mask.rst index 9ca6e43f6f..81e5118b72 100644 --- a/docs/reference/commands/dwi2mask.rst +++ b/docs/reference/commands/dwi2mask.rst @@ -862,22 +862,22 @@ See the Mozilla Public License v. 2.0 for more details. For more details, see http://www.mrtrix.org/. -.. _dwi2mask_synthstrip: +.. _dwi2mask_mtnorm: -dwi2mask synthstrip -=================== +dwi2mask mtnorm +=============== Synopsis -------- -Use the FreeSurfer Synthstrip method on the mean b=0 image +Derives a DWI brain mask by calculating and then thresholding a sum-of-tissue-densities image Usage ----- :: - dwi2mask synthstrip input output [ options ] + dwi2mask mtnorm input output [ options ] - *input*: The input DWI series - *output*: The output mask image @@ -885,23 +885,25 @@ Usage Description ----------- -This algorithm requires that the SynthStrip method be installed and available via PATH; this could be via Freesufer 7.3.0 or later, or the dedicated Singularity container. +This script attempts to identify brain voxels based on the total density of macroscopic tissues as estimated through multi-tissue deconvolution. Following response function estimation and multi-tissue deconvolution, the sum of tissue densities is thresholded at a fixed value (default is 0.5), and subsequent mask image cleaning operations are performed. + +The operation of this script is a subset of that performed by the script "dwibiasnormmask". Many users may find that comprehensive solution preferable; this dwi2mask algorithm is nevertheless provided to demonstrate specifically the mask estimation portion of that command. + +The ODFs estimated within this optimisation procedure are by default of lower maximal spherical harmonic degree than what would be advised for analysis. This is done for computational efficiency. This behaviour can be modified through the -lmax command-line option. Options ------- -Options specific to the 'Synthstrip' algorithm -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -- **-stripped** The output stripped image +Options specific to the "mtnorm" algorithm +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -- **-gpu** Use the GPU +- **-init_mask image** Provide an initial brain mask, which will constrain the response function estimation (if omitted, the default dwi2mask algorithm will be used) -- **-model file** Alternative model weights +- **-lmax values** The maximum spherical harmonic degree for the estimated FODs (see Description); defaults are "4,0,0" for multi-shell and "4,0" for single-shell data) -- **-nocsf** Compute the immediate boundary of brain matter excluding surrounding CSF +- **-threshold value** the threshold on the total tissue density sum image used to derive the brain mask; default is 0.5 -- **-border** Control the boundary distance from the brain +- **-tissuesum image** Export the tissue sum image that was used to generate the mask Options for importing the diffusion gradient table ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -941,7 +943,9 @@ Standard options References ^^^^^^^^^^ -* A. Hoopes, J. S. Mora, A. V. Dalca, B. Fischl, M. Hoffmann. SynthStrip: Skull-Stripping for Any Brain Image. NeuroImage, 2022, 260, 119474 +* Jeurissen, B; Tournier, J-D; Dhollander, T; Connelly, A & Sijbers, J. Multi-tissue constrained spherical deconvolution for improved analysis of multi-shell diffusion MRI data. NeuroImage, 2014, 103, 411-426 + +* Dhollander, T.; Raffelt, D. & Connelly, A. Unsupervised 3-tissue response function estimation from single-shell or multi-shell diffusion MR data without a co-registered T1 image. ISMRM Workshop on Breaking the Barriers of Diffusion MRI, 2016, 5 Tournier, J.-D.; Smith, R. E.; Raffelt, D.; Tabbara, R.; Dhollander, T.; Pietsch, M.; Christiaens, D.; Jeurissen, B.; Yeh, C.-H. & Connelly, A. MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation. NeuroImage, 2019, 202, 116137 @@ -949,7 +953,7 @@ Tournier, J.-D.; Smith, R. E.; Raffelt, D.; Tabbara, R.; Dhollander, T.; Pietsch -**Author:** Ruobing Chen (chrc@student.unimelb.edu.au) +**Author:** Robert E. Smith (robert.smith@florey.edu.au) and Arshiya Sangchooli (asangchooli@student.unimelb.edu.au) **Copyright:** Copyright (c) 2008-2023 the MRtrix3 contributors. @@ -966,22 +970,22 @@ See the Mozilla Public License v. 2.0 for more details. For more details, see http://www.mrtrix.org/. -.. _dwi2mask_mtnorm: +.. _dwi2mask_synthstrip: -dwi2mask mtnorm -=============== +dwi2mask synthstrip +=================== Synopsis -------- -Derives a DWI brain mask by calculating and then thresholding a sum-of-tissue-densities image +Use the FreeSurfer Synthstrip method on the mean b=0 image Usage ----- :: - dwi2mask mtnorm input output [ options ] + dwi2mask synthstrip input output [ options ] - *input*: The input DWI series - *output*: The output mask image @@ -989,25 +993,23 @@ Usage Description ----------- -This script attempts to identify brain voxels based on the total density of macroscopic tissues as estimated through multi-tissue deconvolution. Following response function estimation and multi-tissue deconvolution, the sum of tissue densities is thresholded at a fixed value (default is 0.5), and subsequent mask image cleaning operations are performed. - -The operation of this script is a subset of that performed by the script "dwibiasnormmask". Many users may find that comprehensive solution preferable; this dwi2mask algorithm is nevertheless provided to demonstrate specifically the mask estimation portion of that command. - -The ODFs estimated within this optimisation procedure are by default of lower maximal spherical harmonic degree than what would be advised for analysis. This is done for computational efficiency. This behaviour can be modified through the -lmax command-line option. +This algorithm requires that the SynthStrip method be installed and available via PATH; this could be via Freesufer 7.3.0 or later, or the dedicated Singularity container. Options ------- -Options specific to the "mtnorm" algorithm -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +Options specific to the 'Synthstrip' algorithm +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -- **-init_mask image** Provide an initial brain mask, which will constrain the response function estimation (if omitted, the default dwi2mask algorithm will be used) +- **-stripped** The output stripped image -- **-lmax values** The maximum spherical harmonic degree for the estimated FODs (see Description); defaults are "4,0,0" for multi-shell and "4,0" for single-shell data) +- **-gpu** Use the GPU -- **-threshold value** the threshold on the total tissue density sum image used to derive the brain mask; default is 0.5 +- **-model file** Alternative model weights -- **-tissuesum image** Export the tissue sum image that was used to generate the mask +- **-nocsf** Compute the immediate boundary of brain matter excluding surrounding CSF + +- **-border** Control the boundary distance from the brain Options for importing the diffusion gradient table ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -1047,9 +1049,7 @@ Standard options References ^^^^^^^^^^ -* Jeurissen, B; Tournier, J-D; Dhollander, T; Connelly, A & Sijbers, J. Multi-tissue constrained spherical deconvolution for improved analysis of multi-shell diffusion MRI data. NeuroImage, 2014, 103, 411-426 - -* Dhollander, T.; Raffelt, D. & Connelly, A. Unsupervised 3-tissue response function estimation from single-shell or multi-shell diffusion MR data without a co-registered T1 image. ISMRM Workshop on Breaking the Barriers of Diffusion MRI, 2016, 5 +* A. Hoopes, J. S. Mora, A. V. Dalca, B. Fischl, M. Hoffmann. SynthStrip: Skull-Stripping for Any Brain Image. NeuroImage, 2022, 260, 119474 Tournier, J.-D.; Smith, R. E.; Raffelt, D.; Tabbara, R.; Dhollander, T.; Pietsch, M.; Christiaens, D.; Jeurissen, B.; Yeh, C.-H. & Connelly, A. MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation. NeuroImage, 2019, 202, 116137 @@ -1057,7 +1057,7 @@ Tournier, J.-D.; Smith, R. E.; Raffelt, D.; Tabbara, R.; Dhollander, T.; Pietsch -**Author:** Robert E. Smith (robert.smith@florey.edu.au) and Arshiya Sangchooli (asangchooli@student.unimelb.edu.au) +**Author:** Ruobing Chen (chrc@student.unimelb.edu.au) **Copyright:** Copyright (c) 2008-2023 the MRtrix3 contributors. From 81819d4909b45f653d3ad9a9077bd9ee953d0f73 Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Thu, 16 Mar 2023 11:08:37 +1100 Subject: [PATCH 13/24] dwinormalse: Rename "individual" algorithm to "manual" --- docs/reference/commands/dwinormalise.rst | 12 ++++++------ .../dwinormalise/{individual.py => manual.py} | 4 ++-- 2 files changed, 8 insertions(+), 8 deletions(-) rename lib/mrtrix3/dwinormalise/{individual.py => manual.py} (96%) diff --git a/docs/reference/commands/dwinormalise.rst b/docs/reference/commands/dwinormalise.rst index 6ecfeba737..ef7c320ff4 100644 --- a/docs/reference/commands/dwinormalise.rst +++ b/docs/reference/commands/dwinormalise.rst @@ -15,7 +15,7 @@ Usage dwinormalise algorithm [ options ] ... -- *algorithm*: Select the algorithm to be used to complete the script operation; additional details and options become available once an algorithm is nominated. Options are: group, individual, mtnorm +- *algorithm*: Select the algorithm to be used to complete the script operation; additional details and options become available once an algorithm is nominated. Options are: group, manual, mtnorm Description ----------- @@ -168,22 +168,22 @@ See the Mozilla Public License v. 2.0 for more details. For more details, see http://www.mrtrix.org/. -.. _dwinormalise_individual: +.. _dwinormalise_manual: -dwinormalise individual -======================= +dwinormalise manual +=================== Synopsis -------- -Intensity normalise a DWI series based on the b=0 signal within a supplied mask +Intensity normalise a DWI series based on the b=0 signal within a manually-supplied supplied mask Usage ----- :: - dwinormalise individual input_dwi input_mask output_dwi [ options ] + dwinormalise manual input_dwi input_mask output_dwi [ options ] - *input_dwi*: The input DWI series - *input_mask*: The mask within which a reference b=0 intensity will be sampled diff --git a/lib/mrtrix3/dwinormalise/individual.py b/lib/mrtrix3/dwinormalise/manual.py similarity index 96% rename from lib/mrtrix3/dwinormalise/individual.py rename to lib/mrtrix3/dwinormalise/manual.py index 5ac0500d41..14b290d56b 100644 --- a/lib/mrtrix3/dwinormalise/individual.py +++ b/lib/mrtrix3/dwinormalise/manual.py @@ -23,9 +23,9 @@ def usage(base_parser, subparsers): #pylint: disable=unused-variable - parser = subparsers.add_parser('individual', parents=[base_parser]) + parser = subparsers.add_parser('manual', parents=[base_parser]) parser.set_author('Robert E. Smith (robert.smith@florey.edu.au) and David Raffelt (david.raffelt@florey.edu.au)') - parser.set_synopsis('Intensity normalise a DWI series based on the b=0 signal within a supplied mask') + parser.set_synopsis('Intensity normalise a DWI series based on the b=0 signal within a manually-supplied supplied mask') parser.add_argument('input_dwi', help='The input DWI series') parser.add_argument('input_mask', help='The mask within which a reference b=0 intensity will be sampled') parser.add_argument('output_dwi', help='The output intensity-normalised DWI series') From cb02298395f03f975156b28d89af9cb1dad723a8 Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Thu, 16 Mar 2023 11:54:57 +1100 Subject: [PATCH 14/24] dwibiasnormmask: New config file option "DwibiasnormmaskMaskAlgorithm" This config file option will allow users to specify the default algorithm that will be used wthin the iterative process to derive the brain mask. This is a different step to derivation of an initial brain mask for response function estimation / MSMT CSD / mtnormalise, for which dwi2mask must be used. Rather, during the iterative process, one has access to the ODF tissue sum image, which may be used for brain masking, and so it may be preferable for something other than dwi2mask to be used taking that image as input. --- bin/dwibiasnormmask | 43 +++++++++++++-------- core/file/pyconfig.h | 20 ++++++++++ docs/reference/commands/dwibiasnormmask.rst | 2 +- docs/reference/config_file_options.rst | 22 +++++++++++ 4 files changed, 70 insertions(+), 17 deletions(-) diff --git a/bin/dwibiasnormmask b/bin/dwibiasnormmask index eaca253619..523d2aef4c 100755 --- a/bin/dwibiasnormmask +++ b/bin/dwibiasnormmask @@ -21,7 +21,6 @@ DWIBIASCORRECT_MAX_ITERS = 2 LMAXES_MULTI = [4,0,0] LMAXES_SINGLE = [4,0] REFERENCE_INTENSITY = 1000.0 -# TODO Add config file entry to specify user preferred default MASK_ALGOS = ['dwi2mask', 'fslbet', 'hdbet', 'mrthreshold', 'synthstrip', 'threshold'] MASK_ALGO_DEFAULT = 'threshold' DICE_COEFF_DEFAULT = 1.0 - 1e-3 @@ -65,12 +64,13 @@ def usage(cmdline): #pylint: disable=unused-variable '"dwi2mask": Use the MRtrix3 command dwi2mask on the bias-field-corrected DWI series ' '(ie. do not use the ODF tissue sum image for mask derivation); ' 'the algorithm to be invoked can be controlled by the user via the MRtrix config file entry "Dwi2maskAlgorithm". ' - '"fslbet": Invoke the FSL command "bet". ' - '"hdbet": Invoke the HD-BET command. ' - '"mrthreshold": Invoke the MRtrix3 command "mrthreshold", where an appropriate threshold value will be determined automatically ' - ' (and some heuristic cleanup of the resulting mask will be performed). ' - '"synthstrip": Invoke the FreeSurfer SynthStrip method. ' - '"threshold": Apply a fixed partial volume threshold of 0.5 to the tissue ODF sum image ' + '"fslbet": Invoke the FSL command "bet" on the ODF tissue sum image. ' + '"hdbet": Invoke the HD-BET command on the ODF tissue sum image. ' + '"mrthreshold": Invoke the MRtrix3 command "mrthreshold" on the ODF tissue sum image, ' + 'where an appropriate threshold value will be determined automatically ' + '(and some heuristic cleanup of the resulting mask will be performed). ' + '"synthstrip": Invoke the FreeSurfer SynthStrip method on the ODF tissue sum image. ' + '"threshold": Apply a fixed partial volume threshold of 0.5 to the ODF tissue sum image ' ' (and some heuristic cleanup of the resulting mask will be performed).') cmdline.add_citation('Jeurissen, B; Tournier, J-D; Dhollander, T; Connelly, A & Sijbers, J. ' 'Multi-tissue constrained spherical deconvolution for improved analysis of multi-shell diffusion MRI data. ' @@ -106,7 +106,7 @@ def usage(cmdline): #pylint: disable=unused-variable internal_options.add_argument('-max_iters', type=int, default=DWIBIASCORRECT_MAX_ITERS, metavar='count', help='The maximum number of iterations (see Description); default is ' + str(DWIBIASCORRECT_MAX_ITERS) + '; ' 'set to 0 to proceed until convergence') - internal_options.add_argument('-mask_algo', choices=MASK_ALGOS, default=MASK_ALGO_DEFAULT, metavar='algorithm', + internal_options.add_argument('-mask_algo', choices=MASK_ALGOS, metavar='algorithm', help='The algorithm to use for mask estimation, potentially based on the ODF sum image (see Description); default: ' + MASK_ALGO_DEFAULT) internal_options.add_argument('-lmax', metavar='values', help='The maximum spherical harmonic degree for the estimated FODs (see Description); ' @@ -137,6 +137,18 @@ def execute(): #pylint: disable=unused-variable if app.ARGS.reference <= 0.0: raise MRtrixError('Reference intensity must be positive') + # Check what masking agorithm is going to be used + mask_algo = MASK_ALGO_DEFAULT + if app.ARGS.mask_algo: + mask_algo = app.ARGS.mask_algo + elif 'DwibiasnormmaskMaskAlgorithm' in CONFIG: + mask_algo = CONFIG['DwibiasnormmaskMaskAlgorithm'] + if not mask_algo in MASK_ALGOS: + raise MRtrixError('Invalid masking algorithm selection "%s" in MRtrix config file' % mask_algo) + app.console('"%s" algorithm will be used for brain masking during iteration as specified in config file' % mask_algo) + else: + app.console('Default "%s" algorithm will be used for brain masking during iteration' % MASK_ALGO_DEFAULT) + # Check mask algorithm, including availability of external software if necessary for mask_algo, software, command in [('fslbet', 'FSL', 'bet'), ('hdbet', 'HD-BET', 'hd-bet'), ('synthstrip', 'FreeSurfer', 'mri_synthstrip')]: if app.ARGS.mask_algo == mask_algo and not shutil.which(command): @@ -195,7 +207,6 @@ def execute(): #pylint: disable=unused-variable + CONFIG.get('Dwi2maskAlgorithm', 'legacy') + ' input.mif dwi_mask_init.mif') - # Combined RF estimation / CSD / mtnormalise / mask revision class Tissue(object): #pylint: disable=useless-object-inheritance def __init__(self, name, index): @@ -315,18 +326,18 @@ def execute(): #pylint: disable=unused-variable new_dwi_mask_image = 'dwi_mask' + iter_string + '.mif' tissue_sum_image_nii = None new_dwi_mask_image_nii = None - if app.ARGS.mask_algo in ['fslbet', 'hdbet', 'synthstrip']: + if mask_algo in ['fslbet', 'hdbet', 'synthstrip']: tissue_sum_image_nii = os.path.splitext(tissue_sum_image)[0] + '.nii' run.command('mrconvert ' + tissue_sum_image + ' ' + tissue_sum_image_nii) new_dwi_mask_image_nii = os.path.splitext(new_dwi_mask_image)[0] + '.nii' - if app.ARGS.mask_algo == 'dwi2mask': + if mask_algo == 'dwi2mask': run.command('dwi2mask ' + CONFIG.get('Dwi2maskAlgorithm', 'legacy') + ' ' + new_dwi_image + ' ' + new_dwi_mask_image) - elif app.ARGS.mask_algo == 'fslbet': + elif mask_algo == 'fslbet': run.command('bet ' + tissue_sum_image_nii + ' ' + new_dwi_mask_image_nii + ' -R -m') app.cleanup(fsl.find_image(os.path.splitext(new_dwi_mask_image_nii)[0])) new_dwi_mask_image_nii = fsl.find_image(os.path.splitext(new_dwi_mask_image_nii)[0] + '_mask') run.command('mrcalc ' + new_dwi_mask_image_nii + ' input_pos_mask.mif -mult ' + new_dwi_mask_image) - elif app.ARGS.mask_algo == 'hdbet': + elif mask_algo == 'hdbet': try: run.command('hd-bet -i ' + tissue_sum_image_nii) except run.MRtrixCmdError as e_gpu: @@ -336,8 +347,8 @@ def execute(): #pylint: disable=unused-variable raise run.MRtrixCmdError('hd-bet', 1, e_gpu.stdout + e_cpu.stdout, e_gpu.stderr + e_cpu.stderr) new_dwi_mask_image_nii = os.path.splitext(tissue_sum_image)[0] + '_bet_mask.nii.gz' run.command('mrcalc ' + new_dwi_mask_image_nii + ' input_pos_mask.mif -mult ' + new_dwi_mask_image) - elif app.ARGS.mask_algo in ['mrthreshold', 'threshold']: - mrthreshold_abs_option = ' -abs 0.5' if app.ARGS.mask_algo == 'threshold' else '' + elif mask_algo in ['mrthreshold', 'threshold']: + mrthreshold_abs_option = ' -abs 0.5' if mask_algo == 'threshold' else '' run.command('mrthreshold ' + tissue_sum_image + mrthreshold_abs_option @@ -350,7 +361,7 @@ def execute(): #pylint: disable=unused-variable + ' mrcalc - input_pos_mask.mif -mult ' + new_dwi_mask_image + ' -datatype bit') - elif app.ARGS.mask_algo == 'synthstrip': + elif mask_algo == 'synthstrip': run.command('mri_synthstrip -i ' + tissue_sum_image_nii + ' --mask ' + new_dwi_mask_image_nii) run.command('mrcalc ' + new_dwi_mask_image_nii + ' input_pos_mask.mif -mult ' + new_dwi_mask_image) else: diff --git a/core/file/pyconfig.h b/core/file/pyconfig.h index c0add9c524..a0a1fafd6d 100644 --- a/core/file/pyconfig.h +++ b/core/file/pyconfig.h @@ -71,6 +71,26 @@ //CONF Dwi2maskTemplate is set to "fsl"), specify the configuration //CONF file to be provided to the FSL command fnirt. +//CONF option: DwibiasnormmaskMaskAlgorithm +//CONF default: threshold +//CONF For the "dwibiasnormmask" command, specify the algorithm that +//CONF will be used for brain masking across the iterative process. +//CONF Note that the image data that are used to derive the mask +//CONF may be different between the various options. +//CONF Available options are as follows. +//CONF dwi2mask: Invoke the MRtrix3 command dwi2mask using the +//CONF bias-field-corrected DWI series, using whatever algorithm has +//CONF been specified as the default for that command +//CONF (see config file option Dwi2maskAlgorithm). +//CONF fslbet: Use FSL command "bet" on the ODF tissue sum image. +//CONF hdbet: Use HD-BET on the ODF tissue sum image. +//CONF mrthreshold: Use MRtrix3 command "mrthreshold" on the ODF tissue +//CONF sum image, allowing it to determine the appropriate threshold +//CONF automatically; some mask cleanup operations will additionally be used. +//CONF synthstrip: Use FreeSurfer's SynthStrip on the ODF tissue sum image. +//CONF threshold: Apply a 0.5 threshold to the ODF tissue sum image; +//CONF some mask cleanup operations will additionally be used. + //CONF option: ScriptScratchDir //CONF default: `.` //CONF The location in which to generate the scratch directories to be diff --git a/docs/reference/commands/dwibiasnormmask.rst b/docs/reference/commands/dwibiasnormmask.rst index 3112be9b0f..0df9e89030 100644 --- a/docs/reference/commands/dwibiasnormmask.rst +++ b/docs/reference/commands/dwibiasnormmask.rst @@ -32,7 +32,7 @@ The ODFs estimated within this optimisation procedure are by default of lower ma By default, the optimisation procedure will terminate after only two iterations. This is done because it has been observed for some data / configurations that additional iterations can lead to unstable divergence and erroneous results for bias field estimation and masking. For other configurations, it may be preferable to use a greater number of iterations, and allow the iterative algorithm to converge to a stable solution. This can be controlled via the -max_iters command-line option. -Within the optimisation algorithm, derivation of the mask may potentially be performed differently to a conventional mask derivation that is based on a DWI series (where, in many instances, it is actually only the mean b=0 image that is used). Here, the image corresponding to the sum of tissue signal fractions following spherical deconvolution / bias field correction / intensity normalisation is also available, and this can potentially be used for mask derivation. Available options are as follows. "dwi2mask": Use the MRtrix3 command dwi2mask on the bias-field-corrected DWI series (ie. do not use the ODF tissue sum image for mask derivation); the algorithm to be invoked can be controlled by the user via the MRtrix config file entry "Dwi2maskAlgorithm". "fslbet": Invoke the FSL command "bet". "hdbet": Invoke the HD-BET command. "mrthreshold": Invoke the MRtrix3 command "mrthreshold", where an appropriate threshold value will be determined automatically (and some heuristic cleanup of the resulting mask will be performed). "synthstrip": Invoke the FreeSurfer SynthStrip method. "threshold": Apply a fixed partial volume threshold of 0.5 to the tissue ODF sum image (and some heuristic cleanup of the resulting mask will be performed). +Within the optimisation algorithm, derivation of the mask may potentially be performed differently to a conventional mask derivation that is based on a DWI series (where, in many instances, it is actually only the mean b=0 image that is used). Here, the image corresponding to the sum of tissue signal fractions following spherical deconvolution / bias field correction / intensity normalisation is also available, and this can potentially be used for mask derivation. Available options are as follows. "dwi2mask": Use the MRtrix3 command dwi2mask on the bias-field-corrected DWI series (ie. do not use the ODF tissue sum image for mask derivation); the algorithm to be invoked can be controlled by the user via the MRtrix config file entry "Dwi2maskAlgorithm". "fslbet": Invoke the FSL command "bet" on the ODF tissue sum image. "hdbet": Invoke the HD-BET command on the ODF tissue sum image. "mrthreshold": Invoke the MRtrix3 command "mrthreshold" on the ODF tissue sum image, where an appropriate threshold value will be determined automatically (and some heuristic cleanup of the resulting mask will be performed). "synthstrip": Invoke the FreeSurfer SynthStrip method on the ODF tissue sum image. "threshold": Apply a fixed partial volume threshold of 0.5 to the ODF tissue sum image (and some heuristic cleanup of the resulting mask will be performed). Options ------- diff --git a/docs/reference/config_file_options.rst b/docs/reference/config_file_options.rst index bf53cc4a80..cee14dc07d 100644 --- a/docs/reference/config_file_options.rst +++ b/docs/reference/config_file_options.rst @@ -277,6 +277,28 @@ List of MRtrix3 configuration file options The software to be used for registration and transformation by default within the "dwi2mask template" algorithm. +.. option:: DwibiasnormmaskMaskAlgorithm + + *default: threshold* + + For the "dwibiasnormmask" command, specify the algorithm that + will be used for brain masking across the iterative process. + Note that the image data that are used to derive the mask + may be different between the various options. + Available options are as follows. + dwi2mask: Invoke the MRtrix3 command dwi2mask using the + bias-field-corrected DWI series, using whatever algorithm has + been specified as the default for that command + (see config file option Dwi2maskAlgorithm). + fslbet: Use FSL command "bet" on the ODF tissue sum image. + hdbet: Use HD-BET on the ODF tissue sum image. + mrthreshold: Use MRtrix3 command "mrthreshold" on the ODF tissue + sum image, allowing it to determine the appropriate threshold + automatically; some mask cleanup operations will additionally be used. + synthstrip: Use FreeSurfer's SynthStrip on the ODF tissue sum image. + threshold: Apply a 0.5 threshold to the ODF tissue sum image; + some mask cleanup operations will additionally be used. + .. option:: FailOnWarn *default: 0 (false)* From 1b34d362a5f0296a94f9ad94c2cb2d377b3e620f Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Fri, 17 Mar 2023 13:21:30 +1100 Subject: [PATCH 15/24] Docs: Global intensity normalisation: Reflect command renaming "dwinormalise individual" was renamed to "dwinormalise manual" in 81819d49. --- docs/concepts/global_intensity_normalisation.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/concepts/global_intensity_normalisation.rst b/docs/concepts/global_intensity_normalisation.rst index 31a4e662ba..9a966bda0c 100644 --- a/docs/concepts/global_intensity_normalisation.rst +++ b/docs/concepts/global_intensity_normalisation.rst @@ -46,7 +46,7 @@ folder containing the corresponding whole brain mask images (with the same filen prefix). The script runs by first computing diffusion tensor Fractional Anisotropy (FA) maps, registering these to a groupwise template, then thresholding the template FA map to obtain an *approximate* white matter mask. The mask is then transformed back into the -space of each subject image and used in the :ref:`dwinormalise_individual` command to +space of each subject image and used in the :ref:`dwinormalise_manual` command to normalise the input DW images to have the same b=0 white matter median value. All intensity normalised data will be output in a single folder. As previously mentioned, all DWI data must be bias field corrected *before* applying :ref:`dwinormalise_group`, From 3def23d27e7febc803a95250a3f7482045d3dd42 Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Wed, 22 Mar 2023 09:13:48 +1100 Subject: [PATCH 16/24] dwi2mask: Centralise definition of default algorithm --- bin/dwi2response | 2 +- bin/dwibiascorrect | 2 +- bin/dwibiasnormmask | 2 +- bin/dwifslpreproc | 2 +- bin/dwigradcheck | 2 +- lib/mrtrix3/__init__.py | 7 +++++-- lib/mrtrix3/dwinormalise/mtnorm.py | 2 +- 7 files changed, 11 insertions(+), 8 deletions(-) diff --git a/bin/dwi2response b/bin/dwi2response index 85a18606c5..8461a66378 100755 --- a/bin/dwi2response +++ b/bin/dwi2response @@ -113,7 +113,7 @@ def execute(): #pylint: disable=unused-variable raise MRtrixError('Provided mask image needs to be a 3D image') else: app.console('Computing brain mask (dwi2mask)...') - run.command('dwi2mask ' + CONFIG.get('Dwi2maskAlgorithm', 'legacy') + ' dwi.mif mask.mif', show=False) + run.command('dwi2mask ' + CONFIG['Dwi2maskAlgorithm'] + ' dwi.mif mask.mif', show=False) if not image.statistics('mask.mif', mask='mask.mif').count: raise MRtrixError(('Provided' if app.ARGS.mask else 'Generated') + ' mask image does not contain any voxels') diff --git a/bin/dwibiascorrect b/bin/dwibiascorrect index 990b311207..64387392b9 100755 --- a/bin/dwibiascorrect +++ b/bin/dwibiascorrect @@ -71,7 +71,7 @@ def execute(): #pylint: disable=unused-variable if not image.match('in.mif', 'mask.mif', up_to_dim=3): raise MRtrixError('Provided mask image does not match input DWI') else: - run.command('dwi2mask ' + CONFIG.get('Dwi2maskAlgorithm', 'legacy') + ' in.mif mask.mif') + run.command('dwi2mask ' + CONFIG['Dwi2maskAlgorithm'] + ' in.mif mask.mif') # From here, the script splits depending on what estimation algorithm is being used alg.execute() diff --git a/bin/dwibiasnormmask b/bin/dwibiasnormmask index 523d2aef4c..c98d559dce 100755 --- a/bin/dwibiasnormmask +++ b/bin/dwibiasnormmask @@ -204,7 +204,7 @@ def execute(): #pylint: disable=unused-variable else: app.debug('Performing intial DWI brain masking') run.command('dwi2mask ' - + CONFIG.get('Dwi2maskAlgorithm', 'legacy') + + CONFIG['Dwi2maskAlgorithm'] + ' input.mif dwi_mask_init.mif') # Combined RF estimation / CSD / mtnormalise / mask revision diff --git a/bin/dwifslpreproc b/bin/dwifslpreproc index 9b982834ea..cd6f7f6d7f 100755 --- a/bin/dwifslpreproc +++ b/bin/dwifslpreproc @@ -816,7 +816,7 @@ def execute(): #pylint: disable=unused-variable # Need gradient table if running dwi2mask after applytopup to derive a brain mask for eddy run.command('mrinfo dwi.mif -export_grad_mrtrix grad.b') - dwi2mask_algo = CONFIG.get('Dwi2maskAlgorithm', 'legacy') + dwi2mask_algo = CONFIG['Dwi2maskAlgorithm'] eddy_in_topup_option = '' dwi_post_eddy_crop_option = '' diff --git a/bin/dwigradcheck b/bin/dwigradcheck index 664ef3c644..91161664ae 100755 --- a/bin/dwigradcheck +++ b/bin/dwigradcheck @@ -92,7 +92,7 @@ def execute(): #pylint: disable=unused-variable # Note that gradient table must be explicitly loaded, since there may not # be one in the image header (user may be relying on -grad or -fslgrad input options) if not os.path.exists('mask.mif'): - run.command('dwi2mask ' + CONFIG.get('Dwi2maskAlgorithm', 'legacy') + ' data.mif mask.mif -grad grad.b') + run.command('dwi2mask ' + CONFIG['Dwi2maskAlgorithm'] + ' data.mif mask.mif -grad grad.b') # How many tracks are we going to generate? number_option = ' -select ' + str(app.ARGS.number) diff --git a/lib/mrtrix3/__init__.py b/lib/mrtrix3/__init__.py index f311bda411..a4f29bcd41 100644 --- a/lib/mrtrix3/__init__.py +++ b/lib/mrtrix3/__init__.py @@ -42,8 +42,11 @@ class MRtrixError(MRtrixBaseError): #pylint: disable=unused-variable EXE_LIST = [ os.path.splitext(name)[0] for name in os.listdir(BIN_PATH) ] #pylint: disable=unused-variable -# - 'CONFIG' is a directory containing those entries present in the MRtrix config files -CONFIG = { } +# 'CONFIG' is a dictionary containing those entries present in the MRtrix config files +# Can add default values here that would otherwise appear in multiple locations +CONFIG = { + 'Dwi2maskAlgorithm': 'legacy' +} # Codes for printing information to the terminal diff --git a/lib/mrtrix3/dwinormalise/mtnorm.py b/lib/mrtrix3/dwinormalise/mtnorm.py index b255af84ce..69bec7b683 100644 --- a/lib/mrtrix3/dwinormalise/mtnorm.py +++ b/lib/mrtrix3/dwinormalise/mtnorm.py @@ -112,7 +112,7 @@ def execute(): #pylint: disable=unused-variable if not image.match('input.mif', 'mask.mif', up_to_dim=3): raise MRtrixError('Provided mask image does not match input DWI') else: - run.command('dwi2mask ' + CONFIG.get('Dwi2maskAlgorithm', 'legacy') + ' input.mif mask.mif') + run.command('dwi2mask ' + CONFIG['Dwi2maskAlgorithm'] + ' input.mif mask.mif') # Determine whether we are working with single-shell or multi-shell data bvalues = [ From 854761f8b19a0b50ac5ab89279f08d64146f15ca Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Wed, 22 Mar 2023 09:16:21 +1100 Subject: [PATCH 17/24] Docs: Update DWI masking page - Add description of dwi2mask mtnorm algorithm. - Add description of dwi2mask synthstrip algorithm. - Add description of new command dwibiasnormmask. --- docs/dwi_preprocessing/masking.rst | 198 ++++++++++++++++++++++++++--- 1 file changed, 181 insertions(+), 17 deletions(-) diff --git a/docs/dwi_preprocessing/masking.rst b/docs/dwi_preprocessing/masking.rst index f29d0d6c16..e597a8a5ff 100644 --- a/docs/dwi_preprocessing/masking.rst +++ b/docs/dwi_preprocessing/masking.rst @@ -77,10 +77,10 @@ the core image registration and transformation processes: - ``fsl``: Utilises FSL_ commands as follows: - - flirt_: Initial affine registration; - - fnirt_: Non-linear registration; - - invwarp_: Inversion of warp from subject to template; - - applywarp_: Transform template mask to subject space. + - flirt_: Initial affine registration; + - fnirt_: Non-linear registration; + - invwarp_: Inversion of warp from subject to template; + - applywarp_: Transform template mask to subject space. By default, if no manual selection is made here using either the ``-software`` command-line option or the ``Dwi2maskTemplateSoftware`` configuration file @@ -175,6 +175,51 @@ Operations are as follows: 4. Apply mask cleaning filter to remove small areas only connected to the largest component via thin "bridges". +``dwi2mask mtnorm`` +^^^^^^^^^^^^^^^^^^^ + +This algorithm implements a subset of the functionalities provided in the +``dwibiasnormmask`` script (described in further detail below). +It is based on utilisation of the results generated by the ``mtnormalise`` command. +The basic premise is that, following multi-shell multi-tissue CSD and appropriate +response function bias correction / bias field cocrrection / intensity normalisation, +an image consisting of the sum of all macroscopic tissue ODFs should be approximately +1.0 in brain voxels and 0.0 in non-brain voxels. + +The order of operations is as follows: + +1. If not provided by the user, generate an initial brain mask using the default + ``dwi2mask`` algorithm. + +2. Perform three-tissue response function estimation. + +3. Perform multi-shell multi-tissue CSD + (with all three macroscopic tissues---WM, GM and CSF---if possible, + otherwise only WM and CSF) + +4. Use ``mtnormalise`` to correct: + + 1. Biases in response function magnitudes using ``-balanced`` option + (note that this functionality is *deliberately omitted* from typical + quantitative analysis pipelines as it may regress out effects of interest) + 2. Smoothly-varying bias field + 3. Global intensity scaling + +5. Calculate an image of the sum of tissue ODFs + +6. Apply a threshold to binarize this image + (default threshold is 0.5). + +7. Apply mask cleaning operations (eg. largest connected component). + +``dwi2mask synthstrip`` +^^^^^^^^^^^^^^^^^^^^^^^ + +The SynthStrip_ method is based on a deep learning neural network that has been +trained on a wide range of neuroimaging modalities and data qualities. This +algorithm provides the mean *b*\=0 image to SynthStrip, whether installed as part +of FreeSurfer (version 7.3.0 or later) or as the stand-alone Singularity container. + ``dwi2mask trace`` ^^^^^^^^^^^^^^^^^^ @@ -190,7 +235,7 @@ Its behaviour is as follows: image approximately the same intensity of that of the first shell (this is so that each shell contributes approximately equally toward determination of the mask); - + 3. Calculate the mean trace-weighted image across shells; 4. Automatically determine an intensity threshold for this image to produce @@ -203,17 +248,17 @@ Its behaviour is as follows: 7. If the command-line option ``-iterative`` is *not* used, the algorithm ceases at this point (i.e. the default behaviour); - + 8. For each *b*-value shell, compute the mean and standard deviation of the trace-weighted image intensities inside and outside of the current mask, and use this to derive Cohen's *d* statistic; - + 9. Perform a recombination of the trace-weighted images; but the multiplicative weights applied to each *b*-value shell trace image are, instead of being based on intensity matching as in step 2, the Cohen's *d* statistics calculated in step 8; - -10. Apply a threshold and mask filtering operations as in steps 4-6; + +10. Apply a threshold and mask filtering operations as in steps 4--6; 11. If the resulting mask differs from the previous estimate, go back to step 8; if not, or if a maximum number of iterations is reached, @@ -224,6 +269,8 @@ a hypothetical heuristic, and it is not yet known whether or not its behaviour is reasonable across a range of DWI data; it should therefore be considered entirely experimental. +.. _dwi2mask_algorithm_comparison: + Algorithm comparison -------------------- @@ -236,8 +283,10 @@ Algorithm comparison "``b02template``", "Yes (ANTs_ / FSL_)", "No", "Matches template", "Yes", "No" "``consensus``", "Only if installed", "Yes", "Various", "Various", "No" "``fslbet``", "Yes (FSL_)", "No", "Approx. spherical", "Yes", "No" - "``hdbet``", "Yes (HD-BET_)", "No", "Brain", "Yes", "Yes" + "``hdbet``", "Yes (HD-BET_)", "No", "Human brain", "Yes", "Yes" "``legacy``", "No", "Yes", "Single connected component", "No", "No" + "``mtnorm``", "No", "Yes", "WM / GM / CSF constituency; single connected component", "Yes", "No" + "``synthstrip``", "Yes", "No", "Human brain", "Yes", "Yes" "``trace``", "No", "Yes", "Single connected component", "No", "No" .. _dwi2mask_python: @@ -251,7 +300,7 @@ provided with one explicitly at the command-line) will internally execute the ``dwi2mask`` command. Because it is not possible for the user to manually specify how ``dwi2mask`` -should be utilised in this scenario, there are +should be utilised in this scenario, there are `configuration file options <../reference/config_file_options.html>`_ provided to assist in controlling the behaviour of ``dwi2mask`` in these scenarios (see below). @@ -259,14 +308,127 @@ scenarios (see below). .. csv-table:: :header: "*MRtrix3* Python command", "Purpose of DWI mask" :widths: auto - + + "``dwi2response``", "| Voxels outside of the initial mask are never considered as candidates for response function(s), nor do they contribute to any optimisation of the selection of such." "``dwibiascorrect``", "| Only voxels within the mask are utilised in optimisation of bias field parameters. | For ``ants`` algorithm, field is estimated within the mask but applied to all voxels within the field of view (field basis is extrapolated beyond the extremities of the mask); | for ``fsl`` algorithm, field is both estimated within, and applied to, only those voxels within the mask, producing a discontinuity in image intensity at the outer edge of the mask that can be deleterious for subsequent quantitative analyses." + "``dwibiasnormmask``", "| Determination of an *initial* brain mask by which to constrain the first iteration (see below)." "``dwifslpreproc``", "| Constrains optimisation of distortion parameter estimates in FSL_ ``eddy``. | If performing susceptibility distortion correction, this is applied to the DWI data subsequently to the appplication of FSL_ command ``applytopup``." "``dwigradcheck``", "| Utilised as both seed and mask image for streamlines tractography in the ``tckgen`` command." - "``dwi2response``", "| Voxels outside of the initial mask are never considered as candidates for response function(s), nor do they contribute to any optimisation of the selection of such." + +``dwibiasnormmask`` +------------------- + +This new script is an experimental approach for improving DWI brain mask estimation +(among other things), initially created during development of the MRtrix3_connectome_ +BIDS App. +It is based on the simple observation that the processes of bias field estimation, +intensity normalisation, and brain mask derivation, can have circular dependencies +between one another, and that therefore combining them into a single step may be +beneficial. +It is however noted that the behaviour of this algorithm can vary between different +types of data, and therefore close scruitiny of such is recommended. + +While this script is highly dependent on the operation of the ``mtnormalise`` command +(as was observed to be the case for the ``dwi2mask mtnorm`` algorithm above, +which performs a subset of the functionalities within ``dwibiasnormmask``), +the form of the primary results that it provides are slightly different: + +- *Output intended for usage*: + + With ``mtnormalise``, a set of ODFs are provided as input, and a set of ODFs + are then yielded as output, where the output ODFs have been corrected for + a smoothly-spatially-varying bias field, and global intensity scaling + (and importantly for quantitative applications the same intensity scaling is + applied to all ODFs). + For ``dwibiasnormmask``, the provided input is a DWI series, and the yielded + output is a DWI series, where the output has had the same smoothly-spatially-varying + bias field and global intensity scaling corrections applied. + The process of estimating these corrections is identical; + the only difference is that in ``dwibiasnormmask`` the corrections are *back-projected* + to correct the DWI series from which the ODFs were estimated, + rather than directly utilising the corrected ODFs. + +- *Global intensity normalisation*: + + The topic of :ref:`global-intensity-normalisation` is a long-standing issue + in the domain of CSD analysis. + Unlike other diffusion models, the *b*\=0 intensity of each voxel is not + used as a reference for the modelling of DWIs in that voxel, + and even in the context of multi-tissue CSD the composition of the voxel is + not explicitly forced to be unity. + This does however raise the issue of how to appropriately globally scale the + intensities of the image data in order for observed differences in eg. + Apparent Fibre Density (AFD) to be attributable to the effect of interest + rather then meaningless differential scaling of image intensities between subjects. + + The approach taken by ``mtnormalise`` is to determine the scaling factor that + results in voxels throughout the brain having a sum of tissue densities of + approximately unity. They will not all be exactly unity, even after bias field + correction, but they should be approximately centred around unity. + + ``dwibiasnormmask`` provides a more advanced version of the original proposal + for global intensity normalisation for AFD analysis. + It was first proposed that the *b*\=0 signal intensity in CSF should act as a + reference intensity to normalise across subjects. + However identifying appropriate exemplar voxels to do so can be labour-intensive + and difficult. + In ``dwibiasnormmask``, information from the ``dwi2response dhollander`` response + function estimation algorithm and the ``mtnormalise`` approach are combined in + such a way that the *b*\=0 CSF-like intensity is scaled to a fixed reference intensity: + + - using data from across the entire brain even in the presence of partial volume; + - accounting for potential miscalibrations in response function estimation; + - in conjunction with bias field estimation and correction. + + It is not yet known whether using this approach for global intensity normalisation + may yield greater sensitivity to effects of interest. + It should however be noted that if one were to subsequently execute ``mtnormalise`` + and make use of its output ODFs, + then the effective global intensity normalisation behaviour would revert to + that of ``mtnormalise`` rather than that described above. + +Another key aspect of this algorithm is the data used to derive the brain mask. +Most DWI brain masking approaches base their operation on the mean *b*\=0 image +(see dwi2mask_algorithm_comparison_ above). +In this algorithm, there is an alternative 3D image that can be used to drive brain mask +derivation, being the sum of tissue ODFs. +Depending on the configuration, this image may be used rather than the bias-field-corrected +DWI series to estimate a new brain mask at each iteration; +for instance, duplicating the functionality of the ``dwi2mask mtnorm`` algorithm above. + +The script itself operates as follows: + +1. If no initial mask is provided, then one must be calculated using ``dwi2mask``. + +2. Three-tissue response function estimation using ``dwi2response dhollander``. + +3. Multi-shell multi-tissue CSD, by default using a lower WM *lmax* for computational efficiency. + +4. ``mtnormalise`` to estimate bias field and intensity scaling factors between tissues. + +5. Estimation of a new brain mask, using either the bias-field-corrected DWI series + or the tissue ODF sum image. + +6. Determination of whether to exit, or *loop back to step 2*, based on: + + 1. Adequate similarity of the brain mask between successive iterations; + 2. Masks between successive iterations becoming less similar rather than more similar + (indicating some form of instability or divergence); + 3. Reaching maximal number of iterations. + +While the primary output of the command is the DWI series corrected for bias field and +intensity normalised, the brain mask corresponding to the last stable iteration can be +additionally exported using the ``-output_mask`` option. + +When initially developed, the number of iterations for this approach was fixed at 2, +as the solution was found to erroneously diverge after that in some instances. +It is however possible that for certain data, as well as the subsequent addition of the +explicit check for mask divergence, it is possible to permit a larger number of iterations +and allow the algorithm to converge toward a fixed solution. +Feedback on the success or failure of this experimental script for different data is encouraged. .. _dwi2mask_config: @@ -283,14 +445,14 @@ mentioned here also for discoverability: For those :ref:`dwi2mask_python`, this is the ``dwi2mask`` algorithm that will be invoked. If not explicitly set, the ``legacy`` algorithm will be used. - + .. NOTE:: - + Setting this configuration file option does *not* enable the utilisation of ``dwi2mask`` without manually specifying the algorithm to be used. For manual usage, the algorithm must *always* be specified. This option *only* controls the algorithm that will - be used when ``dwi2mask`` is invoked from inside one of the Python + be used when ``dwi2mask`` is invoked from inside one of the Python scripts provided with *MRtrix3*. - ``Dwi2maskTemplateSoftware`` @@ -305,7 +467,7 @@ mentioned here also for discoverability: - ``Dwi2maskTemplateImage`` and ``Dwi2maskTemplateMask`` This pair of configuration file options allow the user to pre-specify the - filesystem locations of the two images (T2-weighted template and + filesystem locations of the two images (T2-weighted template and corresponding binary mask) to be utilised by the ``dwi2mask ants`` and ``dwi2mask b02template`` algorithms. Note that there is no "default" template to be utilised by these algorithms; so the user *must* either include these @@ -330,6 +492,8 @@ mentioned here also for discoverability: .. _ANTs: http://stnava.github.io/ANTs/ .. _FSL: https://fsl.fmrib.ox.ac.uk/fsl/fslwiki .. _HD-BET: https://github.com/MIC-DKFZ/HD-BET +.. _MRtrix3_connectome: https://github.com/BIDS-Apps/MRtrix3_connectome +.. _SynthStrip: https://surfer.nmr.mgh.harvard.edu/docs/synthstrip/ .. _flirt: https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FLIRT .. _fnirt: https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FNIRT .. _invwarp: https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FNIRT/UserGuide#invwarp From 478f0e420dc330971d6fde6231e4a2a53a79324a Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Mon, 3 Apr 2023 12:53:52 +1000 Subject: [PATCH 18/24] Copyright headers for mtnormalise-related Python scripts --- bin/dwibiasnormmask | 2 +- lib/mrtrix3/dwi2mask/mtnorm.py | 2 +- lib/mrtrix3/dwibiascorrect/mtnorm.py | 2 +- lib/mrtrix3/dwinormalise/mtnorm.py | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/bin/dwibiasnormmask b/bin/dwibiasnormmask index c98d559dce..5baa6a4d64 100755 --- a/bin/dwibiasnormmask +++ b/bin/dwibiasnormmask @@ -1,6 +1,6 @@ #!/usr/bin/python3 -# Copyright (c) 2008-2022 the MRtrix3 contributors. +# Copyright (c) 2008-2023 the MRtrix3 contributors. # # This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this diff --git a/lib/mrtrix3/dwi2mask/mtnorm.py b/lib/mrtrix3/dwi2mask/mtnorm.py index 2c51515581..340e923723 100644 --- a/lib/mrtrix3/dwi2mask/mtnorm.py +++ b/lib/mrtrix3/dwi2mask/mtnorm.py @@ -1,4 +1,4 @@ -# Copyright (c) 2008-2022 the MRtrix3 contributors. +# Copyright (c) 2008-2023 the MRtrix3 contributors. # # This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this diff --git a/lib/mrtrix3/dwibiascorrect/mtnorm.py b/lib/mrtrix3/dwibiascorrect/mtnorm.py index 8d5eb2911e..dbbcd5350c 100644 --- a/lib/mrtrix3/dwibiascorrect/mtnorm.py +++ b/lib/mrtrix3/dwibiascorrect/mtnorm.py @@ -1,4 +1,4 @@ -# Copyright (c) 2008-2022 the MRtrix3 contributors. +# Copyright (c) 2008-2023 the MRtrix3 contributors. # # This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this diff --git a/lib/mrtrix3/dwinormalise/mtnorm.py b/lib/mrtrix3/dwinormalise/mtnorm.py index 69bec7b683..52bfa064c6 100644 --- a/lib/mrtrix3/dwinormalise/mtnorm.py +++ b/lib/mrtrix3/dwinormalise/mtnorm.py @@ -1,4 +1,4 @@ -# Copyright (c) 2008-2022 the MRtrix3 contributors. +# Copyright (c) 2008-2023 the MRtrix3 contributors. # # This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this From 06b4c4084d9e92d2ab5f5579ea5e67e25f11705a Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Mon, 3 Apr 2023 15:07:58 +1000 Subject: [PATCH 19/24] dwibiasnormmask: Multiple fixes - Add code to export bias field image (was in command-line interface but non-functional). - Fix erroneous deletion of tissue sum image prior to export. - Fix operation of -force option for -output_scale. --- bin/dwibiasnormmask | 43 +++++++++++++++++++++++++++---------------- 1 file changed, 27 insertions(+), 16 deletions(-) diff --git a/bin/dwibiasnormmask b/bin/dwibiasnormmask index 5baa6a4d64..d3a18e6dcd 100755 --- a/bin/dwibiasnormmask +++ b/bin/dwibiasnormmask @@ -222,6 +222,7 @@ def execute(): #pylint: disable=unused-variable dwi_image = 'input.mif' dwi_mask_image = 'dwi_mask_init.mif' + bias_field_image = None tissue_sum_image = None iteration = 0 step = 'initialisation' @@ -272,12 +273,13 @@ def execute(): #pylint: disable=unused-variable step = 'mtnormalise' progress.increment() - field_path = 'field' + iter_string + '.mif' + old_bias_field_image = bias_field_image + bias_field_image = 'field' + iter_string + '.mif' factors_path = 'factors' + iter_string + '.txt' run.command('mtnormalise -balanced' + ' -mask ' + eroded_mask - + ' -check_norm ' + field_path + + ' -check_norm ' + bias_field_image + ' -check_factors ' + factors_path + ' ' + ' '.join(tissue.fod_init + ' ' + tissue.fod_norm @@ -297,10 +299,9 @@ def execute(): #pylint: disable=unused-variable (csf_rf_bzero_lzero / csf_balance_factor) new_dwi_image = 'dwi' + iter_string + '.mif' run.command('mrcalc ' + dwi_image + ' ' - + field_path + ' -div ' + + bias_field_image + ' -div ' + str(scale_multiplier) + ' -mult ' + new_dwi_image) - app.cleanup(field_path) old_dwi_image = dwi_image dwi_image = new_dwi_image @@ -366,13 +367,12 @@ def execute(): #pylint: disable=unused-variable run.command('mrcalc ' + new_dwi_mask_image_nii + ' input_pos_mask.mif -mult ' + new_dwi_mask_image) else: assert False - app.cleanup(tissue_sum_image) if tissue_sum_image_nii: app.cleanup(tissue_sum_image_nii) if new_dwi_mask_image_nii: app.cleanup(new_dwi_mask_image_nii) - step = 'mask_comparison' + step = 'mask comparison' progress.increment() dwi_old_mask_count = image.statistics(dwi_mask_image, mask=dwi_mask_image).count @@ -393,11 +393,22 @@ def execute(): #pylint: disable=unused-variable new_dice_coefficient = 2.0 * mask_overlap_count / \ (dwi_old_mask_count + dwi_new_mask_count) + if iteration == app.ARGS.max_iters: + progress.done() + app.console('Terminating due to reaching maximum %d iterations; final Dice coefficient = %f' % (iteration, new_dice_coefficient)) + app.cleanup(old_dwi_image) + app.cleanup(old_dwi_mask_image) + app.cleanup(old_bias_field_image) + app.cleanup(old_tissue_sum_image) + total_scaling_factor *= scale_multiplier + break + if new_dice_coefficient > app.ARGS.dice: progress.done() app.console('Exiting loop after %d iterations due to mask convergence (Dice coefficient = %f)' % (iteration, new_dice_coefficient)) app.cleanup(old_dwi_image) app.cleanup(old_dwi_mask_image) + app.cleanup(old_bias_field_image) app.cleanup(old_tissue_sum_image) total_scaling_factor *= scale_multiplier break @@ -408,21 +419,14 @@ def execute(): #pylint: disable=unused-variable + ' using mask from previous iteration') app.cleanup(dwi_image) app.cleanup(dwi_mask_image) + app.cleanup(bias_field_image) app.cleanup(tissue_sum_image) dwi_image = old_dwi_image dwi_mask_image = old_dwi_mask_image + bias_field_image = old_bias_field_image tissue_sum_image = old_tissue_sum_image break - if iteration == app.ARGS.max_iters: - progress.done() - app.console('Terminating due to reaching maximum %d iterations; final Dice coefficient = %f' % (iteration, new_dice_coefficient)) - app.cleanup(old_dwi_image) - app.cleanup(old_dwi_mask_image) - app.cleanup(old_tissue_sum_image) - total_scaling_factor *= scale_multiplier - break - iteration += 1 prev_dice_coefficient = new_dice_coefficient @@ -431,13 +435,20 @@ def execute(): #pylint: disable=unused-variable mrconvert_keyval=path.from_user(app.ARGS.input, False), force=app.FORCE_OVERWRITE) + if app.ARGS.output_bias: + run.command(['mrconvert', bias_field_image, path.from_user(app.ARGS.output_bias, False)], + mrconvert_keyval=path.from_user(app.ARGS.input, False), + force=app.FORCE_OVERWRITE) + if app.ARGS.output_mask: run.command(['mrconvert', dwi_mask_image, path.from_user(app.ARGS.output_mask, False)], mrconvert_keyval=path.from_user(app.ARGS.input, False), force=app.FORCE_OVERWRITE) if app.ARGS.output_scale: - matrix.save_vector(path.from_user(app.ARGS.output_scale, False), [total_scaling_factor]) + matrix.save_vector(path.from_user(app.ARGS.output_scale, False), + [total_scaling_factor], + force=app.FORCE_OVERWRITE) if app.ARGS.output_tissuesum: run.command(['mrconvert', tissue_sum_image, path.from_user(app.ARGS.output_tissuesum, False)], From 19698ddf0d71a2c1103f9ec06b134359b9c642f8 Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Mon, 3 Apr 2023 15:09:43 +1000 Subject: [PATCH 20/24] dwibiasnormmask: Fix testing and add test data --- testing/scripts/data | 2 +- testing/scripts/tests/dwibiasnormmask | 9 +++++---- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/testing/scripts/data b/testing/scripts/data index 6c6bb12ecb..367e2ddafd 160000 --- a/testing/scripts/data +++ b/testing/scripts/data @@ -1 +1 @@ -Subproject commit 6c6bb12ecb6cb356dec37b97bfa68d2a1ad952cd +Subproject commit 367e2ddafdefa18ab4c63668cbb3c625e9c2d0ea diff --git a/testing/scripts/tests/dwibiasnormmask b/testing/scripts/tests/dwibiasnormmask index 5c5c2e957a..d3936e70ee 100644 --- a/testing/scripts/tests/dwibiasnormmask +++ b/testing/scripts/tests/dwibiasnormmask @@ -1,4 +1,5 @@ -mkdir -p ../tmp/dwibiasnormmask && mrconvert BIDS/sub-01/dwi/sub-01_dwi.nii.gz -fslgrad BIDS/sub-01/dwi/sub-01_dwi.bvec BIDS/sub-01/dwi/sub-01_dwi.bval tmp-sub-01_dwi.mif -force && dwibiasnormmask tmp-sub-01_dwi.mif ../tmp/dwibiasnormmask/default_out.mif -output_mask ../tmp/dwibiasnormmask/default_mask.mif -force -dwibiasnormmask tmp-sub-01_dwi.mif -max_iters 3 ../tmp/dwibiasnormmask/3iters_out.mif -output_mask ../tmp/dwibiasnormmask/3iters_mask.mif -force -dwibiasnormmask tmp-sub-01_dwi.mif -lmax 6,0,0 ../tmp/dwibiasnormmask/lmax600_out.mif -output_mask ../tmp/dwibiasnormmask/lmax600_mask.mif -force - +mkdir -p ../tmp/dwibiasnormmask && mrconvert BIDS/sub-01/dwi/sub-01_dwi.nii.gz -fslgrad BIDS/sub-01/dwi/sub-01_dwi.bvec BIDS/sub-01/dwi/sub-01_dwi.bval tmp-sub-01_dwi.mif -force && dwibiasnormmask tmp-sub-01_dwi.mif ../tmp/dwibiasnormmask/default_out.mif ../tmp/dwibiasnormmask/default_mask.mif -output_bias ../tmp/dwibiasnormmask/default_bias.mif -output_scale ../tmp/dwibiasnormmask/default_scale.txt -output_tissuesum ../tmp/dwibiasnormmask/default_tissuesum.mif -force && testing_diff_image ../tmp/dwibiasnormmask/default_out.mif dwibiasnormmask/default_out.mif.gz -frac 1e-5 && testing_diff_image ../tmp/dwibiasnormmask/default_mask.mif dwibiasnormmask/default_mask.mif.gz && testing_diff_image ../tmp/dwibiasnormmask/default_bias.mif dwibiasnormmask/default_bias.mif.gz -frac 1e-5 && testing_diff_matrix ../tmp/dwibiasnormmask/default_scale.txt dwibiasnormmask/default_scale.txt -frac 1e-5 && testing_diff_image ../tmp/dwibiasnormmask/default_tissuesum.mif dwibiasnormmask/default_tissuesum.mif.gz -abs 1e-5 +dwibiasnormmask tmp-sub-01_dwi.mif -max_iters 3 -dice 1.0 ../tmp/dwibiasnormmask/3iters_out.mif ../tmp/dwibiasnormmask/3iters_mask.mif -output_bias ../tmp/dwibiasnormmask/3iters_bias.mif -output_scale ../tmp/dwibiasnormmask/3iters_scale.txt -output_tissuesum ../tmp/dwibiasnormmask/3iters_tissuesum.mif -force && testing_diff_image ../tmp/dwibiasnormmask/3iters_out.mif dwibiasnormmask/3iters_out.mif.gz -frac 1e-5 && testing_diff_image ../tmp/dwibiasnormmask/3iters_mask.mif dwibiasnormmask/3iters_mask.mif.gz && testing_diff_image ../tmp/dwibiasnormmask/3iters_bias.mif dwibiasnormmask/3iters_bias.mif.gz -frac 1e-5 && testing_diff_matrix ../tmp/dwibiasnormmask/3iters_scale.txt dwibiasnormmask/3iters_scale.txt -frac 1e-5 && testing_diff_image ../tmp/dwibiasnormmask/3iters_tissuesum.mif dwibiasnormmask/3iters_tissuesum.mif.gz -abs 1e-5 +dwibiasnormmask tmp-sub-01_dwi.mif -lmax 6,0,0 ../tmp/dwibiasnormmask/lmax600_out.mif ../tmp/dwibiasnormmask/lmax600_mask.mif -output_bias ../tmp/dwibiasnormmask/lmax600_bias.mif -output_scale ../tmp/dwibiasnormmask/lmax600_scale.txt -output_tissuesum ../tmp/dwibiasnormmask/lmax600_tissuesum.mif -force && testing_diff_image ../tmp/dwibiasnormmask/lmax600_out.mif dwibiasnormmask/lmax600_out.mif.gz -frac 1e-5 && testing_diff_image ../tmp/dwibiasnormmask/lmax600_mask.mif dwibiasnormmask/lmax600_mask.mif.gz && testing_diff_image ../tmp/dwibiasnormmask/lmax600_bias.mif dwibiasnormmask/lmax600_bias.mif.gz -frac 1e-5 && testing_diff_matrix ../tmp/dwibiasnormmask/lmax600_scale.txt dwibiasnormmask/lmax600_scale.txt -frac 1e-5 && testing_diff_image ../tmp/dwibiasnormmask/lmax600_tissuesum.mif dwibiasnormmask/lmax600_tissuesum.mif.gz -abs 1e-5 +dwibiasnormmask tmp-sub-01_dwi.mif -reference 1.0 ../tmp/dwibiasnormmask/reference_out.mif ../tmp/dwibiasnormmask/reference_mask.mif -output_scale ../tmp/dwibiasnormmask/scale.txt -force && mrcalc ../tmp/dwibiasnormmask/reference_out.mif 1000.0 -mult - | testing_diff_image - dwibiasnormmask/default_out.mif.gz -frac 1e-5 +mrcalc tmp-sub-01_dwi.mif 1000.0 -div tmp-sub-01_dwi_scaled.mif -force && dwibiasnormmask tmp-sub-01_dwi_scaled.mif ../tmp/dwibiasnormmask/scaled_out.mif ../tmp/dwibiasnormmask/scaled_mask.mif -output_scale ../tmp/dwibiasnormmask/scaled_scale.txt -force && testing_diff_image ../tmp/dwibiasnormmask/scaled_out.mif dwibiasnormmask/default_out.mif.gz -frac 1e-5 From 151473843ebc0553ff028aaee759bb3ae15454dc Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Mon, 3 Apr 2023 15:57:06 +1000 Subject: [PATCH 21/24] dwinormalise group: Reflect renaming of dwinormalise individual -> manual --- lib/mrtrix3/dwinormalise/group.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/mrtrix3/dwinormalise/group.py b/lib/mrtrix3/dwinormalise/group.py index 3911833978..77260dfbc2 100644 --- a/lib/mrtrix3/dwinormalise/group.py +++ b/lib/mrtrix3/dwinormalise/group.py @@ -111,7 +111,7 @@ def __init__(self, filename, prefix, mask_filename = ''): path.make_dir('wm_mask_warped') for i in input_list: run.command('mrtransform template_wm_mask.mif -interp nearest -warp_full ' + os.path.join('warps', i.prefix + '.mif') + ' ' + os.path.join('wm_mask_warped', i.prefix + '.mif') + ' -from 2 -template ' + os.path.join('fa', i.prefix + '.mif')) - run.command('dwinormalise individual ' + shlex.quote(os.path.join(input_dir, i.filename)) + ' ' + os.path.join('wm_mask_warped', i.prefix + '.mif') + ' temp.mif') + run.command('dwinormalise manual ' + shlex.quote(os.path.join(input_dir, i.filename)) + ' ' + os.path.join('wm_mask_warped', i.prefix + '.mif') + ' temp.mif') run.command('mrconvert temp.mif ' + path.from_user(os.path.join(app.ARGS.output_dir, i.filename)), mrconvert_keyval=path.from_user(os.path.join(input_dir, i.filename), False), force=app.FORCE_OVERWRITE) os.remove('temp.mif') progress.increment() From 79b358c9d80239146db6406ac69819193804e917 Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Mon, 3 Apr 2023 15:57:28 +1000 Subject: [PATCH 22/24] dwinormalise mtnorm: Fix operation of -reference option --- lib/mrtrix3/dwinormalise/mtnorm.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/mrtrix3/dwinormalise/mtnorm.py b/lib/mrtrix3/dwinormalise/mtnorm.py index 52bfa064c6..7f9c5a7589 100644 --- a/lib/mrtrix3/dwinormalise/mtnorm.py +++ b/lib/mrtrix3/dwinormalise/mtnorm.py @@ -164,7 +164,7 @@ def __init__(self, name): balance_factors = matrix.load_vector('factors.txt') app.cleanup('factors.txt') csf_balance_factor = balance_factors[-1] - scale_multiplier = (1000.0 * math.sqrt(4.0*math.pi)) / (csf_rf_bzero_lzero / csf_balance_factor) + scale_multiplier = (app.ARGS.reference * math.sqrt(4.0*math.pi)) / (csf_rf_bzero_lzero / csf_balance_factor) run.command('mrcalc input.mif ' + str(scale_multiplier) + ' -mult - | ' + 'mrconvert - ' + path.from_user(app.ARGS.output), From 528d0c5850de370ebd2b1c2d3259f1b8278e34ef Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Mon, 3 Apr 2023 16:02:55 +1000 Subject: [PATCH 23/24] New tests and data for mtnormalise-related Python script algorithms --- testing/scripts/data | 2 +- testing/scripts/tests/dwi2mask | 3 +++ testing/scripts/tests/dwibiascorrect | 3 +++ testing/scripts/tests/dwinormalise | 9 +++++++-- 4 files changed, 14 insertions(+), 3 deletions(-) diff --git a/testing/scripts/data b/testing/scripts/data index 367e2ddafd..ad1cf773e4 160000 --- a/testing/scripts/data +++ b/testing/scripts/data @@ -1 +1 @@ -Subproject commit 367e2ddafdefa18ab4c63668cbb3c625e9c2d0ea +Subproject commit ad1cf773e49d89ec23c0e2a064672b4a3aa59424 diff --git a/testing/scripts/tests/dwi2mask b/testing/scripts/tests/dwi2mask index 99771a1ee5..9f7cdf3e0a 100644 --- a/testing/scripts/tests/dwi2mask +++ b/testing/scripts/tests/dwi2mask @@ -22,6 +22,9 @@ dwi2mask fslbet tmp-sub-01_dwi.mif ../tmp/dwi2mask/fslbet_rescale.mif -rescale - dwi2mask hdbet tmp-sub-01_dwi.mif ../tmp/dwi2mask/hdbet.mif -force dwi2mask legacy tmp-sub-01_dwi.mif ../tmp/dwi2mask/legacy.mif -force && testing_diff_image ../tmp/dwi2mask/legacy.mif dwi2mask/legacy.mif.gz dwi2mask mean tmp-sub-01_dwi.mif ../tmp/dwi2mask/mean.mif -force && testing_diff_image ../tmp/dwi2mask/mean.mif dwi2mask/mean.mif.gz +dwi2mask mtnorm tmp-sub-01_dwi.mif ../tmp/dwi2mask/mtnorm_default_mask.mif -tissuesum ../tmp/dwi2mask/mtnorm_default_tissuesum.mif -force && testing_diff_image ../tmp/dwi2mask/mtnorm_default_mask.mif dwi2mask/mtnorm_default_mask.mif.gz && testing_diff_image ../tmp/dwi2mask/mtnorm_default_tissuesum.mif dwi2mask/mtnorm_default_tissuesum.mif.gz -abs 1e-5 +dwi2mask mtnorm tmp-sub-01_dwi.mif -init_mask BIDS/sub-01/dwi/sub-01_brainmask.nii.gz ../tmp/dwi2mask/mtnorm_initmask_mask.mif -tissuesum ../tmp/dwi2mask/mtnorm_initmask_tissuesum.mif -force && testing_diff_image ../tmp/dwi2mask/mtnorm_initmask_mask.mif dwi2mask/mtnorm_initmask_mask.mif.gz && testing_diff_image ../tmp/dwi2mask/mtnorm_initmask_tissuesum.mif dwi2mask/mtnorm_initmask_tissuesum.mif.gz -abs 1e-5 +dwi2mask mtnorm tmp-sub-01_dwi.mif -lmax 6,0,0 ../tmp/dwi2mask/mtnorm_lmax600_mask.mif -tissuesum ../tmp/dwi2mask/mtnorm_lmax600_tissuesum.mif -force && testing_diff_image ../tmp/dwi2mask/mtnorm_lmax600_mask.mif dwi2mask/mtnorm_lmax600_mask.mif.gz && testing_diff_image ../tmp/dwi2mask/mtnorm_lmax600_tissuesum.mif dwi2mask/mtnorm_lmax600_tissuesum.mif.gz -abs 1e-5 dwi2mask synthstrip tmp-sub-01_dwi.mif ../tmp/dwi2mask/synthstrip_default.mif -force dwi2mask synthstrip tmp-sub-01_dwi.mif ../tmp/dwi2mask/synthstrip_options.mif -stripped ../tmp/dwi2mask/synthstrip_stripped.mif -gpu -nocsf -border 0 -force dwi2mask trace tmp-sub-01_dwi.mif ../tmp/dwi2mask/trace_default.mif -force && testing_diff_image ../tmp/dwi2mask/trace_default.mif dwi2mask/trace_default.mif.gz diff --git a/testing/scripts/tests/dwibiascorrect b/testing/scripts/tests/dwibiascorrect index 9dd108d370..c75f2f6632 100644 --- a/testing/scripts/tests/dwibiascorrect +++ b/testing/scripts/tests/dwibiascorrect @@ -6,3 +6,6 @@ mkdir -p ../tmp/dwibiascorrect/fsl && dwibiascorrect fsl tmp-sub-01_dwi.mif ../t dwibiascorrect fsl BIDS/sub-01/dwi/sub-01_dwi.nii.gz -fslgrad BIDS/sub-01/dwi/sub-01_dwi.bvec BIDS/sub-01/dwi/sub-01_dwi.bval ../tmp/dwibiascorrect/fsl/fslgrad_out.mif -bias ../tmp/dwibiascorrect/fsl/fslgrad_bias.mif -force && testing_diff_header ../tmp/dwibiascorrect/fsl/fslgrad_out.mif dwibiascorrect/fsl/default_out.mif.gz && testing_diff_header ../tmp/dwibiascorrect/fsl/fslgrad_bias.mif dwibiascorrect/fsl/default_bias.mif.gz dwibiascorrect fsl BIDS/sub-01/dwi/sub-01_dwi.nii.gz -grad tmp-sub-01_dwi.b ../tmp/dwibiascorrect/fsl/grad_out.mif -bias ../tmp/dwibiascorrect/fsl/grad_bias.mif -force && testing_diff_header ../tmp/dwibiascorrect/fsl/grad_out.mif dwibiascorrect/fsl/default_out.mif.gz && testing_diff_header ../tmp/dwibiascorrect/fsl/grad_bias.mif dwibiascorrect/fsl/default_bias.mif.gz dwibiascorrect fsl tmp-sub-01_dwi.mif ../tmp/dwibiascorrect/fsl/masked_out.mif -bias ../tmp/dwibiascorrect/fsl/masked_bias.mif -mask BIDS/sub-01/dwi/sub-01_brainmask.nii.gz -force && testing_diff_header ../tmp/dwibiascorrect/fsl/masked_out.mif dwibiascorrect/fsl/masked_out.mif.gz && testing_diff_header ../tmp/dwibiascorrect/fsl/masked_bias.mif dwibiascorrect/fsl/masked_bias.mif.gz +mkdir -p ../tmp/dwibiascorrect/mtnorm && dwibiascorrect mtnorm tmp-sub-01_dwi.mif ../tmp/dwibiascorrect/mtnorm/default_out.mif -bias ../tmp/dwibiascorrect/mtnorm/default_bias.mif -force && testing_diff_image ../tmp/dwibiascorrect/mtnorm/default_out.mif dwibiascorrect/mtnorm/default_out.mif.gz -frac 1e-5 && testing_diff_image ../tmp/dwibiascorrect/mtnorm/default_bias.mif dwibiascorrect/mtnorm/default_bias.mif.gz -frac 1e-5 +dwibiascorrect mtnorm tmp-sub-01_dwi.mif ../tmp/dwibiascorrect/mtnorm/masked_out.mif -bias ../tmp/dwibiascorrect/mtnorm/masked_bias.mif -mask BIDS/sub-01/dwi/sub-01_brainmask.nii.gz -force && testing_diff_image ../tmp/dwibiascorrect/mtnorm/masked_out.mif dwibiascorrect/mtnorm/masked_out.mif.gz && testing_diff_image ../tmp/dwibiascorrect/mtnorm/masked_bias.mif dwibiascorrect/mtnorm/masked_bias.mif.gz +dwibiascorrect mtnorm tmp-sub-01_dwi.mif -lmax 6,0,0 ../tmp/dwibiascorrect/mtnorm/lmax600_out.mif -bias ../tmp/dwibiascorrect/mtnorm/lmax600_bias.mif -force && testing_diff_image ../tmp/dwibiascorrect/mtnorm/lmax600_out.mif dwibiascorrect/mtnorm/lmax600_out.mif.gz && testing_diff_image ../tmp/dwibiascorrect/mtnorm/lmax600_bias.mif dwibiascorrect/mtnorm/lmax600_bias.mif.gz diff --git a/testing/scripts/tests/dwinormalise b/testing/scripts/tests/dwinormalise index b3e46b60e2..1b66ebfd0e 100644 --- a/testing/scripts/tests/dwinormalise +++ b/testing/scripts/tests/dwinormalise @@ -1,3 +1,8 @@ mkdir -p ../tmp/dwinormalise/group && mkdir -p tmp-dwi && mkdir -p tmp-mask && mrconvert BIDS/sub-02/dwi/sub-02_dwi.nii.gz -fslgrad BIDS/sub-02/dwi/sub-02_dwi.bvec BIDS/sub-02/dwi/sub-02_dwi.bval tmp-dwi/sub-02.mif -force && mrconvert BIDS/sub-02/dwi/sub-02_brainmask.nii.gz tmp-mask/sub-02.mif -force && mrconvert BIDS/sub-03/dwi/sub-03_dwi.nii.gz -fslgrad BIDS/sub-03/dwi/sub-03_dwi.bvec BIDS/sub-03/dwi/sub-03_dwi.bval tmp-dwi/sub-03.mif -force && mrconvert BIDS/sub-03/dwi/sub-03_brainmask.nii.gz tmp-mask/sub-03.mif -force && dwinormalise group tmp-dwi/ tmp-mask/ ../tmp/dwinormalise/group/ ../tmp/dwinormalise/group/fa.mif ../tmp/dwinormalise/group/mask.mif -force && testing_diff_image ../tmp/dwinormalise/group/sub-02.mif dwinormalise/group/sub-02.mif.gz -frac 1e-2 && testing_diff_image ../tmp/dwinormalise/group/sub-03.mif dwinormalise/group/sub-03.mif.gz -frac 1e-2 && testing_diff_image ../tmp/dwinormalise/group/fa.mif dwinormalise/group/fa.mif.gz -abs 1e-3 && testing_diff_image $(mrfilter ../tmp/dwinormalise/group/mask.mif smooth -) $(mrfilter dwinormalise/group/mask.mif.gz smooth -) -abs 0.3 -mkdir ../tmp/dwinormalise/individual && dwinormalise individual BIDS/sub-01/dwi/sub-01_dwi.nii.gz -fslgrad BIDS/sub-01/dwi/sub-01_dwi.bvec BIDS/sub-01/dwi/sub-01_dwi.bval BIDS/sub-01/dwi/sub-01_brainmask.nii.gz ../tmp/dwinormalise/individual/out.mif -force && testing_diff_image ../tmp/dwinormalise/individual/out.mif dwinormalise/individual/out.mif.gz -frac 1e-5 -dwinormalise individual BIDS/sub-01/dwi/sub-01_dwi.nii.gz -fslgrad BIDS/sub-01/dwi/sub-01_dwi.bvec BIDS/sub-01/dwi/sub-01_dwi.bval BIDS/sub-01/dwi/sub-01_brainmask.nii.gz ../tmp/dwinormalise/individual/percentile.mif -percentile 40 -force && testing_diff_image ../tmp/dwinormalise/individual/percentile.mif dwinormalise/individual/percentile.mif.gz -frac 1e-5 +mkdir ../tmp/dwinormalise/manual && dwinormalise manual BIDS/sub-01/dwi/sub-01_dwi.nii.gz -fslgrad BIDS/sub-01/dwi/sub-01_dwi.bvec BIDS/sub-01/dwi/sub-01_dwi.bval BIDS/sub-01/dwi/sub-01_brainmask.nii.gz ../tmp/dwinormalise/manual/out.mif -force && testing_diff_image ../tmp/dwinormalise/manual/out.mif dwinormalise/manual/out.mif.gz -frac 1e-5 +dwinormalise manual BIDS/sub-01/dwi/sub-01_dwi.nii.gz -fslgrad BIDS/sub-01/dwi/sub-01_dwi.bvec BIDS/sub-01/dwi/sub-01_dwi.bval BIDS/sub-01/dwi/sub-01_brainmask.nii.gz ../tmp/dwinormalise/manual/percentile.mif -percentile 40 -force && testing_diff_image ../tmp/dwinormalise/manual/percentile.mif dwinormalise/manual/percentile.mif.gz -frac 1e-5 +mkdir ../tmp/dwinormalise/mtnorm && dwinormalise mtnorm BIDS/sub-01/dwi/sub-01_dwi.nii.gz -fslgrad BIDS/sub-01/dwi/sub-01_dwi.bvec BIDS/sub-01/dwi/sub-01_dwi.bval ../tmp/dwinormalise/mtnorm/default_out.mif -scale ../tmp/dwinormalise/mtnorm/default_scale.txt -force && testing_diff_image ../tmp/dwinormalise/mtnorm/default_out.mif dwinormalise/mtnorm/default_out.mif.gz -frac 1e-5 && testing_diff_matrix ../tmp/dwinormalise/mtnorm/default_scale.txt dwinormalise/mtnorm/default_scale.txt -abs 1e-5 +dwinormalise mtnorm BIDS/sub-01/dwi/sub-01_dwi.nii.gz -mask BIDS/sub-01/dwi/sub-01_brainmask.nii.gz -fslgrad BIDS/sub-01/dwi/sub-01_dwi.bvec BIDS/sub-01/dwi/sub-01_dwi.bval ../tmp/dwinormalise/mtnorm/masked_out.mif -scale ../tmp/dwinormalise/mtnorm/masked_scale.txt -force && testing_diff_image ../tmp/dwinormalise/mtnorm/masked_out.mif dwinormalise/mtnorm/masked_out.mif.gz -frac 1e-5 && testing_diff_matrix ../tmp/dwinormalise/mtnorm/masked_scale.txt dwinormalise/mtnorm/masked_scale.txt -abs 1e-5 +dwinormalise mtnorm BIDS/sub-01/dwi/sub-01_dwi.nii.gz -lmax 6,0,0 -fslgrad BIDS/sub-01/dwi/sub-01_dwi.bvec BIDS/sub-01/dwi/sub-01_dwi.bval ../tmp/dwinormalise/mtnorm/lmax600_out.mif -scale ../tmp/dwinormalise/mtnorm/lmax600_scale.txt -force && testing_diff_image ../tmp/dwinormalise/mtnorm/lmax600_out.mif dwinormalise/mtnorm/lmax600_out.mif.gz -frac 1e-5 && testing_diff_matrix ../tmp/dwinormalise/mtnorm/lmax600_scale.txt dwinormalise/mtnorm/lmax600_scale.txt -abs 1e-5 +mrcalc BIDS/sub-01/dwi/sub-01_dwi.nii.gz 1000.0 -div tmp-sub-01_dwi_scaled.mif -force && dwinormalise mtnorm tmp-sub-01_dwi_scaled.mif -fslgrad BIDS/sub-01/dwi/sub-01_dwi.bvec BIDS/sub-01/dwi/sub-01_dwi.bval ../tmp/dwinormalise/mtnorm/scaled_out.mif -scale ../tmp/dwinormalise/mtnorm/scaled_scale.txt -force && testing_diff_image ../tmp/dwinormalise/mtnorm/scaled_out.mif dwinormalise/mtnorm/default_out.mif.gz -frac 1e-5 +dwinormalise mtnorm BIDS/sub-01/dwi/sub-01_dwi.nii.gz -reference 1.0 -fslgrad BIDS/sub-01/dwi/sub-01_dwi.bvec BIDS/sub-01/dwi/sub-01_dwi.bval ../tmp/dwinormalise/mtnorm/reference_out.mif -scale ../tmp/dwinormalise/mtnorm/reference_scale.txt -force && mrcalc ../tmp/dwinormalise/mtnorm/reference_out.mif 1000.0 -mult - | testing_diff_image - dwinormalise/mtnorm/default_out.mif.gz -frac 1e-5 From 86d6d8976fffe4a1e073a00c395a39218060374e Mon Sep 17 00:00:00 2001 From: Daljit Singh Date: Wed, 21 Jun 2023 08:52:15 +0100 Subject: [PATCH 24/24] Add support for multiple arguments to CXX_COMPILER_LAUNCHER --- configure | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/configure b/configure index 2b1143b79c..5af4790b48 100755 --- a/configure +++ b/configure @@ -1481,8 +1481,9 @@ with open (config_filename, 'w', encoding='utf-8') as config_file: commit (config_file, 'lib_suffix', lib_suffix) if "CXX_COMPILER_LAUNCHER" in os.environ: - print("\n USING COMPILER PREFIX", os.environ["CXX_COMPILER_LAUNCHER"]) - cpp.insert(0, os.environ["CXX_COMPILER_LAUNCHER"]) + cpp_launcher = os.environ["CXX_COMPILER_LAUNCHER"].split() + print("\n USING CXX_COMPILER_LAUNCHER", cpp_launcher) + cpp = cpp_launcher + cpp commit (config_file, 'cpp', cpp) commit (config_file, 'cpp_flags', cpp_flags)