用Mayavi/Python从数据生成3D等高线图

Question

我想使用Mayavi绘制一个三维等值线图，与此页面上的第三个图形（氢电子云模型）完全相同：

http://www.sethanil.com/python-for-reseach/5

我有一组数据点，这些点是我用自己的模型创建的，我想要使用它们。这些数据点存储在一个多维numpy数组中，如下所示：

XYZV = [[1, 2, 3, 4],
        [6, 7, 8, 9],
        ...
        [4, 5, 6, 7]]

这些数据点在XYZ空间中不是均匀分布的，也没有以任何特定的顺序存储。我认为示例使用网格来生成数据点 - 我查阅了一些资料，但完全不理解。任何帮助将不胜感激。

_{(来源：sethanil.com)}

- joshlk

---

展示一下你目前尝试过的内容，我们会从那里开始帮助你。 - Don Question

仅供参考，像这样的问题非常适合[scicomp.SE]。 - David Z

Answer 1

关键是在绘图之前对网格进行插值 - 我会使用 scipy 进行插值。下面的R是一个(500,3)的XYZ数值数组，V是每个XYZ点的"大小"。

from scipy.interpolate import griddata
import numpy as np

# Create some test data, 3D gaussian, 200 points
dx, pts = 2, 100j

N = 500
R = np.random.random((N,3))*2*dx - dx
V = np.exp(-( (R**2).sum(axis=1)) )

# Create the grid to interpolate on
X,Y,Z = np.mgrid[-dx:dx:pts, -dx:dx:pts, -dx:dx:pts]

# Interpolate the data
F = griddata(R, V, (X,Y,Z))

从这里开始，展示我们的数据就很容易：

from mayavi.mlab import *
contour3d(F,contours=8,opacity=.2 )

这将产生一个漂亮（不规则）的高斯分布。

请查看griddata的文档，注意您可以更改插值方法。如果您有更多点（在插值网格和数据集上），插值效果会越来越好，更能代表您试图说明的基础函数。这是10K个点和更细的网格的示例：

- Hooked

---

非常感谢。它像魔法一样运行！只有一个问题：如果我想在拟合网格上将点数加倍，我需要在“dx，pts = 2，100j”的那一行更改什么？ - joshlk

@Josh，你可以将其更改为dx, pts = 2, 200j，但这会使_每个维度_中的点数翻倍，因此您将有8倍的点来进行插值。 dx控制绘图网格的范围。为了更精细的控制，请为每个线性维度使用mgrid。 - Hooked

我发现griddata函数计算起来可能需要很长时间。你知道有什么加速这个过程的技巧吗？ - joshlk

@Josh，你可以尝试将插值方式改为最近邻而不是线性的，我没有进行过时间测试，所以不确定这样做会有多大的差异。此外，你可以将图形分成几个部分——高密度区域可以有更多的点/体素，而低密度区域则可以更稀疏。不过，你需要为此使用单独的“griddata”命令。 - Hooked

Answer 2

你可以使用delaunay3d过滤器从点云中创建单元格。然后你可以为delaunay3d的输出UnstructuredGrid创建iso_surface（）表面。如果你需要ImageData，则可以使用image_data_probe 过滤器。

import numpy as np
from tvtk.api import tvtk
from mayavi import mlab

points = np.random.normal(0, 1, (1000, 3))
ug = tvtk.UnstructuredGrid(points=points)
ug.point_data.scalars = np.sqrt(np.sum(points**2, axis=1))
ug.point_data.scalars.name = "value"
ds = mlab.pipeline.add_dataset(ug)
delaunay = mlab.pipeline.delaunay3d(ds)
iso = mlab.pipeline.iso_surface(delaunay)
iso.actor.property.opacity = 0.1
iso.contour.number_of_contours = 10
mlab.show()

- HYRY

---

嗨，当我使用自己的数据使用这个方法时，会弹出许多错误窗口，最后一个警告是：“警告：在C：\ pisi \ tmp \ VTK-5.6.0-2 \ work \ VTK \ Graphics \ vtkDelaunay3D.cxx，第488行vtkDelaunay3D（03AD2250）：遇到27个退化三角形，网格质量可疑”。你知道发生了什么吗？ - joshlk