interpolarefarav4pct_control.m

%date intrare
clc;clear;
fisier = 'GNSS_ACMTC2.txt';
fisier2 = 'GNSS_ACMTC2.txt';
%sfarsit date intrare
digits(3);
pct=load(fisier);
pctctr=load(fisier2);

x=pct(:,3);
y=pct(:,2);
z=pct(:,4);


Z=pctctr(:,4);
QX(:,1)=pctctr(:,3);
QX(:,2)=pctctr(:,2);

[X1,Y1] = meshgrid(697110:0.2:697160,631708:0.2:631758);
Z1=griddata(x,y,z,X1,Y1,'linear'); 
figure;
surf(X1,Y1,Z1);hold on, plot3(x,y,z,'bo'), hold on;  
title('Metoda de interpolare liniara','fontweight','b'); colorbar; hold off;
xlabel('Est [m]'); ylabel('Nord [m]'); zlabel('Altitudini normale [m]');
QZ1=griddata(X1,Y1,Z1,QX(:,1),QX(:,2),'linear');


[X2,Y2] = meshgrid(697110:0.2:697160,631708:0.2:631758);
Z2=griddata(x,y,z,X2,Y2,'nearest'); 
figure;
surf(X2,Y2,Z2); hold on, plot3(x,y,z,'bo'), hold on;  
title('Metoda de interpolare a vecinului cel mai apropiat','fontweight','b'); colorbar; hold off;
xlabel('Est [m]'); ylabel('Nord [m]'); zlabel('Altitudini normale [m]');
QZ2=griddata(X2,Y2,Z2,QX(:,1),QX(:,2),'nearest');


[X3,Y3] = meshgrid(697110:0.2:697160,631708:0.2:631758);
Z3=griddata(x,y,z,X3,Y3,'cubic');
figure;
surf(X3,Y3,Z3); hold on, plot3(x,y,z,'bo'), hold on;  
title('Metoda de interpolare spline cubica','fontweight','b'); colorbar; hold off;
xlabel('Est [m]'); ylabel('Nord [m]'); zlabel('Altitudini normale [m]');
QZ3=griddata(X3,Y3,Z3,QX(:,1),QX(:,2),'cubic');

[X4,Y4] = meshgrid(697112:0.2:697160,631708:0.2:631758);
Z4  = interp2(X1, Y1, Z1, X4, Y4, 'bicubic');
figure
surf(X4, Y4, Z4); hold on, plot3(x,y,z,'bo'), hold on;  
title('Metoda de interpolare spline bicubica','fontweight','b'); 
xlabel('Est [m]'); ylabel('Nord [m]'); zlabel('Altitudini normale [m]');

% shading flat; 
colormap(jet); view(0, 90)
hold on; 
% plot3(X1, Y1, 10*ones(size(X1)), 'k.', 'MarkerSize', 20)
colorbar
print -dpng -r0.5 Interpolare_Bicubica_Acoperis_Cadastru.png
QZ4=interp2(X4,Y4,Z4,QX(:,1),QX(:,2),'bicubic');


A1=reshape(X1,[],1);
B1=reshape(Y1,[],1);
C1=reshape(Z1,[],1);
D1 = horzcat(A1,B1,C1);

E1=D1(isfinite(D1(:,3)), :);

A2=reshape(X2,[],1);
B2=reshape(Y2,[],1);
C2=reshape(Z2,[],1);
D2 = horzcat(A2,B2,C2);

E2=D2(isfinite(D2(:,3)), :);

A3=reshape(X3,[],1);
B3=reshape(Y3,[],1);
C3=reshape(Z3,[],1);
D3 = horzcat(A3,B3,C3);

E3=D3(isfinite(D3(:,3)), :);

A4=reshape(X4,[],1);
B4=reshape(Y4,[],1);
C4=reshape(Z4,[],1);
D4 = horzcat(A4,B4,C4);

E4=D4(isfinite(D4(:,3)), :);

% 
% figure;
% plot(x,y,'g.','MarkerSize',20);
% hold on
% axis('equal')
% title('Acoperis cladire "CADASTRU-NORD"','fontweight','b')  
% grid
% xlabel('East [m]')
% ylabel('Nord [m]')
% hold on
% plot(X1,Y1,'bx')
% axis('equal')
% hold on;
% 
% raz=0.5;
% 
% for i=1:size(E4)
%          circle(E4(i,1),E4(i,2),raz);
% end

% figure;
% plot3(x,y,z,'ro')
% hold on
% axis('equal')
% title('puncte mun. Iasi')  
% grid

% figure;
% Tri = delaunay(x,y);
% trisurf(Tri,x,y,z,45);

Dif1=Z-QZ1
Dif2=Z-QZ2
Dif3=Z-QZ3
Dif4=Z-QZ4



Dif(:,1)=Dif1;
Dif(:,2)=Dif2;
Dif(:,3)=Dif3;
Dif(:,4)=Dif4;


dlmwrite('diferente_pct_control(0.2)GPS.txt',Dif,'delimiter','\t','precision',3);

% dlmwrite('coord_noduri_gridlidarBICUBIC.txt',E5,'delimiter','\t','precision',9);