[Insight-users] EM-classification

[sami bourouis]

Hi all,
I have a problem with the Expectation Maximization algorithm.
I read the tutorial and looked at the mailing list but there is not a single 
example of how I would read an image, generate its histogram, feed it into 
the EM and get the results.
Below is the code,  EM is terminating with  * termination code=  1 *.

I would appreciate if you could comment on what is wrong.
Thanks a lot for help!




#include "itkScalarImageToHistogramGenerator.h"
#include "itkImage.h"
#include "itkVector.h"
#include "itkListSample.h"
#include "itkGaussianMixtureModelComponent.h"
#include "itkExpectationMaximizationMixtureModelEstimator.h"
// Software Guide : EndCodeSnippet
#include "itkImageFileReader.h"
int main( int argc, char * argv [] )
{
/*
if( argc < 2 )
   {
   std::cerr << "Missing command line arguments" << std::endl;
   std::cerr << "Usage :  ImageHistogram1  inputImageFileName " << 
std::endl;
   return -1;
   }
*/
   typedef unsigned char       PixelType;
const unsigned int          Dimension = 2;
   typedef itk::Image<PixelType, Dimension > ImageType;
   typedef itk::ImageFileReader< ImageType > ReaderType;
   ReaderType::Pointer reader = ReaderType::New();

  reader->SetFileName( "BrainT1Slice.png" );
   try
   {
   reader->Update();
   }
catch( itk::ExceptionObject & excp )
   {
   std::cerr << "Problem encoutered while reading image file : " << argv[1] 
<< std::endl;
   std::cerr << excp << std::endl;
   return -1;
   }
   typedef itk::Statistics::ScalarImageToHistogramGenerator<
                                                   ImageType
                                                         >   
HistogramGeneratorType;
   HistogramGeneratorType::Pointer histogramGenerator = 
HistogramGeneratorType::New();
   histogramGenerator->SetInput(  reader->GetOutput() );
   histogramGenerator->SetNumberOfBins( 255 );
histogramGenerator->SetMarginalScale( 10.0 );
histogramGenerator->Compute();
   typedef HistogramGeneratorType::HistogramType  HistogramType;
   const HistogramType * histogram = histogramGenerator->GetOutput();
   const unsigned int histogramSize = histogram->Size();
   std::cout << "Histogram size " << histogramSize << std::endl;
   unsigned int bin;
for( bin=0; bin < histogramSize; bin++ )
   {
   std::cout << "bin = " << bin << " frequency = ";
   std::cout << histogram->GetFrequency( bin, 0 ) << std::endl;
   }

   //*******************************
//=====end of histogram ========
//======start EM ================
  unsigned int numberOfClasses = 2;
typedef itk::Vector< double, 1 > MeasurementVectorType;
typedef itk::Statistics::ListSample< MeasurementVectorType > SampleType;
SampleType::Pointer sample = SampleType::New();

for ( unsigned int i = 0 ; i < histogramSize ; ++i )
   {
      sample->PushBack( histogram->GetFrequency( i, 0 ) );//0 is for first 
channel
   }
    //initial parameters for EM
typedef itk::Array< double > ParametersType;
ParametersType params( 2 );
   std::vector< ParametersType > initialParameters( numberOfClasses );
params[0] = 50.0;
params[1] = 100.0;
initialParameters[0] = params;
   params[0] = 200.0;
params[1] = 150.0;
initialParameters[1] = params;

// EM will estimate with Gaussians
typedef itk::Statistics::GaussianMixtureModelComponent< SampleType >
   ComponentType;
   std::vector< ComponentType::Pointer > components;
for ( unsigned int i = 0 ; i < numberOfClasses ; i++ )
   {
    components.push_back( ComponentType::New() );
   (components[i])->SetSample( sample );
   (components[i])->SetParameters( initialParameters[i] );
   }
   typedef itk::Statistics::ExpectationMaximizationMixtureModelEstimator<
                          SampleType > EstimatorType;
EstimatorType::Pointer estimator = EstimatorType::New();
   estimator->SetSample( sample );
estimator->SetMaximumIteration( 300 );
   itk::Array< double > initialProportions(numberOfClasses);
initialProportions[0] = 0.5;
initialProportions[1] = 0.5;
   estimator->SetInitialProportions( initialProportions );

for ( unsigned int i = 0 ; i < numberOfClasses ; i++)
   {
   estimator->AddComponent( (ComponentType::Superclass*)
                            (components[i]).GetPointer() );
   }

try
   {
   estimator->Update();
   }
catch( itk::ExceptionObject & excp )
   {
   std::cerr << "Problem encoutered while estimating" << std::endl;
   std::cerr << excp << std::endl;
   return -1;
   }

int termincode = estimator->GetTerminationCode (); // CONVERGED = 0, 
NOT_CONVERGED = 1 }
  std::cout<<"termination code=  "<<termincode<<std::endl;

for ( unsigned int i = 0 ; i < numberOfClasses ; i++ )
   {
   std::cout << "Cluster[" << i << "]" << std::endl;
   std::cout << "    Parameters:" << std::endl;
   std::cout << "         " << (components[i])->GetFullParameters()
             << std::endl;
   std::cout << "    Proportion: ";
   std::cout << "         " << (*estimator->GetProportions())[i] << 
std::endl;
   }
   return 0;

}

_________________________________________________________________
MSN Messenger : discutez en direct avec vos amis ! 
http://www.msn.fr/msger/default.asp