[Insight-users] PCAnalisys

[lucas]

Hi Luis,

I'm trying to perform a PCA on a set of training signed 
distance functions but although they are different from 
each other all the eigenvalues are zero.
I think I might be making a mistake when loading the 
training datasets but I'm not sure.
The following one is the code I'm using: 


#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "vnl/vnl_matrix_fixed.h"
#include "vnl/vnl_math.h"
#include "itkTextOutput.h"
#include "itkVector.h"
#include "itkImagePCAShapeModelEstimator.h"


int main( int argc, char * argv[] )
{
   if( argc < 4 )
     {
     std::cerr << "Usage: " << argv[0];
     std::cerr << std::endl;
     return 1;
     }

   typedef short InputPixelType;
   typedef double   OutputPixelType;

   typedef itk::Image<InputPixelType,2> InputImageType;
   typedef itk::Image<OutputPixelType,2> OutputImageType;

   typedef itk::ImageFileReader< InputImageType  > 
 ReaderType;
   typedef itk::ImageFileWriter< OutputImageType > 
 WriterType;

   ReaderType::Pointer reader = ReaderType::New();
   WriterType::Pointer writer = WriterType::New();

   //----------------------------------------------------------------------
   //----------------------------------------------------------------------
   //Set the image model estimator
   //----------------------------------------------------------------------
   typedef 
itk::ImagePCAShapeModelEstimator<InputImageType, 
OutputImageType>
     ImagePCAShapeModelEstimatorType;

   ImagePCAShapeModelEstimatorType::Pointer
     applyPCAShapeEstimator = 
ImagePCAShapeModelEstimatorType::New();

   //----------------------------------------------------------------------
   //Set the parameters of the clusterer
   //----------------------------------------------------------------------
   applyPCAShapeEstimator->SetNumberOfTrainingImages(4);
   applyPCAShapeEstimator->SetNumberOfPrincipalComponentsRequired(2);

   reader->SetFileName( argv[1] );
   reader->Update();
   applyPCAShapeEstimator->SetInput(0,reader->GetOutput() 
);

   reader->SetFileName( argv[2] );
   reader->Update();
   applyPCAShapeEstimator->SetInput(1,reader->GetOutput());

   reader->SetFileName( argv[3] );
   reader->Update();
   applyPCAShapeEstimator->SetInput(2, 
reader->GetOutput());

   reader->SetFileName( argv[4] );
   reader->Update();
   applyPCAShapeEstimator->SetInput(3, 
reader->GetOutput());

   applyPCAShapeEstimator->Update();

   //----------------------------------------------------------------------
   // Write the mean image
   //----------------------------------------------------------------------
   writer->SetFileName("mean.vtk");
   writer->SetInput(applyPCAShapeEstimator->GetOutput( 0 ) 
);
   try
     {
       writer->Update();
     }
   catch( itk::ExceptionObject exp )
     {
     std::cerr << "Exception caught ! mean writer" << 
std::endl;
     std::cerr <<     exp    << std::endl;
     }

    //----------------------------------------------------------------------
   // Write first mode image
   //----------------------------------------------------------------------
   writer->SetFileName("1mode.vtk");
   writer->SetInput(applyPCAShapeEstimator->GetOutput( 1 ) 
);
   try
     {
       writer->Update();
     }
   catch( itk::ExceptionObject exp )
     {
     std::cerr << "Exception caught ! 1st mode writer" << 
std::endl;
     std::cerr <<     exp    << std::endl;
     }

    //----------------------------------------------------------------------
   // Write second mode image
   //----------------------------------------------------------------------
   writer->SetFileName("2mode.vtk");
   writer->SetInput(applyPCAShapeEstimator->GetOutput( 2 ) 
);
   try
     {
       writer->Update();
     }
   catch( itk::ExceptionObject exp )
     {
     std::cerr << "Exception caught ! 2nd mode writer" << 
std::endl;
     std::cerr <<     exp    << std::endl;
     }

   vnl_vector<double> eigenValues =
     applyPCAShapeEstimator->GetEigenValues();

   for(unsigned int i= 0; i< 4 ; i++ )
     {
     std::cout << eigenValues[ i ] << std::endl;
     }

   return 0;
}

If you have any clue of what am I doing wrong please let 
me know.
Thanks,

Lucas Lorenzo