Why do I get different results, when I export data to *.vtk-files with MATLAB or Python, respectively?

Question

I got some trouble exporting vector-data to *.vtk-file format for later use in ParaView. However, since now I was using MATLAB and especially a script called vtkwrite.m which can be found here. This works fine so far, but I wanted to change to Python using tvtk because of license reasons.

I managed to export my vector-data with tvtk and Python to *.vtk-format, but compared to the MATLAB exported data, the files are quite different! First of all, the MATLAB version is almost twice as big as the Python version (67.2MB to 46.2MB). Also when I visualize streamlines in ParaView, both data look quite different. The MATLAB data is way more smoother than the Python version.

What is the reason for these differences?

Here some coding I used to export the data. Consider vx, vy, vz the 3D vectorial velocity components I want to process.

1) MATLAB:

[x,y,z]=ndgrid(1:size(vx,1),1:size(vx,2),1:size(vx,3));
vtkwrite('/pathToFile/filename.vtk','structured_grid',x,y,z,'vectors','velocity',vx,vy,vz);

2) Python

from tvtk.api import tvtk, write_data 

dim=vx.shape

xx,yy,zz=np.mgrid[0:dim[0],0:dim[1],0:dim[2]]
pts = empty(dim + (3,), dtype=int)
pts[..., 0] = xx
pts[..., 1] = yy
pts[..., 2] = zz

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

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, '/pathToFile/filename.vtk')

As you can see, the workflow in Python is way more difficult, so maybe I made a mistake here?!

Answer 1

If the Python code is working for you, I would wrap it in a function and simplyfy it as follows:

import numpy as np


def vtk_save(
        filepath,
        v_arr,
        x_arr=None,
        label=None,
        ndim=3):
    base_shape = v_arr.shape[:ndim]
    if not isinstance(v_arr, np.ndarray):
        v_arr = np.stack(v_arr[::-1], -1).reshape(-1, ndim)
    if x_arr is None:
        x_arr = np.stack(
            np.mgrid[tuple(slice(0, dim) for dim in v_arr.shape[::-1])], -1) \
            .reshape(-1, ndim)
    elif not isinstance(x_arr, np.ndarray):
        x_arr = np.stack(x_arr[::-1], -1).reshape(-1, ndim)
    sg = tvtk.StructuredGrid(
        dimensions=base_shape, points=x_arr)
    sg.point_data.vectors = v_arr
    sg.point_data.vectors.name = label
    write_data(sg, filepath)

which could be used like:

vtk_save('/pathToFile/filename.vtk', [vx, vy, vz], label='velocity')

this code is modulo silly bugs that happen when writing on-the-fly untested code.