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
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
| # 实现PCA分析和法向量计算，并加载数据集中的文件进行验证
import open3d as o3d
import os
import numpy as npPointCloud
from pyntcloud import PyntCloud
import argparse
# 功能：计算PCA的函数
# 输入：
# data：点云，NX3的矩阵
# correlation：区分np的cov和corrcoef，不输入时默认为False
# sort: 特征值排序，排序是为了其他功能方便使用，不输入时默认为True
# 输出：
# eigenvalues：特征值
# eigenvectors：特征向量
# 计算点云PCA
def PCA(data, correlation=False, sort=True):
N = data.shape[0] # 获取点云数据中的点的数量
X = data.to_numpy() # 将数据转换为numpy格式
# 归一化处理
mu = np.mean(X, axis=0)
X_normalized = X - mu
# 定义函数，这里我们使用暂不使用核函数的方法，而是计算归一化数据的协方差矩阵
func = np.cov if not correlation else np.corrcoef
H = func(X_normalized, rowvar=False, bias=True)
# 计算协方差矩阵的特征值和特征向量,这里需要注意得到的协方差矩阵维数为3
eigenvalues, eigenvectors = np.linalg.eig(H)
# 对特征值进行排序，特征值最大的就是最主要的成分，用于求取PCA；最小的就是最不重要的成分，用于求取法向量；对三个特征值进行排序后得到其由小到大的索引，
# 用该索引对特征向量排序，则[:,0]就是最小特征值对应的特征向量，即法向量；则[:,2]就是最大的特征值对赢得特征向量，即PCA
if sort:
sort = eigenvalues.argsort()[::-1]
eigenvalues = eigenvalues[sort]
eigenvectors = eigenvectors[:, sort]
return eigenvalues, eigenvectors
# 使用open3d可视化PCA
def get_pca_o3d(w, v, points):
# 中心为点云数据的中心
centroid = points.mean()
# 获取尺度
projs = np.dot(points.to_numpy(), v[:, 0])
scale = projs.max() - projs.min()
#
points = centroid.to_numpy() + np.vstack(
(
np.asarray([0.0, 0.0, 0.0]),
scale * v.T
)
).tolist()
lines = [[0, 1], [0, 2], [0, 3]]
colors = np.identity(3).tolist()
pca_o3d = o3d.geometry.LineSet(
points=o3d.utility.Vector3dVector(points),
lines=o3d.utility.Vector2iVector(lines),
)
pca_o3d.colors = o3d.utility.Vector3dVector(colors)
return pca_o3d
def get_surface_normals(pcd, points, knn=5):
# 设定点集搜索方法
pcd_tree = o3d.geometry.KDTreeFlann(pcd)
# 点的数量
N = len(pcd.points)
# 法向量容器
normals = []
# 分别计算每个点的法线
for i in range(N):
# 获取当前点的最临近5个点
[k, idx, _] = pcd_tree.search_knn_vector_3d(pcd.points[i], knn)
# 获取特征值和特征向量
w, v = PCA(points.iloc[idx])
# 得到特征向量的矩阵v，拿到临近点集协方差矩阵的特征值最小的特征向量，即该点的特征向量
normals.append(v[:, 0])
return np.array(normals, dtype=np.float64)
def get_surface_normal_o3d(normals, points, scale=2):
N = points.shape[0]
points = np.vstack(
(points.to_numpy(), points.to_numpy() + scale * normals)
)
lines = [[i, i + N] for i in range(N)]
colors = np.zeros((N, 3)).tolist()
surface_normals_o3d = o3d.geometry.LineSet(
points=o3d.utility.Vector3dVector(points),
lines=o3d.utility.Vector2iVector(lines),
)
surface_normals_o3d.colors = o3d.utility.Vector3dVector(colors)
return surface_normals_o3d
def main(point_cloud_filename):
# 加载原始点云
point_cloud_pynt = PyntCloud.from_file(point_cloud_filename)
point_cloud_o3d = point_cloud_pynt.to_instance("open3d", mesh=False)
# o3d.visualization.draw_geometries([point_cloud_o3d]) # 显示原始点云
# 从点云中获取点，只对点进行处理
points = point_cloud_pynt.points
print('total points number is:', points.shape[0])
# 用PCA分析点云主方向
w, v = PCA(points)
point_cloud_vector = v[:, 2] # 点云主方向对应的向量
print('the main orientation of this pointcloud is: ', point_cloud_vector)
# TODO: 此处只显示了点云，还没有显示PCA
# o3d.visualization.draw_geometries([point_cloud_o3d])
pca_o3d = get_pca_o3d(w, v, points)
# 循环计算每个点的法向量
normals = get_surface_normals(point_cloud_o3d, points)
point_cloud_o3d.normals = o3d.utility.Vector3dVector(normals)
surface_normals_o3d = get_surface_normal_o3d(normals, points)
o3d.visualization.draw_geometries([point_cloud_o3d, pca_o3d,surface_normals_o3d])
def get_arguments():
""" Get command-line arguments
"""
# init parser:
parser = argparse.ArgumentParser("Get PCA and surface normals for given point cloud.")
# add required and optional groups:
required = parser.add_argument_group('Required')
# add required:
required.add_argument(
"-i", dest="input", help="Input path of point cloud in ply format.",
required=True
)
# parse arguments:
return parser.parse_args()
if __name__ == '__main__':
arguments = get_arguments()
main(arguments.input)
|
