Commit e22ff12e authored by Will Schroeder's avatar Will Schroeder

Created generalized point interpolation filter

This is a general framework to interpolate attributes from a data source,
using the input as the set of interpolation points. It requires
specifying an interpolation kernel (several of which have been
created), as well as a locator for fast spatial searching. Although
similar in concept to vtkProbeFilter, it does not interpolate
using a cells interpolation basis, rather it uses a local evaluation
of the specified kernel, applied to the local points (neighborhood)
around each interpolation point, to compute the interpolated data values.
The algorithm is threaded with vtkSMPTools.
parent 7416f8dc
Pipeline #7879 passed with stage
......@@ -36,6 +36,7 @@ vtk_add_test_cxx(${vtk-module}CxxTests tests
TestTreeBFSIterator.cxx
TestTreeDFSIterator.cxx
TestTriangle.cxx
TimePointLocators.cxx
otherCellArray.cxx
otherCellBoundaries.cxx
otherCellPosition.cxx
......
/*=========================================================================
Program: Visualization Toolkit
Module: TestPointLocators.cxx
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notice for more information.
=========================================================================*/
#include "vtkPointLocator.h"
#include "vtkStaticPointLocator.h"
#include "vtkKdTree.h"
#include "vtkKdTreePointLocator.h"
#include "vtkOctreePointLocator.h"
#include "vtkPoints.h"
#include "vtkPolyData.h"
#include "vtkTimerLog.h"
#include "vtkMath.h"
int TimePointLocators(int , char *[])
{
int nPts = 1000000;
int nQ = nPts/10;
int N = 10;
double R = 0.01;
vtkTimerLog* timer = vtkTimerLog::New();
double buildTime[4], cpTime[4], cnpTime[4], crpTime[4];
for (int i=0; i<4; ++i)
{
buildTime[i] = cpTime[i] = cnpTime[i] = crpTime[i] = 0.0;
}
cout << "\nTiming for " << nPts << " points, " << nQ << " queries\n";
// Populate a list of points and query locations
vtkPoints* points = vtkPoints::New();
points->SetDataTypeToDouble();
points->SetNumberOfPoints(nPts);
for (int i=0; i<nPts; ++i)
{
points->SetPoint(i, vtkMath::Random(-1,1), vtkMath::Random(-1,1),
vtkMath::Random(-1,1));
}
vtkPolyData* polydata = vtkPolyData::New();
polydata->SetPoints(points);
points->ComputeBounds();
// Create points array which are positions to probe data with FindClosestPoint().
vtkPoints* qPoints = vtkPoints::New();
qPoints->SetDataTypeToDouble();
qPoints->SetNumberOfPoints(nQ);
vtkMath::RandomSeed(314159);
for (int i=0; i<nQ; ++i)
{
qPoints->SetPoint(i, vtkMath::Random(-1,1), vtkMath::Random(-1,1),
vtkMath::Random(-1,1));
}
vtkIdList *closest = vtkIdList::New();
//---------------------------------------------------------------------------
// The simple uniform binning point locator
vtkPointLocator* uniformLocator = vtkPointLocator::New();
timer->StartTimer();
uniformLocator->SetDataSet(polydata);
uniformLocator->BuildLocator();
timer->StopTimer();
uniformLocator->Delete();
buildTime[0] = timer->GetElapsedTime();
uniformLocator = vtkPointLocator::New();
uniformLocator->SetDataSet(polydata);
uniformLocator->BuildLocator();
timer->StartTimer();
for (int i=0; i<nQ; ++i)
{
uniformLocator->FindClosestPoint(qPoints->GetPoint(i));
}
timer->StopTimer();
cpTime[0] = timer->GetElapsedTime();
timer->StartTimer();
for (int i=0; i<nQ; ++i)
{
uniformLocator->FindClosestNPoints(N, qPoints->GetPoint(i), closest);
}
timer->StopTimer();
cnpTime[0] = timer->GetElapsedTime();
timer->StartTimer();
for (int i=0; i<nQ; ++i)
{
uniformLocator->FindPointsWithinRadius(R, qPoints->GetPoint(i), closest);
}
timer->StopTimer();
crpTime[0] = timer->GetElapsedTime();
uniformLocator->Delete();
//---------------------------------------------------------------------------
// The simple uniform binning point locator, this time built
// statically and threaded (may be much faster on threaded machine).
vtkStaticPointLocator* staticLocator = vtkStaticPointLocator::New();
timer->StartTimer();
staticLocator->SetDataSet(polydata);
staticLocator->BuildLocator();
timer->StopTimer();
staticLocator->Delete();
buildTime[1] = timer->GetElapsedTime();
staticLocator = vtkStaticPointLocator::New();
staticLocator->SetDataSet(polydata);
staticLocator->BuildLocator();
timer->StartTimer();
for (int i=0; i<nQ; ++i)
{
staticLocator->FindClosestPoint(qPoints->GetPoint(i));
}
timer->StopTimer();
cpTime[1] = timer->GetElapsedTime();
timer->StartTimer();
for (int i=0; i<nQ; ++i)
{
staticLocator->FindClosestNPoints(N, qPoints->GetPoint(i), closest);
}
timer->StopTimer();
cnpTime[1] = timer->GetElapsedTime();
timer->StartTimer();
for (int i=0; i<nQ; ++i)
{
staticLocator->FindPointsWithinRadius(R, qPoints->GetPoint(i), closest);
}
timer->StopTimer();
crpTime[1] = timer->GetElapsedTime();
staticLocator->Delete();
//---------------------------------------------------------------------------
// The KD Tree point locator
vtkKdTreePointLocator* kdTreeLocator = vtkKdTreePointLocator::New();
timer->StartTimer();
kdTreeLocator->SetDataSet(polydata);
kdTreeLocator->BuildLocator();
kdTreeLocator->Delete();
timer->StopTimer();
buildTime[2] = timer->GetElapsedTime();
kdTreeLocator = vtkKdTreePointLocator::New();
kdTreeLocator->SetDataSet(polydata);
kdTreeLocator->BuildLocator();
timer->StartTimer();
for (int i=0; i<nQ; ++i)
{
kdTreeLocator->FindClosestPoint(qPoints->GetPoint(i));
}
timer->StopTimer();
cpTime[2] = timer->GetElapsedTime();
timer->StartTimer();
for (int i=0; i<nQ; ++i)
{
kdTreeLocator->FindClosestNPoints(N, qPoints->GetPoint(i), closest);
}
timer->StopTimer();
cnpTime[2] = timer->GetElapsedTime();
timer->StartTimer();
for (int i=0; i<nQ; ++i)
{
kdTreeLocator->FindPointsWithinRadius(R, qPoints->GetPoint(i), closest);
}
timer->StopTimer();
crpTime[2] = timer->GetElapsedTime();
kdTreeLocator->Delete();
//---------------------------------------------------------------------------
// Octree point locator
vtkOctreePointLocator* octreeLocator = vtkOctreePointLocator::New();
timer->StartTimer();
octreeLocator->SetDataSet(polydata);
octreeLocator->BuildLocator();
octreeLocator->Delete();
timer->StopTimer();
buildTime[3] = timer->GetElapsedTime();
octreeLocator = vtkOctreePointLocator::New();
octreeLocator->SetDataSet(polydata);
octreeLocator->BuildLocator();
timer->StartTimer();
for (int i=0; i<nQ; ++i)
{
octreeLocator->FindClosestPoint(qPoints->GetPoint(i));
}
timer->StopTimer();
cpTime[3] = timer->GetElapsedTime();
timer->StartTimer();
for (int i=0; i<nQ; ++i)
{
octreeLocator->FindClosestNPoints(N, qPoints->GetPoint(i), closest);
}
timer->StopTimer();
cnpTime[3] = timer->GetElapsedTime();
timer->StartTimer();
for (int i=0; i<nQ; ++i)
{
octreeLocator->FindPointsWithinRadius(R, qPoints->GetPoint(i), closest);
}
timer->StopTimer();
crpTime[3] = timer->GetElapsedTime();
octreeLocator->Delete();
//---------------------------------------------------------------------------
// Print out the statistics
cout << "Build and delete tree\n";
cout << "\tUniform: " << buildTime[0] << "\n";
cout << "\tStatic: " << buildTime[1] << "\n";
cout << "\tOctree: " << buildTime[2] << "\n";
cout << "\tKD Tree: " << buildTime[3] << "\n";
cout << "Closest point queries\n";
cout << "\tUniform: " << cpTime[0] << "\n";
cout << "\tStatic: " << cpTime[1] << "\n";
cout << "\tOctree: " << cpTime[2] << "\n";
cout << "\tKD Tree: " << cpTime[3] << "\n";
cout << "Closest N points queries\n";
cout << "\tUniform: " << cnpTime[0] << "\n";
cout << "\tStatic: " << cnpTime[1] << "\n";
cout << "\tOctree: " << cnpTime[2] << "\n";
cout << "\tKD Tree: " << cnpTime[3] << "\n";
cout << "Closest points within radius queries\n";
cout << "\tUniform: " << crpTime[0] << "\n";
cout << "\tStatic: " << crpTime[1] << "\n";
cout << "\tOctree: " << crpTime[2] << "\n";
cout << "\tKD Tree: " << crpTime[3] << "\n";
cout << "Total time\n";
cout << "\tUniform: " << (buildTime[0] + cpTime[0] + cnpTime[0] + crpTime[0]) << "\n";
cout << "\tStatic: " << (buildTime[1] + cpTime[1] + cnpTime[1] + crpTime[1]) << "\n";
cout << "\tOctree: " << (buildTime[2] + cpTime[2] + cnpTime[2] + crpTime[2]) << "\n";
cout << "\tKD Tree: " << (buildTime[3] + cpTime[3] + cnpTime[3] + crpTime[3]) << "\n";
timer->Delete();
points->Delete();
qPoints->Delete();
polydata->Delete();
closest->Delete();
// Always return success, although the test infrastructure should catch
// excessive execution times.
return 0;
}
......@@ -28,11 +28,13 @@ set(Module_SRCS
vtkFlyingEdges2D.cxx
vtkFlyingEdges3D.cxx
vtkFlyingEdgesPlaneCutter.cxx
vtkGaussianKernel.cxx
vtkGlyph2D.cxx
vtkGlyph3D.cxx
vtkHedgeHog.cxx
vtkHull.cxx
vtkIdFilter.cxx
vtkInterpolationKernel.cxx
vtkMarchingCubes.cxx
vtkMarchingSquares.cxx
vtkMaskFields.cxx
......@@ -43,6 +45,7 @@ set(Module_SRCS
vtkMergeFields.cxx
vtkMergeFilter.cxx
vtkPointDataToCellData.cxx
vtkPointInterpolator.cxx
vtkPolyDataConnectivityFilter.cxx
vtkPolyDataNormals.cxx
vtkProbeFilter.cxx
......@@ -51,6 +54,8 @@ set(Module_SRCS
vtkRearrangeFields.cxx
vtkResampleToImage.cxx
vtkReverseSense.cxx
vtkSPHKernel.cxx
vtkShepardKernel.cxx
vtkSimpleElevationFilter.cxx
vtkSmoothPolyDataFilter.cxx
vtkStripper.cxx
......@@ -66,6 +71,7 @@ set(Module_SRCS
vtkTubeFilter.cxx
vtkVectorDot.cxx
vtkVectorNorm.cxx
vtkVoronoiKernel.cxx
vtkWindowedSincPolyDataFilter.cxx
vtkCutter.cxx
......@@ -85,6 +91,7 @@ set(Module_SRCS
set_source_files_properties(
vtkEdgeSubdivisionCriterion
vtkInterpolationFunction
vtkStreamerBase
ABSTRACT
)
......
......@@ -22,6 +22,7 @@ vtk_add_test_python(
TestFlyingEdgesPlaneCutter.py
TestGridSynchronizedTemplates3D.py
TestMarchingSquares.py
TestPointInterpolator.py
TestProbeFilterImageInput.py
TestRectilinearSynchronizedTemplates.py
TestScalarTrees.py
......
#!/usr/bin/env python
import vtk
from vtk.test import Testing
from vtk.util.misc import vtkGetDataRoot
VTK_DATA_ROOT = vtkGetDataRoot()
# Parameters for debugging
res = 1000
# create pipeline
#
pl3d = vtk.vtkMultiBlockPLOT3DReader()
pl3d.SetXYZFileName(VTK_DATA_ROOT + "/Data/combxyz.bin")
pl3d.SetQFileName(VTK_DATA_ROOT + "/Data/combq.bin")
pl3d.SetScalarFunctionNumber(100)
pl3d.SetVectorFunctionNumber(202)
pl3d.Update()
output = pl3d.GetOutput().GetBlock(0)
# Create a probe plane
center = output.GetCenter()
plane = vtk.vtkPlaneSource()
plane.SetResolution(res,res)
plane.SetOrigin(0,0,0)
plane.SetPoint1(10,0,0)
plane.SetPoint2(0,10,0)
plane.SetCenter(center)
plane.SetNormal(0,1,0)
# Reuse the locator
locator = vtk.vtkStaticPointLocator()
locator.SetDataSet(output)
locator.BuildLocator()
# Voronoi kernel------------------------------------------------
voronoiKernel = vtk.vtkVoronoiKernel()
interpolator = vtk.vtkPointInterpolator()
interpolator.SetInputConnection(plane.GetOutputPort())
interpolator.SetSourceData(output)
interpolator.SetKernel(voronoiKernel)
interpolator.SetLocator(locator)
# Time execution
timer = vtk.vtkTimerLog()
timer.StartTimer()
interpolator.Update()
timer.StopTimer()
time = timer.GetElapsedTime()
print("Interpolate Points (Voronoi): {0}".format(time))
intMapper = vtk.vtkPolyDataMapper()
intMapper.SetInputConnection(interpolator.GetOutputPort())
intActor = vtk.vtkActor()
intActor.SetMapper(intMapper)
# Create an outline
outline = vtk.vtkStructuredGridOutlineFilter()
outline.SetInputData(output)
outlineMapper = vtk.vtkPolyDataMapper()
outlineMapper.SetInputConnection(outline.GetOutputPort())
outlineActor = vtk.vtkActor()
outlineActor.SetMapper(outlineMapper)
# Gaussian kernel-------------------------------------------------------
gaussianKernel = vtk.vtkGaussianKernel()
gaussianKernel.SetSharpness(4);
interpolator1 = vtk.vtkPointInterpolator()
interpolator1.SetInputConnection(plane.GetOutputPort())
interpolator1.SetSourceData(output)
interpolator1.SetKernel(gaussianKernel)
interpolator1.SetLocator(locator)
interpolator1.SetNullPointsStrategyToNullValue()
# Time execution
timer.StartTimer()
interpolator1.Update()
timer.StopTimer()
time = timer.GetElapsedTime()
print("Interpolate Points (Gaussian): {0}".format(time))
intMapper1 = vtk.vtkPolyDataMapper()
intMapper1.SetInputConnection(interpolator1.GetOutputPort())
intActor1 = vtk.vtkActor()
intActor1.SetMapper(intMapper1)
# Create an outline
outline1 = vtk.vtkStructuredGridOutlineFilter()
outline1.SetInputData(output)
outlineMapper1 = vtk.vtkPolyDataMapper()
outlineMapper1.SetInputConnection(outline1.GetOutputPort())
outlineActor1 = vtk.vtkActor()
outlineActor1.SetMapper(outlineMapper1)
# Shepard kernel-------------------------------------------------------
shepardKernel = vtk.vtkShepardKernel()
shepardKernel.SetPowerParameter(2)
interpolator2 = vtk.vtkPointInterpolator()
interpolator2.SetInputConnection(plane.GetOutputPort())
interpolator2.SetSourceData(output)
interpolator2.SetKernel(shepardKernel)
interpolator2.SetLocator(locator)
interpolator2.SetNullPointsStrategyToMaskPoints()
# Time execution
timer.StartTimer()
interpolator2.Update()
timer.StopTimer()
time = timer.GetElapsedTime()
print("Interpolate Points (Shepard): {0}".format(time))
intMapper2 = vtk.vtkPolyDataMapper()
intMapper2.SetInputConnection(interpolator2.GetOutputPort())
intActor2 = vtk.vtkActor()
intActor2.SetMapper(intMapper2)
# Create an outline
outline2 = vtk.vtkStructuredGridOutlineFilter()
outline2.SetInputData(output)
outlineMapper2 = vtk.vtkPolyDataMapper()
outlineMapper2.SetInputConnection(outline2.GetOutputPort())
outlineActor2 = vtk.vtkActor()
outlineActor2.SetMapper(outlineMapper2)
# SPH kernel-------------------------------------------------------
SPHKernel = vtk.vtkSPHKernel()
interpolator3 = vtk.vtkPointInterpolator()
interpolator3.SetInputConnection(plane.GetOutputPort())
interpolator3.SetSourceData(output)
interpolator3.SetKernel(SPHKernel)
interpolator3.SetLocator(locator)
# Time execution
timer.StartTimer()
interpolator3.Update()
timer.StopTimer()
time = timer.GetElapsedTime()
print("Interpolate Points (SPH): {0}".format(time))
intMapper3 = vtk.vtkPolyDataMapper()
intMapper3.SetInputConnection(interpolator3.GetOutputPort())
intActor3 = vtk.vtkActor()
intActor3.SetMapper(intMapper3)
# Create an outline
outline3 = vtk.vtkStructuredGridOutlineFilter()
outline3.SetInputData(output)
outlineMapper3 = vtk.vtkPolyDataMapper()
outlineMapper3.SetInputConnection(outline3.GetOutputPort())
outlineActor3 = vtk.vtkActor()
outlineActor3.SetMapper(outlineMapper3)
# Create the RenderWindow, Renderer and both Actors
#
ren0 = vtk.vtkRenderer()
ren0.SetViewport(0,0,.5,.5)
ren1 = vtk.vtkRenderer()
ren1.SetViewport(0.5,0,1,.5)
ren2 = vtk.vtkRenderer()
ren2.SetViewport(0,0.5,.5,1)
ren3 = vtk.vtkRenderer()
ren3.SetViewport(0.5,0.5,1,1)
renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(ren0)
renWin.AddRenderer(ren1)
renWin.AddRenderer(ren2)
renWin.AddRenderer(ren3)
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)
# Add the actors to the renderer, set the background and size
#
ren0.AddActor(intActor)
ren0.AddActor(outlineActor)
ren0.SetBackground(0.1, 0.2, 0.4)
ren1.AddActor(intActor1)
ren1.AddActor(outlineActor1)
ren1.SetBackground(0.1, 0.2, 0.4)
ren2.AddActor(intActor2)
ren2.AddActor(outlineActor2)
ren2.SetBackground(0.1, 0.2, 0.4)
ren3.AddActor(intActor3)
ren3.AddActor(outlineActor3)
ren3.SetBackground(0.1, 0.2, 0.4)
renWin.SetSize(500, 500)
cam = ren0.GetActiveCamera()
cam.SetClippingRange(3.95297, 50)
cam.SetFocalPoint(8.88908, 0.595038, 29.3342)
cam.SetPosition(-12.3332, 31.7479, 41.2387)
cam.SetViewUp(0.060772, -0.319905, 0.945498)
ren1.SetActiveCamera(cam)
ren2.SetActiveCamera(cam)
ren3.SetActiveCamera(cam)
iren.Initialize()
# render the image
#
renWin.Render()
iren.Start()
/*=========================================================================
Program: Visualization Toolkit
Module: vtkVoronoiKernel.cxx
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notice for more information.
=========================================================================*/
#include "vtkVoronoiKernel.h"
#include "vtkAbstractPointLocator.h"
//----------------------------------------------------------------------------
vtkVoronoiKernel::vtkVoronoiKernel()
{
}
//----------------------------------------------------------------------------
vtkVoronoiKernel::~vtkVoronoiKernel()
{
}
//----------------------------------------------------------------------------
void vtkVoronoiKernel::
vtkIdType ComputeWeights(double x[3], vtkIdList *pIds,vtkDoubleArray *weights)
{
}
//----------------------------------------------------------------------------
void vtkVoronoiKernel::PrintSelf(ostream& os, vtkIndent indent)
{
this->Superclass::PrintSelf(os,indent);
}
/*=========================================================================
Program: Visualization Toolkit
Module: vtkGaussianKernel.cxx
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notice for more information.
=========================================================================*/
#include "vtkGaussianKernel.h"
#include "vtkAbstractPointLocator.h"
#include "vtkObjectFactory.h"
#include "vtkIdList.h"
#include "vtkDoubleArray.h"
#include "vtkDataSet.h"
#include "vtkPointData.h"
#include "vtkMath.h"
vtkStandardNewMacro(vtkGaussianKernel);
//----------------------------------------------------------------------------
vtkGaussianKernel::vtkGaussianKernel()
{
this->Radius = 1.0;
this->Sharpness = 2.0;
}
//----------------------------------------------------------------------------
vtkGaussianKernel::~vtkGaussianKernel()
{
}
//----------------------------------------------------------------------------
vtkIdType vtkGaussianKernel::
ComputeWeights(double x[3], vtkIdList *pIds, vtkDoubleArray *weights)
{
this->Locator->FindPointsWithinRadius(this->Radius, x, pIds);
vtkIdType numPts = pIds->GetNumberOfIds();
if ( numPts >= 1 ) //Use Gaussian kernel
{
int i;
vtkIdType id;
double d2, y[3], sum = 0.0;
double h2 = (this->Radius*this->Radius) / (this->Sharpness*this->Sharpness);
weights->SetNumberOfTuples(numPts);
double *w = weights->GetPointer(0);
for (i=0; i<numPts; ++i)
{
id = pIds->GetId(i);
this->DataSet->GetPoint(id,y);
d2 = vtkMath::Distance2BetweenPoints(x,y);
if ( d2 == 0.0 ) //precise hit on existing point
{
pIds->SetNumberOfIds(1);
pIds->SetId(0,id);
weights->SetNumberOfTuples(1);
weights->SetValue(0,1.0);
return 1;
}
else
{
w[i] = exp(-(d2/h2));
sum += w[i];
}
}//over all points
// Normalize
for (i=0; i<numPts; ++i)
{
w[i] /= sum;
}
return numPts;
}//using Gaussian Kernel
else //null point
{
return 0;
}
}
//----------------------------------------------------------------------------
void vtkGaussianKernel::PrintSelf(ostream& os, vtkIndent indent)
{
this->Superclass::PrintSelf(os,indent);
}
/*=========================================================================
Program: Visualization Toolkit
Module: vtkGaussianKernel.h
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notice for more information.
=========================================================================*/