1 概述
利用主成分分析(Principal Components Analysis, PCA)方法，可计算待拟合点的法向量，进而得到平面参数。
原理详见参考文献：
Pauly M, Keiser R, Gross M. Multi‐scale feature extraction on point‐sampled surfaces[C]//Computer graphics forum. Oxford, UK: Blackwell Publishing, Inc, 2003, 22(3): 281-289.
2 代码实现
function planes = fitPlane_PCA(data)
% 功能：利用PCA拟合平面
% 输入：data - 原始数据(m*3)
% 输出：planes - 拟合所得平面参数
points = data(:,1:3);
[m,~] = size(points);
% 计算协方差矩阵
M = points-ones(m,1)*(sum(points,1)/m);
C = M.'*M./(m-1);
% 计算特征值与特征向量
[V, D] = eig(C);
% 最小特征值对应的特征向量为法向量
s1 = D(1,1);
s2 = D(2,2);
s3 = D(3,3);
if (s1 <= s2 && s1 <= s3)
normal(1,:) = V(:,1)/norm(V(:,1));
elseif (s2 <= s1 && s2 <= s3)
normal(1,:) = V(:,2)/norm(V(:,2));
elseif (s3 <= s1 && s3 <= s2)
normal(1,:) = V(:,3)/norm(V(:,3));
end
% 平面参数标准化
dtmp = mean(points*normal');
planes(1:3) = normal'*sign(dtmp);
planes(4) = -dtmp*sign(dtmp);
end
3 可视化验证
为了检测平面拟合的效果，采用仿真数据进行验证：
%% 数据准备
% 读取数据
data = load('data.txt');
%% PCA平面拟合
% 平面:Ax+By+Cz+D=0
planes = fitPlane_PCA(data);
A = planes(1);
B = planes(2);
C = planes(3);
D = planes(4);
%% 可视化验证
% 窗口尺寸设置（单位：厘米）
figureUnits = 'centimeters';
figureWidth = 15;
figureHeight = 15;
figureHandle = figure;
set(gcf, 'Units', figureUnits, 'Position', [0 0 figureWidth figureHeight]);
% 原始点云可视化
l = scatter3(data(:,1),data(:,2),data(:,3),15);% 原始点云
hold on
% 拟合平面绘制
xfit = min(data(:,1)):0.01:max(data(:,1));
yfit = min(data(:,2)):0.01:max(data(:,2));
[XFit,YFit]= meshgrid (xfit,yfit);
ZFit = -(D + A * XFit + B * YFit)/C;
s = surf(XFit,YFit,ZFit,'EdgeColor','none');
hTitle = title('基于PCA的平面拟合');
hXLabel = xlabel('XAxis');
hYLabel = ylabel('YAxis');
hZLabel = zlabel('ZAxis');
% 细节优化
colormap(map)
set(l,'Marker','o','MarkerFaceColor',CR,'MarkerEdgeColor',CR)
view(-27.5,46.9);% 视角
set(gca, 'Box', 'on', ... % 边框
'XGrid','on','YGrid','on','ZGrid','on', ... % 网格
'TickDir', 'out', 'TickLength', [0.01 0.01], ... % 刻度
'XMinorTick', 'off', 'YMinorTick', 'off', ... % 小刻度
'XColor', [.1 .1 .1], 'YColor', [.1 .1 .1],... % 坐标轴颜色
'XLim',[-1.4 0.2])
% 字体和字号
set(gca, 'FontName', 'Arial', 'FontSize', 10)
set([hXLabel, hYLabel, hZLabel], 'FontName', 'Arial','FontSize', 11)
set(hTitle, 'FontName', '微软雅黑', 'FontSize', 12, 'FontWeight' , 'bold')
% 背景颜色
set(gcf,'Color',[1 1 1])
%
输出
print('test.png','-r300','-dpng')
其中，为了区分不同对象，从Matlab配色神器TheColor的XKCD和SCI颜色库中选择对比色及渐变色：
% 颜色定义
CR = TheColor('xkcd',154);
map = TheColor('sci',2068);
(点击上图查看TheColor具体功能)
获取方式：公众号（阿昆的科研日常）后台回复 TC
最终结果如下：
以上。
下载方式
原创不易，请按照以下方式获取：
Step1：点击文章最下方的在看按钮
Step2：在本公众号（阿昆的科研日常）后台回复关键字“点云2”下载
转自：“阿昆的科研日常”微信公众号
如有侵权，请联系本站删除！
