mitk: itkBrainMaskExtractionImageFilter.txx Source File

Go to the documentation of this file.
00001 /*=========================================================================
00002 
00003 Program:   Insight Segmentation & Registration Toolkit
00004 Language:  C++
00005 Date:      $Date: 2006-07-19 15:11:41 $
00006 Version:   $Revision: 1.11 $
00007 
00008 Copyright (c) Insight Software Consortium. All rights reserved.
00009 See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
00010 
00011 This software is distributed WITHOUT ANY WARRANTY; without even 
00012 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
00013 PURPOSE.  See the above copyright notices for more information.
00014 
00015 =========================================================================*/
00016 #ifndef __itkBrainMaskExtractionImageFilter_txx
00017 #define __itkBrainMaskExtractionImageFilter_txx
00018 
00019 #include "itkBrainMaskExtractionImageFilter.h"
00020 
00021 #include <itkBinaryThresholdImageFilter.h>
00022 #include <itkBinaryErodeImageFilter.h>
00023 #include <itkBinaryDilateImageFilter.h>
00024 #include <itkVotingBinaryHoleFillingImageFilter.h>
00025 #include <itkVotingBinaryIterativeHoleFillingImageFilter.h>
00026 #include <itkBinaryBallStructuringElement.h>
00027 #include <itkBinaryCrossStructuringElement.h>
00028 #include <itkAndImageFilter.h>
00029 #include <itkRecursiveGaussianImageFilter.h> 
00030 #include <itkMaskImageFilter.h>
00031 
00032 namespace itk {
00033 
00034   template< class TOutputImagePixelType >
00035     BrainMaskExtractionImageFilter< TOutputImagePixelType >
00036     ::BrainMaskExtractionImageFilter()
00037   {
00038     // At least 1 inputs is necessary for a vector image.
00039     // For images added one at a time we need at least six
00040     this->SetNumberOfRequiredInputs( 1 ); 
00041   }
00042 
00043   template< class TOutputImagePixelType >
00044     void BrainMaskExtractionImageFilter< TOutputImagePixelType >
00045     ::GenerateData()
00046   {
00047 
00048     // gaussian filter
00049     typename itk::RecursiveGaussianImageFilter<InputImageType,InputImageType>::Pointer gaussian
00050       = itk::RecursiveGaussianImageFilter<InputImageType,InputImageType>::New();
00051     gaussian->SetInput( this->GetInput(0) );
00052     gaussian->SetSigma( 1.0 );
00053 
00054     try
00055     {
00056       gaussian->Update();
00057     }
00058     catch( itk::ExceptionObject &e)
00059     {
00060       std::cerr << e;
00061       return;
00062     }
00063 
00064     // threshold the image
00065     typename itk::BinaryThresholdImageFilter<InputImageType,OutputImageType>::Pointer threshold = 
00066       itk::BinaryThresholdImageFilter<InputImageType,OutputImageType>::New();
00067 
00068     threshold->SetInput( gaussian->GetOutput() );
00069 
00070     int seuil = static_cast<int>( ComputeHistogram( gaussian->GetOutput() ) );
00071     threshold->SetLowerThreshold( seuil );
00072 
00073     std::cout << "Thresholding..." << std::flush;
00074     try
00075     {
00076       threshold->Update();
00077     }
00078     catch( itk::ExceptionObject &e)
00079     {
00080       std::cerr << e;
00081       return;
00082     }
00083     std::cout << "Done." << std::endl;
00084 
00085 #ifdef DEBUG_ME
00086     {
00087       WriterType::Pointer writer = WriterType::New();
00088       writer->SetInput( threshold->GetOutput() );
00089       writer->SetFileName( "AfterThreshold.hdr" );
00090       writer->Update();
00091     }
00092 #endif
00093 
00094     // erode to remove background noise
00095     typedef itk::BinaryBallStructuringElement<int, 3> StructuralElementType;
00096     StructuralElementType ball;
00097 
00098     typename itk::BinaryErodeImageFilter<OutputImageType,OutputImageType,StructuralElementType>::Pointer erode = 
00099       itk::BinaryErodeImageFilter<OutputImageType,OutputImageType,StructuralElementType>::New();
00100 
00101     ball.SetRadius( 3 );
00102 
00103     erode->SetInput( threshold->GetOutput() );
00104     erode->SetKernel( ball );
00105 
00106     std::cout << "Eroding..." << std::flush;
00107     try
00108     {
00109       erode->Update();
00110     }
00111     catch( itk::ExceptionObject &e)
00112     {
00113       std::cerr << e;
00114       return;
00115     }
00116     std::cout << "Done." << std::endl;
00117 
00118 #ifdef DEBUG_ME
00119     {
00120       WriterType::Pointer writer = WriterType::New();
00121       writer->SetInput( erode->GetOutput() );
00122       writer->SetFileName( "AfterErode.hdr" );
00123       writer->Update();
00124     }
00125 #endif
00126 
00127 
00128     typedef BinaryCrossStructuringElement<int, 3> CrossType;
00129 
00130     typedef BinaryDilateImageFilter<OutputImageType,OutputImageType,CrossType> DilateFilterType;
00131     typedef AndImageFilter<OutputImageType,OutputImageType,OutputImageType> AndFilterType;
00132 
00133     typename OutputImageType::Pointer M0 = threshold->GetOutput();
00134     typename OutputImageType::Pointer Mn = OutputImageType::New();
00135     Mn->SetRegions( M0->GetLargestPossibleRegion() );
00136     Mn->SetSpacing( M0->GetSpacing() );
00137     Mn->SetOrigin( M0->GetOrigin() );
00138     Mn->Allocate();
00139 
00140     typename OutputImageType::Pointer Mnplus1 = erode->GetOutput();
00141 
00142 
00143     CrossType cross;
00144     cross.SetRadius( 1 );
00145     //unsigned long rad[3]={3,3,3};
00146 
00147     //ball2.SetRadius( rad );
00148 
00149     std::cout << "Conditional reconstruction..." << std::flush;
00150     int iter = 0;
00151     do
00152     {
00153       std::cout << "Iteration: " << iter++ << std::endl;
00154       CopyImage( Mn, Mnplus1);
00155 
00156       typename DilateFilterType::Pointer dilater = DilateFilterType::New();
00157       dilater->SetInput( Mn );
00158       dilater->SetKernel( cross );
00159 
00160       try
00161       {
00162         dilater->Update();
00163       }
00164       catch( itk::ExceptionObject &e)
00165       {
00166         std::cerr << e;
00167         return;
00168       }
00169 
00170       typename AndFilterType::Pointer andfilter = AndFilterType::New();
00171       andfilter->SetInput(0, M0);
00172       andfilter->SetInput(1, dilater->GetOutput() );
00173 
00174       try
00175       {
00176         andfilter->Update();
00177       }
00178       catch( itk::ExceptionObject &e)
00179       {
00180         std::cerr << e;
00181         return;
00182       }
00183 
00184       Mnplus1 = andfilter->GetOutput();
00185 
00186       /*
00187       #ifdef DEBUG_ME
00188       {
00189       WriterType::Pointer writer = WriterType::New();
00190       writer->SetInput( andfilter->GetOutput() );
00191       char filename[512];
00192       sprintf( filename, "CondReconstruction_iter_%d.hdr", iter);
00193       writer->SetFileName( filename );
00194       writer->Update();
00195       }
00196       #endif*/
00197 
00198     } while( !CompareImages( Mn, Mnplus1) );
00199 
00200     std::cout << "Done." << std::endl;
00201 
00202 #ifdef DEBUG_ME
00203     {
00204       WriterType::Pointer writer = WriterType::New();
00205       writer->SetInput( Mn );
00206       writer->SetFileName( "AfterCondReconstruction.hdr" );
00207       writer->Update();
00208     }
00209 #endif
00210     // now fill the holes
00211 
00212     typename itk::VotingBinaryIterativeHoleFillingImageFilter< OutputImageType >::Pointer filler =
00213       itk::VotingBinaryIterativeHoleFillingImageFilter< OutputImageType >::New();
00214     filler->SetInput( Mn );
00215     filler->SetMaximumNumberOfIterations (1000);
00216 
00217     std::cout << "Filling the holes..." << std::flush;
00218     try
00219     {
00220       filler->Update();
00221     }
00222     catch( itk::ExceptionObject &e)
00223     {
00224       std::cerr << e;
00225       return;
00226     }
00227     std::cout << "Done." << std::endl;
00228 
00229     typename OutputImageType::Pointer outputImage = 
00230       static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0));
00231 
00232     outputImage->SetSpacing( filler->GetOutput()->GetSpacing() );   // Set the image spacing
00233     outputImage->SetOrigin( filler->GetOutput()->GetOrigin() );     // Set the image origin
00234     outputImage->SetDirection( filler->GetOutput()->GetDirection() );  // Set the image direction
00235     outputImage->SetLargestPossibleRegion( filler->GetOutput()->GetLargestPossibleRegion());
00236     outputImage->SetBufferedRegion( filler->GetOutput()->GetLargestPossibleRegion() );
00237     outputImage->Allocate();
00238 
00239     itk::ImageRegionIterator<OutputImageType> itIn( filler->GetOutput(), filler->GetOutput()->GetLargestPossibleRegion() );
00240     itk::ImageRegionIterator<OutputImageType> itOut( outputImage, outputImage->GetLargestPossibleRegion() );
00241 
00242     while( !itIn.IsAtEnd() )
00243     {
00244         itOut.Set(itIn.Get());
00245         ++itIn;
00246         ++itOut;
00247     }
00248 
00249   }
00250 
00251 
00252   template< class TOutputImagePixelType >
00253     void BrainMaskExtractionImageFilter< TOutputImagePixelType >
00254     ::CopyImage( typename OutputImageType::Pointer target, typename OutputImageType::Pointer source)
00255   {
00256     itk::ImageRegionConstIterator<OutputImageType> itIn( source, source->GetLargestPossibleRegion() );
00257     itk::ImageRegionIterator<OutputImageType> itOut( target, target->GetLargestPossibleRegion() );
00258 
00259     while( !itOut.IsAtEnd() )
00260     {
00261       itOut.Set( itIn.Get() );
00262       ++itIn;
00263       ++itOut;
00264     }
00265   }
00266 
00267 
00268   template< class TOutputImagePixelType >
00269     bool BrainMaskExtractionImageFilter< TOutputImagePixelType >
00270     ::CompareImages( typename OutputImageType::Pointer im1, typename OutputImageType::Pointer im2)
00271   {
00272     itk::ImageRegionConstIterator<OutputImageType> itIn( im1, im1->GetLargestPossibleRegion() );
00273     itk::ImageRegionConstIterator<OutputImageType> itOut( im2, im2->GetLargestPossibleRegion() );
00274 
00275     while( !itOut.IsAtEnd() )
00276     {
00277       if( itOut.Value() != itIn.Value() )
00278       {
00279         return false;
00280       }
00281       ++itOut;
00282       ++itIn;
00283     }
00284 
00285     return true;
00286   }
00287 
00288 
00289   template< class TOutputImagePixelType >
00290     int BrainMaskExtractionImageFilter< TOutputImagePixelType >
00291     ::ComputeHistogram( typename InputImageType::Pointer image)
00292   {
00293     // IMPORTANT: IMAGE MUST BE UNSIGNED SHORT
00294     int N=65535;
00295     int* histogram = new int[N];
00296     for( int i=0; i<N; i++)
00297     {
00298       histogram[i] = 0;
00299     }
00300 
00301     itk::ImageRegionConstIterator<InputImageType> itIn( image, image->GetLargestPossibleRegion() );
00302 
00303     long totVoxels = 0;
00304 
00305     int max = -1;
00306     int min = 9999999;
00307     while( !itIn.IsAtEnd() )
00308     {
00309       histogram[ (int)(itIn.Value()) ]++;
00310 
00311       if( itIn.Value()>max )
00312       {
00313         max = itIn.Value();
00314       }
00315 
00316       if( itIn.Value()<min )
00317       {
00318         min = itIn.Value();
00319       }
00320       ++itIn;
00321       ++totVoxels;
00322     }
00323 
00324 
00325     //int EPS = 1;
00326     int seuil = 0;
00327     //int newseuil = (max + min)/2;
00328 
00329     N = max;
00330 
00331     double V = 0.0;
00332     int S = 0;
00333 
00334     double mean = 0.0;
00335     for( int i=min; i<=max; i++)
00336     {
00337       mean += (double)(i) * (double)histogram[i];
00338     }
00339     mean /= (double)totVoxels;
00340 
00341     //std::cout << "Min/Max/Mean: " << min << "/" << max << "/" << mean << std::endl;
00342 
00343     //while( abs(newseuil - seuil)>EPS )
00344     for( seuil = min; seuil<=max; seuil++)
00345     {
00346       //seuil = newseuil;
00347 
00348       // compute the classes:
00349       double mean0 = 0.0;
00350       double mean1 = 0.0;
00351       //double std0 = 0.0;
00352       //double std1 = 0.0;
00353 
00354       double num0 = 0.0;
00355       double num1 = 0.0;
00356       for( int i=min; i<seuil; i++)
00357       {
00358         //mean0 += (double)(histogram[i])/(double)(totVoxels) * (double)(i);
00359         mean0 += (double)histogram[i] * (double)i;
00360         num0  += (double)histogram[i];
00361       }
00362 
00363       for( int i=seuil; i<max; i++)
00364       {
00365         //mean1 += (double)(histogram[i])/(double)(totVoxels) * (double)(i);
00366         mean1 += (double)histogram[i] * (double)i;
00367         num1  += (double)histogram[i];
00368       }
00369 
00370       if( num0 )
00371       {
00372         mean0/=num0;
00373       }
00374       if( num1 )
00375       {
00376         mean1/=num1;
00377       }
00378 
00379       /*
00380       for( int i=0; i<seuil; i++)
00381       {
00382       std0 += (double)histogram[i]/(double)totVoxels* ((double)i - mean0)*((double)i - mean0);
00383       //std0 += (double)histogram[i]/num0* ((double)i - mean0)*((double)i - mean0);
00384       }
00385 
00386       for( int i=seuil; i<65536; i++)
00387       {
00388       std1 += (double)histogram[i]/(double)totVoxels* ((double)i - mean1)* ((double)i - mean1);
00389       //std1 += (double)histogram[i]/num1* ((double)i - mean1)* ((double)i - mean1);
00390       }
00391       */
00392 
00393 
00394       //std0 = sqrt(std0);
00395       //std1 = sqrt(std1);
00396 
00397       //std::cout << "Classe 0: " << mean0 << " " << std0 << std::endl;
00398       //std::cout << "Classe 1: " << mean1 << " " << std1 << std::endl;
00399 
00400       //std::cout << "Classe 0: " << mean0 << std::endl;
00401       //std::cout << "Classe 1: " << mean1 << std::endl;
00402 
00403       //double alpha = std0*std1/(std0+std1);
00404 
00405       //std::cout << "alpha: " << alpha << " " << std0 << " " << std1 << std::endl;
00406       //newseuil = (int) (alpha*(mean0/std0 + mean1/std1));             
00407 
00408       //std::cout << "New threshold: " << newseuil << std::endl;
00409 
00410       /*
00411       if( newseuil > 65535 || newseuil<0)
00412       {
00413       std::cerr << "Error: threshold is too low or high, exiting" << std::endl;
00414       return -1;
00415       }*/
00416 
00417       double Vm = num0 * (mean0 - mean)*(mean0 - mean) + num1*(mean1 - mean)*(mean1 - mean);
00418       if( Vm > V )
00419       {
00420         V = Vm;
00421         S = seuil;
00422         //std::cout << "New seuil: " << S << std::endl;
00423         //getchar();
00424       }
00425 
00426     }
00427 
00428     delete [] histogram;
00429 
00430     std::cout << "Seuil: " << S << std::endl;
00431 
00432     return S;
00433   }
00434 }
00435 
00436 #endif // __itkBrainMaskExtractionImageFilter_txx