Compute PCA Shape From Training Sample


Wish List Still needs additional work to finish proper creation of example.


Compute a PCA shape model from a training sample.



Help Wanted Implementation of Results for sphinx examples containing this message. Reconfiguration of CMakeList.txt may be necessary. Write An Example <>



Author: Juan Cardelino <juan dot cardelino at gmail dot com>

// ITK
#include "itkBinaryThresholdImageFilter.h"
#include "itkBoundedReciprocalImageFilter.h"
#include "itkChangeInformationImageFilter.h"
#include "itkCommand.h"
#include "itkCurvatureAnisotropicDiffusionImageFilter.h"
#include "itkEuler2DTransform.h"
#include "itkFastMarchingImageFilter.h"
#include "itkGeodesicActiveContourShapePriorLevelSetImageFilter.h"
#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkImagePCAShapeModelEstimator.h"
#include "itkMultiplyImageFilter.h"
#include "itkNumericSeriesFileNames.h"
#include "itkNormalVariateGenerator.h"
#include "itkOnePlusOneEvolutionaryOptimizer.h"
#include "itkPCAShapeSignedDistanceFunction.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkShapePriorMAPCostFunction.h"
#include "itkSpatialFunctionImageEvaluatorFilter.h"

// VNL
#include "vnl/vnl_sample.h"

main(int argc, char * argv[])
  if (argc < 5)
    std::cerr << "Missing Parameters " << std::endl;
    std::cerr << "Usage: " << argv[0];
    std::cerr << " nbTrain  trainFilePattern";
    std::cerr << " nbModes  modeFilePattern";
    std::cerr << std::endl;
    return 1;

  for (int i = 0; i < argc; i++)
    std::cout << "id: " << i << " arg: " << argv[i] << std::endl;
  constexpr unsigned int Dimension = 2;
  using my_PixelType = float;
  using ImageType = itk::Image<my_PixelType, Dimension>;

  using ReaderType = itk::ImageFileReader<ImageType>;
  using WriterType = itk::ImageFileWriter<ImageType>;
  using ScaleType = itk::MultiplyImageFilter<ImageType, ImageType, ImageType>;

  unsigned int nb_train = atoi(argv[1]);

  itk::NumericSeriesFileNames::Pointer fileNamesCreator = itk::NumericSeriesFileNames::New();
  std::vector<ImageType::Pointer>      trainingImages(nb_train);

  fileNamesCreator->SetEndIndex(nb_train - 1);
  const std::vector<std::string> & shapeModeFileNames = fileNamesCreator->GetFileNames();

  for (unsigned int k = 0; k < nb_train; k++)
    ReaderType::Pointer reader = ReaderType::New();
    trainingImages[k] = reader->GetOutput();

  using my_Estimatortype = itk::ImagePCAShapeModelEstimator<ImageType, ImageType>;
  my_Estimatortype::Pointer filter = my_Estimatortype::New();

  for (unsigned int k = 0; k < nb_train; k++)
    filter->SetInput(k, trainingImages[k]);

  unsigned int nb_modes = atoi(argv[3]);

  itk::NumericSeriesFileNames::Pointer fileNamesOutCreator = itk::NumericSeriesFileNames::New();

  fileNamesOutCreator->SetEndIndex(nb_modes - 1);
  const std::vector<std::string> & outFileNames = fileNamesOutCreator->GetFileNames();

  ScaleType::Pointer scaler = ScaleType::New();

  my_Estimatortype::VectorOfDoubleType v = filter->GetEigenValues();
  double                               sv_mean = sqrt(v[0]);

  for (unsigned int k = 0; k < nb_modes; k++)
    double sv = sqrt(v[k]);
    double sv_n = sv / sv_mean;
    // double sv_n=sv;
    std::cout << "writing: " << outFileNames[k] << std::endl;

    std::cout << "svd[" << k << "]: " << sv << " norm: " << sv_n << std::endl;
    WriterType::Pointer writer = WriterType::New();

  return EXIT_SUCCESS;

Classes demonstrated

template<typename TInputImage, typename TOutputImage = Image<double, TInputImage::ImageDimension>>
class ImagePCAShapeModelEstimator : public itk::ImageShapeModelEstimatorBase<TInputImage, TOutputImage>

Base class for ImagePCAShapeModelEstimator object.

itkImagePCAShapeModelEstimator performs a principal component analysis (PCA) on a set of images. The user specifies the number of training images and also the number of desired largest principal components needed. The ITK pipeline mechanism sets up the storage for both input and output images. The number of output images are the user specified number of desired largest principal components plus 1 (for the mean image).

The algorithm uses the VNL library to perform the eigen analysis. To speed the computation of the instead of performing the eigen analysis of the covariance vector A*A’ where A is a matrix with p x t, p = number of pixels or voxels in each images and t = number of training images, we calculate the eigen vectors of the inner product matrix A’*A. The resulting eigen vectors (E) are then multiplied with the the matrix A to get the principal components. The covariance matrix has a dimension of p x p. Since number of pixels in any image being typically very high the eigen decomposition becomes computationally expensive. The inner product on the other hand has the dimension of t x t, where t is typically much smaller that p. Hence the eigen decomposition (most compute intensive part) is an orders of magnitude faster.

The Update() function enables the calculation of the various models, creates the membership function objects and populates them.

ITK Sphinx Examples:

See itk::ImagePCAShapeModelEstimator for additional documentation.