[Insight-users] ImagePCAShapeModelEstimator

[Reinhard Hameeteman]

Hi Siddharth,

Isn't your mean image equal to nr. 56?
If you ask me, this:

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

will set 56 pointers to the same address (the output of the reader),
which will contain the last read image. The result is you only have 1
principle component.

Reinhard

2006/6/17, Siddharth Vikal <siddharthvikal at gmail.com>:
> 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;
> >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<
>