1 /*===================================================================
3 The Medical Imaging Interaction Toolkit (MITK)
5 Copyright (c) German Cancer Research Center,
6 Division of Medical and Biological Informatics.
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
13 See LICENSE.txt or http://www.mitk.org for details.
15 ===================================================================*/
17 /*===================================================================
19 This file is based heavily on a corresponding ITK filter.
21 ===================================================================*/
22 #ifndef __itkBrainMaskExtractionImageFilter_txx
23 #define __itkBrainMaskExtractionImageFilter_txx
25 #include "itkBrainMaskExtractionImageFilter.h"
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>
40 template< class TOutputImagePixelType >
41 BrainMaskExtractionImageFilter< TOutputImagePixelType >
42 ::BrainMaskExtractionImageFilter()
43 : m_MaxNumIterations(itk::NumericTraits<int>::max())
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 );
50 template< class TOutputImagePixelType >
51 void BrainMaskExtractionImageFilter< TOutputImagePixelType >
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 );
65 catch( itk::ExceptionObject &e)
71 // threshold the image
72 typename itk::BinaryThresholdImageFilter<InputImageType,OutputImageType>::Pointer threshold =
73 itk::BinaryThresholdImageFilter<InputImageType,OutputImageType>::New();
75 threshold->SetInput( gaussian->GetOutput() );
77 int seuil = static_cast<int>( ComputeHistogram( gaussian->GetOutput() ) );
78 threshold->SetLowerThreshold( seuil );
80 std::cout << "Thresholding..." << std::flush;
85 catch( itk::ExceptionObject &e)
90 std::cout << "Done." << std::endl;
94 WriterType::Pointer writer = WriterType::New();
95 writer->SetInput( threshold->GetOutput() );
96 writer->SetFileName( "AfterThreshold.hdr" );
101 // erode to remove background noise
102 typedef itk::BinaryBallStructuringElement<int, 3> StructuralElementType;
103 StructuralElementType ball;
105 typename itk::BinaryErodeImageFilter<OutputImageType,OutputImageType,StructuralElementType>::Pointer erode =
106 itk::BinaryErodeImageFilter<OutputImageType,OutputImageType,StructuralElementType>::New();
110 erode->SetInput( threshold->GetOutput() );
111 erode->SetKernel( ball );
113 std::cout << "Eroding..." << std::flush;
118 catch( itk::ExceptionObject &e)
123 std::cout << "Done." << std::endl;
127 WriterType::Pointer writer = WriterType::New();
128 writer->SetInput( erode->GetOutput() );
129 writer->SetFileName( "AfterErode.hdr" );
135 typedef BinaryCrossStructuringElement<int, 3> CrossType;
137 typedef BinaryDilateImageFilter<OutputImageType,OutputImageType,CrossType> DilateFilterType;
138 typedef AndImageFilter<OutputImageType,OutputImageType,OutputImageType> AndFilterType;
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() );
148 typename OutputImageType::Pointer Mnplus1 = erode->GetOutput();
152 cross.SetRadius( 1 );
153 //unsigned long rad[3]={3,3,3};
155 //ball2.SetRadius( rad );
157 std::cout << "Conditional reconstruction..." << std::flush;
161 std::cout << "Iteration: " << iter++ << std::endl;
162 CopyImage( Mn, Mnplus1);
164 typename DilateFilterType::Pointer dilater = DilateFilterType::New();
165 dilater->SetInput( Mn );
166 dilater->SetKernel( cross );
172 catch( itk::ExceptionObject &e)
178 typename AndFilterType::Pointer andfilter = AndFilterType::New();
179 andfilter->SetInput(0, M0);
180 andfilter->SetInput(1, dilater->GetOutput() );
186 catch( itk::ExceptionObject &e)
192 Mnplus1 = andfilter->GetOutput();
197 WriterType::Pointer writer = WriterType::New();
198 writer->SetInput( andfilter->GetOutput() );
200 sprintf( filename, "CondReconstruction_iter_%d.hdr", iter);
201 writer->SetFileName( filename );
206 } while( !CompareImages( Mn, Mnplus1) && iter < m_MaxNumIterations );
208 std::cout << "Done." << std::endl;
212 WriterType::Pointer writer = WriterType::New();
213 writer->SetInput( Mn );
214 writer->SetFileName( "AfterCondReconstruction.hdr" );
218 // now fill the holes
220 typename itk::VotingBinaryIterativeHoleFillingImageFilter< OutputImageType >::Pointer filler =
221 itk::VotingBinaryIterativeHoleFillingImageFilter< OutputImageType >::New();
222 filler->SetInput( Mn );
223 filler->SetMaximumNumberOfIterations (1000);
225 std::cout << "Filling the holes..." << std::flush;
230 catch( itk::ExceptionObject &e)
235 std::cout << "Done." << std::endl;
237 typename OutputImageType::Pointer outputImage =
238 static_cast< OutputImageType * >(this->ProcessObject::GetPrimaryOutput());
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();
247 itk::ImageRegionIterator<OutputImageType> itIn( filler->GetOutput(), filler->GetOutput()->GetLargestPossibleRegion() );
248 itk::ImageRegionIterator<OutputImageType> itOut( outputImage, outputImage->GetLargestPossibleRegion() );
250 while( !itIn.IsAtEnd() )
252 itOut.Set(itIn.Get());
260 template< class TOutputImagePixelType >
261 void BrainMaskExtractionImageFilter< TOutputImagePixelType >
262 ::CopyImage( typename OutputImageType::Pointer target, typename OutputImageType::Pointer source)
264 itk::ImageRegionConstIterator<OutputImageType> itIn( source, source->GetLargestPossibleRegion() );
265 itk::ImageRegionIterator<OutputImageType> itOut( target, target->GetLargestPossibleRegion() );
267 while( !itOut.IsAtEnd() )
269 itOut.Set( itIn.Get() );
276 template< class TOutputImagePixelType >
277 bool BrainMaskExtractionImageFilter< TOutputImagePixelType >
278 ::CompareImages( typename OutputImageType::Pointer im1, typename OutputImageType::Pointer im2)
280 itk::ImageRegionConstIterator<OutputImageType> itIn( im1, im1->GetLargestPossibleRegion() );
281 itk::ImageRegionConstIterator<OutputImageType> itOut( im2, im2->GetLargestPossibleRegion() );
283 while( !itOut.IsAtEnd() )
285 if( itOut.Value() != itIn.Value() )
297 template< class TOutputImagePixelType >
298 int BrainMaskExtractionImageFilter< TOutputImagePixelType >
299 ::ComputeHistogram( typename InputImageType::Pointer image)
301 // IMPORTANT: IMAGE MUST BE UNSIGNED SHORT
303 int* histogram = new int[N];
304 for( int i=0; i<N; i++)
309 itk::ImageRegionConstIterator<InputImageType> itIn( image, image->GetLargestPossibleRegion() );
315 while( !itIn.IsAtEnd() )
317 histogram[ (int)(itIn.Value()) ]++;
319 if( itIn.Value()>max )
324 if( itIn.Value()<min )
335 //int newseuil = (max + min)/2;
343 for( int i=min; i<=max; i++)
345 mean += (double)(i) * (double)histogram[i];
347 mean /= (double)totVoxels;
349 //std::cout << "Min/Max/Mean: " << min << "/" << max << "/" << mean << std::endl;
351 //while( abs(newseuil - seuil)>EPS )
352 for( seuil = min; seuil<=max; seuil++)
356 // compute the classes:
364 for( int i=min; i<seuil; i++)
366 //mean0 += (double)(histogram[i])/(double)(totVoxels) * (double)(i);
367 mean0 += (double)histogram[i] * (double)i;
368 num0 += (double)histogram[i];
371 for( int i=seuil; i<max; i++)
373 //mean1 += (double)(histogram[i])/(double)(totVoxels) * (double)(i);
374 mean1 += (double)histogram[i] * (double)i;
375 num1 += (double)histogram[i];
388 for( int i=0; i<seuil; i++)
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);
394 for( int i=seuil; i<65536; i++)
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);
405 //std::cout << "Classe 0: " << mean0 << " " << std0 << std::endl;
406 //std::cout << "Classe 1: " << mean1 << " " << std1 << std::endl;
408 //std::cout << "Classe 0: " << mean0 << std::endl;
409 //std::cout << "Classe 1: " << mean1 << std::endl;
411 //double alpha = std0*std1/(std0+std1);
413 //std::cout << "alpha: " << alpha << " " << std0 << " " << std1 << std::endl;
414 //newseuil = (int) (alpha*(mean0/std0 + mean1/std1));
416 //std::cout << "New threshold: " << newseuil << std::endl;
419 if( newseuil > 65535 || newseuil<0)
421 std::cerr << "Error: threshold is too low or high, exiting" << std::endl;
425 double Vm = num0 * (mean0 - mean)*(mean0 - mean) + num1*(mean1 - mean)*(mean1 - mean);
430 //std::cout << "New seuil: " << S << std::endl;
438 std::cout << "Seuil: " << S << std::endl;
444 #endif // __itkBrainMaskExtractionImageFilter_txx