Plane Fit and Sphere Fit in Matlab
![](https://iztuts.com/wp-content/uploads/2018/12/Fit-mat-cau.jpg)
Ta có dữ liệu là đám mây điểm, muốn xem bố trí của đám mây điểm đó xấp xỉ một mặt phẳng hoặc mặt cầu chúng ta dùng các hàm sau:
PlaneFit
function C = planefit(x,y,z) xx = x(:); yy = y(:); zz = z(:); N = length(xx); O = ones(N,1); C = [xx yy O]\zz;
SphereFit
function [xc, yc, zc, R] = sumith_fit(x,y,z) [N,dum]=size(x); Sx = sum(x); Sy = sum(y); Sz = sum(z); Sxx = sum(x.*x); Syy = sum(y.*y); Szz = sum(z.*z); Sxy = sum(x.*y); Sxz = sum(x.*z); Syz = sum(y.*z); Sxxx = sum(x.*x.*x); Syyy = sum(y.*y.*y); Szzz = sum(z.*z.*z); Sxyy = sum(x.*y.*y); Sxzz = sum(x.*z.*z); Sxxy = sum(x.*x.*y); Sxxz = sum(x.*x.*z); Syyz =sum(y.*y.*z); Syzz = sum(y.*z.*z); A1 = Sxx +Syy +Szz; a = 2*Sx*Sx-2*N*Sxx; b = 2*Sx*Sy-2*N*Sxy; c = 2*Sx*Sz-2*N*Sxz; d = -N*(Sxxx +Sxyy +Sxzz)+A1*Sx; e = b; %2*Sx*Sy-2*N*Sxy; f = 2*Sy*Sy-2*N*Syy; g = 2*Sy*Sz-2*N*Syz; h = -N*(Sxxy +Syyy +Syzz)+A1*Sy; j = c; %2*Sx*Sz-2*N*Sxz; k = g; %2*Sy*Sz-2*N*Syz; l = 2*Sz*Sz-2*N*Szz; m = -N*(Sxxz +Syyz + Szzz)+A1*Sz; delta = a*(f*l - g*k)-e*(b*l-c*k) + j*(b*g-c*f); xc = (d*(f*l-g*k) -h*(b*l-c*k) +m*(b*g-c*f))/delta; yc = (a*(h*l-m*g) -e*(d*l-m*c) +j*(d*g-h*c))/delta; zc = (a*(f*m-h*k) -e*(b*m-d*k) +j*(b*h-d*f))/delta; R = sqrt(xc^2+yc^2+zc^2 + (A1 -2*(xc*Sx+yc*Sy+zc*Sz))/N); end
Trước khi sử dụng phải biến đổi dũ liệu cho phù hợp:
Ví dụ: có ma trận Height (601x601) điểm:
Chuyển thành 3 biến: x,y,z như sau:
[Sy,Sx] = size(Height);
x = 1:Sx;
y = 1:Sy;
[xx, yy] = meshgrid(x,y);
zz = Height;
xxx = xx(:); //Kích thước 361201
yyy=yy(:);//Kích thước 361201
zzz = zz(:);//Kích thước 361201
Áp dụng hàm Fit mạt cầu:
[xc, yc, zc, R] = sumith_fit(xxx,yyy,zzz)
%Chon khu vuc phang de lam khop
plan = TopoReal(1:100,1:300);
%Tim mat phang khop nhat
[Sy,Sx] = size(plan);
x = 1:Sx;
y = 1:Sy;
[xx, yy] = meshgrid(x,y);
zz = plan;
xxx = xx(:);
yyy = yy(:);
zzz = zz(:);
%Áp d?ng hàm Fit mat cau:
C = planefit(xxx,yyy,zzz);
%Dung lai mat phang do bang kich thuoc voi TopoReal
[Sy,Sx] = size(TopoReal);
x = 1:Sx;
y = 1:Sy;
[xx, yy] = meshgrid(x,y);
zz = TopoReal;
xxx = xx(:);
yyy = yy(:);
zzz = zz(:);
pl = C(1)*xx+C(2)*yy+C(3);
%Day la bien dang 3D khong bi nghieng
XoaNghieng = zz-pl;
%Chon CHom Cau
%Chom = XoaNghieng(178:620,883:1280);
%Chom = XoaNghieng(117:720,796:1280);
Chom = XoaNghieng(87:691,177:747);
% Kich thuoc Pixel laf 240nm
[Sy,Sx] = size(Chom);
x = 1:Sx;
y = 1:Sy;
[xx, yy] = meshgrid(x*238,y*238);
zz = Chom;
%Xoas cac thanh phan khong thuoc mat cau
indices = find(Chom < 10);
zz(indices) = [];
xx(indices) = [];
yy(indices) = [];
xxx = xx(:);
yyy = yy(:);
zzz = zz(:);
%Áp d?ng hàm Fit mat cau:
[xc, yc, zc, R] = sumith_fit(xxx,yyy,zzz);
Rchom = R/1000000z
Bình luận gần đây