From 915b1857a27845aa3cec3a88b52e91c29c358207 Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Sun, 17 Nov 2024 22:09:34 +1100 Subject: [PATCH] New command dwi2noise Closes #3035. --- cmd/dwi2noise.cpp | 194 +++++++++++++++++++++ cmd/dwidenoise.cpp | 13 +- docs/reference/commands/dwi2noise.rst | 123 +++++++++++++ docs/reference/commands/dwidenoise.rst | 68 ++++++-- docs/reference/commands_list.rst | 2 + src/denoise/estimate.cpp | 120 +++++++++++++ src/denoise/{functor.h => estimate.h} | 47 ++--- src/denoise/functor.cpp | 230 ------------------------- src/denoise/recon.cpp | 163 ++++++++++++++++++ src/denoise/recon.h | 67 +++++++ 10 files changed, 749 insertions(+), 278 deletions(-) create mode 100644 cmd/dwi2noise.cpp create mode 100644 docs/reference/commands/dwi2noise.rst create mode 100644 src/denoise/estimate.cpp rename src/denoise/{functor.h => estimate.h} (57%) delete mode 100644 src/denoise/functor.cpp create mode 100644 src/denoise/recon.cpp create mode 100644 src/denoise/recon.h diff --git a/cmd/dwi2noise.cpp b/cmd/dwi2noise.cpp new file mode 100644 index 0000000000..045dd5dd2a --- /dev/null +++ b/cmd/dwi2noise.cpp @@ -0,0 +1,194 @@ +/* Copyright (c) 2008-2024 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/. + */ + +#include + +#include "algo/threaded_loop.h" +#include "command.h" +#include "denoise/estimate.h" +#include "denoise/estimator/estimator.h" +#include "denoise/exports.h" +#include "denoise/kernel/kernel.h" +#include "exception.h" + +using namespace MR; +using namespace App; +using namespace MR::Denoise; + +// clang-format off +void usage() { + + SYNOPSIS = "Noise level estimation using Marchenko-Pastur PCA"; + + DESCRIPTION + + "DWI data noise map estimation" + " by interrogating data redundancy in the PCA domain" + " using the prior knowledge that the eigenspectrum of random covariance matrices" + " is described by the universal Marchenko-Pastur (MP) distribution." + " Fitting the MP distribution to the spectrum of patch-wise signal matrices" + " hence provides an estimator of the noise level 'sigma'." + + + "Unlike the MRtrix3 command dwidenoise," + " this command does not generate a denoised version of the input image series;" + " its primary output is instead a map of the estimated noise level." + " While this can also be obtained from the dwidenoise command using option -noise_out," + " using instead the dwi2noise command gives the ability to obtain a noise map" + " to which filtering can be applied," + " which can then be utilised for the actual image series denoising," + " without generating an unwanted intermiedate denoised image series." + + + "Important note:" + " noise level estimation should only be performed as the first step of an image processing pipeline." + " The routine is invalid if interpolation or smoothing has been applied to the data prior to denoising." + + + "Note that on complex input data," + " the output will be the total noise level across real and imaginary channels," + " so a scale factor sqrt(2) applies." + + + Kernel::shape_description + + + Kernel::size_description; + + AUTHOR = "Daan Christiaens (daan.christiaens@kcl.ac.uk)" + " and Jelle Veraart (jelle.veraart@nyumc.org)" + " and J-Donald Tournier (jdtournier@gmail.com)" + " and Robert E. Smith (robert.smith@florey.edu.au)"; + + REFERENCES + + "Veraart, J.; Fieremans, E. & Novikov, D.S. " // Internal + "Diffusion MRI noise mapping using random matrix theory. " + "Magn. Res. Med., 2016, 76(5), 1582-1593, doi: 10.1002/mrm.26059" + + + "Cordero-Grande, L.; Christiaens, D.; Hutter, J.; Price, A.N.; Hajnal, J.V. " // Internal + "Complex diffusion-weighted image estimation via matrix recovery under general noise models. " + "NeuroImage, 2019, 200, 391-404, doi: 10.1016/j.neuroimage.2019.06.039" + + + "* If using -estimator mrm2022: " + "Olesen, J.L.; Ianus, A.; Ostergaard, L.; Shemesh, N.; Jespersen, S.N. " + "Tensor denoising of multidimensional MRI data. " + "Magnetic Resonance in Medicine, 2022, 89(3), 1160-1172"; + + ARGUMENTS + + Argument("dwi", "the input diffusion-weighted image").type_image_in() + + Argument("noise", "the output estimated noise level map").type_image_out(); + + OPTIONS + + OptionGroup("Options for modifying PCA computations") + + datatype_option + + Estimator::option + + Kernel::options + + // TODO Implement mask option + // Note that behaviour of -mask for dwi2noise may be different to that of dwidenoise + + + OptionGroup("Options for exporting additional data regarding PCA behaviour") + + Option("rank", + "The signal rank estimated for the denoising patch centred at each input image voxel") + + Argument("image").type_image_out() + + OptionGroup("Options for debugging the operation of sliding window kernels") + + Option("max_dist", + "The maximum distance between a voxel and another voxel that was included in the local denoising patch") + + Argument("image").type_image_out() + + Option("voxelcount", + "The number of voxels that contributed to the PCA for processing of each voxel") + + Argument("image").type_image_out(); + + COPYRIGHT = + "Copyright (c) 2016 New York University, University of Antwerp, and the MRtrix3 contributors \n \n" + "Permission is hereby granted, free of charge, to any non-commercial entity ('Recipient') obtaining a copy of " + "this software and " + "associated documentation files (the 'Software'), to the Software solely for non-commercial research, including " + "the rights to " + "use, copy and modify the Software, subject to the following conditions: \n \n" + "\t 1. The above copyright notice and this permission notice shall be included by Recipient in all copies or " + "substantial portions of " + "the Software. \n \n" + "\t 2. THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT " + "LIMITED TO THE WARRANTIES" + "OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR " + "COPYRIGHT HOLDERS BE" + "LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING " + "FROM, OUT OF OR" + "IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. \n \n" + "\t 3. In no event shall NYU be liable for direct, indirect, special, incidental or consequential damages in " + "connection with the Software. " + "Recipient will defend, indemnify and hold NYU harmless from any claims or liability resulting from the use of " + "the Software by recipient. \n \n" + "\t 4. Neither anything contained herein nor the delivery of the Software to recipient shall be deemed to grant " + "the Recipient any right or " + "licenses under any patents or patent application owned by NYU. \n \n" + "\t 5. The Software may only be used for non-commercial research and may not be used for clinical care. \n \n" + "\t 6. Any publication by Recipient of research involving the Software shall cite the references listed below."; +} +// clang-format on + +template +void run(Header &data, + std::shared_ptr kernel, + std::shared_ptr estimator, + Exports &exports) { + auto input = data.get_image().with_direct_io(3); + Image mask; // unused + Estimate func(data, mask, kernel, estimator, exports); + ThreadedLoop("running MP-PCA noise level estimation", data, 0, 3).run(func, input); +} + +void run() { + auto dwi = Header::open(argument[0]); + + if (dwi.ndim() != 4 || dwi.size(3) <= 1) + throw Exception("input image must be 4-dimensional"); + + auto kernel = Kernel::make_kernel(dwi); + assert(kernel); + + auto estimator = Estimator::make_estimator(); + assert(estimator); + + Exports exports(dwi); + exports.set_noise_out(argument[1]); + auto opt = get_options("rank"); + if (!opt.empty()) + exports.set_rank_input(opt[0][0]); + opt = get_options("max_dist"); + if (!opt.empty()) + exports.set_max_dist(opt[0][0]); + opt = get_options("voxelcount"); + if (!opt.empty()) + exports.set_voxelcount(opt[0][0]); + + int prec = get_option_value("datatype", 0); // default: single precision + if (dwi.datatype().is_complex()) + prec += 2; // support complex input data + switch (prec) { + case 0: + INFO("select real float32 for processing"); + run(dwi, kernel, estimator, exports); + break; + case 1: + INFO("select real float64 for processing"); + run(dwi, kernel, estimator, exports); + break; + case 2: + INFO("select complex float32 for processing"); + run(dwi, kernel, estimator, exports); + break; + case 3: + INFO("select complex float64 for processing"); + run(dwi, kernel, estimator, exports); + break; + } +} diff --git a/cmd/dwidenoise.cpp b/cmd/dwidenoise.cpp index d90b0e1c7c..ad55fc03a9 100644 --- a/cmd/dwidenoise.cpp +++ b/cmd/dwidenoise.cpp @@ -14,7 +14,6 @@ * For more details, see http://www.mrtrix.org/. */ -#include #include #include @@ -32,12 +31,12 @@ #include "denoise/estimator/mrm2022.h" #include "denoise/estimator/result.h" #include "denoise/exports.h" -#include "denoise/functor.h" #include "denoise/kernel/cuboid.h" #include "denoise/kernel/data.h" #include "denoise/kernel/kernel.h" #include "denoise/kernel/sphere_radius.h" #include "denoise/kernel/sphere_ratio.h" +#include "denoise/recon.h" using namespace MR; using namespace App; @@ -72,7 +71,7 @@ void usage() { + "By default, optimal value shrinkage based on minimisation of the Frobenius norm " "will be used to attenuate eigenvectors based on the estimated noise level. " - "Hard truncation of sub-threshold components" + "Hard truncation of sub-threshold components and inclusion of supra-threshold components" "---which was the behaviour of the dwidenoise command in version 3.0.x---" "can be activated using -filter truncate." @@ -85,9 +84,9 @@ void usage() { "There are multiple algebraic forms that modulate the weight with which each decomposition " "contributes with greater or lesser strength toward the output image intensities. " "The various options are: " - "'Gaussian': A Gaussian distribution with FWHM equal to twice the voxel size, " + "'gaussian': A Gaussian distribution with FWHM equal to twice the voxel size, " "such that decompisitions centred more closely to the output voxel have greater influence; " - "'invL0': The inverse of the L0 norm (ie. rank) of each decomposition, " + "'invl0': The inverse of the L0 norm (ie. rank) of each decomposition, " "as used in Manjon et al. 2013; " "'rank': The rank of each decomposition, " "such that high-rank decompositions contribute more strongly to the output intensities " @@ -130,10 +129,10 @@ void usage() { + OptionGroup("Options for modifying PCA computations") + datatype_option + Estimator::option - + Kernel::options + OptionGroup("Options that affect reconstruction of the output image series") + // TODO Separate masks for voxels to contribute to patches vs. voxels for which to perform denoising + Option("mask", "Only denoise voxels within the specified binary brain mask image.") + Argument("image").type_image_in() @@ -249,7 +248,7 @@ void run(Header &data, header.datatype() = DataType::from(); auto output = Image::create(output_name, header); // run - Functor func(data, mask, kernel, estimator, filter, aggregator, exports); + Recon func(data, mask, kernel, estimator, filter, aggregator, exports); ThreadedLoop("running MP-PCA denoising", data, 0, 3).run(func, input, output); // Rescale output if performing aggregation if (aggregator == aggregator_type::EXCLUSIVE) diff --git a/docs/reference/commands/dwi2noise.rst b/docs/reference/commands/dwi2noise.rst new file mode 100644 index 0000000000..cf4a1f1383 --- /dev/null +++ b/docs/reference/commands/dwi2noise.rst @@ -0,0 +1,123 @@ +.. _dwi2noise: + +dwi2noise +=================== + +Synopsis +-------- + +Noise level estimation using Marchenko-Pastur PCA + +Usage +-------- + +:: + + dwi2noise [ options ] dwi noise + +- *dwi*: the input diffusion-weighted image +- *noise*: the output estimated noise level map + +Description +----------- + +DWI data noise map estimation by exploiting data redundancy in the PCA domain using the prior knowledge that the eigenspectrum of random covariance matrices is described by the universal Marchenko-Pastur (MP) distribution. Fitting the MP distribution to the spectrum of patch-wise signal matrices hence provides an estimator of the noise level 'sigma'. + +Unlike the MRtrix3 command dwidenoise, this command does not generate a denoised version of the input image series; its primary output is instead a map of the estimated noise level. While this can also be obtained from the dwidenoise command using option -noise_out, using instead the dwi2noise command gives the ability to obtain a noise map to which filtering can be applied, which can then be utilised for the actual image series denoising, without generating an unwanted intermiedate denoised image series. + +Important note: noise level estimation should only be performed as the first step of an image processing pipeline. The routine is invalid if interpolation or smoothing has been applied to the data prior to denoising. + +Note that on complex input data, the output will be the total noise level across real and imaginary channels, so a scale factor sqrt(2) applies. + +The sliding spatial window behaves differently at the edges of the image FoV depending on the shape / size selected for that window. The default behaviour is to use a spherical kernel centred at the voxel of interest, whose size is some multiple of the number of input volumes; where some such voxels lie outside of the image FoV, the radius of the kernel will be increased until the requisite number of voxels are used. For a spherical kernel of a fixed radius, no such expansion will occur, and so for voxels near the image edge a reduced number of voxels will be present in the kernel. For a cuboid kernel, the centre of the kernel will be offset from the voxel being processed such that the entire volume of the kernel resides within the image FoV. + +The size of the default spherical kernel is set to select a number of voxels that is 1.0 / 0.85 ~ 1.18 times the number of volumes in the input series. If a cuboid kernel is requested, but the -extent option is not specified, the command will select the smallest isotropic patch size that exceeds the number of DW images in the input data; e.g., 5x5x5 for data with <= 125 DWI volumes, 7x7x7 for data with <= 343 DWI volumes, etc. + +Options +------- + +Options for modifying PCA computations +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-datatype float32/float64** Datatype for the eigenvalue decomposition (single or double precision). For complex input data, this will select complex float32 or complex float64 datatypes. + +- **-estimator algorithm** Select the noise level estimator (default = Exp2), either: |br| + * Exp1: the original estimator used in Veraart et al. (2016); |br| + * Exp2: the improved estimator introduced in Cordero-Grande et al. (2019); |br| + * MRM2022: the alternative estimator introduced in Olesen et al. (2022). + +Options for controlling the sliding spatial window kernel +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-shape choice** Set the shape of the sliding spatial window. Options are: cuboid,sphere; default: sphere + +- **-radius_mm value** Set an absolute spherical kernel radius in mm + +- **-radius_ratio value** Set the spherical kernel size as a ratio of number of voxels to number of input volumes (default: 1.0/0.85 ~= 1.18) + +- **-extent window** Set the patch size of the cuboid kernel; can be either a single odd integer or a comma-separated triplet of odd integers + +Options for exporting additional data regarding PCA behaviour +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-rank image** The signal rank estimated for the denoising patch centred at each input image voxel + +Options for debugging the operation of sliding window kernels +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-max_dist image** The maximum distance between a voxel and another voxel that was included in the local denoising patch + +- **-voxelcount image** The number of voxels that contributed to the PCA for processing of each voxel + +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 (caution: using the same file as input and output might cause unexpected behaviour). + +- **-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 +^^^^^^^^^^ + +Veraart, J.; Fieremans, E. & Novikov, D.S. Diffusion MRI noise mapping using random matrix theory. Magn. Res. Med., 2016, 76(5), 1582-1593, doi: 10.1002/mrm.26059 + +Cordero-Grande, L.; Christiaens, D.; Hutter, J.; Price, A.N.; Hajnal, J.V. Complex diffusion-weighted image estimation via matrix recovery under general noise models. NeuroImage, 2019, 200, 391-404, doi: 10.1016/j.neuroimage.2019.06.039 + +* If using -estimator mrm2022: Olesen, J.L.; Ianus, A.; Ostergaard, L.; Shemesh, N.; Jespersen, S.N. Tensor denoising of multidimensional MRI data. Magnetic Resonance in Medicine, 2022, 89(3), 1160-1172 + +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:** Daan Christiaens (daan.christiaens@kcl.ac.uk) and Jelle Veraart (jelle.veraart@nyumc.org) and J-Donald Tournier (jdtournier@gmail.com) and Robert E. Smith (robert.smith@florey.edu.au) + +**Copyright:** Copyright (c) 2016 New York University, University of Antwerp, and the MRtrix3 contributors + +Permission is hereby granted, free of charge, to any non-commercial entity ('Recipient') obtaining a copy of this software and associated documentation files (the 'Software'), to the Software solely for non-commercial research, including the rights to use, copy and modify the Software, subject to the following conditions: + + 1. The above copyright notice and this permission notice shall be included by Recipient in all copies or substantial portions of the Software. + + 2. THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIESOF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BELIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF ORIN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + 3. In no event shall NYU be liable for direct, indirect, special, incidental or consequential damages in connection with the Software. Recipient will defend, indemnify and hold NYU harmless from any claims or liability resulting from the use of the Software by recipient. + + 4. Neither anything contained herein nor the delivery of the Software to recipient shall be deemed to grant the Recipient any right or licenses under any patents or patent application owned by NYU. + + 5. The Software may only be used for non-commercial research and may not be used for clinical care. + + 6. Any publication by Recipient of research involving the Software shall cite the references listed below. + diff --git a/docs/reference/commands/dwidenoise.rst b/docs/reference/commands/dwidenoise.rst index a1a4660b08..c769924859 100644 --- a/docs/reference/commands/dwidenoise.rst +++ b/docs/reference/commands/dwidenoise.rst @@ -21,28 +21,72 @@ Usage Description ----------- -DWI data denoising and noise map estimation by exploiting data redundancy in the PCA domain using the prior knowledge that the eigenspectrum of random covariance matrices is described by the universal Marchenko-Pastur (MP) distribution. Fitting the MP distribution to the spectrum of patch-wise signal matrices hence provides an estimator of the noise level 'sigma', as was first shown in Veraart et al. (2016) and later improved in Cordero-Grande et al. (2019). This noise level estimate then determines the optimal cut-off for PCA denoising. +DWI data denoising and noise map estimation by exploiting data redundancy in the PCA domain using the prior knowledge that the eigenspectrum of random covariance matrices is described by the universal Marchenko-Pastur (MP) distribution. Fitting the MP distribution to the spectrum of patch-wise signal matrices hence provides an estimator of the noise level 'sigma'; this noise level estimate then determines the optimal cut-off for PCA denoising. Important note: image denoising must be performed as the first step of the image processing pipeline. The routine will fail if interpolation or smoothing has been applied to the data prior to denoising. Note that this function does not correct for non-Gaussian noise biases present in magnitude-reconstructed MRI images. If available, including the MRI phase data can reduce such non-Gaussian biases, and the command now supports complex input data. +The sliding spatial window behaves differently at the edges of the image FoV depending on the shape / size selected for that window. The default behaviour is to use a spherical kernel centred at the voxel of interest, whose size is some multiple of the number of input volumes; where some such voxels lie outside of the image FoV, the radius of the kernel will be increased until the requisite number of voxels are used. For a spherical kernel of a fixed radius, no such expansion will occur, and so for voxels near the image edge a reduced number of voxels will be present in the kernel. For a cuboid kernel, the centre of the kernel will be offset from the voxel being processed such that the entire volume of the kernel resides within the image FoV. + +The size of the default spherical kernel is set to select a number of voxels that is 1.0 / 0.85 ~ 1.18 times the number of volumes in the input series. If a cuboid kernel is requested, but the -extent option is not specified, the command will select the smallest isotropic patch size that exceeds the number of DW images in the input data; e.g., 5x5x5 for data with <= 125 DWI volumes, 7x7x7 for data with <= 343 DWI volumes, etc. + +By default, optimal value shrinkage based on minimisation of the Frobenius norm will be used to attenuate eigenvectors based on the estimated noise level. Hard truncation of sub-threshold components and inclusion of supra-threshold components---which was the behaviour of the dwidenoise command in version 3.0.x---can be activated using -filter truncate. + +-aggregation exclusive corresponds to the behaviour of the dwidenoise command in version 3.0.x, where the output intensities for a given image voxel are determined exclusively from the PCA decomposition where the sliding spatial window is centred at that voxel. In all other use cases, so-called "overcomplete local PCA" is performed, where the intensities for an output image voxel are some combination of all PCA decompositions for which that voxel is included in the local spatial kernel. There are multiple algebraic forms that modulate the weight with which each decomposition contributes with greater or lesser strength toward the output image intensities. The various options are: 'gaussian': A Gaussian distribution with FWHM equal to twice the voxel size, such that decompisitions centred more closely to the output voxel have greater influence; 'invl0': The inverse of the L0 norm (ie. rank) of each decomposition, as used in Manjon et al. 2013; 'rank': The rank of each decomposition, such that high-rank decompositions contribute more strongly to the output intensities regardless of distance between the output voxel and the centre of the decomposition kernel; 'uniform': All decompositions that include the output voxel in the sliding spatial window contribute equally. + Options ------- -- **-mask image** Only process voxels within the specified binary brain mask image. +Options for modifying PCA computations +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -- **-extent window** Set the patch size of the denoising filter. By default, the command will select the smallest isotropic patch size that exceeds the number of DW images in the input data, e.g., 5x5x5 for data with <= 125 DWI volumes, 7x7x7 for data with <= 343 DWI volumes, etc. +- **-datatype float32/float64** Datatype for the eigenvalue decomposition (single or double precision). For complex input data, this will select complex float32 or complex float64 datatypes. -- **-noise level** The output noise map, i.e., the estimated noise level 'sigma' in the data.Note that on complex input data, this will be the total noise level across real and imaginary channels, so a scale factor sqrt(2) applies. +- **-estimator algorithm** Select the noise level estimator (default = Exp2), either: |br| + * Exp1: the original estimator used in Veraart et al. (2016); |br| + * Exp2: the improved estimator introduced in Cordero-Grande et al. (2019); |br| + * MRM2022: the alternative estimator introduced in Olesen et al. (2022). -- **-rank cutoff** The selected signal rank of the output denoised image. +Options for controlling the sliding spatial window kernel +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -- **-datatype float32/float64** Datatype for the eigenvalue decomposition (single or double precision). For complex input data, this will select complex float32 or complex float64 datatypes. +- **-shape choice** Set the shape of the sliding spatial window. Options are: cuboid,sphere; default: sphere + +- **-radius_mm value** Set an absolute spherical kernel radius in mm + +- **-radius_ratio value** Set the spherical kernel size as a ratio of number of voxels to number of input volumes (default: 1.0/0.85 ~= 1.18) + +- **-extent window** Set the patch size of the cuboid kernel; can be either a single odd integer or a comma-separated triplet of odd integers + +Options that affect reconstruction of the output image series +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-mask image** Only denoise voxels within the specified binary brain mask image. + +- **-filter choice** Modulate how component contributions are filtered based on the cumulative eigenvalues relative to the noise level; options are: truncate,frobenius; default: frobenius (Optimal Shrinkage based on minimisation of the Frobenius norm) -- **-estimator Exp1/Exp2** Select the noise level estimator (default = Exp2), either: |br| - * Exp1: the original estimator used in Veraart et al. (2016), or |br| - * Exp2: the improved estimator introduced in Cordero-Grande et al. (2019). +- **-aggregator choice** Select how the outcomes of multiple PCA outcomes centred at different voxels contribute to the reconstructed DWI signal in each voxel; options are: exclusive,gaussian,invl0,rank,uniform; default: Gaussian + +Options for exporting additional data regarding PCA behaviour +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-noise_out image** The output noise map, i.e., the estimated noise level 'sigma' in the data. Note that on complex input data, this will be the total noise level across real and imaginary channels, so a scale factor sqrt(2) applies. + +- **-rank_input image** The signal rank estimated for the denoising patch centred at each input image voxel + +- **-rank_output image** An estimated rank for the output image data, accounting for multi-patch aggregation + +Options for debugging the operation of sliding window kernels +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-max_dist image** The maximum distance between a voxel and another voxel that was included in the local denoising patch + +- **-voxelcount image** The number of voxels that contributed to the PCA for processing of each voxel + +- **-sum_aggregation image** The sum of aggregation weights of those patches contributing to each output voxel + +- **-sum_optshrink image** the sum of eigenvector weights computed for the denoising patch centred at each voxel as a result of performing optimal shrinkage Standard options ^^^^^^^^^^^^^^^^ @@ -72,13 +116,17 @@ Veraart, J.; Fieremans, E. & Novikov, D.S. Diffusion MRI noise mapping using ran Cordero-Grande, L.; Christiaens, D.; Hutter, J.; Price, A.N.; Hajnal, J.V. Complex diffusion-weighted image estimation via matrix recovery under general noise models. NeuroImage, 2019, 200, 391-404, doi: 10.1016/j.neuroimage.2019.06.039 +* If using -estimator mrm2022: Olesen, J.L.; Ianus, A.; Ostergaard, L.; Shemesh, N.; Jespersen, S.N. Tensor denoising of multidimensional MRI data. Magnetic Resonance in Medicine, 2022, 89(3), 1160-1172 + +* If using anything other than -aggregation exclusive: Manjon, J.V.; Coupe, P.; Concha, L.; Buades, A.; D. Collins, D.L.; Robles, M. Diffusion Weighted Image Denoising Using Overcomplete Local PCA. PLoS ONE, 2013, 8(9), e73021 + 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:** Daan Christiaens (daan.christiaens@kcl.ac.uk) and Jelle Veraart (jelle.veraart@nyumc.org) and J-Donald Tournier (jdtournier@gmail.com) +**Author:** Daan Christiaens (daan.christiaens@kcl.ac.uk) and Jelle Veraart (jelle.veraart@nyumc.org) and J-Donald Tournier (jdtournier@gmail.com) and Robert E. Smith (robert.smith@florey.edu.au) **Copyright:** Copyright (c) 2016 New York University, University of Antwerp, and the MRtrix3 contributors diff --git a/docs/reference/commands_list.rst b/docs/reference/commands_list.rst index 78cfaa506e..b32465b80e 100644 --- a/docs/reference/commands_list.rst +++ b/docs/reference/commands_list.rst @@ -34,6 +34,7 @@ List of MRtrix3 commands commands/dwi2adc commands/dwi2fod commands/dwi2mask + commands/dwi2noise commands/dwi2response commands/dwi2tensor commands/dwibiascorrect @@ -167,6 +168,7 @@ List of MRtrix3 commands |cpp.png|, :ref:`dwi2adc`, "Calculate ADC and/or IVIM parameters." |cpp.png|, :ref:`dwi2fod`, "Estimate fibre orientation distributions from diffusion data using spherical deconvolution" |python.png|, :ref:`dwi2mask`, "Generate a binary mask from DWI data" + |cpp.png|, :ref:`dwi2noise`, "Noise level estimation using Marchenko-Pastur PCA" |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" diff --git a/src/denoise/estimate.cpp b/src/denoise/estimate.cpp new file mode 100644 index 0000000000..c2c19a269a --- /dev/null +++ b/src/denoise/estimate.cpp @@ -0,0 +1,120 @@ +/* Copyright (c) 2008-2024 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/. + */ + +#include "denoise/estimate.h" + +#include + +#include "math/math.h" + +namespace MR::Denoise { + +template +Estimate::Estimate(const Header &header, + Image &mask, + std::shared_ptr kernel, + std::shared_ptr estimator, + Exports &exports) + : m(header.size(3)), + mask(mask), + kernel(kernel), + estimator(estimator), + X(m, kernel->estimated_size()), + XtX(std::min(m, kernel->estimated_size()), std::min(m, kernel->estimated_size())), + eig(std::min(m, kernel->estimated_size())), + s(std::min(m, kernel->estimated_size())), + exports(exports) {} + +template void Estimate::operator()(Image &dwi) { + + // Process voxels in mask only + if (mask.valid()) { + assign_pos_of(dwi, 0, 3).to(mask); + if (!mask.value()) + return; + } + + // Load list of voxels from which to load data + neighbourhood = (*kernel)({dwi.index(0), dwi.index(1), dwi.index(2)}); + const ssize_t n = neighbourhood.voxels.size(); + const ssize_t r = std::min(m, n); + const ssize_t q = std::max(m, n); + + // Expand local storage if necessary + if (n > X.cols()) { + DEBUG("Expanding data matrix storage from " + str(m) + "x" + str(X.cols()) + " to " + str(m) + "x" + str(n)); + X.resize(m, n); + } + if (r > XtX.cols()) { + DEBUG("Expanding decomposition matrix storage from " + str(X.rows()) + " to " + str(r)); + XtX.resize(r, r); + s.resize(r); + } + + // Fill matrices with NaN when in debug mode; + // make sure results from one voxel are not creeping into another + // due to use of block oberations to prevent memory re-allocation + // in the presence of variation in kernel sizes +#ifndef NDEBUG + X.fill(std::numeric_limits::signaling_NaN()); + XtX.fill(std::numeric_limits::signaling_NaN()); + s.fill(std::numeric_limits::signaling_NaN()); +#endif + + load_data(dwi, neighbourhood.voxels); + + // Compute Eigendecomposition: + if (m <= n) + XtX.topLeftCorner(r, r).template triangularView() = X.leftCols(n) * X.leftCols(n).adjoint(); + else + XtX.topLeftCorner(r, r).template triangularView() = X.leftCols(n).adjoint() * X.leftCols(n); + eig.compute(XtX.topLeftCorner(r, r)); + // eigenvalues sorted in increasing order: + s.head(r) = eig.eigenvalues().template cast(); + + // Marchenko-Pastur optimal threshold determination + threshold = (*estimator)(s, m, n); + const ssize_t in_rank = r - threshold.cutoff_p; + + // Store additional output maps if requested + if (exports.noise_out.valid()) { + assign_pos_of(dwi, 0, 3).to(exports.noise_out); + exports.noise_out.value() = float(std::sqrt(threshold.sigma2)); + } + if (exports.rank_input.valid()) { + assign_pos_of(dwi, 0, 3).to(exports.rank_input); + exports.rank_input.value() = in_rank; + } + if (exports.max_dist.valid()) { + assign_pos_of(dwi, 0, 3).to(exports.max_dist); + exports.max_dist.value() = neighbourhood.max_distance; + } + if (exports.voxelcount.valid()) { + assign_pos_of(dwi, 0, 3).to(exports.voxelcount); + exports.voxelcount.value() = n; + } +} + +template void Estimate::load_data(Image &image, const std::vector &voxels) { + const Kernel::Voxel::index_type pos({image.index(0), image.index(1), image.index(2)}); + for (ssize_t i = 0; i != voxels.size(); ++i) { + assign_pos_of(voxels[i].index, 0, 3).to(image); + X.col(i) = image.row(3); + } + assign_pos_of(pos, 0, 3).to(image); +} + +} // namespace MR::Denoise diff --git a/src/denoise/functor.h b/src/denoise/estimate.h similarity index 57% rename from src/denoise/functor.h rename to src/denoise/estimate.h index 11cf94c42a..babbabfb60 100644 --- a/src/denoise/functor.h +++ b/src/denoise/estimate.h @@ -16,9 +16,7 @@ #pragma once -#include #include -#include #include @@ -34,58 +32,45 @@ namespace MR::Denoise { -template class Functor { +template class Estimate { public: using MatrixType = Eigen::Matrix; - Functor(const Header &header, - Image &mask, - std::shared_ptr kernel, - std::shared_ptr estimator, - filter_type filter, - aggregator_type aggregator, - Exports &exports); + Estimate(const Header &header, + Image &mask, + std::shared_ptr kernel, + std::shared_ptr estimator, + Exports &exports); - void operator()(Image &dwi, Image &out); + void operator()(Image &dwi); -private: - // Denoising configuration +protected: const ssize_t m; + + // Denoising configuration Image mask; std::shared_ptr kernel; std::shared_ptr estimator; - filter_type filter; - aggregator_type aggregator; - double gaussian_multiplier; // Reusable memory + Kernel::Data neighbourhood; MatrixType X; MatrixType XtX; Eigen::SelfAdjointEigenSolver eig; eigenvalues_type s; + Estimator::Result threshold; vector_type clam; - vector_type w; // Export images Exports exports; - // Data that can only be written in a thread-safe manner - // Note that this applies not just to this scratch buffer, but also the output image - // (while it would be thread-safe to create a full copy of the output image for each thread - // and combine them only at destruction time, - // this runs the risk of becoming prohibitively large) - // Not placing this within a MutexProtexted<> as the image type is still templated - static std::mutex mutex_aggregator; - void load_data(Image &image, const std::vector &voxels); }; -template std::mutex Functor::mutex_aggregator; - -template class Functor; -template class Functor; -template class Functor; -template class Functor; +template class Estimate; +template class Estimate; +template class Estimate; +template class Estimate; } // namespace MR::Denoise diff --git a/src/denoise/functor.cpp b/src/denoise/functor.cpp deleted file mode 100644 index 4fab62d51f..0000000000 --- a/src/denoise/functor.cpp +++ /dev/null @@ -1,230 +0,0 @@ -/* Copyright (c) 2008-2024 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/. - */ - -#include "denoise/functor.h" - -#include "math/math.h" - -namespace MR::Denoise { - -template -Functor::Functor(const Header &header, - Image &mask, - std::shared_ptr kernel, - std::shared_ptr estimator, - filter_type filter, - aggregator_type aggregator, - Exports &exports) - : m(header.size(3)), - mask(mask), - kernel(kernel), - estimator(estimator), - filter(filter), - aggregator(aggregator), - // FWHM = 2 x cube root of voxel spacings - gaussian_multiplier(-std::log(2.0) / - Math::pow2(std::cbrt(header.spacing(0) * header.spacing(1) * header.spacing(2)))), - X(m, kernel->estimated_size()), - XtX(std::min(m, kernel->estimated_size()), std::min(m, kernel->estimated_size())), - eig(std::min(m, kernel->estimated_size())), - s(std::min(m, kernel->estimated_size())), - clam(std::min(m, kernel->estimated_size())), - w(std::min(m, kernel->estimated_size())), - exports(exports) {} - -template void Functor::operator()(Image &dwi, Image &out) { - // Process voxels in mask only - if (mask.valid()) { - assign_pos_of(dwi, 0, 3).to(mask); - if (!mask.value()) - return; - } - - // Load list of voxels from which to load data - const Kernel::Data neighbourhood = (*kernel)({dwi.index(0), dwi.index(1), dwi.index(2)}); - const ssize_t n = neighbourhood.voxels.size(); - const ssize_t r = std::min(m, n); - const ssize_t q = std::max(m, n); - - // Expand local storage if necessary - if (n > X.cols()) { - DEBUG("Expanding data matrix storage from " + str(m) + "x" + str(X.cols()) + " to " + str(m) + "x" + str(n)); - X.resize(m, n); - } - if (r > XtX.cols()) { - DEBUG("Expanding decomposition matrix storage from " + str(X.rows()) + " to " + str(r)); - XtX.resize(r, r); - s.resize(r); - clam.resize(r); - w.resize(r); - } - - // Fill matrices with NaN when in debug mode; - // make sure results from one voxel are not creeping into another - // due to use of block oberations to prevent memory re-allocation - // in the presence of variation in kernel sizes -#ifndef NDEBUG - X.fill(std::numeric_limits::signaling_NaN()); - XtX.fill(std::numeric_limits::signaling_NaN()); - s.fill(std::numeric_limits::signaling_NaN()); - clam.fill(std::numeric_limits::signaling_NaN()); - w.fill(std::numeric_limits::signaling_NaN()); -#endif - - load_data(dwi, neighbourhood.voxels); - - // Compute Eigendecomposition: - if (m <= n) - XtX.topLeftCorner(r, r).template triangularView() = X.leftCols(n) * X.leftCols(n).adjoint(); - else - XtX.topLeftCorner(r, r).template triangularView() = X.leftCols(n).adjoint() * X.leftCols(n); - eig.compute(XtX.topLeftCorner(r, r)); - // eigenvalues sorted in increasing order: - s.head(r) = eig.eigenvalues().template cast(); - - // Marchenko-Pastur optimal threshold determination - const Estimator::Result threshold = (*estimator)(s, m, n); - - // Generate weights vector - double sum_weights = 0.0; - ssize_t out_rank = 0; - switch (filter) { - case filter_type::TRUNCATE: - out_rank = r - threshold.cutoff_p; - w.head(threshold.cutoff_p).setZero(); - w.segment(threshold.cutoff_p, r - threshold.cutoff_p).setOnes(); - sum_weights = double(out_rank); - break; - case filter_type::FROBENIUS: { - const double beta = r / q; - const double transition = 1.0 + std::sqrt(beta); - double clam = 0.0; - for (ssize_t i = 0; i != r; ++i) { - const double lam = std::max(s[i], 0.0) / q; - clam += lam; - const double y = clam / (threshold.sigma2 * (i + 1)); - double nu = 0.0; - if (y > transition) { - nu = std::sqrt(Math::pow2(Math::pow2(y) - beta - 1.0) - (4.0 * beta)) / y; - ++out_rank; - } - w[i] = nu / y; - sum_weights += w[i]; - } - } break; - default: - assert(false); - } - - // recombine data using only eigenvectors above threshold - // If only the data computed when this voxel was the centre of the patch - // is to be used for synthesis of the output image, - // then only that individual column needs to be reconstructed; - // if however the result from this patch is to contribute to the synthesized image - // for all voxels that were utilised within this patch, - // then we need to instead compute the full projection - switch (aggregator) { - case aggregator_type::EXCLUSIVE: - if (m <= n) - X.col(neighbourhood.centre_index) = - eig.eigenvectors() * (w.head(r).cast().matrix().asDiagonal() * - (eig.eigenvectors().adjoint() * X.col(neighbourhood.centre_index))); - else - X.col(neighbourhood.centre_index) = - X.leftCols(n) * (eig.eigenvectors() * (w.head(r).cast().matrix().asDiagonal() * - eig.eigenvectors().adjoint().col(neighbourhood.centre_index))); - assign_pos_of(dwi).to(out); - out.row(3) = X.col(neighbourhood.centre_index); - if (exports.sum_aggregation.valid()) { - assign_pos_of(dwi, 0, 3).to(exports.sum_aggregation); - exports.sum_aggregation.value() = 1.0; - } - if (exports.rank_output.valid()) { - assign_pos_of(dwi, 0, 3).to(exports.rank_output); - exports.rank_output.value() = out_rank; - } - break; - default: { - if (m <= n) - X = eig.eigenvectors() * (w.head(r).cast().matrix().asDiagonal() * (eig.eigenvectors().adjoint() * X)); - else - X.leftCols(n) = X.leftCols(n) * - (eig.eigenvectors() * (w.head(r).cast().matrix().asDiagonal() * eig.eigenvectors().adjoint())); - std::lock_guard lock(mutex_aggregator); - for (size_t voxel_index = 0; voxel_index != neighbourhood.voxels.size(); ++voxel_index) { - assign_pos_of(neighbourhood.voxels[voxel_index].index, 0, 3).to(out); - assign_pos_of(neighbourhood.voxels[voxel_index].index).to(exports.sum_aggregation); - double weight = std::numeric_limits::signaling_NaN(); - switch (aggregator) { - case aggregator_type::EXCLUSIVE: - assert(false); - break; - case aggregator_type::GAUSSIAN: - weight = std::exp(gaussian_multiplier * neighbourhood.voxels[voxel_index].sq_distance); - break; - case aggregator_type::INVL0: - weight = 1.0 / (1 + out_rank); - break; - case aggregator_type::RANK: - weight = out_rank; - break; - case aggregator_type::UNIFORM: - weight = 1.0; - break; - } - out.row(3) += weight * X.col(voxel_index); - exports.sum_aggregation.value() += weight; - if (exports.rank_output.valid()) { - assign_pos_of(neighbourhood.voxels[voxel_index].index, 0, 3).to(exports.rank_output); - exports.rank_output.value() += weight * out_rank; - } - } - } break; - } - - // Store additional output maps if requested - if (exports.noise_out.valid()) { - assign_pos_of(dwi, 0, 3).to(exports.noise_out); - exports.noise_out.value() = float(std::sqrt(threshold.sigma2)); - } - if (exports.rank_input.valid()) { - assign_pos_of(dwi, 0, 3).to(exports.rank_input); - exports.rank_input.value() = out_rank; - } - if (exports.sum_optshrink.valid()) { - assign_pos_of(dwi, 0, 3).to(exports.sum_optshrink); - exports.sum_optshrink.value() = sum_weights; - } - if (exports.max_dist.valid()) { - assign_pos_of(dwi, 0, 3).to(exports.max_dist); - exports.max_dist.value() = neighbourhood.max_distance; - } - if (exports.voxelcount.valid()) { - assign_pos_of(dwi, 0, 3).to(exports.voxelcount); - exports.voxelcount.value() = n; - } -} - -template void Functor::load_data(Image &image, const std::vector &voxels) { - const Kernel::Voxel::index_type pos({image.index(0), image.index(1), image.index(2)}); - for (ssize_t i = 0; i != voxels.size(); ++i) { - assign_pos_of(voxels[i].index, 0, 3).to(image); - X.col(i) = image.row(3); - } - assign_pos_of(pos, 0, 3).to(image); -} - -} // namespace MR::Denoise diff --git a/src/denoise/recon.cpp b/src/denoise/recon.cpp new file mode 100644 index 0000000000..4304e31aa8 --- /dev/null +++ b/src/denoise/recon.cpp @@ -0,0 +1,163 @@ +/* Copyright (c) 2008-2024 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/. + */ + +#include "denoise/recon.h" + +#include "math/math.h" + +namespace MR::Denoise { + +template +Recon::Recon(const Header &header, + Image &mask, + std::shared_ptr kernel, + std::shared_ptr estimator, + filter_type filter, + aggregator_type aggregator, + Exports &exports) + : Estimate(header, mask, kernel, estimator, exports), + filter(filter), + aggregator(aggregator), + // FWHM = 2 x cube root of voxel spacings + gaussian_multiplier(-std::log(2.0) / + Math::pow2(std::cbrt(header.spacing(0) * header.spacing(1) * header.spacing(2)))), + w(std::min(Estimate::m, kernel->estimated_size())) {} + +template void Recon::operator()(Image &dwi, Image &out) { + + Estimate (*this)(dwi); + + const ssize_t n = Estimate::neighbourhood.voxels.size(); + const ssize_t r = std::min(Estimate::m, n); + const ssize_t q = std::max(Estimate::m, n); + const ssize_t in_rank = r - Estimate::threshold.cutoff_p; + + if (r > w.size()) + w.resize(r); +#ifndef NDEBUG + w.fill(std::numeric_limits::signaling_NaN()); +#endif + + // Generate weights vector + double sum_weights = 0.0; + ssize_t out_rank = 0; + switch (filter) { + case filter_type::TRUNCATE: + out_rank = in_rank; + w.head(Estimate::threshold.cutoff_p).setZero(); + w.segment(Estimate::threshold.cutoff_p, in_rank).setOnes(); + sum_weights = double(out_rank); + break; + case filter_type::FROBENIUS: { + const double beta = r / q; + const double transition = 1.0 + std::sqrt(beta); + double clam = 0.0; + for (ssize_t i = 0; i != r; ++i) { + const double lam = std::max(Estimate::s[i], 0.0) / q; + clam += lam; + const double y = clam / (Estimate::threshold.sigma2 * (i + 1)); + double nu = 0.0; + if (y > transition) { + nu = std::sqrt(Math::pow2(Math::pow2(y) - beta - 1.0) - (4.0 * beta)) / y; + ++out_rank; + } + w[i] = nu / y; + sum_weights += w[i]; + } + } break; + default: + assert(false); + } + + // recombine data using only eigenvectors above threshold + // If only the data computed when this voxel was the centre of the patch + // is to be used for synthesis of the output image, + // then only that individual column needs to be reconstructed; + // if however the result from this patch is to contribute to the synthesized image + // for all voxels that were utilised within this patch, + // then we need to instead compute the full projection + // TODO Use a new data member local to Recon<> + switch (aggregator) { + case aggregator_type::EXCLUSIVE: + if (Estimate::m <= n) + Estimate::X.col(Estimate::neighbourhood.centre_index) = + Estimate::eig.eigenvectors() * + (w.head(r).cast().matrix().asDiagonal() * + (Estimate::eig.eigenvectors().adjoint() * Estimate::X.col(Estimate::neighbourhood.centre_index))); + else + Estimate::X.col(Estimate::neighbourhood.centre_index) = + Estimate::X.leftCols(n) * + (Estimate::eig.eigenvectors() * + (w.head(r).cast().matrix().asDiagonal() * + Estimate::eig.eigenvectors().adjoint().col(Estimate::neighbourhood.centre_index))); + assign_pos_of(dwi).to(out); + out.row(3) = Estimate::X.col(Estimate::neighbourhood.centre_index); + if (Estimate::exports.sum_aggregation.valid()) { + assign_pos_of(dwi, 0, 3).to(Estimate::exports.sum_aggregation); + Estimate::exports.sum_aggregation.value() = 1.0; + } + if (Estimate::exports.rank_output.valid()) { + assign_pos_of(dwi, 0, 3).to(Estimate::exports.rank_output); + Estimate::exports.rank_output.value() = out_rank; + } + break; + default: { + if (Estimate::m <= n) + Estimate::X = Estimate::eig.eigenvectors() * (w.head(r).cast().matrix().asDiagonal() * + (Estimate::eig.eigenvectors().adjoint() * Estimate::X)); + else + Estimate::X.leftCols(n) = + Estimate::X.leftCols(n) * (Estimate::eig.eigenvectors() * (w.head(r).cast().matrix().asDiagonal() * + Estimate::eig.eigenvectors().adjoint())); + std::lock_guard lock(mutex_aggregator); + for (size_t voxel_index = 0; voxel_index != Estimate::neighbourhood.voxels.size(); ++voxel_index) { + assign_pos_of(Estimate::neighbourhood.voxels[voxel_index].index, 0, 3).to(out); + assign_pos_of(Estimate::neighbourhood.voxels[voxel_index].index).to(Estimate::exports.sum_aggregation); + double weight = std::numeric_limits::signaling_NaN(); + switch (aggregator) { + case aggregator_type::EXCLUSIVE: + assert(false); + break; + case aggregator_type::GAUSSIAN: + weight = std::exp(gaussian_multiplier * Estimate::neighbourhood.voxels[voxel_index].sq_distance); + break; + case aggregator_type::INVL0: + weight = 1.0 / (1 + out_rank); + break; + case aggregator_type::RANK: + weight = out_rank; + break; + case aggregator_type::UNIFORM: + weight = 1.0; + break; + } + out.row(3) += weight * Estimate::X.col(voxel_index); + Estimate::exports.sum_aggregation.value() += weight; + if (Estimate::exports.rank_output.valid()) { + assign_pos_of(Estimate::neighbourhood.voxels[voxel_index].index, 0, 3).to(Estimate::exports.rank_output); + Estimate::exports.rank_output.value() += weight * out_rank; + } + } + } break; + } + + if (Estimate::exports.sum_optshrink.valid()) { + assign_pos_of(dwi, 0, 3).to(Estimate::exports.sum_optshrink); + Estimate::exports.sum_optshrink.value() = sum_weights; + } +} + +} // namespace MR::Denoise diff --git a/src/denoise/recon.h b/src/denoise/recon.h new file mode 100644 index 0000000000..946d01d711 --- /dev/null +++ b/src/denoise/recon.h @@ -0,0 +1,67 @@ +/* Copyright (c) 2008-2024 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/. + */ + +#pragma once + +#include +#include +#include + +#include + +#include "denoise/estimate.h" +#include "denoise/estimator/base.h" +#include "denoise/exports.h" +#include "denoise/kernel/base.h" +#include "header.h" +#include "image.h" + +namespace MR::Denoise { + +template class Recon : public Estimate { + +public: + Recon(const Header &header, + Image &mask, + std::shared_ptr kernel, + std::shared_ptr estimator, + filter_type filter, + aggregator_type aggregator, + Exports &exports); + + void operator()(Image &dwi, Image &out); + +protected: + // Denoising configuration + filter_type filter; + aggregator_type aggregator; + double gaussian_multiplier; + + // Reusable memory + vector_type w; + + // Some data can only be written in a thread-safe manner + static std::mutex mutex_aggregator; +}; + +template std::mutex Recon::mutex_aggregator; + +template class Recon; +template class Recon; +template class Recon; +template class Recon; + +} // namespace MR::Denoise