ITK  4.4.0
Insight Segmentation and Registration Toolkit
itkPatchBasedDenoisingImageFilter.h
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright Insight Software Consortium
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  * http://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 
21 #include "itkPatchBasedDenoisingBaseImageFilter.h"
22 #include "itkImageRegionIterator.h"
23 #include "itkNeighborhoodAlgorithm.h"
24 #include "itkVector.h"
25 #include "itkVectorImage.h"
26 #include "itkRGBPixel.h"
27 #include "itkRGBAPixel.h"
28 #include "itkDiffusionTensor3D.h"
29 #include "itkFixedArray.h"
30 #include "itkMatrix.h"
31 #include "itkRegionConstrainedSubsampler.h"
32 #include "itkEnableIf.h"
33 #include "itkIsSame.h"
34 
35 #include <vector>
36 
37 namespace itk
38 {
39 
62 template <class TInputImage, class TOutputImage>
63 class ITK_EXPORT PatchBasedDenoisingImageFilter :
64  public PatchBasedDenoisingBaseImageFilter<TInputImage, TOutputImage>
65 {
66 public:
68  typedef PatchBasedDenoisingImageFilter Self;
69  typedef PatchBasedDenoisingBaseImageFilter<TInputImage, TOutputImage> Superclass;
70  typedef SmartPointer<Self> Pointer;
71  typedef SmartPointer<const Self> ConstPointer;
72  typedef typename Superclass::OutputImagePointer OutputImagePointer;
73 
75  itkNewMacro(Self);
76 
78  itkTypeMacro(PatchBasedDenoisingImageFilter, ImageToImageFilter);
79 
81  typedef typename Superclass::InputImageType InputImageType;
82  typedef typename Superclass::OutputImageType OutputImageType;
83 
85  itkStaticConstMacro(ImageDimension, unsigned int,
86  Superclass::ImageDimension);
87 
89  typedef typename InputImageType::RegionType InputImageRegionType;
90 
92  typedef ImageRegionIterator<OutputImageType> OutputImageRegionIteratorType;
93  typedef ImageRegionConstIterator<InputImageType> InputImageRegionConstIteratorType;
94 
97  typedef typename Superclass::PixelType PixelType;
98  typedef typename Superclass::PixelValueType PixelValueType;
99 
100  typedef typename NumericTraits< PixelType >::RealType RealType;
101  typedef typename NumericTraits< PixelValueType >::RealType RealValueType;
102  typedef Array<PixelValueType> PixelArrayType;
103  typedef Array<RealValueType> RealArrayType;
104  typedef Array<unsigned short> ShortArrayType;
105 
107  typedef typename Superclass::ListAdaptorType ListAdaptorType;
108  typedef typename Superclass::PatchRadiusType PatchRadiusType;
109  typedef typename Superclass::InputImagePatchIterator InputImagePatchIterator;
110  typedef ListAdaptorType PatchSampleType;
111  typedef typename Superclass::PatchWeightsType PatchWeightsType;
112 
114  typedef itk::Statistics::RegionConstrainedSubsampler<
115  PatchSampleType, InputImageRegionType > BaseSamplerType;
116  typedef typename BaseSamplerType::Pointer BaseSamplerPointer;
117  typedef typename BaseSamplerType::InstanceIdentifier InstanceIdentifier;
118 
124  typedef FixedArray< PixelValueType, 3 > EigenValuesArrayType;
125  typedef Matrix< PixelValueType, 3, 3 > EigenVectorsMatrixType;
126  typedef std::vector<EigenValuesArrayType*> EigenValuesCacheType;
127  typedef std::vector<EigenVectorsMatrixType*> EigenVectorsCacheType;
128 
129  struct ThreadDataStruct
130  {
131  ShortArrayType validDerivatives;
132  RealArrayType entropyFirstDerivative;
133  RealArrayType entropySecondDerivative;
134  ShortArrayType validNorms;
135  RealArrayType minNorm;
136  RealArrayType maxNorm;
137  BaseSamplerPointer sampler;
138  EigenValuesCacheType eigenValsCache;
139  EigenVectorsCacheType eigenVecsCache;
140  };
141 
146  itkSetMacro(UseSmoothDiscPatchWeights, bool);
147  itkBooleanMacro(UseSmoothDiscPatchWeights);
148  itkGetConstMacro(UseSmoothDiscPatchWeights, bool);
150 
155  void SetKernelBandwidthSigma(const RealArrayType& kernelSigma);
156  itkGetConstMacro(KernelBandwidthSigma, RealArrayType);
158 
163  itkSetClampMacro(KernelBandwidthFractionPixelsForEstimation, double, 0.01, 1.0);
164  itkGetConstReferenceMacro(KernelBandwidthFractionPixelsForEstimation, double);
166 
169  itkSetMacro(ComputeConditionalDerivatives, bool);
170  itkBooleanMacro(ComputeConditionalDerivatives);
171  itkGetConstMacro(ComputeConditionalDerivatives, bool);
173 
187  itkSetMacro(UseFastTensorComputations, bool);
188  itkBooleanMacro(UseFastTensorComputations);
189  itkGetConstMacro(UseFastTensorComputations, bool);
191 
193  itkStaticConstMacro(MaxSigmaUpdateIterations, unsigned int,
194  20);
195 
202  itkSetClampMacro(KernelBandwidthMultiplicationFactor, double, 0.01, 100);
203  itkGetConstReferenceMacro(KernelBandwidthMultiplicationFactor, double);
205 
209  void SetNoiseSigma(const RealType& sigma);
210 
211  itkGetConstMacro(NoiseSigma, RealType);
212 
214  itkSetObjectMacro(Sampler, BaseSamplerType);
215  itkGetModifiableObjectMacro(Sampler, BaseSamplerType);
217 
218 protected:
219  PatchBasedDenoisingImageFilter();
220  ~PatchBasedDenoisingImageFilter();
221  virtual void PrintSelf(std::ostream& os, Indent indent) const;
222 
224  virtual void EmptyCaches();
225 
227  virtual void AllocateUpdateBuffer();
228 
229  virtual void CopyInputToOutput();
230 
231  virtual void GenerateInputRequestedRegion();
232 
239  template< typename T>
240  typename EnableIfC<
241  IsSame<T, typename NumericTraits<T>::ValueType>::Value,
242  T >::Type
243  GetComponent(const T pix,
244  unsigned int itkNotUsed( idx ) ) const
245  {
246  // The enable if idiom is used to overload this method for both
247  // scalars and multi-component types. By exploiting that
248  // NumericTraits' ValueType typedef (defines the per-element type
249  // for multi-component types ) is different then the parameterize
250  // type, the bracket operator is used only for multi-component
251  // types.
252  return pix;
253  }
254 
255  template< typename T>
256  typename DisableIfC<
257  IsSame<T, typename NumericTraits<T>::ValueType>::Value,
258  typename NumericTraits<T>::ValueType>::Type
259  GetComponent(const T& pix,
260  unsigned int idx ) const
261  {
262  return pix[idx];
263  }
264 
266  template< typename T >
267  void
268  SetComponent( T &pix,
269  unsigned int itkNotUsed( idx ),
270  typename EnableIfC< IsSame<T,
271  typename NumericTraits<T>::ValueType>::Value, RealValueType>::Type val) const
272  {
273  pix = val;
274  }
275 
276  template< typename T >
277  void
278  SetComponent( T &pix,
279  unsigned int idx,
280  typename DisableIfC< IsSame<T,
281  typename NumericTraits<T>::ValueType>::Value,
282  RealValueType>::Type val) const
283  {
284  pix[idx] = val;
285  }
286 
290  void ComputeMinMax(const Image< DiffusionTensor3D<PixelValueType> , ImageDimension>* img)
291  {
292  if (this->m_ComponentSpace == Superclass::RIEMANNIAN)
293  {
294  DispatchedRiemannianMinMax(img);
295  }
296  else
297  {
298  DispatchedArrayMinMax(img);
299  }
300  }
301 
302  template< typename TImageType>
303  typename EnableIfC<
304  IsSame<typename TImageType::PixelType, typename NumericTraits<typename TImageType::PixelType>::ValueType>::Value
305  >::Type
306  ComputeMinMax(const TImageType* img)
307  {
308  DispatchedMinMax(img);
309  }
310 
311  template< typename TImageType>
312  typename DisableIfC<
313  IsSame<typename TImageType::PixelType, typename NumericTraits<typename TImageType::PixelType>::ValueType>::Value
314  >::Type
315  ComputeMinMax(const TImageType* img)
316  {
317  DispatchedArrayMinMax(img);
318  }
319 
328  void ComputeDifferenceAndWeightedSquaredNorm(const DiffusionTensor3D<PixelValueType>& a,
329  const DiffusionTensor3D<PixelValueType>& b,
330  const RealArrayType& weight,
331  bool useCachedComputations,
332  SizeValueType cacheIndex,
333  EigenValuesCacheType& eigenValsCache,
334  EigenVectorsCacheType& eigenVecsCache,
335  RealType& diff, RealArrayType& norm)
336  {
337  if (this->m_ComponentSpace == Superclass::RIEMANNIAN)
338  {
339  ComputeLogMapAndWeightedSquaredGeodesicDifference(a, b, weight,
340  useCachedComputations, cacheIndex,
341  eigenValsCache, eigenVecsCache,
342  diff, norm);
343  }
344  else
345  {
346  ComputeSignedEuclideanDifferenceAndWeightedSquaredNorm(a, b, weight,
347  useCachedComputations, cacheIndex,
348  eigenValsCache, eigenVecsCache,
349  diff, norm);
350  }
351  }
353 
354  template <class PixelT>
355  void ComputeDifferenceAndWeightedSquaredNorm(const PixelT& a,
356  const PixelT& b,
357  const RealArrayType& weight,
358  bool useCachedComputations,
359  SizeValueType cacheIndex,
360  EigenValuesCacheType& eigenValsCache,
361  EigenVectorsCacheType& eigenVecsCache,
362  RealType& diff, RealArrayType& norm)
363  {
364  ComputeSignedEuclideanDifferenceAndWeightedSquaredNorm(a, b, weight,
365  useCachedComputations, cacheIndex,
366  eigenValsCache, eigenVecsCache,
367  diff, norm);
368  }
369 
373  RealType AddUpdate(const DiffusionTensor3D<RealValueType>& a,
374  const RealType& b)
375  {
376  if (this->m_ComponentSpace == Superclass::RIEMANNIAN)
377  {
378  return this->AddExponentialMapUpdate(a, b);
379  }
380  else
381  {
382  return this->AddEuclideanUpdate(a, b);
383  }
384  }
386 
387  template <class RealT>
388  RealType AddUpdate(const RealT& a,
389  const RealType& b)
390  {
391  return this->AddEuclideanUpdate(a, b);
392  }
393 
394  virtual void EnforceConstraints();
395 
396  virtual void Initialize();
397 
398  virtual void InitializeKernelSigma();
399 
400  virtual void InitializePatchWeights();
401 
402  virtual void InitializePatchWeightsSmoothDisc();
403 
404  virtual void InitializeIteration();
405 
406  virtual void ComputeKernelBandwidthUpdate(); // derived from base class;
407 
408  // define here
409 
410  virtual ThreadDataStruct ThreadedComputeSigmaUpdate(const InputImageRegionType& regionToProcess,
411  const int itkNotUsed(threadId),
412  ThreadDataStruct threadData);
413 
414  virtual RealArrayType ResolveSigmaUpdate();
415 
416  virtual void ComputeImageUpdate();
417 
418  virtual ThreadDataStruct ThreadedComputeImageUpdate(const InputImageRegionType& regionToProcess,
419  const int threadId,
420  ThreadDataStruct threadData);
421 
422  virtual RealType ComputeGradientJointEntropy(InstanceIdentifier id,
423  typename ListAdaptorType::Pointer& inList,
424  BaseSamplerPointer& sampler,
425  ThreadDataStruct& threadData);
426 
427  virtual void ApplyUpdate();
428 
429  virtual void ThreadedApplyUpdate(const InputImageRegionType& regionToProcess,
430  const int itkNotUsed(threadId) );
431 
432  virtual void PostProcessOutput();
433 
434  virtual void SetThreadData(int threadId, const ThreadDataStruct& data);
435 
436  virtual ThreadDataStruct GetThreadData(int threadId);
437 
438  std::vector<ThreadDataStruct> m_ThreadData;
439 
441  typename OutputImageType::Pointer m_UpdateBuffer;
442 
443  unsigned int m_NumPixelComponents;
444  unsigned int m_NumIndependentComponents;
445  unsigned int m_TotalNumberPixels;
446  //
447  bool m_UseSmoothDiscPatchWeights;
448  //
449  bool m_UseFastTensorComputations;
450  //
451  RealArrayType m_KernelBandwidthSigma;
452  bool m_KernelBandwidthSigmaIsSet;
453  RealArrayType m_IntensityRescaleInvFactor;
454  PixelType m_ZeroPixel;
455  PixelArrayType m_ImageMin;
456  PixelArrayType m_ImageMax;
457  double m_KernelBandwidthFractionPixelsForEstimation;
458  bool m_ComputeConditionalDerivatives;
459  double m_MinSigma;
460  double m_MinProbability;
461  unsigned int m_SigmaUpdateDecimationFactor;
462  double m_SigmaUpdateConvergenceTolerance;
463  ShortArrayType m_SigmaConverged;
464  double m_KernelBandwidthMultiplicationFactor;
465  //
466  RealType m_NoiseSigma;
467  RealType m_NoiseSigmaSquared;
468  bool m_NoiseSigmaIsSet;
469  //
470  BaseSamplerPointer m_Sampler;
471  typename ListAdaptorType::Pointer m_SearchSpaceList;
472 
473 private:
474  PatchBasedDenoisingImageFilter(const Self&); // purposely not implemented
475  void operator=(const Self&); // purposely not implemented
476 
479  static ITK_THREAD_RETURN_TYPE ComputeSigmaUpdateThreaderCallback( void *arg );
480 
483  static ITK_THREAD_RETURN_TYPE ComputeImageUpdateThreaderCallback( void *arg );
484 
487  static ITK_THREAD_RETURN_TYPE ApplyUpdateThreaderCallback( void *arg );
488 
489  template <typename TInputImageType>
490  void DispatchedMinMax(const TInputImageType* img);
491 
492  template <typename TInputImageType>
493  void DispatchedArrayMinMax(const TInputImageType* img);
494 
495  template <typename TInputImageType>
496  void DispatchedVectorMinMax(const TInputImageType* img);
497 
498  template <typename TInputImageType>
499  void DispatchedRiemannianMinMax(const TInputImageType* img);
500 
503  static ITK_THREAD_RETURN_TYPE RiemannianMinMaxThreaderCallback( void *arg );
504 
505  ThreadDataStruct ThreadedRiemannianMinMax(const InputImageRegionType& regionToProcess,
506  const int itkNotUsed(threadId),
507  const InputImageType* img,
508  ThreadDataStruct threadData);
509 
510  virtual void ResolveRiemannianMinMax();
511 
512  void ComputeSignedEuclideanDifferenceAndWeightedSquaredNorm(const PixelType& a, const PixelType& b,
513  const RealArrayType& weight,
514  bool useCachedComputations,
515  SizeValueType cacheIndex,
516  EigenValuesCacheType& eigenValsCache,
517  EigenVectorsCacheType& eigenVecsCache,
518  RealType& diff, RealArrayType& norm);
519 
521  void ComputeLogMapAndWeightedSquaredGeodesicDifference(const DiffusionTensor3D<PixelValueType>& spdMatrixA,
522  const DiffusionTensor3D<PixelValueType>& spdMatrixB,
523  const RealArrayType& weight,
524  bool useCachedComputations,
525  SizeValueType cacheIndex,
526  EigenValuesCacheType& eigenValsCache,
527  EigenVectorsCacheType& eigenVecsCache,
528  RealType& symMatrixLogMap, RealArrayType& geodesicDist);
529 
530  template <typename TensorValueT>
531  void Compute3x3EigenAnalysis(const DiffusionTensor3D<TensorValueT>& spdMatrix,
532  FixedArray< TensorValueT, 3 >& eigenVals,
533  Matrix< TensorValueT, 3, 3 >& eigenVecs);
534 
535  RealType AddEuclideanUpdate(const RealType& a, const RealType& b);
536 
538  RealType AddExponentialMapUpdate(const DiffusionTensor3D<RealValueType>& spdMatrix,
539  const DiffusionTensor3D<RealValueType>& symMatrix);
540 
541  struct ThreadFilterStruct
542  {
543  PatchBasedDenoisingImageFilter *Filter;
544  InputImageType *Img;
545  };
546 
547 }; // end class PatchBasedDenoisingImageFilter
548 
549 } // end namespace itk
550 
551 #ifndef ITK_MANUAL_INSTANTIATION
552 # include "itkPatchBasedDenoisingImageFilter.hxx"
553 #endif
554 
555 #endif
556 
Tarballs of the nightly generated Doxygen documentation are available for the html, xml, and tag file.
Generated on Mon May 13 2013 01:03:58 for ITK by   doxygen 1.8.3.1