%% LoadGoogleEarthImages.m

% This script loads Google Earth aerial images and applies a geometric
% transformation to obtain the 44x49.5m focus area in orthogonal view. It
% furthermore computes the (relative) VDVI-values of this matrix.

% Date: 19 December 2019
% 
% Dataset creator and contact person:
% Roeland C. van de Vijsel
% roeland.van.de.vijsel[at]nioz.nl
% NIOZ Royal Netherlands Institute for Sea Research (www.nioz.nl)
% Yerseke, The Netherlands
% 
% This data belongs to the dataset that can be downloaded from the NIOZ
% Data Repository (https://dataverse.nioz.nl) and/or can be found as 
% DOI: 10.25850/nioz/7b.b.m (https://doi.org/10.25850/nioz/7b.b.m).
% 
% This dataset is part of the following publication:
% van de Vijsel, R.C., van Belzen, J., Bouma, T.J., van der Wal, D.,
% Cusseddu, V., Purkis, S.J., Rietkerk, M. & van de Koppel, J. (2019).
% Estuarine biofilm patterns: Modern analogues for Precambrian
% self-organization. Earth Surface Processes and Landforms. 
% DOI: 10.1002/esp.4783 (https://doi.org/10.1002/esp.4783).
% 
% Whenever any of these data (scripts, figures, raw and processed data) are
% used elsewhere, please refer correctly to this publication and its
% authors and please cite the dataset. Options for citing the dataset can
% be found on https://dataverse.nioz.nl. Dataset and publication published
% under the Creative Commons Attribution License. 

clear all
close all
clc

on = 1; off = 0; % Define on/off switches

SAVE = off; % Whether or not to save output files and plots

%% 1.1 Load image and manually select the cornerpoints (with the cursor, export cursor data)
Image0 = imread('input/2004_with_49.5m_long_legend.jpg');
figure(1),
imagesc(Image0);colorbar; axis image
disp('1.1 done')

%% 1.2 Register the image onto a 44 x 49.5m "orthogonal reference matrix"
moving = Image0; % The unregistered image that has to be registered
fixed = zeros(440,495); % The registered "reference" matrix (44x49.5m with 0.10m grid spacing)

cpselect(moving,fixed) % Use control point selection to do the registering: manually select middle of corner "bullets" in "moving" and outermost corners in "fixed"
disp('1.2 done')

%% 1.3 Store control points
%%% "movingPoints" = [3863.37433337891,2503.03459592507;3828.25220560031,39.0657204961003;1065.61691508273,82.3590865581841;1104.69179133103,2540.42194092827]
%%% "fixedPoints" = [495.500000000000,440.500000000000;495.500000000000,0.500000000000000;0.500000000000000,0.500000000000000;0.500000000000000,440.500000000000]
disp('1.3 done')

%% 1.4 Perform the actual image transformation
unregistered = Image0; % Unregistered aerial image
ortho = fixed; % Registered reference matrix
t_concord = fitgeotrans(movingPoints,fixedPoints,'projective'); % Coordinate transformation (assuming projective distortion) to register the unregister drone image into the same coordinate system as the ortho (DEM)
Rfixed = imref2d(size(ortho));
registered = imwarp(unregistered,t_concord,'OutputView',Rfixed); % 'OutputView' ensures the registered image will be aligned for elementwise comparison with the orthophoto

Image0reg = registered;

disp('1.4 done')

%% 1.5 Check output
figure(2)
subplot(1,2,1)
imagesc(Image0); axis image

subplot(1,2,2)
imagesc(Image0reg); axis image

size(Image0reg)
disp('1.5 done')

%% 1.6 Apply this transformation/selection procedure to the "real" images (without labels) and oompute VDVI-matrix

yylist = [2004,2009,2012,2013,2015:2017];

for yy = yylist
    
    Image{yy} = imread(['input/',num2str(yy),'.jpg']); % Full photo for this year
    
    % Apply transformation
        unregistered = Image{yy}; % Unregistered aerial image
        ortho = fixed; % Registered reference matrix
        t_concord = fitgeotrans(movingPoints,fixedPoints,'projective'); % Coordinate transformation (assuming projective distortion) to register the unregister drone image into the same coordinate system as the ortho (DEM)
        Rfixed = imref2d(size(ortho));
        registered = imwarp(unregistered,t_concord,'OutputView',Rfixed); % 'OutputView' ensures the registered image will be aligned for elementwise comparison with the orthophoto

    ImReg{yy} = registered; % Registered and cropped image
    
    VDVI{yy} = (2*cast(ImReg{yy}(:,:,2),'double') ... % Converted to VDVI
                - cast(ImReg{yy}(:,:,1),'double')...
                - cast(ImReg{yy}(:,:,3),'double'))...
                 ./(2*cast(ImReg{yy}(:,:,2),'double') ...
                    + cast(ImReg{yy}(:,:,1),'double')...
                    + cast(ImReg{yy}(:,:,3),'double'));
    
    % Save output image
    figure(yy)
    imagesc(ImReg{yy});axis image
    set(gca,'XTick',[],'YTick',[])
    if SAVE == on
        saveas(gcf,[output/num2str(yy),'Registered.jpg'])
        saveas(gcf,[output/num2str(yy),'Registered.pdf'])
        saveas(gcf,[output/num2str(yy),'Registered.png'])
    end
    
disp(['Registering ',num2str(yy),' done'])
end
disp('1.6 done')

%% 1.7 Determine relative VDVI, i.e. subtract background VDVI

for yy = yylist

    Z = VDVI{yy}; % Total VDVI
    [row col]=size(Z);
    xlin=linspace(0,1,col);
    ylin=linspace(0,1,row);
    [X Y]=meshgrid(xlin,ylin);
    x=X(:);y=Y(:);
    z=Z(:);
 
        xid=[ones(size(x)) x y x.^2 x.*y y.^2 x.^3 (x.^2).*y x.*(y.^2) y.^3 x.^4 (x.^3).*y (x.^2).*(y.^2) x.*(y.^3) y.^4]; % Fourth-order 2D-regression (incl. cross-terms)
        [b,bint,r,rint,stats] = regress(z(:),xid); % Find best fit (= background VDVI)
        T4=b(1)+b(2)*X+b(3)*Y+b(4)*X.^2+b(5)*X.*Y+b(6)*Y.^2+b(7)*X.^3+b(8)*(X.^2).*Y+b(9)*X.*(Y.^2)+b(10)*Y.^3+b(11)*X.^4+b(12)*(X.^3).*Y+b(13)*(X.^2).*(Y.^2)+b(14)*X.*(Y.^3)+b(15)*Y.^4; % Background VDVI
        Zdet4=Z-T4; % "Detrended", i.e. relative VDVI
        stats4 = stats;
        Rsquared4 = stats4(1) % Show R-squared

    VDVIdet{yy} = Zdet4; % Relative ("detrended") VDVI
    VDVIT{yy}   = T4; % Trend in VDVI

end
disp('1.7 done')

%% 1.8 Plot total and relative VDVI

for yy = yylist
    
    figure(yy + 1000)

    subplot(2,2,1)
    imagesc(VDVI{yy});
    axis image
    colorbar
    caxis([-0.1 0.1])

    subplot(2,2,2)
    imagesc(VDVIT{yy});
    axis image
    colorbar
    caxis([-0.1 0.1])

    subplot(2,2,3)
    imagesc(VDVIdet{yy});
    axis image
    colorbar
    caxis([-0.1 0.1])
    
end
disp('1.8 done')

%% 1.9 Save to output file

if SAVE == on
    
    for yy = yylist
        
        VI = VDVI{yy};
        VIT = VDVIT{yy};
        VId = VDVIdet{yy};
        
        save(['output/VDVI',num2str(yy),'.txt'],'VI','-ascii')
        save(['output/VDVI_Trend',num2str(yy),'.txt'],'VIT','-ascii')
        save(['output/VDVIdet',num2str(yy),'.txt'],'VId','-ascii')
    end
    
end
disp('1.9 done')
