Sunday, August 8, 2010

Reading a vtk file with multiple scalar fields

Sometimes I come across data files in .vtk format. One of these files stores multiple scalar values in a single file, where the scalar are stored as structured points. For example I received a file containing quantA, quantB and quantC. I would like to read in these files and convert them to something else. So the main goal was to obtain arrays with the 3 scalar fields. The code is not yet complete but the reading is quite involved and it took me a lot of time to patch the required information together from all different examples and the DOxygen online manual of VTK (which I have come to dislike to a great extend). Perhaps the VTK books would have helped me out here but I am not willing to buy them just for this single application of VTK.

Okay here we go:
//First you have to create a reader object:
vtkSmartPointer < vtkstructuredpointsreader> reader = vtkSmartPointer < vtkstructuredpointsreader>::New();
reader->SetFileName(inputFilename.c_str());
//Now we can get the number of scalars in the file
vtkIdType numScal=reader->GetNumberOfScalarsInFile();
std::cout << "number of scalars in file: " << numScal << std::endl;
//Show which scalars we have
int i;
for(i=0;i< numScalComp;i++){
  std::cout << "Scalar "<< i <<": " << reader->GetScalarsNameInFile(i)<< std::endl;
}
//To read a certain scalar field we have to tell the reader which one we want
//we do that by name in this example we take the third one
reader->SetScalarsName(reader->GetScalarsNameInFile(2));
reader->Update();  //I think this actually makes the reader do something
//the rest is what vtk people call setting up the pipeline
//To get to the point data we need a few intermediate steps and objects
vtkStructuredPoints* structuredPoints = reader->GetOutput();
vtkPointData *pd=structuredPoints->GetPointData();
pd->Update();
vtkDataArray *scalars ;
scalars=pd->GetScalars(reader->GetScalarsNameInFile(2));
//here the same name as selected in the reader must be used!!
//I tried this line with the SetScalarsName call to the reader and
//that doesn't work.
//Now we can access the points for example as
vtkIdType numPoints = structuredPoints->GetNumberOfPoints();
for(i=0;i< numPoints;i++){
   std::cout<< i<< ": " << scalars->GetComponent(i,0)<< std::endl;
}
The next step is to write a .vtk file from my own data. The mayavi2 and ParaView programs look quite nice for imaging of 3d data.

1 comment:

  1. Nope, no... The VTK books are just the doxygen output printed and bound. Absolutely appalling boatload of koolaid drinkers regarding the "rightness" of their documentation scheme amongst that crowd.

    ReplyDelete