点云（网格）PCA及其存在的问题

🕗 发布于 2024-11-20 11:08 计算机视觉 3d

PCA基础

步骤

计算样本的各变量的均值，然后将样本减去均值进行去中心化，相当于将原点移动到样本的中心。
计算去中心化后的协方差矩阵（去中心化后很容易计算，把样本矩阵（m行n列代表m个样本n维变量）的转置和样本矩阵自身相乘就可以得到）
计算协方差矩阵的特征值和特征向量（python和matlab都用现成的库可以计算）
最后将数据投影到新的坐标轴，也就是特征向量上，具体实现是坐标与特征向量做内积。如果是降维的话就根据特征值大小来选择前百分之多少的主成分，然后将数据投影到选出来的主成分上达到降维的目的。（根据线性代数和统计学理论，协方差矩阵的特征值一定是非负的）

特征值与特征向量的含义

协方差矩阵的特征值，据StatQuest中说是数据投影到某个特征向量轴上的投影点距原点距离的平方。其实反映的是该特征值对应的某一主成分的贡献，特征值越大，则贡献越高，具体还可以通过计算某一特征值与所有特征值之和的比值来量化这个贡献的多少。
协方差矩阵的特征向量其实是一个混合了多种变量的主成分，假如一个特征向量是 $[4, 1, 2]$ ，那么代表新的一个主成分有4份的变量1、1份的变量2以及2份的变量3组成。注意对于实对称矩阵而言，特征值不同的特征向量是两两正交的。
实际上把原始数据与特征向量矩阵做乘积就达到了原始数据投影在特征向量所组成的新坐标轴的目的。

点云（网格）处理的PCA

步骤

通过体积矩计算点云的质心，具体计算可参考[4]和[5]，然后将点云坐标减去质心坐标进行去中心化。
对去中心化的点云计算二阶矩构成协方差矩阵（具体计算参考论文中）。
对二阶矩构成的协方差矩阵进行特征向量分解得到特征值和特征向量（按照特征值从大到小进行排序）。
将原始点云数据与得到的特征向量矩阵进行相乘得到旋转的点云（可通过制定特征向量列的顺序来制定某个主轴，也就是主方向是和X轴、Y轴还是Z轴对齐）

实现代码如下：

import open3d as o3d
import numpy as np
from Plyobject import *

# o3d读取
model_name = "bunny.ply"
o3d_model = o3d.io.read_triangle_mesh("./"+model_name)
model = Plyobject(np.asarray(o3d_model.vertices), np.asarray(o3d_model.triangles))

# 计算质心去中心化
center = model.get_center_by_vm()
print("center of mass:", center)

center_model = Plyobject(model.pts-center, model.face)

# 计算二阶矩及特征向量（算的是去中心化后的矩）
M1 = center_model.get_1ord_moments()
M2 = center_model.get_2ord_moments()
# print("归一化前的矩：", M1, '\n', M2)
eigenvalues, eigenvectors = np.linalg.eig(M2) #每一列是一个特征向量

print(eigenvalues, "\n", eigenvectors)

# 进行PCA旋转
eigenvectors = eigenvectors[:,[2,1,0]] # 最大主轴和Z轴对齐，最小轴和X轴对齐
pca_points = np.matmul(center_model.pts, eigenvectors)

# 测试一下归一化后的一阶矩和部分二阶矩是否为零
pca_model = Plyobject(pca_points, center_model.face)
pca_M1 = pca_model.get_1ord_moments()
pca_M2 = pca_model.get_2ord_moments()
pca_M3 = pca_model.get_3ord_moments()
print("归一化后的矩：", pca_M1, '\n', pca_M2, '\n', pca_M3)


# o3d写出
pca_mesh = o3d.geometry.TriangleMesh()
pca_mesh.vertices = o3d.utility.Vector3dVector(pca_points)
pca_mesh.triangles = o3d.utility.Vector3iVector(center_model.face)
o3d.io.write_triangle_mesh('PCA_'+ model_name, pca_mesh, write_ascii=True)

Plyobject是自己定义的一个包，有质心和不同阶矩的计算方法：

class Plyobject:
    # 构造方法
    def __init__(self, pts, face):
        self.pts = pts
        self.face = face
        self.area = self.get_area()

    # 计算三角形面的面积
    def get_face_area(self, v1, v2, v3):
        v1v2 = v2 - v1
        v1v3 = v3 - v1
        vc = np.cross(v1v2, v1v3)
        area = np.linalg.norm(vc) * 0.5
        return area

    # 计算物体表面积
    def get_area(self):
        sum_area = 0
        for f in self.face:
            v1 = self.pts[f[0]]
            v2 = self.pts[f[1]]
            v3 = self.pts[f[2]]
            area = self.get_face_area(v1, v2, v3)
            sum_area += area
        self.area = sum_area
        return sum_area

    # 计算面积矩的物体中心
    def get_center_by_area(self):
        sum = 0
        for f in self.face:
            v1 = self.pts[f[0]]
            v2 = self.pts[f[1]]
            v3 = self.pts[f[2]]
            area = self.get_face_area(v1, v2, v3)
            sum += area * 1 / 3 * (v1 + v2 + v3)
        center = sum / self.area
        self.center = center
        return center

    def get_center_by_avg(self):
        xb = np.mean(self.pts[:, 0])
        yb = np.mean(self.pts[:, 1])
        zb = np.mean(self.pts[:, 2])
        return [xb, yb, zb]
    def get_center_by_vm(self):
        M100, M010, M001, M000 = 0, 0, 0, 0
        for i in self.face:
            index_v1 = i[0]
            index_v2 = i[1]
            index_v3 = i[2]
            t_M000 = 1 / 6 * (- self.pts[index_v3][0] * self.pts[index_v2][1] * self.pts[index_v1][2]
                              + self.pts[index_v2][0] * self.pts[index_v3][1] * self.pts[index_v1][2]
                              + self.pts[index_v3][0] * self.pts[index_v1][1] * self.pts[index_v2][2]
                              - self.pts[index_v1][0] * self.pts[index_v3][1] * self.pts[index_v2][2]
                              - self.pts[index_v2][0] * self.pts[index_v1][1] * self.pts[index_v3][2]
                              + self.pts[index_v1][0] * self.pts[index_v2][1] * self.pts[index_v3][2])
            M100 += 1 / 4 * (self.pts[index_v1][0] + self.pts[index_v2][0] + self.pts[index_v3][0]) * t_M000
            M010 += 1 / 4 * (self.pts[index_v1][1] + self.pts[index_v2][1] + self.pts[index_v3][1]) * t_M000
            M001 += 1 / 4 * (self.pts[index_v1][2] + self.pts[index_v2][2] + self.pts[index_v3][2]) * t_M000
            M000 += t_M000

        xb = M100 / M000
        yb = M010 / M000
        zb = M001 / M000
        return [xb, yb, zb]
    
    def get_3ord_moments(self):
        M300, M030, M003, M000 = 0,0,0,0
        for i in self.face:
            index_v1 = i[0]
            index_v2 = i[1]
            index_v3 = i[2]
            t_M000 = 1 / 6 * (- self.pts[index_v3][0] * self.pts[index_v2][1] * self.pts[index_v1][2]
                            + self.pts[index_v2][0] * self.pts[index_v3][1] * self.pts[index_v1][2]
                            + self.pts[index_v3][0] * self.pts[index_v1][1] * self.pts[index_v2][2]
                            - self.pts[index_v1][0] * self.pts[index_v3][1] * self.pts[index_v2][2]
                            - self.pts[index_v2][0] * self.pts[index_v1][1] * self.pts[index_v3][2]
                            + self.pts[index_v1][0] * self.pts[index_v2][1] * self.pts[index_v3][2])
            M000 += t_M000
            M300 = 1 / 20 *(self.pts[index_v1][0] ** 3 + self.pts[index_v2][0] ** 3 + self.pts[index_v3][0] ** 3 +
                            self.pts[index_v1][0] ** 2 * (self.pts[index_v2][0] + self.pts[index_v3][0])+
                            self.pts[index_v2][0] ** 2 * (self.pts[index_v1][0] + self.pts[index_v3][0])+
                            self.pts[index_v3][0] ** 2 * (self.pts[index_v1][0] + self.pts[index_v2][0])+
                            self.pts[index_v1][0] * self.pts[index_v2][0] * self.pts[index_v3][0]) * M000
            M030 = 1 / 20 *(self.pts[index_v1][1] ** 3 + self.pts[index_v2][1] ** 3 + self.pts[index_v3][1] ** 3 +
                            self.pts[index_v1][1] ** 2 * (self.pts[index_v2][1] + self.pts[index_v3][1])+
                            self.pts[index_v2][1] ** 2 * (self.pts[index_v1][1] + self.pts[index_v3][1])+
                            self.pts[index_v3][1] ** 2 * (self.pts[index_v1][1] + self.pts[index_v2][1])+
                            self.pts[index_v1][1] * self.pts[index_v2][2] * self.pts[index_v3][2]) * M000
            M003 = 1 / 20 *(self.pts[index_v1][2] ** 3 + self.pts[index_v2][2] ** 3 + self.pts[index_v3][2] ** 3 +
                            self.pts[index_v1][2] ** 2 * (self.pts[index_v2][2] + self.pts[index_v3][2])+
                            self.pts[index_v2][2] ** 2 * (self.pts[index_v1][2] + self.pts[index_v3][2])+
                            self.pts[index_v3][2] ** 2 * (self.pts[index_v1][2] + self.pts[index_v2][2])+
                            self.pts[index_v1][2] * self.pts[index_v2][2] * self.pts[index_v3][2]) * M000
        return [M300, M030, M003]

    def get_2ord_moments(self):
        M200,M020,M002,M110,M101,M011,M100,M010,M001,M000 = 0,0,0,0,0,0,0,0,0,0
        for i in self.face:
            index_v1 = i[0]
            index_v2 = i[1]
            index_v3 = i[2]
            t_M000 = 1 / 6 * (- self.pts[index_v3][0] * self.pts[index_v2][1] * self.pts[index_v1][2]
                            + self.pts[index_v2][0] * self.pts[index_v3][1] * self.pts[index_v1][2]
                            + self.pts[index_v3][0] * self.pts[index_v1][1] * self.pts[index_v2][2]
                            - self.pts[index_v1][0] * self.pts[index_v3][1] * self.pts[index_v2][2]
                            - self.pts[index_v2][0] * self.pts[index_v1][1] * self.pts[index_v3][2]
                            + self.pts[index_v1][0] * self.pts[index_v2][1] * self.pts[index_v3][2])
            M100 += 1 / 4 * (self.pts[index_v1][0] + self.pts[index_v2][0] + self.pts[index_v3][0]) * t_M000
            M010 += 1 / 4 * (self.pts[index_v1][1] + self.pts[index_v2][1] + self.pts[index_v3][1]) * t_M000
            M001 += 1 / 4 * (self.pts[index_v1][2] + self.pts[index_v2][2] + self.pts[index_v3][2]) * t_M000
            M000 += t_M000



            M200 += 1 / 10 * (
                self.pts[index_v1][0] ** 2 + self.pts[index_v2][0] ** 2 + self.pts[index_v3][0] ** 2
                + self.pts[index_v1][0] * self.pts[index_v2][0]
                + self.pts[index_v2][0] * self.pts[index_v3][0]
                + self.pts[index_v1][0] * self.pts[index_v3][0]) * t_M000
            M020 += 1 / 10 * (
                    self.pts[index_v1][1] ** 2 + self.pts[index_v2][1] ** 2 + self.pts[index_v3][1] ** 2
                    + self.pts[index_v1][1] * self.pts[index_v2][1]
                    + self.pts[index_v2][1] * self.pts[index_v3][1]
                    + self.pts[index_v1][1] * self.pts[index_v3][1]) * t_M000
            M002 += 1 / 10 * (
                    self.pts[index_v1][2] ** 2 + self.pts[index_v2][2] ** 2 + self.pts[index_v3][2] ** 2
                    + self.pts[index_v1][2] * self.pts[index_v2][2]
                    + self.pts[index_v2][2] * self.pts[index_v3][2]
                    + self.pts[index_v1][2] * self.pts[index_v3][2]) * t_M000
            M110 += 1 / 20 * (2 * (self.pts[index_v1][0] * self.pts[index_v1][1] +
                                self.pts[index_v2][0] * self.pts[index_v2][1] +
                                self.pts[index_v3][0] * self.pts[index_v3][1]) +
                            self.pts[index_v1][0] * self.pts[index_v2][1] +
                            self.pts[index_v2][0] * self.pts[index_v1][1] +
                            self.pts[index_v1][0] * self.pts[index_v3][1] +
                            self.pts[index_v3][0] * self.pts[index_v1][1] +
                            self.pts[index_v2][0] * self.pts[index_v3][1] +
                            self.pts[index_v3][0] * self.pts[index_v2][1]
                            ) * t_M000

            M101 += 1 / 20 * (2 * (self.pts[index_v1][0] * self.pts[index_v1][2] +
                                self.pts[index_v2][0] * self.pts[index_v2][2] +
                                self.pts[index_v3][0] * self.pts[index_v3][2]) +
                            self.pts[index_v1][0] * self.pts[index_v2][2] +
                            self.pts[index_v2][0] * self.pts[index_v1][2] +
                            self.pts[index_v1][0] * self.pts[index_v3][2] +
                            self.pts[index_v3][0] * self.pts[index_v1][2] +
                            self.pts[index_v2][0] * self.pts[index_v3][2] +
                            self.pts[index_v3][0] * self.pts[index_v2][2]
                            ) * t_M000

            M011 += 1 / 20 * (2 * (self.pts[index_v1][1] * self.pts[index_v1][2] +
                                self.pts[index_v2][1] * self.pts[index_v2][2] +
                                self.pts[index_v3][1] * self.pts[index_v3][2]) +
                            self.pts[index_v1][1] * self.pts[index_v2][2] +
                            self.pts[index_v2][1] * self.pts[index_v1][2] +
                            self.pts[index_v1][1] * self.pts[index_v3][2] +
                            self.pts[index_v3][1] * self.pts[index_v1][2] +
                            self.pts[index_v2][1] * self.pts[index_v3][2] +
                            self.pts[index_v3][1] * self.pts[index_v2][2]
                            ) * t_M000
            M2 = np.array([[M200, M110, M101],
                            [M110, M020, M011],
                            [M101, M011, M002]])
            return M2
    
    def get_1ord_moments(self):
        M100,M010,M001,M000 = 0,0,0,0
        for i in self.face:
            index_v1 = i[0]
            index_v2 = i[1]
            index_v3 = i[2]
            t_M000 = 1 / 6 * (- self.pts[index_v3][0] * self.pts[index_v2][1] * self.pts[index_v1][2]
                            + self.pts[index_v2][0] * self.pts[index_v3][1] * self.pts[index_v1][2]
                            + self.pts[index_v3][0] * self.pts[index_v1][1] * self.pts[index_v2][2]
                            - self.pts[index_v1][0] * self.pts[index_v3][1] * self.pts[index_v2][2]
                            - self.pts[index_v2][0] * self.pts[index_v1][1] * self.pts[index_v3][2]
                            + self.pts[index_v1][0] * self.pts[index_v2][1] * self.pts[index_v3][2])
            M100 += 1 / 4 * (self.pts[index_v1][0] + self.pts[index_v2][0] + self.pts[index_v3][0]) * t_M000
            M010 += 1 / 4 * (self.pts[index_v1][1] + self.pts[index_v2][1] + self.pts[index_v3][1]) * t_M000
            M001 += 1 / 4 * (self.pts[index_v1][2] + self.pts[index_v2][2] + self.pts[index_v3][2]) * t_M000
            M000 += t_M000
        return [M100, M010, M001]

PCA点云配准存在的问题

当模型旋转角度过大时或者模型噪声较多或采样不规则时，PCA会存在配准结果不一的情况，具体可参加参考[3]论文中的分析。

参考文献/链接

[1] https://www.bilibili.com/video/BV1X54y1R7g7/?spm_id_from=333.337.search-card.all.click&vd_source=bac8ddf04ec0b6386d58110f67353bc7

[2] https://www.bilibili.com/video/BV1AZ4y1w7Y3/?spm_id_from=333.337.search-card.all.click&vd_source=bac8ddf04ec0b6386d58110f67353bc7

[3] 2006_What is wrong with mesh PCA in coordinate direction normalization

[4] 2001_Efficient feature extraction for 2d3d objects IN MESH REPRESENTATION

[5] 2011_Robust and blind mesh watermarking based on volume moments

原文地址：https://blog.csdn.net/qq_45607390/article/details/143870161

免责声明：本站文章内容转载自网络资源，如本站内容侵犯了原著者的合法权益，可联系本站删除。更多内容请关注自学内容网（zxcms.com）！

上一篇：数据分析-51-时间序列分解之局部均值分解LMD
下一篇：torch.__version__的torch版本和conda list的torch版本不一致

STM32:Bootloader(AB备份，自动回滚)
STM32的OTA升级，Bootloader项目级代码
阅读更多2024-12-04
微信小程序px和rpx单位互转方法
【代码】微信小程序px和rpx单位互转方法。
阅读更多2024-12-04
vue-注册全局组件并使用
vue全局组件如何通过方法调用
阅读更多2024-12-04
反向代理模块
反向代理是指以代理服务器来接收客户端的请求，然后将请求转发给内部网络上的服务器，将从服务器上得到的结果返回给客户端，此时代理服务器对外表现为一个反向代理服务器。不适用代理的情况下，外网不能访问企业内网
阅读更多2024-12-04
【k8s 深入学习之 event 聚合】event count累记聚合（采用 Patch），Message 聚合形成聚合 event（采用Create）
【代码】【k8s 深入学习之 event 聚合】event count累记聚合（采用 Patch），Message 聚合形成聚合 event（采用Create）
阅读更多2024-12-04
adb导出系统apk
3.导出apk到电脑使用adb pull命令将APK文件导出到本地计算机： adb pull <路径>==<apk名称>2.查找路径 adb shell pm path &
阅读更多2024-12-04
【adb】AndroidStudio调试
1：配置系统adb环境变量，找到 Android SDK下的 platform-tools目录，复制路径，添加到系统环境变量中的Path中。3.在adb目录，运行cmd，输入命令 adb connec
阅读更多2024-12-04
python笔记3
简单（直接）输出print （）输出到指定文件fp.close()调用日期模块输出import datetime # 调用日期模块datetimeprint('现在的年份是' + str(date.y
阅读更多2024-12-04
什么语言适合做 Serverless 开发？
随着云计算的普及，**无服务器架构（Serverless Architecture）**成为一种流行的开发模式，它使得开发者无需管理服务器基础设施，专注于编写应用逻辑。无服务器架构通常按需提供计算资源
阅读更多2024-12-04
Qt Designer Ui设计
【代码】Qt Designer Ui设计。
阅读更多2024-12-04