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
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