当我分别使用 MATLAB 或 Python 将数据导出到 *.vtk 文件时,为什么会得到不同的结果?

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

我在将矢量数据导出为 *.vtk 文件格式以便稍后在 ParaView 中使用时遇到了一些问题。然而,从现在开始我一直在使用 MATLAB,尤其是一个名为 vtkwrite.m 的脚本,可以在 here 中找到它。到目前为止这工作正常,但由于许可原因,我想使用 tvtk 更改为 Python。

我成功地将tvtk和Python的矢量数据导出为*.vtk格式,但与MATLAB导出的数据相比,文件有很大不同!首先,MATLAB 版本几乎是 Python 版本的两倍(67.2MB 到 46.2MB)。此外,当我在 ParaView 中可视化流线时,两个数据看起来完全不同。 MATLAB 数据比 Python 版本更平滑。

造成这些差异的原因是什么?

这里是我用来导出数据的一些编码。考虑 vx, vy, vz 我要处理的 3D 矢量速度分量。

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')

如您所见,Python 中的工作流程要困难得多,所以我可能在这里犯了错误?!

如果 Python 代码对你有用,我会将其包装在一个函数中并按如下方式简单化:

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)

可以这样使用:

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

此代码是 modulo 编写即时未经测试的代码时发生的愚蠢错误。