<DIV id=RTEContent>Hi all, </DIV> <DIV> </DIV> <DIV>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! </DIV> <DIV>Esen</DIV> <DIV> </DIV> <DIV><STRONG>Results:</STRONG></DIV> <DIV>Cluster[0]<BR> Parameters:<BR> [0, 0]<BR>
Proportion: 1<BR>Cluster[1]<BR> Parameters:<BR> [0, 0]<BR> Proportion: 2.22156e-053<BR><STRONG>Code:</STRONG></DIV> <DIV> </DIV> <DIV>#include "itkScalarImageToHistogramGenerator.h"<BR>#include "itkImage.h"<BR>#include "itkVector.h"<BR>#include "itkListSample.h"<BR>#include "itkGaussianMixtureModelComponent.h"<BR>#include "itkExpectationMaximizationMixtureModelEstimator.h"<BR>// Software Guide : EndCodeSnippet</DIV> <DIV>#include "itkImageFileReader.h"</DIV> <DIV>int main( int argc, char * argv [] )<BR>{</DIV> <DIV>/*<BR> if( argc < 2 )<BR> {<BR> std::cerr << "Missing command line arguments" << std::endl;<BR> std::cerr << "Usage : ImageHistogram1 inputImageFileName " <<
std::endl;<BR> return -1;<BR> }<BR>*/</DIV> <DIV> typedef unsigned char PixelType;<BR> const unsigned int Dimension = 2;</DIV> <DIV> typedef itk::Image<PixelType, Dimension > ImageType;</DIV> <DIV> typedef itk::ImageFileReader< ImageType > ReaderType;</DIV> <DIV> ReaderType::Pointer reader = ReaderType::New();</DIV> <DIV> // reader->SetFileName( argv[1] );<BR> reader->SetFileName( "Simulatedchecker.png" );</DIV> <DIV> try<BR> {<BR> reader->Update();<BR> }<BR> catch( itk::ExceptionObject & excp )<BR> {<BR> std::cerr << "Problem encoutered while reading image file : " << argv[1] << std::endl;<BR> std::cerr << excp <<
std::endl;<BR> return -1;<BR> }</DIV> <DIV> typedef itk::Statistics::ScalarImageToHistogramGenerator< <BR> ImageType <BR> > HistogramGeneratorType;</DIV> <DIV> HistogramGeneratorType::Pointer histogramGenerator = HistogramGeneratorType::New();</DIV> <DIV>
histogramGenerator->SetInput( reader->GetOutput() );</DIV> <DIV> histogramGenerator->SetNumberOfBins( 255 );<BR> histogramGenerator->SetMarginalScale( 10.0 );<BR> histogramGenerator->Compute();</DIV> <DIV> typedef HistogramGeneratorType::HistogramType HistogramType;</DIV> <DIV> const HistogramType * histogram = histogramGenerator->GetOutput();</DIV> <DIV> const unsigned int histogramSize = histogram->Size();</DIV> <DIV> std::cout << "Histogram size " << histogramSize << std::endl;</DIV> <DIV> unsigned int bin;<BR> for( bin=0; bin < histogramSize; bin++ )<BR> {<BR> std::cout << "bin = " << bin << " frequency = ";<BR> std::cout << histogram->GetFrequency( bin, 0 ) << std::endl;<BR> }</DIV> <DIV> </DIV> <DIV> //*******************************<BR>
//=====end of histogram ========<BR> //======start EM ================</DIV> <DIV> unsigned int numberOfClasses = 2;<BR> typedef itk::Vector< double, 1 > MeasurementVectorType;<BR> typedef itk::Statistics::ListSample< MeasurementVectorType > SampleType;<BR> SampleType::Pointer sample = SampleType::New();<BR> <BR> for ( unsigned int i = 0 ; i < histogramSize ; ++i )<BR> {</DIV> <DIV> sample->PushBack( histogram->GetFrequency( bin, 0 ) );//0 is for first channel<BR> }</DIV> <DIV> //initial parameters for EM<BR> typedef itk::Array< double > ParametersType;<BR> ParametersType params( 2 );</DIV> <DIV> std::vector< ParametersType > initialParameters( numberOfClasses );<BR> params[0] = 50.0;<BR> params[1] = 100.0;<BR> initialParameters[0] = params;</DIV> <DIV> params[0] = 200.0;<BR> params[1] =
150.0;<BR> initialParameters[1] = params;<BR> <BR> // EM will estimate with Gaussians<BR> typedef itk::Statistics::GaussianMixtureModelComponent< SampleType > <BR> ComponentType;</DIV> <DIV> std::vector< ComponentType::Pointer > components;<BR> for ( unsigned int i = 0 ; i < numberOfClasses ; i++ )<BR> {<BR> components.push_back( ComponentType::New() );<BR> (components[i])->SetSample( sample );<BR> (components[i])->SetParameters( initialParameters[i] );<BR> }</DIV> <DIV> typedef itk::Statistics::ExpectationMaximizationMixtureModelEstimator< <BR> SampleType > EstimatorType;<BR> EstimatorType::Pointer estimator = EstimatorType::New();</DIV>
<DIV> estimator->SetSample( sample );<BR> estimator->SetMaximumIteration( 300 );</DIV> <DIV> itk::Array< double > initialProportions(numberOfClasses);<BR> initialProportions[0] = 0.5;<BR> initialProportions[1] = 0.5;</DIV> <DIV> estimator->SetInitialProportions( initialProportions );<BR> <BR> for ( unsigned int i = 0 ; i < numberOfClasses ; i++)<BR> {<BR> estimator->AddComponent( (ComponentType::Superclass*)<BR> (components[i]).GetPointer() );<BR> }<BR> <BR> try<BR> {<BR> estimator->Update();<BR> }<BR> catch( itk::ExceptionObject & excp )<BR> {<BR> std::cerr << "Problem
encoutered while estimating" << std::endl;<BR> std::cerr << excp << std::endl;<BR> return -1;<BR> }</DIV> <DIV><BR> int termincode = estimator->GetTerminationCode (); // CONVERGED = 0, NOT_CONVERGED = 1 }<BR> std::cout<<"termination code= "<<termincode<<std::endl;</DIV> <DIV><BR> for ( unsigned int i = 0 ; i < numberOfClasses ; i++ )<BR> {<BR> std::cout << "Cluster[" << i << "]" << std::endl;<BR> std::cout << " Parameters:" << std::endl;<BR> std::cout << " " << (components[i])->GetFullParameters() <BR> << std::endl;<BR> std::cout << " Proportion:
";<BR> std::cout << " " << (*estimator->GetProportions())[i] << std::endl;<BR> }</DIV> <DIV> return 0;<BR> <BR>}</DIV><p>
<hr size=1> <a href="http://us.lrd.yahoo.com/_ylc=X3oDMTFqODRtdXQ4BF9TAzMyOTc1MDIEX3MDOTY2ODgxNjkEcG9zAzEEc2VjA21haWwtZm9vdGVyBHNsawNmYw--/SIG=110oav78o/**http%3a//farechase.yahoo.com/">Yahoo! FareChase - Search multiple travel sites in one click.</a>