PCA主成分分析法在点云曲面拟合及曲面法向量计算中的应用
PCA主成分分析法在点云曲面拟合及曲面法向量计算中的应用
PCA（主成分分析）是一种数学变换方法，利用降维的思想在变换中保持变量的总方差不变，将给定的一组变量线性变换为另一组不相关的变量，并且使变换后的第一变量的方差最大，即第一主成分，其他分量的方差依次递减。以二维数据为例：
图1 二维数据PCA分析
在这个图中，原始数据是二维的，所以PCA的目标是找到两个新的坐标轴（称为主成分），使得这些新坐标轴上的数据方差最大化。图中的第一条绿色直线（y’轴）表示第一个主成分，它与原始x轴有一个角度。这个角度表示第一个主成分的方向，它捕捉了数据的最大方差。第二条绿色直线（x’轴）垂直于第一条直线，表示第二个主成分，它捕捉了剩余的方差。蓝色的点在原始坐标系（x, y）中分布不均匀，但在新的坐标系（x’, y’）中看起来更加均匀。这说明PCA成功地将数据转换到了一个新的空间，使得数据在新的坐标轴上的分布更加均匀，更容易进行后续的分析或可视化。
在点云数据中，变量为三维点坐标的集合，其变量为X、Y、Z三个坐标值，则经过变换后，应有三个主成分。对于一个空间平面，在平行于平面的方向上点集分布最为离散，方差最大，在垂直于平面的方向上，点集分布最为集中，方差最小，即空间平面的第三主成分为垂直于空间平面的向量即曲面的法向量。由于平面拟合最关键的为法向量的拟合，利用PCA得到点集的第三主成分，即能进一步拟合出平面方程，如图2所示。
图2 点云图
拟合得到的曲面如图3所示：
图3 三次样条插值法拟合曲面
% 指定球面的参数
R = 500; % 半径
theta0=acos((7.5/2)/R);
phi0=asin((6.5/2)/R);
theta_range = [theta0, pi-theta0]; % 极角范围
phi_range = [0, phi0]; % 方位角范围
% 生成球面的参数
N = 20; % 离散点的数量
[theta, phi] = meshgrid(linspace(theta_range(1), theta_range(2), N), linspace(phi_range(1), phi_range(2), N));
% 计算球面上的点的坐标
x = R * sin(theta) .* cos(phi);
y = R * sin(theta) .* sin(phi);
z = R * cos(theta);
% 添加噪声，模拟测量不准
max_noise_amplitude = 0.001;
noise = max_noise_amplitude * randn(size(x));
x_with_noise = x + noise;
% 绘制球面部分
surf(x, y, z);
axis equal;
xlabel('X');
ylabel('Y');
zlabel('Z');
title('截取球面部分');
% 合并为二维数组
data = [reshape(z, [], 1), reshape(y, [], 1), reshape(x_with_noise, [], 1)];
%模拟吊装不平，使之旋转一定角度
xita = 10 * pi / 180; % 定义旋转角度
% 计算绕 x 轴旋转的旋转矩阵
Rx = [1, 0, 0;
0, cos(xita), -sin(xita);
0, sin(xita), cos(xita)];
data = (Rx * data')'; % 将数据点乘以旋转矩阵
final_data = rotateplane_PCA( data );%获取垂直于z轴的曲面
% 对点云进行密集三次样条插值，构成曲面
% 为了进行三次样条插值，我们需要一个规则的网格
[xq, yq] = meshgrid(linspace(min(final_data(:,1)), max(final_data(:,1)), 50), linspace(min(final_data(:,2)), max(final_data(:,2)), 50));
% 使用scatteredInterpolant进行插值
F = scatteredInterpolant(final_data(:,1), final_data(:,2), final_data(:,3), 'natural', 'linear');
zq = F(xq, yq);
% 绘制插值后的点云
figure;
surf(xq, yq, zq, 'EdgeColor', 'none'); % 绘制曲面
hold on;
% scatter3(final_data(:,1), final_data(:,2), final_data(:,3), 50, final_data(:,3), 'filled'); % 绘制原始点云
colorbar; % 添加颜色条
title('Interpolated Point Cloud with Color Mapping by Z-axis');
xlabel('X'); ylabel('Y'); zlabel('Z');
axis equal;
view(3);
% 映射颜色到z轴
colormap jet; % 使用 jet 颜色图
caxis([min(final_data(:,3)), max(final_data(:,3))]); % 设置颜色轴范围
hold off;
function [ final_data ] = rotateplane_PCA( data )
%UNTITLED4 此处显示有关此函数的摘要
% 该函数用PCA主成分分析法生成旋转矩阵，使曲面平铺在yoz平面上
%输入为三维散点数据，输出为旋转后的散点数据
% 将数据中心化（去均值）
mean_points = mean(data);
centered_points = data - mean_points;
%计算协方差矩阵
cov_matrix = cov(centered_points);
% 进行特征值分解
[V, D] = eig(cov_matrix);
%按特征值从大到小排序
[~, order] = sort(diag(D), 'descend');
V = V(:, order);
%旋转点云数据，使其在新坐标系下的z轴方向上的点的方差最小
final_data = centered_points * V;
%绘制旋转前后的点云
figure;
subplot(1,2,1);
scatter3(data(:,1), data(:,2), data(:,3), 'filled','r');
title('Original Point Cloud');
xlabel('X'); ylabel('Y'); zlabel('Z');
axis equal;
subplot(1,2,2);
scatter3(final_data(:,1), final_data(:,2), final_data(:,3), 'filled');
title('Rotated Point Cloud');
xlabel('X'); ylabel('Y'); zlabel('Z');
axis equal;
% % 7. 旋转矩阵
% rotation_matrix = V;
% disp('Rotation Matrix:');
% disp(rotation_matrix);
end
其实从本质上来说，这里所用到的插值法拟合曲面并非是真正的拟合曲面，而是在原始点云中插入了许多点，很密集看起来像是一个曲面了。还有就是代码中并没有写关于法向量的求取，但是可以根据旋转矩阵V计算出来。例如原来曲面的法向量是L,旋转之后的曲面平行于XOY，所以新曲面的法向量就是L1=（0,0,1）.因此LV=L1,所以L=L1inv(V)。
本文如有不当之处，欢迎各位批评指正。
本文参考了两篇文章：
https://blog.csdn.net/weixin_45710335/article/details/136559483
https://zhuanlan.zhihu.com/p/668272208