深蓝学院-----3D点云课程第一章笔记

weixin_47343723

已于 2022-08-07 18:18:11 修改

阅读量930

点赞数 2

分类专栏：点云学习文章标签：深度学习

于 2022-08-06 13:11:38 首次发布

本文链接：https://blog.csdn.net/weixin_47343723/article/details/126192990

版权

点云学习专栏收录该内容

5 篇文章 1 订阅

订阅专栏

主成分分析（PCA）

1.作用

PCA用于找到某一坐标系，使得点云在坐标系下的投影数据最多，即：降维后尽可能小的减少数据损失。
由上述介绍可知：PCA可以应用在数据降维、压缩等方面。

2.操作步骤

计算PCA的流程如下：

计算PCA前需要了解一些知识 （主要是第四条），此处指路b站：用最直观的方式告诉你：什么是主成分分析PCA 介绍的很详细：

D’=RSD
其中，D’是拉伸旋转后的点云数据，S是拉伸矩阵 S=S^T ，R是旋转矩阵 R^-1=R^T。我们希望拉伸旋转后的点云数据在坐标轴上的投影数据损失最小。
协方差矩阵C= DD^T/(n-1)
奇异值分解（SVD）：将一个矩阵分解为三个矩阵的乘积，即：M=ULV^T。
其中，V是M^TM的特征向量，U是MM^T的特征向量，L是MM^T或者M^TM的奇异值（也就是特征值开根号）。
SVD与PCA的关系：V=R。由此可知，计算R时，可以不必求原始数据D的协方差矩阵，而通过计算D^TD的特征向量、特征值求得R、S，大大减少了计算量。注意：特征向量相同，特征值差n-1倍。（此处我感觉是这样，由线代矩阵扩大k倍，特征值扩大k，特征向量不变的定理得出。）

计算主成分有以下两种方法：
1.直接计算D的协方差矩阵C，求R和S。
2.中心化之后，求D^TD，求特征向量即为R，特征值/(n-1)即为S。

【例题】
这里用Matlab测试C与D^TD的关系

# 方法一
a = [2.5,0.5,2.2,1.9,3.1,2.3,2,1,1.5,1.1]; 
b = [2.4,0.7,2.9,2.2,3,2.7,1.6,1.1,1.6,0.9];
cov(a,b)
# 输出
ans =
    0.6166    0.6154
    0.6154    0.7166

# 方法二
# 1.先中心化 （没有中心化得不到上述结果）
a1 = a-mean(a);
b1 = b-mean(b); 
# 2. 拼接a1,b1
a1 = a';
b1 = b';
x = [a1,b1]

x =
    0.6900    0.4900
   -1.3100   -1.2100
    0.3900    0.9900
    0.0900    0.2900
    1.2900    1.0900
    0.4900    0.7900
    0.1900   -0.3100
   -0.8100   -0.8100
   -0.3100   -0.3100
   -0.7100   -1.0100
# 3. 求(xT*x)/(n-1)
x'*x/9
# 输出
ans =
    0.6166    0.6154
    0.6154    0.7166

此处需要注意n的含义 ，n是样本的个数。
计算出主成分之后，降至n维：
取前n个最大的特征值对应的特征向量组成矩阵P，D*P即为降维后的结果（第一列x，第二列y，使用matplotlib画二维散点图）。

【例题】
将1个35的向量降维至32，即：5个特征降到2个。res = np.matmul(nor,v)即为降维。

import numpy as np

mat = [[-1,-1,0,2,1],[2,0,0,-1,-1],[2,0,1,1,0]]
mat = np.array(mat)
avg = np.mean(mat,axis=0)  # 0：列求平均
nor = mat-avg
n = nor.shape[0]
eigenvalues, eigenvectors = np.linalg.eig(nor.T.dot(nor)/(n-1))
sort = eigenvalues.argsort()[::-1]  #表示对特征值eigenvalues 从大到小排序，返回索引值
eigenvalues = eigenvalues[sort]
eigenvectors = eigenvectors[:, sort]

print(eigenvectors.shape)  # (5, 5)
v = eigenvectors[:,0:2]    # (5, 2) 特征向量是按列排

res = np.matmul(nor,v) #！！！！！降维降维！！！！！！
print(res.shape) # (3, 2)
print(res)
#[[ 2.6838453  -0.36098161]
 [-2.09303664 -0.78689112]
 [-0.59080867  1.14787272]]

3.作业

此处有坑！！小白我想了一天才发现，传入的data是个n*3的矩阵，不能加颜色那3个维度，也就是如下：

points = points.loc[:,["x","y","z"]]  #不加颜色那3个维度
w, v, mid = PCA(points)
p = v[:,0:2]
res = np.matmul(mid,p)  # 降维！！
res = np.asarray(res)
# 
plt.scatter(res.loc[:,[0]],res.loc[:,[1]])
plt.show()

def PCA(data, correlation=False, sort=True):
    # 屏蔽开始
    # data的数据类型是pandas,因此加减乘除都需要符合pandas		
    avg = np.mean(data,axis=0)
    # 中心化
    nor = data.sub(avg,axis=1)
    #SVD
    n = nor.shape[0]
    eigenvalues, eigenvectors = np.linalg.eig(nor.T.dot(nor)/(n-1))  #自带函数求解特征值特征向量
    # 屏蔽结束
    if sort:
        sort = eigenvalues.argsort()[::-1]
        eigenvalues = eigenvalues[sort]
        eigenvectors = eigenvectors[:, sort]

    return eigenvalues, eigenvectors，nor

结果展示如下:

在这里插入图片描述

法向量

用open3d包计算并显示法向量

pyntcloud类型转为open3d类型

p3d = point_cloud_pynt.to_instance("open3d", mesh=False)

多种方法计算法向量

p3d.estimate_normals(
        search_param=o3d.geometry.KDTreeSearchParamKNN(knn=15))  # 当前点与15个近邻点构成平面的法向量
p3d.estimate_normals(
        search_param=o3d.geometry.KDTreeSearchParamRadius(radius=0.01))  # 搜索半径1cm
p3d.estimate_normals(
        search_param=o3d.geometry.KDTreeSearchParamHybrid(radius=0.01, max_nn=20))  # 搜索半径1cm，只考虑邻域内的20个点
# p3d应该有默认的法向量，使用上述方法会更改默认法向量        
normals = np.asarray(point_cloud_o3d.normals)	# 将法向量转为numpy类型用于打印和显示

显示法向量*

normals = np.array(normals, dtype=np.float64)
p3d.normals = o3d.utility.Vector3dVector(normals)
o3d.visualization.draw_geometries([point_cloud_o3d],point_show_normal=True) 
# point_show_normal=True用于显示法向量

point_show_normal=False时，显示原始点云

point_show_normal=True时，显示包含法向量的点云 在这里插入图片描述

pyntcloud官网介绍

中心/随机体素滤波(下采样)

在这里插入图片描述

随机滤波

利用如上公式求得每个样本点对应的网格索引值index（每个样本对应一个网格索引值），随机选取相同的网格索引值下的一个点（样本索引值）加入filtered_points列表中，并再列表l（l中存储样本索引值）中删除属于该索引的所有元素，重复上述过程至l清空。根据索引列表filtered_points获得点云的6维特征（使用.select_by_index()函数），并显示。

# 实现voxel滤波，并加载数据集中的文件进行验证

import open3d as o3d
import numpy as np
from pyntcloud import PyntCloud
import random

def voxel_filter(point_cloud, leaf_size):
    filtered_points = []
    # 作业3
    # 屏蔽开始
    uppers = np.max(point_cloud, axis=0)
    lowers = np.min(point_cloud, axis=0)
    # 每个格子的长宽高
    D = (uppers - lowers) / leaf_size

    # 所有点的索引值
    H = (point_cloud - lowers) // leaf_size
    H = np.asarray(H)

    index = H[:, 0] + H[:, 1] * D['x'] + H[:, 2] * D['x'] * D['y']  # 每个样本的索引值

    ind = index.argsort()  # 存放排序后的索引
    l = list(ind)

    while len(l) != 0:
        i = 0
        find_n = np.where(index == index[l[i]])

        # 随机抽取元素
        ran = random.choice(find_n[0])
        filtered_points.append(ran)
        l = list(set(l) - set(find_n[0]))


    # 屏蔽结束

    # 把点云格式改成array，并对外返回
    # filtered_points = np.array(filtered_points, dtype=np.float64)
    # print(filtered_points)
    return filtered_points


def main():
    data_dir = 'D:/BaiduNetdiskDownload/modelnet40_normal_resampled/modelnet40_normal_resampled/'
    with open(data_dir + 'modelnet40_shape_names.txt') as f:
        a = f.readlines()
   
    point_cloud_pynt = PyntCloud.from_file(
        data_dir + 'airplane/airplane_0002.txt', sep=",",
        names=["x", "y", "z", "nx", "ny", "nz"])

    point = point_cloud_pynt.points
    point = point.loc[:, ["x", "y", "z"]]

    # 转成open3d能识别的格式
    point_cloud_o3d = point_cloud_pynt.to_instance("open3d", mesh=False)
    # o3d.visualization.draw_geometries([point_cloud_o3d]) # 显示原始点云

    # 调用voxel滤波函数，实现滤波
    filtered_cloud = voxel_filter(point, 0.1)  # 100

    for_show = point_cloud_o3d.select_by_index(filtered_cloud)


    # 显示滤波后的点云
    o3d.visualization.draw_geometries([for_show])
    o3d.visualization.draw_geometries([point_cloud_o3d])


if __name__ == '__main__':
    main()

中心滤波

# voxel_filter.py 放了全部代码
# 实现voxel滤波，并加载数据集中的文件进行验证

import open3d as o3d 
import os
import numpy as np
from pyntcloud import PyntCloud
import random
import pandas as pd

def voxel_filter(point_cloud, leaf_size):
    point = point_cloud.loc[:, ["x", "y", "z"]]

    filtered_points = pd.DataFrame()
    # 作业3
    # 屏蔽开始
    uppers = np.max(point, axis=0)
    lowers = np.min(point, axis=0)
    #每个格子的长宽高
    D = (uppers-lowers)/leaf_size

    #所有点的索引值
    H = (point-lowers)//leaf_size
    H = np.asarray(H)

    index = H[:,0]+H[:,1]*D['x']+H[:,2]*D['x']*D['y'] #每个样本的索引值

    ind = index.argsort()  #存放排序后的索引
    l = list(ind)

    while len(l)!=0:
        find_n = np.where(index==index[l[0]])

        selected = point_cloud.iloc[find_n[0]]
        p = np.mean(selected)
        filtered_points = pd.concat([filtered_points,p],axis=1)
        l = list(set(l)-set(find_n[0]))

    # 屏蔽结束

    # 把点云格式改成array，并对外返回
    # filtered_points = np.array(filtered_points, dtype=np.float64)
    # print(filtered_points)
    filtered_points = filtered_points.T

    return filtered_points

def main():
    data_dir='D:/BaiduNetdiskDownload/modelnet40_normal_resampled/modelnet40_normal_resampled/'

    point_cloud_pynt = PyntCloud.from_file(
            data_dir+'cup/cup_0002.txt', sep=",",
            names=["x", "y", "z", "nx", "ny", "nz"])

    point = point_cloud_pynt.points


    # 转成open3d能识别的格式
    # point_cloud_o3d = point_cloud_pynt.to_instance("open3d", mesh=False)
    # o3d.visualization.draw_geometries([point_cloud_o3d]) # 显示原始点云

    # 调用voxel滤波函数，实现滤波
    filtered_cloud = voxel_filter(point, 0.05)  #100

    point_cloud_pynt.points = filtered_cloud
    point_cloud_o3d = point_cloud_pynt.to_instance("open3d", mesh=False)
    # point_cloud_o3d.points = o3d.utility.Vector3dVector(for_show)
    # 显示滤波后的点云
    o3d.visualization.draw_geometries([point_cloud_o3d])


if __name__ == '__main__':
    main()

两种方式的对比结果如下：

在这里插入图片描述
特别注意: 两种方法的代码其函数voxel_filter()的输入输出有所不同，因为随机滤波只需要x,y,z三个特征，所以输入只需要前三个特征，输出为下采样后样本的索引值。中心滤波需要求全部特征的平均值，所以对6个特征求平均值，并存储至DataFrame中进行输出，返回结果就是下采样后的点云。