python 科学计算三维可视化笔记 第三周 高级进阶

python 科学计算三维可视化笔记 第三周 高级进阶一、Mayavi 入门1. Mayavi 库的基本元素2. 快速绘制实例3. Mayavi 管线二、Mlab 基础（一）基于 numpy 数组的绘图函数1. 0D 数据：``points3d()`` 函数2. 1D 数据：``plot3d()`` 函数3. 2D 数据：``imshow()`` 函数4. 2D 数据：``surf()`` 函数5. 2D 数据：``contour_surf()`` 函数6. 2D 数据：其他函数7. 3D 数据：``contour3d()`` 函数8. 3D 数据：``quiver3d()`` 函数9. 3D 数据：``flow()`` 函数（二）改变物体的颜色（三）mlab 控制函数（四）鼠标选取交互操作1. ``on_mouse_pick()`` 函数2. 选取红色小球问题分析3. 代码4. 程序框架的优化5. 优化代码6. 结果（五）mlab 管线控制函数（六）实例1：标量数据可视化1. 绘制等值面2. 绘制切平面3. 复合观测方法（七）实例2：矢量数据可视化1. 使用 ``quiver3d()`` 函数2. masking vector 采样3. Cut Plane 切面4. 级数的等值面5. 流线 flow 的可视化6. 复合观测方法

内容来自中国大学MOOC，北京理工大学，python数据分析与展示课程，侵删。

如有错误，烦请指出。

python 科学计算三维可视化笔记 第三周 高级进阶

一、Mayavi 入门

1. Mayavi 库的基本元素

Mayavi.mlab 模块：

Mayavi API：

2. 快速绘制实例

(1) 实例1

'''快速绘制实例1'''from mayavi import mlab# 建立数据x = [[-1,1,1,-1,-1], [-1,1,1,-1,-1]]y = [[-1,-1,-1,-1,-1], [1,1,1,1,1]]z = [[1,1,-1,-1,1], [1,1,-1,-1,1]]# 对数据进行三维可视化s = mlab.mesh(x,y,z)mlab.show()

(2) 实例2

'''快速绘制实例2'''from numpy import pi, sin, cos, mgridfrom mayavi import mlab# 建立数据dphi, dtheta = pi/250.0, pi/250.0[phi,theta] = mgrid[0:pi+dphi*1.5:dphi,0:2*pi+dtheta*1.5:dtheta]m0 = 4; m1 = 3; m2 = 2; m3 = 3; m4 = 6; m5 = 2; m6 = 6; m7 = 4;r = sin(m0*phi)**m1 + cos(m2*phi)**m3 + sin(m4*theta)**m5 + cos(m6*theta)**m7x = r*sin(phi)*cos(theta)y = r*cos(phi)z = r*sin(phi)*sin(theta)# 对数据进行三维可视化s = mlab.mesh(x, y, z)# s = mlab.mesh(x, y, z, representation='wireframe', line_width=1.0)mlab.show()

mlab.mesh()函数默认为 surface（左图），可以更改为 wireframe（右图）：

可以用键盘和鼠标对场景进行操作：

旋转场景：左键拖动或键盘的方向键平移场景：按住 shift 键并使用左键拖动，shift+方向键盘缩放场景：鼠标右键上下拖动或使用 + 和 - 按键滚动相机：按住 ctrl 键并用左键拖动工具栏：从坐标轴 6 个方向观察场景、等角投影、切换平行透视和成角透视等

3. Mayavi 管线

engine：建立和销毁 scenesscenes：多个数据集合 sourcesfilters：对数据进行变换module manager：控制颜色，Colors and legendsmodules：最终数据的表示，如线条、平面等

点击 Mayavi pipeline 或使用mlab.show_pipeline()打开 mayavi 管线 对话框：

(1) 管线中的对象 scene

Mayavi Scene：处于树的最顶层的对象，表示场景

使用s = mlab.gcf()获取当前场景

使用print(s.scene.background)输出当前场景的背景色

(2) 管线中的对象 GridSource

GridSource 是网格数据源，每一项是标量数据的名称

使用source = s.children[0]获取 GridSource 对象

使用print(source.name)输出名称

使用print(repr(source.data.points))输出坐标

使用print(repr(source.data.point_data.scalars))输出每个点对应的标量数组

(3) 管线中的对象 PolyDataNormals

PolyDataNormals：数据源的法向量

(4) 管线中的对象 Colors and legends

更改 lut mode 颜色为蓝色，勾选 show legend 添加颜色条显示颜色和标量之间的关系：

colors = manager.children[0]colors.scalar_lut_manager.lut_mode = 'Blues'colors.scalar_lut_manager.show_legend = 'True'

结果：

(5) 管线中的对象 Surface

将 GridSource 中输出的 PolyData 数据转换为最终在场景中显示的三维实体，将 representation 改为 wireframe，将 line width 设置为 0，将 opacity (透明度) 设置为 0.6：

surface = colors.children[0]surface.actor.property.representation = 'wireframe'surface.actor.property.opacity = 0.6

结果：

程序配置属性的步骤：

获取场景对象，mlab.gcf()通过 chirldren 属性，在管线中找到需要修改的对象配置窗口有多个选项卡，属性需要一级一级获得通过工具栏 “红色原点” 菜单，可以记录配置过程

二、Mlab 基础

（一）基于 numpy 数组的绘图函数

mlab 对 numpy 建立可视化过程：

建立数据源使用 Filter（可选）添加可视化模块

3D 绘图函数，包括points3d()和plot3d()：

1. 0D 数据：points3d()函数

(1) 函数形式

points3d(x, y, z ...)points3d(x, y, z, s, ...)points3d(x, y, z, f, ...)

x, y, z 表示 numpy 数组、列表或者其他形式的点三维坐标s 表示在该坐标点处的标量值f 表示通过函数f(x, y, z)返回的标量值

(2) 函数参数

(3) 实例

'''points3d()'''import numpy as npfrom mayavi import mlab# 建立数据t = np.linspace(0, 4*np.pi, 20)x = np.sin(2*t)y = np.cos(t)z = np.cos(2*t)s = 2 + np.sin(t)# 对数据进行可视化points = mlab.points3d(x, y, z, s, colormap='Reds', scale_factor=.25)# points = mlab.points3d(x, y, z, s, colormap='Greens', scale_factor=.25)mlab.show()

2. 1D 数据：plot3d()函数

(1) 函数形式

plod3d(x, y, z ...)plot3d(x, y, z, s, ...)

x, y, z 表示 numpy 数组、列表或者其他形式的点三维坐标s 表示在该坐标点处的标量值

(2) 函数参数

(3) 实例

'''plot3d'''import numpy as npfrom mayavi import mlab# 建立数据n_mer, n_long = 6, 11dphi = np.pi / 1000.0phi = np.arange(0.0, 2*np.pi+0.5*dphi, dphi)mu = phi * n_merx = np.cos(mu) * (1+np.cos(n_long*mu/n_mer)*0.5)y = np.sin(mu) * (1+np.cos(n_long*mu/n_mer)*0.5)z = np.sin(n_long*mu/n_mer) * 0.5# 对数据进行可视化l = mlab.plot3d(x, y, z, np.sin(mu), tube_radius=0.025, colormap='Spectral')mlab.show()

3. 2D 数据：imshow()函数

将二维数组可视化为一张图像

(1) 函数形式

imshow(s, ...)

s 是一个二维数组，s 的值使用 colormap 被映射为颜色

(2) 函数参数

(3) 实例

'''imshow()'''import numpy as npfrom mayavi import mlab# 建立数据s = np.random.random((10, 10))# 对数据进行可视化img = mlab.imshow(s, colormap='gist_earth')mlab.show()

4. 2D 数据：surf()函数

将二维数组可视化为一个平面，z 轴描述了数组点的高度

(1) 函数形式

surf(s, ...)surf(x, y, s, ...)surf(x, y, f, ...)

x, y 是一维或二维数组，由np.ogrid或np.mgrid得到s 是一个高程矩阵，用二维数组表示

(2) 实例

'''surf()'''import numpy as npfrom mayavi import mlab# 建立数据def f(x,y):return np.sin(x-y) + np.cos(x+y)x,y = np.mgrid[-7.:7.05:0.1, -5:5.05:0.05]# 对数据进行可视化s = mlab.surf(x, y, f)mlab.show()

5. 2D 数据：contour_surf()函数

将二维数组可视化为等高线，高度值由数组点的值来确定

(1) 函数形式

contour_surf(s, ...)contour_surf(x, y, s, ...)contour_surf(x, y, f, ...)

x, y 是一维或二维数组，由np.ogrid或np.mgrid得到s 是一个高程矩阵，用二维数组表示

(2) 实例

'''contour_surf()'''import numpy as npfrom mayavi import mlab# 建立数据def f(x,y):return np.sin(x-y) + np.cos(x+y)x,y = np.mgrid[-7.:7.05:0.1, -5:5.05:0.05]# 对数据进行可视化s = mlab.contour_surf(x, y, f)mlab.show()

6. 2D 数据：其他函数

mesh()函数绘制由三个二维数组 x、y、z 描述坐标点的网格平面barchart()函数根据二维、三维或者点云数据绘制的三位柱状图triangular_mesh()函数绘制由 x、y、z 坐标点描述的三角网格面

7. 3D 数据：contour3d()函数

三维数组定义的体数据的等值面可视化

(1) 函数形式

contour3d(scalars, ...)contour3d(x, y, z, scalars, ...)

scalars 网格上的数据，用三维 numpy 数组表示x、y、z 是三维空间坐标

(2) 函数参数

contours：定义等值面的数量transparent：True 表示透明

(3) 实例

'''contour3d()'''import numpy as npfrom mayavi import mlab# 建立数据x, y, z = np.ogrid[-5:5:64j, -5:5:64j, -5:5:64j]scalars = x*x + y*y + z*z# 对数据进行可视化obj = mlab.contour3d(scalars, contours=8, transparent=True)mlab.show()

8. 3D 数据：quiver3d()函数

三维矢量数组的可视化，箭头表示在该点的矢量数据

(1) 函数形式

quiver3d(u, v, w, ...)quiver3d(x, y, z, u, v, w, ...)quiver3d(x, y, z, f, ...)

u、v、w 是用 numpy 数组表示的向量x、y、z 表示箭头的位置f 需要返回在给定位置 (x, y, z) 的 (u, v, w) 矢量

(2) 实例

'''quiver3d()'''import numpy as npfrom mayavi import mlab# 建立数据x,y,z = np.mgrid[-2:3, -2:3, -2:3]r = np.sqrt(x**2 + y**2 + z**4)u = y*np.sin(r) / (r+0.001)v = -x*np.sin(r) / (r+0.001)w = np.zeros_like(z)# 对数据进行可视化obj = mlab.quiver3d(x, y, z, u, v, w, line_width=3, scale_factor=1)mlab.show()

9. 3D 数据：flow()函数

绘制三维数组描述的向量场的粒子轨迹

（二）改变物体的颜色

colormap 定义的颜色，也叫 LUT (Look Up Table)LUT 是一个 255×4 的数组，列向量表示RGBA，每个值的范围从 0 到 255

'''colormaps'''import numpy as npfrom mayavi import mlab# 建立数据x, y = np.mgrid[-10:10:200j, -10:10:200j]z = 100 * np.sin(x*y) / (x*y)# 对数据进行可视化mlab.figure(bgcolor=(1,1,1))surf = mlab.surf(z, colormap='cool')mlab.show()

改变 colormap（上面代码，左图）；改变 LUT（下面代码，右图）：

'''LUT'''import numpy as npfrom mayavi import mlab# 建立数据x, y = np.mgrid[-10:10:200j, -10:10:200j]z = 100 * np.sin(x*y) / (x*y)# 对数据进行可视化mlab.figure(bgcolor=(1,1,1))surf = mlab.surf(z, colormap='cool')# 访问surf对象的LUTlut = surf.module_manager.scalar_lut_manager.lut.table.to_array()# 增加透明梯度，修改alpha通道lut[:, -1] = np.linspace(0, 255, 256)surf.module_manager.scalar_lut_manager.lut.table = lutmlab.show()

（三）mlab 控制函数

图像控制函数

图像装饰函数

相机控制函数

其他控制函数

（四）鼠标选取交互操作

选取一个物体，查看数据选取物体上一点，查看数据

1.on_mouse_pick()函数

on_mouse_pick(callback, type='point', Button='Left', Remove=False)

type：'point'、'cell'、'world'Button：'left'、'Middle'、'Right'Remove：如果值为 True，则 callback 函数不起作用返回：一个 vtk picker 对象

2. 选取红色小球问题分析

建立一个 figure随机生成红、白小球初始化红色小球选取外框鼠标选取任意红色小球，外框移动到该小球上 (callback)建立 on_mouse_pick() 响应机制

3. 代码

'''鼠标选取'''import numpy as npfrom mayavi import mlab####### 场景初始化 #######figure = mlab.gcf()# 用mlab.points3d建立红色小球的集合x1, y1, z1 = np.random.random((3, 10))red_glyphs = mlab.points3d(x1, y1, z1, color=(1, 0, 0), resolution=10)# 用mlab.points3d建立白色小球的集合x2, y2, z2 = np.random.random((3, 10))white_glyphs = mlab.points3d(x2, y2, z2, color=(0.9, 0.9, 0.9), resolution=10)# 绘制选取框，并放在第一个小球上outline = mlab.outline(line_width=3)outline.outline_mode = 'cornered'outline.bounds = (x1[0] - 0.1, x1[0] + 0.1,y1[0] - 0.1, y1[0] + 0.1,z1[0] - 0.1, z1[0] + 0.1)####### 处理选取事件 ######## 获取构成一个红色小球的顶点列表glyph_points = red_glyphs.glyph.glyph_source.glyph_source.output.points.to_array()# 处理选取时间def picker_callback(picker):if picker.actor in red_glyphs.actor.actors:# 计算被选取的小球的ID号point_id = int(picker.point_id / glyph_points.shape[0]) # int向下取整 if point_id != -1: # 如果没有小球被选取，则point_id = -1# 计算与此红色小球相关的坐标x, y, z = x1[point_id], y1[point_id], z1[point_id]# 将外框移到小球上outline.bounds = (x - 0.1, x + 0.1,y - 0.1, y + 0.1,z - 0.1, z + 0.1)# 建立响应机制picker = figure.on_mouse_pick(picker_callback)# 设置窗口的标题并显示窗口mlab.title('Click on red balls')mlab.show()

4. 程序框架的优化

程序运行有两个问题：小球初始速度太慢、鼠标选取不精确加入以下三条语句：

figure.scene.disable_render = Truefigure.scene.disable_render = FalsePicker.tolerance = 0.01

5. 优化代码

'''鼠标选取 优化'''import numpy as npfrom mayavi import mlab####### 场景初始化 #######figure = mlab.gcf()figure.scene.disable_render = True# 用mlab.points3d建立红色小球的集合x1, y1, z1 = np.random.random((3, 10))red_glyphs = mlab.points3d(x1, y1, z1, color=(1, 0, 0), resolution=10)# 用mlab.points3d建立白色小球的集合x2, y2, z2 = np.random.random((3, 10))white_glyphs = mlab.points3d(x2, y2, z2, color=(0.9, 0.9, 0.9), resolution=10)figure.scene.disable_render = False# 绘制选取框，并放在第一个小球上outline = mlab.outline(line_width=3)outline.outline_mode = 'cornered'outline.bounds = (x1[0] - 0.1, x1[0] + 0.1,y1[0] - 0.1, y1[0] + 0.1,z1[0] - 0.1, z1[0] + 0.1)####### 处理选取事件 ######## 获取构成一个红色小球的顶点列表glyph_points = red_glyphs.glyph.glyph_source.glyph_source.output.points.to_array()# 处理选取时间def picker_callback(picker):if picker.actor in red_glyphs.actor.actors:# 计算被选取的小球的ID号point_id = int(picker.point_id / glyph_points.shape[0]) # int向下取整 if point_id != -1: # 如果没有小球被选取，则point_id = -1# 计算与此红色小球相关的坐标x, y, z = x1[point_id], y1[point_id], z1[point_id]# 将外框移到小球上outline.bounds = (x - 0.1, x + 0.1,y - 0.1, y + 0.1,z - 0.1, z + 0.1)# 建立响应机制picker = figure.on_mouse_pick(picker_callback)# 设置 tolerance 参数，提高选取精度Picker.tolerance = 0.01# 设置窗口的标题并显示窗口mlab.title('Click on red balls')mlab.show()

6. 结果

（五）mlab 管线控制函数

mlab.pipeline.function()

Sources：数据源Filters：用来数据变换Modules：用来实现可视化

Sources：

Filters：

Modules：

（六）实例1：标量数据可视化

1. 绘制等值面

使用mlab.contour3d(s)绘制等值面，只能看到最外层的等值面（左图）

'''实例1：标量数据可视化 等值面绘制'''import numpy as npfrom mayavi import mlab# 生成标量数据x, y, z = np.ogrid[-10:10:20j, -10:10:20j, -10:10:20j]s = np.sin(x*y*z)/(x*y*z)# 对数据进行可视化mlab.contour3d(s)mlab.show()

2. 绘制切平面

'''实例1：标量数据可视化 切平面绘制'''from mayavi import mlabfrom mayavi.tools import pipelinemlab.pipeline.image_plane_widget(mlab.pipeline.scalar_field(s),plane_orientation='x_axes',slice_index=10,)mlab.pipeline.image_plane_widget(mlab.pipeline.scalar_field(s),plane_orientation='y_axes',slice_index=10,)mlab.outline()mlab.show()

3. 复合观测方法

'''实例1：标量数据可视化 复合观测方法'''from mayavi import mlabfrom mayavi.tools import pipelinesrc = mlab.pipeline.scalar_field(s)mlab.pipeline.iso_surface(src, contours=[s.min()+0.1*s.ptp(), ], opacity=0.1)mlab.pipeline.iso_surface(src, contours=[s.max()-0.1*s.ptp(), ])mlab.pipeline.image_plane_widget(src,plane_orientation='z_axes',slice_index=10,)mlab.show()

（七）实例2：矢量数据可视化

1. 使用quiver3d()函数

简单，但结果非常密集

'''实例2：矢量数据可视化 quiver3d()'''import numpy as npfrom mayavi import mlabx, y, z = np.mgrid[0:1:20j, 0:1:20j, 0:1:20j]u = np.sin(np.pi*x) * np.cos(np.pi*z)v = -2*np.sin(np.pi*y) * np.cos(2*np.pi*z)w = np.cos(np.pi*x)*np.sin(np.pi*z) + np.cos(np.pi*y)*np.sin(2*np.pi*z)mlab.quiver3d(u,v,w)mlab.outline()mlab.show()

2. masking vector 采样

降低采样密度

'''实例2：矢量数据可视化 masking vector'''import numpy as npfrom mayavi import mlabx, y, z = np.mgrid[0:1:20j, 0:1:20j, 0:1:20j]u = np.sin(np.pi*x) * np.cos(np.pi*z)v = -2*np.sin(np.pi*y) * np.cos(2*np.pi*z)w = np.cos(np.pi*x)*np.sin(np.pi*z) + np.cos(np.pi*y)*np.sin(2*np.pi*z)src = mlab.pipeline.vector_field(u, v, w)mlab.pipeline.vectors(src, mask_points=10, scale_factor=2.0)mlab.show()

3. Cut Plane 切面

'''实例2：矢量数据可视化 cut plane'''import numpy as npfrom mayavi import mlabx, y, z = np.mgrid[0:1:20j, 0:1:20j, 0:1:20j]u = np.sin(np.pi*x) * np.cos(np.pi*z)v = -2*np.sin(np.pi*y) * np.cos(2*np.pi*z)w = np.cos(np.pi*x)*np.sin(np.pi*z) + np.cos(np.pi*y)*np.sin(2*np.pi*z)src = mlab.pipeline.vector_field(u, v, w)mlab.pipeline.vector_cut_plane(src, mask_points=10, scale_factor=2.0)mlab.show()

4. 级数的等值面

显示矢量的法向等值面

'''实例2：矢量数据可视化 级数的等值面'''import numpy as npfrom mayavi import mlabx, y, z = np.mgrid[0:1:20j, 0:1:20j, 0:1:20j]u = np.sin(np.pi*x) * np.cos(np.pi*z)v = -2*np.sin(np.pi*y) * np.cos(2*np.pi*z)w = np.cos(np.pi*x)*np.sin(np.pi*z) + np.cos(np.pi*y)*np.sin(2*np.pi*z)src = mlab.pipeline.vector_field(u, v, w)magnitude = mlab.pipeline.extract_vector_norm(src)mlab.pipeline.iso_surface(magnitude, contours=[2.0,0.5])mlab.outline()mlab.show()

5. 流线 flow 的可视化

'''实例2：矢量数据可视化 流线flow的可视化'''import numpy as npfrom mayavi import mlabx, y, z = np.mgrid[0:1:20j, 0:1:20j, 0:1:20j]u = np.sin(np.pi*x) * np.cos(np.pi*z)v = -2*np.sin(np.pi*y) * np.cos(2*np.pi*z)w = np.cos(np.pi*x)*np.sin(np.pi*z) + np.cos(np.pi*y)*np.sin(2*np.pi*z)flow = mlab.flow(u, v, w, seed_scale=1, seed_resolution=5, integration_direction='both')mlab.outline()mlab.show()

6. 复合观测方法

'''实例2：矢量数据可视化 复合观测方法'''import numpy as npfrom mayavi import mlabx, y, z = np.mgrid[0:1:20j, 0:1:20j, 0:1:20j]u = np.sin(np.pi*x) * np.cos(np.pi*z)v = -2*np.sin(np.pi*y) * np.cos(2*np.pi*z)w = np.cos(np.pi*x)*np.sin(np.pi*z) + np.cos(np.pi*y)*np.sin(2*np.pi*z)src = mlab.pipeline.vector_field(u, v, w)magnitude = mlab.pipeline.extract_vector_norm(src)iso = mlab.pipeline.iso_surface(magnitude, contours=[2.0, ], opacity=0.3)vec = mlab.pipeline.vectors(magnitude, mask_points=40, line_width=1,color=(0.8,0.8,0.8),scale_factor=4.)flow = mlab.pipeline.streamline(magnitude, seedtype='plane',seed_visible=False,seed_scale=0.5,seed_resolution=1,linetype='ribbon')vcp = mlab.pipeline.vector_cut_plane(magnitude, mask_points=2, scale_factor=4,colormap='jet',plane_orientation='x_axes')mlab.outline()mlab.show()