
close all
clearvars
addgeopaths




%% --- SETTINGS -----------------------------------------------
%
%
settings=[];
settings.regionname='ZA79'
settings.timerange=[365-16*2 365+16*2];
%settings.timerange=[365*.9 365*3+16];%365*1.5]; %min and max temporal range between images
%settings.timerange=[365-16 365+16];%365*1.5]; %min and max temporal range between images
settings.maxvelocity=7; %[m/day] used to calculate a search window size. Better too large than too small
settings.chipsize=100; %*15m
settings.roi=[79.7068 -20.8740; 79.5486 -23.2031; 79.5963 -24.9170; 79.5167 -27.2900; 78.2380 -26.4111; 78.5430 -23.1152; 78.6991 -21.0718; 78.5518 -21.0938; 78.7206 -19.3579; 79.1713 -19.0942; 79.1009 -20.3467; 79.2577 -21.7090; 79.4041 -20.5225; 79.3474 -19.2920; 79.7695 -18.9844; 79.7577 -19.6216; 79.9205 -19.6216; 79.9052 -20.5664];

%use measures grid!!!
measures=[];
[measures.Vx,measures.mx,measures.my,Im]=geoimread('c:\users\ag\hugedata\measures\greenland_vel_mosaic500_2015_2016_vx_v2_comp.tif',settings.roi(:,1),settings.roi(:,2));
[measures.Vy,measures.mx,measures.my,Im]=geoimread('c:\users\ag\hugedata\measures\greenland_vel_mosaic500_2015_2016_vx_v2_comp.tif',settings.roi(:,1),settings.roi(:,2));
[MX,MY]=meshgrid(measures.mx,measures.my);


%settings.filepatterns={'_00[789]003_' '_01[01]003_'};
settings.filepatterns={'_00[89]003_' '_01[01]003_'}; %ignore 7's here at first

settings.trackgrid=[];
[settings.trackgrid.lat,settings.trackgrid.lon]=mapproj(Im,MX,MY);


sz=3;
gf=zeros(sz);gf([1 end])=[1 -1];
% gf=conv2(gf,fspecial('disk',sz/sqrt(2)));
gf=gf+1i*rot90(gf);
fang=@(A)exp(1i*atan2(imag(A),real(A)));
myfilter=@(A)fang(imfilter(A,gf,'replicate'));
%  myfilter=@(A)exp(1i*atan2(imfilter(A,gf,'replicate'),imfilter(A,rot90(gf),'replicate')));

outputfolder='output';
mkdir(outputfolder)

%%------------------------------------------------------

%find all images....

% if false
%     r=recursivedir('D:\hagen\extractc1\','*.tif');
%     
%     images=landsatparse(r{1});
%     for ii=2:length(r)
%         imgA=landsatparse(r{ii});
%     end
% else
    try
        images=landsatindex(settings.filepatterns);
        save('output\images.mat','images')
    catch
        load('output\images.mat','images');
        warning('fallback image list')
    end
    
% end

[~,six]=unique([[images.path]' [images.row]' [images.date]'],'rows'); %unique and sorting.
images=images(six);

Nimages=length(images);

dt=diff([images.date]');dt(dt<0)=nan;
dispmtx([[images.path]' [images.row]'  [dt;nan]],'%7.0f',{'path' 'row' 'dt'},cellstr(datestr([images.date]')))

for kk=1:Nimages-1
    for ii=1:Nimages-kk
        jj=min(ii+kk,Nimages);
        imgA=images(ii);
        imgB=images(jj);
        fA=imgA.filename; fB=imgB.filename;
        mtla=cachedcall(@MTL_parser,{regexprep(fA,'_B\d+.TIF','_MTL.txt')});

        tA=imgA.date; tB=imgB.date;
        dt=tB-tA;
        if ~((imgA.path==imgB.path)&&(imgA.row==imgB.row))
            continue
        end
        if ~(all(imgA.sat==imgB.sat))
            continue
        end
        if ~( (dt>=min(settings.timerange)) && (dt<=max(settings.timerange)))
            continue
        end
        
        %pairid=sprintf('%s_%s-%s',settings.regionname,datestr(tA,'yyyymmdd'),datestr(tB,'yyyymmdd'));
        [ya,ma,da]=datevec(tA);jda=tA-datenum(ya,1,0);
        [yb,mb,db]=datevec(tB);jdb=tB-datenum(yb,1,0);
        
        if ismember(ma,6:8) || ismember(mb,6:8) %skip fastest flow in mid july... %TODO later
            continue
        end
        
        pairid=sprintf('%s_%04i%03i-%04i%03i_%03i%03i',settings.regionname,ya,jda,yb,jdb,imgA.path,imgA.row);
        
        foutdata=fullfile(outputfolder,sprintf('V_%s.mat',pairid));
        fprintf('\n\nPair: %s   dt=%.0f\n',pairid, dt)
        
        
        
        if ~isempty(dir(foutdata)) %&& (~strcmpi(pairid,'KNSSMALL005015_2013313-2014044'))
            fprintf('already processed - skipping!\n')
            continue
        end
        
        %--------------------------- Do the tracking --------------
        
        
        disp('loading');drawnow;
        try
            mtla=cachedcall(@MTL_parser,{regexprep(fA,'_B\d+.TIF','_MTL.txt')});
            mtlb=cachedcall(@MTL_parser,{regexprep(fB,'_B\d+.TIF','_MTL.txt')});
            thres=15;
            if mtla.IMAGE_ATTRIBUTES.CLOUD_COVER_LAND>thres || mtlb.IMAGE_ATTRIBUTES.CLOUD_COVER_LAND>thres
                continue
            end
            [A,xa,ya,Ia]=cachedcall(@geoimread,{fA,settings.roi(:,1),settings.roi(:,2),4000});
            [B,xb,yb,Ib]=cachedcall(@geoimread,{fB,settings.roi(:,1),settings.roi(:,2),4000});
            if (xa(2)-xa(1))>(xb(2)-xb(1))
                B=imfilter(B,ones(2)/4);
                B=B(1:2:end,1:2:end);
                xb=(xb(1:2:end)+xb(2:2:end))/2;
                yb=(yb(1:2:end)+yb(2:2:end))/2;
            elseif (xa(2)-xa(1))<(xb(2)-xb(1))
                A=imfilter(A,ones(2)/4);
                A=A(1:2:end,1:2:end);
                xa=(xa(1:2:end)+xa(2:2:end))/2;
                ya=(ya(1:2:end)+ya(2:2:end))/2;
            end
        catch
            warning('skipping - cant read file')
            continue
        end
        
        if ~strcmpi(Ia.Projection,Ib.Projection)
            disp('not same projection')
            continue
        end
        
        fig('thumbnail');clf
        subaxis(1,2,1,'s',0,'m',0)
        imshow(imresize(A,.1));
        subaxis(2);
        imshow(imresize(B,.1));
        
        deltax=xa(2)-xa(1); deltay=ya(2)-ya(1);
        A=im2single(A);
        B=im2single(B);
        mapoffset=[(xa(1)-xb(1))/deltax (ya(1)-yb(1))/deltay];
        
        
        
        BG=A;BG=BG/max(BG(:));BG=im2uint8(BG);
        
        disp('filtering');drawnow;
        A(A==0)=nan;B(B==0)=nan;
        A=imnoisefill(A,35);
        B=imnoisefill(B,35);
        A=myfilter(A); B=myfilter(B);
        A=A./max(abs(A(:))); B=B./max(abs(B(:)));
        A(isnan(A))=0;B(isnan(B))=0;
        
        
        
        %chipsz=settings.chipsize;
        chipsz=settings.chipsize*15/deltax;
        super=1; %no need to super sample
        
        
        
        sz=size(A);
        
        [pu,pv]=projfwd(Ia,settings.trackgrid.lat,settings.trackgrid.lon);
        pu=round(interp1(xa,1:length(xa),pu,'linear','extrap'));
        pv=round(interp1(ya,1:length(ya),pv,'linear','extrap'));
        
        %REDUCE GRID
        %             pu=pu(1:5:end,1:5:end);
        %             pv=pv(1:5:end,1:5:end);
        
        dpu=pu*0;dpv=pv*0; %initial guess for displacement
        
        px=interp1(xa,pu,'linear','extrap');
        py=interp1(ya,pv,'linear','extrap');
        if ~isempty(settings.roi)
            [roi,roi(:,2)]=projfwd(Ia,settings.roi(:,1),settings.roi(:,2));
            dpu(~inpolygon( px , py ,roi(:,1),roi(:,2)))=nan;
            dpu(isnan(interp2(A,pu,pv)))=nan;
        end
        
        
        
        maxdisplacement=round(settings.maxvelocity*dt/deltax);
        
        
        %         rockpoly=[568004 7134316;569572 7131690;572550 7130632;572746 7129574;574275 7129261;578899 7133885;583759 7135805;582818 7138588;580898 7141449;577018 7139019;572315 7135805;568592 7135257];
        %         inrock=inpolygon(px,py,rockpoly(:,1),rockpoly(:,2));
        %         rockpoly=[556874 7135845;561185 7132866;563576 7130319;564869 7127144;567103 7123774;573452 7117621;569493 7115230;563889 7119776;556678 7115387];
        %         inrock=inrock | inpolygon(px,py,rockpoly(:,1),rockpoly(:,2));
        
        options=[];
        options.pu=pu;
        options.pv=pv;
        options.showprogress=true;%{datestr(tA) datestr(tB)}; %true
        options.searchwidth=chipsz+maxdisplacement*2; %reduced search because of initial guess
        options.templatewidth=chipsz;
        options.supersample=1; %supersample the input images
        options.initialdu=dpu+mapoffset(1);
        options.initialdv=dpv+mapoffset(2);
        if ~isreal(A)
            options.method='CCF';
        end
        
        
        
        
        
        disp('tracking');drawnow;
        [dpu,dpv,C,Cnoise,pu,pv]=templatematch(A,B,options);
        
        dx=(dpu-mapoffset(1))*deltax/dt;
        dy=(dpv-mapoffset(2))*deltay/dt;
        
        dv=sqrt(dx.^2+dy.^2);
        
        snr=(C./filter2(ones(3)/9,Cnoise));
        snr(isnan(snr))=C(isnan(snr))./Cnoise(isnan(snr));
        
        
        
        disp('saving')
        save(foutdata,'img*','MX','MY','Im','dx','dy','dt','t*','C','Cnoise','snr','settings','options')
        
        
        
        fig('B');clf
        axes('position',[0 0 1 1])
        showimg(xa,ya,BG);
        xlim([xa(1) xa(end)]);
        ylim([min(ya) max(ya)]);
        text(1,0,sprintf('%s   dt=%.0f ',pairid, dt),'interpreter','none','backgroundcolor',[1 1 1],'verticalalignment','bottom','horizontalalignment','right','units','normalized','margin',0.001);
        %title(sprintf('%s   dt=%.0f',pairid, dt),'interpreter','none')
        xl=xlim;
        hold on;
        
        
        w=3;
        inconsistent=sqrt((dx-medfilt2(dx,[w w])).^2+(dy-medfilt2(dy,[w w])).^2);
        keep=(snr>1) & (inconsistent<settings.maxvelocity*.7) & ~isnan(dpu); %change the thresholds
        
        %         rockpoly=[71.11113 -26.47324;71.12332 -26.45414;71.14027 -26.44013;71.14985 -26.44649;71.15546 -26.47205;71.15499 -26.49667;71.09739 -26.55081;71.09387 -26.50072;];
        %         [rockpoly(:,1),rockpoly(:,2)]=projfwd(Ia,rockpoly(:,1),rockpoly(:,2));
        %         inrock=inpolygon(px,py,rockpoly(:,1),rockpoly(:,2));
        %         dxrock=nanmedian(dx(keep&inrock));
        %         dyrock=nanmedian(dy(keep&inrock));
        %         dx=dx-dxrock;
        %         dy=dy-dyrock;
        dv=sqrt(dx.^2+dy.^2);
        
        
        if sum(keep(:))==0
            q=[];
        end
        
        % keep=(snr>8)&(abs(dv)<25); %change the thresholds
        %fastscatter(pu(keep),pv(keep),dv(keep),'markersize',7);
        h=alphawarp(px,py,dv,keep*.5+(~isnan(dv))*.2);
        
        
        
        %plot(px(keep),py(keep),'k.')
        quiver(px(keep),py(keep),dx(keep),dy(keep),0.5,'k')
        hslcolormap('bR',1,[.3 .85])
        
        clim=[0 settings.maxvelocity];
        caxis(clim)
        
        
        axes('position',[0.05 0.01 0.2 0.025])
        xx=linspace(clim(1),clim(2),100);
        imagesc(xx,0.5,xx)
        caxis(clim);
        set(gca,'ytick',[],'xaxislocation','top')
        
        
        sz=[5 4];
        savefigure(fullfile(outputfolder,sprintf('V_%s.png',pairid)),'s',sz,'dpi',200);
        %set(gcf,'PaperUnits','inches','PaperPositionMode','manual','PaperSize',sz,'PaperPosition',[0 0 sz]);
        %print('-dpng',fullfile(outputfolder,sprintf('V_%s.png',pairid)),'-r200')
        
    end
end



