Compute PCA Shape From Training Sample#

Note

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

Synopsis#

Compute a PCA shape model from a training sample.

Results#

Note

Help Wanted Implementation of Results for sphinx examples containing this message. Reconfiguration of CMakeList.txt may be necessary. Write An Example <https://itk.org/ITKExamples/Documentation/Contribute/WriteANewExample.html>

Code#

C++#

/*
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"

int
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 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->SetStartIndex(0);
  fileNamesCreator->SetEndIndex(nb_train - 1);
  fileNamesCreator->SetSeriesFormat(argv[2]);
  const std::vector<std::string> & shapeModeFileNames = fileNamesCreator->GetFileNames();

  for (unsigned int k = 0; k < nb_train; ++k)
  {
    trainingImages[k] = itk::ReadImage<ImageType>(shapeModeFileNames[k]);
  }

  using my_Estimatortype = itk::ImagePCAShapeModelEstimator<ImageType, ImageType>;
  auto filter = my_Estimatortype::New();
  filter->SetNumberOfTrainingImages(nb_train);
  filter->SetNumberOfPrincipalComponentsRequired(2);

  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->SetStartIndex(0);
  fileNamesOutCreator->SetEndIndex(nb_modes - 1);
  fileNamesOutCreator->SetSeriesFormat(argv[4]);
  const std::vector<std::string> & outFileNames = fileNamesOutCreator->GetFileNames();

  auto scaler = ScaleType::New();

  filter->Update();
  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;

    scaler->SetInput(filter->GetOutput(k));
    scaler->SetConstant(sv_n);

    itk::WriteImage(scaler->GetOutput(), outFileNames[k]);
  }

  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:

All ITK Sphinx Examples

Compute PCA Shape From Training Sample

See itk::ImagePCAShapeModelEstimator for additional documentation.