# InterpolateTerrain

vtk-examples/Cxx/PolyData/InterpolateTerrain

### Description¶

This example samples a "terrain map" using two approaches.

• Creates an image from the x,y points of the terrain and then uses vtkProbeFilter to interpolate the heights.
• Uses vtkCellLocator to directly interpolate the terrain polydata.

Note that the results differ when the point is not one of the original terrain points. This is because the image has quadrilateral elements and the polydata has triangles. As the resolution of the grid increases, the two results converge.

Question

### Code¶

InterpolateTerrain.cxx

#include <vtkCellArray.h>
#include <vtkCellLocator.h>
#include <vtkDelaunay2D.h>
#include <vtkDoubleArray.h>
#include <vtkImageData.h>
#include <vtkLine.h>
#include <vtkMath.h>
#include <vtkMinimalStandardRandomSequence.h>
#include <vtkNew.h>
#include <vtkPointData.h>
#include <vtkPoints.h>
#include <vtkPolyData.h>
#include <vtkProbeFilter.h>
#include <vtkTriangle.h>
#include <vtkXMLPolyDataWriter.h>

int main(int, char*[])
{
vtkNew<vtkImageData> image;
image->SetExtent(0, 9, 0, 9, 0, 0);
image->AllocateScalars(VTK_DOUBLE, 1);

// Create a random set of heights on a grid. This is often called a
//"terrain map"
vtkNew<vtkPoints> points;

vtkNew<vtkMinimalStandardRandomSequence> randomSequence;
randomSequence->SetSeed(8775070);
unsigned int GridSize = 10;
for (unsigned int x = 0; x < GridSize; x++)
{
for (unsigned int y = 0; y < GridSize; y++)
{
double val = randomSequence->GetRangeValue(-1, 1);
randomSequence->Next();
points->InsertNextPoint(x, y, val);
image->SetScalarComponentFromDouble(x, y, 0, 0, val);
}
}

// add the grid points to a polydata object
vtkNew<vtkPolyData> polydata;
polydata->SetPoints(points);

// triangulate the grid points
vtkNew<vtkDelaunay2D> delaunay;
delaunay->SetInputData(polydata);
delaunay->Update();

vtkNew<vtkXMLPolyDataWriter> writer;
writer->SetFileName("surface.vtp");
writer->SetInputConnection(delaunay->GetOutputPort());
writer->Write();

// Add some points to interpolate
vtkNew<vtkPoints> probePoints;
probePoints->InsertNextPoint(5.2, 3.2, 0);
probePoints->InsertNextPoint(5.0, 3.0, 0);
probePoints->InsertNextPoint(0.0, 0.0, 0);

vtkNew<vtkPolyData> probePolyData;
probePolyData->SetPoints(probePoints);

vtkNew<vtkProbeFilter> probe;
probe->SetSourceData(image);
probe->SetInputData(probePolyData);
probe->Update();

vtkDataArray* data = probe->GetOutput()->GetPointData()->GetScalars();
vtkDoubleArray* doubleData = dynamic_cast<vtkDoubleArray*>(data);

for (int i = 0; i < doubleData->GetNumberOfTuples(); i++)
{
double val = doubleData->GetValue(i);
cout << "Interpolation using ProbeFilter ";
cout << "doubleData->GetValue(" << i << "): " << val << endl;
}

// Now find the elevation with a CellLocator
vtkNew<vtkCellLocator> cellLocator;
cellLocator->SetDataSet(delaunay->GetOutput());
cellLocator->BuildLocator();

for (int i = 0; i < doubleData->GetNumberOfTuples(); i++)
{
int subId;
double t, xyz[3], pcoords[3];
double rayStart[3], rayEnd[3];
probePoints->GetPoint(i, rayStart);
rayStart[2] += 1000.0;
probePoints->GetPoint(i, rayEnd);
rayEnd[2] -= 1000.0;

if (cellLocator->IntersectWithLine(rayStart, rayEnd, 0.0001, t, xyz,
pcoords, subId))
{
cout << "Interpolation using CellLocator ";
cout << "Elevation at " << rayStart[0] << ", " << rayStart[1] << " is "
<< xyz[2] << endl;
}
}
return EXIT_SUCCESS;
}


### CMakeLists.txt¶

cmake_minimum_required(VERSION 3.12 FATAL_ERROR)

project(InterpolateTerrain)

find_package(VTK COMPONENTS
CommonCore
CommonDataModel
FiltersCore
IOXML
)

if (NOT VTK_FOUND)
message(FATAL_ERROR "InterpolateTerrain: Unable to find the VTK build folder.")
endif()

# Prevent a "command line is too long" failure in Windows.
set(CMAKE_NINJA_FORCE_RESPONSE_FILE "ON" CACHE BOOL "Force Ninja to use response files.")
target_link_libraries(InterpolateTerrain PRIVATE ${VTK_LIBRARIES} ) # vtk_module_autoinit is needed vtk_module_autoinit( TARGETS InterpolateTerrain MODULES${VTK_LIBRARIES}
)


cd InterpolateTerrain/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:

./InterpolateTerrain


WINDOWS USERS

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