[Insight-users] ImagePCAShapeModelEstimator

[Siddharth Vikal]

Hello Zach and ITK,

I'm traying to generate a PCA shape model from a set of binary
training images, 56 in number, with noticeable variation in shapes
from nearly circular to super-elliptical type shapes.

I'm able to get the mean shape correctly(at visibly it looks so), the
problem is that all the principal components are same. There is no
variation in 1st principal component to last principal component,
though it should be since the training shapes are varying. Initially I
thought that I was not looking enough principal components, but later
I set number of principal component images to be computed to number of
training images. Still, all the principal component images are same.

Please find below inline the code that I have written and used. I'm
using ParaView to view the output results. Can someone please tell me
what is it that I'm doing wrong or what more do I need to do to get it
right?

waiting for responses.

thanks

regards
siddharth

>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  const    unsigned int    Dimension = 2;
  typedef  float    PixelType;


  //define input training image type and the output distance map type image
  typedef itk::Image <PixelType, Dimension> InternalImageType;


  // image reader
  typedef itk::ImageFileReader< InternalImageType  > ImageReaderType;
  ImageReaderType::Pointer  ImageReader  = ImageReaderType::New();

  //define the distance map generator filter
  typedef itk::SignedDanielssonDistanceMapImageFilter<InternalImageType,
InternalImageType> DistanceMapFilterType;
  DistanceMapFilterType::Pointer filterDaniel = DistanceMapFilterType::New();


  //define the PCA Shape Model estimator
  typedef itk::ImagePCAShapeModelEstimator<InternalImageType,
InternalImageType> ShapeEstimatorType;
  ShapeEstimatorType::Pointer shapeEstimator = ShapeEstimatorType::New();
  shapeEstimator->SetNumberOfTrainingImages(56);
  shapeEstimator->SetNumberOfPrincipalComponentsRequired(56);


//first read the training images
  const std::string FilePrefix = "D:\\TestData\\AlignedShapes\\Shape_";
  char *Suffix = new char[20];
  std::string Filename;

  int ImageNumber = 1;

  for (ImageNumber = 1; ImageNumber <=56; ImageNumber++)
  {	
    //form the file name
	sprintf(Suffix, "%03d.tiff", ImageNumber);
	Filename = FilePrefix + Suffix;
	ImageReader->SetFileName( Filename.c_str() );
	filterDaniel->SetInput(ImageReader->GetOutput());
	filterDaniel->Update();
       	shapeEstimator->SetInput(ImageNumber-1, filterDaniel->GetOutput());
  }

//run the shape estimator
  try
  {
	  shapeEstimator->Update();
  }
  catch( itk::ExceptionObject & err )
  {
    std::cerr << "ExceptionObject caught !" << std::endl;
    std::cerr << err << std::endl;
    return -1;
  }

//get the results
//mean image
  InternalImageType::Pointer MeanImage = InternalImageType::New();
  MeanImage = shapeEstimator->GetOutput(0);

  writer->SetFileName( "D:\\MeanShape.mha" );
  writer->SetInput( MeanImage);
  writer->Update();

// principal component images
// first collect eigen values to get standard deviations

  vnl_vector<double> eigenValues(56);

  eigenValues = shapeEstimator->GetEigenValues();

  typedef itk::ImageFileWriter< InternalImageType >  WriterType;


  typedef itk::ShiftScaleImageFilter<InternalImageType,
InternalImageType> DivideFilterType;


  WriterType::Pointer      writer =  WriterType::New();
  DivideFilterType::Pointer divider = DivideFilterType::New();


  InternalImageType::Pointer EigenImage = InternalImageType::New();


  const std::string WriteFilePrefix = "D:\\EigenShape_";
  for (ImageNumber = 1; ImageNumber <=56; ImageNumber++)
  {	
    //form the file name
	sprintf(Suffix, "%02d.mha", ImageNumber);
	Filename = WriteFilePrefix + Suffix;

	EigenImage = shapeEstimator->GetOutput(ImageNumber);
	
	writer->SetFileName( Filename.c_str());
	divider->SetInput(EigenImage);
//divide by standard deviation for normalization
	divider->SetScale(1/sqrt(eigenValues[ImageNumber-1]));
	writer->SetInput( divider->GetOutput());
	writer->Update();
	
  }


  return 0;
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<