我想编写一个 python 脚本来生成一个类似于 Paraview 的下一个屏幕截图中右侧所示的绘图:
我有一系列使用命令生成的文件foamToVTK
:
VTK 中是否有类似于PlotOverLine
Paraview 方法的功能?
ParaViewPlotOverLine
在 vtk 中为vtkProbeFilter
.
最小的工作示例:
import vtk
# Original data
source = vtk.vtkRTAnalyticSource()
# the line to plot over
line = vtk.vtkLineSource()
# filter
probeFilter = vtk.vtkProbeFilter()
probeFilter.SetInputConnection(line.GetOutputPort())
probeFilter.SetSourceConnection(source.GetOutputPort())
# rendering
plot = vtk.vtkXYPlotActor()
plot.AddDataSetInputConnection(probeFilter.GetOutputPort())
plot.SetTitle("My plot")
window = vtk.vtkRenderWindow()
interactor = vtk.vtkRenderWindowInteractor()
interactor.SetRenderWindow(window)
renderer = vtk.vtkRenderer()
renderer.AddActor2D(plot)
window.AddRenderer(renderer)
window.Render()
interactor.Start()
您可以在此处找到更复杂的(多图、颜色 ...)c++ 示例:https ://lorensen.github.io/VTKExamples/site/Cxx/Annotation/XYPlot/ python API 是相同的。
使用vtkplotter的解决方案:
from vtkplotter import *
img = load(datadir+'vase.vti').imagedata()
# write a test file, return a vtkMultiBlockData
mblock = write([img], "multiblock.vtm")
# load back from file into a list of meshes/volumes
mobjs = load("multiblock.vtm")
p1, p2 = (10,10,10), (80,80,80)
pl = probeLine(mobjs[0], p1, p2).lw(3)
vals = pl.getPointArray()
xvals = pl.points()[:,0]
plt = plotxy([xvals, vals],
lc="r", # line color
marker="*", # marker style
mc="dr", # marker color
ms=0.6, # marker color
)
show([(img, pl), plt], N=2, sharecam=False, axes=1, bg='w')
我想出了解决这个问题的办法。不过,这不太可能是最佳选择。对于这个解决方案,我首先将网格转换为 a vtkUnstructuredGrid
(在本例中使用 256 点的分辨率。)
下面我包括我用来执行此操作的代码:
import matplotlib.pyplot as plt
from scipy.interpolate import griddata
import numpy as np
import vtk
from vtk.util.numpy_support import vtk_to_numpy
from os import walk, path, system
import pandas as pd
### Create the VTK files
system("foamToVTK")
### Initialization of variables
cnt=1
fig= plt.figure()
npts = 256 #dimensions of the grid
### Get the file names of each step of the simulation
(dirpath, dirnames, filenames) = next(walk('VTK'))
ids=[]
for dir in dirnames:
ids.append(int(dir.split("_")[1]))
ids = sorted(ids)
basename = dirnames[0].split("_")[0]
### Iteration of time steps
for id in ids[1:]:
### Read values from the file of this time step
filename = "%s/%s_%d/internal.vtu" % (dirpath, basename, id)
reader = vtk.vtkXMLUnstructuredGridReader()
reader.SetFileName(filename)
reader.Update()
### Get the coordinates of nodes in the mesh using VTK methods
nodes_vtk_array= reader.GetOutput().GetPoints().GetData()
vtk_array = reader.GetOutput().GetPointData().GetArray('U') #Velocity (3 dimensions)
numpy_array = vtk_to_numpy(vtk_array)
nodes_nummpy_array = vtk_to_numpy(nodes_vtk_array)
x,y,z= nodes_nummpy_array[:,0] , nodes_nummpy_array[:,1] , nodes_nummpy_array[:,2]
xmin, xmax = min(x), max(x)
ymin, ymax = min(y), max(y)
### Define grid
xi = np.linspace(xmin, xmax, npts)
yi = np.linspace(ymin, ymax, npts)
### Grid the data
interpolated = griddata((x, y), numpy_array, (xi[None,:], yi[:,None]), method='cubic')
### Create the list of points
plotOverLine=[]
for point in range(len(interpolated[0])):
plotOverLine.append(interpolated[127][point])
### Update and plot the chart for this time step
df = pd.DataFrame(plotOverLine, columns=['X', 'Y', 'Z'])
plt.clf()
plt.title('Frame %d' % cnt)
plt.plot(df)
plt.legend(df.columns)
axes = plt.gca()
axes.set_ylim([-15,10])
plt.draw()
plt.pause(.05)
对于每个时间步,它都会更新并显示如下图:
受@Nico Vuaille 和https://stackoverflow.com/a/21990296/4286662答案的启发,我遇到了另一种基于探针过滤器的不同解决方案。该解决方案忠实于 Paraview 结果。这是工作代码:
from vtk.util import numpy_support as VN
from matplotlib import pyplot as plt
import vtk
import numpy as np
from os import walk, path, system
import pandas as pd
def getFilenames():
### Get the file names of each step of the simulation
(dirpath, dirnames, filenames) = next(walk('VTK'))
ids=[]
for dir in dirnames:
ids.append(int(dir.split("_")[1]))
ids = sorted(ids)
basename = dirnames[0].split("_")[0]
return ids, basename, dirpath
def createLine(p1, p2):
# Create the line along which you want to sample
line = vtk.vtkLineSource()
line.SetResolution(numPoints)
line.SetPoint1(p1)
line.SetPoint2(p2)
line.Update()
return line
def probeOverLine(line,reader):
### Interpolate the data from the VTK-file on the created line.
probe = vtk.vtkProbeFilter()
probe.SetInputConnection(line.GetOutputPort())
probe.SetSourceConnection(reader.GetOutputPort())
probe.Update()
### Get the velocity from the probe
return VN.vtk_to_numpy(probe.GetOutput().GetPointData().GetArray('U'))
### Initialization of variables
cnt=1
fig= plt.figure()
numPoints=100
ids, basename, dirpath = getFilenames()
### Iteration of time steps
for id in ids[1:]:
### Read values from the file of this time step
filename = "%s/%s_%d/internal.vtu" % (dirpath, basename, id)
reader = vtk.vtkXMLUnstructuredGridReader()
reader.SetFileName(filename)
reader.Update()
if cnt==1:
### Get points for a line in the center of the object
bounds = reader.GetOutput().GetBounds()
p1 = [(bounds[1] - bounds[0]) / 2.0 + bounds[0], bounds[2], 0]
p2 = [(bounds[3] - bounds[2]) / 2.0 + bounds[2], bounds[3], 0]
line = createLine(p1, p2)
plotOverLine = probeOverLine(line, reader)
df = pd.DataFrame(plotOverLine, columns=['X', 'Y', 'Z'])
plt.clf()
plt.title('Frame %d' % cnt)
plt.plot(df)
plt.legend(df.columns)
axes = plt.gca()
plt.draw()
plt.pause(.05)
cnt+=1
plt.show()
结果如下: