diff --git a/Modules/Filtering/AdaptiveDenoising/CMakeLists.txt b/Modules/Filtering/AdaptiveDenoising/CMakeLists.txt new file mode 100644 index 00000000000..e0387995bf4 --- /dev/null +++ b/Modules/Filtering/AdaptiveDenoising/CMakeLists.txt @@ -0,0 +1,5 @@ +project(AdaptiveDenoising) + +set(AdaptiveDenoising_LIBRARIES AdaptiveDenoising) + +itk_module_impl() diff --git a/Modules/Filtering/AdaptiveDenoising/LICENSE b/Modules/Filtering/AdaptiveDenoising/LICENSE new file mode 100644 index 00000000000..d6456956733 --- /dev/null +++ b/Modules/Filtering/AdaptiveDenoising/LICENSE @@ -0,0 +1,202 @@ + + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + + APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "[]" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + + Copyright [yyyy] [name of copyright owner] + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. diff --git a/Modules/Filtering/AdaptiveDenoising/README.md b/Modules/Filtering/AdaptiveDenoising/README.md new file mode 100644 index 00000000000..e93d97bee27 --- /dev/null +++ b/Modules/Filtering/AdaptiveDenoising/README.md @@ -0,0 +1,31 @@ +# AdaptiveDenoising + +In-tree ITK module providing patch-based adaptive denoising filters for +scalar and vector images. The flagship class is +`itk::AdaptiveNonLocalMeansDenoisingImageFilter`, which implements a +non-local-means denoiser whose patch-similarity radius and smoothing +strength adapt to local image content. Supporting infrastructure +includes `itk::NonLocalPatchBasedImageFilter` (the base class for +non-local-means style filters) and `itk::VarianceImageFilter` (a +neighborhood variance estimator used to drive the adaptive parameter +selection). + +## Origin + +Ingested from the standalone remote module +[**ntustison/ITKAdaptiveDenoising**](https://github.com/ntustison/ITKAdaptiveDenoising) +on 2026-05-08, following the v4 ingestion guidelines defined in +InsightSoftwareConsortium/ITK#6204. The upstream repository will be +archived read-only after this PR merges; it remains reachable at the +URL above for historical reference (notably the `examples/` +directory, which was intentionally left in the upstream archive). + +## References + +- Manjón J.V., Coupé P., Martí-Bonmatí L., Collins D.L., Robles M. + *Adaptive non-local means denoising of MR images with spatially + varying noise levels.* Journal of Magnetic Resonance Imaging. + 2010;31(1):192-203. +- Tustison N., Manjón J.V. *Adaptive Non-Local Means Filtering for + Image Denoising.* The Insight Journal. 2020. + diff --git a/Modules/Filtering/AdaptiveDenoising/include/itkAdaptiveNonLocalMeansDenoisingImageFilter.h b/Modules/Filtering/AdaptiveDenoising/include/itkAdaptiveNonLocalMeansDenoisingImageFilter.h new file mode 100644 index 00000000000..55b826aa038 --- /dev/null +++ b/Modules/Filtering/AdaptiveDenoising/include/itkAdaptiveNonLocalMeansDenoisingImageFilter.h @@ -0,0 +1,222 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkAdaptiveNonLocalMeansDenoisingImageFilter_h +#define itkAdaptiveNonLocalMeansDenoisingImageFilter_h + +#include "itkNonLocalPatchBasedImageFilter.h" + +#include "itkConstNeighborhoodIterator.h" +#include "itkGaussianOperator.h" + +namespace itk +{ + +/** + * \class AdaptiveNonLocalMeansDenoisingImageFilter + * \brief Implementation of a denoising image filter. + * + * \author Jose V. Manjon with ITK porting by Nick Tustison + * + * Contributed by + * + * \par REFERENCE + * + * J. V. Manjon, P. Coupe, Luis Marti-Bonmati, D. L. Collins, + * and M. Robles. "Adaptive Non-Local Means Denoising of MR Images With + * Spatially Varying Noise Levels, Journal of Magnetic Resonance Imaging, + * 31:192-203, June 2010. + * + * \ingroup AdaptiveDenoising + */ + +template > +class ITK_TEMPLATE_EXPORT AdaptiveNonLocalMeansDenoisingImageFilter final + : public NonLocalPatchBasedImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(AdaptiveNonLocalMeansDenoisingImageFilter); + + /** Standard class typedefs. */ + using Self = AdaptiveNonLocalMeansDenoisingImageFilter; + using Superclass = NonLocalPatchBasedImageFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Runtime information support. */ + itkOverrideGetNameOfClassMacro(AdaptiveNonLocalMeansDenoisingImageFilter); + + /** Standard New method. */ + itkNewMacro(Self); + + /** ImageDimension constants */ + static constexpr unsigned int ImageDimension = TInputImage::ImageDimension; + + /** Some convenient typedefs. */ + using InputImageType = TInputImage; + using InputPixelType = typename InputImageType::PixelType; + using OutputImageType = TOutputImage; + using RegionType = typename Superclass::RegionType; + + using MaskImageType = TMaskImage; + using MaskPixelType = typename MaskImageType::PixelType; + using LabelType = typename MaskImageType::PixelType; + + using RealType = typename Superclass::RealType; + using RealImageType = typename Superclass::RealImageType; + using RealImagePointer = typename Superclass::RealImagePointer; + using IndexType = typename Superclass::IndexType; + + using ConstNeighborhoodIteratorType = typename Superclass::ConstNeighborhoodIteratorType; + using NeighborhoodRadiusType = typename Superclass::NeighborhoodRadiusType; + using NeighborhoodOffsetType = typename Superclass::NeighborhoodOffsetType; + using NeighborhoodOffsetListType = typename Superclass::NeighborhoodOffsetListType; + + using ModifiedBesselCalculatorType = GaussianOperator; + + /** + * The image expected for input for noise correction. + */ + void + SetInput1(const InputImageType * image) + { + this->SetInput(image); + } + + /** + * Set mask image function. If a binary mask image is specified, only + * those input image voxels corresponding with the mask image. + */ + void + SetMaskImage(const MaskImageType * mask) + { + this->SetNthInput(1, const_cast(mask)); + } + void + SetInput2(const MaskImageType * mask) + { + this->SetMaskImage(mask); + } + + /** + * Get mask image function. If a binary mask image is specified, only + * those input image voxels corresponding with the mask image. + */ + const MaskImageType * + GetMaskImage() const + { + return static_cast(this->ProcessObject::GetInput(1)); + } + + /** + * Employ Rician noise model. Otherwise use a Gaussian noise model. + * Default = true. + */ + itkSetMacro(UseRicianNoiseModel, bool); + itkGetConstMacro(UseRicianNoiseModel, bool); + itkBooleanMacro(UseRicianNoiseModel); + + /** + * Smoothing factor for noise. Default = 1.0. + */ + itkSetMacro(SmoothingFactor, RealType); + itkGetConstMacro(SmoothingFactor, RealType); + + /** + * Smoothing variance for Rician noise. Default = 2.0. + */ + itkSetMacro(SmoothingVariance, RealType); + itkGetConstMacro(SmoothingVariance, RealType); + + /** + * Epsilon for minimum value of mean and variance at a pixel. + * Default = 0.00001. + */ + itkSetMacro(Epsilon, RealType); + itkGetConstMacro(Epsilon, RealType); + + /** + * Mean threshold. + * Default = 0.95. + */ + itkSetMacro(MeanThreshold, RealType); + itkGetConstMacro(MeanThreshold, RealType); + + /** + * Variance threshold. + * Default = 0.5. + */ + itkSetMacro(VarianceThreshold, RealType); + itkGetConstMacro(VarianceThreshold, RealType); + + /** + * Neighborhood for computing local mean and variance images. + * Default = 1x1x... + */ + itkSetMacro(NeighborhoodRadiusForLocalMeanAndVariance, NeighborhoodRadiusType); + itkGetConstMacro(NeighborhoodRadiusForLocalMeanAndVariance, NeighborhoodRadiusType); + +protected: + AdaptiveNonLocalMeansDenoisingImageFilter(); + ~AdaptiveNonLocalMeansDenoisingImageFilter() override = default; + + void + PrintSelf(std::ostream & os, Indent indent) const override; + + void + ThreadedGenerateData(const RegionType &, ThreadIdType) override; + + void + BeforeThreadedGenerateData() override; + + void + AfterThreadedGenerateData() override; + +private: + RealType CalculateCorrectionFactor(RealType); + + bool m_UseRicianNoiseModel; + + ModifiedBesselCalculatorType m_ModifiedBesselCalculator; + + RealType m_Epsilon; + RealType m_MeanThreshold; + RealType m_VarianceThreshold; + RealType m_SmoothingFactor; + RealType m_SmoothingVariance; + + RealType m_MaximumInputPixelIntensity; + RealType m_MinimumInputPixelIntensity; + + RealImagePointer m_MeanImage; + RealImagePointer m_RicianBiasImage; + RealImagePointer m_VarianceImage; + RealImagePointer m_ThreadContributionCountImage; + RealImagePointer m_IntensitySquaredDistanceImage; + + NeighborhoodRadiusType m_NeighborhoodRadiusForLocalMeanAndVariance; +}; + +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkAdaptiveNonLocalMeansDenoisingImageFilter.hxx" +#endif + +#endif diff --git a/Modules/Filtering/AdaptiveDenoising/include/itkAdaptiveNonLocalMeansDenoisingImageFilter.hxx b/Modules/Filtering/AdaptiveDenoising/include/itkAdaptiveNonLocalMeansDenoisingImageFilter.hxx new file mode 100644 index 00000000000..f0f76c58073 --- /dev/null +++ b/Modules/Filtering/AdaptiveDenoising/include/itkAdaptiveNonLocalMeansDenoisingImageFilter.hxx @@ -0,0 +1,524 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkAdaptiveNonLocalMeansDenoisingImageFilter_hxx +#define itkAdaptiveNonLocalMeansDenoisingImageFilter_hxx + + +#include "itkArray.h" +#include "itkDiscreteGaussianImageFilter.h" +#include "itkImageRegionConstIterator.h" +#include "itkImageRegionIterator.h" +#include "itkImageRegionIteratorWithIndex.h" +#include "itkMath.h" +#include "itkMeanImageFilter.h" +#include "itkNeighborhoodIterator.h" +#include "itkProgressReporter.h" +#include "itkStatisticsImageFilter.h" +#include "itkVarianceImageFilter.h" + +#include + +namespace itk +{ + +template +AdaptiveNonLocalMeansDenoisingImageFilter:: + AdaptiveNonLocalMeansDenoisingImageFilter() + : m_UseRicianNoiseModel(true) + , m_Epsilon(0.00001) + , m_MeanThreshold(0.95) + , m_VarianceThreshold(0.5) + , m_SmoothingFactor(1.0) + , m_SmoothingVariance(2.0) + , m_MaximumInputPixelIntensity(NumericTraits::NonpositiveMin()) + , m_MinimumInputPixelIntensity(NumericTraits::max()) +{ + this->SetNumberOfRequiredInputs(1); + + this->m_MeanImage = nullptr; + this->m_VarianceImage = nullptr; + this->m_IntensitySquaredDistanceImage = nullptr; + this->m_ThreadContributionCountImage = nullptr; + + this->m_RicianBiasImage = nullptr; + + this->m_NeighborhoodRadiusForLocalMeanAndVariance.Fill(1); + this->DynamicMultiThreadingOff(); +} + +template +void +AdaptiveNonLocalMeansDenoisingImageFilter::BeforeThreadedGenerateData() +{ + Superclass::BeforeThreadedGenerateData(); + + const InputImageType * inputImage = this->GetInput(); + + typedef MeanImageFilter MeanImageFilterType; + typename MeanImageFilterType::Pointer meanImageFilter = MeanImageFilterType::New(); + meanImageFilter->SetInput(inputImage); + meanImageFilter->SetRadius(this->GetNeighborhoodRadiusForLocalMeanAndVariance()); + + this->m_MeanImage = meanImageFilter->GetOutput(); + this->m_MeanImage->Update(); + this->m_MeanImage->DisconnectPipeline(); + + typedef VarianceImageFilter VarianceImageFilterType; + typename VarianceImageFilterType::Pointer varianceImageFilter = VarianceImageFilterType::New(); + varianceImageFilter->SetInput(inputImage); + varianceImageFilter->SetRadius(this->GetNeighborhoodRadiusForLocalMeanAndVariance()); + + this->m_VarianceImage = varianceImageFilter->GetOutput(); + this->m_VarianceImage->Update(); + this->m_VarianceImage->DisconnectPipeline(); + + typedef StatisticsImageFilter StatsFilterType; + typename StatsFilterType::Pointer statsFilter = StatsFilterType::New(); + statsFilter->SetInput(inputImage); + statsFilter->Update(); + + this->m_MaximumInputPixelIntensity = static_cast(statsFilter->GetMaximum()); + this->m_MinimumInputPixelIntensity = static_cast(statsFilter->GetMinimum()); + + this->m_ThreadContributionCountImage = RealImageType::New(); + this->m_ThreadContributionCountImage->CopyInformation(inputImage); + this->m_ThreadContributionCountImage->SetRegions(inputImage->GetRequestedRegion()); + this->m_ThreadContributionCountImage->Allocate(true); + + if (this->m_UseRicianNoiseModel) + { + this->m_RicianBiasImage = RealImageType::New(); + this->m_RicianBiasImage->CopyInformation(inputImage); + this->m_RicianBiasImage->SetRegions(inputImage->GetRequestedRegion()); + this->m_RicianBiasImage->Allocate(true); + } + + this->AllocateOutputs(); + // Output buffer needs to be zero initialized + this->GetOutput()->FillBuffer(0.0); +} + +template +void +AdaptiveNonLocalMeansDenoisingImageFilter::ThreadedGenerateData( + const RegionType & region, + ThreadIdType threadId) +{ + ProgressReporter progress(this, threadId, region.GetNumberOfPixels(), 100); + + const InputImageType * inputImage = this->GetInput(); + const MaskImageType * maskImage = this->GetMaskImage(); + + OutputImageType * outputImage = this->GetOutput(); + RegionType targetImageRegion = this->GetTargetImageRegion(); + + NeighborhoodOffsetListType neighborhoodPatchOffsetList = this->GetNeighborhoodPatchOffsetList(); + + NeighborhoodRadiusType neighborhoodSearchRadius = this->GetNeighborhoodSearchRadius(); + + ConstNeighborhoodIterator ItV(neighborhoodSearchRadius, this->m_VarianceImage, region); + ConstNeighborhoodIterator ItM(neighborhoodSearchRadius, this->m_MeanImage, region); + + const unsigned int neighborhoodSearchSize = this->GetNeighborhoodSearchSize(); + const unsigned int neighborhoodPatchSize = this->GetNeighborhoodPatchSize(); + + Array weightedAverageIntensities(neighborhoodPatchSize); + + ItM.GoToBegin(); + ItV.GoToBegin(); + + while (!ItM.IsAtEnd()) + { + typename InputImageType::IndexType centerIndex = ItM.GetIndex(); + + InputPixelType inputCenterPixel = inputImage->GetPixel(centerIndex); + RealType meanCenterPixel = this->m_MeanImage->GetPixel(centerIndex); + RealType varianceCenterPixel = this->m_VarianceImage->GetPixel(centerIndex); + + RealType maxWeight = NumericTraits::ZeroValue(); + RealType sumOfWeights = NumericTraits::ZeroValue(); + + weightedAverageIntensities.Fill(NumericTraits::ZeroValue()); + + RealType meanNeighborhoodPixel = NumericTraits::ZeroValue(); + RealType varianceNeighborhoodPixel = NumericTraits::ZeroValue(); + + if (inputCenterPixel > 0 && meanCenterPixel > this->m_Epsilon && varianceCenterPixel > this->m_Epsilon && + (!maskImage || maskImage->GetPixel(centerIndex) != NumericTraits::ZeroValue())) + { + // Calculate the minimum distance + + RealType minimumDistance = NumericTraits::max(); + for (unsigned int m = 0; m < neighborhoodSearchSize; m++) + { + if (!ItM.IndexInBounds(m) || m == static_cast(0.5 * neighborhoodSearchSize)) + { + continue; + } + + IndexType neighborhoodIndex = ItM.GetIndex(m); + + if (inputImage->GetPixel(neighborhoodIndex) <= 0) + { + continue; + } + + meanNeighborhoodPixel = this->m_MeanImage->GetPixel(neighborhoodIndex); + varianceNeighborhoodPixel = this->m_VarianceImage->GetPixel(neighborhoodIndex); + + if (meanNeighborhoodPixel <= this->m_Epsilon || varianceNeighborhoodPixel <= this->m_Epsilon) + { + continue; + } + + const RealType meanRatio = meanCenterPixel / meanNeighborhoodPixel; + const RealType meanRatioInverse = (this->m_MaximumInputPixelIntensity - meanCenterPixel) / + (this->m_MaximumInputPixelIntensity - meanNeighborhoodPixel); + + const RealType varianceRatio = varianceCenterPixel / varianceNeighborhoodPixel; + + if (((meanRatio > this->m_MeanThreshold && + meanRatio < itk::NumericTraits::OneValue() / this->m_MeanThreshold) || + (meanRatioInverse > this->m_MeanThreshold && + meanRatioInverse < itk::NumericTraits::OneValue() / this->m_MeanThreshold)) && + varianceRatio > this->m_VarianceThreshold && + varianceRatio < itk::NumericTraits::OneValue() / this->m_VarianceThreshold) + { + + RealType averageDistance = itk::NumericTraits::ZeroValue(); + RealType count = itk::NumericTraits::ZeroValue(); + + for (unsigned int n = 0; n < neighborhoodPatchSize; n++) + { + IndexType neighborhoodPatchIndex = neighborhoodIndex + neighborhoodPatchOffsetList[n]; + + if (!targetImageRegion.IsInside(neighborhoodPatchIndex)) + { + continue; + } + RealType neighborhoodInputImagePixel = static_cast(inputImage->GetPixel(neighborhoodPatchIndex)); + RealType neighborhoodMeanImagePixel = this->m_MeanImage->GetPixel(neighborhoodPatchIndex); + averageDistance += itk::Math::sqr(neighborhoodInputImagePixel - neighborhoodMeanImagePixel); + + count += itk::NumericTraits::OneValue(); + } + averageDistance /= count; + minimumDistance = std::min(averageDistance, minimumDistance); + } + } + + if (itk::Math::AlmostEquals(minimumDistance, NumericTraits::ZeroValue())) + { + minimumDistance = NumericTraits::OneValue(); + } + + // Rician correction + + if (this->m_UseRicianNoiseModel) + { + for (unsigned int n = 0; n < neighborhoodPatchSize; n++) + { + IndexType neighborhoodPatchIndex = centerIndex + neighborhoodPatchOffsetList[n]; + if (!targetImageRegion.IsInside(neighborhoodPatchIndex)) + { + continue; + } + + if (itk::Math::AlmostEquals(minimumDistance, NumericTraits::max())) + { + this->m_RicianBiasImage->SetPixel(neighborhoodPatchIndex, 0.0); + } + else + { + this->m_RicianBiasImage->SetPixel(neighborhoodPatchIndex, minimumDistance); + } + } + } + + // Patch filtering + + for (unsigned int m = 0; m < neighborhoodSearchSize; m++) + { + if (!ItM.IndexInBounds(m) || m == static_cast(0.5 * neighborhoodSearchSize)) + { + continue; + } + + IndexType neighborhoodIndex = ItM.GetIndex(m); + + if (inputImage->GetPixel(neighborhoodIndex) <= 0) + { + continue; + } + + meanNeighborhoodPixel = this->m_MeanImage->GetPixel(neighborhoodIndex); + varianceNeighborhoodPixel = this->m_VarianceImage->GetPixel(neighborhoodIndex); + + if (meanNeighborhoodPixel <= this->m_Epsilon || varianceNeighborhoodPixel <= this->m_Epsilon) + { + continue; + } + + const RealType meanRatio = meanCenterPixel / meanNeighborhoodPixel; + const RealType meanRatioInverse = (this->m_MaximumInputPixelIntensity - meanCenterPixel) / + (this->m_MaximumInputPixelIntensity - meanNeighborhoodPixel); + + const RealType varianceRatio = varianceCenterPixel / varianceNeighborhoodPixel; + + if (((meanRatio > this->m_MeanThreshold && + meanRatio < itk::NumericTraits::OneValue() / this->m_MeanThreshold) || + (meanRatioInverse > this->m_MeanThreshold && + meanRatioInverse < itk::NumericTraits::OneValue() / this->m_MeanThreshold)) && + varianceRatio > this->m_VarianceThreshold && + varianceRatio < itk::NumericTraits::OneValue() / this->m_VarianceThreshold) + { + + RealType averageDistance = 0.0; + RealType count = 0.0; + for (unsigned int n = 0; n < neighborhoodPatchSize; n++) + { + IndexType searchNeighborhoodPatchIndex = neighborhoodIndex + neighborhoodPatchOffsetList[n]; + IndexType centerNeighborhoodPatchIndex = centerIndex + neighborhoodPatchOffsetList[n]; + if (!targetImageRegion.IsInside(searchNeighborhoodPatchIndex) || + !targetImageRegion.IsInside(centerNeighborhoodPatchIndex)) + { + continue; + } + RealType distance1 = inputImage->GetPixel(searchNeighborhoodPatchIndex) - + this->m_MeanImage->GetPixel(searchNeighborhoodPatchIndex); + RealType distance2 = inputImage->GetPixel(centerNeighborhoodPatchIndex) - + this->m_MeanImage->GetPixel(centerNeighborhoodPatchIndex); + averageDistance += itk::Math::sqr(distance1 - distance2); + count += itk::NumericTraits::OneValue(); + } + averageDistance /= count; + + RealType weight = itk::NumericTraits::ZeroValue(); + if (averageDistance <= static_cast(3.0) * minimumDistance) + { + weight = std::exp(-averageDistance / minimumDistance); + } + if (weight > maxWeight) + { + maxWeight = weight; + } + + if (weight > itk::NumericTraits::ZeroValue()) + { + for (unsigned int n = 0; n < neighborhoodPatchSize; n++) + { + IndexType neighborhoodPatchIndex = neighborhoodIndex + neighborhoodPatchOffsetList[n]; + if (!targetImageRegion.IsInside(neighborhoodPatchIndex)) + { + continue; + } + if (this->m_UseRicianNoiseModel) + { + weightedAverageIntensities[n] += weight * itk::Math::sqr(inputImage->GetPixel(neighborhoodPatchIndex)); + } + else + { + weightedAverageIntensities[n] += weight * inputImage->GetPixel(neighborhoodPatchIndex); + } + } + sumOfWeights += weight; + } + } + } + + if (itk::Math::AlmostEquals(maxWeight, NumericTraits::ZeroValue())) + { + maxWeight = NumericTraits::OneValue(); + } + } + else + { + maxWeight = NumericTraits::OneValue(); + } + + for (unsigned int n = 0; n < neighborhoodPatchSize; n++) + { + IndexType neighborhoodPatchIndex = centerIndex + neighborhoodPatchOffsetList[n]; + if (!targetImageRegion.IsInside(neighborhoodPatchIndex)) + { + continue; + } + if (this->m_UseRicianNoiseModel) + { + weightedAverageIntensities[n] += maxWeight * itk::Math::sqr(inputImage->GetPixel(neighborhoodPatchIndex)); + } + else + { + weightedAverageIntensities[n] += maxWeight * inputImage->GetPixel(neighborhoodPatchIndex); + } + } + sumOfWeights += maxWeight; + + if (sumOfWeights > itk::NumericTraits::ZeroValue()) + { + for (unsigned int n = 0; n < neighborhoodPatchSize; n++) + { + IndexType neighborhoodPatchIndex = centerIndex + neighborhoodPatchOffsetList[n]; + if (!targetImageRegion.IsInside(neighborhoodPatchIndex)) + { + continue; + } + typename OutputImageType::PixelType estimate = outputImage->GetPixel(neighborhoodPatchIndex); + estimate += (weightedAverageIntensities[n] / sumOfWeights); + + outputImage->SetPixel(neighborhoodPatchIndex, estimate); + this->m_ThreadContributionCountImage->SetPixel( + neighborhoodPatchIndex, this->m_ThreadContributionCountImage->GetPixel(neighborhoodPatchIndex) + 1); + } + } + + ++ItM; + ++ItV; + + progress.CompletedPixel(); + } +} + +template +void +AdaptiveNonLocalMeansDenoisingImageFilter::AfterThreadedGenerateData() +{ + const MaskImageType * maskImage = this->GetMaskImage(); + + if (this->m_UseRicianNoiseModel) + { + typedef DiscreteGaussianImageFilter SmootherType; + typename SmootherType::Pointer smoother = SmootherType::New(); + smoother->SetInput(this->m_RicianBiasImage); + smoother->SetVariance(this->m_SmoothingVariance); + smoother->SetUseImageSpacing(true); + smoother->Update(); + + ImageRegionConstIterator ItS(smoother->GetOutput(), smoother->GetOutput()->GetRequestedRegion()); + ImageRegionConstIteratorWithIndex ItM(this->m_MeanImage, this->m_MeanImage->GetRequestedRegion()); + ImageRegionIterator ItB(this->m_RicianBiasImage, this->m_RicianBiasImage->GetRequestedRegion()); + ItS.GoToBegin(); + ItM.GoToBegin(); + ItB.GoToBegin(); + + while (!ItS.IsAtEnd()) + { + if (ItS.Get() > itk::NumericTraits::ZeroValue() && + (!maskImage || maskImage->GetPixel(ItM.GetIndex()) != NumericTraits::ZeroValue())) + { + const RealType snr = ItM.Get() / std::sqrt(ItS.Get()); + + RealType bias = static_cast(2.0) * ItS.Get() / this->CalculateCorrectionFactor(snr); + + if (std::isnan(bias) || std::isinf(bias)) + { + bias = itk::NumericTraits::ZeroValue(); + } + ItB.Set(bias); + } + + ++ItS; + ++ItM; + ++ItB; + } + } + + ImageRegionIteratorWithIndex ItO(this->GetOutput(), this->GetOutput()->GetRequestedRegion()); + ImageRegionConstIterator ItL(this->m_ThreadContributionCountImage, + this->m_ThreadContributionCountImage->GetRequestedRegion()); + + for (ItO.GoToBegin(), ItL.GoToBegin(); !ItO.IsAtEnd(); ++ItO, ++ItL) + { + RealType estimate = ItO.Get(); + + if (itk::Math::FloatAlmostEqual(ItL.Get(), itk::NumericTraits::ZeroValue())) + { + continue; + } + + estimate /= ItL.Get(); + + if (this->m_UseRicianNoiseModel) + { + RealType bias = this->m_RicianBiasImage->GetPixel(ItO.GetIndex()); + + estimate -= bias; + if (estimate < itk::NumericTraits::ZeroValue()) + { + estimate = itk::NumericTraits::ZeroValue(); + } + estimate = std::sqrt(estimate); + } + + ItO.Set(estimate); + } +} + +template +typename AdaptiveNonLocalMeansDenoisingImageFilter::RealType +AdaptiveNonLocalMeansDenoisingImageFilter::CalculateCorrectionFactor( + RealType snr) +{ + const RealType snrSquared = itk::Math::sqr(snr); + + RealType value = + static_cast(2.0) + snrSquared - + static_cast(0.125) * static_cast(Math::pi) * + static_cast(std::exp(static_cast(-0.5) * snrSquared)) * + itk::Math::sqr((static_cast(2.0) + snrSquared) * + static_cast( + this->m_ModifiedBesselCalculator.ModifiedBesselI0(static_cast(0.25) * snrSquared)) + + snrSquared * static_cast(this->m_ModifiedBesselCalculator.ModifiedBesselI1( + static_cast(0.25) * snrSquared))); + + if (value < static_cast(0.001) || value > static_cast(10.0)) + { + value = itk::NumericTraits::OneValue(); + } + return value; +} + +template +void +AdaptiveNonLocalMeansDenoisingImageFilter::PrintSelf(std::ostream & os, + Indent indent) const +{ + Superclass::PrintSelf(os, indent); + + if (this->m_UseRicianNoiseModel) + { + os << indent << "Using Rician noise model." << std::endl; + } + else + { + os << indent << "Using Gaussian noise model." << std::endl; + } + + os << indent << "Epsilon = " << this->m_Epsilon << std::endl; + os << indent << "Mean threshold = " << this->m_MeanThreshold << std::endl; + os << indent << "Variance threshold = " << this->m_VarianceThreshold << std::endl; + os << indent << "Smoothing variance = " << this->m_SmoothingVariance << std::endl; + + os << indent + << "Neighborhood radius for local mean and variance = " << this->m_NeighborhoodRadiusForLocalMeanAndVariance + << std::endl; +} + +} // end namespace itk + +#endif diff --git a/Modules/Filtering/AdaptiveDenoising/include/itkNonLocalPatchBasedImageFilter.h b/Modules/Filtering/AdaptiveDenoising/include/itkNonLocalPatchBasedImageFilter.h new file mode 100644 index 00000000000..90b00ebaec6 --- /dev/null +++ b/Modules/Filtering/AdaptiveDenoising/include/itkNonLocalPatchBasedImageFilter.h @@ -0,0 +1,217 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkNonLocalPatchBasedImageFilter_h +#define itkNonLocalPatchBasedImageFilter_h + +#include +#include "itkImageToImageFilter.h" +#include "AdaptiveDenoisingExport.h" + +#include "itkConstNeighborhoodIterator.h" + +namespace itk +{ + +/** + * \class NonLocalPatchBasedImageFilterEnums + * \brief Non-local patch-based image filter enum classes. + * \ingroup AdaptiveDenoising + */ + +class NonLocalPatchBasedImageFilterEnums +{ +public: + /**\class SimilarityMetric + * \brief Neighborhood patch similarity metric enumerated type. + * \ingroup AdaptiveDenoising + */ + enum class SimilarityMetric : uint8_t + { + PEARSON_CORRELATION = 0, + MEAN_SQUARES = 1 + }; +}; + +extern AdaptiveDenoising_EXPORT std::ostream & +operator<<(std::ostream & out, const NonLocalPatchBasedImageFilterEnums::SimilarityMetric value); + + +/** + * \class NonLocalPatchBasedImageFilter + * \brief Implementation of a non-local image filter. + * + * \ingroup AdaptiveDenoising + */ + +template +class ITK_TEMPLATE_EXPORT NonLocalPatchBasedImageFilter : public ImageToImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(NonLocalPatchBasedImageFilter); + + /** Standard class typedefs. */ + using Self = NonLocalPatchBasedImageFilter; + using Superclass = ImageToImageFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Runtime information support. */ + itkOverrideGetNameOfClassMacro(NonLocalPatchBasedImageFilter); + + /** Standard New method. */ + itkNewMacro(Self); + + /** ImageDimension constants */ + static constexpr unsigned int ImageDimension = TInputImage::ImageDimension; + + /** Some convenient typedefs. */ + using InputImageType = TInputImage; + using InputPixelType = typename InputImageType::PixelType; + using InputImagePointer = typename InputImageType::Pointer; + using InputImageList = std::vector; + using InputImageSetList = std::vector; + using RegionType = typename InputImageType::RegionType; + + using OutputImageType = TOutputImage; + using OutputPixelType = typename OutputImageType::PixelType; + + using InputImagePixelVectorType = std::vector; + + using RealType = float; + using RealImageType = Image; + using RealImagePointer = typename RealImageType::Pointer; + using IndexType = typename RealImageType::IndexType; + + using NeighborhoodType = Neighborhood; + using NeighborhoodSizeType = SizeValueType; + + using ConstNeighborhoodIteratorType = ConstNeighborhoodIterator; + using NeighborhoodRadiusType = typename ConstNeighborhoodIteratorType::RadiusType; + using NeighborhoodOffsetType = typename ConstNeighborhoodIteratorType::OffsetType; + + using NeighborhoodOffsetListType = std::vector; + + using SimilarityMetricEnum = NonLocalPatchBasedImageFilterEnums::SimilarityMetric; +#if !defined(ITK_LEGACY_REMOVE) + using SimilarityMetricType = SimilarityMetricEnum; + static constexpr SimilarityMetricType PEARSON_CORRELATION = SimilarityMetricEnum::PEARSON_CORRELATION; + static constexpr SimilarityMetricType MEAN_SQUARES = SimilarityMetricEnum::MEAN_SQUARES; +#endif + + /** + * Get/set neighborhood search radius. + * Default = 3x3x... + */ + itkSetMacro(NeighborhoodSearchRadius, NeighborhoodRadiusType); + itkGetConstMacro(NeighborhoodSearchRadius, NeighborhoodRadiusType); + + /** + * Get/set neighborhood search size. + */ + itkSetMacro(NeighborhoodSearchSize, NeighborhoodSizeType); + itkGetConstMacro(NeighborhoodSearchSize, NeighborhoodSizeType); + + /** + * Get/set neighborhood search offset list. + */ + virtual void + SetNeighborhoodSearchOffsetList(const NeighborhoodOffsetListType list) + { + this->m_NeighborhoodSearchOffsetList = list; + this->Modified(); + } + itkGetConstMacro(NeighborhoodSearchOffsetList, NeighborhoodOffsetListType); + + /** + * Get/set neighborhood patch radius. + * Default = 1x1x... + */ + itkSetMacro(NeighborhoodPatchRadius, NeighborhoodRadiusType); + itkGetConstMacro(NeighborhoodPatchRadius, NeighborhoodRadiusType); + + /** + * Get/set neighborhood patch size. + */ + itkSetMacro(NeighborhoodPatchSize, NeighborhoodSizeType); + itkGetConstMacro(NeighborhoodPatchSize, NeighborhoodSizeType); + + /** + * Get/set neighborhood patch offset list. + */ + virtual void + SetNeighborhoodPatchOffsetList(const NeighborhoodOffsetListType list) + { + this->m_NeighborhoodPatchOffsetList = list; + this->Modified(); + } + itkGetConstMacro(NeighborhoodPatchOffsetList, NeighborhoodOffsetListType); + + /** + * Enumerated type for neighborhood similarity. Default = MEAN_SQUARES + */ + itkSetMacro(SimilarityMetric, SimilarityMetricEnum); + itkGetConstMacro(SimilarityMetric, SimilarityMetricEnum); + +protected: + NonLocalPatchBasedImageFilter(); + ~NonLocalPatchBasedImageFilter() override = default; + + void + BeforeThreadedGenerateData() override; + + void + PrintSelf(std::ostream & os, Indent indent) const override; + + RealType + ComputeNeighborhoodPatchSimilarity(const InputImageList &, + const IndexType, + const InputImagePixelVectorType &, + const bool); + + InputImagePixelVectorType + VectorizeImageListPatch(const InputImageList &, const IndexType, const bool); + + InputImagePixelVectorType + VectorizeImagePatch(const InputImagePointer, const IndexType, const bool); + + void + GetMeanAndStandardDeviationOfVectorizedImagePatch(const InputImagePixelVectorType &, RealType &, RealType &); + + itkSetMacro(TargetImageRegion, RegionType); + itkGetConstMacro(TargetImageRegion, RegionType); + + SimilarityMetricEnum m_SimilarityMetric; + + SizeValueType m_NeighborhoodSearchSize; + NeighborhoodRadiusType m_NeighborhoodSearchRadius; + NeighborhoodOffsetListType m_NeighborhoodSearchOffsetList; + + SizeValueType m_NeighborhoodPatchSize; + NeighborhoodRadiusType m_NeighborhoodPatchRadius; + NeighborhoodOffsetListType m_NeighborhoodPatchOffsetList; + + RegionType m_TargetImageRegion; +}; + +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkNonLocalPatchBasedImageFilter.hxx" +#endif + +#endif diff --git a/Modules/Filtering/AdaptiveDenoising/include/itkNonLocalPatchBasedImageFilter.hxx b/Modules/Filtering/AdaptiveDenoising/include/itkNonLocalPatchBasedImageFilter.hxx new file mode 100644 index 00000000000..c1661890440 --- /dev/null +++ b/Modules/Filtering/AdaptiveDenoising/include/itkNonLocalPatchBasedImageFilter.hxx @@ -0,0 +1,260 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkNonLocalPatchBasedImageFilter_hxx +#define itkNonLocalPatchBasedImageFilter_hxx + + +#include "itkNeighborhood.h" + +namespace itk +{ + +template +NonLocalPatchBasedImageFilter::NonLocalPatchBasedImageFilter() +{ + this->m_SimilarityMetric = SimilarityMetricEnum::MEAN_SQUARES; + + this->m_NeighborhoodPatchRadius.Fill(1); + this->m_NeighborhoodPatchOffsetList.clear(); + + this->m_NeighborhoodSearchRadius.Fill(3); + this->m_NeighborhoodSearchOffsetList.clear(); +} + +template +void +NonLocalPatchBasedImageFilter::BeforeThreadedGenerateData() +{ + // Set up the search neighborhood parameters + + this->m_NeighborhoodSearchOffsetList.clear(); + + NeighborhoodType searchNeighborhood; + searchNeighborhood.SetRadius(this->m_NeighborhoodSearchRadius); + + this->m_NeighborhoodSearchSize = searchNeighborhood.Size(); + for (unsigned int n = 0; n < this->m_NeighborhoodSearchSize; n++) + { + this->m_NeighborhoodSearchOffsetList.push_back(searchNeighborhood.GetOffset(n)); + } + + // Set up the patch neighborhood parameters + + this->m_NeighborhoodPatchOffsetList.clear(); + + NeighborhoodType patchNeighborhood; + patchNeighborhood.SetRadius(this->m_NeighborhoodPatchRadius); + + this->m_NeighborhoodPatchSize = patchNeighborhood.Size(); + for (unsigned int n = 0; n < this->m_NeighborhoodPatchSize; n++) + { + this->m_NeighborhoodPatchOffsetList.push_back(patchNeighborhood.GetOffset(n)); + } + + this->m_TargetImageRegion = this->GetInput()->GetRequestedRegion(); +} + +template +typename NonLocalPatchBasedImageFilter::InputImagePixelVectorType +NonLocalPatchBasedImageFilter::VectorizeImageListPatch(const InputImageList & imageList, + const IndexType index, + const bool normalize) +{ + InputImagePixelVectorType patchVector(this->m_NeighborhoodPatchSize * imageList.size()); + for (unsigned int i = 0; i < imageList.size(); i++) + { + InputImagePixelVectorType patchVectorPerModality = this->VectorizeImagePatch(imageList[i], index, normalize); + for (unsigned int j = 0; j < this->m_NeighborhoodPatchSize; j++) + { + patchVector[i * this->m_NeighborhoodPatchSize + j] = patchVectorPerModality[j]; + } + } + return patchVector; +} + +template +typename NonLocalPatchBasedImageFilter::InputImagePixelVectorType +NonLocalPatchBasedImageFilter::VectorizeImagePatch(const InputImagePointer image, + const IndexType index, + const bool normalize) +{ + InputImagePixelVectorType patchVector(this->m_NeighborhoodPatchSize); + for (SizeValueType i = 0; i < this->m_NeighborhoodPatchSize; i++) + { + IndexType neighborhoodIndex = index + this->m_NeighborhoodPatchOffsetList[i]; + + bool isInBounds = this->m_TargetImageRegion.IsInside(neighborhoodIndex); + if (isInBounds) + { + InputPixelType pixel = image->GetPixel(neighborhoodIndex); + patchVector[i] = pixel; + } + else + { + patchVector[i] = std::numeric_limits::quiet_NaN(); + } + } + + if (normalize) + { + // Convert the current image patch to a z-score. + RealType mean = 0.0; + RealType standardDeviation = 0.0; + this->GetMeanAndStandardDeviationOfVectorizedImagePatch(patchVector, mean, standardDeviation); + if (standardDeviation < NumericTraits::epsilon()) + { + for (auto & it : patchVector) + { + it = 0.; + } + } + else + { + const auto inv_standardDeviation = 1. / standardDeviation; + for (auto & it : patchVector) + { + it = (it - mean) * inv_standardDeviation; + } + } + } + return patchVector; +} + +template +void +NonLocalPatchBasedImageFilter::GetMeanAndStandardDeviationOfVectorizedImagePatch( + const InputImagePixelVectorType & patchVector, + RealType & mean, + RealType & standardDeviation) +{ + RealType sum = 0.0; + RealType sumOfSquares = 0.0; + RealType count = 0.0; + + for (const auto it : patchVector) + { + if (std::isfinite(it)) // Silently skip non-finite values used to indicate + // out-of-bounds. + { + sum += it; + sumOfSquares += itk::Math::sqr(it); + count += itk::NumericTraits::OneValue(); + } + } + mean = sum / count; + standardDeviation = + std::sqrt((sumOfSquares - count * itk::Math::sqr(mean)) / (count - itk::NumericTraits::OneValue())); +} + +template +typename NonLocalPatchBasedImageFilter::RealType +NonLocalPatchBasedImageFilter::ComputeNeighborhoodPatchSimilarity( + const InputImageList & imageList, + const IndexType index, + const InputImagePixelVectorType & patchVectorY, + const bool useOnlyFirstImage) +{ + unsigned int numberOfImagesToUse = imageList.size(); + if (useOnlyFirstImage) + { + numberOfImagesToUse = 1; + } + + RealType sumX = 0.0; + RealType sumOfSquaresX = 0.0; + RealType sumOfSquaredDifferencesXY = 0.0; + RealType sumXY = 0.0; + RealType N = 0.0; + + SizeValueType count = 0; + for (SizeValueType i = 0; i < numberOfImagesToUse; i++) + { + for (SizeValueType j = 0; j < this->m_NeighborhoodPatchSize; j++) + { + IndexType neighborhoodIndex = index + this->m_NeighborhoodPatchOffsetList[j]; + + bool isInBounds = this->m_TargetImageRegion.IsInside(neighborhoodIndex); + if (isInBounds && std::isfinite(patchVectorY[count])) + { + auto x = static_cast(imageList[i]->GetPixel(neighborhoodIndex)); + auto y = static_cast(patchVectorY[count]); + + sumX += x; + sumOfSquaresX += itk::Math::sqr(x); + sumXY += (x * y); + + sumOfSquaredDifferencesXY += itk::Math::sqr(y - x); + N += itk::NumericTraits::OneValue(); + } + ++count; + } + } + + // If we are on the boundary, a neighborhood patch might not overlap + // with the image. If we have 2 voxels or less for a neighborhood patch + // we don't consider it to be a suitable match. + if (N < static_cast(3.0)) + { + return NumericTraits::max(); + } + if (this->m_SimilarityMetric == SimilarityMetricEnum::PEARSON_CORRELATION) + { + const RealType varianceX = sumOfSquaresX - itk::Math::sqr(sumX) / N; + const RealType measure = + (varianceX > std::numeric_limits::epsilon()) ? itk::Math::sqr(sumXY) / varianceX : 0.; + if (sumXY > 0) + { + return -measure; + } + else + { + return measure; + } + } + else if (this->m_SimilarityMetric == SimilarityMetricEnum::MEAN_SQUARES) + { + return (sumOfSquaredDifferencesXY / N); + } + else + { + itkExceptionMacro("Unrecognized similarity metric."); + } +} + +template +void +NonLocalPatchBasedImageFilter::PrintSelf(std::ostream & os, Indent indent) const +{ + Superclass::PrintSelf(os, indent); + + if (this->m_SimilarityMetric == SimilarityMetricEnum::PEARSON_CORRELATION) + { + os << "Using Pearson correlation to measure the patch similarity." << std::endl; + } + else if (this->m_SimilarityMetric == SimilarityMetricEnum::MEAN_SQUARES) + { + os << "Using mean squares to measure the patch similarity." << std::endl; + } + + os << indent << "Neighborhood search radius = " << this->m_NeighborhoodSearchRadius << std::endl; + os << indent << "Neighborhood patch radius = " << this->m_NeighborhoodPatchRadius << std::endl; +} + +} // end namespace itk + +#endif diff --git a/Modules/Filtering/AdaptiveDenoising/include/itkVarianceImageFilter.h b/Modules/Filtering/AdaptiveDenoising/include/itkVarianceImageFilter.h new file mode 100644 index 00000000000..cf5133b1521 --- /dev/null +++ b/Modules/Filtering/AdaptiveDenoising/include/itkVarianceImageFilter.h @@ -0,0 +1,102 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkVarianceImageFilter_h +#define itkVarianceImageFilter_h + +#include "itkBoxImageFilter.h" +#include "itkImage.h" +#include "itkNumericTraits.h" + +namespace itk +{ +/** \class VarianceImageFilter + * \brief Applies a variance filter to an image + * + * Computes an image where a given pixel is the sample variance of the + * pixels in a neighborhood about the corresponding input pixel. + * + * A variance filter is a nonlinear (quadratic) neighborhood filter. + * + * \ingroup AdaptiveDenoising + */ +template +class ITK_TEMPLATE_EXPORT VarianceImageFilter final : public BoxImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(VarianceImageFilter); + + /** Extract dimension from input and output image. */ + static constexpr unsigned int InputImageDimension = TInputImage::ImageDimension; + static constexpr unsigned int OutputImageDimension = TOutputImage::ImageDimension; + + /** Convenient typedefs for simplifying declarations. */ + using InputImageType = TInputImage; + using OutputImageType = TOutputImage; + + /** Standard class typedefs. */ + using Self = VarianceImageFilter; + using Superclass = BoxImageFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkOverrideGetNameOfClassMacro(VarianceImageFilter); + + /** Image typedef support. */ + using InputPixelType = typename InputImageType::PixelType; + using OutputPixelType = typename OutputImageType::PixelType; + using InputRealType = typename NumericTraits::RealType; + + using InputImageRegionType = typename InputImageType::RegionType; + using OutputImageRegionType = typename OutputImageType::RegionType; + + using InputSizeType = typename InputImageType::SizeType; + +#ifdef ITK_USE_CONCEPT_CHECKING + // Begin concept checking + itkConceptMacro(InputHasNumericTraitsCheck, (Concept::HasNumericTraits)); + // End concept checking +#endif + +protected: + VarianceImageFilter(); + ~VarianceImageFilter() override = default; + + /** VarianceImageFilter can be implemented as a multithreaded filter. + * Therefore, this implementation provides a ThreadedGenerateData() + * routine which is called for each processing thread. The output + * image data is allocated automatically by the superclass prior to + * calling ThreadedGenerateData(). ThreadedGenerateData can only + * write to the portion of the output image specified by the + * parameter "outputRegionForThread" + * + * \sa BoxImageFilter::ThreadedGenerateData(), + * BoxImageFilter::GenerateData() */ + void + ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, ThreadIdType threadId) override; +}; +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkVarianceImageFilter.hxx" +#endif + +#endif diff --git a/Modules/Filtering/AdaptiveDenoising/include/itkVarianceImageFilter.hxx b/Modules/Filtering/AdaptiveDenoising/include/itkVarianceImageFilter.hxx new file mode 100644 index 00000000000..ddb67125485 --- /dev/null +++ b/Modules/Filtering/AdaptiveDenoising/include/itkVarianceImageFilter.hxx @@ -0,0 +1,101 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkVarianceImageFilter_hxx +#define itkVarianceImageFilter_hxx + +#include "itkConstNeighborhoodIterator.h" +#include "itkNeighborhoodInnerProduct.h" +#include "itkImageRegionIterator.h" +#include "itkNeighborhoodAlgorithm.h" +#include "itkOffset.h" +#include "itkProgressReporter.h" + +namespace itk +{ +template +VarianceImageFilter::VarianceImageFilter() +{ + this->DynamicMultiThreadingOff(); +} + +template +void +VarianceImageFilter::ThreadedGenerateData( + const OutputImageRegionType & outputRegionForThread, + ThreadIdType threadId) +{ + unsigned int i; + + ZeroFluxNeumannBoundaryCondition nbc; + + ConstNeighborhoodIterator bit; + ImageRegionIterator it; + + // Allocate output + typename OutputImageType::Pointer output = this->GetOutput(); + typename InputImageType::ConstPointer input = this->GetInput(); + + // Find the data-set boundary "faces" + typename NeighborhoodAlgorithm::ImageBoundaryFacesCalculator::FaceListType faceList; + NeighborhoodAlgorithm::ImageBoundaryFacesCalculator bC; + faceList = bC(input, outputRegionForThread, this->GetRadius()); + + typename NeighborhoodAlgorithm::ImageBoundaryFacesCalculator::FaceListType::iterator fit; + + // support progress methods/callbacks + ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels()); + + InputRealType sum; + InputRealType sumOfSquares; + + // Process each of the boundary faces. These are N-d regions which border + // the edge of the buffer. + for (fit = faceList.begin(); fit != faceList.end(); ++fit) + { + bit = ConstNeighborhoodIterator(this->GetRadius(), input, *fit); + unsigned int neighborhoodSize = bit.Size(); + it = ImageRegionIterator(output, *fit); + bit.OverrideBoundaryCondition(&nbc); + bit.GoToBegin(); + + while (!bit.IsAtEnd()) + { + sum = NumericTraits::ZeroValue(); + sumOfSquares = NumericTraits::ZeroValue(); + + for (i = 0; i < neighborhoodSize; ++i) + { + sum += static_cast(bit.GetPixel(i)); + sumOfSquares += itk::Math::sqr(static_cast(bit.GetPixel(i))); + } + + // get the variance value + const double num = static_cast(neighborhoodSize); + OutputPixelType var = (sumOfSquares - (itk::Math::sqr(sum) / num)) / (num - 1.0); + + it.Set(var); + + ++bit; + ++it; + progress.CompletedPixel(); + } + } +} +} // end namespace itk + +#endif diff --git a/Modules/Filtering/AdaptiveDenoising/itk-module.cmake b/Modules/Filtering/AdaptiveDenoising/itk-module.cmake new file mode 100644 index 00000000000..2f10fe55fea --- /dev/null +++ b/Modules/Filtering/AdaptiveDenoising/itk-module.cmake @@ -0,0 +1,26 @@ +# the top-level README is used for describing this module, just +# re-used it for documentation here +# itk_module() defines the module dependencies in AdaptiveDenoising +# AdaptiveDenoising depends on ITKCommon +# The testing module in AdaptiveDenoising depends on ITKTestKernel +# and ITKMetaIO(besides AdaptiveDenoising and ITKCore) +# By convention those modules outside of ITK are not prefixed with +# ITK. + +# define the dependencies of the include module and the tests +itk_module( + AdaptiveDenoising + DEPENDS + ITKCommon + ITKSmoothing + ITKStatistics + ITKImageStatistics + COMPILE_DEPENDS + ITKImageSources + TEST_DEPENDS + ITKTestKernel + ITKMetaIO + DESCRIPTION "Module ingested from upstream." + EXCLUDE_FROM_DEFAULT + ENABLE_SHARED +) diff --git a/Modules/Filtering/AdaptiveDenoising/src/CMakeLists.txt b/Modules/Filtering/AdaptiveDenoising/src/CMakeLists.txt new file mode 100644 index 00000000000..b8b206626a5 --- /dev/null +++ b/Modules/Filtering/AdaptiveDenoising/src/CMakeLists.txt @@ -0,0 +1,3 @@ +set(AdaptiveDenoising_SRCS itkNonLocalPatchBasedImageFilterEnums.cxx) + +itk_module_add_library(AdaptiveDenoising ${AdaptiveDenoising_SRCS}) diff --git a/Modules/Filtering/AdaptiveDenoising/src/itkNonLocalPatchBasedImageFilterEnums.cxx b/Modules/Filtering/AdaptiveDenoising/src/itkNonLocalPatchBasedImageFilterEnums.cxx new file mode 100644 index 00000000000..931c9f32037 --- /dev/null +++ b/Modules/Filtering/AdaptiveDenoising/src/itkNonLocalPatchBasedImageFilterEnums.cxx @@ -0,0 +1,42 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#include "itkNonLocalPatchBasedImageFilter.h" + + +namespace itk +{ + +/** Print enum values */ +std::ostream & +operator<<(std::ostream & out, const NonLocalPatchBasedImageFilterEnums::SimilarityMetric value) +{ + return out << [value] { + switch (value) + { + case NonLocalPatchBasedImageFilterEnums::SimilarityMetric::PEARSON_CORRELATION: + return "itk::NonLocalPatchBasedImageFilterEnums::SimilarityMetric::PEARSON_CORRELATION"; + case NonLocalPatchBasedImageFilterEnums::SimilarityMetric::MEAN_SQUARES: + return "itk::NonLocalPatchBasedImageFilterEnums::SimilarityMetric::MEAN_SQUARES"; + default: + return "INVALID VALUE FOR itk::NonLocalPatchBasedImageFilterEnums::SimilarityMetric"; + } + }(); +} + +} // end namespace itk diff --git a/Modules/Filtering/AdaptiveDenoising/test/Baseline/r16denoised_mean_squares.nrrd.cid b/Modules/Filtering/AdaptiveDenoising/test/Baseline/r16denoised_mean_squares.nrrd.cid new file mode 100644 index 00000000000..2d29f602e36 --- /dev/null +++ b/Modules/Filtering/AdaptiveDenoising/test/Baseline/r16denoised_mean_squares.nrrd.cid @@ -0,0 +1 @@ +bafkreigltzq46r5ojiiyel37egvlosnxxgyth6ppwckzjtuxyywtlijoii diff --git a/Modules/Filtering/AdaptiveDenoising/test/Baseline/r16denoised_pearson_correlation.nrrd.cid b/Modules/Filtering/AdaptiveDenoising/test/Baseline/r16denoised_pearson_correlation.nrrd.cid new file mode 100644 index 00000000000..dfcb5cf3f88 --- /dev/null +++ b/Modules/Filtering/AdaptiveDenoising/test/Baseline/r16denoised_pearson_correlation.nrrd.cid @@ -0,0 +1 @@ +bafkreickjyvdwmq6eyz76e7cukkegm4dd75raimxj4pag2o6bkanvbofwm diff --git a/Modules/Filtering/AdaptiveDenoising/test/CMakeLists.txt b/Modules/Filtering/AdaptiveDenoising/test/CMakeLists.txt new file mode 100644 index 00000000000..6bdf0be53c7 --- /dev/null +++ b/Modules/Filtering/AdaptiveDenoising/test/CMakeLists.txt @@ -0,0 +1,31 @@ +itk_module_test() + +set(AdaptiveDenoisingTests itkAdaptiveNonLocalMeansDenoisingImageFilterTest.cxx) + +createtestdriver(AdaptiveDenoising "${AdaptiveDenoising-Test_LIBRARIES}" "${AdaptiveDenoisingTests}") + +itk_add_test( + NAME AdaptiveNonLocalMeansDenoisingImageFilterTest1 + COMMAND + AdaptiveDenoisingTestDriver + --compare + DATA{Baseline/r16denoised_pearson_correlation.nrrd} + ${ITK_TEST_OUTPUT_DIR}/r16denoised_pearson_correlation.nrrd + itkAdaptiveNonLocalMeansDenoisingImageFilterTest + DATA{Input/r16slice.nrrd} + ${ITK_TEST_OUTPUT_DIR}/r16denoised_pearson_correlation.nrrd + 0 +) + +itk_add_test( + NAME AdaptiveNonLocalMeansDenoisingImageFilterTest2 + COMMAND + AdaptiveDenoisingTestDriver + --compare + DATA{Baseline/r16denoised_mean_squares.nrrd} + ${ITK_TEST_OUTPUT_DIR}/r16denoised_mean_squares.nrrd + itkAdaptiveNonLocalMeansDenoisingImageFilterTest + DATA{Input/r16slice.nrrd} + ${ITK_TEST_OUTPUT_DIR}/r16denoised_mean_squares.nrrd + 1 +) diff --git a/Modules/Filtering/AdaptiveDenoising/test/Input/r16slice.nrrd.cid b/Modules/Filtering/AdaptiveDenoising/test/Input/r16slice.nrrd.cid new file mode 100644 index 00000000000..e03b4e30a39 --- /dev/null +++ b/Modules/Filtering/AdaptiveDenoising/test/Input/r16slice.nrrd.cid @@ -0,0 +1 @@ +bafkreiafwnstu4vzkrl7gpt73tl3q53ywpdylbs2d4emezgzutj2thoyfy diff --git a/Modules/Filtering/AdaptiveDenoising/test/itkAdaptiveNonLocalMeansDenoisingImageFilterTest.cxx b/Modules/Filtering/AdaptiveDenoising/test/itkAdaptiveNonLocalMeansDenoisingImageFilterTest.cxx new file mode 100644 index 00000000000..078c731ea5a --- /dev/null +++ b/Modules/Filtering/AdaptiveDenoising/test/itkAdaptiveNonLocalMeansDenoisingImageFilterTest.cxx @@ -0,0 +1,181 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#include +#include "itkAdaptiveNonLocalMeansDenoisingImageFilter.h" + +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" +#include "itkTestingMacros.h" +#include "itkMath.h" + +template +class CommandProgressUpdate : public itk::Command +{ +public: + using Self = CommandProgressUpdate; + using Superclass = itk::Command; + using Pointer = itk::SmartPointer>; + itkNewMacro(CommandProgressUpdate); + +protected: + CommandProgressUpdate() = default; + + using FilterType = TFilter; + + unsigned int m_CurrentProgress{ 0 }; + +public: + void + Execute(itk::Object * caller, const itk::EventObject & event) override + { + auto * po = dynamic_cast(caller); + if (!po) + return; + if (typeid(event) == typeid(itk::ProgressEvent)) + { + if (this->m_CurrentProgress < 99) + { + this->m_CurrentProgress++; + if (this->m_CurrentProgress % 10 == 0) + { + std::cout << this->m_CurrentProgress << std::flush; + } + else + { + std::cout << "*" << std::flush; + } + } + } + } + + void + Execute(const itk::Object * object, const itk::EventObject & event) override + { + auto * po = dynamic_cast(const_cast(object)); + if (!po) + return; + + if (typeid(event) == typeid(itk::ProgressEvent)) + { + if (this->m_CurrentProgress < 99) + { + this->m_CurrentProgress++; + if (this->m_CurrentProgress % 10 == 0) + { + std::cout << this->m_CurrentProgress << std::flush; + } + else + { + std::cout << "*" << std::flush; + } + } + } + } +}; + +int +itkAdaptiveNonLocalMeansDenoisingImageFilterTest(int argc, char * argv[]) +{ + if (argc < 4) + { + std::cerr << "Missing parameters." << std::endl; + std::cerr << "Usage: " << itkNameOfTestExecutableMacro(argv) << " inputImage" + << " outputImage" + << " similarityMetric (0: PEARSON_CORRELATION; 1: MEAN_SQUARES)" << std::endl; + return EXIT_FAILURE; + } + + constexpr unsigned int Dimension = 2; + using PixelType = float; + using ImageType = itk::Image; + + using ReaderType = itk::ImageFileReader; + ReaderType::Pointer reader = ReaderType::New(); + reader->SetFileName(argv[1]); + + using DenoiserType = itk::AdaptiveNonLocalMeansDenoisingImageFilter; + DenoiserType::Pointer filter = DenoiserType::New(); + + ITK_EXERCISE_BASIC_OBJECT_METHODS(filter, AdaptiveNonLocalMeansDenoisingImageFilter, NonLocalPatchBasedImageFilter); + + filter->SetInput(reader->GetOutput()); + filter->SetUseRicianNoiseModel(false); + ITK_TEST_SET_GET_VALUE(false, filter->GetUseRicianNoiseModel()); + + DenoiserType::NeighborhoodRadiusType neighborhoodPatchRadius; + DenoiserType::NeighborhoodRadiusType neighborhoodSearchRadius; + + neighborhoodPatchRadius.Fill(1); + neighborhoodSearchRadius.Fill(2); + + filter->SetNeighborhoodSearchRadius(neighborhoodSearchRadius); + filter->SetNeighborhoodPatchRadius(neighborhoodPatchRadius); + + DenoiserType::NeighborhoodRadiusType neighborhoodRadiusForLocalMeanAndVariance; + neighborhoodRadiusForLocalMeanAndVariance.Fill(1); + + filter->SetNeighborhoodRadiusForLocalMeanAndVariance(neighborhoodRadiusForLocalMeanAndVariance); + + filter->SetEpsilon(0.00001f); + ITK_TEST_EXPECT_TRUE(itk::Math::FloatAlmostEqual(0.00001f, filter->GetEpsilon())); + + filter->SetMeanThreshold(0.95f); + ITK_TEST_EXPECT_TRUE(itk::Math::FloatAlmostEqual(0.95f, filter->GetMeanThreshold())); + + filter->SetVarianceThreshold(0.5f); + ITK_TEST_EXPECT_TRUE(itk::Math::FloatAlmostEqual(0.5f, filter->GetVarianceThreshold())); + + filter->SetSmoothingFactor(1.0f); + ITK_TEST_EXPECT_TRUE(itk::Math::FloatAlmostEqual(1.0f, filter->GetSmoothingFactor())); + + filter->SetSmoothingVariance(2.0f); + ITK_TEST_EXPECT_TRUE(itk::Math::FloatAlmostEqual(2.0f, filter->GetSmoothingVariance())); + + auto similarityMetric = static_cast(std::atoi(argv[3])); + filter->SetSimilarityMetric(similarityMetric); + ITK_TEST_SET_GET_VALUE(similarityMetric, filter->GetSimilarityMetric()); + + using CommandType = CommandProgressUpdate; + CommandType::Pointer observer = CommandType::New(); + filter->AddObserver(itk::ProgressEvent(), observer); + + ITK_TRY_EXPECT_NO_EXCEPTION(filter->Update()); + + using WriterType = itk::ImageFileWriter; + WriterType::Pointer writer = WriterType::New(); + writer->SetFileName(argv[2]); + writer->SetInput(filter->GetOutput()); + + ITK_TRY_EXPECT_NO_EXCEPTION(writer->Update()); + + + // Test streaming enumeration for NonLocalPatchBasedImageFilterEnums::SimilarityMetric elements + const std::set allSimilarityMetric{ + itk::NonLocalPatchBasedImageFilterEnums::SimilarityMetric::PEARSON_CORRELATION, + itk::NonLocalPatchBasedImageFilterEnums::SimilarityMetric::MEAN_SQUARES + }; + for (const auto & ee : allSimilarityMetric) + { + std::cout << "STREAMED ENUM VALUE NonLocalPatchBasedImageFilterEnums::SimilarityMetric: " << ee << std::endl; + } + + + std::cout << "Test finished" << std::endl; + return EXIT_SUCCESS; +} diff --git a/Modules/Filtering/AdaptiveDenoising/wrapping/CMakeLists.txt b/Modules/Filtering/AdaptiveDenoising/wrapping/CMakeLists.txt new file mode 100644 index 00000000000..c2e65ecce77 --- /dev/null +++ b/Modules/Filtering/AdaptiveDenoising/wrapping/CMakeLists.txt @@ -0,0 +1,11 @@ +itk_wrap_module(AdaptiveDenoising) + +set( + WRAPPER_SUBMODULE_ORDER + itkVarianceImageFilter + itkNonLocalPatchBasedImageFilter + itkAdaptiveNonLocalMeansDenoisingImageFilter +) + +itk_auto_load_submodules() +itk_end_wrap_module() diff --git a/Modules/Filtering/AdaptiveDenoising/wrapping/itkAdaptiveNonLocalMeansDenoisingImageFilter.wrap b/Modules/Filtering/AdaptiveDenoising/wrapping/itkAdaptiveNonLocalMeansDenoisingImageFilter.wrap new file mode 100644 index 00000000000..4ecd0e25c81 --- /dev/null +++ b/Modules/Filtering/AdaptiveDenoising/wrapping/itkAdaptiveNonLocalMeansDenoisingImageFilter.wrap @@ -0,0 +1,3 @@ +itk_wrap_class("itk::AdaptiveNonLocalMeansDenoisingImageFilter" POINTER) +itk_wrap_image_filter("${WRAP_ITK_SCALAR}" 2) +itk_end_wrap_class() diff --git a/Modules/Filtering/AdaptiveDenoising/wrapping/itkNonLocalPatchBasedImageFilter.wrap b/Modules/Filtering/AdaptiveDenoising/wrapping/itkNonLocalPatchBasedImageFilter.wrap new file mode 100644 index 00000000000..859f4e064c0 --- /dev/null +++ b/Modules/Filtering/AdaptiveDenoising/wrapping/itkNonLocalPatchBasedImageFilter.wrap @@ -0,0 +1,8 @@ +set(WRAPPER_AUTO_INCLUDE_HEADERS OFF) +itk_wrap_include("itkNonLocalPatchBasedImageFilter.h") + +itk_wrap_simple_class("itk::NonLocalPatchBasedImageFilterEnums" ENUM) + +itk_wrap_class("itk::NonLocalPatchBasedImageFilter" POINTER) +itk_wrap_image_filter("${WRAP_ITK_SCALAR}" 2) +itk_end_wrap_class() diff --git a/Modules/Filtering/AdaptiveDenoising/wrapping/itkVarianceImageFilter.wrap b/Modules/Filtering/AdaptiveDenoising/wrapping/itkVarianceImageFilter.wrap new file mode 100644 index 00000000000..ffc8554bdd4 --- /dev/null +++ b/Modules/Filtering/AdaptiveDenoising/wrapping/itkVarianceImageFilter.wrap @@ -0,0 +1,3 @@ +itk_wrap_class("itk::VarianceImageFilter" POINTER) +itk_wrap_image_filter("${WRAP_ITK_SCALAR}" 2) +itk_end_wrap_class() diff --git a/Modules/Remote/AdaptiveDenoising.remote.cmake b/Modules/Remote/AdaptiveDenoising.remote.cmake deleted file mode 100644 index 289dbdfa9ff..00000000000 --- a/Modules/Remote/AdaptiveDenoising.remote.cmake +++ /dev/null @@ -1,61 +0,0 @@ -#-- # Grading Level Criteria Report -#-- EVALUATION DATE: 09/28/2020 -#-- EVALUATORS: Nick Tustison -#-- -#-- ## Compliance level 5 star (AKA ITK main modules, or remote modules that could become core modules) -#-- - [ ] Widespread community dependance -#-- - [ ] Above 90% code coverage -#-- - [ ] CI dashboards and testing monitored rigorously -#-- - [ ] Key API features are exposed in wrapping interface -#-- - [ ] All requirements of Levels 4,3,2,1 -#-- -#-- ## Compliance Level 4 star (Very high-quality code, perhaps small community dependance) -#-- - [ ] Meets all ITK code style standards -#-- - [ X ] No external requirements beyond those needed by ITK proper -#-- - [ ] Builds and passes tests on all supported platforms within 1 month of each core tagged release -#-- - [ ] Windows Shared Library Build with Visual Studio -#-- - [ ] Mac with clang compiller -#-- - [ ] Linux with gcc compiler -#-- - [ X ] Active developer community dedicated to maintaining code-base -#-- - [ ] 75% code coverage demonstrated for testing suite -#-- - [ X ] Continuous integration testing performed -#-- - [ ] All requirements of Levels 3,2,1 -#-- -#-- ## Compliance Level 3 star (Quality beta code) -#-- - [ X ] API | executable interface is considered mostly stable and feature complete -#-- - [ X ] 10% C0-code coverage demonstrated for testing suite -#-- - [ X ] Some tests exist and pass on at least some platform -#-- - [ X ] All requirements of Levels 2,1 -#-- -#-- ## Compliance Level 2 star (Alpha code feature API development or niche community/execution environment dependance ) -#-- - [ X ] Compiles for at least 1 niche set of execution envirionments, and perhaps others -#-- (may depend on specific external tools like a java environment, or specific external libraries to work ) -#-- - [ X ] All requirements of Levels 1 -#-- -#-- ## Compliance Level 1 star (Pre-alpha features under development and code of unknown quality) -#-- - [ X ] Code complies on at least 1 platform -#-- -#-- ## Compliance Level 0 star ( Code/Feature of known poor-quality or deprecated status ) -#-- - [ ] Code reviewed and explicitly identified as not recommended for use -#-- -#-- ### Please document here any justification for the criteria above -# Code style enforced by clang-format on 2020-02-19, and clang-tidy modernizations completed - -itk_fetch_module( - AdaptiveDenoising - "A spatially adaptive denoising image filter using non-local means. - A patch-based framework for new ITK functionality: Joint fusion, denoising, - and non-local super-resolution. - Tustison N., Avants B., Wang H., Xie L., Coupe P., Yushkevich P., Manjon J. - Insight Journal. - https://doi.org/10.54294/ywuouh. - Two Luis Miguel fans walk into a bar in Nagoya ---> (yada, yada, yada) - ---> an ITK-implementation of a popular patch-based denoising filter. - Tustison N., Manjon J.V. - Insight Journal. - https://doi.org/10.54294/9f5wt3. - " - MODULE_COMPLIANCE_LEVEL 3 - GIT_REPOSITORY https://github.com/ntustison/ITKAdaptiveDenoising.git - GIT_TAG 99db4507f05ea933600c41ccf0d53b4202f4cab2 - ) diff --git a/pyproject.toml b/pyproject.toml index c73ae38c83d..3f7fc02a387 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -52,6 +52,7 @@ cmd = '''cmake \ -DITK_USE_CCACHE:BOOL=ON \ -DCMAKE_C_COMPILER_LAUNCHER:STRING=ccache \ -DCMAKE_CXX_COMPILER_LAUNCHER:STRING=ccache \ + -DModule_AdaptiveDenoising:BOOL=ON \ -DModule_AnisotropicDiffusionLBR:BOOL=ON \ -DModule_Cuberille:BOOL=ON \ -DModule_Montage:BOOL=ON \