rust的学习笔记

电气化、自动化、数字化、智能化、智慧化

0%

ImagePy解析: 28 -- 三维可视化

本文解析一下ImagePy的三维画布。
以如下例子入手:
demo
首先,原始图像是一个5乘5的方形图像,其中间是4乘4的白色,周围是一圈黑色。
由这张原始图根据距离变换得到右上角的高程图,继而对该高程图做三维可视化。

渲染插件

二维平面的三维可视化插件是这样写的:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
class Surface2D(Simple):
title = '2D Surface'
note = ['8-bit', '16-bit', 'float']
para = {'name':'undifine', 'sample':2, 'sigma':2,'h':0.3, 'cm':'gray'}
view = [(str, 'name', 'Name', ''),
(int, 'sample', (1,10), 0, 'down sample', 'pix'),
(int, 'sigma', (0,30), 0, 'sigma', ''),
(float, 'h', (0.1,10), 1, 'scale z', ''),
('cmap', 'cm', 'color map')]

def run(self, ips, imgs, para = None):
ds, sigma, cm = para['sample'], para['sigma'], ColorManager.get(para['cm'])
mesh = Surface2d(ips.img, sample=ds, sigma=sigma, k=para['h'], cmap=cm)
self.app.show_mesh(mesh, para['name'])

其界面为:
interface
即设定名字、下采样率、平滑率和z轴伸缩率,以及渲染所用的colormap。(在该demo中,就按如图中的参数进行设置)
然后将这些参数传给Surface2d这个ImagePy定义的Mesh对象。
最后调用show_mesh方法将其呈现出来。
下面是一步步分析这个Mesh对象及其绘制方法。

Mesh对象

如前所述,二维高程图传给了Surface2d这一类,具体看一下其代码实现:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
class Surface2d(Mesh):
def __init__(self, img=None, sample=1, sigma=0, k=0.3, **key):
self.img, self.sample, self.sigma, self.k = img, sample, sigma, k
Mesh.__init__(self, **key)
self.set_data(img, sample, sigma, k)

def set_data(self, img=None, sample=None, sigma=None, k=None, **key):
if not img is None: self.img = img
if not sample is None: self.sample = sample
if not sigma is None: self.sigma = sigma
if not k is None: self.k = k
if sum([not i is None for i in (img, sample, sigma, k)])>0:
from ..util import meshutil
vert, fs = meshutil.create_surface2d(self.img, self.sample, self.sigma, self.k)
Mesh.set_data(self, verts=vert, faces=fs.astype(np.uint32), colors=vert[:,2], **key)
else: Mesh.set_data(self, **key)

可以看到,在它的初始化函数中调用了set_data方法。进一步地,在该方法中有两个核心方法:将图像转化为顶点和面,然后再转为Mesh对象。

位图提取格点坐标和像素值

即如下方法:

1
vert, fs = meshutil.create_surface2d(self.img, self.sample, self.sigma, self.k)

源码及注释为:
1
2
3
4
5
6
7
8
9
10
11
12
13
def create_surface2d(img, sample=1, sigma=0, k=0.3):
from scipy.ndimage import gaussian_filter
#start = time()
# 以采样率为步长进行图像的重新提取
img = img[::sample, ::sample].astype(np.float32)
# 如果指定了平滑率,则使用高斯滤波进行平滑
if sigma>0: img = gaussian_filter(img, sigma)
# 根据采样后的图像形状生成网格格点
xs, ys = np.mgrid[:img.shape[0],:img.shape[1]]
# 根据采样率,将格点范围伸缩到之前的大小
xs *= sample; ys *= sample
# 将图像像素值乘以伸缩大小,作为z轴的值,与格点坐标xy传入下面的方法
return create_grid_mesh(xs, ys, img*k)

在此例中,依照上面的参数,来看一下各个中间结果:
首先高程图在降采样后,图像矩阵为:
1
2
3
[[0. 0. 0.]
[0. 2. 0.]
[0. 0. 0.]]

然后在$\sigma=1$的高斯滤波后,图像矩阵为:
1
2
3
[[0.17534617 0.2415005  0.17534617]
[0.2415005 0.3326134 0.2415005 ]
[0.17534617 0.2415005 0.17534617]]

其再经过k倍的伸缩,变为:
1
2
3
[[0.8767308 1.2075025 0.8767308]
[1.2075025 1.6630671 1.2075025]
[0.8767308 1.2075025 0.8767308]]

同时xs和ys即网格格点坐标,也经过了降采样,以及范围伸缩,变为:
1
2
3
[[0 0 0]
[2 2 2]
[4 4 4]]

获取格点和面的信息

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
def create_grid_mesh(xs, ys, zs):
h, w = xs.shape
# 将xy坐标位置和z值合并起来
vts = np.array([xs, ys, zs], dtype=np.float32)
# 这一步是定义以某格点为参考点的坐标系下它与哪些点形成面
# 在局部坐标系下,参考点索引为0,那么它所构成的面有两个
# 分别与(1, 1+w)这两个点构成一个面,与(1+w, w)这两个点构成一个面
# 比如在此例下,did的值就是[[0 1 4 0 4 3]]
did = np.array([[0, 1, 1+w, 0, 1+w, w]], dtype=np.uint32)
# rcs由两部分构成
# 第一部分是获取全局坐标系下每一排的第一个元素的索引,所以是以w为步长
# 注意排除最后一排,即是w*h还要减去w,因为最后一排元素所参与的面可以通过倒数第二排来获得
# 对于此例,就是[[0], [3]],注意这里使用None来增加一个维度
# 第二部分是获取全局坐标系下每一列的索引,所以是以1为步长
# 注意是排除最后一列,即w要减去1,这也是因为最后一列参与的面可以通过前一列得到
# 对于此例,就是[0, 1]
# 最终rcs就是numpy数组的[[0],[3]]+[0, 1]
# 这里用到了numpy的广播: https://www.rust-cn.com/2019/10/20/python-indexing/
# [[0, 0],[3, 3]] + [[0, 1], [0, 1]] = [[0, 1], [3, 4]]
# 代表的意思就是在全局坐标系中,第一排取索引为0和1的格点,在第二排取索引为3和4的格点
rcs = np.arange(0,w*h-w,w)[:,None] + np.arange(0,w-1,1)
# 接下来就是根据上面取得的格点,得到每个格点上所形成的面
# 首先第一部分是将rcs拉直为[[0], [1], [3], [4]],即这四个格点索引拉平到一个维度上
# 然后加上上面的局部坐标系下形成面的格点索引[[0 1 4 0 4 3]]
# 同样根据广播原则,就得到了每个格点与相邻点所形成的面
# 结果为:[[0 1 4 0 4 3], [1 2 5 1 5 4], [3 4 7 3 7 6], [4 5 8 4 8 7]]
# 即每个格点上都参与形成两个面,具体每一个面的格点组成看上面的序列
faces = rcs.reshape(-1,1) + did
# 返回值是两个
# 第一个就是格点坐标及其上面的值,并按这三个值合并起来算一个重新改变形状
# 第二个就是由格点所形成的面的信息
return vts.reshape(3,-1).T.copy(), faces.reshape(-1,3)

具体的解析过程见上面源码。
最后说一下最终返回的格点信息和面信息,分别是:

1
2
3
4
5
6
7
8
9
[[0.        0.        0.8767308]
[0. 2. 1.2075025]
[0. 4. 0.8767308]
[2. 0. 1.2075025]
[2. 2. 1.6630671]
[2. 4. 1.2075025]
[4. 0. 0.8767308]
[4. 2. 1.2075025]
[4. 4. 0.8767308]]

以第一个格点为例,它是在(0, 0)坐标,同时上面的值是0.8767308。
以及面信息:
1
2
3
4
5
6
7
8
[[0 1 4]
[0 4 3]
[1 2 5]
[1 5 4]
[3 4 7]
[3 7 6]
[4 5 8]
[4 8 7]]

以第一个面为例,它由(0, 1, 4)号格点组成。

构建Mesh对象

1
2
# 传入上面的格点、面、颜色(这里取的是格点上的z值)以及cmap
Mesh.set_data(self, verts=vert, faces=fs.astype(np.uint32), colors=vert[:,2], **key)

上述代码是调用了Mesh对象的set_data方法。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
class Mesh:
# 在初始化函数中传入一个Mesh对象所需要的信息
def __init__(self, verts=None, faces=None, colors=None, cmap=None, **key):
# 如果有格点信息,但没有面信息
if faces is None and not verts is None:
# 则直接按格点个数-1生成面,即两个相邻格点相连,就成为面
faces = np.arange(len(verts), dtype=np.uint32)
# 传入格点信息
self.verts = verts.astype(np.float32, copy=False) if not verts is None else None
# 传入面信息
self.faces = faces.astype(np.uint32, copy=False) if not faces is None else None
# 传入颜色信息
self.colors = colors
# 设置模式、可见性和dirty属性等
self.mode, self.visible, self.dirty = 'mesh', True, 'geom'
# 设置alpha透明度和边信息
self.alpha = 1; self.edges = None
# 设置高光、colormap
self.high_light = False; self.cmap = 'gray' if cmap is None else cmap
# 调用set_data方法
self.set_data(**key)

def set_data(self, verts=None, faces=None, colors=None, **key):
# 同上面的初始化功能近似,区别是可以直接调用它来配置信息
if faces is None and not verts is None:
faces = np.arange(len(verts), dtype=np.uint32)
if not verts is None: self.verts = verts.astype(np.float32, copy=False)
if not faces is None: self.faces = faces.astype(np.uint32, copy=False)
if not colors is None: self.colors = colors
if not faces is None: self.edge = None
if sum([i is None for i in [verts, faces, colors]])<3: self.dirty = 'geom'
if not self.faces is None and self.faces.ndim==1: key['mode'] = 'points'
elif not self.faces is None and self.faces.shape[1]==2:
if key.get('mode', self.mode)=='mesh': key['mode'] = 'grid'
if key.get('mode', self.mode) != self.mode: self.dirty = 'geom'
self.mode = key.get('mode', self.mode)
self.visible = key.get('visible', self.visible)
self.alpha = key.get('alpha', self.alpha)
self.high_light = key.get('high_light', False)
self.cmap = key.get('cmap', self.cmap)
self.dirty = self.dirty or True

可视化Mesh

即将Mesh对象通过三维画布展示出来:

1
self.app.show_mesh(mesh, para['name'])

这里就是调用了app的show_mesh方法。
ImagePy的三维画布是基于VisPy的,同时又进行了封装,最底层的是如下这个类:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
class Canvas3D(scene.SceneCanvas):
def __new__(cls, parent, scene3d=None):
self = super().__new__(cls)
scene.SceneCanvas.__init__(self, app="wx", parent=parent, keys='interactive', show=True, dpi=150)
canvas = parent.GetChildren()[-1]
self.unfreeze()
self.canvas = weakref.ref(canvas)
self.view = self.central_widget.add_view()
self.set_scene(scene3d or Scene())
self.visuals = {}
self.curobj = None
self.freeze()
canvas.Bind(wx.EVT_IDLE, self.on_idle)
canvas.tool = None
canvas.camera = scene.cameras.TurntableCamera(parent=self.view.scene, fov=45, name='Turntable')
canvas.set_camera = self.set_camera
canvas.fit = lambda : self.set_camera(auto=True)
canvas.at = self.at
self.view.camera = canvas.camera
return canvas

VisPy的教程略微有点少,留坑待填。