clear all 
close all

%cette opération permet d´avoir une matrice vide de la taille  
fid = fopen('mnt_aragon_2007_utm.asc', 'r'); 
mnt = textscan(fid, '%f %f %f', 'headerlines', 0, 'delimiter',',','CollectOutput', 0); 
fclose(fid); 

% Dans MNT se retrouvent les valeurs xyz dans les colonnes 1,2 et 3.

%Pour avoir un MNT régulier, il est utile d´avoir un grid régulier avec la
%fonction meshgrid.A partir des valeurs maximales et minimales,  meshgrid construit 
%une matrice qui reproduit a chaque endroit les valeurs x et y correspondant.       

minx=min(mnt{1}); 
maxx=max(mnt{1}); 
miny=min(mnt{2}); 
maxy=max(mnt{2}); 

% pour avoir un Mnt régulier avec des nombres entiers, il faut arondir les
% valeurs extremes.    
minx=floor(minx); 
miny=floor(miny); 
maxx=ceil(maxx); 
maxy=ceil(maxy); 

% stockage des trois colonnes spéarément
X=mnt{1};
Y=mnt{2};
Z=mnt{3};

% il faut déterminer le nb lignes/ nb colonnes. Lire simultanément les
% colonnes X,Y. =OK
dx=maxx-minx
dy=maxy-miny
rows=ceil(dy/100)
cols=ceil(dx/100)

%***********************************************************************

I=zeros(rows,cols);
for i=1:size(X)
I(ceil((Y(i)-miny)/100), ceil((X(i)-minx)/100))=Z(i);
end

% Création de l'image initiale avec un pas de 100

%Pour chq x,y faire le calcul pour revenir en coord image?????????????????????? 
%(on arrondit au plus proche voisin) et affecter au
% pixel correspondant la valeur du z à partir de la troisième colonne.

% affichage image
I=double(I);

figure,imagesc(I),axis image,colormap('jet'),title('image initale res=100');

%Création image interpolée
%demander à l'utilisateur le pas d'échantillonage
%res = résolution souhaité
res=50
s=100/res
ro=s*rows
co=s*cols
Ifinal=zeros(ro,co);
for i=1:1:ro  %correction NC :l'ordre des boucles
    for j=1:1:co
       x=(1/double(s))*j;
       y=(1/double(s))*i;
       Ifinal(i,j)=pv(I,x,y,'linear'); % correction NC: utilisation de la fonction pv 
    end
end
% correction NC: le pb à l'affichage, il faut passer l'image en double pour
% la voir en couleur 
Ifinal=double(Ifinal);
figure,imagesc(Ifinal),axis image,colormap('jet'),title('image interpolée');

%pour chaque pt de l'image interpolée -> calcul pour retrouver le point(en
%double) x,y dans l'image initiale



% --> interpolation ( 3 types: plus proche voisin, interpolation bilinéaire
% et inverse distance.



%******************************************************
% Exmeple d'utilisation de la fonction d'interpolation de matlab
%******************************************************
% Construction d´une grille des coordonnées avec un intervalle de 20 mètres.    
% coordonnées interpolées
%[xmesh, ymesh]=meshgrid(minx:20:maxx, miny:20:maxy); 

%ZI = griddata(X,Y,Z,xmesh,ymesh,'cubic');
%figure,mesh(xmesh,ymesh,ZI),hold on;
%plot3(X,Y,Z,'o'), hold off
%figure, surf(xmesh,ymesh,ZI);

%*******************************************************