[Insight-users] Help required on ITK Expectation Maximization
Algorithm - Aplogies for the resend
Prashanth Ravindran (BioImagene)
prashanth at bioimagene.com
Fri Dec 28 02:42:56 EST 2007
Sorry for the resend.. the last time I sent I noticed that only the code was visible on the mailing list and the question was not. The code is now inline
Hi insight-users-
I am trying to use ITK’s Expectation Maximization algorithm with measurement vectors that are 2-D (itk::Vector<double,2>). I have modified the 1-D example provided to support 2-D. When trying to initialize the algorithm I use the following initialization condition
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
// Isotropic gaussians assumed here
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 );
However when I run the code the output that I get is :
C:\InsightToolkit-3.4.0\Utilities\vxl\core\vnl/algo/vnl_svd.txx: suspicious return value (2) from SVDC
C:\InsightToolkit-3.4.0\Utilities\vxl\core\vnl/algo/vnl_svd.txx: M is 2x2
M = [ ...
0 0.0000000000000
0.0000000000000 0 ]
Cluster[0]
Parameters:
[1.25925e-302, 1.09379e-303]
Proportion: 1.#QNAN
Cluster[1]
Parameters:
[1.25925e-302, 1.09379e-303]
Proportion: 1.#QNAN
Cluster[2]
Parameters:
[1.25925e-302, 1.09379e-303]
Proportion: 1.#QNAN
My hunch is that the way I am setting the initial parameters is wrong and causing the algorithm to not converge. The whole source is attached to the mail. I would be very grateful for any help with this.
Thanks,
Prashanth Ravindran
BioImagene India
Code:
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;
}
More information about the Insight-users
mailing list