[Insight-users] problem solved--- cannot get proper classification using itk:GaussianDensityFunction after EM clustering
Luis Ibanez
luis.ibanez at kitware.com
Tue Mar 17 18:29:13 EDT 2009
Hi Baoyun,
Thanks for letting us know that you found the solution to the problem.
Indeed, when you are using the Gaussian mixture modeling estimator,
the classification stage you should do it by using the
itkSampleClassifier class.
The code that you posted in your previous email, is basically what
the classifier contains inside, but with the added flexibility of
being reusable and allowing for different types of decision rules.
Looking at your previous code in its current form, it is hard to
tell what was going wrong.
However, some suspicious elements are:
A) OutputPixelType classlabel;
is uninitialized by the time you start iterating over pixels.
B) The actual classification code was conditioned to the
value of "acumulate1". Why was such thing needed ?
A simple way to find the problem in your code, will be to recompile
it with a higher settings of warnings (e.g. -Wall if you use gcc).
Chances are that the culprit of the problem is already begin flagged
out by the compiler.
BTW: You may find interesting to start using the new Statistics
Framework, that is available at:
http://www.na-mic.org/svn/NAMICSandBox/trunk/ITKStatisticsPipelineRefactoring/
You will find a description of what was changed (and why)
in the following wiki pages:
http://www.itk.org/Wiki/Proposals:Refactoring_Statistics_Framework_2007
Regards,
Luis
------------------
Baoyun Li wrote:
> Dear All:
>
> I just figure out the problem. I need to use
> itk::GaussianMixtureModelComponent->Evaluate to make decision. The
> result is ok.
>
> but i am still not clear what happened if I used
> itk::GaussianDensityFunction->Evaluate() to get decision. Maybe I did
> not use it properly. Still interested to know what happened.
>
> Can somebody tell me whether my code is wrong or I did not used
> itk::GaussianDensityFunction->Evaluate() properly?
>
> Thanks
>
> Baoyun
>
> ------------------------------------------------------------------------
> *From:* Baoyun Li <baoyun_li123 at yahoo.com>
> *To:* Luis Ibanez <luis.ibanez at kitware.com>
> *Cc:* insight-users at itk.org
> *Sent:* Tuesday, March 17, 2009 12:24:04 PM
> *Subject:* cannot get proper classification using
> itk:GaussianDensityFunction after EM clustering
>
> Dear ALL:
>
> I am using EM gassian mixutrue estimator to estimate the mean and
> covariance of vector image. I first adapt the image to sample, and the
> EM works well.
> It produced mean value for Gray matter , whithe matter and CSF, the mean
> value looks ok to me.
>
> But when I try to assign the label to the image, I cannot get the proper
> labled image. I got the estimated all the voxle to the one class. And
> the pdf of GassianDensityFunction always have the biggest value for CSF,
> that should not happend.
>
> Below is part of my code.
>
> 1, I first get the pararameter of each component by :
> finalparameters=(components[i])->GetFullParameters();
> 2, Then I decompose the parameters to estimatedmean and estimatedcovariance
> 3, set the mean for GaussianDensityFunction by :
> membershipFunctions[i]->SetMean(&estimatedmean );
> membershipFunctions[i]->SetCovariance( &estimatedcovariance);
>
> 4, In the lableing program, I used
> eveluateresult[i]=(membershipFunctions[i]->Evaluate(mv)) to get pdf for
> each class, mv is the measurement vector obtained of images by itk iterator
> 5, then I find the class with biggest pdf, assign the label .
>
> The code is listed below, can some body help me to find what is wrong.
>
>
> Thanks
>
> Baoyun
>
>
> ///////give the mean value and covariance to GaussianDensityFunction
> typedef
> itk::Statistics::GaussianDensityFunction<ComponentType::MeasurementVectorType
> > MembershipFunctionType;
> std::vector< MembershipFunctionType::Pointer >
> membershipFunctions(numberOfClasses);
> unsigned int paramIndex;
> ParametersType
> finalparameters((NumberOfComponents+1)*NumberOfComponents);
> typedef MembershipFunctionType::MeanType MeanType;
> typedef MembershipFunctionType::CovarianceType CovarianceType;
> MeanType estimatedmean(NumberOfComponents);
> CovarianceType estimatedcovariance;
> typedef ComponentType::MeasurementVectorSizeType
> MeasurementVectorSizeType;
> MeasurementVectorSizeType measurementVectorSize=NumberOfComponents;
> estimatedcovariance.SetSize( measurementVectorSize,
> measurementVectorSize );
> //estimatedmean.SetSize( measurementVectorSize, 1);
> // estimatedmean.Fill(0.0) ;
> for ( unsigned int i = 0 ; i < numberOfClasses ; i++ )
> {
> // membershipFunctions.push_back(MembershipFunctionType::New());
> membershipFunctions[i]=MembershipFunctionType::New();
> std::cout << " Size of Gaussian : ";
> std::cout << " " <<
> (membershipFunctions[i])->GetMeasurementVectorSize() << std::endl;
>
> finalparameters=(components[i])->GetFullParameters();
> paramIndex=0;
> for (unsigned int i1 = 0 ; i1 < measurementVectorSize ; i1++)
> {
>
> estimatedmean[i1] = finalparameters[paramIndex] ;
> ++paramIndex ;
> }
>
> for (unsigned int i1 = 0 ; i1 < measurementVectorSize ; i1++ )
> {
> for (unsigned int j1 = 0 ; j1 < measurementVectorSize; j1++ )
> {
>
> estimatedcovariance.GetVnlMatrix().put(i1, j1,
> finalparameters[paramIndex]) ;
> ++paramIndex ;
> }
> }
>
> // finalProportions[i]=(*estimator->GetProportions())[i];
> std::cout << "********************************" << std::endl;
> std::cout << "class[" << i << "]" << std::endl;
> std::cout << " Estimated Mean:" << std::endl;
> std::cout << " " << estimatedmean
> << std::endl;
> std::cout << " Estimated Covariance ";
> std::cout << " " << estimatedcovariance << std::endl;
>
>
> membershipFunctions[i]->SetMean(&estimatedmean );
> membershipFunctions[i]->SetCovariance( &estimatedcovariance);
> std::cout << "***************************" << std::endl;
> std::cout << "class[" << i << "]" << std::endl;
> std::cout << " Gaussain Mean:" << std::endl;
> std::cout << " "
> <<*(membershipFunctions[i]->GetMean()+0)<< std::endl;
> // std::cout << " "
> <<*(membershipFunctions[i]->GetMean()+1)<< std::endl;
> std::cout << " Gaussian Covariance ";
> std::cout << " " <<
> *(membershipFunctions[i]->GetCovariance()+0) << std::endl;
> // std::cout << " " <<
> *(membershipFunctions[i]->GetCovariance()+1) << std::endl;
> // std::cout << " " <<
> *(membershipFunctions[i]->GetCovariance()+3) << std::endl;
> }
> /////to lable the image
> typedef unsigned char OutputPixelType;
> typedef itk::Image< OutputPixelType, 3 > OutputImageType;
>
> OutputImageType::Pointer outputPtr=reader1->GetOutput();
> typedef itk::ImageRegionIterator< OutputImageType > ImageIterator;
> typedef VectorImageType::RegionType RegionType;
> RegionType region = outputPtr->GetBufferedRegion();
> //
> //
> ImageIterator pixel( outputPtr, region );
> typedef itk::Vector< double, numberOfClasses > EvaluationVectorType;
> EvaluationVectorType eveluateresult;
>
>
> //pixel.GoToBegin();
> // itk::Array< double > finalProportions(numberOfClasses);
> for(unsigned char i=0;i<numberOfClasses;i++)
> {
> finalProportions[i]=(*estimator->GetProportions())[i];
> }
>
> double maxvalue;
> OutputPixelType classlabel;
> for ( pixel.GoToBegin(),constIterator.GoToBegin(); !pixel.IsAtEnd();
> ++pixel,++constIterator )
> {
> mv=constIterator.Get();
> accumulate1=0;
> maxvalue=0;
> for( unsigned int i=0;i<NumberOfComponents;i++)
> {
> accumulate1=accumulate1+mv[i];
> }
> if(accumulate1>20)
> {
>
> for(unsigned char i=0;i<numberOfClasses;i++)
> {
>
> eveluateresult[i]=(membershipFunctions[i]->Evaluate(mv))*finalProportions[i];
> }
>
> for(unsigned char i=0;i<numberOfClasses;i++)
> {
> if(eveluateresult[i]>=maxvalue)
> {
> maxvalue=eveluateresult[i];
> classlabel=i+1;
> }
> // if(mv[0]==94&&mv[1]==77)
> // {
> // std::cout << "evaluation result" <<
> eveluateresult << "]" << std::endl;
> // }
> }
> pixel.Set(classlabel);
> }
> else
> {
> pixel.Set(0);
> }
>
> }
>
>
More information about the Insight-users
mailing list