vtk value at x,y,z point

Question

I have a vtk file which maps temperature in 3 dimensions. I would like to determine the temperature at a given x, y, z point. I will use the following code in order to load the vtk file (Reading .vtk file):

int main(int argc, char *argv[]) 
{
  // simply set filename here (oh static joy)
  std::string inputFilename = "setYourPathToVtkFileHere";

  // Get all data from the file
  vtkSmartPointer<vtkGenericDataObjectReader> reader =
        vtkSmartPointer<vtkGenericDataObjectReader>::New();
  reader->SetFileName(inputFilename.c_str());
  reader->Update();

  // All of the standard data types can be checked and obtained like this:
  if (reader->IsFilePolyData()) 
  {
    std::cout << "output is a polydata" << std::endl;
    vtkPolyData* output = reader->GetPolyDataOutput();
    std::cout << "output has " << output->GetNumberOfPoints() << "     points." << std::endl;
  }

  return EXIT_SUCCESS;
}

However, when searching through the extensive list of methods in the vtk library I cant find the appropriate function to extract a value at a specific location. Any suggestions?

Answer 1

The proper way to retrieve a scalar value at a given position depends on two questions:

1. How is your data laid out and
2. From which position do you want to retrieve an attribute

Concerning the data layout there are two main layouts:

• Structured: The data resides in a uniform grid
• Unstructured: The point samples are arbitrary

Concerning the position you can have two situations:

• Query at a sample position: You ask for a point that is directly a sample in your dataset
• Query at an arbitrary position: You ask for a point that is somewhere in your domain but does not neccessarily coincide with a sample of your data.

Independent of the data layout, to retrieve the data at a sample position (i.e. a sample of your original data set) you can use the vtkPointLocator class. Use the class as follows (untested):

// Build locator object
vtkSmartPointer<vtkPointLocator> locator = vtkPointLocator::New();
locator->SetDataSet(polyData);
locator->BuildLocator();
// Define query position
double pt[3] = {0.1, 0.2, 0.3};
// Get the ID of the point that is closest to the query position
vtkIdType id = locator->FindClosestPoint(pt);
// Retrieve the first attribute value from this point
double value = polyData->GetPointData()->GetScalars()->GetTuple(id, 0);

This will give you the point value of the closest sample of your data. Note that this will not give you the explicit position of the point in the data set as this is implicitly encoded in the variable id. To retrieve the actual position of the closest point, you can write:

double *dataPt = polyData->GetPoint(id);

If you want to retrieve data at an arbitrary position of your domain you will need some way of interpolation. Here, the data layout is important.

• For structured data you can first convert your data to a vtkImage and then perform queries on it. If you want to retrieve an interpolated attribute using linear or cubic schemes you can add a vtkImageInterpolator in your filter chain and then retrieve a point using the GetScalarComponentAsDouble method.
• For unstructured data you should first decide on an interpolation scheme. vtk has various filters to reconstruct continuous data from data samples. Options include Delaunay triangulation/tetrahedrization (vtkDelaunay2D, vtkDelaunay3D) , as well as the Shepard method (vtkShepardMethod). Either method will give you a new data set with possibilities to query arbitrary points. If you want to retrieve the scalar attributes for a (batch of) points without actually reconstructing a complete data set you can also look at the vtkProbeFilter.

Answer 2

You need to first extract the polyData from the reader. Then, store the points via vtkPolyData::getPoints into a vtksmartPointer<vtkPoints>. To finish, create a std::vector of a custom struct and store them while iterating over your vtkPoints.

Here's some code to illustrate :

#include <vtkDataArray.h>
#include <vtkDataSet.h>
#include <vtkGenericDataObjectReader.h>
#include <vtkPointLocator.h>
#include <vtkPointData.h>
#include <vtkPolyData.h>
#include <vtkSmartPointer.h>
#include <vtkStructuredGrid.h>
#include <string>


struct Pnt {
    double x_, y_, z_;
    Pnt(double x, double y, double z) : x_(x), y_(y), z_(z) {}
};

    int main ( int argc, char *argv[] )
{
  // Ensure a filename was specified
  if(argc != 2)
    {
    std::cerr << "Usage: " << argv[0] << " InputFilename" << endl;
    return EXIT_FAILURE;
    }

  // simply set filename here (oh static joy)
  std::string inputFilename = "setYourPathToVtkFileHere";

  // Get all data from the file
  vtkSmartPointer<vtkGenericDataObjectReader> reader = vtkSmartPointer<vtkGenericDataObjectReader>::New();
  reader->SetFileName(inputFilename.c_str());
  reader->Update();

  vtkSmartPointer<vtkPolyData> polydata = reader->GetPolyDataOutput();
  vtkSmartPointer<vtkPoints> vtk_points = polydata->GetPoints();

  std::vector<Pnt> my_points;
  for (int i = 0; i < vtk_points->GetNumberOfPoints(); i++){
      const auto pnt = vtk_points->GetPoint(i);
      my_points.emplace_back(pnt[0], pnt[1], pnt[2]);
  }

  return EXIT_SUCCESS;
}

And here is the version with the vtkPointLocator mentionned in QnD's answer :

int main(int argc, char *argv[])
{
    // Ensure a filename was specified
    if (argc != 2)
    {
        std::cerr << "Usage: " << argv[0] << " InputFilename" << endl;
        return EXIT_FAILURE;
    }

    // simply set filename here (oh static joy)
    std::string inputFilename = "setYourPathToVtkFileHere";

    // Get all data from the file
    vtkSmartPointer<vtkGenericDataObjectReader> reader = vtkSmartPointer<vtkGenericDataObjectReader>::New();
    reader->SetFileName(inputFilename.c_str());
    reader->Update();

    vtkSmartPointer<vtkPolyData> polydata = reader->GetPolyDataOutput();

    //Building locator
    vtkSmartPointer<vtkPointLocator> locator = vtkPointLocator::New();
    locator->SetDataSet(polydata);
    locator->BuildLocator();

    //Finding point
    const double pt[3] = { 0.1, 0.2, 0.3 };
    vtkIdType id = locator->FindClosestPoint(pt);
    double pnt_found[3];
    polydata->GetPointData()->GetScalars()->GetTuple(id, pnt_found);

    return EXIT_SUCCESS;
}

And the CMakeLists.txt

cmake_minimum_required(VERSION 2.8)

PROJECT(GenericDataObjectReader)

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

add_executable(GenericDataObjectReader MACOSX_BUNDLE GenericDataObjectReader)

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

Please let me know if I need to do any specific action to correctly credit the other answer. I'm new and don't know all the specifics yet.