Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
itkNonLocalMeansDenoisingFilter.txx
Go to the documentation of this file.
1 /*===================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center,
6 Division of Medical and Biological Informatics.
7 All rights reserved.
8 
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
11 A PARTICULAR PURPOSE.
12 
13 See LICENSE.txt or http://www.mitk.org for details.
14 
15 ===================================================================*/
16 
17 #ifndef __itkNonLocalMeansDenoisingFilter_txx
18 #define __itkNonLocalMeansDenoisingFilter_txx
19 
20 #include <time.h>
21 #include <stdio.h>
22 #include <stdlib.h>
23 
24 #define _USE_MATH_DEFINES
25 #include <math.h>
26 
27 #include "itkImageRegionIterator.h"
28 #include "itkNeighborhoodIterator.h"
29 #include <itkImageRegionIteratorWithIndex.h>
30 #include <vector>
31 
32 namespace itk {
33 
34 
35 template< class TPixelType >
36 NonLocalMeansDenoisingFilter< TPixelType >
37 ::NonLocalMeansDenoisingFilter()
38  : m_SearchRadius(4),
39  m_ComparisonRadius(1),
40  m_UseJointInformation(false),
41  m_UseRicianAdaption(false),
42  m_Variance(1),
43  m_Mask(NULL)
44 {
45  this->SetNumberOfRequiredInputs( 1 );
46 }
47 
48 template< class TPixelType >
49 void
50 NonLocalMeansDenoisingFilter< TPixelType >
51 ::BeforeThreadedGenerateData()
52 {
53 
54  MITK_INFO << "SearchRadius: " << m_SearchRadius;
55  MITK_INFO << "ComparisonRadius: " << m_ComparisonRadius;
56  MITK_INFO << "Noisevariance: " << m_Variance;
57  MITK_INFO << "Use Rician Adaption: " << std::boolalpha << m_UseRicianAdaption;
58  MITK_INFO << "Use Joint Information: " << std::boolalpha << m_UseJointInformation;
59 
60 
61  typename InputImageType::Pointer inputImagePointer = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) );
62  if (m_Mask.IsNull())
63  {
64  // If no mask is used generate a mask of the complete image
65 
66  m_Mask = MaskImageType::New();
67  m_Mask->SetRegions(inputImagePointer->GetLargestPossibleRegion());
68  m_Mask->Allocate();
69  m_Mask->FillBuffer(1);
70  }
71  else
72  {
73  // Calculation of the smallest masked region
74 
75  typename OutputImageType::Pointer outputImage =
76  static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0));
77  ImageRegionIterator< OutputImageType > oit(outputImage, inputImagePointer->GetLargestPossibleRegion());
78  oit.GoToBegin();
79  ImageRegionIterator< MaskImageType > mit(m_Mask, m_Mask->GetLargestPossibleRegion());
80  mit.GoToBegin();
81  typename MaskImageType::IndexType minIndex;
82  typename MaskImageType::IndexType maxIndex;
83  minIndex.Fill(10000);
84  maxIndex.Fill(0);
85  typename OutputImageType::PixelType outpix;
86  outpix.SetSize(inputImagePointer->GetVectorLength());
87  outpix.Fill(0);
88  while (!mit.IsAtEnd())
89  {
90 
91  if (mit.Get())
92  {
93  // calculation of the start & end index of the smallest masked region
94  minIndex[0] = minIndex[0] < mit.GetIndex()[0] ? minIndex[0] : mit.GetIndex()[0];
95  minIndex[1] = minIndex[1] < mit.GetIndex()[1] ? minIndex[1] : mit.GetIndex()[1];
96  minIndex[2] = minIndex[2] < mit.GetIndex()[2] ? minIndex[2] : mit.GetIndex()[2];
97 
98  maxIndex[0] = maxIndex[0] > mit.GetIndex()[0] ? maxIndex[0] : mit.GetIndex()[0];
99  maxIndex[1] = maxIndex[1] > mit.GetIndex()[1] ? maxIndex[1] : mit.GetIndex()[1];
100  maxIndex[2] = maxIndex[2] > mit.GetIndex()[2] ? maxIndex[2] : mit.GetIndex()[2];
101  }
102  else
103  {
104  oit.Set(outpix);
105  }
106  ++mit;
107  ++oit;
108  }
109 
110 
111  // calculation of the masked region
112  typename OutputImageType::SizeType size;
113  size[0] = maxIndex[0] - minIndex[0] + 1;
114  size[1] = maxIndex[1] - minIndex[1] + 1;
115  size[2] = maxIndex[2] - minIndex[2] + 1;
116 
117  typename OutputImageType::RegionType region (minIndex, size);
118 
119  outputImage->SetRequestedRegion(region);
120  }
121 
122  m_CurrentVoxelCount = 0;
123 }
124 
125 template< class TPixelType >
126 void
127 NonLocalMeansDenoisingFilter< TPixelType >
128 ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, ThreadIdType )
129 {
130 
131 
132  // initialize iterators
133  typename OutputImageType::Pointer outputImage =
134  static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0));
135 
136  ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread);
137  oit.GoToBegin();
138 
139  ImageRegionIterator< MaskImageType > mit(m_Mask, outputRegionForThread);
140  mit.GoToBegin();
141 
142 
143 
144  typedef ImageRegionIteratorWithIndex <InputImageType> InputIteratorType;
145  typename InputImageType::Pointer inputImagePointer = NULL;
146  inputImagePointer = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) );
147 
148  InputIteratorType git(inputImagePointer, outputRegionForThread );
149  git.GoToBegin();
150 
151  // iterate over complete image region
152  while( !git.IsAtEnd() )
153  {
154  typename OutputImageType::PixelType outpix;
155  outpix.SetSize (inputImagePointer->GetVectorLength());
156 
157  if (mit.Get() != 0 && !this->GetAbortGenerateData())
158  {
159  if(!m_UseJointInformation)
160  {
161  for (int i = 0; i < (int)inputImagePointer->GetVectorLength(); ++i)
162  {
163  double summw = 0;
164  double sumj = 0;
165  double w = 0;
166  std::vector<double> wj;
167  std::vector<double> p;
168  typename InputIteratorType::IndexType index;
169  index = git.GetIndex();
170 
171  for (int x = index.GetElement(0) - m_SearchRadius; x <= index.GetElement(0) + m_SearchRadius; ++x)
172  {
173  for (int y = index.GetElement(1) - m_SearchRadius; y <= index.GetElement(1) + m_SearchRadius; ++y)
174  {
175  for (int z = index.GetElement(2) - m_SearchRadius; z <= index.GetElement(2) + m_SearchRadius; ++z)
176  {
177  typename InputIteratorType::IndexType indexV;
178  indexV.SetElement(0, x);
179  indexV.SetElement(1, y);
180  indexV.SetElement(2, z);
181  if (inputImagePointer->GetLargestPossibleRegion().IsInside(indexV))
182  {
183  TPixelType pixelJ = inputImagePointer->GetPixel(indexV)[i];
184  double sumk = 0;
185  double size = 0;
186  for (int xi = index.GetElement(0) - m_ComparisonRadius, xj = x - m_ComparisonRadius; xi <= index.GetElement(0) + m_ComparisonRadius; ++xi, ++xj)
187  {
188  for (int yi = index.GetElement(1) - m_ComparisonRadius, yj = y - m_ComparisonRadius; yi <= index.GetElement(1) + m_ComparisonRadius; ++yi, ++yj)
189  {
190  for (int zi = index.GetElement(2) - m_ComparisonRadius, zj = z - m_ComparisonRadius; zi <= index.GetElement(2) + m_ComparisonRadius; ++zi, ++zj)
191  {
192  typename InputIteratorType::IndexType indexI, indexJ;
193  indexI.SetElement(0, xi);
194  indexI.SetElement(1, yi);
195  indexI.SetElement(2, zi);
196  indexJ.SetElement(0, xj);
197  indexJ.SetElement(1, yj);
198  indexJ.SetElement(2, zj);
199 
200 
201  // Compare neighborhoods ni & nj
202  if (inputImagePointer->GetLargestPossibleRegion().IsInside(indexI) && inputImagePointer->GetLargestPossibleRegion().IsInside(indexJ))
203  {
204  int diff = inputImagePointer->GetPixel(indexI)[i] - inputImagePointer->GetPixel(indexJ)[i];
205  sumk += (double)(diff*diff);
206  ++size;
207  }
208  }
209  }
210  }
211  // weight all neighborhoods
212  w = std::exp( - sumk / size / m_Variance);
213  wj.push_back(w);
214  if (m_UseRicianAdaption)
215  {
216  p.push_back((double)(pixelJ*pixelJ));
217  }
218  else
219  {
220  p.push_back((double)(pixelJ));
221  }
222  summw += w;
223  }
224  }
225  }
226  }
227  for (unsigned int n = 0; n < wj.size(); ++n)
228  {
229  sumj += (wj[n]/summw) * p[n];
230  }
231  if (m_UseRicianAdaption)
232  {
233  sumj -=2 * m_Variance;
234  }
235 
236  if (sumj < 0)
237  {
238  sumj = 0;
239  }
240 
241  TPixelType outval;
242  if (m_UseRicianAdaption)
243  {
244  outval = std::floor(std::sqrt(sumj) + 0.5);
245  }
246  else
247  {
248  outval = std::floor(sumj + 0.5);
249  }
250  outpix.SetElement(i, outval);
251  }
252  }
253 
254  else
255  {
256  // same procedure for vektoranalysis
257 
258  double Z = 0;
259  itk::VariableLengthVector<double> sumj;
260  sumj.SetSize(inputImagePointer->GetVectorLength());
261  sumj.Fill(0.0);
262  double w = 0;
263  std::vector<double> wj;
264  std::vector<itk::VariableLengthVector <double> > p;
265  typename InputIteratorType::IndexType index;
266  index = git.GetIndex();
267 
268  for (int x = index.GetElement(0) - m_SearchRadius; x <= index.GetElement(0) + m_SearchRadius; ++x)
269  {
270  for (int y = index.GetElement(1) - m_SearchRadius; y <= index.GetElement(1) + m_SearchRadius; ++y)
271  {
272  for (int z = index.GetElement(2) - m_SearchRadius; z <= index.GetElement(2) + m_SearchRadius; ++z)
273  {
274  typename InputIteratorType::IndexType indexV;
275  indexV.SetElement(0, x);
276  indexV.SetElement(1, y);
277  indexV.SetElement(2, z);
278  if (inputImagePointer->GetLargestPossibleRegion().IsInside(indexV))
279  {
280  typename InputImageType::PixelType pixelJ = inputImagePointer->GetPixel(indexV);
281  double sumk = 0;
282  double size = 0;
283  for (int xi = index.GetElement(0) - m_ComparisonRadius, xj = x - m_ComparisonRadius; xi <= index.GetElement(0) + m_ComparisonRadius; ++xi, ++xj)
284  {
285  for (int yi = index.GetElement(1) - m_ComparisonRadius, yj = y - m_ComparisonRadius; yi <= index.GetElement(1) + m_ComparisonRadius; ++yi, ++yj)
286  {
287  for (int zi = index.GetElement(2) - m_ComparisonRadius, zj = z - m_ComparisonRadius; zi <= index.GetElement(2) + m_ComparisonRadius; ++zi, ++zj)
288  {
289  typename InputIteratorType::IndexType indexI, indexJ;
290  indexI.SetElement(0, xi);
291  indexI.SetElement(1, yi);
292  indexI.SetElement(2, zi);
293  indexJ.SetElement(0, xj);
294  indexJ.SetElement(1, yj);
295  indexJ.SetElement(2, zj);
296  // Compare neighborhoods ni & nj
297  if (inputImagePointer->GetLargestPossibleRegion().IsInside(indexI) && inputImagePointer->GetLargestPossibleRegion().IsInside(indexJ))
298  {
299  typename InputImageType::PixelType diff = inputImagePointer->GetPixel(indexI) - inputImagePointer->GetPixel(indexJ);
300  sumk += (double)(diff.GetSquaredNorm());
301  ++size;
302  }
303  }
304  }
305  }
306  // weight all neighborhoods
307  size *= inputImagePointer->GetVectorLength() + 1;
308  w = std::exp( - (sumk / size) / m_Variance);
309  wj.push_back(w);
310  if (m_UseRicianAdaption)
311  {
312  itk::VariableLengthVector <double> m;
313  m.SetSize(inputImagePointer->GetVectorLength());
314  for (unsigned int i = 0; i < inputImagePointer->GetVectorLength(); ++i)
315  {
316  double sp = (double)(pixelJ.GetElement(i) * pixelJ.GetElement(i));
317  m.SetElement(i,sp);
318  }
319  p.push_back(m);
320  }
321  else
322  ++size;
323  {
324  p.push_back(pixelJ);
325  }
326  Z += w;
327  }
328  }
329  }
330  }
331  for (unsigned int n = 0; n < wj.size(); ++n)
332  {
333  sumj = sumj + ((wj[n]/Z) * p[n]);
334  }
335  if (m_UseRicianAdaption)
336  {
337  sumj = sumj - (2 * m_Variance);
338  }
339 
340  for (unsigned int i = 0; i < inputImagePointer->GetVectorLength(); ++i)
341  {
342  double a = sumj.GetElement(i);
343  if (a < 0)
344  {
345  a = 0;
346  }
347  TPixelType outval;
348  if (m_UseRicianAdaption)
349  {
350  outval = std::floor(std::sqrt(a) + 0.5);
351  }
352  else
353  {
354  outval = std::floor(a + 0.5);
355  }
356  outpix.SetElement(i, outval);
357  }
358  }
359  }
360  else
361  {
362  outpix.Fill(0);
363  }
364 
365  oit.Set(outpix);
366  ++oit;
367  ++m_CurrentVoxelCount;
368  ++git;
369  ++mit;
370  }
371 
372  MITK_INFO << "One Thread finished calculation";
373 }
374 
375 template< class TPixelType >
376 void NonLocalMeansDenoisingFilter< TPixelType >::SetInputImage(const InputImageType* image)
377 {
378  this->SetNthInput(0, const_cast< InputImageType* >(image));
379 }
380 
381 template< class TPixelType >
382 void NonLocalMeansDenoisingFilter< TPixelType >::SetInputMask(MaskImageType* mask)
383 {
384  m_Mask = mask;
385 }
386 }
387 
388 #endif // __itkNonLocalMeansDenoisingFilter_txx