function?[abcd]=PlaneFitting（s）?%ABCD為平面方程系數(shù)，s為測量點(diǎn)陣求平面方程平面方程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的互協(xié)方差
oxyz=[oxoxoxoyoxoz;oxoyoyoyoyoz;oxozoyozozoz];%xyz的協(xié)方差陣
%?[VD]=eig（oxyz）%進(jìn)行特征分解
%?[~t]=min（[abs（[D（11）D（22）D（33）]）]）%最小特征值所對應(yīng)的特征向量就是平面方程系數(shù)
%?a=V（1t）;
%?b=V（2t）;
%?c=V（3t）;
[USV]=svd（oxyz）;%?協(xié)方差矩陣的SVD變換中，最小奇異值對應(yīng)的奇異向量就是平面的方向
a=V（13）;
b=V（23）;
c=V（33）;
d=a*xo+b*yo+c*zo;%輸出平面方程系數(shù)

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