[Insight-users] help needed for Expectation Maximization algorithm

[Seniha Esen Yuksel]

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 (although all of them exist and work seperately, my combination is resisting to work). Below is the code,  EM is terminating with  * termination code=  1 *. And the clusters are getting *0* mean and variance values. I would appreciate if you could comment on what is wrong. The image I am using has two nicely gaussian classes, not to be cut, I will send it in a the next email. Thanks a lot for help! 
  Esen
   
  Results:
  Cluster[0]
    Parameters:
         [0, 0]
    Proportion:          1
Cluster[1]
    Parameters:
         [0, 0]
    Proportion:          2.22156e-053
Code:
   
  #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( argv[1] );
   reader->SetFileName( "Simulatedchecker.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( bin, 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;
  
}

		
---------------------------------
 Yahoo! FareChase - Search multiple travel sites in one click.  
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20051119/41ff99c0/attachment.html