[Insight-users] EM-algorithm

[sami bourouis]

>From: "sami bourouis" <sami.bourouis at hotmail.com>
>To: insight-users-request at itk.org, insight-users at itk.org
>Subject: [Insight-users] EM-algorithm
>Date: Wed, 06 Dec 2006 09:56:17 +0000
>
>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/
>
>_______________________________________________
>Insight-users mailing list
>Insight-users at itk.org
>http://www.itk.org/mailman/listinfo/insight-users

_________________________________________________________________
MSN Hotmail sur i-mode™ : envoyez et recevez des e-mails depuis votre 
téléphone portable ! http://www.msn.fr/hotmailimode/