[Insight-users] cannot get proper classification using itk:GaussianDensityFunction after EM clustering

[Baoyun Li]

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);
           }
     
     }



      
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20090317/78abefc4/attachment-0001.htm>