CODEWITHSUNDEEP
pythonnumpyvtk
Phtagen
Phtagen

Reputation: 121

How to export 3D vector field from numpy array to *.vtk-file using pyvtk?

I'm struggling to export some 3D vector-arrays (numpy arrays) from python to a *.vtk-file for later use in ParaView.

I have three 3D MR-Velocimetry Images, each of 100x100x200 voxels containing the velocity components in x, y and z. What I want is to export this vector field to a *.vtk-file using the pyvtk-module from here.

Unfortunately I don't understand how it works :-(

What I tried so far:

from pyvtk import *

vectors = [flow['vx'], flow['vy'], flow['vz']]

dim=flow['vx'].shape

pointdata=PointData(Vectors([flow['vx'],flow['vy'],flow['vz']]))

vtk=VtkData(StructuredPoints(dim[0],dim[1],dim[2]), pointdata)

Where flow['...'] contains the vector-components. I got the following error:

ValueError: DataSet (size=100) and PointData (size=3) have different sizes

Yeeees, what is it trying to tell me? Ok, guess something like dimension mismatch, but how do I set up my input properly?

Any help would be appreciated. Thanks in advance

Upvotes: 3

Views: 2028

Answers (1)

Phtagen
Phtagen

Reputation: 121

I found a proper solution for my problem, using TVTK instead of PyVTK. So everyone who is interested, a possible workaround could be as follows:

from tvtk.api import tvtk, write_data
# Unpack velocity information

vx=flow['vx']
vy=flow['vy']
vz=flow['vz']

dim=vx.shape

# Generate the grid
xx,yy,zz=np.mgrid[0:dim[0],0:dim[1],0:dim[2]]
pts = empty(vx.shape + (3,), dtype=int)
pts[..., 0] = xx
pts[..., 1] = yy
pts[..., 2] = zz

vectors = empty(vx.shape + (3,), dtype=float)
vectors[..., 0] = vx
vectors[..., 1] = vy
vectors[..., 2] = vz

# We reorder the points and vectors so this is as per VTK's
# requirement of x first, y next and z last.
pts = pts.transpose(2, 1, 0, 3).copy()
pts.shape = pts.size // 3, 3

vectors = vectors.transpose(2, 1, 0, 3).copy()
vectors.shape = vectors.size // 3, 3

sg = tvtk.StructuredGrid(dimensions=xx.shape, points=pts)

sg.point_data.vectors = vectors
sg.point_data.vectors.name = 'velocity'

write_data(sg, 'vtktest.vtk')

Greetings

Phtagen

Upvotes: 3

Related Questions