Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
itkBrainMaskExtractionImageFilter.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 /*===================================================================
18 
19 This file is based heavily on a corresponding ITK filter.
20 
21 ===================================================================*/
22 #ifndef __itkBrainMaskExtractionImageFilter_txx
23 #define __itkBrainMaskExtractionImageFilter_txx
24 
25 #include "itkBrainMaskExtractionImageFilter.h"
26 
27 #include <itkBinaryThresholdImageFilter.h>
28 #include <itkBinaryErodeImageFilter.h>
29 #include <itkBinaryDilateImageFilter.h>
30 #include <itkVotingBinaryHoleFillingImageFilter.h>
31 #include <itkVotingBinaryIterativeHoleFillingImageFilter.h>
32 #include <itkBinaryBallStructuringElement.h>
33 #include <itkBinaryCrossStructuringElement.h>
34 #include <itkAndImageFilter.h>
35 #include <itkRecursiveGaussianImageFilter.h>
36 #include <itkMaskImageFilter.h>
37 
38 namespace itk {
39 
40  template< class TOutputImagePixelType >
41  BrainMaskExtractionImageFilter< TOutputImagePixelType >
42  ::BrainMaskExtractionImageFilter()
43  : m_MaxNumIterations(itk::NumericTraits<int>::max())
44  {
45  // At least 1 inputs is necessary for a vector image.
46  // For images added one at a time we need at least six
47  this->SetNumberOfRequiredInputs( 1 );
48  }
49 
50  template< class TOutputImagePixelType >
51  void BrainMaskExtractionImageFilter< TOutputImagePixelType >
52  ::GenerateData()
53  {
54 
55  // gaussian filter
56  typename itk::RecursiveGaussianImageFilter<InputImageType,InputImageType>::Pointer gaussian
57  = itk::RecursiveGaussianImageFilter<InputImageType,InputImageType>::New();
58  gaussian->SetInput( this->GetInput(0) );
59  gaussian->SetSigma( 1.0 );
60 
61  try
62  {
63  gaussian->Update();
64  }
65  catch( itk::ExceptionObject &e)
66  {
67  std::cerr << e;
68  return;
69  }
70 
71  // threshold the image
72  typename itk::BinaryThresholdImageFilter<InputImageType,OutputImageType>::Pointer threshold =
73  itk::BinaryThresholdImageFilter<InputImageType,OutputImageType>::New();
74 
75  threshold->SetInput( gaussian->GetOutput() );
76 
77  int seuil = static_cast<int>( ComputeHistogram( gaussian->GetOutput() ) );
78  threshold->SetLowerThreshold( seuil );
79 
80  std::cout << "Thresholding..." << std::flush;
81  try
82  {
83  threshold->Update();
84  }
85  catch( itk::ExceptionObject &e)
86  {
87  std::cerr << e;
88  return;
89  }
90  std::cout << "Done." << std::endl;
91 
92 #ifdef DEBUG_ME
93  {
94  WriterType::Pointer writer = WriterType::New();
95  writer->SetInput( threshold->GetOutput() );
96  writer->SetFileName( "AfterThreshold.hdr" );
97  writer->Update();
98  }
99 #endif
100 
101  // erode to remove background noise
102  typedef itk::BinaryBallStructuringElement<int, 3> StructuralElementType;
103  StructuralElementType ball;
104 
105  typename itk::BinaryErodeImageFilter<OutputImageType,OutputImageType,StructuralElementType>::Pointer erode =
106  itk::BinaryErodeImageFilter<OutputImageType,OutputImageType,StructuralElementType>::New();
107 
108  ball.SetRadius( 3 );
109 
110  erode->SetInput( threshold->GetOutput() );
111  erode->SetKernel( ball );
112 
113  std::cout << "Eroding..." << std::flush;
114  try
115  {
116  erode->Update();
117  }
118  catch( itk::ExceptionObject &e)
119  {
120  std::cerr << e;
121  return;
122  }
123  std::cout << "Done." << std::endl;
124 
125 #ifdef DEBUG_ME
126  {
127  WriterType::Pointer writer = WriterType::New();
128  writer->SetInput( erode->GetOutput() );
129  writer->SetFileName( "AfterErode.hdr" );
130  writer->Update();
131  }
132 #endif
133 
134 
135  typedef BinaryCrossStructuringElement<int, 3> CrossType;
136 
137  typedef BinaryDilateImageFilter<OutputImageType,OutputImageType,CrossType> DilateFilterType;
138  typedef AndImageFilter<OutputImageType,OutputImageType,OutputImageType> AndFilterType;
139 
140  typename OutputImageType::Pointer M0 = threshold->GetOutput();
141  typename OutputImageType::Pointer Mn = OutputImageType::New();
142  Mn->SetRegions( M0->GetLargestPossibleRegion() );
143  Mn->SetSpacing( M0->GetSpacing() );
144  Mn->SetOrigin( M0->GetOrigin() );
145  Mn->SetDirection( M0->GetDirection() );
146  Mn->Allocate();
147 
148  typename OutputImageType::Pointer Mnplus1 = erode->GetOutput();
149 
150 
151  CrossType cross;
152  cross.SetRadius( 1 );
153  //unsigned long rad[3]={3,3,3};
154 
155  //ball2.SetRadius( rad );
156 
157  std::cout << "Conditional reconstruction..." << std::flush;
158  int iter = 0;
159  do
160  {
161  std::cout << "Iteration: " << iter++ << std::endl;
162  CopyImage( Mn, Mnplus1);
163 
164  typename DilateFilterType::Pointer dilater = DilateFilterType::New();
165  dilater->SetInput( Mn );
166  dilater->SetKernel( cross );
167 
168  try
169  {
170  dilater->Update();
171  }
172  catch( itk::ExceptionObject &e)
173  {
174  std::cerr << e;
175  return;
176  }
177 
178  typename AndFilterType::Pointer andfilter = AndFilterType::New();
179  andfilter->SetInput(0, M0);
180  andfilter->SetInput(1, dilater->GetOutput() );
181 
182  try
183  {
184  andfilter->Update();
185  }
186  catch( itk::ExceptionObject &e)
187  {
188  std::cerr << e;
189  return;
190  }
191 
192  Mnplus1 = andfilter->GetOutput();
193 
194  /*
195  #ifdef DEBUG_ME
196  {
197  WriterType::Pointer writer = WriterType::New();
198  writer->SetInput( andfilter->GetOutput() );
199  char filename[512];
200  sprintf( filename, "CondReconstruction_iter_%d.hdr", iter);
201  writer->SetFileName( filename );
202  writer->Update();
203  }
204  #endif*/
205 
206  } while( !CompareImages( Mn, Mnplus1) && iter < m_MaxNumIterations );
207 
208  std::cout << "Done." << std::endl;
209 
210 #ifdef DEBUG_ME
211  {
212  WriterType::Pointer writer = WriterType::New();
213  writer->SetInput( Mn );
214  writer->SetFileName( "AfterCondReconstruction.hdr" );
215  writer->Update();
216  }
217 #endif
218  // now fill the holes
219 
220  typename itk::VotingBinaryIterativeHoleFillingImageFilter< OutputImageType >::Pointer filler =
221  itk::VotingBinaryIterativeHoleFillingImageFilter< OutputImageType >::New();
222  filler->SetInput( Mn );
223  filler->SetMaximumNumberOfIterations (1000);
224 
225  std::cout << "Filling the holes..." << std::flush;
226  try
227  {
228  filler->Update();
229  }
230  catch( itk::ExceptionObject &e)
231  {
232  std::cerr << e;
233  return;
234  }
235  std::cout << "Done." << std::endl;
236 
237  typename OutputImageType::Pointer outputImage =
238  static_cast< OutputImageType * >(this->ProcessObject::GetPrimaryOutput());
239 
240  outputImage->SetSpacing( filler->GetOutput()->GetSpacing() ); // Set the image spacing
241  outputImage->SetOrigin( filler->GetOutput()->GetOrigin() ); // Set the image origin
242  outputImage->SetLargestPossibleRegion( filler->GetOutput()->GetLargestPossibleRegion());
243  outputImage->SetBufferedRegion( filler->GetOutput()->GetLargestPossibleRegion() );
244  outputImage->SetDirection( filler->GetOutput()->GetDirection() );
245  outputImage->Allocate();
246 
247  itk::ImageRegionIterator<OutputImageType> itIn( filler->GetOutput(), filler->GetOutput()->GetLargestPossibleRegion() );
248  itk::ImageRegionIterator<OutputImageType> itOut( outputImage, outputImage->GetLargestPossibleRegion() );
249 
250  while( !itIn.IsAtEnd() )
251  {
252  itOut.Set(itIn.Get());
253  ++itIn;
254  ++itOut;
255  }
256 
257  }
258 
259 
260  template< class TOutputImagePixelType >
261  void BrainMaskExtractionImageFilter< TOutputImagePixelType >
262  ::CopyImage( typename OutputImageType::Pointer target, typename OutputImageType::Pointer source)
263  {
264  itk::ImageRegionConstIterator<OutputImageType> itIn( source, source->GetLargestPossibleRegion() );
265  itk::ImageRegionIterator<OutputImageType> itOut( target, target->GetLargestPossibleRegion() );
266 
267  while( !itOut.IsAtEnd() )
268  {
269  itOut.Set( itIn.Get() );
270  ++itIn;
271  ++itOut;
272  }
273  }
274 
275 
276  template< class TOutputImagePixelType >
277  bool BrainMaskExtractionImageFilter< TOutputImagePixelType >
278  ::CompareImages( typename OutputImageType::Pointer im1, typename OutputImageType::Pointer im2)
279  {
280  itk::ImageRegionConstIterator<OutputImageType> itIn( im1, im1->GetLargestPossibleRegion() );
281  itk::ImageRegionConstIterator<OutputImageType> itOut( im2, im2->GetLargestPossibleRegion() );
282 
283  while( !itOut.IsAtEnd() )
284  {
285  if( itOut.Value() != itIn.Value() )
286  {
287  return false;
288  }
289  ++itOut;
290  ++itIn;
291  }
292 
293  return true;
294  }
295 
296 
297  template< class TOutputImagePixelType >
298  int BrainMaskExtractionImageFilter< TOutputImagePixelType >
299  ::ComputeHistogram( typename InputImageType::Pointer image)
300  {
301  // IMPORTANT: IMAGE MUST BE UNSIGNED SHORT
302  int N=65535;
303  int* histogram = new int[N];
304  for( int i=0; i<N; i++)
305  {
306  histogram[i] = 0;
307  }
308 
309  itk::ImageRegionConstIterator<InputImageType> itIn( image, image->GetLargestPossibleRegion() );
310 
311  long totVoxels = 0;
312 
313  int max = -1;
314  int min = 9999999;
315  while( !itIn.IsAtEnd() )
316  {
317  histogram[ (int)(itIn.Value()) ]++;
318 
319  if( itIn.Value()>max )
320  {
321  max = itIn.Value();
322  }
323 
324  if( itIn.Value()<min )
325  {
326  min = itIn.Value();
327  }
328  ++itIn;
329  ++totVoxels;
330  }
331 
332 
333  //int EPS = 1;
334  int seuil = 0;
335  //int newseuil = (max + min)/2;
336 
337  N = max;
338 
339  double V = 0.0;
340  int S = 0;
341 
342  double mean = 0.0;
343  for( int i=min; i<=max; i++)
344  {
345  mean += (double)(i) * (double)histogram[i];
346  }
347  mean /= (double)totVoxels;
348 
349  //std::cout << "Min/Max/Mean: " << min << "/" << max << "/" << mean << std::endl;
350 
351  //while( abs(newseuil - seuil)>EPS )
352  for( seuil = min; seuil<=max; seuil++)
353  {
354  //seuil = newseuil;
355 
356  // compute the classes:
357  double mean0 = 0.0;
358  double mean1 = 0.0;
359  //double std0 = 0.0;
360  //double std1 = 0.0;
361 
362  double num0 = 0.0;
363  double num1 = 0.0;
364  for( int i=min; i<seuil; i++)
365  {
366  //mean0 += (double)(histogram[i])/(double)(totVoxels) * (double)(i);
367  mean0 += (double)histogram[i] * (double)i;
368  num0 += (double)histogram[i];
369  }
370 
371  for( int i=seuil; i<max; i++)
372  {
373  //mean1 += (double)(histogram[i])/(double)(totVoxels) * (double)(i);
374  mean1 += (double)histogram[i] * (double)i;
375  num1 += (double)histogram[i];
376  }
377 
378  if( num0 )
379  {
380  mean0/=num0;
381  }
382  if( num1 )
383  {
384  mean1/=num1;
385  }
386 
387  /*
388  for( int i=0; i<seuil; i++)
389  {
390  std0 += (double)histogram[i]/(double)totVoxels* ((double)i - mean0)*((double)i - mean0);
391  //std0 += (double)histogram[i]/num0* ((double)i - mean0)*((double)i - mean0);
392  }
393 
394  for( int i=seuil; i<65536; i++)
395  {
396  std1 += (double)histogram[i]/(double)totVoxels* ((double)i - mean1)* ((double)i - mean1);
397  //std1 += (double)histogram[i]/num1* ((double)i - mean1)* ((double)i - mean1);
398  }
399  */
400 
401 
402  //std0 = sqrt(std0);
403  //std1 = sqrt(std1);
404 
405  //std::cout << "Classe 0: " << mean0 << " " << std0 << std::endl;
406  //std::cout << "Classe 1: " << mean1 << " " << std1 << std::endl;
407 
408  //std::cout << "Classe 0: " << mean0 << std::endl;
409  //std::cout << "Classe 1: " << mean1 << std::endl;
410 
411  //double alpha = std0*std1/(std0+std1);
412 
413  //std::cout << "alpha: " << alpha << " " << std0 << " " << std1 << std::endl;
414  //newseuil = (int) (alpha*(mean0/std0 + mean1/std1));
415 
416  //std::cout << "New threshold: " << newseuil << std::endl;
417 
418  /*
419  if( newseuil > 65535 || newseuil<0)
420  {
421  std::cerr << "Error: threshold is too low or high, exiting" << std::endl;
422  return -1;
423  }*/
424 
425  double Vm = num0 * (mean0 - mean)*(mean0 - mean) + num1*(mean1 - mean)*(mean1 - mean);
426  if( Vm > V )
427  {
428  V = Vm;
429  S = seuil;
430  //std::cout << "New seuil: " << S << std::endl;
431  //getchar();
432  }
433 
434  }
435 
436  delete [] histogram;
437 
438  std::cout << "Seuil: " << S << std::endl;
439 
440  return S;
441  }
442 }
443 
444 #endif // __itkBrainMaskExtractionImageFilter_txx