[vtkusers] PCA in VTK

[Jonathan Bailleul]

#include "vtkPCAAnalysisFilter.h"

It is simple to use and features many things you might need...
It takes in input polydatas and can generate shape instances with
coefficients you provide, or estimate coefficients of an unknown
polydata.

(of course, polydata designate landmarked data).

Here is old code playing around with it (not reference code ok?), just
showing how to use some features...


So, you managed to implement the Frangi method, didn't you? ;)



#include "vtkRenderer.h"
#include "vtkRenderWindow.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkPolyData.h"
#include "vtkPoints.h"
#include "vtkCellArray.h"
#include "vtkFloatArray.h"
#include "vtkPolyDataMapper.h"
#include "vtkActor.h"
#include "vtkCamera.h"

#include "vtkPolyDataNormals.h"
#include "vtkPCAAnalysisFilter.h"
#include "vtkProcrustesAlignmentFilter.h"
#include "vtkPolyDataReader.h"
#include "vtkPolyDataWriter.h"
#include "vtkLandmarkTransform.h"

#include "vtkPointSetToPointSetFilter.h"
#include "vtkDataSetWriter.h"

#include <assert.h>


/* GLOS (Graphic Library in Open Source), an ANSI Common Lisp OpenGL
subset.
   Copyright (C) 2001 the GLOS development team
(http://glos.sourceforge.net) */



static void
usage(int argc, char **argv) 
{
  if (argc - 1 <= 2)
    {
      printf("Displays PCA mode of given already procrustes-aligned vtk
files\n");
      printf("Usage: %s <mode number> <output instances number> <input
vtk landmark files> \n", argv[0]);
      printf("\ncreate shape files deviated from mean shape along
selected mode. Instances number: e.g 1 means instances for phi = {-3.0,
0, -3.0}\n");
      
      exit(1);
      
    }  
}


 

/* mode: index of mode from 0*/
void 
CreateShapeInstance(vtkPCAAnalysisFilter *pca_filter, char* filename, 
		    int mode, float variance_fact, vtkPointSet *shape_container)
{
  vtkFloatArray *params = vtkFloatArray::New();
  params -> SetNumberOfComponents(1); // 1 per tuple
  params -> SetNumberOfTuples(1 + mode); // 1 instance 


  // float val[1]; val[0] = variance_fact; params -> SetTuple(0, val);
  for (int i = 0; i < mode; i++) {params -> SetTuple1(i, 0.0);}
  params -> SetTuple1(mode, variance_fact);
  
  pca_filter -> GetParameterisedShape(params, shape_container);

  vtkDataSetWriter *writer = vtkDataSetWriter::New();

  writer -> SetFileName(filename);
  writer -> SetInput(shape_container);
  writer -> Update();
}



int 
main( int argc, char *argv[] )
{
  usage(argc, argv);

  char* filename_in;
  char* filename_out;
  int i;
  int inputs = argc - 3;
  int mode_number = atoi(argv[1]);
  int mode_instances = atoi(argv[2]); // shape instances in range ]0;
+-3]


  //
----------------------------------------------------------------------
  /* STEP 1: PCA computation */
  //
----------------------------------------------------------------------


  // allocate one reader per file since input pca filter source
  // are read in parallel when filter is executed
  vtkPolyDataReader **readers;
  assert(readers = (vtkPolyDataReader**)calloc(inputs,
sizeof(vtkPolyDataReader*))); 
  for (i = 0; i < inputs; i++)
    readers[i] = vtkPolyDataReader::New();
  

  vtkPCAAnalysisFilter *pca = vtkPCAAnalysisFilter::New();
  pca -> SetNumberOfInputs(inputs);
  

  cout << "Inputfiles: " << inputs << endl;
  
  // assign pca inputs with landmark file readers  
  for (i = 0; i < inputs; i++)
    {
      filename_in = argv[i + 3];
      cout << "infile: " << filename_in << endl;
      
      readers[i] -> SetFileName(filename_in);
      readers[i] -> Update();
      pca -> SetInput(i, readers[i] -> GetOutput());
    }
  
  
  // execute pca
  pca -> Update();   
  pca -> PrintSelf(cout, 0);

  

  //
----------------------------------------------------------------------
  /* STEP 2 */
  //
----------------------------------------------------------------------

  
  // eigenmodes display  
  cout << endl << "eigenvalues:";

  inputs = pca -> GetEvals() -> GetNumberOfTuples();
  for (i = 0; i < inputs; i++)
    cout << " " << pca -> GetEvals() -> GetValue(i);
  cout << endl;
  

  float fact_delta = 3.0 / mode_instances;

  float fact;
  char* outfile;
  assert(outfile = (char *) calloc (150, sizeof(char)));


  sprintf(outfile, "./file_out_meanshape.vtk");//, mode_number, 0.0);
  CreateShapeInstance(pca, outfile, mode_number, 0.0, readers[0] ->
GetOutput());

  for (i = 0, fact = fact_delta; i < mode_instances; i++, fact +=
fact_delta)
    {
      sprintf(outfile, "./file_out_mode_%d_varpos_%d.vtk", mode_number,
i);
      CreateShapeInstance(pca, outfile, mode_number, fact, readers[0] ->
GetOutput());
    }
  
  for (i = 0, fact = -fact_delta; i < mode_instances; i++, fact -=
fact_delta)
    {
      sprintf(outfile, "./file_out_mode_%d_varneg_%d.vtk", mode_number,
i);
      CreateShapeInstance(pca, outfile, mode_number, fact, readers[0] ->
GetOutput());
    }


  cout << endl;
  return 0;
}




-- 
-----------------------------------
Jonathan BAILLEUL, Doctorant
GREYC Image - Université de Caen 
http://www.greyc.ismra.fr/~bailleul