VTK/Examples/Cxx/Utilities/PCADemo

PCADemo.cxx

#include <vtkVersion.h>
#include "vtkSmartPointer.h"
#include "vtkLine.h"
#include "vtkLineSource.h"
#include "vtkTransform.h"
#include "vtkVertexGlyphFilter.h"
#include "vtkMath.h"
#include "vtkDoubleArray.h"
#include "vtkMultiBlockDataSet.h"
#include "vtkPCAStatistics.h"
#include "vtkStringArray.h"
#include "vtkTable.h"
#include "vtkTransformPolyDataFilter.h"
#include <vtkActor.h>
#include <vtkPoints.h>
#include <vtkPolyData.h>
#include <vtkPolyDataMapper.h>
#include <vtkProperty.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>
#include <vtkVertexGlyphFilter.h>

int main(int, char*[])
{
  vtkMath::RandomSeed(0);
  
  vtkSmartPointer<vtkPoints> points =
    vtkSmartPointer<vtkPoints>::New();
  for(unsigned int i = 0; i < 200; i++)
    {
    double x = vtkMath::Gaussian(0, 2);
    double y = vtkMath::Gaussian(0, 5);
    points->InsertNextPoint(x,y,0);
    }
  vtkSmartPointer<vtkPolyData> polydata =
    vtkSmartPointer<vtkPolyData>::New();
  polydata->SetPoints(points);

  vtkSmartPointer<vtkTransform> transform =
    vtkSmartPointer<vtkTransform>::New();
  transform->RotateZ(30);
  
  vtkSmartPointer<vtkTransformPolyDataFilter> transformFilter =
    vtkSmartPointer<vtkTransformPolyDataFilter>::New();
  transformFilter->SetTransform(transform);
#if VTK_MAJOR_VERSION <= 5
  transformFilter->SetInputConnection(polydata->GetProducerPort());
#else
  transformFilter->SetInputData(polydata);
#endif
  transformFilter->Update();

  // These would be all of your "x" values.
  vtkSmartPointer<vtkDoubleArray> xArray =
    vtkSmartPointer<vtkDoubleArray>::New();
  xArray->SetNumberOfComponents(1);
  xArray->SetName("x");
  
  // These would be all of your "y" values.
  vtkSmartPointer<vtkDoubleArray> yArray =
    vtkSmartPointer<vtkDoubleArray>::New();
  yArray->SetNumberOfComponents(1);
  yArray->SetName("y");

  for(vtkIdType i = 0; i < polydata->GetNumberOfPoints(); i++)
    {
    double p[3];
    transformFilter->GetOutput()->GetPoint(i,p);
    xArray->InsertNextValue(p[0]);
    yArray->InsertNextValue(p[1]);
    }
    
  vtkSmartPointer<vtkTable> datasetTable =
    vtkSmartPointer<vtkTable>::New();
  datasetTable->AddColumn(xArray);
  datasetTable->AddColumn(yArray);

  vtkSmartPointer<vtkPCAStatistics> pcaStatistics =
    vtkSmartPointer<vtkPCAStatistics>::New();
#if VTK_MAJOR_VERSION <= 5
  pcaStatistics->SetInput( vtkStatisticsAlgorithm::INPUT_DATA, datasetTable );
#else
  pcaStatistics->SetInputData( vtkStatisticsAlgorithm::INPUT_DATA, datasetTable );
#endif

  pcaStatistics->SetColumnStatus("x", 1 );
  pcaStatistics->SetColumnStatus("y", 1 );

  pcaStatistics->RequestSelectedColumns();
  pcaStatistics->SetDeriveOption(true);
  pcaStatistics->Update();

  ///////// Eigenvalues ////////////
  vtkSmartPointer<vtkDoubleArray> eigenvalues =
    vtkSmartPointer<vtkDoubleArray>::New();
  pcaStatistics->GetEigenvalues(eigenvalues);
  for(vtkIdType i = 0; i < eigenvalues->GetNumberOfTuples(); i++)
    {
    std::cout << "Eigenvalue " << i << " = " << eigenvalues->GetValue(i) << std::endl;
    }

  ///////// Eigenvectors ////////////
  vtkSmartPointer<vtkDoubleArray> eigenvectors =
    vtkSmartPointer<vtkDoubleArray>::New();

  pcaStatistics->GetEigenvectors(eigenvectors);

  vtkSmartPointer<vtkDoubleArray> evec1 =
    vtkSmartPointer<vtkDoubleArray>::New();
  pcaStatistics->GetEigenvector(0, evec1);

  vtkSmartPointer<vtkDoubleArray> evec2 =
    vtkSmartPointer<vtkDoubleArray>::New();
  pcaStatistics->GetEigenvector(1, evec2);

  double scale = 3.0;
  vtkSmartPointer<vtkLineSource> vector1Source =
    vtkSmartPointer<vtkLineSource>::New();
  vector1Source->SetPoint1(0,0,0);
  vector1Source->SetPoint1(scale * evec1->GetValue(0), scale* evec1->GetValue(1),0);

  vtkSmartPointer<vtkPolyDataMapper> vec1Mapper =
    vtkSmartPointer<vtkPolyDataMapper>::New();
  vec1Mapper->SetInputConnection(vector1Source->GetOutputPort());

  vtkSmartPointer<vtkActor> vector1Actor =
    vtkSmartPointer<vtkActor>::New();
  vector1Actor->SetMapper(vec1Mapper);
  vector1Actor->GetProperty()->SetColor(0,1,0);
  vector1Actor->GetProperty()->SetLineWidth(3);
  
  vtkSmartPointer<vtkLineSource> vector2Source =
    vtkSmartPointer<vtkLineSource>::New();
  vector2Source->SetPoint1(0,0,0);
  vector2Source->SetPoint1(scale* evec2->GetValue(0), scale*evec2->GetValue(1),0);


  vtkSmartPointer<vtkPolyDataMapper> vec2Mapper =
    vtkSmartPointer<vtkPolyDataMapper>::New();
  vec2Mapper->SetInputConnection(vector2Source->GetOutputPort());

  vtkSmartPointer<vtkActor> vector2Actor =
    vtkSmartPointer<vtkActor>::New();
  vector2Actor->SetMapper(vec2Mapper);
  vector2Actor->GetProperty()->SetColor(1,0,0);
  vector2Actor->GetProperty()->SetLineWidth(3);

  // Project all of the points onto the eigenvector with
  // the largest eigenvalues
  double p0[2];
  p0[0] = -100*evec1->GetValue(0);
  p0[1] = -100*evec1->GetValue(1);
  double p1[2];
  p1[0] = 100*evec1->GetValue(0);
  p1[1] = 100*evec1->GetValue(1);

  vtkSmartPointer<vtkPoints> projectedPoints =
    vtkSmartPointer<vtkPoints>::New();
  for(vtkIdType i = 0; i < polydata->GetNumberOfPoints(); i++)
    {
    double p[3];
    transformFilter->GetOutput()->GetPoint(i,p);
    double t;
    double closestPoint[3];
    vtkLine::DistanceToLine(p, p0, p1, t, closestPoint);
    double newP[3];
    double v[3];
    vtkMath::Subtract(p1, p0, v);
    vtkMath::Normalize(v);
    vtkMath::MultiplyScalar(v, t);
    vtkMath::Add(p0, v, newP);
    projectedPoints->InsertNextPoint(t, 0, 0);
    }
    
  vtkSmartPointer<vtkPolyData> projectedPolyData =
    vtkSmartPointer<vtkPolyData>::New();
  projectedPolyData->SetPoints(projectedPoints);

  vtkSmartPointer<vtkVertexGlyphFilter> projectedGlyphFilter =
    vtkSmartPointer<vtkVertexGlyphFilter>::New();
#if VTK_MAJOR_VERSION <= 5
  projectedGlyphFilter->SetInputConnection(projectedPolyData->GetProducerPort());
#else
  projectedGlyphFilter->SetInputData(projectedPolyData);
#endif
  projectedGlyphFilter->Update();

  vtkSmartPointer<vtkPolyDataMapper> projectedMapper =
    vtkSmartPointer<vtkPolyDataMapper>::New();
  projectedMapper->SetInputConnection(projectedGlyphFilter->GetOutputPort());
  vtkSmartPointer<vtkActor> projectedActor =
    vtkSmartPointer<vtkActor>::New();
  projectedActor->SetMapper(projectedMapper);
  projectedActor->GetProperty()->SetPointSize(3);
  
  // There will be one render window
  vtkSmartPointer<vtkRenderWindow> renderWindow =
    vtkSmartPointer<vtkRenderWindow>::New();
  renderWindow->SetSize(600, 300);

  // And one interactor
  vtkSmartPointer<vtkRenderWindowInteractor> interactor =
    vtkSmartPointer<vtkRenderWindowInteractor>::New();
  interactor->SetRenderWindow(renderWindow);

  // Define viewport ranges
  // (xmin, ymin, xmax, ymax)
  double leftViewport[4] = {0.0, 0.0, 0.5, 1.0};
  double rightViewport[4] = {0.5, 0.0, 1.0, 1.0};

  // Setup both renderers
  vtkSmartPointer<vtkRenderer> leftRenderer =
    vtkSmartPointer<vtkRenderer>::New();
  renderWindow->AddRenderer(leftRenderer);
  leftRenderer->SetViewport(leftViewport);
  leftRenderer->SetBackground(.6, .5, .4);

  vtkSmartPointer<vtkRenderer> rightRenderer =
    vtkSmartPointer<vtkRenderer>::New();
  renderWindow->AddRenderer(rightRenderer);
  rightRenderer->SetViewport(rightViewport);
  rightRenderer->SetBackground(.4, .5, .6);

  vtkSmartPointer<vtkVertexGlyphFilter> glyphFilter =
    vtkSmartPointer<vtkVertexGlyphFilter>::New();
  glyphFilter->SetInputConnection(transformFilter->GetOutputPort());
  glyphFilter->Update();
  
  vtkSmartPointer<vtkPolyDataMapper> originalMapper =
    vtkSmartPointer<vtkPolyDataMapper>::New();
  originalMapper->SetInputConnection(glyphFilter->GetOutputPort());
  vtkSmartPointer<vtkActor> originalActor =
    vtkSmartPointer<vtkActor>::New();
  originalActor->SetMapper(originalMapper);
  originalActor->GetProperty()->SetPointSize(3);
  
  leftRenderer->AddActor(originalActor);
  leftRenderer->AddActor(vector1Actor);
  leftRenderer->AddActor(vector2Actor);

  rightRenderer->AddActor(projectedActor);

  leftRenderer->ResetCamera();
  rightRenderer->ResetCamera();

  renderWindow->Render();
  interactor->Start();

  return EXIT_SUCCESS;
}

Please try the new VTKExamples website.

CMakeLists.txt

cmake_minimum_required(VERSION 2.8)

PROJECT(PCADemo)

find_package(VTK REQUIRED)
include(${VTK_USE_FILE})

add_executable(PCADemo MACOSX_BUNDLE PCADemo.cxx)

if(VTK_LIBRARIES)
  target_link_libraries(PCADemo ${VTK_LIBRARIES})
else()
  target_link_libraries(PCADemo vtkHybrid vtkWidgets)
endif()

Download and Build PCADemo

Click here to download PCADemo. and its CMakeLists.txt file.

Once the tarball PCADemo.tar has been downloaded and extracted,

cd PCADemo/build 

• If VTK is installed:

cmake ..

• If VTK is not installed but compiled on your system, you will need to specify the path to your VTK build:

cmake -DVTK_DIR:PATH=/home/me/vtk_build ..

Build the project:

make

and run it:

./PCADemo

WINDOWS USERS PLEASE NOTE: Be sure to add the VTK bin directory to your path. This will resolve the VTK dll's at run time.