基于open3d重建点云表面

1
import matplotlib.pyplot as plt
2
from aicsimageio import AICSImage
3
from glob import glob
4
import numpy as np
5
import open3d as o3d
6
from sklearn.cluster import DBSCAN
7

8

9
fps = glob("*.czi")
10
a = AICSImage(fps[1])
11
data = a.get_image_data("ZYX", C=0)
12
zs, ys, xs = np.where(data>0)
13
locs = np.vstack([xs, ys, zs]).T
14

15

16
# 第一次聚类用于去噪
17
clustering = DBSCAN(eps=3, min_samples=3).fit(locs)
18
inds = np.where(clustering.labels_>-1)
19
locs2 = locs[inds]
20
# 第二次聚类用于分开不同的大cluster
21
clustering2 = DBSCAN(eps=10, min_samples=3).fit(locs2)
22

23

24
# 取其中一个点云cluster进行尝试
25
inds3 = np.where(clustering2.labels_==3)
26
locs3 = locs2[inds3]
27

28

29
# 转换为o3d的点云对象
30
pt = o3d.geometry.PointCloud()
31
pt.points = o3d.utility.Vector3dVector(locs3)
32

33

34
# 使用 alphashape 计算点云的表面
35
mesh1 = o3d.geometry.TriangleMesh.create_from_point_cloud_alpha_shape(pt, alpha=2)
36
# alpha参数越大，表面越光滑，但可能舍弃的点越多
37
mesh1.compute_vertex_normals()
38
# 查看 alphashape 方法得到的 mesh
39
# o3d.visualization.draw_geometries([pt, mesh1],  width=400, height=300,  mesh_show_back_face=True)
40

41

42
# 基于泊松方法重建点云表面
43
pt2 = mesh1.sample_points_poisson_disk(len(locs3))
44
# disk可以输入一个整数，代表采样密度，如果一个单位体积中点密度高，这个数值就写高点
45
pt2.estimate_normals()
46
mesh2, densities = o3d.geometry.TriangleMesh.create_from_point_cloud_poisson(pt2, depth=3)
47
mesh2.compute_vertex_normals()
48
# 按照密度上伪彩
49
d = np.asarray(densities)
50
dc = plt.get_cmap('plasma')((d-d.min())/(d.max()-d.min()))
51
dc = dc[:,:3]
52
mesh2.vertex_colors = o3d.utility.Vector3dVector(dc)
53
# 查看泊松重建表面效果
54
o3d.visualization.draw_geometries([pt, mesh2], width=400, height=300, mesh_show_back_face=True)
55
# 保存网格为3d文件，支持其它三维作图软件查看
56
o3d.io.write_triangle_mesh("psf.gltf", mesh2)