Merge branch 'dev' into clangformat

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')

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()

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 = ''

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)

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)

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

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`,

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,22 +300,135 @@ 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).
 
 .. 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