[Insight-users] Help required on ITK Expectation Maximization Algorithm

[Prashanth Ravindran]

Skipped content of type multipart/alternative-------------- next part --------------
# include "itkRescaleIntensityImageFilter.h"
#include "itkVector.h"
#include "itkListSample.h"
#include "itkGaussianMixtureModelComponent.h"
#include "itkExpectationMaximizationMixtureModelEstimator.h"
#include "itkNormalVariateGenerator.h"
# include "itkImageToListGenerator.h"
# include "itkImageRegionIterator.h"
# include "itkImageRegionConstIterator.h"
# include <ctime>

//typedef itk::Vector<double,2> PixelType;
typedef float PixelType;
typedef itk::Image<PixelType,2> ImageType;
typedef unsigned char LabelType;
typedef itk::Image<unsigned char,2> LabelImageType;
typedef itk::ImageFileReader<ImageType> ReaderType;
typedef itk::BayesianClassifierInitializationImageFilter< ImageType > 
                                                BayesianInitializerType;



int main(int argc, char * argv[] )
{

// EM- test

  typedef itk::Vector<double,2> PixelType;
  typedef itk::Image<PixelType,2> ImageType;
  typedef itk::ImageFileReader<ImageType> ReaderType;
  ReaderType::Pointer reader = ReaderType::New();
  reader->SetFileName(argv[1]);
  try{
	  reader->Update();
  }
  catch(itk::ExceptionObject & excep)
  {
	  std::cerr<<excep;
	  return -1;
  }
   
  unsigned int numberOfClasses = 3;
  typedef itk::Vector< double, 2 > MeasurementVectorType;
  //typedef PixelType MeasurementVectorType;
  typedef itk::Statistics::ListSample< MeasurementVectorType > SampleType;
  typedef itk::ImageRegionConstIterator<ImageType> ConstIteratorType;
  ConstIteratorType sourceit(reader->GetOutput(), (reader->GetOutput())->GetLargestPossibleRegion());
  sourceit.GoToBegin();
  SampleType::Pointer sample = SampleType::New();
  sample->SetMeasurementVectorSize( 2 ); // length of measurement vectors
  //   
                                    // in the sample.
  //std::clock_t start = clock();
  MeasurementVectorType mv;
    
  typedef itk::Statistics::NormalVariateGenerator NormalGeneratorType;
  NormalGeneratorType::Pointer normalGenerator1 = NormalGeneratorType::New();
  NormalGeneratorType::Pointer normalGenerator2 = NormalGeneratorType::New();

  normalGenerator1->Initialize( 101 );
  normalGenerator2->Initialize( 110 );

  //MeasurementVectorType mv;
  double mean1 = 100;
  double mean2 = 100;
  double standardDeviation1 = 30;
  double standardDeviation2 = 30;
  for ( unsigned int i = 0 ; i < 100 ; ++i )
    {
    mv[0] = ( normalGenerator1->GetVariate() * standardDeviation1 ) + mean1;
	mv[1] = ( normalGenerator2->GetVariate() * standardDeviation2 ) + mean2;
	sample->PushBack( mv );
    }

  normalGenerator1->Initialize( 3024 );
  normalGenerator2->Initialize(324 );
  mean1 = 200;
  mean2 = 200;
  standardDeviation1 = 30;
  standardDeviation2 = 30;
  for ( unsigned int i = 0 ; i < 100 ; ++i )
    {
    mv[0] = ( normalGenerator1->GetVariate() * standardDeviation1 ) + mean1;
	mv[1] = ( normalGenerator2->GetVariate() * standardDeviation2 ) + mean2;
    sample->PushBack( mv );
    }
    normalGenerator1->Initialize( 7654 );
	normalGenerator2->Initialize( 5467 );
    mean1 = 300;
    mean2 = 300;
	standardDeviation1 = 60;
	standardDeviation2 = 60;
  for ( unsigned int i = 0 ; i < 100 ; ++i )
    {
    mv[0] = ( normalGenerator1->GetVariate() * standardDeviation1 ) + mean1;
	mv[1] = ( normalGenerator2->GetVariate() * standardDeviation2 ) + mean2;

    sample->PushBack( mv );
    }
  typedef itk::Array< double > ParametersType;
  ParametersType params( 4 );

  std::vector< ParametersType > initialParameters( numberOfClasses );
  params[0] = 100.0;      // Mean of Component 1, gaussian 1 - Estimate
  params[1] = 100.0;      // Mean of Component 2, gaussian 1 - Estimate
  params[2] = 800.0;      // Variance of Component 1, gaussian 1 - Estimate
  params[3] = 800.0;      // Variance of  Componet 2, gaussian 1 - Estimate

  initialParameters[0] = params;

  params[0] = 200.0;
  params[1] = 200.0;
  params[2] = 850.0;
  params[3] = 850.0;
  initialParameters[1] = params;

  params[0] = 300.0;
  params[1] = 300.0;
  params[2] = 1500.0;
  params[3] = 1500.0;

  initialParameters[2] = params;

  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( 1200 );
  

  itk::Array< double > initialProportions(numberOfClasses);
  initialProportions[0] = 0.4;
  initialProportions[1] = 0.4;
  initialProportions[2] = 0.2;	

  estimator->SetInitialProportions( initialProportions );
  
  for ( unsigned int i = 0 ; i < numberOfClasses ; i++)
    {
    estimator->AddComponent( (ComponentType::Superclass*)
                             (components[i]).GetPointer() );
    }
  
  estimator->Update();
  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;
}