1. 程式人生 > >根據空間點座標擬合平面和直線(線性迴歸和svd分解)

根據空間點座標擬合平面和直線(線性迴歸和svd分解)

根據一組點的座標擬合空間平面,有兩種方法
第一種:如果在測量得到的資料中,x,y值都是確認沒有誤差的,而誤差只是出現在z值上,則可以使用線性迴歸的方法,此方法最小二乘的目標是在z方向上de殘差
Matlab 程式碼
  1. % 隨機生成一組(x,y,z),這些點的座標離一個空間平面比較近
  2. x0=1,L1=2;
  3. y0=1,L2=2;
  4. x=x0+rand(20,1)*L1;
  5. y=y0+rand(20,1)*L2;
  6. z=1+2*x+3*y;
  7. scatter3(x,y,z,'filled')
  8. hold on;
  9. X = [ones(length(x),1) x y];
  10. % 擬合,其實是線性迴歸,但可以用來擬合平面
  11. % 輸出為 b = [b(1) b(2) b(3)] 表示 z = b(1) + b(2)*x + b(3)*y 是擬合出來的平面的方程
  12. [b,bint,r,rint,stats] = regress(z,X,95);
  13. % 圖形繪製
  14. xfit = min(x):0.1:max(x);
  15. yfit = min(y):0.1:max(y);
  16. [XFIT,YFIT]= meshgrid (xfit,yfit);
  17. ZFIT = b(1) + b(2) * XFIT + b(3) * YFIT;
  18. mesh(XFIT,YFIT,ZFIT);
複製程式碼輸出結果如圖1

第二中: 如果在測量得到的資料中,x,y,z都存在誤差,則最小化的目標應該是測量點到平面距離的殘差
Matlab程式碼
  1. % 隨機生成一組(x,y,z),這些點的座標離一個空間平面比較近
  2. x0=1,L1=2;
  3. y0=1,L2=2;
  4. x=x0+rand(20,1)*L1;
  5. y=y0+rand(20,1)*L2;
  6. z=1+2*x+3*y;
  7. scatter3(x,y,z,'filled')
  8. hold on;
  9. planeData=[x,y,z];
  10. % 協方差矩陣的SVD變換中,最小奇異值對應的奇異向量就是平面的方向
  11. xyz0=mean(planeData,1);
  12. centeredPlane=bsxfun(@minus,planeData,xyz0);
  13. [U,S,V]=svd(centeredPlane);
  14. a=V(1,3);
  15. b=V(2,3);
  16. c=V(3,3);
  17. d=-dot([a b c],xyz0);
  18. % 圖形繪製
  19. xfit = min(x):0.1:max(x);
  20. yfit = min(y):0.1:max(y);
  21. [XFIT,YFIT]= meshgrid (xfit,yfit);
  22. ZFIT = -(d + a * XFIT + b * YFIT)/c;
  23. mesh(XFIT,YFIT,ZFIT);
複製程式碼輸出結果如圖2

而根據空間點擬合一條空間直線的思路比較直接,就是最小化這些點到直線的距離
  1. % 隨機生成一組點,這寫點距離直線l比較近,l的過點[1,1,1],方向向量為[1,2,3]
  2. lineData=bsxfun(@plus, [1,1,1], (-1:.1:1).'*[1,2,3]);
  3. Noise=rand(size(lineData))*.1;
  4. lineData=lineData+Noise;
  5. scatter3(lineData(:,1), lineData(:,2), lineData(:,3),'filled')
  6. hold on;
  7. % 擬合的直線必過所有座標的算數平均值
  8. xyz0=mean(lineData,1),
  9. % 協方差矩陣奇異變換,與擬合平面不同的是
  10. % 所得直線的方向實際上與最大奇異值對應的奇異向量相同
  11. centeredLine=bsxfun(@minus,lineData,xyz0);
  12. [U,S,V]=svd(centeredLine);
  13. direction=V(:,1);
  14. % 畫圖
  15. t=-8:0.1:8;
  16. xx=xyz0(1)+direction(1)*t;
  17. yy=xyz0(2)+direction(2)*t;
  18. zz=xyz0(3)+direction(3)*t;
  19. plot3(xx,yy,zz)
複製程式碼輸出結果如圖