diff --git a/cmd/dwi2noise.cpp b/cmd/dwi2noise.cpp new file mode 100644 index 0000000000..a6d6be0158 --- /dev/null +++ b/cmd/dwi2noise.cpp @@ -0,0 +1,247 @@ +/* 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 "axes.h" +#include "command.h" +#include "denoise/denoise.h" +#include "denoise/estimate.h" +#include "denoise/estimator/estimator.h" +#include "denoise/exports.h" +#include "denoise/kernel/kernel.h" +#include "denoise/precondition.h" +#include "denoise/subsample.h" +#include "dwi/gradient.h" +#include "exception.h" +#include "filter/demodulate.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." + + + demodulation_description + + + Kernel::shape_description + + + Kernel::default_size_description + + + Kernel::cuboid_size_description; + + AUTHOR = "Robert E. Smith (robert.smith@florey.edu.au)" + " and Daan Christiaens (daan.christiaens@kcl.ac.uk)" + " and Jelle Veraart (jelle.veraart@nyumc.org)" + " and J-Donald Tournier (jdtournier@gmail.com)"; + + 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" + + + "* If using -estimator med: " + "Gavish, M.; Donoho, D.L. " + "The Optimal Hard Threshold for Singular Values is 4/sqrt(3). " + "IEEE Transactions on Information Theory, 2014, 60(8), 5040-5053."; + + 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::estimator_option + + Kernel::options + + subsample_option + + precondition_options(false) + + + DWI::GradImportOptions() + + DWI::GradExportOptions() + + + OptionGroup("Options for exporting additional data regarding PCA behaviour") + + Option("rank", + "The signal rank estimated for each denoising patch") + + Argument("image").type_image_out() + + OptionGroup("Options for debugging the operation of sliding window kernels") + + Option("max_dist", + "The maximum distance between the centre of the patch and a voxel that was included within that patch") + + Argument("image").type_image_out() + + Option("voxelcount", + "The number of voxels that contributed to the PCA for processing of each patch") + + Argument("image").type_image_out() + + Option("patchcount", + "The number of unique patches to which an input image voxel contributes") + + Argument("image").type_image_out(); + +} +// clang-format on + +template +void run(Image &input, + std::shared_ptr subsample, + std::shared_ptr kernel, + std::shared_ptr estimator, + Exports &exports) { + Estimate func(input, subsample, kernel, estimator, exports); + ThreadedLoop("running MP-PCA noise level estimation", input, 0, 3).run(func, input); +} + +template +void run(Header &dwi, + const Demodulation &demodulation, + const demean_type demean, + Image &vst_noise_image, + std::shared_ptr subsample, + std::shared_ptr kernel, + std::shared_ptr estimator, + Exports &exports) { + auto opt_preconditioned = get_options("preconditioned"); + if (!demodulation && demean == demean_type::NONE && !vst_noise_image.valid()) { + if (!opt_preconditioned.empty()) { + WARN("-preconditioned option ignored: no preconditioning taking place"); + } + Image input = dwi.get_image().with_direct_io(3); + run(input, subsample, kernel, estimator, exports); + return; + } + Image input(dwi.get_image()); + const Precondition preconditioner(input, demodulation, demean, vst_noise_image); + Header H_preconditioned(input); + Stride::set(H_preconditioned, Stride::contiguous_along_axis(3, input)); + Image input_preconditioned; + input_preconditioned = opt_preconditioned.empty() + ? Image::scratch(H_preconditioned, "Preconditioned version of \"" + input.name() + "\"") + : Image::create(opt_preconditioned[0][0], H_preconditioned); + preconditioner(input, input_preconditioned, false); + run(input_preconditioned, subsample, kernel, estimator, exports); + if (vst_noise_image.valid()) { + Interp::Cubic> vst(vst_noise_image); + const Transform transform(exports.noise_out); + for (auto l = Loop(exports.noise_out)(exports.noise_out); l; ++l) { + vst.scanner(transform.voxel2scanner * Eigen::Vector3d({default_type(exports.noise_out.index(0)), + default_type(exports.noise_out.index(1)), + default_type(exports.noise_out.index(2))})); + exports.noise_out.value() *= vst.value(); + } + } + if (preconditioner.rank() == 1 && exports.rank_input.valid()) { + for (auto l = Loop(exports.rank_input)(exports.rank_input); l; ++l) + exports.rank_input.value() = + std::max(uint16_t(exports.rank_input.value()) + uint16_t(1), uint16_t(dwi.size(3))); + } +} + +void run() { + auto dwi = Header::open(argument[0]); + if (dwi.ndim() != 4 || dwi.size(3) <= 1) + throw Exception("input image must be 4-dimensional"); + bool complex = dwi.datatype().is_complex(); + + const Demodulation demodulation = select_demodulation(dwi); + const demean_type demean = select_demean(dwi); + Image vst_noise_image; + auto opt = get_options("vst"); + if (!opt.empty()) + vst_noise_image = Image::open(opt[0][0]); + + auto subsample = Subsample::make(dwi, Denoise::default_subsample_ratio); + assert(subsample); + + auto kernel = Kernel::make_kernel(dwi, subsample->get_factors()); + assert(kernel); + + auto estimator = Estimator::make_estimator(vst_noise_image, false); + assert(estimator); + + Exports exports(dwi, subsample->header()); + exports.set_noise_out(argument[1]); + 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]); + opt = get_options("patchcount"); + if (!opt.empty()) + exports.set_patchcount(opt[0][0]); + + int prec = get_option_value("datatype", 0); // default: single precision + if (complex) + prec += 2; // support complex input data + switch (prec) { + case 0: + assert(demodulation.axes.empty()); + INFO("select real float32 for processing"); + run(dwi, demodulation, demean, vst_noise_image, subsample, kernel, estimator, exports); + break; + case 1: + assert(demodulation.axes.empty()); + INFO("select real float64 for processing"); + run(dwi, demodulation, demean, vst_noise_image, subsample, kernel, estimator, exports); + break; + case 2: + INFO("select complex float32 for processing"); + run(dwi, demodulation, demean, vst_noise_image, subsample, kernel, estimator, exports); + break; + case 3: + INFO("select complex float64 for processing"); + run(dwi, demodulation, demean, vst_noise_image, subsample, kernel, estimator, exports); + break; + } +} diff --git a/cmd/dwidenoise.cpp b/cmd/dwidenoise.cpp index 20d9242cd4..8e64b3f618 100644 --- a/cmd/dwidenoise.cpp +++ b/cmd/dwidenoise.cpp @@ -15,6 +15,7 @@ */ #include "command.h" +#include "denoise/denoise.h" #include "image.h" #include @@ -32,7 +33,7 @@ void usage() { SYNOPSIS = "dMRI noise level estimation and denoising using Marchenko-Pastur PCA"; DESCRIPTION - + "DWI data denoising and noise map estimation" + + "DWI data noise map estimation and denoising" " 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." @@ -42,14 +43,9 @@ void usage() { " and later improved in Cordero-Grande et al. (2019)." " 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." + + Denoise::first_step_description - + "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."; + + Denoise::non_gaussian_noise_description; AUTHOR = "Daan Christiaens (daan.christiaens@kcl.ac.uk)" " and Jelle Veraart (jelle.veraart@nyumc.org)" diff --git a/cmd/dwidenoise2.cpp b/cmd/dwidenoise2.cpp new file mode 100644 index 0000000000..72e2789556 --- /dev/null +++ b/cmd/dwidenoise2.cpp @@ -0,0 +1,477 @@ +/* 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 "command.h" +#include "filter/demodulate.h" +#include "header.h" +#include "image.h" +#include "stride.h" + +#include +#include + +#include "denoise/denoise.h" +#include "denoise/estimate.h" +#include "denoise/estimator/base.h" +#include "denoise/estimator/estimator.h" +#include "denoise/estimator/exp.h" +#include "denoise/estimator/mrm2022.h" +#include "denoise/estimator/rank.h" +#include "denoise/estimator/result.h" +#include "denoise/exports.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/precondition.h" +#include "denoise/recon.h" +#include "denoise/subsample.h" + +using namespace MR; +using namespace App; +using namespace MR::Denoise; + +// clang-format off +void usage() { + + SYNOPSIS = "Improved dMRI denoising using Marchenko-Pastur PCA"; + + AUTHOR = "Robert E. Smith (robert.smith@florey.edu.au)" + " and Daan Christiaens (daan.christiaens@kcl.ac.uk)" + " and Jelle Veraart (jelle.veraart@nyumc.org)" + " and J-Donald Tournier (jdtournier@gmail.com)"; + + DESCRIPTION + + "This command performs DWI data denoising," + " additionally with data-driven noise map estimation if not provided explicitly." + " The output denoised DWI data is formed based on filtering of eigenvectors of PCA decompositions:" + " for a set of patches each of which consists of a set of proximal voxels," + " the PCA decomposition is applied," + " and the DWI signal for those voxels is reconstructed" + " where the contribution of each eigenvector is modulated based on its classification as noise." + " In many use cases," + " a threshold that classifies eigenvalues as belonging to signal of interest vs. random thermal noise" + " is based on the prior knowledge that the eigenspectrum of random covariance matrices" + " is described by the universal Marchenko-Pastur (MP) distribution." + + + "This command includes many capabilities absent from the original dwidenoise command. " + "These include:" + " - Multiple sliding window kernel shapes," + " including a spherical kernel that dilates at image edges to preserve aspect ratio;" + " - A greater number of mechanisms for noise level estimation," + " including taking a pre-estimated noise map as input;" + " - Preconditioning, including (per-shell) demeaning," + " phase demodulation (linear or nonlinear)," + " and variance-stabilising transform to compensate for within-patch heteroscedasticity;" + " - Overcomplete local PCA;" + " - Subsampling (performing fewer PCAs than there are input voxels);" + " - Optimal shrinkage of eigenvalues." + + + Denoise::first_step_description + + + Denoise::non_gaussian_noise_description + + + demodulation_description + + + Kernel::shape_description + + + Kernel::default_size_description + + + Kernel::cuboid_size_description + + + Denoise::filter_description + + + Denoise::aggregation_description; + + EXAMPLES + + Example("To approximately replicate the behaviour of the original dwidenoise command", + "dwidenoise2 DWI.mif out.mif -shape cuboid -subsample 1 -demodulate none -demean none -filter truncate -aggregator exclusive", + "While this is neither guaranteed to match exactly the output of the original dwidenoise command" + " nor is it a recommended use case," + " it may nevertheless be informative in demonstrating those advanced features of dwidenoise2 active by default" + " that must be explicitly disabled in order to approximate that behaviour.") + + + Example("Denoising multi-echo fMRI data", + "mrcat original_TE1.mif original_TE2.mif original_TE3.mif original_all.mif -axis 4;" + " dwidenoise2 original_all.mif denoised_all.mif;" + " for_each 1 2 3 : mrconvert denoised_all.mif denoised_TEIN.mif -coord 4 IN -axes 0,1,2,3", + "By concatenating echoes along a new fifth image axis," + " the program will by default calculate the mean intensity per echo " + " and regress this from the data prior to PCA; " + " this should reduce the signal rank and increase PCA precision."); + + REFERENCES + + "Veraart, J.; Novikov, D.S.; Christiaens, D.; Ades-aron, B.; Sijbers, J. & Fieremans, E. " // Internal + "Denoising of diffusion MRI using random matrix theory. " + "NeuroImage, 2016, 142, 394-406, doi: 10.1016/j.neuroimage.2016.08.016" + + + "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" + + + "* 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" + + + "* If using -estimator med or -filter optthresh: " + "Gavish, M.; Donoho, D.L." + "The Optimal Hard Threshold for Singular Values is 4/sqrt(3). " + "IEEE Transactions on Information Theory, 2014, 60(8), 5040-5053."; + + ARGUMENTS + + Argument("dwi", "the input diffusion-weighted image.").type_image_in() + + Argument("out", "the output denoised DWI image.").type_image_out(); + + OPTIONS + + OptionGroup("Options for modifying PCA computations") + + datatype_option + + Estimator::estimator_denoise_options + + Kernel::options + + subsample_option + + precondition_options(true) + + + OptionGroup("Options that affect reconstruction of the output image series") + + Option("filter", + "Modulate how component contributions are filtered " + "based on the cumulative eigenvalues relative to the noise level; " + "options are: " + join(filters, ",") + "; " + "default: optshrink (Optimal Shrinkage based on minimisation of the Frobenius norm)") + + Argument("choice").type_choice(filters) + + Option("aggregator", + "Select how the outcomes of multiple PCA outcomes centred at different voxels " + "contribute to the reconstructed DWI signal in each voxel; " + "options are: " + join(aggregators, ",") + "; default: Gaussian") + + Argument("choice").type_choice(aggregators) + // TODO For specifically the Gaussian aggregator, + // should ideally be possible to select the FWHM of the aggregator + + + OptionGroup("Options for exporting additional data regarding PCA behaviour") + + Option("noise_out", + "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.") + + Argument("image").type_image_out() + + Option("rank_input", + "The signal rank estimated for each denoising patch") + + Argument("image").type_image_out() + + Option("rank_output", + "An estimated rank for the output image data, accounting for multi-patch aggregation") + + 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() + + Option("patchcount", + "The number of unique patches to which an image voxel contributes") + + Argument("image").type_image_out() + + Option("sum_aggregation", + "The sum of aggregation weights of those patches contributing to each output voxel") + + Argument("image").type_image_out() + + Option("sum_optshrink", + "the sum of eigenvector weights computed for the denoising patch centred at each voxel " + "as a result of performing optimal shrinkage") + + Argument("image").type_image_out(); + +} +// clang-format on + +// Necessary to allow normalisation by sum of aggregation weights +// where the image type is cdouble, but aggregation weights are float +// (operations combining complex & real types not allowed to be of different precision) +std::complex operator/(const std::complex &c, const float n) { return c / double(n); } + +template +void run(Image &input, + std::shared_ptr subsample, + std::shared_ptr kernel, + std::shared_ptr estimator, + filter_type filter, + aggregator_type aggregator, + Image &output, + Exports &exports) { + Recon func(input, subsample, kernel, estimator, filter, aggregator, exports); + ThreadedLoop("running MP-PCA denoising", input, 0, 3).run(func, input, output); + // Rescale output if aggregation was performed + if (aggregator == aggregator_type::EXCLUSIVE) + return; + for (auto l_voxel = Loop(exports.sum_aggregation)(output, exports.sum_aggregation); l_voxel; ++l_voxel) { + for (auto l_volume = Loop(3)(output); l_volume; ++l_volume) + output.value() /= float(exports.sum_aggregation.value()); + } + if (exports.rank_output.valid()) { + for (auto l = Loop(exports.sum_aggregation)(exports.rank_output, exports.sum_aggregation); l; ++l) + exports.rank_output.value() /= exports.sum_aggregation.value(); + } +} + +template +void run(Header &dwi, + const Demodulation &demodulation, + const demean_type demean, + Image &vst_noise_image, + std::shared_ptr subsample, + std::shared_ptr kernel, + std::shared_ptr estimator, + filter_type filter, + aggregator_type aggregator, + const std::string &output_name, + Exports &exports) { + auto opt_preconditioned_input = get_options("preconditioned_input"); + auto opt_preconditioned_output = get_options("preconditioned_output"); + if (!demodulation && demean == demean_type::NONE && !vst_noise_image.valid()) { + if (!opt_preconditioned_input.empty()) { + WARN("-preconditioned_input option ignored: no preconditioning taking place"); + } + if (!opt_preconditioned_output.empty()) { + WARN("-preconditioned_output option ignored: no preconditioning taking place"); + } + auto input = dwi.get_image().with_direct_io(3); + Header H(dwi); + H.datatype() = DataType::from(); + auto output = Image::create(output_name, H); + run(input, subsample, kernel, estimator, filter, aggregator, output, exports); + return; + } + auto input = dwi.get_image(); + // perform preconditioning + const Precondition preconditioner(input, demodulation, demean, vst_noise_image); + Image input_preconditioned; + input_preconditioned = + opt_preconditioned_input.empty() + ? Image::scratch(preconditioner.header(), "Preconditioned version of \"" + dwi.name() + "\"") + : Image::create(opt_preconditioned_input[0][0], preconditioner.header()); + preconditioner(input, input_preconditioned, false); + Image output = Image::create(output_name, dwi); + Image output_preconditioned; + if (!opt_preconditioned_output.empty()) + output_preconditioned = Image::create(opt_preconditioned_output[0][0], preconditioner.header()); + // If we can make use of the output image for preconditioned denoised data + // rather than explicitly allocating an intermediate scratch buffer for that purpose, + // then do so + else if (DataType::from() == dwi.datatype() && dwi.ndim() == 4) + output_preconditioned = output; + else + output_preconditioned = Image::scratch( // + preconditioner.header(), // + "Scratch buffer for denoised data before undoing preconditioning"); // + // do the denoising itself + run(input_preconditioned, subsample, kernel, estimator, filter, aggregator, output_preconditioned, exports); + // reverse effects of preconditioning + preconditioner(output_preconditioned, output, true); + // compensate for effects of preconditioning where relevant + if (exports.noise_out.valid() && vst_noise_image.valid()) { + Interp::Cubic> vst(vst_noise_image); + const Transform transform(exports.noise_out); + for (auto l = Loop(exports.noise_out)(exports.noise_out); l; ++l) { + vst.scanner(transform.voxel2scanner * Eigen::Vector3d{default_type(exports.noise_out.index(0)), + default_type(exports.noise_out.index(1)), + default_type(exports.noise_out.index(2))}); + exports.noise_out.value() *= vst.value(); + } + } + if (preconditioner.rank() == 1) { + if (exports.rank_input.valid()) { + for (auto l = Loop(exports.rank_input)(exports.rank_input); l; ++l) { + if (exports.rank_input.value() > 0) + exports.rank_input.value() = + std::min(uint16_t(exports.rank_input.value()) + uint16_t(1), uint16_t(dwi.size(3))); + } + } + if (exports.rank_output.valid()) { + for (auto l = Loop(exports.rank_output)(exports.rank_output); l; ++l) + exports.rank_output.value() = std::min(float(exports.rank_output.value()) + 1.0f, float(dwi.size(3))); + } + } +} + +void run() { + auto dwi = Header::open(argument[0]); + if (dwi.ndim() < 4) + throw Exception("input image must be at least 4-dimensional"); + ssize_t intensities_per_voxel = dwi.size(3); + for (ssize_t axis = 4; axis != dwi.ndim(); ++axis) + intensities_per_voxel *= dwi.size(axis); + if (intensities_per_voxel == 1) + throw Exception("input image must be non-singleton across non-spatial dimensions"); + + const Demodulation demodulation = select_demodulation(dwi); + const demean_type demean = select_demean(dwi); + Image vst_noise_image; + auto opt = get_options("vst"); + if (!opt.empty()) { + vst_noise_image = Image::open(opt[0][0]); + if (vst_noise_image.ndim() != 3) + throw Exception("Variance-stabilising noise image must be 3D"); + } + + aggregator_type aggregator = aggregator_type::GAUSSIAN; + opt = get_options("aggregator"); + if (!opt.empty()) + aggregator = aggregator_type(int(opt[0][0])); + + auto subsample = Subsample::make(dwi, aggregator == aggregator_type::EXCLUSIVE ? 1 : default_subsample_ratio); + assert(subsample); + if (aggregator == aggregator_type::EXCLUSIVE && + *std::max_element(subsample->get_factors().begin(), subsample->get_factors().end()) > 1) { + WARN("Utilising subsampling in conjunction with -aggregator exclusive " + "will result in an output image with holes, " + "as not all input voxels will have their own patch"); + } + + auto kernel = Kernel::make_kernel(dwi, subsample->get_factors()); + assert(kernel); + + auto estimator = Estimator::make_estimator(vst_noise_image, true); + assert(estimator); + + filter_type filter = get_options("fixed_rank").empty() ? filter_type::OPTSHRINK : filter_type::TRUNCATE; + opt = get_options("filter"); + if (!opt.empty()) + filter = filter_type(int(opt[0][0])); + + Exports exports(dwi, subsample->header()); + opt = get_options("noise_out"); + if (!opt.empty()) + exports.set_noise_out(opt[0][0]); + opt = get_options("rank_input"); + if (!opt.empty()) + exports.set_rank_input(opt[0][0]); + opt = get_options("rank_output"); + if (!opt.empty()) { + if (aggregator == aggregator_type::EXCLUSIVE && filter == filter_type::TRUNCATE) { + WARN("When using -aggregator exclusive and -filter truncate, " + "the output of -rank_output will be identical to the output of -rank_input, " + "as there is no aggregation of multiple patches per output voxel " + "and no optimal shrinkage to reduce output rank relative to estimated input rank"); + } + exports.set_rank_output(opt[0][0]); + } + opt = get_options("sum_optshrink"); + if (!opt.empty()) { + if (filter == filter_type::TRUNCATE) { + WARN("Note that with a truncation filter, " + "output image from -sumweights option will be equivalent to rank_input"); + } + exports.set_sum_optshrink(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]); + opt = get_options("patchcount"); + if (!opt.empty()) + exports.set_patchcount(opt[0][0]); + + opt = get_options("sum_aggregation"); + if (!opt.empty()) { + if (aggregator == aggregator_type::EXCLUSIVE) { + WARN("Output from -sum_aggregation will just contain 1 for every voxel processed: " + "no patch aggregation takes place when output series comes exclusively from central patch"); + } + exports.set_sum_aggregation(opt[0][0]); + } else if (aggregator != aggregator_type::EXCLUSIVE) { + exports.set_sum_aggregation(""); + } + + 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: + assert(demodulation.axes.empty()); + INFO("select real float32 for processing"); + run( // + dwi, // + demodulation, // + demean, // + vst_noise_image, // + subsample, // + kernel, // + estimator, // + filter, // + aggregator, // + argument[1], // + exports); // + break; + case 1: + assert(demodulation.axes.empty()); + INFO("select real float64 for processing"); + run( // + dwi, // + demodulation, // + demean, // + vst_noise_image, // + subsample, // + kernel, // + estimator, // + filter, // + aggregator, // + argument[1], // + exports); // + break; + case 2: + INFO("select complex float32 for processing"); + run( // + dwi, // + demodulation, // + demean, // + vst_noise_image, // + subsample, // + kernel, // + estimator, // + filter, // + aggregator, // + argument[1], // + exports); // + break; + case 3: + INFO("select complex float64 for processing"); + run( // + dwi, // + demodulation, // + demean, // + vst_noise_image, // + subsample, // + kernel, // + estimator, // + filter, // + aggregator, // + argument[1], // + exports); // + break; + } +} diff --git a/cmd/mrfilter.cpp b/cmd/mrfilter.cpp index e97e2882a4..be2713d0c0 100644 --- a/cmd/mrfilter.cpp +++ b/cmd/mrfilter.cpp @@ -15,28 +15,39 @@ */ #include +#include #include "command.h" #include "filter/base.h" +#include "filter/demodulate.h" #include "filter/gradient.h" +#include "filter/kspace.h" #include "filter/median.h" #include "filter/normalise.h" #include "filter/smooth.h" #include "filter/zclean.h" #include "image.h" +#include "interp/cubic.h" #include "math/fft.h" using namespace MR; using namespace App; -const std::vector filters = {"fft", "gradient", "median", "smooth", "normalise", "zclean"}; +const std::vector filters = { + "demodulate", "fft", "gradient", "kspace", "median", "smooth", "normalise", "zclean"}; // clang-format off -const OptionGroup FFTOption = OptionGroup ("Options for FFT filter") +const OptionGroup FFTAxesOption = OptionGroup ("Options applicable to demodulate / FFT / k-space filters") + Option ("axes", "the axes along which to apply the Fourier Transform." " By default, the transform is applied along the three spatial axes." " Provide as a comma-separate list of axis indices.") - + Argument ("list").type_sequence_int() + + Argument ("list").type_sequence_int(); + +const OptionGroup DemodulateOption = OptionGroup ("Options applicable to demodulate filter") + + Option ("linear", "only demodulate based on a linear phase ramp, " + "rather than a filtered k-space"); + +const OptionGroup FFTOption = OptionGroup ("Options for FFT filter") + Option ("inverse", "apply the inverse FFT") + Option ("magnitude", "output a magnitude image rather than a complex-valued image") + Option ("rescale", "rescale values so that inverse FFT recovers original values") @@ -57,6 +68,15 @@ const OptionGroup GradientOption = OptionGroup ("Options for gradient filter") + Option ("scanner", "define the gradient with respect to" " the scanner coordinate frame of reference."); +const OptionGroup KSpaceOption = OptionGroup ("Options for k-space filtering") + + Option ("window", "specify the shape of the k-space window filter; " + "options are: " + join(Filter::kspace_window_choices, ",") + " " + "(no default; must be specified for \"kspace\" operation)") + + Argument("name").type_choice(Filter::kspace_window_choices) + + Option ("strength", "modulate the strength of the chosen filter " + "(exact interpretation & defaultmay depend on the exact filter chosen)") + + Argument("value").type_float(0.0, 1.0); + const OptionGroup MedianOption = OptionGroup ("Options for median filter") + Option ("extent", "specify extent of median filtering neighbourhood in voxels." " This can be specified either as a single value to be used for all 3 axes," @@ -116,7 +136,7 @@ void usage() { DESCRIPTION + "The available filters are:" - " fft, gradient, median, smooth, normalise, zclean." + " demodulate, fft, gradient, median, smooth, normalise, zclean." + "Each filter has its own unique set of optional parameters." + "For 4D images, each 3D volume is processed independently."; @@ -126,8 +146,11 @@ void usage() { + Argument ("output", "the output image.").type_image_out (); OPTIONS + + FFTAxesOption + + DemodulateOption + FFTOption + GradientOption + + KSpaceOption + MedianOption + NormaliseOption + SmoothOption @@ -137,44 +160,70 @@ void usage() { } // clang-format on +// TODO Use presence of SliceEncodingDirection to select defaults +std::vector get_axes(const Header &H, const std::vector &default_axes) { + auto opt = get_options("axes"); + std::vector axes = default_axes; + if (!opt.empty()) { + axes = parse_ints(opt[0][0]); + for (const auto axis : axes) + if (axis >= H.ndim()) + throw Exception("axis provided with -axes option is out of range"); + if (std::set(axes.begin(), axes.end()).size() != axes.size()) + throw Exception("axis indices must not contain duplicates"); + } + return axes; +} + void run() { const size_t filter_index = argument[1]; switch (filter_index) { - // FFT + // Phase demodulation case 0: { - // FIXME Had to use cdouble throughout; seems to fail at compile time even trying to - // convert between cfloat and cdouble... - auto input = Image::open(argument[0]); + Header H_in = Header::open(argument[0]); + if (!H_in.datatype().is_complex()) + throw Exception("demodulation filter only applicable for complex image data"); + auto input = H_in.get_image(); - std::vector axes = {0, 1, 2}; - auto opt = get_options("axes"); - if (!opt.empty()) { - axes = parse_ints(opt[0][0]); - for (const auto axis : axes) - if (axis >= input.ndim()) - throw Exception("axis provided with -axes option is out of range"); - } + const std::vector inner_axes = get_axes(H_in, {0, 1}); + + Filter::Demodulate filter(input, inner_axes, !get_options("linear").empty()); + + Header H_out(H_in); + Stride::set_from_command_line(H_out); + auto output = Image::create(argument[2], H_out); + + filter(input, output); + } break; + + // FFT + case 1: { + Header H_in = Header::open(argument[0]); + std::vector axes = get_axes(H_in, {0, 1, 2}); const int direction = get_options("inverse").empty() ? FFTW_FORWARD : FFTW_BACKWARD; const bool centre_FFT = !get_options("centre_zero").empty(); const bool magnitude = !get_options("magnitude").empty(); - Header header = input; - Stride::set_from_command_line(header); - header.datatype() = magnitude ? DataType::Float32 : DataType::CFloat64; - auto output = Image::create(argument[2], header); + Header H_out(H_in); + Stride::set_from_command_line(H_out); + H_out.datatype() = magnitude ? DataType::Float32 : DataType::CFloat64; + auto output = Image::create(argument[2], H_out); double scale = 1.0; + // FIXME Had to use cdouble throughout; seems to fail at compile time even trying to + // convert between cfloat and cdouble... + auto input = H_in.get_image(); + Image in(input), out; for (size_t n = 0; n < axes.size(); ++n) { scale *= in.size(axes[n]); - if (n >= (axes.size() - 1) && !magnitude) { + if (n == (axes.size() - 1) && !magnitude) { out = output; - } else { - if (!out.valid()) - out = Image::scratch(input); + } else if (!out.valid()) { + out = Image::scratch(H_in); } Math::FFT(in, out, axes[n], direction, centre_FFT); @@ -187,15 +236,15 @@ void run() { [](decltype(out) &a, decltype(output) &b) { a.value() = abs(cdouble(b.value())); }, output, out); } if (!get_options("rescale").empty()) { - scale = std::sqrt(scale); - ThreadedLoop(out).run([&scale](decltype(out) &a) { a.value() /= scale; }, output); + scale = 1.0 / std::sqrt(scale); + ThreadedLoop(out).run([&scale](decltype(out) &a) { a.value() *= scale; }, output); } break; } // Gradient - case 1: { + case 2: { auto input = Image::open(argument[0]); Filter::Gradient filter(input, !get_options("magnitude").empty()); @@ -223,8 +272,42 @@ void run() { break; } + // k-space filtering + case 3: { + auto opt_window = get_options("window"); + if (opt_window.empty()) + throw Exception("-window option is compulsory for k-space filtering"); + + Header H_in = Header::open(argument[0]); + const std::vector axes = get_axes(H_in, {0, 1, 2}); + const bool is_complex = H_in.datatype().is_complex(); + auto input = H_in.get_image(); + + Image window; + switch (Filter::kspace_windowfn_t(int(opt_window[0][0]))) { + case Filter::kspace_windowfn_t::TUKEY: + window = Filter::KSpace::window_tukey(H_in, axes, get_option_value("strength", Filter::default_tukey_width)); + break; + default: + assert(false); + } + Filter::KSpace filter(H_in, window); + Header H_out(H_in); + + if (is_complex) { + auto output = Image::create(argument[2], H_out); + filter(input, output); + } else { + H_out.datatype() = DataType::Float32; + H_out.datatype().set_byte_order_native(); + auto output = Image::create(argument[2], H_out); + filter(input, output); + } + break; + } + // Median - case 2: { + case 4: { auto input = Image::open(argument[0]); Filter::Median filter(input); @@ -241,7 +324,7 @@ void run() { } // Smooth - case 3: { + case 5: { auto input = Image::open(argument[0]); Filter::Smooth filter(input); @@ -272,7 +355,7 @@ void run() { } // Normalisation - case 4: { + case 6: { auto input = Image::open(argument[0]); Filter::Normalise filter(input); @@ -289,7 +372,7 @@ void run() { } // Zclean - case 5: { + case 7: { auto input = Image::open(argument[0]); Filter::ZClean filter(input); diff --git a/core/filter/demodulate.h b/core/filter/demodulate.h new file mode 100644 index 0000000000..aaa366d3ef --- /dev/null +++ b/core/filter/demodulate.h @@ -0,0 +1,260 @@ +/* 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 "algo/copy.h" +#include "algo/loop.h" +#include "filter/base.h" +#include "filter/kspace.h" +#include "filter/smooth.h" +#include "image.h" +#include "interp/cubic.h" +#include "math/fft.h" +#include "progressbar.h" + +namespace MR::Filter { +/** \addtogroup Filters +@{ */ + +// From Manzano-Patron et al. 2024 +constexpr default_type default_tukey_FWHM_demodulate = 0.58; +// TODO Ideally do more experimentation to figure out a reasonable default here +// Too high and everything ends up in the real axis; +// too low and disjointed phase cound drive up signal rank +constexpr default_type default_tukey_alpha_demodulate = 2.0 * (1.0 - default_tukey_FWHM_demodulate); + +/*! Estimate a linear phase ramp of a complex image and demodulate by such + */ +class Demodulate : public Base { +public: + template + Demodulate(ImageType &in, const std::vector &inner_axes, const bool linear) + : Base(in), + phase(Image::scratch(in, + std::string("Scratch image storing ") // + + (linear ? "linear" : "nonlinear") // + + " phase for demodulator")) { // + + using value_type = typename ImageType::value_type; + using real_type = typename ImageType::value_type::value_type; + + ImageType input(in); + // if (!in.datatype().is_complex()) + // throw Exception("demodulation filter only applicable for complex image data"); + + std::vector outer_axes; + std::vector::const_iterator it = inner_axes.begin(); + for (size_t axis = 0; axis != input.ndim(); ++axis) { + if (it != inner_axes.end() && *it == axis) + ++it; + else + outer_axes.push_back(axis); + } + + ProgressBar progress(std::string("estimating ") + (linear ? "linear" : "nonlinear") + " phase modulation", + inner_axes.size() + 1 + (linear ? 0 : inner_axes.size())); + + // FFT currently hard-wired to double precision; + // have to manually load into cdouble memory + Image kspace; + { + Image temp; + for (ssize_t n = 0; n != inner_axes.size(); ++n) { + switch (n) { + case 0: + // TODO For now doing a straight copy; + // would be preferable if, if the type of the input image is cdouble, + // to just make temp a view of input + // In retrospect it's probably necessary to template the whole k-space derivation + temp = Image::scratch(in, "Scratch k-space for phase demodulator (1)"); + copy(input, temp); + kspace = Image::scratch(in, "Scratch k-space for phase demodulator (2)"); + break; + // TODO If it's possible to prevent the Algo::copy() call above, + // then the below code block can be re-introduced + // case 1: + // temp = kspace; + // kspace = Image::scratch(in, "Scratch k-space for phase demodulator (2)"); + // break; + default: + std::swap(temp, kspace); + break; + } + // Centred FFT if linear (so we can do peak-finding without wraparound); + // not centred if non-linear (for compatibility with window filter generation function) + Math::FFT(temp, kspace, inner_axes[n], FFTW_FORWARD, linear); + ++progress; + } + } + + // TODO Can likely remove "input" from here + auto gen_linear_phase = + [&](Image &input, Image &kspace, Image &phase, const std::vector &axes) { + std::array axis_mask({false, false, false}); + for (auto axis : axes) { + if (axis > 2) + throw Exception("Linear phase demodulator does not support non-spatial inner axes"); + axis_mask[axis] = true; + } + + std::vector index_of_max(kspace.ndim()); + std::vector::const_iterator it = axes.begin(); + for (ssize_t axis = 0; axis != kspace.ndim(); ++axis) { + if (it != axes.end() && *it == axis) { + index_of_max[axis] = -1; + ++it; + } else { + index_of_max[axis] = kspace.index(axis); + } + } + double max_magnitude = -std::numeric_limits::infinity(); + cdouble value_at_max = cdouble(std::numeric_limits::signaling_NaN(), // + std::numeric_limits::signaling_NaN()); // + for (auto l_inner = Loop(axes)(kspace); l_inner; ++l_inner) { + const cdouble value = kspace.value(); + const double magnitude = std::abs(value); + if (magnitude > max_magnitude) { + max_magnitude = magnitude; + value_at_max = value; + for (auto axis : axes) + index_of_max[axis] = kspace.index(axis); + } + } + + using pos_type = Eigen::Matrix; + using gradient_type = Eigen::Matrix; + pos_type pos({real_type(index_of_max[0]), real_type(index_of_max[1]), real_type(index_of_max[2])}); + gradient_type gradient; + cdouble value; + Interp::SplineInterp, + Math::UniformBSpline, + Math::SplineProcessingType::ValueAndDerivative> + interp(kspace); + assign_pos_of(index_of_max).to(interp); + interp.voxel(pos.template cast()); + interp.value_and_gradient(value, gradient); + const double mag = std::abs(value); + for (ssize_t iter = 0; iter != 10; ++iter) { + pos_type grad_mag({0.0, 0.0, 0.0}); // Gradient of the magnitude of the complex k-space data + for (ssize_t axis = 0; axis != 3; ++axis) { + if (axis_mask[axis]) + grad_mag[axis] = (value.real() * gradient[axis].real() + value.imag() * gradient[axis].imag()) / mag; + } + grad_mag /= max_magnitude; + pos += 1.0 * grad_mag; + interp.value_and_gradient(value, gradient); + } + value_at_max = value; + const real_type phase_at_max = std::arg(value_at_max); + max_magnitude = std::abs(value_at_max); + + // Determine direction and frequency of harmonic + pos_type kspace_origin; + for (ssize_t axis = 0; axis != 3; ++axis) + // Do integer division, then convert to floating-point + kspace_origin[axis] = axis_mask[axis] ? real_type((kspace.size(axis) - 1) / 2) : real_type(0); + const pos_type kspace_offset({axis_mask[0] ? (pos[0] - kspace_origin[0]) : real_type(0), + axis_mask[1] ? (pos[1] - kspace_origin[1]) : real_type(0), + axis_mask[2] ? (pos[2] - kspace_origin[2]) : real_type(0)}); + const pos_type cycle_voxel({axis_mask[0] ? (input.size(0) / kspace_offset[0]) : real_type(0), + axis_mask[1] ? (input.size(1) / kspace_offset[1]) : real_type(0), + axis_mask[2] ? (input.size(2) / kspace_offset[2]) : real_type(0)}); + const pos_type voxel_origin({real_type(0.5 * (input.size(0) - 1.0)), + real_type(0.5 * (input.size(1) - 1.0)), + real_type(0.5 * (input.size(2) - 1.0))}); + for (auto l = Loop(axes)(input, phase); l; ++l) { + pos = {real_type(input.index(0)), real_type(input.index(1)), real_type(input.index(2))}; + const real_type voxel_phase = + phase_at_max + (2.0 * Math::pi * (pos - voxel_origin).dot(cycle_voxel) / cycle_voxel.squaredNorm()); + const value_type modulator(std::cos(voxel_phase), std::sin(voxel_phase)); + phase.value() = cfloat(modulator); + } + }; + + auto gen_nonlinear_phase = + [&](Image &input, Image &kspace, Image &phase, const std::vector &axes) { + Image window = Filter::KSpace::window_tukey(input, axes, default_tukey_alpha_demodulate); + Adapter::Replicate> replicating_window(window, in); + for (auto l = Loop(kspace)(kspace, replicating_window); l; ++l) + kspace.value() *= double(replicating_window.value()); + for (auto axis : axes) { + Math::FFT(kspace, kspace, axis, FFTW_BACKWARD, false); + ++progress; + } + for (auto l = Loop(kspace)(kspace, phase); l; ++l) { + cdouble value = cdouble(kspace.value()); + value /= std::sqrt(norm(value)); + phase.value() = {float(value.real()), float(value.imag())}; + } + }; + + if (linear) { + if (outer_axes.size()) { + // TODO Multi-thread + for (auto l_outer = Loop(outer_axes)(input, kspace, phase); l_outer; ++l_outer) + gen_linear_phase(input, kspace, phase, inner_axes); + } else { + gen_linear_phase(input, kspace, phase, inner_axes); + } + } else { + gen_nonlinear_phase(input, kspace, phase, inner_axes); + } + } + + template + void operator()(InputImageType &in, OutputImageType &out, const bool forward = false) { + // if (!out.datatype().is_complex()) + // throw Exception("Cannot modulate phase for output image that is not complex"); + if (forward) { + for (auto l = Loop("re-applying phase modulation")(in, phase, out); l; ++l) + out.value() = typename OutputImageType::value_type(typename InputImageType::value_type(in.value()) * + typename InputImageType::value_type(cfloat(phase.value()))); + } else { + for (auto l = Loop("performing phase demodulation")(in, phase, out); l; ++l) + out.value() = + typename OutputImageType::value_type(typename InputImageType::value_type(in.value()) * + typename InputImageType::value_type(std::conj(cfloat(phase.value())))); + } + } + + template void operator()(ImageType &image, const bool forward = false) { + // if (!image.datatype().is_complex()) + // throw Exception("Cannot modulate phase for output image that is not complex"); + if (forward) { + for (auto l = Loop("re-applying phase modulation")(image, phase); l; ++l) + image.value() = typename ImageType::value_type(typename ImageType::value_type(image.value()) * + typename ImageType::value_type(cfloat(phase.value()))); + } else { + for (auto l = Loop("performing phase demodulation")(image, phase); l; ++l) + image.value() = + typename ImageType::value_type(typename ImageType::value_type(image.value()) * + typename ImageType::value_type(std::conj(cfloat(phase.value())))); + } + } + + Image operator()() const { return phase; } + +protected: + // TODO Change to Image; can produce complex value at processing time + Image phase; + + // TODO Define a protected class that can be utilised to generate the phase image in a multi-threaded manner +}; +//! @} +} // namespace MR::Filter diff --git a/core/filter/kspace.h b/core/filter/kspace.h new file mode 100644 index 0000000000..dd5cb33d78 --- /dev/null +++ b/core/filter/kspace.h @@ -0,0 +1,200 @@ +/* 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 "adapter/replicate.h" +#include "filter/base.h" +#include "header.h" +#include "image.h" +#include "math/fft.h" +#include "types.h" + +namespace MR::Filter { + +const std::vector kspace_window_choices({"tukey"}); +enum class kspace_windowfn_t { TUKEY }; +constexpr default_type default_tukey_width = 0.5; + +/** \addtogroup Filters +@{ */ + +/*! Apply (or reverse) k-space filtering + */ +// TODO Add ability to undo prior filtering +class KSpace : public Base { +public: + KSpace(const Header &H, Image &window) : Base(H), window(window) { + for (size_t axis = 0; axis != H.ndim(); ++axis) { + if (axis < window.ndim() && window.size(axis) > 1) + inner_axes.push_back(axis); + else + outer_axes.push_back(axis); + } + } + + template + typename std::enable_if>::value, void>::type + operator()(InputImageType &in, OutputImageType &out) { + Image kspace; + Image temp; + { + for (ssize_t n = 0; n != inner_axes.size(); ++n) { + switch (n) { + case 0: + kspace = Image::scratch(in, + "Scratch k-space for \"" + in.name() + "\"" // + + (inner_axes.size() > 1 ? " (1 of 2)" : "")); // + Math::FFT(in, kspace, inner_axes[0], FFTW_FORWARD, false); + break; + case 1: + temp = kspace; + kspace = Image::scratch(in, + "Scratch k-space for \"" + in.name() + "\" " // + + "(2 of 2)"); // + Math::FFT(temp, kspace, inner_axes[n], FFTW_FORWARD, false); + break; + default: + std::swap(temp, kspace); + Math::FFT(temp, kspace, inner_axes[n], FFTW_FORWARD, false); + break; + } + } + } + + if (outer_axes.size()) { + Adapter::Replicate> replicating_window(window, in); + for (auto l = Loop(in)(kspace, replicating_window); l; ++l) + kspace.value() *= double(replicating_window.value()); + } else { + for (auto l = Loop(in)(kspace, window); l; ++l) + kspace.value() *= double(window.value()); + } + + for (ssize_t n = 0; n != inner_axes.size(); ++n) { + if (n == inner_axes.size() - 1) { + // Final FFT: + // use output image if applicable; + // perform amplitude transform if necessary + do_final_fft(out, kspace, temp, inner_axes[n]); + } else { + Math::FFT(kspace, temp, inner_axes[n], FFTW_BACKWARD, false); + std::swap(kspace, temp); + } + } + } + + template + typename std::enable_if::value, void>::type + operator()(InputImageType &in, OutputImageType &out) { + Image temp = Image::scratch(in, "Scratch \"" + in.name() + "\" converted to cdouble for FFT"); + for (auto l = Loop(in)(in, temp); l; ++l) + temp.value() = {double(in.value()), double(0)}; + (*this)(temp, out); + } + + // TODO Static functions to create different windows + // (different windows may require different input parameters) + // TODO Currently extending 1D window functions to ND via outer product; + // this is however not a unique choice + static Image window_tukey(const Header &header, // + const std::vector &inner_axes, // + const default_type cosine_frac) { // + assert(cosine_frac >= 0.0 && cosine_frac <= 1.0); + Image window = Image::scratch(make_window_header(header, inner_axes), + "Scratch Tukey filter window with alpha=" + str(cosine_frac)); + for (auto l = Loop(window)(window); l; ++l) + window.value() = 1.0; + for (auto axis : inner_axes) { + const size_t N = header.size(axis); + Eigen::Array window1d = Eigen::Array::Ones(N); + const default_type transition_lower = 0.5 - 0.5 * cosine_frac; + const default_type transition_upper = 0.5 + 0.5 * cosine_frac; + // Beware of FFT being non-centred + for (size_t n = 0; n != N; ++n) { + const default_type pos = default_type(n) / default_type(N); + if (pos > transition_lower && pos < transition_upper) + window1d[n] = 0.5 + 0.5 * std::cos(2.0 * Math::pi * (pos - transition_lower) / cosine_frac); + } + window1d *= 1.0 / double(N); + // Need to loop over all inner axes other than the current one + std::vector inner_excluding_axis; + for (auto a : inner_axes) { + if (a != axis) + inner_excluding_axis.push_back(a); + } + window.reset(); + for (auto l = Loop(inner_excluding_axis)(window); l; ++l) { + for (size_t n = 0; n != N; ++n) { + window.index(axis) = n; + window.value() *= window1d[n]; + } + } + } + return window; + } + +protected: + Image window; + std::vector inner_axes; + std::vector outer_axes; + + template + typename std::enable_if::value, void>::type + do_final_fft(ImageType &out, Image &kspace, Image &scratch, const size_t axis) { + TRACE; + Math::FFT(kspace, out, axis, FFTW_BACKWARD, false); + } + template + typename std::enable_if::value, void>::type + do_final_fft(ImageType &out, Image &kspace, Image &scratch, const size_t axis) { + TRACE; + assert(scratch.valid()); + Math::FFT(kspace, scratch, axis, FFTW_BACKWARD, false); + for (auto l = Loop(out)(scratch, out); l; ++l) + out.value() = {float(cdouble(scratch.value()).real()), float(cdouble(scratch.value()).imag())}; + } + template + typename std::enable_if::value, void>::type + do_final_fft(ImageType &out, Image &kspace, Image &scratch, const size_t axis) { + TRACE; + assert(scratch.valid()); + Math::FFT(kspace, scratch, axis, FFTW_BACKWARD, false); + for (auto l = Loop(out)(scratch, out); l; ++l) + out.value() = std::abs(cdouble(scratch.value())); + } + + static Header make_window_header(const Header &header, const std::vector &inner_axes) { + Header H(header); + H.datatype() = DataType::Float64; + H.datatype().set_byte_order_native(); + std::vector::const_iterator it = inner_axes.begin(); + for (size_t axis = 0; axis != header.ndim(); ++axis) { + if (it != inner_axes.end() && *it == axis) + ++it; + else + H.size(axis) = 1; + } + squeeze_dim(H); + return H; + } +}; +//! @} + +} // namespace MR::Filter diff --git a/docs/reference/commands/dwi2noise.rst b/docs/reference/commands/dwi2noise.rst new file mode 100644 index 0000000000..6235ff54f7 --- /dev/null +++ b/docs/reference/commands/dwi2noise.rst @@ -0,0 +1,162 @@ +.. _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 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. + +If the input data are of complex type, then a smooth non-linear phase will be demodulated removed from each k-space prior to PCA. In the absence of metadata indicating otherwise, it is inferred that the first two axes correspond to acquired slices, and different slices / volumes will be demodulated individually; this behaviour can be modified using the -demod_axes option. A strictly linear phase term can instead be regressed from each k-space, similarly to performed in Cordero-Grande et al. 2019, by specifying -demodulate linear. + +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. + +Permissible sizes for the cuboid kernel depend on the subsampling factor. If no subsampling is performed, or the subsampling factor is odd, then the extent(s) of the kernel must be odd, such that a unique voxel lies at the very centre of each kernel. If however an even subsampling factor is used, then the extent(s) of the kernel must be even, reflecting the fact that it is a voxel corner that resides at the centre of the kernel.In either case, if the extent is specified manually, the user can either provide a single integer---which will determine the number of voxels in the kernel across all three spatial axes---or a comma-separated list of three integers,individually defining the number of voxels in the kernel for all three spatial axes. + +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| + * Med: estimate based on the median eigenvalue as in Gavish and Donohue (2014); |br| + * MRM2022: the alternative estimator introduced in Olesen et al. (2022). |br| + Operation will be bypassed if -noise_in or -fixed_rank are specified + +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 integer or a comma-separated triplet of integers (see Description) + +- **-subsample factor** reduce the number of PCA kernels relative to the number of image voxels; can provide either an integer subsampling factor, or a comma-separated list of three factors; default: 2 + +Options for preconditioning data prior to PCA +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-demodulate mode** select form of phase demodulation; options are: none,linear,nonlinear (default: nonlinear) + +- **-demod_axes axes** comma-separated list of axis indices along which FFT can be applied for phase demodulation + +- **-demean mode** select method of demeaning prior to PCA; options are: none,shells,all (default: 'shells' if DWI gradient table available, 'all' otherwise) + +- **-vst image** apply a within-patch variance-stabilising transformation based on a pre-estimated noise level map + +- **-preconditioned image** export the preconditioned version of the input image that is the input to PCA + +DW gradient table import options +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-grad file** Provide the diffusion-weighted gradient scheme used in the acquisition in a text file. This should be supplied as a 4xN text file with each line in the format [ X Y Z b ], where [ X Y Z ] describe the direction of the applied gradient, and b gives the b-value in units of s/mm^2. If a diffusion gradient scheme is present in the input image header, the data provided with this option will be instead used. + +- **-fslgrad bvecs bvals** Provide the diffusion-weighted gradient scheme used in the acquisition in FSL bvecs/bvals format files. If a diffusion gradient scheme is present in the input image header, the data provided with this option will be instead used. + +DW gradient table export options +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-export_grad_mrtrix path** export the diffusion-weighted gradient table to file in MRtrix format + +- **-export_grad_fsl bvecs_path bvals_path** export the diffusion-weighted gradient table to files in FSL (bvecs / bvals) format + +Options for exporting additional data regarding PCA behaviour +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-rank image** The signal rank estimated for each denoising patch + +Options for debugging the operation of sliding window kernels +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-max_dist image** The maximum distance between the centre of the patch and a voxel that was included within that patch + +- **-voxelcount image** The number of voxels that contributed to the PCA for processing of each patch + +- **-patchcount image** The number of unique patches to which an input image voxel contributes + +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 + +* If using -estimator med: Gavish, M.; Donoho, D.L. The Optimal Hard Threshold for Singular Values is 4/sqrt(3). IEEE Transactions on Information Theory, 2014, 60(8), 5040-5053. + +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..68dda49771 100644 --- a/docs/reference/commands/dwidenoise.rst +++ b/docs/reference/commands/dwidenoise.rst @@ -21,11 +21,11 @@ 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 noise map estimation and denoising 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. -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. +Important note: image denoising must be performed as the first step of the image processing pipeline. The routine will not operate correctly 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. +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 as part of a complex input image can reduce such non-Gaussian biases. Options ------- diff --git a/docs/reference/commands/dwidenoise2.rst b/docs/reference/commands/dwidenoise2.rst new file mode 100644 index 0000000000..52589ab671 --- /dev/null +++ b/docs/reference/commands/dwidenoise2.rst @@ -0,0 +1,187 @@ +.. _dwidenoise2: + +dwidenoise2 +=================== + +Synopsis +-------- + +Improved dMRI denoising using Marchenko-Pastur PCA + +Usage +-------- + +:: + + dwidenoise2 [ options ] dwi out + +- *dwi*: the input diffusion-weighted image. +- *out*: the output denoised DWI image. + +Description +----------- + +This command performs DWI data denoising, additionally with data-driven noise map estimation if not provided explicitly. The output denoised DWI data is formed based on filtering of eigenvectors of PCA decompositions: for a set of patches each of which consists of a set of proximal voxels, the PCA decomposition is applied, and the DWI signal for those voxels is reconstructed where the contribution of each eigenvector is modulated based on its classification as noise. In many use cases, a threshold that classifies eigenvalues as belonging to signal of interest vs. random thermal noise is based on the prior knowledge that the eigenspectrum of random covariance matrices is described by the universal Marchenko-Pastur (MP) distribution. + +This command includes many capabilities absent from the original dwidenoise command. These include: - Multiple sliding window kernel shapes, including a spherical kernel that dilates at image edges to preserve aspect ratio; - A greater number of mechanisms for noise level estimation, including taking a pre-estimated noise map as input; - Preconditioning, including (per-shell) demeaning, phase demodulation (linear or nonlinear), and variance-stabilising transform to compensate for within-patch heteroscedasticity; - Overcomplete local PCA; - Subsampling (performing fewer PCAs than there are input voxels); - Optimal shrinkage of eigenvalues. + +Important note: image denoising must be performed as the first step of the image processing pipeline. The routine will not operate correctly 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 as part of a complex input image can reduce such non-Gaussian biases. + +If the input data are of complex type, then a smooth non-linear phase will be demodulated removed from each k-space prior to PCA. In the absence of metadata indicating otherwise, it is inferred that the first two axes correspond to acquired slices, and different slices / volumes will be demodulated individually; this behaviour can be modified using the -demod_axes option. A strictly linear phase term can instead be regressed from each k-space, similarly to performed in Cordero-Grande et al. 2019, by specifying -demodulate linear. + +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. + +Permissible sizes for the cuboid kernel depend on the subsampling factor. If no subsampling is performed, or the subsampling factor is odd, then the extent(s) of the kernel must be odd, such that a unique voxel lies at the very centre of each kernel. If however an even subsampling factor is used, then the extent(s) of the kernel must be even, reflecting the fact that it is a voxel corner that resides at the centre of the kernel.In either case, if the extent is specified manually, the user can either provide a single integer---which will determine the number of voxels in the kernel across all three spatial axes---or a comma-separated list of three integers,individually defining the number of voxels in the kernel for all three spatial axes. + +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.Alternatively, optimal truncation as described in Gavish and Donoho 2014 can be utilised by specifying -filter optthresh. + +-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. + +Example usages +-------------- + +- *To approximately replicate the behaviour of the original dwidenoise command*:: + + $ dwidenoise2 DWI.mif out.mif -shape cuboid -subsample 1 -demodulate none -demean none -filter truncate -aggregator exclusive + + While this is neither guaranteed to match exactly the output of the original dwidenoise command nor is it a recommended use case, it may nevertheless be informative in demonstrating those advanced features of dwidenoise2 active by default that must be explicitly disabled in order to approximate that behaviour. + +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. + +Options relating to signal / noise level estimation for denoising +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-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| + * Med: estimate based on the median eigenvalue as in Gavish and Donohue (2014); |br| + * MRM2022: the alternative estimator introduced in Olesen et al. (2022). |br| + Operation will be bypassed if -noise_in or -fixed_rank are specified + +- **-noise_in image** import a pre-estimated noise level map for denoising rather than estimating this level from data + +- **-fixed_rank value** set a fixed input signal rank rather than estimating the noise level from the data + +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 integer or a comma-separated triplet of integers (see Description) + +- **-subsample factor** reduce the number of PCA kernels relative to the number of image voxels; can provide either an integer subsampling factor, or a comma-separated list of three factors; default: 2 + +Options for preconditioning data prior to PCA +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-demodulate mode** select form of phase demodulation; options are: none,linear,nonlinear (default: nonlinear) + +- **-demod_axes axes** comma-separated list of axis indices along which FFT can be applied for phase demodulation + +- **-demean mode** select method of demeaning prior to PCA; options are: none,shells,all (default: 'shells' if DWI gradient table available, 'all' otherwise) + +- **-vst image** apply a within-patch variance-stabilising transformation based on a pre-estimated noise level map + +- **-preconditioned image** export the preconditioned version of the input image that is the input to PCA + +Options that affect reconstruction of the output image series +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-filter choice** Modulate how component contributions are filtered based on the cumulative eigenvalues relative to the noise level; options are: optshrink,optthresh,truncate; default: optshrink (Optimal Shrinkage based on minimisation of the Frobenius norm) + +- **-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 each denoising patch + +- **-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 + +- **-patchcount image** The number of unique patches to which an image voxel contributes + +- **-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 +^^^^^^^^^^^^^^^^ + +- **-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.; Novikov, D.S.; Christiaens, D.; Ades-aron, B.; Sijbers, J. & Fieremans, E. Denoising of diffusion MRI using random matrix theory. NeuroImage, 2016, 142, 394-406, doi: 10.1016/j.neuroimage.2016.08.016 + +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 + +* 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 + +* If using -estimator med or -filter optthresh: Gavish, M.; Donoho, D.L.The Optimal Hard Threshold for Singular Values is 4/sqrt(3). IEEE Transactions on Information Theory, 2014, 60(8), 5040-5053. + +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 Daan Christiaens (daan.christiaens@kcl.ac.uk) and Jelle Veraart (jelle.veraart@nyumc.org) and J-Donald Tournier (jdtournier@gmail.com) + +**Copyright:** 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/. + + diff --git a/docs/reference/commands/mrfilter.rst b/docs/reference/commands/mrfilter.rst index 5a1a7da432..153324171c 100644 --- a/docs/reference/commands/mrfilter.rst +++ b/docs/reference/commands/mrfilter.rst @@ -22,7 +22,7 @@ Usage Description ----------- -The available filters are: fft, gradient, median, smooth, normalise, zclean. +The available filters are: demodulate, fft, gradient, median, smooth, normalise, zclean. Each filter has its own unique set of optional parameters. @@ -31,11 +31,19 @@ For 4D images, each 3D volume is processed independently. Options ------- -Options for FFT filter -^^^^^^^^^^^^^^^^^^^^^^ +Options applicable to demodulate / FFT / k-space filters +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - **-axes list** the axes along which to apply the Fourier Transform. By default, the transform is applied along the three spatial axes. Provide as a comma-separate list of axis indices. +Options applicable to demodulate filter +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-linear** only demodulate based on a linear phase ramp, rather than a filtered k-space + +Options for FFT filter +^^^^^^^^^^^^^^^^^^^^^^ + - **-inverse** apply the inverse FFT - **-magnitude** output a magnitude image rather than a complex-valued image @@ -53,6 +61,13 @@ Options for gradient filter - **-scanner** define the gradient with respect to the scanner coordinate frame of reference. +Options for k-space filtering +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **-window name** specify the shape of the k-space window filter; options are: tukey (no default; must be specified for "kspace" operation) + +- **-strength value** modulate the strength of the chosen filter (exact interpretation & defaultmay depend on the exact filter chosen) + Options for median filter ^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/docs/reference/commands_list.rst b/docs/reference/commands_list.rst index 78cfaa506e..110fa5cd86 100644 --- a/docs/reference/commands_list.rst +++ b/docs/reference/commands_list.rst @@ -34,12 +34,14 @@ List of MRtrix3 commands commands/dwi2adc commands/dwi2fod commands/dwi2mask + commands/dwi2noise commands/dwi2response commands/dwi2tensor commands/dwibiascorrect commands/dwibiasnormmask commands/dwicat commands/dwidenoise + commands/dwidenoise2 commands/dwiextract commands/dwifslpreproc commands/dwigradcheck @@ -167,12 +169,14 @@ 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" |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:`dwidenoise2`, "Improved dMRI denoising using Marchenko-Pastur PCA" |cpp.png|, :ref:`dwiextract`, "Extract diffusion-weighted volumes, b=0 volumes, or certain shells from a DWI dataset" |python.png|, :ref:`dwifslpreproc`, "Perform diffusion image pre-processing using FSL's eddy tool; including inhomogeneity distortion correction using FSL's topup tool if possible" |python.png|, :ref:`dwigradcheck`, "Check the orientation of the diffusion gradient table" diff --git a/src/denoise/denoise.cpp b/src/denoise/denoise.cpp new file mode 100644 index 0000000000..ad93366697 --- /dev/null +++ b/src/denoise/denoise.cpp @@ -0,0 +1,72 @@ +/* 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/denoise.h" + +#include "axes.h" + +namespace MR::Denoise { + +using namespace App; + +const char *first_step_description = + "Important note:" + " image denoising must be performed as the first step of the image processing pipeline." + " The routine will not operate correctly if interpolation or smoothing" + " has been applied to the data prior to denoising."; + +const char *non_gaussian_noise_description = + "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 as part of a complex input image" + " can reduce such non-Gaussian biases."; + +const char *filter_description = + "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." + "Alternatively, optimal truncation as described in Gavish and Donoho 2014 " + "can be utilised by specifying -filter optthresh."; + +const char *aggregation_description = + "-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."; + +const Option datatype_option = Option("datatype", + "Datatype for the eigenvalue decomposition" + " (single or double precision). " + "For complex input data," + " this will select complex float32 or complex float64 datatypes.") + + Argument("float32/float64").type_choice(dtypes); + +} // namespace MR::Denoise diff --git a/src/denoise/denoise.h b/src/denoise/denoise.h new file mode 100644 index 0000000000..fb609e532c --- /dev/null +++ b/src/denoise/denoise.h @@ -0,0 +1,47 @@ +/* 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 "app.h" +#include "header.h" + +namespace MR::Denoise { + +constexpr ssize_t default_subsample_ratio = 2; + +using eigenvalues_type = Eigen::Matrix; +using vector_type = Eigen::Array; + +extern const char *first_step_description; +extern const char *non_gaussian_noise_description; +extern const char *filter_description; +extern const char *aggregation_description; + +const std::vector dtypes = {"float32", "float64"}; +extern const App::Option datatype_option; + +const std::vector filters = {"optshrink", "optthresh", "truncate"}; +enum class filter_type { OPTSHRINK, OPTTHRESH, TRUNCATE }; + +const std::vector aggregators = {"exclusive", "gaussian", "invl0", "rank", "uniform"}; +enum class aggregator_type { EXCLUSIVE, GAUSSIAN, INVL0, RANK, UNIFORM }; + +} // namespace MR::Denoise diff --git a/src/denoise/estimate.cpp b/src/denoise/estimate.cpp new file mode 100644 index 0000000000..9c3b6e4ef3 --- /dev/null +++ b/src/denoise/estimate.cpp @@ -0,0 +1,208 @@ +/* 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 "interp/cubic.h" +#include "math/math.h" + +namespace MR::Denoise { + +template +Estimate::Estimate(const Header &header, + std::shared_ptr subsample, + std::shared_ptr kernel, + std::shared_ptr estimator, + Exports &exports, + const bool enable_recon) + : m(header.size(3)), + subsample(subsample), + kernel(kernel), + estimator(estimator), +#ifdef DWIDENOISE2_USE_BDCSVD + svd_saves_uv(enable_recon), +#endif + X(m, kernel->estimated_size()), +#ifdef DWIDENOISE2_USE_BDCSVD + SVD(m, + kernel->estimated_size(), + svd_saves_uv ? (Eigen::ComputeThinU | Eigen::ComputeThinV) : Eigen::EigenvaluesOnly), +#else + XtX(std::min(m, kernel->estimated_size()), std::min(m, kernel->estimated_size())), + eig(std::min(m, kernel->estimated_size())), +#endif + s(std::min(m, kernel->estimated_size())), + exports(exports) { +} + +template +Estimate::Estimate(const Estimate &that) + : m(that.m), + subsample(that.subsample), + kernel(that.kernel), + estimator(that.estimator), +#ifdef DWIDENOISE2_USE_BDCSVD + svd_saves_uv(that.svd_saves_uv), +#endif + X(m, kernel->estimated_size()), +#ifndef DWIDENOISE2_USE_BDCSVD + XtX(std::min(m, kernel->estimated_size()), std::min(m, kernel->estimated_size())), + eig(std::min(m, kernel->estimated_size())), +#endif + s(std::min(m, kernel->estimated_size())), + exports(that.exports) { +} + +template void Estimate::operator()(Image &dwi) { + + // There are two options here for looping in the presence of subsampling: + // 1. Loop over the input image + // Skip voxels that don't lie at the centre of a patch + // Have to transform input image voxel indices to subsampled image voxel indices for some optional outputs + // 2. Loop over the subsampled image + // In some use cases there may not be any image created that conforms to this voxel grid + // Have to transform the subsampled voxel index into an input image voxel index for the centre of the patch + // Going to go with 1. for now, as for 2. may not have a suitable image over which to loop + Kernel::Voxel::index_type voxel({dwi.index(0), dwi.index(1), dwi.index(2)}); + if (!subsample->process(voxel)) + return; + + // Load list of voxels from which to import data + patch = (*kernel)(voxel); + const ssize_t n = patch.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); + } +#ifndef DWIDENOISE2_USE_BDCSVD + if (r > XtX.cols()) { + DEBUG("Expanding decomposition matrix storage from " + str(X.rows()) + " to " + str(r)); + XtX.resize(r, r); + } +#endif + if (r > s.size()) { + DEBUG("Expanding eigenvalue storage from " + str(s.size()) + " to " + str(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()); +#ifndef DWIDENOISE2_USE_BDCSVD + XtX.fill(std::numeric_limits::signaling_NaN()); +#endif + s.fill(std::numeric_limits::signaling_NaN()); +#endif + + load_data(dwi); + assert(X.leftCols(n).allFinite()); + + // Compute Eigendecomposition +#ifdef DWIDENOISE2_USE_BDCSVD + SVD.compute(X.leftCols(n), svd_saves_uv ? (Eigen::ComputeThinU | Eigen::ComputeThinV) : Eigen::EigenvaluesOnly); + bool successful_decomposition = SVD.info() == Eigen::Success; + if (successful_decomposition) { + // eigenvalues sorted in increasing order: + s.head(r) = SVD.singularValues().array().reverse().square().template cast(); +#else + 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)); + bool successful_decomposition = eig.info() == Eigen::Success; + if (successful_decomposition) { + // eigenvalues sorted in increasing order: + s.head(r) = eig.eigenvalues().template cast(); +#endif + +#ifndef NDEBUG + // Make sure that eigenvalues are in fact sorted; + // may be some instances in which decomposition fails to flag an issue, + // but subsequent threhsold determination is predicated on these data being sorted + double prev = s[0]; + bool sorted = true; + for (ssize_t i = 1; i != r; ++i) { + if (s[i] < prev) { + sorted = false; + break; + } + prev = s[i]; + } + if (!sorted) + throw Exception("Decomposition reports success but eigenvalues are not sorted"); +#endif + + // Threshold determination, possibly via Marchenko-Pastur + threshold = (*estimator)(s, m, n, patch.centre_realspace); + } else { + s.head(r).fill(std::numeric_limits::signaling_NaN()); + threshold = Estimator::Result(); + } + + // Store additional output maps if requested + auto ss_index = subsample->in2ss(voxel); + if (exports.noise_out.valid()) { + assign_pos_of(ss_index).to(exports.noise_out); + exports.noise_out.value() = bool(threshold) // + ? float(std::sqrt(threshold.sigma2)) // + : std::numeric_limits::quiet_NaN(); // + } + if (exports.rank_input.valid()) { + assign_pos_of(ss_index).to(exports.rank_input); + if (!successful_decomposition) + exports.rank_input.value() = 0; + else if (bool(threshold)) + exports.rank_input.value() = r - threshold.cutoff_p; + else + exports.rank_input.value() = r; + } + if (exports.max_dist.valid()) { + assign_pos_of(ss_index).to(exports.max_dist); + exports.max_dist.value() = patch.max_distance; + } + if (exports.voxelcount.valid()) { + assign_pos_of(ss_index).to(exports.voxelcount); + exports.voxelcount.value() = n; + } + if (exports.patchcount.valid()) { + std::lock_guard lock(Estimate::mutex); + for (const auto &v : patch.voxels) { + assign_pos_of(v.index).to(exports.patchcount); + exports.patchcount.value() = exports.patchcount.value() + 1; + } + } +} + +template void Estimate::load_data(Image &image) { + const Kernel::Voxel::index_type pos({image.index(0), image.index(1), image.index(2)}); + for (ssize_t i = 0; i != patch.voxels.size(); ++i) { + assign_pos_of(patch.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/estimate.h b/src/denoise/estimate.h new file mode 100644 index 0000000000..2e53c23038 --- /dev/null +++ b/src/denoise/estimate.h @@ -0,0 +1,115 @@ +/* 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 + +// Switch from SelfAdjointEigenSolver to BDCSVD +#define DWIDENOISE2_USE_BDCSVD + +#include +#include + +#include +#ifdef DWIDENOISE2_USE_BDCSVD +#include +#else +#include +#endif + +#include "denoise/denoise.h" +#include "denoise/estimator/base.h" +#include "denoise/estimator/result.h" +#include "denoise/exports.h" +#include "denoise/kernel/base.h" +#include "denoise/kernel/data.h" +#include "denoise/kernel/voxel.h" +#include "denoise/subsample.h" +#include "header.h" +#include "image.h" +#include "transform.h" + +namespace MR::Denoise { + +template class Estimate { + +public: + using MatrixType = Eigen::Matrix; + + Estimate(const Header &header, + std::shared_ptr subsample, + std::shared_ptr kernel, + std::shared_ptr estimator, + Exports &exports, + const bool enable_recon = false); + + Estimate(const Estimate &); + + void operator()(Image &dwi); + +protected: + const ssize_t m; + + // Denoising configuration + std::shared_ptr subsample; + std::shared_ptr kernel; + std::shared_ptr estimator; +#ifdef DWIDENOISE2_USE_BDCSVD + bool svd_saves_uv; +#endif + + // Reusable memory + Kernel::Data patch; + MatrixType X; + // TODO For both BDCSVD and SelfAdjointEigenSolver, + // the template type is MatrixType, + // and it doesn't seem to be possible to define an Eigen::Block as this template type; + // as such, most likely in both circumstances it is actually constructing a MatrixType from Eigen::Block + // in order to construct the decomposition + // What could conceivably be done instead, + // given that these matrices are relatively small + // and the number of unique patch sizes is small (though not necessarily one), + // would be to construct a std::map<> from patch size to PCA memory; + // each processing thread would allocate new memory for new patch sizes not yet encountered by it, + // but the total memory consumption should still be relatively small; + // note that "X" would be subsumed within such a mechanism also +#ifdef DWIDENOISE2_USE_BDCSVD + Eigen::BDCSVD SVD; +#else + MatrixType XtX; + Eigen::SelfAdjointEigenSolver eig; +#endif + eigenvalues_type s; + Estimator::Result threshold; + + // Export images + // Note: One instance created per thread, + // so that when possible output image data can be written without mutex-locking + Exports exports; + + // Some data can only be written in a thread-safe manner + static std::mutex mutex; + + void load_data(Image &image); +}; + +template std::mutex Estimate::mutex; + +template class Estimate; +template class Estimate; +template class Estimate; +template class Estimate; + +} // namespace MR::Denoise diff --git a/src/denoise/estimator/base.h b/src/denoise/estimator/base.h new file mode 100644 index 0000000000..efa6190a8f --- /dev/null +++ b/src/denoise/estimator/base.h @@ -0,0 +1,34 @@ +/* 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 "denoise/denoise.h" +#include "denoise/estimator/result.h" + +namespace MR::Denoise::Estimator { + +class Base { +public: + Base() = default; + Base(const Base &) = delete; + virtual Result operator()(const eigenvalues_type &eigenvalues, + const ssize_t m, + const ssize_t n, + const Eigen::Vector3d &pos) const = 0; +}; + +} // namespace MR::Denoise::Estimator diff --git a/src/denoise/estimator/estimator.cpp b/src/denoise/estimator/estimator.cpp new file mode 100644 index 0000000000..bf66aaffa5 --- /dev/null +++ b/src/denoise/estimator/estimator.cpp @@ -0,0 +1,90 @@ +/* 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/estimator/estimator.h" + +#include "denoise/estimator/base.h" +#include "denoise/estimator/exp.h" +#include "denoise/estimator/import.h" +#include "denoise/estimator/med.h" +#include "denoise/estimator/mrm2022.h" +#include "denoise/estimator/rank.h" + +namespace MR::Denoise::Estimator { + +using namespace App; + +// clang-format off +const Option estimator_option = + Option("estimator", + "Select the noise level estimator" + " (default = Exp2)," + " either: \n" + "* Exp1: the original estimator used in Veraart et al. (2016); \n" + "* Exp2: the improved estimator introduced in Cordero-Grande et al. (2019); \n" + "* Med: estimate based on the median eigenvalue as in Gavish and Donohue (2014); \n" + "* MRM2022: the alternative estimator introduced in Olesen et al. (2022). \n" + "Operation will be bypassed if -noise_in or -fixed_rank are specified") + + Argument("algorithm").type_choice(estimators); + +const OptionGroup estimator_denoise_options = + OptionGroup("Options relating to signal / noise level estimation for denoising") + + + estimator_option + + + Option("noise_in", + "import a pre-estimated noise level map for denoising rather than estimating this level from data") + + Argument("image").type_image_in() + + + Option("fixed_rank", + "set a fixed input signal rank rather than estimating the noise level from the data") + + Argument("value").type_integer(1); + +std::shared_ptr make_estimator(Image &vst_noise_in, const bool permit_bypass) { + auto opt = get_options("estimator"); + if (permit_bypass) { + auto noise_in = get_options("noise_in"); + auto fixed_rank = get_options("fixed_rank"); + if (!noise_in.empty()) { + if (!opt.empty()) + throw Exception("Cannot both provide an input noise level image and specify a noise level estimator"); + if (!fixed_rank.empty()) + throw Exception("Cannot both provide an input noise level image and request a fixed signal rank"); + return std::make_shared(noise_in[0][0], vst_noise_in); + } + if (!fixed_rank.empty()) { + if (!opt.empty()) + throw Exception("Cannot both provide an input signal rank and specify a noise level estimator"); + return std::make_shared(fixed_rank[0][0]); + } + } + const estimator_type est = opt.empty() ? estimator_type::EXP2 : estimator_type((int)(opt[0][0])); + switch (est) { + case estimator_type::EXP1: + return std::make_shared>(); + case estimator_type::EXP2: + return std::make_shared>(); + case estimator_type::MED: + return std::make_shared(); + case estimator_type::MRM2022: + return std::make_shared(); + default: + assert(false); + } + return nullptr; +} + +} // namespace MR::Denoise::Estimator diff --git a/src/denoise/estimator/estimator.h b/src/denoise/estimator/estimator.h new file mode 100644 index 0000000000..1384b157ee --- /dev/null +++ b/src/denoise/estimator/estimator.h @@ -0,0 +1,36 @@ +/* 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 "app.h" +#include "image.h" + +namespace MR::Denoise::Estimator { + +class Base; + +extern const App::Option estimator_option; +extern const App::OptionGroup estimator_denoise_options; +const std::vector estimators = {"exp1", "exp2", "med", "mrm2022"}; +enum class estimator_type { EXP1, EXP2, MED, MRM2022 }; +std::shared_ptr make_estimator(Image &vst_noise_in, const bool permit_bypass); + +} // namespace MR::Denoise::Estimator diff --git a/src/denoise/estimator/exp.cpp b/src/denoise/estimator/exp.cpp new file mode 100644 index 0000000000..3e09e5bf99 --- /dev/null +++ b/src/denoise/estimator/exp.cpp @@ -0,0 +1,62 @@ +/* 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/estimator/exp.h" + +namespace MR::Denoise::Estimator { + +template +Result Exp::operator()(const eigenvalues_type &s, + const ssize_t m, + const ssize_t n, + const Eigen::Vector3d & /*unused*/) const { + Result result; + const ssize_t r = std::min(m, n); + const ssize_t q = std::max(m, n); + const double lam_r = std::max(s[0], 0.0) / q; + double clam = 0.0; + for (ssize_t p = 0; p < r; ++p) // p+1 is the number of noise components + { // (as opposed to the paper where p is defined as the number of signal components) + const double lam = std::max(s[p], 0.0) / q; + clam += lam; + double denominator = std::numeric_limits::signaling_NaN(); + switch (version) { + case 1: + denominator = q; + break; + case 2: + denominator = q - (r - p - 1); + break; + default: + assert(false); + } + const double gam = double(p + 1) / denominator; + const double sigsq1 = clam / double(p + 1); + const double sigsq2 = (lam - lam_r) / (4.0 * std::sqrt(gam)); + // sigsq2 > sigsq1 if signal else noise + if (sigsq2 < sigsq1) { + result.sigma2 = sigsq1; + result.cutoff_p = p + 1; + result.lamplus = lam; + } + } + return result; +} + +template class Exp<1>; +template class Exp<2>; + +} // namespace MR::Denoise::Estimator diff --git a/src/denoise/estimator/exp.h b/src/denoise/estimator/exp.h new file mode 100644 index 0000000000..2e1815dcdf --- /dev/null +++ b/src/denoise/estimator/exp.h @@ -0,0 +1,34 @@ +/* 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 "denoise/estimator/base.h" +#include "denoise/estimator/result.h" + +namespace MR::Denoise::Estimator { + +template class Exp : public Base { +public: + Exp() {} + ~Exp() {} + Result operator()(const eigenvalues_type &s, + const ssize_t m, + const ssize_t n, + const Eigen::Vector3d & /*unused*/) const final; +}; + +} // namespace MR::Denoise::Estimator diff --git a/src/denoise/estimator/import.h b/src/denoise/estimator/import.h new file mode 100644 index 0000000000..fde0801bba --- /dev/null +++ b/src/denoise/estimator/import.h @@ -0,0 +1,87 @@ +/* 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 "denoise/denoise.h" +#include "denoise/estimator/base.h" +#include "denoise/estimator/result.h" +#include "image.h" +#include "interp/cubic.h" + +namespace MR::Denoise::Estimator { + +class Import : public Base { +public: + Import(const std::string &path, Image &vst_noise_in) // + : noise_image(Image::open(path)), // + vst_noise_image(vst_noise_in) {} // + Result operator()(const eigenvalues_type &s, // + const ssize_t m, // + const ssize_t n, // + const Eigen::Vector3d &pos) const final { // + Result result; + const ssize_t r = std::min(m, n); + const ssize_t q = std::max(m, n); + { + // Construct on each call to preserve const-ness & thread-safety + Interp::Cubic> interp(noise_image); + // TODO This will cause issues at the edge of the image FoV + // Addressing this may require integration of the mrfilter changes + // that provide wrappers for various handling of FoV edges + // For now, just expect that denoising won't do anything + // where the patch centre is too close to the image edge for cubic interpolation + if (!interp.scanner(pos)) + return result; + // If the data have been preconditioned at input based on a pre-estimated noise level, + // then we need to rescale the threshold that we load from this image + // based on knowledge of that rescaling + if (vst_noise_image.valid()) { + Interp::Cubic> vst_interp(vst_noise_image); + if (!vst_interp.scanner(pos)) + return result; + result.sigma2 = Math::pow2(interp.value() / vst_interp.value()); + } else { + result.sigma2 = Math::pow2(interp.value()); + } + } + // From this noise level, + // estimate the upper bound of the MP distribution and rank of signal + // given the ordered list of eigenvalues + double cumulative_lambda = 0.0; + for (ssize_t p = 0; p != r; ++p) { + const double lambda = std::max(s[p], 0.0) / q; + cumulative_lambda += lambda; + const double sigma_sq = cumulative_lambda / (p + 1); + if (sigma_sq < result.sigma2) { + result.cutoff_p = p; + result.lamplus = lambda; + } + } + // TODO It would be nice if the upper bound, lambda_plus, + // could be yielded at a higher precision than the discrete eigenvalues, + // as optimal shrinkage / optimal thresholding could make use of this precision if available + return result; + } + +private: + Image noise_image; + Image vst_noise_image; +}; + +} // namespace MR::Denoise::Estimator diff --git a/src/denoise/estimator/med.h b/src/denoise/estimator/med.h new file mode 100644 index 0000000000..c508fc4964 --- /dev/null +++ b/src/denoise/estimator/med.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 "denoise/estimator/base.h" +#include "denoise/estimator/result.h" +#include "math/math.h" +#include "math/median.h" + +namespace MR::Denoise::Estimator { + +class Med : public Base { +public: + Med() = default; + Result operator()(const eigenvalues_type &s, + const ssize_t m, + const ssize_t n, + const Eigen::Vector3d & /*unused*/) const final { + Result result; + const ssize_t r = std::min(m, n); + const ssize_t q = std::max(m, n); + const ssize_t beta = double(r) / double(q); + // Eigenvalues should already be sorted; + // no need to execute a sort for median calculation + const double ymed = s.size() & 1 ? s[s.size() / 2] : (0.5 * (s[s.size() / 2 - 1] + s[s.size() / 2])); + result.lamplus = ymed / (q * mu(beta)); + // Mechanism intrinsically assumes half rank + result.cutoff_p = r / 2; + // Calculate noise level based on MP distribution + result.sigma2 = 2.0 * s.head(s.size() / 2).sum() / (q * s.size()); + return result; + } + +protected: + // Coefficients as provided in Gavish and Donohue 2014 + // double omega(const double beta) const { + // const double betasq = Math::pow2(beta); + // return (0.56*beta*betasq - 0.95*betasq + 1.82*beta + 1.43); + // } + // Median of Marcenko-Pastur distribution + // Third-order polynomial fit to data generated using Matlab code supplementary to Gavish and Donohue 2014 + double mu(const double beta) const { + const double betasq = Math::pow2(beta); + return ((-0.005882794526340723 * betasq * beta) // + - (0.007508551496715836 * betasq) // + - (0.3338169644754149 * beta) // + + 1.0); // + } +}; + +} // namespace MR::Denoise::Estimator diff --git a/src/denoise/estimator/mrm2022.h b/src/denoise/estimator/mrm2022.h new file mode 100644 index 0000000000..f78f3146ff --- /dev/null +++ b/src/denoise/estimator/mrm2022.h @@ -0,0 +1,62 @@ +/* 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 "denoise/estimator/base.h" +#include "denoise/estimator/result.h" +#include "math/math.h" + +namespace MR::Denoise::Estimator { + +class MRM2022 : public Base { +public: + MRM2022() = default; + Result operator()(const eigenvalues_type &s, + const ssize_t m, + const ssize_t n, + const Eigen::Vector3d & /*unused*/) const final { + Result result; + const ssize_t mprime = std::min(m, n); + const ssize_t nprime = std::max(m, n); + const double sigmasq_to_lamplus = Math::pow2(std::sqrt(nprime) + std::sqrt(mprime)); + double clam = 0.0; + for (ssize_t i = 0; i != mprime; ++i) + clam += std::max(s[i], 0.0); + clam /= nprime; + // Unlike Exp# code, + // MRM2022 article uses p to index number of signal components, + // and here doing a direct translation of the manuscript content to code + double lamplusprev = -std::numeric_limits::infinity(); + for (ssize_t p = 0; p < mprime; ++p) { + const ssize_t i = mprime - 1 - p; + const double lam = std::max(s[i], 0.0) / nprime; + if (lam < lamplusprev) + return result; + clam -= lam; + const double sigmasq = clam / ((mprime - p) * (nprime - p)); + lamplusprev = sigmasq * sigmasq_to_lamplus; + result.cutoff_p = i; + result.sigma2 = sigmasq; + result.lamplus = lamplusprev; + } + return result; + } +}; + +} // namespace MR::Denoise::Estimator diff --git a/src/denoise/estimator/rank.h b/src/denoise/estimator/rank.h new file mode 100644 index 0000000000..69a73f36ca --- /dev/null +++ b/src/denoise/estimator/rank.h @@ -0,0 +1,48 @@ +/* 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 "denoise/estimator/base.h" +#include "denoise/estimator/result.h" + +namespace MR::Denoise::Estimator { + +class Rank : public Base { +public: + Rank(const ssize_t r) : rank(r) {} + Result operator()(const eigenvalues_type &s, + const ssize_t m, + const ssize_t n, + const Eigen::Vector3d & /*unused*/) const final { + Result result; + const ssize_t r = std::min(m, n); + const ssize_t q = std::max(m, n); + result.cutoff_p = r - rank; + double clam = 0.0; + for (ssize_t p = 0; p <= result.cutoff_p; ++p) + clam += std::max(s[p], 0.0); + clam /= q; + result.sigma2 = clam / (result.cutoff_p + 1); + result.lamplus = std::max(s[result.cutoff_p], 0.0); + return result; + } + +protected: + const ssize_t rank; +}; + +} // namespace MR::Denoise::Estimator diff --git a/src/denoise/estimator/result.h b/src/denoise/estimator/result.h new file mode 100644 index 0000000000..3e6423f6bc --- /dev/null +++ b/src/denoise/estimator/result.h @@ -0,0 +1,34 @@ +/* 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 + +namespace MR::Denoise::Estimator { + +class Result { +public: + Result() + : cutoff_p(-1), + sigma2(std::numeric_limits::signaling_NaN()), + lamplus(std::numeric_limits::signaling_NaN()) {} + operator bool() const { return std::isfinite(sigma2); } + bool operator!() const { return !bool(*this); } + ssize_t cutoff_p; + double sigma2; + double lamplus; +}; + +} // namespace MR::Denoise::Estimator diff --git a/src/denoise/exports.cpp b/src/denoise/exports.cpp new file mode 100644 index 0000000000..fa442cfb05 --- /dev/null +++ b/src/denoise/exports.cpp @@ -0,0 +1,69 @@ +/* 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/exports.h" + +#include "dwi/gradient.h" + +namespace MR::Denoise { + +Exports::Exports(const Header &in, const Header &ss) // + : H_in(std::make_shared
(in)), // + H_ss(std::make_shared
(ss)) { // + DWI::clear_DW_scheme(*H_in); + H_in->ndim() = 3; + H_in->reset_intensity_scaling(); + H_in->datatype() = DataType::Float32; + H_in->datatype().set_byte_order_native(); +} + +void Exports::set_noise_out(const std::string &path) { noise_out = Image::create(path, *H_ss); } + +void Exports::set_rank_input(const std::string &path) { + Header H(*H_ss); + H.datatype() = DataType::UInt16; + H.datatype().set_byte_order_native(); + rank_input = Image::create(path, H); +} + +void Exports::set_rank_output(const std::string &path) { rank_output = Image::create(path, *H_in); } + +void Exports::set_sum_optshrink(const std::string &path) { sum_optshrink = Image::create(path, *H_ss); } + +void Exports::set_max_dist(const std::string &path) { max_dist = Image::create(path, *H_ss); } + +void Exports::set_voxelcount(const std::string &path) { + Header H(*H_ss); + H.datatype() = DataType::UInt16; + H.datatype().set_byte_order_native(); + voxelcount = Image::create(path, H); +} + +void Exports::set_patchcount(const std::string &path) { + Header H(*H_in); + H.datatype() = DataType::UInt16; + H.datatype().set_byte_order_native(); + patchcount = Image::create(path, H); +} + +void Exports::set_sum_aggregation(const std::string &path) { + if (path.empty()) + sum_aggregation = Image::scratch(*H_in, "Scratch image for patch aggregation sums"); + else + sum_aggregation = Image::create(path, *H_in); +} + +} // namespace MR::Denoise diff --git a/src/denoise/exports.h b/src/denoise/exports.h new file mode 100644 index 0000000000..250875c335 --- /dev/null +++ b/src/denoise/exports.h @@ -0,0 +1,54 @@ +/* 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 "header.h" +#include "image.h" + +namespace MR::Denoise { + +class Exports { +public: + Exports(const Header &in, const Header &ss); + Exports(const Exports &that) = default; + + void set_noise_out(const std::string &path); + void set_rank_input(const std::string &path); + void set_rank_output(const std::string &path); + void set_sum_optshrink(const std::string &path); + void set_max_dist(const std::string &path); + void set_voxelcount(const std::string &path); + void set_patchcount(const std::string &path); + void set_sum_aggregation(const std::string &path); + + Image noise_out; + Image rank_input; + Image rank_output; + Image sum_optshrink; + Image max_dist; + Image voxelcount; + Image patchcount; + Image sum_aggregation; + +protected: + std::shared_ptr
H_in; + std::shared_ptr
H_ss; +}; + +} // namespace MR::Denoise diff --git a/src/denoise/kernel/base.h b/src/denoise/kernel/base.h new file mode 100644 index 0000000000..490bad437c --- /dev/null +++ b/src/denoise/kernel/base.h @@ -0,0 +1,59 @@ +/* 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 "denoise/kernel/data.h" +#include "denoise/kernel/voxel.h" +#include "header.h" +#include "transform.h" + +namespace MR::Denoise::Kernel { + +class Base { +public: + Base(const Header &H, const std::array &subsample_factors) + : H(H), + transform(H), + halfvoxel_offsets({subsample_factors[0] & 1 ? 0.0 : 0.5, + subsample_factors[1] & 1 ? 0.0 : 0.5, + subsample_factors[2] & 1 ? 0.0 : 0.5}) {} + Base(const Base &) = default; + virtual ~Base() = default; + // This is just for pre-allocating matrices + virtual ssize_t estimated_size() const = 0; + // This is the interface that kernels must provide + virtual Data operator()(const Voxel::index_type &) const = 0; + +protected: + const Header H; + const Transform transform; + std::array halfvoxel_offsets; + + // For translating the index of a processed voxel + // into a realspace position corresponding to the centre of the patch, + // accounting for the fact that subsampling may be introducing an offset + // such that the actual centre of the patch is not at the centre of this voxel + Eigen::Vector3d voxel2real(const Kernel::Voxel::index_type &pos) const { + return ( // + transform.voxel2scanner * // + Eigen::Vector3d({pos[0] + halfvoxel_offsets[0], // + pos[1] + halfvoxel_offsets[1], // + pos[2] + halfvoxel_offsets[2]})); // + } +}; + +} // namespace MR::Denoise::Kernel diff --git a/src/denoise/kernel/cuboid.cpp b/src/denoise/kernel/cuboid.cpp new file mode 100644 index 0000000000..06aefa30d8 --- /dev/null +++ b/src/denoise/kernel/cuboid.cpp @@ -0,0 +1,82 @@ +/* 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/kernel/cuboid.h" + +namespace MR::Denoise::Kernel { + +Cuboid::Cuboid(const Header &header, + const std::array &subsample_factors, + const std::array &extent) + : Base(header, subsample_factors), + size(extent[0] * extent[1] * extent[2]), + // Only sensible if no subsampling is performed, + // and every single DWI voxel is reconstructed from a patch centred at that voxel, + // with no overcomplete local PCA (aggregator == EXCLUSIVE) + centre_index(subsample_factors == std::array({1, 1, 1}) ? (size / 2) : -1) { + for (ssize_t axis = 0; axis != 3; ++axis) { + if (subsample_factors[axis] % 2) { + if (!(extent[axis] % 2)) + throw Exception("For no subsampling / subsampling by an odd number, " + "size of cubic kernel must be an odd integer"); + bounding_box(axis, 0) = -extent[axis] / 2; + bounding_box(axis, 1) = extent[axis] / 2; + } else { + if (extent[axis] % 2) + throw Exception("For subsampling by an even number, " + "size of cubic kernel must be an even integer"); + bounding_box(axis, 0) = 1 - extent[axis] / 2; + bounding_box(axis, 1) = extent[axis] / 2; + } + } +} + +namespace { +// patch handling at image edges +inline ssize_t wrapindex(int p, int r, int bbminus, int bbplus, int max) { + int rr = p + r; + if (rr < 0) + rr = bbplus - r; + if (rr >= max) + rr = (max - 1) + bbminus - r; + return rr; +} +} // namespace + +Data Cuboid::operator()(const Voxel::index_type &pos) const { + Data result(voxel2real(pos), centre_index); + Voxel::index_type voxel; + Offset::index_type offset; + for (offset[2] = bounding_box(2, 0); offset[2] <= bounding_box(2, 1); ++offset[2]) { + voxel[2] = wrapindex(pos[2], offset[2], bounding_box(2, 0), bounding_box(2, 1), H.size(2)); + for (offset[1] = bounding_box(1, 0); offset[1] <= bounding_box(1, 1); ++offset[1]) { + voxel[1] = wrapindex(pos[1], offset[1], bounding_box(1, 0), bounding_box(1, 1), H.size(1)); + for (offset[0] = bounding_box(0, 0); offset[0] <= bounding_box(0, 1); ++offset[0]) { + voxel[0] = wrapindex(pos[0], offset[0], bounding_box(0, 0), bounding_box(0, 1), H.size(0)); + assert(!is_out_of_bounds(H, voxel, 0, 3)); + const default_type sq_distance = Math::pow2(pos[0] + halfvoxel_offsets[0] - voxel[0]) * H.spacing(0) + + Math::pow2(pos[1] + halfvoxel_offsets[1] - voxel[1]) * H.spacing(1) + + Math::pow2(pos[2] + halfvoxel_offsets[2] - voxel[2]) * H.spacing(2); + result.voxels.push_back(Voxel(voxel, sq_distance)); + result.max_distance = std::max(result.max_distance, sq_distance); + } + } + } + result.max_distance = std::sqrt(result.max_distance); + return result; +} + +} // namespace MR::Denoise::Kernel diff --git a/src/denoise/kernel/cuboid.h b/src/denoise/kernel/cuboid.h new file mode 100644 index 0000000000..c90d915efd --- /dev/null +++ b/src/denoise/kernel/cuboid.h @@ -0,0 +1,42 @@ +/* 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 "denoise/kernel/base.h" +#include "denoise/kernel/data.h" +#include "header.h" + +namespace MR::Denoise::Kernel { + +class Cuboid : public Base { + +public: + Cuboid(const Header &header, const std::array &subsample_factors, const std::array &extent); + Cuboid(const Cuboid &) = default; + ~Cuboid() final = default; + Data operator()(const Voxel::index_type &pos) const override; + ssize_t estimated_size() const override { return size; } + +private: + Eigen::Array bounding_box; + const ssize_t size; + const ssize_t centre_index; +}; + +} // namespace MR::Denoise::Kernel diff --git a/src/denoise/kernel/data.h b/src/denoise/kernel/data.h new file mode 100644 index 0000000000..f4659d21a6 --- /dev/null +++ b/src/denoise/kernel/data.h @@ -0,0 +1,52 @@ +/* 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 "denoise/kernel/voxel.h" +#include "types.h" + +namespace MR::Denoise::Kernel { + +class Data { +public: + Data(const Eigen::Vector3d &pos, const ssize_t i) + : centre_realspace(pos), + centre_index(i), + max_distance(-std::numeric_limits::infinity()), + centre_noise(std::numeric_limits::signaling_NaN()) {} + Data(const ssize_t i) + : centre_realspace(Eigen::Vector3d::Constant(std::numeric_limits::signaling_NaN())), + centre_index(i), + max_distance(-std::numeric_limits::infinity()), + centre_noise(std::numeric_limits::signaling_NaN()) {} + Data() + : centre_realspace(Eigen::Vector3d::Constant(std::numeric_limits::signaling_NaN())), + centre_index(-1), + max_distance(-std::numeric_limits::infinity()), + centre_noise(std::numeric_limits::signaling_NaN()) {} + Eigen::Vector3d centre_realspace; + std::vector voxels; + ssize_t centre_index; + default_type max_distance; + double centre_noise; +}; + +} // namespace MR::Denoise::Kernel diff --git a/src/denoise/kernel/kernel.cpp b/src/denoise/kernel/kernel.cpp new file mode 100644 index 0000000000..f935bf13b6 --- /dev/null +++ b/src/denoise/kernel/kernel.cpp @@ -0,0 +1,153 @@ +/* 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/kernel/kernel.h" + +#include "denoise/kernel/base.h" +#include "denoise/kernel/cuboid.h" +#include "denoise/kernel/sphere_radius.h" +#include "denoise/kernel/sphere_ratio.h" +#include "math/math.h" + +namespace MR::Denoise::Kernel { + +using namespace App; + +const char *const shape_description = + "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."; + +const char *const default_size_description = + "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."; + +const char *const cuboid_size_description = + "Permissible sizes for the cuboid kernel depend on the subsampling factor. " + "If no subsampling is performed, or the subsampling factor is odd, " + "then the extent(s) of the kernel must be odd, " + "such that a unique voxel lies at the very centre of each kernel. " + "If however an even subsampling factor is used, " + "then the extent(s) of the kernel must be even, " + "reflecting the fact that it is a voxel corner that resides at the centre of the kernel." + "In either case, if the extent is specified manually, " + "the user can either provide a single integer---" + "which will determine the number of voxels in the kernel across all three spatial axes---" + "or a comma-separated list of three integers," + "individually defining the number of voxels in the kernel for all three spatial axes."; + +// clang-format off +const OptionGroup options = OptionGroup("Options for controlling the sliding spatial window kernel") ++ Option("shape", + "Set the shape of the sliding spatial window. " + "Options are: " + join(shapes, ",") + "; default: sphere") + + Argument("choice").type_choice(shapes) ++ Option("radius_mm", "Set an absolute spherical kernel radius in mm") + + Argument("value").type_float(0.0) ++ Option("radius_ratio", + "Set the spherical kernel size as a ratio of number of voxels to number of input volumes " + "(default: 1.0/0.85 ~= 1.18)") + + Argument("value").type_float(0.0) +// TODO Command-line option that allows user to specify minimum absolute number of voxels in kernel ++ Option("extent", + "Set the patch size of the cuboid kernel; " + "can be either a single integer or a comma-separated triplet of integers (see Description)") + + Argument("window").type_sequence_int(); +// clang-format on + +std::shared_ptr make_kernel(const Header &H, const std::array &subsample_factors) { + auto opt = App::get_options("shape"); + const Kernel::shape_type shape = opt.empty() ? Kernel::shape_type::SPHERE : Kernel::shape_type((int)(opt[0][0])); + std::shared_ptr kernel; + + switch (shape) { + case Kernel::shape_type::SPHERE: { + // TODO Could infer that user wants a cuboid kernel if -extent is used, even if -shape is not + if (!get_options("extent").empty()) + throw Exception("-extent option does not apply to spherical kernel"); + opt = get_options("radius_mm"); + if (opt.empty()) + return std::make_shared( + H, subsample_factors, get_option_value("radius_ratio", sphere_multiplier_default)); + return std::make_shared(H, subsample_factors, opt[0][0]); + } + case Kernel::shape_type::CUBOID: { + if (!get_options("radius_mm").empty() || !get_options("radius_ratio").empty()) + throw Exception("-radius_* options are inapplicable if cuboid kernel shape is selected"); + opt = get_options("extent"); + std::array extent; + if (!opt.empty()) { + auto userinput = parse_ints(opt[0][0]); + if (userinput.size() == 1) + extent = {userinput[0], userinput[0], userinput[0]}; + else if (userinput.size() == 3) + extent = {userinput[0], userinput[1], userinput[2]}; + else + throw Exception("-extent must be either a scalar or a list of length 3"); + for (ssize_t axis = 0; axis < 3; ++axis) { + if (extent[axis] > H.size(axis)) + throw Exception("-extent must not exceed the image dimensions"); + if ((subsample_factors[axis] & 1) != (extent[axis] & 1)) + throw Exception("-extent must match subsampling factors " + "(odd for no subsampling or subsampling by an odd factor; " + "even for subsampling by an even factor)"); + } + } else { + extent = {subsample_factors[0] & 1 ? 3 : 2, subsample_factors[1] & 1 ? 3 : 2, subsample_factors[2] & 1 ? 3 : 2}; + ssize_t prev_num_voxels = 0; // Exit loop below if maximum achievable extent is reached + while (extent[0] * extent[1] * extent[2] < std::max(H.size(3), prev_num_voxels)) { + prev_num_voxels = extent[0] * extent[1] * extent[2]; + // If multiple axes are tied for spatial extent in mm, increment all of them + const default_type min_length = + std::min({extent[0] * H.spacing(0), extent[1] * H.spacing(1), extent[2] * H.spacing(2)}); + for (ssize_t axis = 0; axis != 3; ++axis) { + if (extent[axis] * H.spacing(axis) == min_length && extent[axis] + 2 <= H.size(axis)) + extent[axis] += 2; + } + } + } + INFO("selected cuboid patch size: " + str(extent[0]) + " x " + str(extent[1]) + " x " + str(extent[2])); + + if (std::min(H.size(3), extent[0] * extent[1] * extent[2]) < 15) { + WARN("The number of volumes or the patch size is small; " + "this may lead to discretisation effects in the noise level " + "and cause inconsistent denoising between adjacent voxels"); + } + + return std::make_shared(H, subsample_factors, extent); + } break; + default: + assert(false); + } + return nullptr; +} + +} // namespace MR::Denoise::Kernel diff --git a/src/denoise/kernel/kernel.h b/src/denoise/kernel/kernel.h new file mode 100644 index 0000000000..e091074648 --- /dev/null +++ b/src/denoise/kernel/kernel.h @@ -0,0 +1,42 @@ +/* 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 "app.h" +#include "header.h" +#include "types.h" + +namespace MR::Denoise::Kernel { + +class Base; + +extern const char *const shape_description; +extern const char *const default_size_description; +extern const char *const cuboid_size_description; + +const std::vector shapes = {"cuboid", "sphere"}; +enum class shape_type { CUBOID, SPHERE }; +extern const App::OptionGroup options; + +std::shared_ptr make_kernel(const Header &H, const std::array &subsample_factors); + +} // namespace MR::Denoise::Kernel diff --git a/src/denoise/kernel/sphere_base.cpp b/src/denoise/kernel/sphere_base.cpp new file mode 100644 index 0000000000..11216a994d --- /dev/null +++ b/src/denoise/kernel/sphere_base.cpp @@ -0,0 +1,76 @@ +/* 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/kernel/sphere_base.h" + +#include "math/math.h" + +namespace MR::Denoise::Kernel { + +SphereBase::Shared::Shared(const Header &voxel_grid, + const std::array &subsample_factors, + const std::array &halfvoxel_offsets, + const default_type max_radius) { + const default_type max_radius_sq = Math::pow2(max_radius); + Eigen::Array bounding_box; + for (ssize_t axis = 0; axis != 3; ++axis) { + if (subsample_factors[axis] % 2) { + bounding_box(axis, 1) = int(std::ceil(max_radius / voxel_grid.spacing(axis))); + bounding_box(axis, 0) = -bounding_box(axis, 1); + } else { + bounding_box(axis, 0) = -int(std::ceil((max_radius / voxel_grid.spacing(axis)) - 0.5)); + bounding_box(axis, 1) = int(std::ceil((max_radius / voxel_grid.spacing(axis)) + 0.5)); + } + } + // Build the searchlight + data.reserve(size_t(bounding_box(0, 1) + 1 - bounding_box(0, 0)) * // + size_t(bounding_box(1, 1) + 1 - bounding_box(1, 0)) * // + size_t(bounding_box(2, 1) + 1 - bounding_box(2, 0))); // + Offset::index_type offset({0, 0, 0}); + for (offset[2] = bounding_box(2, 0); offset[2] <= bounding_box(2, 1); ++offset[2]) { + for (offset[1] = bounding_box(1, 0); offset[1] <= bounding_box(1, 1); ++offset[1]) { + for (offset[0] = bounding_box(0, 0); offset[0] <= bounding_box(0, 1); ++offset[0]) { + const default_type squared_distance = + Math::pow2((offset[0] - halfvoxel_offsets[0]) * voxel_grid.spacing(0)) // + + Math::pow2((offset[1] - halfvoxel_offsets[1]) * voxel_grid.spacing(1)) // + + Math::pow2((offset[2] - halfvoxel_offsets[2]) * voxel_grid.spacing(2)); // + if (squared_distance <= max_radius_sq) + data.emplace_back(Offset(offset, squared_distance)); + } + } + } + std::sort(data.begin(), data.end()); + DEBUG("Spherical searchlight construction:"); + DEBUG("Voxel spacing: [" // + + str(voxel_grid.spacing(0)) + "," // + + str(voxel_grid.spacing(1)) + "," // + + str(voxel_grid.spacing(2)) + "]"); // + DEBUG("Maximum nominated radius: " + str(max_radius)); + DEBUG("Halfvoxel offsets: [" // + + str(halfvoxel_offsets[0]) + "," // + + str(halfvoxel_offsets[1]) + "," // + + str(halfvoxel_offsets[2]) + "]"); // + DEBUG("Bounding box for search: [" // + "[" + + str(bounding_box(0, 0)) + " " + str(bounding_box(0, 1)) + "] " + // + "[" + str(bounding_box(1, 0)) + " " + str(bounding_box(1, 1)) + "] " + // + "[" + str(bounding_box(2, 0)) + " " + str(bounding_box(2, 1)) + "]]"); // + DEBUG("First element: " + str(data.front().index.transpose()) + " @ " + str(data.front().distance())); + DEBUG("Last element: " + str(data.back().index.transpose()) + " @ " + str(data.back().distance())); + DEBUG("Number of elements: " + str(data.size())); +} + +} // namespace MR::Denoise::Kernel diff --git a/src/denoise/kernel/sphere_base.h b/src/denoise/kernel/sphere_base.h new file mode 100644 index 0000000000..6d8dc5ace1 --- /dev/null +++ b/src/denoise/kernel/sphere_base.h @@ -0,0 +1,61 @@ +/* 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 "denoise/kernel/base.h" +#include "denoise/kernel/kernel.h" +#include "denoise/kernel/voxel.h" +#include "header.h" + +namespace MR::Denoise::Kernel { + +class SphereBase : public Base { + +public: + SphereBase(const Header &voxel_grid, const std::array &subsample_factors, const default_type max_radius) + : Base(voxel_grid, subsample_factors), + shared(new Shared(voxel_grid, subsample_factors, halfvoxel_offsets, max_radius)), + centre_index(subsample_factors == std::array({1, 1, 1}) ? 0 : -1) {} + + SphereBase(const SphereBase &) = default; + + virtual ~SphereBase() override {} + +protected: + class Shared { + public: + using TableType = std::vector; + Shared(const Header &voxel_grid, + const std::array &subsample_factors, + const std::array &halfvoxel_offsets, + const default_type max_radius); + TableType::const_iterator begin() const { return data.begin(); } + TableType::const_iterator end() const { return data.end(); } + + private: + TableType data; + }; + + std::shared_ptr shared; + const ssize_t centre_index; +}; + +} // namespace MR::Denoise::Kernel diff --git a/src/denoise/kernel/sphere_radius.cpp b/src/denoise/kernel/sphere_radius.cpp new file mode 100644 index 0000000000..be5e4cc1e5 --- /dev/null +++ b/src/denoise/kernel/sphere_radius.cpp @@ -0,0 +1,37 @@ +/* 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/kernel/sphere_radius.h" + +namespace MR::Denoise::Kernel { + +Data SphereFixedRadius::operator()(const Voxel::index_type &pos) const { + Data result(voxel2real(pos), centre_index); + result.voxels.reserve(maximum_size); + for (auto map_it = shared->begin(); map_it != shared->end(); ++map_it) { + const Voxel::index_type voxel({pos[0] + map_it->index[0], // + pos[1] + map_it->index[1], // + pos[2] + map_it->index[2]}); // + if (!is_out_of_bounds(H, voxel, 0, 3)) { + result.voxels.push_back(Voxel(voxel, map_it->sq_distance)); + result.max_distance = map_it->sq_distance; + } + } + result.max_distance = std::sqrt(result.max_distance); + return result; +} + +} // namespace MR::Denoise::Kernel diff --git a/src/denoise/kernel/sphere_radius.h b/src/denoise/kernel/sphere_radius.h new file mode 100644 index 0000000000..1ebb0607c9 --- /dev/null +++ b/src/denoise/kernel/sphere_radius.h @@ -0,0 +1,44 @@ +/* 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 "denoise/kernel/data.h" +#include "denoise/kernel/sphere_base.h" +#include "header.h" +#include "types.h" + +namespace MR::Denoise::Kernel { + +class SphereFixedRadius : public SphereBase { +public: + SphereFixedRadius(const Header &voxel_grid, // + const std::array &subsample_factors, // + const default_type radius) // + : SphereBase(voxel_grid, subsample_factors, radius), // + maximum_size(std::distance(shared->begin(), shared->end())) { // + INFO("Maximum number of voxels in " + str(radius) + "mm fixed-radius kernel is " + str(maximum_size)); + } + SphereFixedRadius(const SphereFixedRadius &) = default; + ~SphereFixedRadius() final = default; + Data operator()(const Voxel::index_type &pos) const; + ssize_t estimated_size() const override { return maximum_size; } + +private: + const ssize_t maximum_size; +}; + +} // namespace MR::Denoise::Kernel diff --git a/src/denoise/kernel/sphere_ratio.cpp b/src/denoise/kernel/sphere_ratio.cpp new file mode 100644 index 0000000000..f4532fc002 --- /dev/null +++ b/src/denoise/kernel/sphere_ratio.cpp @@ -0,0 +1,62 @@ +/* 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/kernel/sphere_ratio.h" + +namespace MR::Denoise::Kernel { + +Data SphereRatio::operator()(const Voxel::index_type &pos) const { + Data result(voxel2real(pos), centre_index); + auto table_it = shared->begin(); + while (table_it != shared->end()) { + // If there's a tie in distances, want to include all such offsets in the kernel, + // even if the size of the utilised kernel extends beyond the minimum size + if (result.voxels.size() >= min_size && table_it->sq_distance != result.max_distance) + break; + const Voxel::index_type voxel({pos[0] + table_it->index[0], // + pos[1] + table_it->index[1], // + pos[2] + table_it->index[2]}); // + if (!is_out_of_bounds(H, voxel, 0, 3)) { + result.voxels.push_back(Voxel(voxel, table_it->sq_distance)); + result.max_distance = table_it->sq_distance; + } + ++table_it; + } + if (table_it == shared->end()) { + throw Exception( // + std::string("Inadequate spherical kernel initialisation ") // + + "(lookup table " + str(std::distance(shared->begin(), shared->end())) + "; " // + + "min size " + str(min_size) + "; " // + + "read size " + str(result.voxels.size()) + ")"); // + } + result.max_distance = std::sqrt(result.max_distance); + return result; +} + +default_type SphereRatio::compute_max_radius(const Header &voxel_grid, const default_type min_ratio) const { + const size_t num_volumes = voxel_grid.size(3); + const default_type voxel_volume = voxel_grid.spacing(0) * voxel_grid.spacing(1) * voxel_grid.spacing(2); + const default_type sphere_volume = 8.0 * num_volumes * min_ratio * voxel_volume; + const default_type approx_radius = std::sqrt(sphere_volume * 0.75 / Math::pi); + const Voxel::index_type half_extents({ssize_t(std::ceil(approx_radius / voxel_grid.spacing(0))), // + ssize_t(std::ceil(approx_radius / voxel_grid.spacing(1))), // + ssize_t(std::ceil(approx_radius / voxel_grid.spacing(2)))}); // + return std::max({half_extents[0] * voxel_grid.spacing(0), + half_extents[1] * voxel_grid.spacing(1), + half_extents[2] * voxel_grid.spacing(2)}); +} + +} // namespace MR::Denoise::Kernel diff --git a/src/denoise/kernel/sphere_ratio.h b/src/denoise/kernel/sphere_ratio.h new file mode 100644 index 0000000000..36778ac2e4 --- /dev/null +++ b/src/denoise/kernel/sphere_ratio.h @@ -0,0 +1,51 @@ +/* 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 "denoise/kernel/data.h" +#include "denoise/kernel/sphere_base.h" +#include "header.h" + +namespace MR::Denoise::Kernel { + +constexpr default_type sphere_multiplier_default = 1.0 / 0.85; + +class SphereRatio : public SphereBase { + +public: + SphereRatio(const Header &voxel_grid, const std::array &subsample_factors, const default_type min_ratio) + : SphereBase(voxel_grid, subsample_factors, compute_max_radius(voxel_grid, min_ratio)), + min_size(std::ceil(voxel_grid.size(3) * min_ratio)) {} + + SphereRatio(const SphereRatio &) = default; + + ~SphereRatio() final = default; + + Data operator()(const Voxel::index_type &pos) const override; + + ssize_t estimated_size() const override { return min_size; } + +private: + ssize_t min_size; + + // Determine an appropriate bounding box from which to generate the search table + // Find the radius for which 7/8 of the sphere will contain the minimum number of voxels, then round up + // This is only for setting the maximal radius for generation of the lookup table + default_type compute_max_radius(const Header &voxel_grid, const default_type min_ratio) const; +}; + +} // namespace MR::Denoise::Kernel diff --git a/src/denoise/kernel/voxel.h b/src/denoise/kernel/voxel.h new file mode 100644 index 0000000000..0570d2c8d7 --- /dev/null +++ b/src/denoise/kernel/voxel.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 "types.h" + +namespace MR::Denoise::Kernel { + +template class VoxelBase { +public: + using index_type = Eigen::Array; + VoxelBase(const index_type &index, const default_type sq_distance) : index(index), sq_distance(sq_distance) {} + VoxelBase(const VoxelBase &) = default; + VoxelBase(VoxelBase &&) = default; + ~VoxelBase() {} + VoxelBase &operator=(const VoxelBase &that) { + index = that.index; + sq_distance = that.sq_distance; + return *this; + } + VoxelBase &operator=(VoxelBase &&that) noexcept { + index = that.index; + sq_distance = that.sq_distance; + return *this; + } + bool operator<(const VoxelBase &that) const { return sq_distance < that.sq_distance; } + default_type distance() const { return std::sqrt(sq_distance); } + + index_type index; + default_type sq_distance; +}; + +// Need signed integer to represent offsets from the centre of the kernel; +// however absolute voxel indices should be unsigned +// For nonstationarity correction, +// "Voxel" also needs a noise level estimate per voxel in the patch, +// as the scaling on insertion of data into the matrix +// needs to be reversed when reconstructing the denoised signal from the eigenvectors +class Voxel : public VoxelBase { +public: + using index_type = VoxelBase::index_type; + Voxel(const index_type &index, const default_type sq_distance, const default_type noise_level) + : VoxelBase(index, sq_distance), noise_level(noise_level) {} + Voxel(const index_type &index, const default_type sq_distance) + : VoxelBase(index, sq_distance), noise_level(std::numeric_limits::signaling_NaN()) {} + default_type noise_level; +}; + +using Offset = VoxelBase; + +} // namespace MR::Denoise::Kernel diff --git a/src/denoise/precondition.cpp b/src/denoise/precondition.cpp new file mode 100644 index 0000000000..7ba883bf69 --- /dev/null +++ b/src/denoise/precondition.cpp @@ -0,0 +1,526 @@ +/* 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/precondition.h" + +#include + +#include "algo/copy.h" +#include "app.h" +#include "axes.h" +#include "dwi/gradient.h" +#include "dwi/shells.h" +#include "transform.h" + +using namespace MR::App; + +namespace MR::Denoise { + +const char *const demodulation_description = + "If the input data are of complex type, " + "then a smooth non-linear phase will be demodulated removed from each k-space prior to PCA. " + "In the absence of metadata indicating otherwise, " + "it is inferred that the first two axes correspond to acquired slices, " + "and different slices / volumes will be demodulated individually; " + "this behaviour can be modified using the -demod_axes option. " + "A strictly linear phase term can instead be regressed from each k-space, " + "similarly to performed in Cordero-Grande et al. 2019, " + "by specifying -demodulate linear."; + +// clang-format off +OptionGroup precondition_options(const bool include_output) +{ + OptionGroup result ("Options for preconditioning data prior to PCA"); + result + + Option("demodulate", + "select form of phase demodulation; " + "options are: " + join(demodulation_choices, ",") + " " + "(default: nonlinear)") + + Argument("mode").type_choice(demodulation_choices); + + Option("demod_axes", + "comma-separated list of axis indices along which FFT can be applied for phase demodulation") + + Argument("axes").type_sequence_int(); + + Option("demean", + "select method of demeaning prior to PCA; " + "options are: " + join(demean_choices, ",") + " " + "(default: 'shells' if DWI gradient table available; 'volume_groups' if volume groups present; 'all' otherwise)") + + Argument("mode").type_choice(demean_choices); + + Option("vst", + "apply a within-patch variance-stabilising transformation based on a pre-estimated noise level map") + + Argument("image").type_image_in(); + if (include_output) { + result + + Option("preconditioned_input", + "export the preconditioned version of the input image that is the input to PCA") + + Argument("image").type_image_out() + + Option("preconditioned_output", + "export the denoised data prior to reversal of preconditioning") + + Argument("image").type_image_out(); + } else { + result + + Option("preconditioned", + "export the preconditioned version of the input image that is the input to PCA") + + Argument("image").type_image_out(); + } + return result; +} +// clang-format on + +Demodulation select_demodulation(const Header &H) { + const bool complex = H.datatype().is_complex(); + auto opt_mode = get_options("demodulate"); + auto opt_axes = get_options("demod_axes"); + Demodulation result; + if (opt_mode.empty()) { + if (complex) { + result.mode = demodulation_t::NONLINEAR; + } else { + if (!opt_axes.empty()) { + throw Exception("Option -demod_axes cannot be specified: " + "no phase demodulation of magnitude data"); + } + } + } else { + result.mode = demodulation_t(int(opt_mode[0][0])); + if (!complex) { + switch (result.mode) { + case demodulation_t::NONE: + WARN("Specifying -demodulate none is redundant: " + "never any phase demodulation for magnitude input data"); + break; + default: + throw Exception("Phase modulation cannot be utilised for magnitude-only input data"); + } + } + } + if (!complex) + return result; + if (opt_axes.empty()) { + auto slice_encoding_it = H.keyval().find("SliceEncodingDirection"); + if (slice_encoding_it == H.keyval().end()) { + // TODO Ideally this would be the first two axes *on disk*, + // not following transform realignment + INFO("No header information on slice encoding; " + "assuming first two axes are within-slice"); + result.axes = {0, 1}; + } else { + auto dir = Axes::id2dir(slice_encoding_it->second); + for (size_t axis = 0; axis != 3; ++axis) { + if (!dir[axis]) + result.axes.push_back(axis); + } + INFO("For header SliceEncodingDirection=\"" + slice_encoding_it->second + + "\", " + "chose demodulation axes: " + + join(result.axes, ",")); + } + } else { + result.axes = parse_ints(opt_axes[0][0]); + for (auto axis : result.axes) { + if (axis > 2) + throw Exception("Phase demodulation implementation not yet robust to non-spatial axes"); + } + } + return result; +} + +demean_type select_demean(const Header &H) { + bool shells_available = false; + try { + auto grad = DWI::get_DW_scheme(H); + auto shells = DWI::Shells(grad); + shells_available = true; + } catch (Exception &) { + } + const bool volume_groups_available = H.ndim() > 4; + auto opt = get_options("demean"); + if (opt.empty()) { + if (shells_available) { + if (volume_groups_available) + throw Exception("Cannot automatically determine how to demean, " + "as input series is >4D and also has a diffusion gradient table"); + INFO("Automatically demeaning per b-value shell based on input gradient table"); + return demean_type::SHELLS; + } + if (volume_groups_available) { + INFO("Automatically demeaning by volume groups"); + return demean_type::VOLUME_GROUPS; + } + INFO("Automatically demeaning across all volumes"); + return demean_type::ALL; + } + const demean_type user_selection = demean_type(int(opt[0][0])); + if (user_selection == demean_type::SHELLS && !shells_available) + throw Exception("Cannot demean by b-value shells as shell structure could not be inferred"); + if (user_selection == demean_type::VOLUME_GROUPS && !volume_groups_available) + throw Exception("Cannot demean by volume groups as image does not possess volume groups"); + return user_selection; +} + +template +Precondition::Precondition(Image &image, + const Demodulation &demodulation, + const demean_type demean, + Image &vst_image) + : H_in(image), // + H_out(image), // + num_volume_groups(1), // + vst_image(vst_image) { // + + for (ssize_t axis = 4; axis != H_in.ndim(); ++axis) { + num_volume_groups *= H_in.size(axis); + H_out.size(3) *= H_in.size(axis); + } + H_out.ndim() = 4; + Stride::set(H_out, Stride::contiguous_along_axis(3)); + H_out.datatype() = DataType::from(); + H_out.datatype().set_byte_order_native(); + + if (H_in.ndim() > 4) { + Header H_serialise(H_in); + for (ssize_t axis = 3; axis != H_in.ndim(); ++axis) { + H_serialise.size(axis - 3) = H_in.size(axis); + H_serialise.stride(axis - 3) = axis - 2; + H_serialise.spacing(axis - 3) = 1.0; + } + for (ssize_t axis = H_in.ndim() - 3; axis != 3; ++axis) { + H_serialise.size(axis) = 1; + H_serialise.stride(axis) = axis + 1; + H_serialise.spacing(axis) = 1.0; + } + H_serialise.ndim() = std::max(size_t(3), H_in.ndim() - 3); + H_serialise.datatype() = DataType::from(); + H_serialise.datatype().set_byte_order_native(); + H_serialise.transform().setIdentity(); + serialise_image = Image::scratch( // + H_serialise, // + "Scratch image for serialising non-spatial indices into Casorati matrix"); // + + ssize_t output_index = 0; + for (auto l = Loop(serialise_image)(serialise_image); l; ++l) + serialise_image.value() = output_index++; + serialise_image.reset(); + assert(output_index == H_out.size(3)); + } + + // Step 1: Phase demodulation + Image dephased; + if (demodulation.mode == demodulation_t::NONE) { + dephased = image; + } else { + typename DemodulatorSelector::type demodulator(image, // + demodulation.axes, // + demodulation.mode == demodulation_t::LINEAR); // + phase_image = demodulator(); + // Only actually perform the dephasing of the input image within this constructor + // if that result needs to be utilised in calculation of the mean + if (demean != demean_type::NONE) { + dephased = Image::scratch(H_in, "Scratch dephased version of \"" + image.name() + "\" for mean calculation"); + demodulator(image, dephased, false); + } + } + + // Step 2: Demeaning + Header H_mean(H_out); + switch (demean) { + case demean_type::NONE: + break; + case demean_type::VOLUME_GROUPS: { + assert(serialise_image.valid()); + if (H_in.ndim() < 5) + throw Exception("Cannot demean by volume groups if input image is <= 4D"); + index2group.resize(H_out.size(3)); + ssize_t group_index = 0; + for (auto l_group = Loop(serialise_image, 1)(serialise_image); l_group; ++l_group, ++group_index) { + for (auto l_volumes = Loop(serialise_image, 0, 1)(serialise_image); l_volumes; ++l_volumes) + index2group[serialise_image.value()] = group_index; + } + serialise_image.reset(); + assert(group_index == num_volume_groups); + H_mean.size(3) = num_volume_groups; + mean_image = Image::scratch(H_mean, "Scratch image for per-volume-group mean intensity"); + for (auto l_outer = Loop("Computing mean intensity image per volume group", H_in, 3)(dephased); // + l_outer; // + ++l_outer) { // + // Which volume group does this volume belong to? + for (ssize_t axis = 3; axis != H_in.ndim(); ++axis) + serialise_image.index(axis - 3) = dephased.index(axis); + mean_image.index(3) = index2group[serialise_image.value()]; + // Add the values within this volume to the mean intensity of the respective volume group + for (auto l_voxel = Loop(H_in, 0, 3)(dephased, mean_image); l_voxel; ++l_voxel) + mean_image.value() += dephased.value(); + } + const default_type multiplier = 1.0 / H_in.size(3); + for (auto l = Loop(H_mean)(mean_image); l; ++l) + mean_image.value() *= multiplier; + } break; + case demean_type::SHELLS: { + Eigen::Matrix grad; + try { + grad = DWI::get_DW_scheme(H_mean); + } catch (Exception &e) { + throw Exception(e, "Cannot demean by shells as unable to obtain valid gradient table"); + } + try { + DWI::Shells shells(grad); + index2shell.resize(image.size(3), -1); + for (ssize_t shell_idx = 0; shell_idx != shells.count(); ++shell_idx) { + for (auto v : shells[shell_idx].get_volumes()) + index2shell[v] = shell_idx; + } + assert(*std::min_element(index2shell.begin(), index2shell.end()) == 0); + H_mean.size(3) = shells.count(); + DWI::stash_DW_scheme(H_mean, grad); + mean_image = Image::scratch(H_mean, "Scratch image for per-shell mean intensity"); + for (auto l_voxel = Loop("Computing mean intensities within shells", H_mean, 0, 3)(dephased, mean_image); // + l_voxel; // + ++l_voxel) { // + for (ssize_t volume_idx = 0; volume_idx != H_in.size(3); ++volume_idx) { + dephased.index(3) = volume_idx; + mean_image.index(3) = index2shell[volume_idx]; + mean_image.value() += dephased.value(); + } + for (ssize_t shell_idx = 0; shell_idx != shells.count(); ++shell_idx) { + mean_image.index(3) = shell_idx; + mean_image.value() /= T(shells[shell_idx].count()); + } + } + } catch (Exception &e) { + throw Exception(e, "Cannot demean by shells as unable to establish b-value shell structure"); + } + } break; + case demean_type::ALL: { + H_mean.ndim() = 3; + DWI::clear_DW_scheme(H_mean); + mean_image = Image::scratch(H_mean, "Scratch image for mean intensity across all volumes"); + const T multiplier = T(1) / T(H_out.size(3)); + for (auto l_voxel = Loop("Computing mean intensity across all volumes", H_mean)(dephased, mean_image); // + l_voxel; // + ++l_voxel) { // + T mean(T(0)); + for (auto l_volume = Loop(3)(dephased); l_volume; ++l_volume) + mean += T(dephased.value()); + mean_image.value() = multiplier * mean; + } + } break; + } + + // Step 3: Variance-stabilising transform + // Image vst is already set within constructor definition; + // no preparation work to do here +} + +namespace { +// Private functions to prevent compiler attempting to create complex functions for real types +template +typename std::enable_if::value, T>::type demodulate(const cfloat in, const cfloat phase) { + return in * std::conj(phase); +} +template +typename std::enable_if::value, T>::type demodulate(const cdouble in, const cfloat phase) { + return in * std::conj(cdouble(phase)); +} +template +typename std::enable_if::value, T>::type demodulate(const T in, const cfloat phase) { + assert(false); + return in; +} +template +typename std::enable_if::value, T>::type modulate(const cfloat in, const cfloat phase) { + return in * phase; +} +template +typename std::enable_if::value, T>::type modulate(const cdouble in, const cfloat phase) { + return in * cdouble(phase); +} +template typename std::enable_if::value, T>::type modulate(const T in, const cfloat phase) { + assert(false); + return in; +} +} // namespace + +template void Precondition::operator()(Image input, Image output, const bool inverse) const { + + // For thread-safety / const-ness + const Transform transform(input); + Image serialise(serialise_image); + Image phase(phase_image); + Image mean(mean_image); + std::unique_ptr>> vst; + if (vst_image.valid()) + vst.reset(new Interp::Cubic>(vst_image)); + + Eigen::Array data(H_out.size(3)); + if (inverse) { + + assert(dimensions_match(H_out, input)); + assert(dimensions_match(H_in, output)); + + for (auto l_voxel = Loop("Reversing data preconditioning", H_in, 0, 3)(input, output); l_voxel; ++l_voxel) { + + // Step 3: Reverse variance-stabilising transform + if (vst) { + vst->scanner(transform.voxel2scanner * // + Eigen::Vector3d({default_type(input.index(0)), // + default_type(input.index(1)), // + default_type(input.index(2))})); // + const T multiplier = T(vst->value()); + for (ssize_t v = 0; v != H_out.size(3); ++v) { + input.index(3) = v; + data[v] = T(input.value()) * multiplier; + } + } else { + for (ssize_t v = 0; v != H_out.size(3); ++v) { + input.index(3) = v; + data[v] = input.value(); + } + } + + // Step 2: Reverse demeaning + if (mean.valid()) { + assign_pos_of(input, 0, 3).to(mean); + if (mean.ndim() == 3) { + const T mean_value = mean.value(); + data += mean_value; + } else if (!index2shell.empty()) { + for (ssize_t v = 0; v != H_out.size(3); ++v) { + mean.index(3) = index2shell[v]; + data[v] += T(mean.value()); + } + } else if (!index2group.empty()) { + for (ssize_t v = 0; v != H_out.size(3); ++v) { + mean.index(3) = index2group[v]; + data[v] += T(mean.value()); + } + } else { + assert(false); + data.fill(std::numeric_limits::signaling_NaN()); + } + } + + // Step 1: Reverse phase demodulation + if (phase.valid()) { + assign_pos_of(input, 0, 3).to(phase); + if (serialise.valid()) { + for (auto l = Loop(H_in, 3)(phase); l; ++l) { + for (ssize_t axis = 3; axis != H_in.ndim(); ++axis) + serialise.index(axis - 3) = phase.index(axis); + data[serialise.value()] = modulate(data[serialise.value()], phase.value()); + } + } else { + for (ssize_t v = 0; v != H_out.size(3); ++v) { + phase.index(3) = v; + data[v] = modulate(data[v], phase.value()); + } + } + } + + // Write to output + if (serialise.valid()) { + for (auto l = Loop(H_in, 3)(output); l; ++l) { + for (ssize_t axis = 3; axis != H_in.ndim(); ++axis) + serialise.index(axis - 3) = output.index(axis); + output.value() = data[serialise.value()]; + } + } else { + for (ssize_t v = 0; v != H_out.size(3); ++v) { + output.index(3) = v; + output.value() = data[v]; + } + } + } + return; + } + + assert(dimensions_match(H_in, input)); + assert(dimensions_match(H_out, output)); + + // Applying forward preconditioning + for (auto l_voxel = Loop("Applying data preconditioning", H_in, 0, 3)(input, output); l_voxel; ++l_voxel) { + + // Serialise all data within this voxel into "data" + if (H_in.ndim() == 4) { + for (ssize_t v = 0; v != H_out.size(3); ++v) { + input.index(3) = v; + data[v] = input.value(); + } + } else { + for (auto l = Loop(H_in, 3)(input); l; ++l) { + for (ssize_t axis = 3; axis != H_in.ndim(); ++axis) + serialise.index(axis - 3) = input.index(axis); + data[serialise.value()] = input.value(); + } + } + + // Step 1: Phase demodulation + if (phase.valid()) { + assign_pos_of(input, 0, 3).to(phase); + if (H_in.ndim() == 4) { + for (ssize_t v = 0; v != H_out.size(3); ++v) { + phase.index(3) = v; + data[v] = demodulate(data[v], phase.value()); + } + } else { + for (auto l = Loop(H_in, 3)(phase); l; ++l) { + for (ssize_t axis = 3; axis != H_in.ndim(); ++axis) + serialise.index(axis - 3) = phase.index(axis); + data[serialise.value()] = demodulate(data[serialise.value()], phase.value()); + } + } + } + + // Step 2: Demeaning + if (mean.valid()) { + assign_pos_of(input, 0, 3).to(mean); + if (mean.ndim() == 3) { + const T mean_value = mean.value(); + for (ssize_t v = 0; v != H_out.size(3); ++v) + data[v] -= mean_value; + } else if (!index2shell.empty()) { + for (ssize_t v = 0; v != H_out.size(3); ++v) { + mean.index(3) = index2shell[v]; + data[v] -= T(mean.value()); + } + } else if (!index2group.empty()) { + for (ssize_t v = 0; v != H_out.size(3); ++v) { + mean.index(3) = index2group[v]; + data[v] -= T(mean.value()); + } + } else { + assert(false); + data.fill(std::numeric_limits::signaling_NaN()); + } + } + + // Step 3: Variance-stabilising transform + if (vst) { + vst->scanner(transform.voxel2scanner // + * Eigen::Vector3d({default_type(input.index(0)), // + default_type(input.index(1)), // + default_type(input.index(2))})); // + const default_type multiplier = 1.0 / vst->value(); + data *= multiplier; + } + + // Write to output + for (ssize_t v = 0; v != H_out.size(3); ++v) { + output.index(3) = v; + output.value() = data[v]; + } + } +} + +} // namespace MR::Denoise diff --git a/src/denoise/precondition.h b/src/denoise/precondition.h new file mode 100644 index 0000000000..4f179e17f3 --- /dev/null +++ b/src/denoise/precondition.h @@ -0,0 +1,101 @@ +/* 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 "app.h" +#include "denoise/kernel/voxel.h" +#include "filter/demodulate.h" +#include "header.h" +#include "image.h" +#include "interp/cubic.h" +#include "types.h" + +namespace MR::Denoise { + +extern const char *const demodulation_description; + +const std::vector demodulation_choices({"none", "linear", "nonlinear"}); +enum class demodulation_t { NONE, LINEAR, NONLINEAR }; + +const std::vector demean_choices = {"none", "volume_groups", "shells", "all"}; +enum class demean_type { NONE, VOLUME_GROUPS, SHELLS, ALL }; + +App::OptionGroup precondition_options(const bool include_output); + +class Demodulation { +public: + Demodulation(demodulation_t mode) : mode(mode) {} + Demodulation() : mode(demodulation_t::NONE) {} + explicit operator bool() const { return mode != demodulation_t::NONE; } + bool operator!() const { return mode == demodulation_t::NONE; } + demodulation_t mode; + std::vector axes; +}; +Demodulation select_demodulation(const Header &); + +demean_type select_demean(const Header &); + +// Need to SFINAE define the demodulator type, +// so that it does not attempt to compile the demodulation filter for non-complex types +class DummyDemodulator { +public: + template DummyDemodulator(ImageType &, const std::vector &, const bool) {} + template + void operator()(InputImageType &, OutputImageType &, const bool) { + assert(false); + } + Image operator()() { return Image(); } +}; +template struct DemodulatorSelector { + using type = DummyDemodulator; +}; +template struct DemodulatorSelector> { + using type = Filter::Demodulate; +}; + +template class Precondition { +public: + Precondition(Image &image, const Demodulation &demodulation, const demean_type demean, Image &vst); + Precondition(Precondition &) = default; + void operator()(Image input, Image output, const bool inverse = false) const; + ssize_t rank() const { return phase_image.valid() || mean_image.valid() ? 1 : 0; } + const Header &header() const { return H_out; } + +private: + const Header H_in; + Header H_out; + // For serialisation of >4D images + ssize_t num_volume_groups; + Image serialise_image; + // First step: Phase demodulation + Image phase_image; + // Second step: Demeaning + std::vector index2shell; + std::vector index2group; + Image mean_image; + // Third step: Variance-stabilising transform + Image vst_image; +}; +template class Precondition; +template class Precondition; +template class Precondition; +template class Precondition; + +} // namespace MR::Denoise diff --git a/src/denoise/recon.cpp b/src/denoise/recon.cpp new file mode 100644 index 0000000000..36b65d21e6 --- /dev/null +++ b/src/denoise/recon.cpp @@ -0,0 +1,260 @@ +/* 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, + std::shared_ptr subsample, + std::shared_ptr kernel, + std::shared_ptr estimator, + filter_type filter, + aggregator_type aggregator, + Exports &exports) + : Estimate(header, subsample, kernel, estimator, exports, true), + filter(filter), + aggregator(aggregator), + // FWHM = 2 x cube root of spacings between kernels + gaussian_multiplier(-std::log(2.0) / // + Math::pow2(std::cbrt(subsample->get_factors()[0] * header.spacing(0) // + * subsample->get_factors()[1] * header.spacing(1) // + * subsample->get_factors()[2] * header.spacing(2)))), // + w(std::min(Estimate::m, kernel->estimated_size())), + Xr(Estimate::m, aggregator == aggregator_type::EXCLUSIVE ? 1 : kernel->estimated_size()) {} + +template void Recon::operator()(Image &dwi, Image &out) { + + if (!Estimate::subsample->process({dwi.index(0), dwi.index(1), dwi.index(2)})) + return; + + Estimate::operator()(dwi); + + const ssize_t n = Estimate::patch.voxels.size(); + const ssize_t r = std::min(Estimate::m, n); + const ssize_t q = std::max(Estimate::m, n); + const double beta = double(r) / double(q); + const ssize_t in_rank = bool(Estimate::threshold) // + ? (r - Estimate::threshold.cutoff_p) // + : -1; // + + if (r > w.size()) + w.resize(r); + if (aggregator != aggregator_type::EXCLUSIVE && n > Xr.cols()) + Xr.resize(Estimate::m, n); +#ifndef NDEBUG + w.fill(std::numeric_limits::signaling_NaN()); + Xr.fill(std::numeric_limits::signaling_NaN()); +#endif + + // Generate weights vector + double sum_weights = 0.0; + ssize_t out_rank = 0; + if (bool(Estimate::threshold)) { + switch (filter) { + case filter_type::OPTSHRINK: { + const double transition = 1.0 + std::sqrt(beta); + for (ssize_t i = 0; i != r; ++i) { + const double lam = std::max(Estimate::s[i], 0.0) / q; + const double y = lam / Estimate::threshold.sigma2; + double nu = 0.0; + if (y > transition) { + // Occasionally floating-point precision will drive this calculation to fractionally greater than y, + // which will erroneously yield a weight fractionally greater than 1.0 + nu = std::min(y, std::sqrt(Math::pow2(Math::pow2(y) - beta - 1.0) - (4.0 * beta)) / y); + ++out_rank; + } + w[i] = lam > 0.0 ? (nu / y) : 0.0; + assert(w[i] >= 0.0 && w[i] <= 1.0); + sum_weights += w[i]; + } + } break; + case filter_type::OPTTHRESH: { + const std::map::const_iterator it = beta2lambdastar.find(beta); + double lambda_star = 0.0; + if (it == beta2lambdastar.end()) { + lambda_star = + sqrt(2.0 * (beta + 1.0) + ((8.0 * beta) / (beta + 1.0 + std::sqrt(Math::pow2(beta) + 14.0 * beta + 1.0)))); + beta2lambdastar[beta] = lambda_star; + } else { + lambda_star = it->second; + } + const double tau_star = lambda_star * std::sqrt(q) * std::sqrt(Estimate::threshold.sigma2); + // TODO Unexpected requisite square applied to q here + const double threshold = tau_star * Math::pow2(q); + for (ssize_t i = 0; i != r; ++i) { + if (Estimate::s[i] >= threshold) { + w[i] = 1.0; + ++out_rank; + } else { + w[i] = 0.0; + } + } + sum_weights = out_rank; + } break; + 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; + default: + assert(false); + } + assert(std::isfinite(sum_weights)); + } else { // Threshold for this patch is invalid + // Erring on the conservative side: + // If the decomposition fails, or a threshold can't be found, + // copy the input data to the output data as-is, + // regardless of whether performing overcomplete local PCA + w.head(r).setOnes(); + out_rank = r; + sum_weights = r; + } + assert(w.head(r).allFinite()); + + // 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 +#ifdef DWIDENOISE2_USE_BDCSVD + const Eigen::Matrix wrev = w.head(r).reverse().cast(); +#endif + switch (aggregator) { + case aggregator_type::EXCLUSIVE: { + assert(Estimate::patch.centre_index >= 0); + if (bool(Estimate::threshold)) { +#ifdef DWIDENOISE2_USE_BDCSVD + assert(Estimate::SVD.matrixU().allFinite()); + assert(Estimate::SVD.matrixV().allFinite()); + assert(wrev.allFinite()); + assert(Estimate::SVD.singularValues().allFinite()); + // TODO Re-try reconstruction without use of V: + // https://github.com/MRtrix3/mrtrix3/pull/2906/commits/eb34f3c57dd460d2b3bd86b9653066be15e916c6 + // It might be that in the case of anything other than EXCLUSIVE, + // computing V is no more expensive than doing the full patch reconstruction in its absence, + // whereas for EXCLUSIVE since only a small portion of V is used it's worthwhile + Xr.noalias() = Estimate::SVD.matrixU() * // + (wrev.array() * Estimate::SVD.singularValues().array()).matrix().asDiagonal() * // + Estimate::SVD.matrixV().row(Estimate::patch.centre_index).adjoint(); // +#else + if (Estimate::m <= n) + Xr.noalias() = // + Estimate::eig.eigenvectors() * // + (w.head(r).cast().matrix().asDiagonal() * // + (Estimate::eig.eigenvectors().adjoint() * // + Estimate::X.col(Estimate::patch.centre_index))); // + else + Xr.noalias() = // + Estimate::X.leftCols(n) * // + (Estimate::eig.eigenvectors() * // + (w.head(r).cast().matrix().asDiagonal() * // + Estimate::eig.eigenvectors().adjoint().col(Estimate::patch.centre_index))); // +#endif + assert(Xr.allFinite()); + } else { + // In the case of -aggregator exclusive, + // where a decomposition fails or we can't find a threshold, + // we simply copy the input data into the output image + Xr.noalias() = Estimate::X.col(Estimate::patch.centre_index); + } + assign_pos_of(dwi).to(out); + out.row(3) = Xr.col(0); + 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: { // All aggregators other than EXCLUSIVE + if (!Estimate::threshold) { + Xr.leftCols(n).noalias() = Estimate::X.leftCols(n); +#ifdef DWIDENOISE2_USE_BDCSVD + } else { + Xr.leftCols(n).noalias() = Estimate::SVD.matrixU() * // + (wrev.array() * Estimate::SVD.singularValues().array()).matrix().asDiagonal() * // + Estimate::SVD.matrixV().adjoint(); // +#else + } else if (Estimate::m <= n) { + Xr.leftCols(n).noalias() = // + Estimate::eig.eigenvectors() * // + (w.head(r).cast().matrix().asDiagonal() * // + (Estimate::eig.eigenvectors().adjoint() * // + Estimate::X.leftCols(n))); // + } else { + Xr.leftCols(n).noalias() = // + Estimate::X.leftCols(n) * // + (Estimate::eig.eigenvectors() * // + (w.head(r).cast().matrix().asDiagonal() * // + Estimate::eig.eigenvectors().adjoint())); // +#endif + } + assert(Xr.leftCols(n).allFinite()); + // Undo prior within-patch variance-stabilising transform + if (std::isfinite(Estimate::patch.centre_noise)) { + for (ssize_t i = 0; i != n; ++i) + if (Estimate::patch.voxels[i].noise_level > 0.0) + Xr.col(i) *= Estimate::patch.voxels[i].noise_level / Estimate::patch.centre_noise; + } + std::lock_guard lock(Estimate::mutex); + for (size_t voxel_index = 0; voxel_index != Estimate::patch.voxels.size(); ++voxel_index) { + assign_pos_of(Estimate::patch.voxels[voxel_index].index, 0, 3).to(out); + assign_pos_of(Estimate::patch.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::patch.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 * Xr.col(voxel_index); + Estimate::exports.sum_aggregation.value() += weight; + if (Estimate::exports.rank_output.valid()) { + assign_pos_of(Estimate::patch.voxels[voxel_index].index, 0, 3).to(Estimate::exports.rank_output); + Estimate::exports.rank_output.value() += weight * out_rank; + } + } + } break; + } + + auto ss_index = Estimate::subsample->in2ss({dwi.index(0), dwi.index(1), dwi.index(2)}); + if (Estimate::exports.sum_optshrink.valid()) { + assign_pos_of(ss_index, 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..6db9e40d59 --- /dev/null +++ b/src/denoise/recon.h @@ -0,0 +1,64 @@ +/* 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, + std::shared_ptr subsample, + 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; + typename Estimate::MatrixType Xr; + std::map beta2lambdastar; +}; + +template class Recon; +template class Recon; +template class Recon; +template class Recon; + +} // namespace MR::Denoise diff --git a/src/denoise/subsample.cpp b/src/denoise/subsample.cpp new file mode 100644 index 0000000000..ae9e91c824 --- /dev/null +++ b/src/denoise/subsample.cpp @@ -0,0 +1,108 @@ +/* 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/subsample.h" + +namespace MR::Denoise { + +const App::Option subsample_option = + App::Option("subsample", + "reduce the number of PCA kernels relative to the number of image voxels; " + "can provide either an integer subsampling factor, " + "or a comma-separated list of three factors; " + "default: 2") + + App::Argument("factor").type_integer(1); + +Subsample::Subsample(const Header &in, const std::array &factors) + : H_in(make_input_header(in)), + factors(factors), + size({(in.size(0) + factors[0] - 1) / factors[0], + (in.size(1) + factors[1] - 1) / factors[1], + (in.size(2) + factors[2] - 1) / factors[2]}), + origin({(in.size(0) - factors[0] * (size[0] - 1) - 1) / 2, + (in.size(1) - factors[1] * (size[1] - 1) - 1) / 2, + (in.size(2) - factors[2] * (size[2] - 1) - 1) / 2}), + H_ss(make_subsample_header()) {} + +bool Subsample::process(const Kernel::Voxel::index_type &pos) const { + for (ssize_t axis = 0; axis != 3; ++axis) { + if (pos[axis] % factors[axis] != origin[axis]) + return false; + } + return true; +} + +std::array Subsample::in2ss(const Kernel::Voxel::index_type &pos) const { + // Do not attempt to map an unprocessed voxel to a voxel index in subsampled space + assert(process(pos)); + assert(!is_out_of_bounds(H_in, pos, 0, 3)); + return std::array({(pos[0] - origin[0]) / factors[0], // + (pos[1] - origin[1]) / factors[1], // + (pos[2] - origin[2]) / factors[2]}); // +} + +std::array Subsample::ss2in(const Kernel::Voxel::index_type &pos) const { + assert(!is_out_of_bounds(H_ss, pos)); + return std::array({pos[0] * factors[0] + origin[0], // + pos[1] * factors[1] + origin[1], // + pos[2] * factors[2] + origin[2]}); // +} + +std::shared_ptr Subsample::make(const Header &in, const ssize_t default_ratio) { + auto opt = App::get_options("subsample"); + std::array factors({default_ratio, default_ratio, default_ratio}); + if (!opt.empty()) { + const std::vector userinput = parse_ints(opt[0][0]); + if (userinput.size() == 1) + factors = {userinput[0], userinput[0], userinput[0]}; + else if (userinput.size() == 3) + factors = {userinput[0], userinput[1], userinput[2]}; + else + throw Exception("Subsampling factor must be either a single positive integer, " + "or a comma-separated list of three positive integers"); + } + return std::make_shared(in, factors); +} + +Header Subsample::make_input_header(const Header &H_in) const { + Header H(H_in); + H.ndim() = 3; + H.reset_intensity_scaling(); + H.datatype() = DataType::Float32; + H.datatype().set_byte_order_native(); + return H; +} + +Header Subsample::make_subsample_header() const { + Header H(H_in); + H.ndim() = 3; + H.reset_intensity_scaling(); + H.datatype() = DataType::Float32; + H.datatype().set_byte_order_native(); + std::array halfvoxel_offsets; + for (ssize_t axis = 0; axis != 3; ++axis) { + H.size(axis) = size[axis]; + H.spacing(axis) *= factors[axis]; + halfvoxel_offsets[axis] = factors[axis] & 1 ? 0.0 : 0.5; + } + H.transform().translation() = + H_in.transform() * Eigen::Matrix({(origin[0] + halfvoxel_offsets[0]) * H_in.spacing(0), // + (origin[1] + halfvoxel_offsets[1]) * H_in.spacing(1), // + (origin[2] + halfvoxel_offsets[2]) * H_in.spacing(2)}); // + return H; +} + +} // namespace MR::Denoise diff --git a/src/denoise/subsample.h b/src/denoise/subsample.h new file mode 100644 index 0000000000..7333ed0491 --- /dev/null +++ b/src/denoise/subsample.h @@ -0,0 +1,56 @@ +/* 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 "app.h" +#include "denoise/denoise.h" +#include "denoise/kernel/voxel.h" +#include "header.h" + +namespace MR::Denoise { + +extern const App::Option subsample_option; + +class Subsample { +public: + Subsample(const Header &in, const std::array &factors); + + const Header &header() const { return H_ss; } + + // TODO May want to move definition of Kernel::Voxel out of Kernel namespace + bool process(const Kernel::Voxel::index_type &pos) const; + std::array in2ss(const Kernel::Voxel::index_type &pos) const; + std::array ss2in(const Kernel::Voxel::index_type &pos) const; + const std::array &get_factors() const { return factors; } + + static std::shared_ptr make(const Header &in, const ssize_t default_ratio); + +protected: + const Header H_in; + const std::array factors; + const std::array size; + const std::array origin; + const Header H_ss; + + Header make_input_header(const Header &) const; + Header make_subsample_header() const; +}; + +} // namespace MR::Denoise