当前位置： 首页 > article >正文

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

article 2024/11/20 7:44:18

PCA基础

步骤
  1. 计算样本的各变量的均值，然后将样本减去均值进行去中心化，相当于将原点移动到样本的中心。
  2. 计算去中心化后的协方差矩阵（去中心化后很容易计算，把样本矩阵（m行n列代表m个样本n维变量）的转置和样本矩阵自身相乘就可以得到）
  3. 计算协方差矩阵的特征值和特征向量（python和matlab都用现成的库可以计算）
  4. 最后将数据投影到新的坐标轴，也就是特征向量上，具体实现是坐标与特征向量做内积。如果是降维的话就根据特征值大小来选择前百分之多少的主成分，然后将数据投影到选出来的主成分上达到降维的目的。（根据线性代数和统计学理论，协方差矩阵的特征值一定是非负的）
特征值与特征向量的含义
  • 协方差矩阵的特征值，据StatQuest中说是数据投影到某个特征向量轴上的投影点距原点距离的平方。其实反映的是该特征值对应的某一主成分的贡献，特征值越大，则贡献越高，具体还可以通过计算某一特征值与所有特征值之和的比值来量化这个贡献的多少。
  • 协方差矩阵的特征向量其实是一个混合了多种变量的主成分，假如一个特征向量是 [ 4 , 1 , 2 ] [4, 1, 2] [4,1,2]，那么代表新的一个主成分有4份的变量1、1份的变量2以及2份的变量3组成。注意对于实对称矩阵而言，特征值不同的特征向量是两两正交的。
  • 实际上把原始数据与特征向量矩阵做乘积就达到了原始数据投影在特征向量所组成的新坐标轴的目的。

点云（网格）处理的PCA

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


http://www.kler.cn/a/401223.html

相关文章：

  • 【第七节】在RadAsm中使用OllyDBG调试器
  • 任务管理功能拆解——如何高效管理项目任务？
  • Android集成FCM（Firebace Cloud Messaging ）
  • AWD脚本编写_1
  • android 如何获取当前 Activity 的类名和包名
  • 【HCIP]——OSPF综合实验
  • 「三」体验HarmonyOS端云一体化开发模板——使用DevEco Studio直接创建端云一体化工程
  • 计算机视觉：赋予机器“看”的能力
  • element表单校验
  • 高级java每日一道面试题-2024年11月12日-框架篇[SpringBoot篇]-SpringBoot中的监视器是什么?
  • 【mongo8社区办】mongosh MongoServerSelectionError 超时
  • uniapp实现中英文切换
  • 15分钟学 Go 第 59 天 ：更高级的Go话题——接触微服务
  • 力扣题目总结
  • 解析 Apache SeaTunnel 的任务运行过程
  • 第7章硬件测试-7.4 专业实验
  • 内容分发网络CDN、动态内容缓存简介
  • Vite 基础理解及应用
  • ThreadLocal 和 Caffeine 缓存是两种不同的缓存机制，它们在用途和实现上有明显的区别
  • hash表和B树
  • 《人工智能深度学习的基本路线图》
  • DevOps-Jenkins-新手入门级
  • IndentationError: unindent does not match any outer indentation level
  • 汽车资讯新高度：Spring Boot技术飞跃
  • 百度智能云 VectorDB 优势数量 TOP 1
  • 湘潭大学软件工程算法设计与分析考试复习笔记（三）