clear all; close all;

s1 = load('geo_090908-091024_test.txt');
s2 = load('geo_100727-101027_test.txt');
s3 = load('geo_080721-090908_test.txt');
t = load('dem.txt');

[x,y] = size(s1); nPix = x.*y;

ss1 = reshape(s1,nPix,1);
ss2 = reshape(s2,nPix,1);
ss3 = reshape(s3,nPix,1);
tt = reshape(t,nPix,1);

%figure; plot(ss1,tt,'ro');

%figure; subplot(1,4,1); imagesc(s1,[-15 15]);
%subplot(1,4,2); imagesc(s2,[-15 15]);
%subplot(1,4,3); imagesc(s3,[-15 15]);
%subplot(1,4,4); imagesc(t);

s = s3;
t = s2;

y1 = fft2(detrend(s));    %s in freq
y2 = fft2(detrend(t));    %t in freq


%--coherence
COR12 = ((y1).*conj((y2)));
COR11 = ((y1).*conj((y1)));
COR22 = ((y2).*conj((y2)));

filtSize = 10; 
myFilt = ones(filtSize) ./ filtSize.^2;

COR12 = imfilter(COR12, myFilt);
COR11 = imfilter(COR11, myFilt);
COR22 = imfilter(COR22, myFilt);
RCOR = real(COR12)./sqrt(COR11.*COR22);
ICOR = imag(COR12)./sqrt(COR11.*COR22);
COR = RCOR.^2+ICOR.^2;
COH = ifft2((COR12)./sqrt(COR11.*COR22),'symmetric');
%COH2 = abs(COR12./sqrt(COR11.*COR22));
%COH3 = atan(ICOR./RCOR);
%COH4 = sqrt(ICOR.^2+RCOR.^2);
%COR4 = ifft2(COH4,'symmetric');

filt = 1 - COR;
 
%inverse FFT
meanST = (y1+y2)/2;
fs = real(ifft2(y1.*COR));
ft = real(ifft2(y2.*filt));

%residual
ress = s - fs;
rest = t - ft;

figure; subplot(2,3,1); imagesc(s,[-15 15]);
subplot(2,3,2); imagesc(fs,[-15 15]); 
subplot(2,3,3); imagesc(ress,[-15 15]);
subplot(2,3,4); imagesc(t,[-15 15]); 
subplot(2,3,5); imagesc(ft,[-15 15]); 
subplot(2,3,6); imagesc(rest,[-15 15]);

%{
%%%%%%%%cut profile
%first point of the profile
X1 = 48;       
Y1 = 54;
%second point of the profile
X2 = 321;       
Y2 = 195;
%the lateral distance (m)
d = 10;

fig1 = profile_rot(s,X1,X2,Y1,Y2,d);
fig2 = profile_rot(fs,X1,X2,Y1,Y2,d);


%plot fft
f3=fft2(s3);
imagesc(fftshift(log(f3.*conj(f3))));
%}








