PCADemo
VTKExamples/Cxx/Utilities/PCADemo
Code¶
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; }
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) target_link_libraries(PCADemo ${VTK_LIBRARIES})
Download and Build PCADemo¶
Danger
The generation of tar files has not been ported to the new VTKExamples. Some tarballs may be missing or out-of-date.
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.