三维点云学习(1)上-PCA主成分分析 法向量估计
三维点云学习(1)上
环境安装
1.系统环境 win10 或者 ubuntu
2. Anaconda3+python3.6
使用Anaconda创建的conda虚拟环境进行python的编写
环境安装主要参考如下网址
安装Anaconda3
Anconda3 安装 open3d
3. 使用conda install 或者 pip install 下载需要的py模块
open3d numpy matplotlib pandas plyfile pyntcloud
#offto_ply.py
import os
import numpy as np
from plyfile import PlyData
from plyfile import PlyElement
#pca_normal.py
import open3d as o3d
import os
import numpy as np
import matplotlib.pyplot as plt
#voxel_filter.py
import open3d as o3d
import os
import numpy as np
from pyntcloud import PyntCloud
4.数据集下载
为40种物体的三维点云数据集
链接:https://pan.baidu.com/s/1LX9xeiXJ0t-Fne8BCGSjlQ<br> 提取码:es14
5.单独读取数据集,并在open3d中显示原点云图
注意:要把对应的数据集文件 如:"sofa_0001.txt"放在对应的代码路径下
读取方法一(这个读法忽略了后三个为法向量,因而处理错误,推荐方法二)
import open3d as o3d #导入open3d
import numpy as np
import matplotlib as plt
raw_point_cloud_matrix = np.genfromtxt(r"sofa_0001.txt", delimiter=",").reshape((-1,3))
# print(raw_point_cloud_matrix)
pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(raw_point_cloud_matrix)
print(pcd)
o3d.visualization.draw_geometries([pcd])
原始点云数据在open3d下运行结果如下所示:
读取方法二
import open3d as o3d
import os
import numpy as np
import matplotlib.pyplot as plt
from pandas import DataFrame
from pyntcloud import PyntCloudpoint_cloud_raw = np.genfromtxt(r"plant_0001.txt", delimiter=",") #为 xyz的 N*3矩阵
point_cloud_raw = DataFrame(point_cloud_raw[:, 0:3]) # 选取每一列 的 第0个元素到第二个元素 [0,3)
point_cloud_raw.columns = ['x', 'y', 'z'] # 给选取到的数据 附上标题
point_cloud_pynt = PyntCloud(point_cloud_raw) # 将points的数据 存到结构体中point_cloud_o3d = point_cloud_pynt.to_instance("open3d", mesh=False) # 实例化o3d.visualization.draw_geometries([point_cloud_o3d]) # 显示原始点云
6.便捷的显示模块摘自enginelong的博客代码片
# matplotlib显示点云函数
def Point_Cloud_Show(points):fig = plt.figure(dpi=150)ax = fig.add_subplot(111, projection='3d')ax.scatter(points[:, 0], points[:, 1], points[:, 2], cmap='spectral', s=2, linewidths=0, alpha=1, marker=".")plt.title('Point Cloud')ax.set_xlabel('x')ax.set_ylabel('y')ax.set_zlabel('z')plt.show()
# 二维点云显示函数
def Point_Show(pca_point_cloud):x = []y = []pca_point_cloud = np.asarray(pca_point_cloud)for i in range(10000):x.append(pca_point_cloud[i][0])y.append(pca_point_cloud[i][1])plt.scatter(x, y)plt.show()
7.PCA 主成分分析法
参考公式网址如下所示
三维点云处理学习笔记
PCA原理解释
对协方差矩阵的通俗理解
代码参考网址如下所示
点云学习(1)参考网址1
点云学习(1)参考网址2
PCA算法的优化目标
1.降维后同一纬度的方差最大
2.不同维度之间的相关性为0
PCA函数编写
主要流程
1.取均值,去中心化
average_data = np.mean(data,axis=0) #求 NX3 向量的均值decentration_matrix = data - average_data #去中心化
2.求取协方差矩阵H,并用SVD奇异值分解,求解出相应的特征值、特征向量
H = np.dot(decentration_matrix.T,decentration_matrix) #求解协方差矩阵 Heigenvectors,eigenvalues,eigenvectors_T = np.linalg.svd(H) # SVD求解特征值、特征向量
3.对2步求解出的特征值特征向量进行降序排序,并存放到列表中
if sort:sort = eigenvalues.argsort()[::-1] #降序排列eigenvalues = eigenvalues[sort] #索引eigenvectors = eigenvectors[:, sort]
PCA处理整体代码块:
# 功能:计算PCA的函数
# 输入:
# data:点云,NX3的矩阵
# correlation:区分np的cov和corrcoef,不输入时默认为False
# sort: 特征值排序,排序是为了其他功能方便使用,不输入时默认为True
# 输出:
# eigenvalues:特征值
# eigenvectors:特征向量
def PCA(data, correlation=False, sort=True):# 作业1# 屏蔽开始average_data = np.mean(data,axis=0) #求 NX3 向量的均值decentration_matrix = data - average_data #去中心化H = np.dot(decentration_matrix.T,decentration_matrix) #求解协方差矩阵 Heigenvectors,eigenvalues,eigenvectors_T = np.linalg.svd(H) # SVD求解特征值、特征向量# 屏蔽结束if sort:sort = eigenvalues.argsort()[::-1] #降序排列eigenvalues = eigenvalues[sort] #索引eigenvectors = eigenvectors[:, sort]return eigenvalues, eigenvectors
4.调用结果
通过调用PCA算法,并显示两个主成分方向,如下图所示黑色线为第一主成分,红色线为第二主成分
5.PCA的应用
降维(Encoder)
#将原数据进行降维度处理point_cloud_encode = (np.dot(point_cloud_vector.T,point_cloud_raw.T)).T #主成分的转置 dot 原数据Point_Show(point_cloud_encode)
效果如下所示:
降维后效果图
升维(Decoder)
#使用主方向进行升维point_cloud_decode = (np.dot(point_cloud_vector,point_cloud_encode.T)).TPoint_Cloud_Show(point_cloud_decode)
效果如下所示:
使用两个主方向再次升维后的结果
8.法向量估计
法向量估计
思想:选取点云中每一点,对其进行临近点的搜索,将包含该点的临近点拟合成曲面,对曲面中的点进行PCA主成分分析,查找特征值最小的对应的特征向量,该特征向量即为该拟合曲面的法向量(摘自:秦乐乐CSDN博客)
1.编码流程
2.代码展示
# 循环计算每个点的法向量pcd_tree = o3d.geometry.KDTreeFlann(point_cloud_o3d) #将原始点云数据输入到KD,进行近邻取点normals = [] #储存曲面的法向量# 作业2# 屏蔽开始print(points.shape[0]) #打印当前点数 20000个点for i in range(points.shape[0]):# search_knn_vector_3d函数 , 输入值[每一点,x] 返回值 [int, open3d.utility.IntVector, open3d.utility.DoubleVector][
,idx,] = pcd_tree.search_knn_vector_3d(point_cloud_o3d.points[i],10) #取10个临近点进行曲线拟合# asarray和array 一样 但是array会copy出一个副本,asarray不会,节省内存k_nearest_point = np.asarray(point_cloud_o3d.points)[idx, :] # 找出每一点的10个临近点,类似于拟合成曲面,然后进行PCA找到特征向量最小的值,作为法向量w, v = PCA(k_nearest_point)normals.append(v[:, 2])# 屏蔽结束normals = np.array(normals, dtype=np.float64)# TODO: 此处把法向量存放在了normals中point_cloud_o3d.normals = o3d.utility.Vector3dVector(normals)o3d.visualization.draw_geometries([point_cloud_o3d])
3.效果展示
如下图1所示为图2的法向量加粗展示,
法线向量显示:在显示窗口按n 可按 + - 更改点的大小(o3d)
完整代码(包含PCA算法、PCA应用(升降维)、法向量估计展示)
# 实现PCA分析和法向量计算,并加载数据集中的文件进行验证import open3d as o3d
import os
import numpy as np
import matplotlib.pyplot as plt
from pandas import DataFrame
from pyntcloud import PyntCloud# matplotlib显示点云函数
def Point_Cloud_Show(points):fig = plt.figure(dpi=150)ax = fig.add_subplot(111, projection='3d')ax.scatter(points[:, 0], points[:, 1], points[:, 2], cmap='spectral', s=2, linewidths=0, alpha=1, marker=".")plt.title('Point Cloud')ax.set_xlabel('x')ax.set_ylabel('y')ax.set_zlabel('z')plt.show()# 二维点云显示函数
def Point_Show(pca_point_cloud):x = []y = []pca_point_cloud = np.asarray(pca_point_cloud)for i in range(10000):x.append(pca_point_cloud[i][0])y.append(pca_point_cloud[i][1])plt.scatter(x, y)plt.show()# 功能:计算PCA的函数
# 输入:
# data:点云,NX3的矩阵
# correlation:区分np的cov和corrcoef,不输入时默认为False
# sort: 特征值排序,排序是为了其他功能方便使用,不输入时默认为True
# 输出:
# eigenvalues:特征值
# eigenvectors:特征向量
def PCA(data, correlation=False, sort=True):# 作业1# 屏蔽开始average_data = np.mean(data,axis=0) #求 NX3 向量的均值decentration_matrix = data - average_data #去中心化H = np.dot(decentration_matrix.T,decentration_matrix) #求解协方差矩阵 Heigenvectors,eigenvalues,eigenvectors_T = np.linalg.svd(H) # SVD求解特征值、特征向量# 屏蔽结束if sort:sort = eigenvalues.argsort()[::-1] #降序排列eigenvalues = eigenvalues[sort] #索引eigenvectors = eigenvectors[:, sort]return eigenvalues, eigenvectorsdef main():# 指定点云路径# cat_index = 10 # 物体编号,范围是0-39,即对应数据集中40个物体# root_dir = '/Users/renqian/cloud_lesson/ModelNet40/ply_data_points' # 数据集路径# cat = os.listdir(root_dir)# filename = os.path.join(root_dir, cat[cat_index],'train', cat[cat_index]+'_0001.ply') # 默认使用第一个点云# 加载原始点云,txt处理point_cloud_raw = np.genfromtxt(r"/home/renzhanqi/code/pythonCode/shenLanLidarPrcess/data/modelnet40_normal_resampled/plant/plant_0001.txt", delimiter=",") #为 xyz的 N*3矩阵point_cloud_raw = DataFrame(point_cloud_raw[:, 0:3]) # 选取每一列 的 第0个元素到第二个元素 [0,3)point_cloud_raw.columns = ['x', 'y', 'z'] # 给选取到的数据 附上标题point_cloud_pynt = PyntCloud(point_cloud_raw) # 将points的数据 存到结构体中point_cloud_o3d = point_cloud_pynt.to_instance("open3d", mesh=False) # 实例化o3d.visualization.draw_geometries([point_cloud_o3d]) # 显示原始点云# 从点云中获取点,只对点进行处理print(point_cloud_o3d) #打印点数# 用PCA分析点云主方向w, v = PCA(point_cloud_raw) # w为特征值 v为主方向point_cloud_vector1 = v[:, 0] #点云主方向对应的向量,第一主成分point_cloud_vector2 = v[:, 1] # 点云主方向对应的向量,第二主成分point_cloud_vector = v[:,0:2] # 点云主方向与次方向print('the main orientation of this pointcloud is: ', point_cloud_vector1)print('the main orientation of this pointcloud is: ', point_cloud_vector2)#在原点云中画图point = [[0,0,0],point_cloud_vector1,point_cloud_vector2] #画点:原点、第一主成分、第二主成分lines = [[0,1],[0,2]] #画出三点之间两两连线colors = [[1,0,0],[0,0,0]]#构造open3d中的LineSet对象,用于主成分显示line_set = o3d.geometry.LineSet()line_set.lines = o3d.utility.Vector2iVector(lines)line_set.colors = o3d.utility.Vector3dVector(colors)line_set.points =o3d.utility.Vector3dVector(point)o3d.visualization.draw_geometries([point_cloud_o3d,line_set]) # 显示原始点云和PCA后的连线#将原数据进行降维度处理point_cloud_encode = (np.dot(point_cloud_vector.T,point_cloud_raw.T)).T #主成分的转置 dot 原数据Point_Show(point_cloud_encode)#使用主方向进行升维point_cloud_decode = (np.dot(point_cloud_vector,point_cloud_encode.T)).TPoint_Cloud_Show(point_cloud_decode)# 循环计算每个点的法向量pcd_tree = o3d.geometry.KDTreeFlann(point_cloud_o3d) #将原始点云数据输入到KD,进行近邻取点normals = [] #储存曲面的法向量# 作业2# 屏蔽开始print(point_cloud_raw.shape[0]) #打印当前点数 20000个点for i in range(point_cloud_raw.shape[0]):# search_knn_vector_3d函数 , 输入值[每一点,x] 返回值 [int, open3d.utility.IntVector, open3d.utility.DoubleVector][_,idx,_] = pcd_tree.search_knn_vector_3d(point_cloud_o3d.points[i],10) #取10个临近点进行曲线拟合# asarray和array 一样 但是array会copy出一个副本,asarray不会,节省内存k_nearest_point = np.asarray(point_cloud_o3d.points)[idx, :] # 找出每一点的10个临近点,类似于拟合成曲面,然后进行PCA找到特征向量最小的值,作为法向量w, v = PCA(k_nearest_point)normals.append(v[:, 2])# 屏蔽结束normals = np.array(normals, dtype=np.float64)# TODO: 此处把法向量存放在了normals中point_cloud_o3d.normals = o3d.utility.Vector3dVector(normals)o3d.visualization.draw_geometries([point_cloud_o3d])if __name__ == '__main__':main()
三维点云学习(1)上-PCA主成分分析 法向量估计相关推荐
- 三维点云学习(2)上- 二叉树实现K-NN Radius-NN Search
三维点云学习(2)上 二叉树实现K-NN Radius-NN Search 代码来自 黎老师github 个人心得 二叉树的搜寻方法 正如老师课堂所说,实现二叉树的搜寻有两种方法,一种是递归,一种是循 ...
- 三维点云学习(5)4-实现Deeplearning-PointNet-1-数据集的批量读取
三维点云学习(5)4-实现Deeplearning-PointNet-1-数据集的批量读取 Github PointNet源码 数据集下载:为40种物体的三维点云数据集 提取码:es14 因为本人初次 ...
- 三维点云学习(4)4-Hough Transform
三维点云学习(4)4-Hough Transform 霍夫变换的理论通俗理解 Hough Tansform 霍夫变换 核心思想: a.原始空间中得点->参数空间中的线 b.原始空间得线-> ...
- 三维点云学习(4)7-ransac 地面分割+ DBSCAN聚类比较
三维点云学习(4)7-ransac 地面分割+ DBSCAN聚类比较 回顾: 实现ransac地面分割 DBSCNA python 复现-1- 距离矩阵法 DBSCNA python 复现-2-kd- ...
- 三维点云学习(4)6-ransac 地面分割
三维点云学习(4)6-ransac 地面分割 ransac课堂笔记 git大神参考代码 ransac代码主要参考如下知乎大佬的ransac的线性拟合 ransac的线性拟合 使用ransac进行地面分 ...
- 三维点云学习(4)5-DBSCNA python 复现-2-kd-_tree加速
三维点云学习(4)5-DBSCNA python 复现-2-kd-tree加速 因为在上一章DBSCAN在构建距离矩阵时,需要构建一个NN的距离矩阵,严重占用资源,古采用kd_tree搜索进行进一步的 ...
- 三维点云学习(3)7- 实现GMM
三维点云学习(3)7- 实现GMM github大神参考代码 高斯混合模型的通俗理解 GMM课程个人总结笔记 最终效果图 原图 进行高斯聚类后的图 代码编写流程 1.输入数据集x1 x2 -xn,和K ...
- 三维点云学习(3)4-Expectation-Maximization (EM)
三维点云学习(3)4-Expectation-Maximization (EM) EM算法的概述 EM算法中的M-step不断优化Q函数 使用两种方法推导出GMM 1.最大似然法 2.EM算法 K-M ...
- 三维点云学习(3)2- K-Means
三维点云学习(3)2- K-Means K-Means的具体构建流程 将输入的N个数据点,分为N个类 1.随机选取K个中心点 2.E-Step(expectation):N个点.K个中心,求N个点到K ...
最新文章
- 【超值干货】10个案例告诉你,数据如何驱动产品设计
- Python matplotlib 和PIL
- Python3.4 django使用mysql
- Asp.Net MVC4入门指南(9):查询详细信息和删除记录
- 【1】TCP三次握手的第三次的 ack包丢失会怎样?
- 如何开始第一个开源项目?
- 明明管理失败,跟距离远有什么关系?
- 人工智能的未来是强化学习_多主体强化学习与AI的未来
- 企业如何培养新型员工队伍
- 关于Eclipse配置Tomcat8的问题
- 灵格斯与众多常用软件的冲突问题
- Windows开发入门:工具-WinDbg的安装和使用教程
- Nginx基础之错误页面配置/流量控制/访问控制/变量/监控/HTTPS配置/性能优化
- [转载]刘光斗-刘晚苍系武学传人概况
- 阿里云调用api配置access_key
- 数据可视化总结——matplotlib、seaborn
- 美国洛杉矶时间转 格林威治时间
- vscode代码格式管理插件prettier-Code formatter安装和设置
- inter处理器 至强e5与i7有什么区别?Inter(R) Xeon(R) CPU E5-2660 v4@2.0GHz 2.00GHz
- 数据挖掘—逻辑回归分类—信用卡欺诈分析
热门文章
- vue项目打包部署-----解决打包后访问资源失败问题
- VS2017 启动调试出现 无法启动程序“http://localhost:15613” 操作在当前状态中是非法的。 同时附加进程也是错误的解决方法
- mysql5.7 解压版 中文乱码_MySQL 5.7解压版安装、卸载及乱码问题的图文解决方法...
- Aubo i5真机 ros - melodic 版驱动下载 [ 驱动下载 ]
- svn安装使用subversion
- Android:安卓线性布局(属性)
- c java 字节流_Java 学习(23)---(IO流之字节流)
- 案例-图片缩放(CSS3)
- ES6 数组高频使用方法
- mysql最多多少个索引_在一个球的周围,最多能摆放多少个相同尺寸的球在它周围?...