[Insight-users] EM-algorithm

[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 Hotmail sur i-mode™ : envoyez et recevez des e-mails depuis votre 
téléphone portable ! http://www.msn.fr/hotmailimode/