基于最小二乘发的平面拟合

function [abcd]=PlaneFitting（s） %ABCD为平面方程系数，s为测量点阵求平面方程平面方程D=Ax+By+Cz
x=s（:1）;%s的第1列x
y=s（:2）;%s的第2列
z=s（:3）;%s的第3列
so=mean（s1）;%对xyz求均值
xo=so（:1）;%x的均值
yo=so（:2）;%y的均值
zo=so（:3）;%z的均值
os=bsxfun（@minussso）;%xyz与其均值之差
ox=os（:1）;%x与其均值之差
oy=os（:2）;%y与其均值之差
oz=os（:3）;%z与其均值之差
oxox=sum（bsxfun（@timesoxox））;
oyoy=sum（bsxfun（@timesoyoy））;
ozoz=sum（bsxfun（@timesozoz））;%xyz的方差
oxoy=sum（bsxfun（@timesoxoy））;
oxoz=sum（bsxfun（@timesoxoz））;
oyoz=sum（bsxfun（@timesoyoz））;%xyz的互协方差
oxyz=[oxoxoxoyoxoz;oxoyoyoyoyoz;oxozoyozozoz];%xyz的协方差阵
% [VD]=eig（oxyz）%进行特征分解
% [~t]=min（[abs（[D（11）D（22）D（33）]）]）%最小特征值所对应的特征向量就是平面方程系数
% a=V（1t）;
% b=V（2t）;
% c=V（3t）;
[USV]=svd（oxyz）;% 协方差矩阵的SVD变换中，最小奇异值对应的奇异向量就是平面的方向
a=V（13）;
b=V（23）;
c=V（33）;
d=a*xo+b*yo+c*zo;%输出平面方程系数

% % % n=size（s1）;
% % % xx=x.*x;
% % % yy=y.*y;
% % % zz=z.*z;
% % % xy=x.*y;
% % % xz=x