Showing posts with label vtk. Show all posts
Showing posts with label vtk. Show all posts

2012-06-30

VTK Polygons and other cells as vtkCellArray in Python

After hours of googling and playing with iPython, I finally figured out the way to access polygons or other cells in VTK files using Python.

The Python binding of vtk library seems missing a very important function for users to access cells/polygons in Python, as mentioned in a VTK mailing list post back to 2002 and another post in 2011.

A workaround I just found out is as follows. Suppose you have a VTK file called test.vtk containing the following data.

# vtk DataFile Version 2.0
Cube example
ASCII

DATASET POLYDATA
POINTS 8 float
0.0 0.0 0.0
1.0 0.0 0.0
1.0 1.0 0.0
0.0 1.0 0.0
0.0 0.0 1.0
1.0 0.0 1.0
1.0 1.0 1.0
0.0 1.0 1.0

POLYGONS 3 12
3 0 1 2 
3 4 5 6
3 7 4 2

Now I use interactions on iPython to demonstrate the accessing to POLYGONS.

First, we prepare the accessing.
In [1]: import vtk

In [2]: Reader = vtk.vtkDataSetReader()

In [3]: Reader.SetFileName('test.vtk')

In [4]: Reader.Update()

In [5]: Data = Reader.GetOutput()

In [6]: CellArray = Data.GetPolys()

In [7]: Polygons = CellArray.GetData()

Now check the number of cells/polygons and number of points in cells/polygons
In [8]: CellArray.GetNumberOfCells()
Out[8]: 3L

In [9]: Polygons.GetNumberOfTuples()
Out[9]: 12L

All cells/polygons can be accessed like this:
In [10]: for i in xrange(0,  Polygons.GetNumberOfTuples()):
   ....:         print Polygons.GetValue(i)
   ....: 
3
0
1
2
3
4
5
6
3
7
4
2
Please note that the numbers (3's here) indicating sizes of cells (i.e., the numbers at the beginning of every line in Cell/Polygon segment in a VTK file) are also retrieved and printed.

If all your cells/polygons are of the same size, e.g., all triangles, here is an easy way.
In [11]: for i in xrange(0,  CellArray.GetNumberOfCells()):
   ....:         print [Polygons.GetValue(j) for j in xrange(i*3+1, i*3+4) ]
   ....: 
[0L, 1L, 2L]
[3L, 4L, 5L]
[6L, 3L, 7L]

2012-01-07

Double and Float types of data in VTK files

by Forrest Sheng Bao http://fsbao.net

I just realized that they are different, at least to python-vtk, the official Python wrapper for VTK library. If a scalar is of type Float, python-vtk will pad many digits after the decimal point. For example, 205187 in a pure ASCII VTK file becomes 0.20518699288368225. All Float data are padded to the 16-th digit after decimal point. However, if the type is Double, no such a problem.

I don't know what VTK's semantics for Double and Float are. But it should be related.

I am using python-vtk and Python2.7 that comes with Ubuntu Linux 11.10 64bit.

2012-01-04

Adding two SCALARS in POINTDATA for one vtkPolyData object in VTK

by Forrest Sheng Bao http://fsbao.net

In VTK (either the file format or the library), we sometimes associate more than one scalars to points. I just figured out how to do this in VTK (in C++, similarly in its Python, Tcl or Java wrapper).

Suppose I have a vtkPolyData pointer
vtkPolyData* mesh;
and two vtkDoubleArray (you can consider a vtkDoubleArray as a list of scalars) pointers
vtkDoubleArray* depth; 
 vtkDoubleArray* curv; 

This is how I do it:

depth->SetName("Depth");
mesh->GetPointData()->SetScalars(depth);
curv->SetName("Curvature");
mesh->GetPointData()->AddArray(curv);

You may test this by writing mesh into a VTK-format file:

vtkPolyDataWriter* writer=vtkPolyDataWriter::New();
 writer->SetFileName("test_dump.vtk");
 
 writer->SetInput(mesh);
 writer->Update();
 writer->Write();
 writer->Delete();

My only question is whether the use of AddArray() function is correct. What is I wanna set a Normal? or a Tensor?

2011-12-28

An incompatible problem between python-vtk and pyvtk

I am using python-vtk (official Python binding of VTK) and pyvtk (only for VTK-format file I/O) these days. It's not a pleasant experience. It took me a few hours to debug my code and found out the problem at somewhere I never expected - the incompatibility of the libraries I am using. On top of that, one of them has compatibility issue with Python2.7.

I used python-vtk to read from VTK files (because I couldn't google out examples using pyvtk to do so). And then I used pyvtk to write into VTK files (because I couldn't google out examples using python-vtk to create scalar POINTDATA).

The problem is that in python-vtk, vtk.vtkDataSetReader().GetOutput().GetPoint() returns a 3-tuple, which is the X-, Y- and Z-coordinates of a POINT in DATASET POLYDATA block, whereas in pyvtk, pyvtk.PolyData() takes in coordinates as a list, not a 3-tuple.

I don't understand how this could cause problems, because lists and tuples are very similar in Python. After mandatory type conversion, the problem was gone.

pyvtk also has problem with Python2.7. I haven't debugged out the cause. But I just know it does not work with Python2.7.

PS: There is a great lack of resources on these two libraries. I use one library for reading VTK files and the other for writing VTK files, though each of them has the ability to do both. This "complimentary" combination is because I couldn't google out examples on using the opposite library for the opposite function. I really don't understand this. Without document, a software package is nothing.

I didn't find any documentation on pyvtk's website. For python-vtk, released with VTK (C++ mainly, along with Tcl, Python and Java bindings), I still didn't find docs for Python. I had to use iPython to learn.

Does anyone know why?

2011-12-04

Reading VTK files in Python via python-vtk

by Forrest Sheng Bao http://fsbao.net

Update 2012-06-29: I just figured out a way to access polygons/cell using VTK's python binding. Check here: http://forrestbao.blogspot.com/2012/06/vtk-polygons-and-other-cells-as.html

Surprisingly, I noticed that there isn't a good document covering frequently-used functions of the Python module vtk (provided via python-vtk on Debian/Ubuntu Linux systems). So I decide to write a very simple one here, covering all functions that I have used.

Note:
  • This tutorial is for Python 2.X though slight changes can make it work with Python 3.X.
  • I assume you are familiar with VTK data format, thus you know what header, DATASET and POINT_DATA are.
  • Notations like In [123] or Out[123] are prompts in iPython (an interactive Python shell). They show a line of code and its output/effect, respectively. They should NOT appear in your code. And, when you use other Python interpreter/shell, you may not see it.
  • I did not use print function in iPython when showing something. But when you write a Python program, you need print to display.
Step 1: Set up the reader.
import vtk
reader = vtk.vtkDataSetReader()
reader.SetFileName("lh.sulc.fundi.from.pits.pial.vtk")
reader.ReadAllScalarsOn()  # Activate the reading of all scalars
reader.Update()

The reader is the top level object we use to access a VTK file. For example, you can know the header of a VTK file by
In [119]: reader.GetHeader()
Out[119]: 'Created by Mindboggle'

Step 2: Get data in DATASET block.

Again, reader is the top level. To get data in DATASET block, we need to call a method on reader.
data=reader.GetOutput()

We can know how many (not ``much'' here) data of each type are in the DATASET by the function GetNumberOf{VTK_Data_Type}(). For example, to know how may Points are there, we use:
In [11]: data.GetNumberOfPoints()
Out[11]: 128895

I used tab-completion function of iPython to find out all GetNumberOf{VTK_Data_Type}() for data:
data.GetNumberOfCells   data.GetNumberOfPolys
data.GetNumberOfLines   data.GetNumberOfStrips
data.GetNumberOfPieces  data.GetNumberOfVerts
data.GetNumberOfPoints  

After knowing the size of data of each type, we can access them. For example, to get the first point's coordinate, we can do this:
In [136]: data.GetPoint(0)
Out[136]: (-11.605026245117188, -97.47259521484375, 3.8222298622131348)

To access vertexes in the VERTICES block, it is a little different.
In [15]: data.GetNumberOfVerts()
Out[15]: 1L
In [18]: vt = data.GetVerts()
In [20]: vt.GetSize()
Out[20]: 56094L
In [23]: vt.GetData().GetValue(2)
Out[23]: 1091L
In [24]: vt.GetData().GetValue(3)
Out[24]: 1108L
Please note that some people, like myself, prefer to put all vertexes in one line, so the result of data.GetNumberOfVerts() and data.GetVerts().GetSize() look different. But it shouldn't bother you from using data.GetVerts().GetData().GetValue(x) to access them.

Step 3: get data in POINT_DATA/CELL_DATA block.

Data in POINT_DATA and CELL_DATA can be loaded via function GetPointData() and GetCellData(). For example,
d=data.GetPointData()

We will use loading a scalar array from POINT_DATA as an example below.
Vectors and tensors can be accessed in similar way by different functions.

All scalar arrays in a file can be viewed by:
In [140]: reader.GetNumberOfScalarsInFile() # get number of scalars
Out[140]: 5
In [141]: reader.GetScalarsNameInFile(1) # get scalar name string
Out[140]: curv

A scalar can be accessed by its name:
In [38]: array=d.GetArray('curv') # 'curv' is the scalar name
In [50]: array.GetValue(260957-260761+1)
Out[50]: 0.043181419372558594

In the end, you should see the following variables:
Variable   Type         Data/Info
---------------------------------
curv       vtkobject    vtkFloatArray (0x9e5cd38)<...>)\n  Array: 0x9f1a9a0\n\n
d          vtkobject    vtkPointData (0x9789b50)\<...>  PedigreeIds: (none)\n\n
data       vtkobject    vtkPolyData (0x9e6c650)\n<...>: 0\n  Ghost Level: 0\n\n
reader     vtkobject    vtkDataSetReader (0x9c431<...> InputStringLength: 0\n\n
vtk        module       < module 'vtk' from '/usr/<...>hon2.6/vtk/__init__.pyc'>

The end. Comments and questions are welcomed.

References: