Plane Fit and Sphere Fit in Matlab

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