ITK  6.0.0
Insight Toolkit
itkPatchBasedDenoisingImageFilter.h
Go to the documentation of this file.
1/*=========================================================================
2 *
3 * Copyright NumFOCUS
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License");
6 * you may not use this file except in compliance with the License.
7 * You may obtain a copy of the License at
8 *
9 * https://www.apache.org/licenses/LICENSE-2.0.txt
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 *
17 *=========================================================================*/
18#ifndef itkPatchBasedDenoisingImageFilter_h
19#define itkPatchBasedDenoisingImageFilter_h
20
24#include "itkVector.h"
25#include "itkVectorImage.h"
26#include "itkRGBPixel.h"
27#include "itkRGBAPixel.h"
29#include "itkFixedArray.h"
30#include "itkMatrix.h"
32#include <type_traits>
33
34#include <vector>
35#include "ITKDenoisingExport.h"
36
37namespace itk
38{
61template <typename TInputImage, typename TOutputImage>
62class ITK_TEMPLATE_EXPORT PatchBasedDenoisingImageFilter
63 : public PatchBasedDenoisingBaseImageFilter<TInputImage, TOutputImage>
64{
65public:
66 ITK_DISALLOW_COPY_AND_MOVE(PatchBasedDenoisingImageFilter);
74 using typename Superclass::OutputImagePointer;
75
77 itkNewMacro(Self);
78
80 itkOverrideGetNameOfClassMacro(PatchBasedDenoisingImageFilter);
81
83 using typename Superclass::InputImageType;
84 using typename Superclass::OutputImageType;
85
87 static constexpr unsigned int ImageDimension = Superclass::ImageDimension;
88
91
95
98 using typename Superclass::PixelType;
99 using typename Superclass::PixelValueType;
100
106
108 using typename Superclass::ListAdaptorType;
109 using typename Superclass::PatchRadiusType;
110 using typename Superclass::InputImagePatchIterator;
112 using typename Superclass::PatchWeightsType;
113
118
126 using EigenValuesCacheType = std::vector<EigenValuesArrayType>;
127 using EigenVectorsCacheType = std::vector<EigenVectorsMatrixType>;
128
130 {
140 };
141
146 itkSetMacro(UseSmoothDiscPatchWeights, bool);
147 itkBooleanMacro(UseSmoothDiscPatchWeights);
148 itkGetConstMacro(UseSmoothDiscPatchWeights, bool);
155 void
157 itkGetConstMacro(KernelBandwidthSigma, RealArrayType);
164 itkSetClampMacro(KernelBandwidthFractionPixelsForEstimation, double, 0.01, 1.0);
165 itkGetConstReferenceMacro(KernelBandwidthFractionPixelsForEstimation, double);
170 itkSetMacro(ComputeConditionalDerivatives, bool);
171 itkBooleanMacro(ComputeConditionalDerivatives);
172 itkGetConstMacro(ComputeConditionalDerivatives, bool);
188 itkSetMacro(UseFastTensorComputations, bool);
189 itkBooleanMacro(UseFastTensorComputations);
190 itkGetConstMacro(UseFastTensorComputations, bool);
194 static constexpr unsigned int MaxSigmaUpdateIterations = 20;
195
202 itkSetClampMacro(KernelBandwidthMultiplicationFactor, double, 0.01, 100);
203 itkGetConstReferenceMacro(KernelBandwidthMultiplicationFactor, double);
209 void
210 SetNoiseSigma(const RealType & sigma);
211
212 itkGetConstMacro(NoiseSigma, RealType);
213
215 itkSetObjectMacro(Sampler, BaseSamplerType);
216 itkGetModifiableObjectMacro(Sampler, BaseSamplerType);
220 itkGetConstMacro(NumIndependentComponents, unsigned int);
221
222protected:
225 void
226 PrintSelf(std::ostream & os, Indent indent) const override;
227
229 virtual void
231
233 void
235
236 void
238
239 void
241
242 template <typename T, typename U = void>
243 using DisableIfMultiComponent = typename std::enable_if<std::is_same_v<T, typename NumericTraits<T>::ValueType>, U>;
244
245 template <typename T, typename U = void>
246 using EnableIfMultiComponent = typename std::enable_if<!std::is_same_v<T, typename NumericTraits<T>::ValueType>, U>;
247
248
255 template <typename T>
257 GetComponent(const T pix, unsigned int itkNotUsed(idx)) const
258 {
259 // The enable if idiom is used to overload this method for both
260 // scalars and multi-component types. By exploiting that
261 // NumericTraits' ValueType type alias (defines the per-element type
262 // for multi-component types ) is different then the parameterize
263 // type, the bracket operator is used only for multi-component
264 // types.
265 return pix;
266 }
267
268 template <typename T>
269 typename EnableIfMultiComponent<T, typename NumericTraits<T>::ValueType>::type
270 GetComponent(const T & pix, unsigned int idx) const
271 {
272 return pix[idx];
273 }
274
276 template <typename T>
277 void
279 unsigned int itkNotUsed(idx),
281 {
282 pix = val;
283 }
284
285 template <typename T>
286 void
287 SetComponent(T & pix, unsigned int idx, typename EnableIfMultiComponent<T, RealValueType>::type val) const
288 {
289 pix[idx] = val;
290 }
291
295 void
297 {
298 if (this->GetComponentSpace() == Superclass::ComponentSpaceEnum::RIEMANNIAN)
299 {
300 DispatchedRiemannianMinMax(img);
301 }
302 else
303 {
304 DispatchedArrayMinMax(img);
305 }
306 }
307
308 template <typename TImageType>
309 typename DisableIfMultiComponent<typename TImageType::PixelType>::type
310 ComputeMinMax(const TImageType * img)
311 {
312 DispatchedMinMax(img);
313 }
314
315 template <typename TImageType>
316 typename EnableIfMultiComponent<typename TImageType::PixelType>::type
317 ComputeMinMax(const TImageType * img)
318 {
319 DispatchedArrayMinMax(img);
320 }
321
330 void
333 const RealArrayType & weight,
334 bool useCachedComputations,
335 SizeValueType cacheIndex,
336 EigenValuesCacheType & eigenValsCache,
337 EigenVectorsCacheType & eigenVecsCache,
338 RealType & diff,
339 RealArrayType & norm)
340 {
341 if (this->GetComponentSpace() == Superclass::ComponentSpaceEnum::RIEMANNIAN)
342 {
343 ComputeLogMapAndWeightedSquaredGeodesicDifference(
344 a, b, weight, useCachedComputations, cacheIndex, eigenValsCache, eigenVecsCache, diff, norm);
345 }
346 else
347 {
348 ComputeSignedEuclideanDifferenceAndWeightedSquaredNorm(
349 a, b, weight, useCachedComputations, cacheIndex, eigenValsCache, eigenVecsCache, diff, norm);
350 }
351 }
354 template <typename PixelT>
355 void
357 const PixelT & b,
358 const RealArrayType & weight,
359 bool useCachedComputations,
360 SizeValueType cacheIndex,
361 EigenValuesCacheType & eigenValsCache,
362 EigenVectorsCacheType & eigenVecsCache,
363 RealType & diff,
364 RealArrayType & norm)
365 {
366 ComputeSignedEuclideanDifferenceAndWeightedSquaredNorm(
367 a, b, weight, useCachedComputations, cacheIndex, eigenValsCache, eigenVecsCache, diff, norm);
368 }
369
373 RealType
375 {
376 if (this->GetComponentSpace() == Superclass::ComponentSpaceEnum::RIEMANNIAN)
377 {
378 return this->AddExponentialMapUpdate(a, b);
379 }
380 else
381 {
382 return this->AddEuclideanUpdate(a, b);
383 }
384 }
387 template <typename RealT>
388 RealType
389 AddUpdate(const RealT & a, const RealType & b)
390 {
391 return this->AddEuclideanUpdate(a, b);
392 }
393
394 virtual void
396
397 void
398 Initialize() override;
399
400 virtual void
402
403 void
405
406 virtual void
408
409 void
411
412 void
413 ComputeKernelBandwidthUpdate() override; // derived from base class;
414
415 // define here
416
417 virtual ThreadDataStruct
419 const int itkNotUsed(threadId),
420 ThreadDataStruct threadData);
421
422 virtual RealArrayType
424
425 void
427
428 virtual ThreadDataStruct
430 const int threadId,
431 ThreadDataStruct threadData);
432
433 virtual RealType
435 typename ListAdaptorType::Pointer & inList,
436 BaseSamplerPointer & sampler,
437 ThreadDataStruct & threadData);
438
439 void
440 ApplyUpdate() override;
441
442 virtual void
443 ThreadedApplyUpdate(const InputImageRegionType & regionToProcess, const int itkNotUsed(threadId));
444
445 void
447
448 virtual void
449 SetThreadData(int threadId, const ThreadDataStruct & data);
450
451 virtual ThreadDataStruct
452 GetThreadData(int threadId);
453
454private:
459
464
469
470 template <typename TInputImageType>
471 void
472 DispatchedMinMax(const TInputImageType * img);
473
474 template <typename TInputImageType>
475 void
476 DispatchedArrayMinMax(const TInputImageType * img);
477
478 template <typename TInputImageType>
479 void
480 DispatchedVectorMinMax(const TInputImageType * img);
481
482 template <typename TInputImageType>
483 void
484 DispatchedRiemannianMinMax(const TInputImageType * img);
485
490
493 const int itkNotUsed(threadId),
494 const InputImageType * img,
495 ThreadDataStruct threadData);
496
497 virtual void
499
500 void
502 const PixelType & b,
503 const RealArrayType & weight,
504 bool useCachedComputations,
505 SizeValueType cacheIndex,
506 EigenValuesCacheType & eigenValsCache,
507 EigenVectorsCacheType & eigenVecsCache,
508 RealType & diff,
509 RealArrayType & norm);
510
512 void
514 const DiffusionTensor3D<PixelValueType> & spdMatrixB,
515 const RealArrayType & weight,
516 bool useCachedComputations,
517 SizeValueType cacheIndex,
518 EigenValuesCacheType & eigenValsCache,
519 EigenVectorsCacheType & eigenVecsCache,
520 RealType & symMatrixLogMap,
521 RealArrayType & geodesicDist);
522
523 template <typename TensorValueT>
524 void
526 FixedArray<TensorValueT, 3> & eigenVals,
527 Matrix<TensorValueT, 3, 3> & eigenVecs);
528
530 AddEuclideanUpdate(const RealType & a, const RealType & b);
531
535 const DiffusionTensor3D<RealValueType> & symMatrix);
536
538 {
541 };
542
543 std::vector<ThreadDataStruct> m_ThreadData{};
544
546 typename OutputImageType::Pointer m_UpdateBuffer{};
547
548 unsigned int m_NumPixelComponents{ 0 };
549 unsigned int m_NumIndependentComponents{ 0 };
550 unsigned int m_TotalNumberPixels{ 0 };
551
552 bool m_UseSmoothDiscPatchWeights{ true };
553
554 bool m_UseFastTensorComputations{ true };
555
556 RealArrayType m_KernelBandwidthSigma{};
557 bool m_KernelBandwidthSigmaIsSet{ false };
558 RealArrayType m_IntensityRescaleInvFactor{};
559 PixelType m_ZeroPixel{};
560 PixelArrayType m_ImageMin{};
561 PixelArrayType m_ImageMax{};
562 double m_KernelBandwidthFractionPixelsForEstimation{ 0.20 };
563 bool m_ComputeConditionalDerivatives{ false };
564 double m_MinSigma{};
565 double m_MinProbability{};
566 unsigned int m_SigmaUpdateDecimationFactor{};
567 double m_SigmaUpdateConvergenceTolerance{ 0.01 };
568 ShortArrayType m_SigmaConverged{};
569 double m_KernelBandwidthMultiplicationFactor{ 1.0 };
570
571 RealType m_NoiseSigma{};
572 RealType m_NoiseSigmaSquared{};
573 bool m_NoiseSigmaIsSet{ false };
574
576 typename ListAdaptorType::Pointer m_SearchSpaceList{};
577};
578} // end namespace itk
579
580#ifndef ITK_MANUAL_INSTANTIATION
581# include "itkPatchBasedDenoisingImageFilter.hxx"
582#endif
583
584#endif
Represent a diffusion tensor as used in DTI images.
Simulate a standard C array with copy semantics.
Definition: itkFixedArray.h:54
A multi-dimensional iterator templated over image type that walks a region of pixels.
A multi-dimensional iterator templated over image type that walks a region of pixels.
Base class for all process objects that output image data.
typename InputImageType::RegionType InputImageRegionType
Templated n-dimensional image class.
Definition: itkImage.h:89
Control indentation during Print() invocation.
Definition: itkIndent.h:50
A templated class holding a M x N size Matrix.
Definition: itkMatrix.h:53
Define additional traits for native types such as int or float.
Base class for patch-based denoising algorithms.
typename itk::Statistics::ImageToNeighborhoodSampleAdaptor< OutputImageType, BoundaryConditionType > ListAdaptorType
Derived class implementing a specific patch-based denoising algorithm, as detailed below.
void DispatchedArrayMinMax(const TInputImageType *img)
static ITK_THREAD_RETURN_FUNCTION_CALL_CONVENTION RiemannianMinMaxThreaderCallback(void *arg)
void Compute3x3EigenAnalysis(const DiffusionTensor3D< TensorValueT > &spdMatrix, FixedArray< TensorValueT, 3 > &eigenVals, Matrix< TensorValueT, 3, 3 > &eigenVecs)
void PrintSelf(std::ostream &os, Indent indent) const override
static ITK_THREAD_RETURN_FUNCTION_CALL_CONVENTION ComputeImageUpdateThreaderCallback(void *arg)
void DispatchedMinMax(const TInputImageType *img)
DisableIfMultiComponent< typenameTImageType::PixelType >::type ComputeMinMax(const TImageType *img)
void ComputeDifferenceAndWeightedSquaredNorm(const DiffusionTensor3D< PixelValueType > &a, const DiffusionTensor3D< PixelValueType > &b, const RealArrayType &weight, bool useCachedComputations, SizeValueType cacheIndex, EigenValuesCacheType &eigenValsCache, EigenVectorsCacheType &eigenVecsCache, RealType &diff, RealArrayType &norm)
void ComputeMinMax(const Image< DiffusionTensor3D< PixelValueType >, ImageDimension > *img)
virtual ThreadDataStruct ThreadedComputeSigmaUpdate(const InputImageRegionType &regionToProcess, const int, ThreadDataStruct threadData)
static ITK_THREAD_RETURN_FUNCTION_CALL_CONVENTION ApplyUpdateThreaderCallback(void *arg)
DisableIfMultiComponent< T, T >::type GetComponent(const T pix, unsigned int) const
A method to generically get a component.
virtual ThreadDataStruct ThreadedComputeImageUpdate(const InputImageRegionType &regionToProcess, const int threadId, ThreadDataStruct threadData)
virtual void SetThreadData(int threadId, const ThreadDataStruct &data)
void ComputeLogMapAndWeightedSquaredGeodesicDifference(const DiffusionTensor3D< PixelValueType > &spdMatrixA, const DiffusionTensor3D< PixelValueType > &spdMatrixB, const RealArrayType &weight, bool useCachedComputations, SizeValueType cacheIndex, EigenValuesCacheType &eigenValsCache, EigenVectorsCacheType &eigenVecsCache, RealType &symMatrixLogMap, RealArrayType &geodesicDist)
std::vector< EigenVectorsMatrixType > EigenVectorsCacheType
virtual void InitializePatchWeightsSmoothDisc()
typename std::enable_if<!std::is_same_v< T, typename NumericTraits< T >::ValueType >, U > EnableIfMultiComponent
RealType AddUpdate(const RealT &a, const RealType &b)
typename BaseSamplerType::Pointer BaseSamplerPointer
virtual ThreadDataStruct GetThreadData(int threadId)
void ComputeKernelBandwidthUpdate() override
void SetKernelBandwidthSigma(const RealArrayType &kernelSigma)
RealType AddExponentialMapUpdate(const DiffusionTensor3D< RealValueType > &spdMatrix, const DiffusionTensor3D< RealValueType > &symMatrix)
void ComputeDifferenceAndWeightedSquaredNorm(const PixelT &a, const PixelT &b, const RealArrayType &weight, bool useCachedComputations, SizeValueType cacheIndex, EigenValuesCacheType &eigenValsCache, EigenVectorsCacheType &eigenVecsCache, RealType &diff, RealArrayType &norm)
RealType AddEuclideanUpdate(const RealType &a, const RealType &b)
void ComputeSignedEuclideanDifferenceAndWeightedSquaredNorm(const PixelType &a, const PixelType &b, const RealArrayType &weight, bool useCachedComputations, SizeValueType cacheIndex, EigenValuesCacheType &eigenValsCache, EigenVectorsCacheType &eigenVecsCache, RealType &diff, RealArrayType &norm)
typename std::enable_if< std::is_same_v< T, typename NumericTraits< T >::ValueType >, U > DisableIfMultiComponent
virtual void ThreadedApplyUpdate(const InputImageRegionType &regionToProcess, const int)
static ITK_THREAD_RETURN_FUNCTION_CALL_CONVENTION ComputeSigmaUpdateThreaderCallback(void *arg)
void SetComponent(T &pix, unsigned int idx, typename EnableIfMultiComponent< T, RealValueType >::type val) const
EnableIfMultiComponent< T, typenameNumericTraits< T >::ValueType >::type GetComponent(const T &pix, unsigned int idx) const
typename NumericTraits< PixelValueType >::RealType RealValueType
EnableIfMultiComponent< typenameTImageType::PixelType >::type ComputeMinMax(const TImageType *img)
void SetNoiseSigma(const RealType &sigma)
typename NumericTraits< PixelType >::RealType RealType
std::vector< EigenValuesArrayType > EigenValuesCacheType
typename BaseSamplerType::InstanceIdentifier InstanceIdentifier
void SetComponent(T &pix, unsigned int, typename DisableIfMultiComponent< T, RealValueType >::type val) const
A method to generically set a component.
virtual RealArrayType ResolveSigmaUpdate()
void DispatchedVectorMinMax(const TInputImageType *img)
void DispatchedRiemannianMinMax(const TInputImageType *img)
virtual RealType ComputeGradientJointEntropy(InstanceIdentifier id, typename ListAdaptorType::Pointer &inList, BaseSamplerPointer &sampler, ThreadDataStruct &threadData)
void GenerateInputRequestedRegion() override
ThreadDataStruct ThreadedRiemannianMinMax(const InputImageRegionType &regionToProcess, const int, const InputImageType *img, ThreadDataStruct threadData)
RealType AddUpdate(const DiffusionTensor3D< RealValueType > &a, const RealType &b)
The base class for all process objects (source, filters, mappers) in the Insight data processing pipe...
This an abstract subsampler that constrains subsamples to be contained within a given image region.
typename TSample::InstanceIdentifier InstanceIdentifier
SmartPointer< Self > Pointer
The "itk" namespace contains all Insight Segmentation and Registration Toolkit (ITK) classes....
unsigned long SizeValueType
Definition: itkIntTypes.h:86
itk::ITK_THREAD_RETURN_TYPE ITK_THREAD_RETURN_FUNCTION_CALL_CONVENTION