Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
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