ITK  4.3.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  itkSetMacro(GaussianKernelSigma, RealArrayType);
156  itkGetConstMacro(GaussianKernelSigma, RealArrayType);
158 
163  itkSetClampMacro(FractionPixelsForSigmaUpdate, double, 0.01, 1.0);
164  itkGetConstReferenceMacro(FractionPixelsForSigmaUpdate, 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 
201  itkSetClampMacro(SigmaMultiplicationFactor, double, 0.01, 100);
202  itkGetConstReferenceMacro(SigmaMultiplicationFactor, double);
204 
208  void SetNoiseSigma(const RealType& sigma);
209 
210  itkGetConstMacro(NoiseSigma, RealType);
211 
213  itkSetObjectMacro(Sampler, BaseSamplerType);
214  itkGetObjectMacro(Sampler, BaseSamplerType);
216 
217 protected:
218  PatchBasedDenoisingImageFilter();
219  ~PatchBasedDenoisingImageFilter();
220  virtual void PrintSelf(std::ostream& os, Indent indent) const;
221 
223  virtual void EmptyCaches();
224 
226  virtual void AllocateUpdateBuffer();
227 
228  virtual void CopyInputToOutput();
229 
230  virtual void GenerateInputRequestedRegion();
231 
238  template< typename T>
239  typename EnableIfC<
240  IsSame<T, typename NumericTraits<T>::ValueType>::Value,
241  T >::Type
242  GetComponent(const T pix,
243  unsigned int itkNotUsed( idx ) ) const
244  {
245  // The enable if idiom is used to overload this method for both
246  // scalars and multi-component types. By exploiting that
247  // NumericTraits' ValueType typedef (defines the per-element type
248  // for multi-component types ) is different then the parameterize
249  // type, the bracket operator is used only for multi-component
250  // types.
251  return pix;
252  }
253 
254  template< typename T>
255  typename DisableIfC<
256  IsSame<T, typename NumericTraits<T>::ValueType>::Value,
257  typename NumericTraits<T>::ValueType>::Type
258  GetComponent(const T& pix,
259  unsigned int idx ) const
260  {
261  return pix[idx];
262  }
263 
265  template< typename T >
266  void
267  SetComponent( T &pix,
268  unsigned int itkNotUsed( idx ),
269  typename EnableIfC< IsSame<T,
270  typename NumericTraits<T>::ValueType>::Value, RealValueType>::Type val) const
271  {
272  pix = val;
273  }
274 
275  template< typename T >
276  void
277  SetComponent( T &pix,
278  unsigned int idx,
279  typename DisableIfC< IsSame<T,
280  typename NumericTraits<T>::ValueType>::Value,
281  RealValueType>::Type val) const
282  {
283  pix[idx] = val;
284  }
285 
289  void ComputeMinMax(const Image< DiffusionTensor3D<PixelValueType> , ImageDimension>* img)
290  {
291  if (this->m_ComponentSpace == Superclass::RIEMANNIAN)
292  {
293  DispatchedRiemannianMinMax(img);
294  }
295  else
296  {
297  DispatchedArrayMinMax(img);
298  }
299  }
300 
301  template< typename TImageType>
302  typename EnableIfC<
303  IsSame<typename TImageType::PixelType, typename NumericTraits<typename TImageType::PixelType>::ValueType>::Value
304  >::Type
305  ComputeMinMax(const TImageType* img)
306  {
307  DispatchedMinMax(img);
308  }
309 
310  template< typename TImageType>
311  typename DisableIfC<
312  IsSame<typename TImageType::PixelType, typename NumericTraits<typename TImageType::PixelType>::ValueType>::Value
313  >::Type
314  ComputeMinMax(const TImageType* img)
315  {
316  DispatchedArrayMinMax(img);
317  }
318 
327  void ComputeDifferenceAndWeightedSquaredNorm(const DiffusionTensor3D<PixelValueType>& a,
328  const DiffusionTensor3D<PixelValueType>& b,
329  const RealArrayType& weight,
330  bool useCachedComputations,
331  SizeValueType cacheIndex,
332  EigenValuesCacheType& eigenValsCache,
333  EigenVectorsCacheType& eigenVecsCache,
334  RealType& diff, RealArrayType& norm)
335  {
336  if (this->m_ComponentSpace == Superclass::RIEMANNIAN)
337  {
338  ComputeLogMapAndWeightedSquaredGeodesicDifference(a, b, weight,
339  useCachedComputations, cacheIndex,
340  eigenValsCache, eigenVecsCache,
341  diff, norm);
342  }
343  else
344  {
345  ComputeSignedEuclideanDifferenceAndWeightedSquaredNorm(a, b, weight,
346  useCachedComputations, cacheIndex,
347  eigenValsCache, eigenVecsCache,
348  diff, norm);
349  }
350  }
352 
353  template <class PixelT>
354  void ComputeDifferenceAndWeightedSquaredNorm(const PixelT& a,
355  const PixelT& b,
356  const RealArrayType& weight,
357  bool useCachedComputations,
358  SizeValueType cacheIndex,
359  EigenValuesCacheType& eigenValsCache,
360  EigenVectorsCacheType& eigenVecsCache,
361  RealType& diff, RealArrayType& norm)
362  {
363  ComputeSignedEuclideanDifferenceAndWeightedSquaredNorm(a, b, weight,
364  useCachedComputations, cacheIndex,
365  eigenValsCache, eigenVecsCache,
366  diff, norm);
367  }
368 
372  RealType AddUpdate(const DiffusionTensor3D<RealValueType>& a,
373  const RealType& b)
374  {
375  if (this->m_ComponentSpace == Superclass::RIEMANNIAN)
376  {
377  return this->AddExponentialMapUpdate(a, b);
378  }
379  else
380  {
381  return this->AddEuclideanUpdate(a, b);
382  }
383  }
385 
386  template <class RealT>
387  RealType AddUpdate(const RealT& a,
388  const RealType& b)
389  {
390  return this->AddEuclideanUpdate(a, b);
391  }
392 
393  virtual void EnforceConstraints();
394 
395  virtual void Initialize();
396 
397  virtual void InitializeKernelSigma();
398 
399  virtual void InitializePatchWeights();
400 
401  virtual void InitializePatchWeightsSmoothDisc();
402 
403  virtual void InitializeIteration();
404 
405  virtual void ComputeKernelBandwidthUpdate(); // derived from base class;
406 
407  // define here
408 
409  virtual ThreadDataStruct ThreadedComputeSigmaUpdate(const InputImageRegionType& regionToProcess,
410  const int itkNotUsed(threadId),
411  ThreadDataStruct threadData);
412 
413  virtual RealArrayType ResolveSigmaUpdate();
414 
415  virtual void ComputeImageUpdate();
416 
417  virtual ThreadDataStruct ThreadedComputeImageUpdate(const InputImageRegionType& regionToProcess,
418  const int threadId,
419  ThreadDataStruct threadData);
420 
421  virtual RealType ComputeGradientJointEntropy(InstanceIdentifier id,
422  typename ListAdaptorType::Pointer& inList,
423  BaseSamplerPointer& sampler,
424  ThreadDataStruct& threadData);
425 
426  virtual void ApplyUpdate();
427 
428  virtual void ThreadedApplyUpdate(const InputImageRegionType& regionToProcess,
429  const int itkNotUsed(threadId) );
430 
431  virtual void PostProcessOutput();
432 
433  virtual void SetThreadData(int threadId, const ThreadDataStruct& data);
434 
435  virtual ThreadDataStruct GetThreadData(int threadId);
436 
437  std::vector<ThreadDataStruct> m_ThreadData;
438 
440  typename OutputImageType::Pointer m_UpdateBuffer;
441 
442  unsigned int m_NumPixelComponents;
443  unsigned int m_NumIndependentComponents;
444  unsigned int m_TotalNumberPixels;
445  //
446  bool m_UseSmoothDiscPatchWeights;
447  //
448  bool m_UseFastTensorComputations;
449  //
450  RealArrayType m_GaussianKernelSigma;
451  RealArrayType m_IntensityRescaleInvFactor;
452  PixelType m_ZeroPixel;
453  PixelArrayType m_ImageMin;
454  PixelArrayType m_ImageMax;
455  double m_FractionPixelsForSigmaUpdate;
456  bool m_ComputeConditionalDerivatives;
457  double m_MinSigma;
458  double m_MinProbability;
459  unsigned int m_SigmaUpdateDecimationFactor;
460  double m_SigmaUpdateConvergenceTolerance;
461  ShortArrayType m_SigmaConverged;
462  double m_SigmaMultiplicationFactor;
463  //
464  RealType m_NoiseSigma;
465  RealType m_NoiseSigmaSquared;
466  bool m_NoiseSigmaIsSet;
467  //
468  BaseSamplerPointer m_Sampler;
469  typename ListAdaptorType::Pointer m_SearchSpaceList;
470 
471 private:
472  PatchBasedDenoisingImageFilter(const Self&); // purposely not implemented
473  void operator=(const Self&); // purposely not implemented
474 
477  static ITK_THREAD_RETURN_TYPE ComputeSigmaUpdateThreaderCallback( void *arg );
478 
481  static ITK_THREAD_RETURN_TYPE ComputeImageUpdateThreaderCallback( void *arg );
482 
485  static ITK_THREAD_RETURN_TYPE ApplyUpdateThreaderCallback( void *arg );
486 
487  template <typename TInputImageType>
488  void DispatchedMinMax(const TInputImageType* img);
489 
490  template <typename TInputImageType>
491  void DispatchedArrayMinMax(const TInputImageType* img);
492 
493  template <typename TInputImageType>
494  void DispatchedVectorMinMax(const TInputImageType* img);
495 
496  template <typename TInputImageType>
497  void DispatchedRiemannianMinMax(const TInputImageType* img);
498 
501  static ITK_THREAD_RETURN_TYPE RiemannianMinMaxThreaderCallback( void *arg );
502 
503  ThreadDataStruct ThreadedRiemannianMinMax(const InputImageRegionType& regionToProcess,
504  const int itkNotUsed(threadId),
505  const InputImageType* img,
506  ThreadDataStruct threadData);
507 
508  virtual void ResolveRiemannianMinMax();
509 
510  void ComputeSignedEuclideanDifferenceAndWeightedSquaredNorm(const PixelType& a, const PixelType& b,
511  const RealArrayType& weight,
512  bool useCachedComputations,
513  SizeValueType cacheIndex,
514  EigenValuesCacheType& eigenValsCache,
515  EigenVectorsCacheType& eigenVecsCache,
516  RealType& diff, RealArrayType& norm);
517 
519  void ComputeLogMapAndWeightedSquaredGeodesicDifference(const DiffusionTensor3D<PixelValueType>& spdMatrixA,
520  const DiffusionTensor3D<PixelValueType>& spdMatrixB,
521  const RealArrayType& weight,
522  bool useCachedComputations,
523  SizeValueType cacheIndex,
524  EigenValuesCacheType& eigenValsCache,
525  EigenVectorsCacheType& eigenVecsCache,
526  RealType& symMatrixLogMap, RealArrayType& geodesicDist);
527 
528  template <typename TensorValueT>
529  void Compute3x3EigenAnalysis(const DiffusionTensor3D<TensorValueT>& spdMatrix,
530  FixedArray< TensorValueT, 3 >& eigenVals,
531  Matrix< TensorValueT, 3, 3 >& eigenVecs);
532 
533  RealType AddEuclideanUpdate(const RealType& a, const RealType& b);
534 
536  RealType AddExponentialMapUpdate(const DiffusionTensor3D<RealValueType>& spdMatrix,
537  const DiffusionTensor3D<RealValueType>& symMatrix);
538 
539  struct ThreadFilterStruct
540  {
541  PatchBasedDenoisingImageFilter *Filter;
542  InputImageType *Img;
543  };
544 
545 }; // end class PatchBasedDenoisingImageFilter
546 
547 } // end namespace itk
548 
549 #ifndef ITK_MANUAL_INSTANTIATION
550 # include "itkPatchBasedDenoisingImageFilter.hxx"
551 #endif
552 
553 #endif
554 
Tarballs of the nightly generated Doxygen documentation are available for the html, xml, and tag file.
Generated on Sun Dec 9 2012 01:22:14 for ITK by   doxygen 1.8.2