% Haerer, Bernhardt, Schulz and Corripio (2012)
% "PRACTISE - Photo Rectification And ClassificaTIon SoftwarE"
%
    % incorporating adapted versions of
    %
    % Corripio, 2004 
    % "Snow surface albedo estimation using terrestrial photography"
    % (Int J Rem Sens)
    % 
    % Wang, Robinson and White, 2000 
    % "Generating Viewsheds without Using Sightlines" (PE&RS)
    % 
    % Tolson and Shoemaker, 2007 
    % "Dynamically dimensioned search algorithm for computationally
    % efficient watershed model calibration" (WRR)
%
%   written by
%   Stefan Härer (LMU Munich)
%   08/2012
%   contact: s.haerer@iggf.geo.uni-muenchen.de
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Name:       PRACTISE
%   Purpose:    Main file of PRACTISE
%   Comment:    This file starts and controls the complete software tool,
%               be careful with any changes.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% input %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; close all;
disp('Starting PRACTISE - Photo Rectification And ClassificaTIon SoftwarE')
%
disp('Using the following options:')
tic
% load input configuration
Input_PRACTISE
% display chosen switches
%   viewshed
if vs == false
    disp('- use existing viewshed (load from ASCII-file [cf. ESRI ArcGIS])');
elseif vs == true
    disp('- generate viewshed with specified camera angles');
elseif ~exist('vs', 'var')
    error('No viewshed switch specified, check parameter file');
else
    error('Viewshed switch option not available, check parameter file');
end
%   ground control points       
if gs == false
    disp('- w/o ground control points');
elseif gs == true
    disp('- w ground control points');
elseif ~exist('gs', 'var')
    error('No ground control points switch specified, check parameter file');
else
    error('Ground control points switch option not available, check parameter file');
end
%   optimization
if os == false
    disp('- w/o DDS optimization');
elseif os == true & gs == true
    disp('- w DDS optimization');
elseif os == true & gs == false
    error('Optimization switch on, but ground control points switch off');
elseif ~exist('os', 'var')
    error('No optimization switch specified, check parameter file');
else
    error('Optimization switch option not avaialable, check parameter file');
end
%   image
if is == false
    disp('- classify a single image');
elseif is == true & os == false
    disp('- classify all images in a specific folder');
elseif is == true & os == true
    error('Image switch on and optimization switch on, this does not work');
elseif ~exist('is', 'var')
    error('No image switch specified, check parameter file');
else
    error('Image switch option not available, check parameter file');    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% check %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Camera positions (viewpoint and targetpoint) defined?
if ~exist('cam', 'var')
    error('Camera viewpoint: longitude and latitude position is not defined.');
end
if size(cam,2)<2
    error('Camera targetpoint: longitude and latitude position is not defined.');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% load %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Loading and preparing input')
% load
%   DEM
fid=fopen(fin_demW);
if fid==-1
    error('DEM file could not be found.');
end
for i=1:6
    dummy=fgetl(fid);
    header_W(i,:)=cellstr(dummy);
%     headern_W(i,:)=cellstr(dummy(1:14));
    headerv_W(i,:)=str2num(dummy(15:end));
end
clear dummy i
demWrcz = fscanf(fid, '%f', [headerv_W(1,1),inf])';
demWrcz(demWrcz==headerv_W(6,1))=NaN;
fclose(fid);
clear ans fid
%   Viewshed (if ASCII-file provided)
if vs==false
    fid=fopen(fin_viewW);
    if fid==-1
        error('Viewshed file could not be found.');
    end
    for i=1:6
        dummy=fgetl(fid);
        if dummy==char(header_W(i,:))
            clear headerdummy
        else
            error('Viewshed and DEM header not equal, please check.');
        end
    end
    clear dummy i
    viewW = fscanf(fid, '%f', [headerv_W(1,1),inf])';
    fclose(fid);
    clear ans fid
end
%   Ground control points (if ASCII-table provided)
if gs==true
    fid=fopen(fin_gcpW);
    if fid==-1
        error('Ground control points file could not be found.');
    end
    cell_gcpW=textscan(fid, '%f %f %f %u %u %s', N_gcpW+1, 'headerlines', 1);
    for i=1:5 % longitude, latitude, altitude, photo rows and photo cols
        gcpW(i,:)=cell_gcpW{1,i}';
    end
    name_gcpW(1,:)=cell_gcpW{1,6}';
    clear i %cell_gcpW
    fclose(fid);
    clear ans fid 
end
%   Photograph (RGB)
%       single photo
if  is==false
    if exist('fin_image', 'var') 
        photo = imread(fin_image); 
        N_images=1;
    else
        error('Photo file could not be found.')
    end
%       multiple photos
%           loading while loop in projection part 
elseif is==true
    if exist('fin_imformat', 'var') & exist('fin_imfolder', 'var')
        list=dir(strcat(fin_imfolder,'*',fin_imformat));
        for i=1:length(list)
            fin_images{i,1}=list(i).name;
        end
        clear i list
        N_images=length(fin_images);
    else
        error('Photo folder and format could not be found.')
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%% DDS optimization %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% create DDS parameters
%   seed (X0)
X0=[cam(1,1), cam(2,1), cam_off, cam_rol, cam(1,2), cam(2,2), cam_foc];
if os==true | gs==true
% save intermediate file for RMSE function
save('DDSintermediatedelete.mat', 'gcpW', 'pix_c', 'pix_r', 'cam_wid', ...
    'cam_hei') 
end
if os==true
    toc
    disp('Starting DDS optimization')
%   upper boundary (UB) using UBD
    UB=X0+UBD;
%   lower boundary (LB) using LBD
    LB=X0+LBD;
% plot photograph and hold on for visual investigation of the optimization
    figure1=figure(1);
    %
    set(figure1, 'Color',[1 1 1], 'InvertHardcopy','off', ...
      'PaperPosition',[0 0 50 35], 'PaperSize',[50 35], 'PaperType','<custom>');
    axes1 = axes('FontSize',15.91, 'FontWeight','demi', 'Layer','top',...
      'TickDir','out', 'XAxisLocation','top', 'YDir','reverse', ...
      'DataAspectRatio', [1,1,1], 'Parent',figure1);
    axis(axes1,[0.5 pix_c+0.5 0.5 pix_r+0.5]);
    xlabel(axes1,'Pixel column','FontSize',18.185);
    ylabel(axes1,'Pixel row','FontSize',18.185);
    box(axes1,'on');
    hold(axes1,'all');
    %
    image(photo, 'Parent',axes1);                                                    
    hold on        
    plot(gcpW(4,:)+0.5, gcpW(5,:)+0.5, 'gx', 'LineWidth', 2, 'MarkerSize',10)
    %
    clear figure1 axes1
% call DDS optimization
[Xopt, Fopt]=DDS_PRACTISE('DDS_RMSE_PRACTISE', X0, LB, UB, DDS_R, DDS_MaxEval, ...
    fout_folder, demWrcz, headerv_W);
% save original camera parameters (incl. altitude)
[cam_orig, cam_rc_orig]=Cam_PRACTISE(cam, cam_off, demWrcz, headerv_W);
% forward new parameters to existing variables
cam(1,1)=Xopt(1,1);
cam(2,1)=Xopt(1,2);
cam_off=Xopt(1,3);
cam_rol=Xopt(1,4);
cam(1,2)=Xopt(1,5);
cam(2,2)=Xopt(1,6);
cam_foc=Xopt(1,7);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% viewshed %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculate camera parameters (raster position and altitude)
[cam, cam_rc]=Cam_PRACTISE(cam, cam_off, demWrcz, headerv_W);
% Specified camera angles viewshed
if vs==true
    toc
%   Normalize DEM to viewpoint altitude
    demWvs=(demWrcz/headerv_W(5,1))-cam_rc(3,1);
    demWvs(cam_rc(2,1),cam_rc(1,1))=0;
%   Cam properties & viewshed raster initialization
    disp('Start calculation of specified camera angles (+ 0.1rad) viewshed')
%       field of view (use maximum range, if cam_rol~=0 + 0.1rad)
    roldummy=cam_rol*pi/180; % deg to rad
    FOV_vert=atan((0.5*cam_hei*abs(cos(roldummy))+0.5*cam_wid*abs(sin(roldummy))) ...
                    /cam_foc)+0.1;
    FOV_hor=atan((0.5*cam_wid*abs(cos(roldummy))+0.5*cam_hei*abs(sin(roldummy))) ...
                    /cam_foc)+0.1;
    clear roldummy
%       calculate azimuth values (clockwise from Azi(1,1) to Azi(2,1),
%       with N=0, E=pi/2, S=pi, W=3*pi/2)
%           delta DEM in the viewing direction [DEM pixel] -> Cam target-Cam position
    N0_rc(1,1)=cam_rc(1,2)-cam_rc(1,1); % delta longitude (=delta cols)
    N0_rc(2,1)=cam_rc(2,1)-cam_rc(2,2); % delta latitude (=-delta rows)        
    N0_rc(3,1)=cam_rc(3,2)-cam_rc(3,1); % delta altitude
%           calculate angles in viewing direction
    if N0_rc(2,1)>0 && N0_rc(1,1)>0 % delta latitude,delta longitude > 0 = NE
        N0_rc_hor=pi/2-atan(N0_rc(2,1)/N0_rc(1,1)); 
    elseif N0_rc(2,1)<0 && N0_rc(1,1)>0 % delta latitude < 0 & delta longitude > 0 = SE
        N0_rc_hor=pi/2+atan(N0_rc(2,1)/N0_rc(1,1)); 
    elseif N0_rc(2,1)<0 && N0_rc(1,1)<0 % delta latitude,delta longitude < 0 = SW
        N0_rc_hor=3*pi/2-atan(N0_rc(2,1)/N0_rc(1,1));
    elseif N0_rc(2,1)>0 && N0_rc(1,1)<0 % delta latitude > 0 & delta longitude < 0 = NW
        N0_rc_hor=3*pi/2+atan(N0_rc(2,1)/N0_rc(1,1)); 
    elseif N0_rc(2,1)==0 && N0_rc(1,1)==0
        error(strcat('Camera position and target at the same xy coordinate,', ... 
            ' observing the sky or the soil below?'));
    elseif N0_rc(2,1)==0
        if N0_rc(1,1)>0
            N0_rc_hor=pi/2;
        elseif N0_rc(1,1)<0
            N0_rc_hor=3*pi/2;
        end
    elseif N0_rc(1,1)==0
        if N0_rc(2,1)>0
            N0_rc_hor=0;
        elseif N0_rc(2,1)<0
            N0_rc_hor=pi;
        end
    end
%           calculate azimuth boundaries (with N=0, E=pi/2, ...)
    Azi(1,1)=N0_rc_hor-FOV_hor;
    if Azi(1,1)<0
        Azi(1,1)=2*pi+Azi(1,1);
    end
    Azi(2,1)=N0_rc_hor+FOV_hor;
    if Azi(2,1)>2*pi
        Azi(2,1)=Azi(2,1)-2*pi;
    end
%           define viewshed sectors with NE=1, SE=2, SW=3, NW=4
    for i=1:2
        if Azi(i,1)>=0 & Azi(i,1)<pi/2 %N-E
            Azi(i,2)=1;
        elseif Azi(i,1)>=pi/2 & Azi(i,1)<pi %E-S
            Azi(i,2)=2;
        elseif Azi(i,1)>=pi & Azi(i,1)<3*pi/2 %S-W
            Azi(i,2)=3;
        elseif Azi(i,1)>=3*pi/2 & Azi(i,1)<2*pi %W-N
            Azi(i,2)=4;
        end
    end
    if Azi(1,1)>Azi(2,1) % Azimuth angle including N (2*pi to 0)
        azi_sec=[Azi(2,2):4, 1:Azi(1,2)]
    else
        azi_sec=[Azi(1,2):Azi(2,2)]
    end
    disp('The above sector(s) with NE=1, SE=2, SW=3, NW=4 will be calculated');
%       vertical angle (from Vert(1,1) to Vert(2,1); upwards to downwards)
%           upwards=+, horizontal=0, downwards=-
    N0_rc_vert=atan(N0_rc(3,1)/norm(N0_rc(1:2,1))); % norm(N0_rc(1:2,1)) is never 0
    Vert(1,1)=N0_rc_vert+FOV_vert;
    Vert(2,1)=N0_rc_vert-FOV_vert;
    if Vert(1,1)>pi/2
        Vert(1,1)=pi/2;
        disp('sky observation angle adjusted to pi/2 as higher values are not usable,')
        disp('-> if wanted, use 2 different viewshed calculations and intersect them!')
    elseif Vert(2,1)<-pi/2
        Vert(2,1)=-pi/2;
        disp('surface observation angle adjusted to -pi/2 as lower values are not usable,')
        disp('-> if wanted, use 2 different viewshed calculations and intersect them!')
    end
    clear  FOV_hor FOV_vert  
%       viewshed raster initialization
    viewW=-1*ones(headerv_W(2,1),headerv_W(1,1));
    viewW(isnan(demWvs))=NaN;
% initialization of temporary rasters and variables
%   R=maximal normalized altitude of visibile pixel
%   index=number-indexed raster (row 1 & col 1 = 1, row 2 & col 1 = 2, ...) 
    R=viewW; 
    index=reshape([1:headerv_W(2,1)*headerv_W(1,1)],headerv_W(2,1),[]); %
    r=cam_rc(2,1);
    c=cam_rc(1,1);
    cam_rc_index=index(r,c);
% viewshed 
%   2 assumptions for specified camera angles computation: 
%       1) algorithm assumes that all pixels in the chosen sectors
%          that lie inbetween the vertical angles are visible
%       2) but azimuth angles are saved in case of vertical visibility and
%          taken out later if they do not fulfill the azimuth conditions
%   calculation of the viewpoint and the 8 neighbouring pixels 
    viewW(r-1,c+1)=pi/4; %NE
    viewW(r+1,c+1)=3/4*pi; %SE
    viewW(r+1,c-1)=pi+pi/4; %SW
    viewW(r-1,c-1)=pi+3/4*pi; %NW
    viewW(r-1,c)=0; %N
    viewW(r+1,c)=pi; %S
    viewW(r,c-1)=3*pi/2; %W
    viewW(r,c+1)=pi/2; %E
    for i=-1:1
        for j=-1:1
            if i~=0 | j~=0
                vertdummy=atan(demWvs(r+i,c+j)/sqrt(i^2+j^2));
                if vertdummy>Vert(1,1) | vertdummy<Vert(2,1)
                    viewW(r+i,c+j)=-1;
                end
            else
                if Vert(2,1)>-pi/2
                    viewW(r,c)=-1;
                else
                    viewW(r,c)=N0_rc_hor; 
                end
            end
        end
    end
    R(r-1:r+1,c-1:c+1)=demWvs(r-1:r+1,c-1:c+1);
%   calculation of the N, E, S and W directions, seen from the viewpoint
%   (using proxy vectors)
    for dir=1:4
        if dir==1 % N
%       flipud to S direction (=from viewpoint to the south)
            demdummy=flipud(demWvs(1:r-1,c));
            viewdummy=flipud(viewW(1:r-1,c));
            Rdummy=flipud(R(1:r-1,c));
            indexdummy=flipud(index(1:r-1,c));
            azidummy=0;
        elseif dir==2 % E
%       transpose to S direction
            demdummy=demWvs(r,c+1:headerv_W(1,1))'; 
            viewdummy=viewW(r,c+1:headerv_W(1,1))';
            Rdummy=R(r,c+1:headerv_W(1,1))';
            indexdummy=index(r,c+1:headerv_W(1,1))';
            azidummy=pi/2;
        elseif dir==3 % S
            demdummy=demWvs(r+1:headerv_W(2,1),c);
            viewdummy=viewW(r+1:headerv_W(2,1),c);
            Rdummy=R(r+1:headerv_W(2,1),c);
            indexdummy=index(r+1:headerv_W(2,1),c);
            azidummy=pi;
        elseif dir==4 % W
%       transpose and flipud to S direction
            demdummy=flipud(demWvs(r,1:c-1)');
            viewdummy=flipud(viewW(r,1:c-1)');
            Rdummy=flipud(R(r,1:c-1)');
            indexdummy=flipud(index(r,1:c-1)');
            azidummy=3*pi/2;
        end
%       S direction or rotated/mirrored other directions to S direction
        for l=2:size(demdummy,1)
            Z=Rdummy(l-1,1)*l/(l-1);
            if isfinite(demdummy(l,1)) 
                if demdummy(l,1)>Z
                    vertdummy=atan(demdummy(l)/l);
                    if vertdummy>Vert(1,1) | vertdummy<Vert(2,1)
                        viewdummy(l,1)=-1; % invisible
                    else
                        viewdummy(l,1)=azidummy; % vertically visible
                    end
                    Rdummy(l,1)=demdummy(l,1);
                else
                    viewdummy(l,1)=-1; % invisible
                    Rdummy(l,1)=Z;
                end
            else
                viewdummy(l,1)=NaN; % like invisible, but NaN
                Rdummy(l,1)=Z;
            end
        end
        viewW(indexdummy)=viewdummy;
        R(indexdummy)=Rdummy;
        clear demdummy viewdummy Rdummy indexdummy
    end
    clear l dir azidummy 
%   calculation of the chosen sectors (using proxy arrays)
%       (W-NW + NW-N = W-N, S-SW + SW-W = S-W,
%       N-NE + NE-E = N-E, E-SE + SE-S = E-S)
    for dir=azi_sec(1):azi_sec(end)
        if dir==1 % NE-E & N-NE
%       leftrright array (NE flipped to NW)
            demdummy=fliplr(demWvs(1:r,c:headerv_W(1,1)));
            viewdummy=fliplr(viewW(1:r,c:headerv_W(1,1)));
            Rdummy=fliplr(R(1:r,c:headerv_W(1,1)));
            indexdummy=fliplr(index(1:r,c:headerv_W(1,1)));
            secdummy=-pi/2; % angle clockwise from N=0
        elseif dir==2 % E-SE & SE-S
%       upsidedown + leftrright array (SE flipped to NW)
            demdummy=flipud(fliplr(demWvs(r:headerv_W(2,1),c:headerv_W(1,1))));
            viewdummy=flipud(fliplr(viewW(r:headerv_W(2,1),c:headerv_W(1,1))));
            Rdummy=flipud(fliplr(R(r:headerv_W(2,1),c:headerv_W(1,1))));
            indexdummy=flipud(fliplr(index(r:headerv_W(2,1),c:headerv_W(1,1))));
            secdummy=pi/2; 
        elseif dir==3 % W-SW & S-SW 
%       upsidedown array (SW flipped to NW)
            demdummy=flipud(demWvs(r:headerv_W(2,1),1:c));
            viewdummy=flipud(viewW(r:headerv_W(2,1),1:c));
            Rdummy=flipud(R(r:headerv_W(2,1),1:c));
            indexdummy=flipud(index(r:headerv_W(2,1),1:c));
            secdummy=-3*pi/2;
        elseif dir==4 % W-NW & NW-N
            demdummy=demWvs(1:r,1:c);
            viewdummy=viewW(1:r,1:c);
            Rdummy=R(1:r,1:c);
            indexdummy=index(1:r,1:c);
            secdummy=3*pi/2; 
        end
        rr=size(demdummy,1);
        cc=size(demdummy,2);
        m=rr-2;
        n=cc-1;
%       W-NW and NW-N = NW sector or rotated/mirrored other sectors to NW sector
        while ( n>0 ) % cols
            if ( m>0 ) % rows
                if ( rr-m<cc-n ) % "south of NW-diagonal"
                    Z=(rr-m)*(Rdummy(m,n+1)-Rdummy(m+1,n+1)) ...
                        +(cc-n)*((rr-m)*(Rdummy(m+1,n+1)-Rdummy(m,n+1))+Rdummy(m,n+1)) ...
                          /(cc-n-1);
                elseif ( rr-m>cc-n ) % "north of NW-diagonal"
                    Z=(cc-n)*(Rdummy(m+1,n)-Rdummy(m+1,n+1)) ...
                        +(rr-m)*((cc-n)*(Rdummy(m+1,n+1)-Rdummy(m+1,n))+Rdummy(m+1,n)) ...
                          /(rr-m-1);
                elseif ( rr-m==cc-n ) % "NW-diagonal"
                    Z=Rdummy(m+1,n+1)*(rr-m)/(rr-m-1);
                end
                if isfinite(demdummy(m,n))
                    if demdummy(m,n)>Z
                        azidummy=abs(atan((rr-m)/(cc-n))+secdummy);
                        vertdummy=atan(demdummy(m,n)/sqrt((rr-m)^2+(cc-n)^2));
                        if vertdummy>Vert(1,1) | vertdummy<Vert(2,1)
                            viewdummy(m,n)=-1; % invisible
                        else
                            viewdummy(m,n)=azidummy; % vertically visible
                        end
                        Rdummy(m,n)=demdummy(m,n);
                    else
                        viewdummy(m,n)=-1; % invisible
                        Rdummy(m,n)=Z;
                    end
                else
                    viewdummy(m,n)=NaN; % like invisible, but NaN
                    Rdummy(m,n)=Z;
                end
                m=m-1;
            else
                m=rr-1;
                n=n-1;
            end
        end
        clear m n
        viewW(indexdummy)=viewdummy;
        R(indexdummy)=Rdummy;
        clear demdummy viewdummy Rdummy indexdummy 
    end
    clear rr cc secdummy azidummy vertdummy Z dir
%   visibility decision: azimuth of vertically visible pixels
    if Azi(1,1)>Azi(2,1)
        viewW(viewW<=Azi(1,1))=1;
        viewW(viewW>=Azi(2,1))=1;
    else
        viewW(viewW<Azi(1,1))=0;
        viewW(viewW>Azi(2,1))=0;
        viewW(viewW>0)=1;
    end
    viewW(viewW<0)=0;
    viewW(isnan(viewW))=-9999;
end % if vs==true
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% projection %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
toc
disp('Starting projection')
% prepare DEM dependent on the viewshed
%   initialization 
demW=NaN(6,headerv_W(2,1)*headerv_W(1,1)); % if row, col included: 6
k=0;
%   DEM grid point coordinates (middle of pixel)
for i=1:headerv_W(1,1) % ncols
    for j=1:headerv_W(2,1) % nrows
        if viewW(j,i)==1 % visible
            k=k+1;
            demW(1,k)=headerv_W(3,1)-(0.5*headerv_W(5,1))+(i*headerv_W(5,1)); % longitude
            demW(2,k)=headerv_W(4,1)+(headerv_W(2,1)*headerv_W(5,1)) ...
                +(0.5*headerv_W(5,1))-(j*headerv_W(5,1)); % latitude
            demW(3,k)=demWrcz(j,i); %altitude
            demW(4,k)=1;
            demW(5,k)=j; %ascii row
            demW(6,k)=i; %ascii col
        end
    end
end
demW=demW(1:6,isfinite(demW(1,:))); % if row, col included: 1:6
clear i j k 
% prepare projection parameters
%   in case that optimization off
if os==false
    Xopt=X0;
end
% Project and scale the world coordinates (xyzW in m) to the photo plane 
% (crP in pixels)
%   1) translate world reference system W to reference system T 
%      -> set coordinate system origin = viewpoint 
%   2) rotate reference system T to camera reference system C 
%      -> set coordinate system axis in the directions right (U), up (V) 
%      and depth (N) of image
%   3) project and scale camera reference system C to perspective
%   projection system P
%      -> project to metric cam image using central projection 
%      (at depth of focal length)
%      -> scale metric to pixel based image size
%      -> translate to the origin of the 2D coordinate system of
%      the photograph
%      -> ceil to integer pixel values
%   4) extract RGB pixel values and classify
%      -> keep rows and columns of DEM (ASCII file [cf. ESRI ArcGIS]) 
%      -> keep pixel rows and columns of photograph
%      -> extract RGB values for pixel rows and columns
%      -> use a threshold (e.g. [R,G,B] > [150,150,150] = snow)
%      -> use delta approach to figure out bright stone
%      (e.g. [max(R,G,B)-min(R,G,B)] < 10 = snow )
%
%   call projection for steps 1) to 3)
%       for DEM
[demP]=Proj_PRACTISE(demW, Xopt, pix_c, pix_r, cam_wid, cam_hei, ...
    demWrcz, headerv_W);
%       for GCP
if gs==true
%           original    
    [gcpRMSE_orig, gcpP]=DDS_RMSE_PRACTISE(X0, demWrcz, headerv_W);
%           optimal 
%               in case of optimization
    if os==true
        [gcpRMSE_opt, gcpP]=DDS_RMSE_PRACTISE(Xopt, demWrcz, headerv_W);
        clear Fopt
        %
        plot(gcpP(1,:)+0.5, gcpP(2,:)+0.5, '.', 'MarkerEdgeColor', [1 0 0], 'MarkerSize', 10)
    end
end
%%%%% loop for multiple photograph classification %%%%%
for photoloop=1:N_images
    if is==true
%       define the respective photo file
        fin_image=fin_images(photoloop);
        photo = imread(strcat(fin_imfolder,char(fin_image))); 
%       define the respective output files
%           classification (ASCII-file [cf. ESRI ArcGIS])
        fin_imageDot=cell2mat(strfind(fin_image,'.'));
        fout_classW_asc=strcat(fin_image{1,1}(1:fin_imageDot-1),'_class.asc');
%           matlab (MAT-file)
        fout_mat=strcat(fin_image{1,1}(1:fin_imageDot-1),'_proj.mat');
    end
%
disp('Starting classification of the photograph:')
disp(fin_image)
%   initialize for 4)
rgbimage=NaN(8,size(demW,2));
%   4) extract RGB pixel values and classify
% extract pixel positions
i=find(demP(1,:)<pix_c & demP(1,:)>0 & demP(2,:)<pix_r & demP(2,:)>0);
rgbimage(1:2,i)=demW(5:6,i); % DEM: row col
rgbimage(3:4,i)=demP(:,i); % image: col row
% delete pixels (outside of the photograph boundaries and temporary data)
demW=demW(1:3,i);
demP=demP(1:2,i);
% extract RGB values
for j=1:length(i)
    rgbimage(5:7,i(j))=photo(rgbimage(4,i(j)),rgbimage(3,i(j)),1:3); % RGB values
end
% delete NaNs
i=find(isnan(rgbimage(1,:)));
rgbimage(:,i)=[];
clear i j
% classify
if cs==false
    disp('All visible pixels with RGB values above the threshold')
    disp(thres_rgb) % classification rule 1
    disp('and below the delta of max(RGB) to min(RGB)')
    disp(delta_rgb) % classification rule 2
    disp('are classified as snow.')
    %
    rgbimage(8,:)=0; % no snow
    % 2 classification rules
    i=find(rgbimage(5,:)>=thres_rgb(1,1) & rgbimage(6,:)>=thres_rgb(2,1) & ...
                    rgbimage(7,:)>=thres_rgb(3,1));
    j=find((max(rgbimage(5:7,i))-min(rgbimage(5:7,i)))<=delta_rgb);
    % 
    rgbimage(8,i(j))=1; % snow
elseif cs==true
    % RGB histogram investigation (Salvatori, 2011)
    %   initialization
    rgbimage(8,:)=0; % no snow
    %   blue band
    rgbhist=hist(rgbimage(7,:), 0:255);
    rgbhistmean=tsmovavg(rgbhist,'s',movavgwindow);
    rgbhistmean=rgbhistmean([(ceil(movavgwindow/2)):end,1:(ceil(movavgwindow/2))-1]);
    %   find local minimum
    rgbhistmeandiff=[diff(rgbhistmean),NaN];
    i=find(rgbhistmeandiff(thres_b_orig:max(rgbimage(7,:)))<0);
    i=i(1)+thres_b_orig-1;
    j=find(rgbhistmeandiff(i:max(rgbimage(7,:)))>0);
    if j<max(rgbimage(7,:))
        thres_b_image=j(1)+i-1-1;
        disp('The classification threshold in the blue band is set to')
        disp(thres_b_image)
    else
        thres_b_image=127;
        disp('No local minima in the smoothed DN histogram could be found beyond DN=126')
        disp('The classification threshold in the blue band is set to')
        disp(thres_b_image)
    end
    clear i j rgbhistmeandiff
    %   classification
    rgbimage(8,rgbimage(7,:)>=thres_b_image)=1; % snow
    %   plot DN frequency histogram
    figure2=figure(2);
    %
    set(figure2, 'Color',[1 1 1], 'InvertHardcopy','off', ...
        'PaperPosition',[0 0 50 35], 'PaperSize',[50 35], 'PaperType','<custom>');
    axes1 = axes('FontSize',15.76, 'FontWeight','demi', 'Layer','top',...
        'TickDir','out', 'XTick',[0 25 50 75 100 125 150 175 200 225 250],...
        'Parent',figure2);
    xlabel(axes1,'Digital Number','FontSize',18.011);
    ylabel(axes1,'Number of pixels','FontSize',18.011);
    box(axes1,'on');
    hold(axes1,'all');
    hist(rgbimage(7,:),0:255)
    figure2hist = findobj(gca,'Type','patch');
    set(figure2hist, 'EdgeColor', 'none', 'FaceColor', [0 0.749 0.749])
    axis1=axis;
    axis([0, 255, 0, axis1(4)])
    %
    plot(0:255,rgbhistmean, 'LineWidth', 3, 'Color', [0.0784 0.1686 0.5490])
    plot([thres_b_orig,thres_b_orig], [0,axis1(4)], 'r', 'LineWidth', 2)
    plot([thres_b_image,thres_b_image], [0,axis1(4)], 'g', 'LineWidth', 2)
    %
    clear figure2 axes1 figure2hist axis1 
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% plot %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
toc
disp('start plotting, saving and writing files')
% plot photograph, classification (projection) and GCP
figure3=figure(3);
%
set(figure3, 'Color',[1 1 1], 'InvertHardcopy','off', ...
  'PaperPosition',[0 0 50 35], 'PaperSize',[50 35], 'PaperType','<custom>');
axes1 = axes('FontSize',15.91, 'FontWeight','demi', 'Layer','top',...
  'TickDir','out', 'XAxisLocation','top', 'YDir','reverse', ...
  'DataAspectRatio', [1,1,1], 'Parent',figure3);
axis(axes1,[0.5 pix_c+0.5 0.5 pix_r+0.5]);
xlabel(axes1,'Pixel column','FontSize',18.185);
ylabel(axes1,'Pixel row','FontSize',18.185);
box(axes1,'on');
hold(axes1,'all');
%
image(photo, 'Parent',axes1);                                                    
plot(rgbimage(3,(rgbimage(8,:)==0)), rgbimage(4,(rgbimage(8,:)==0)), 'b.','MarkerSize',6)
plot(rgbimage(3,(rgbimage(8,:)==1)), rgbimage(4,(rgbimage(8,:)==1)), 'r.','MarkerSize',6)
%
clear figure3 axes1
%
% plot(demP(1,:)+0.5, demP(2,:)+0.5, 'r.','MarkerSize',2)
% plot(camP(1,1)+0.5, camP(2,1)+0.5,'k.', 'MarkerSize',6)
% plot(gcpP(1,:)+0.5, gcpP(2,:)+0.5, 'ro', 'MarkerSize',10)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% save %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~exist('fout_mat')
    error('Matlab output file name is not defined.')
end
save(strcat(fout_folder,fout_mat), 'cam', 'cam_rc', 'cam_off', 'cam_rol', 'cam_foc', ...
     'header_W', 'headerv_W', 'cam_hei', 'cam_wid', ...
     'pix_r', 'pix_c', 'fin_demW',  'demW', 'demP', ...
     'fin_image', 'fout_classW_asc', 'rgbimage', ...
     'vs', 'gs', 'os', 'cs', 'is', 'fout_folder');
if vs==false
    save(strcat(fout_folder,fout_mat), 'fin_viewW', '-append');
else
    save(strcat(fout_folder,fout_mat), 'fout_viewW_asc', 'Azi', 'Vert', '-append'); 
end
if gs==true
    save(strcat(fout_folder,fout_mat), 'N_gcpW', 'name_gcpW', 'gcpW', 'gcpP', ...
        'gcpRMSE_orig', 'fin_gcpW', '-append');
end
if cs==false
    save(strcat(fout_folder,fout_mat), 'thres_rgb', 'delta_rgb', '-append');
elseif cs==true
    save(strcat(fout_folder,fout_mat), 'thres_b_orig', 'movavgwindow', ...
        'thres_b_image', 'rgbhist', 'rgbhistmean', '-append');
end
if os==true
    save(strcat(fout_folder,fout_mat), 'X0', 'UB', 'LB', 'DDS_MaxEval', 'DDS_R', ...
        'Xopt', 'gcpRMSE_opt', 'cam_orig', 'cam_rc_orig', '-append');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% write %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% write GIS ASCII file: classification
classW=NaN(headerv_W(2,1),headerv_W(1,1));
classW(:,:)=-9999;
for i=1:size(rgbimage,2)
    classW(rgbimage(1,i),rgbimage(2,i))=rgbimage(8,i);
end
fid=fopen(strcat(fout_folder,fout_classW_asc), 'w');
if fid==-1
    print=('Classification output file could not be created.');
end
for i=1:6
    fprintf(fid, '%s\n', char(header_W(i,:)));
end
fclose(fid);
clear i fid 
dlmwrite(strcat(fout_folder,fout_classW_asc), classW, 'delimiter', ' ', '-append')
% write GIS ASCII file: viewshed
if vs==true
    viewW(classW(:,:)==-9999)=-9999;
    fid=fopen(strcat(fout_folder,fout_viewW_asc), 'w');
    if fid==-1
        print=('Viewshed output file could not be created.');
    end
    for i=1:6
        fprintf(fid, '%s\n', char(header_W(i,:)));
    end
    fclose(fid);
    clear i fid 
    dlmwrite(strcat(fout_folder,fout_viewW_asc), viewW, 'delimiter', ' ', '-append')
end
%
disp('Load workspace with the filename:')
disp(fout_mat) 
% write Matlab figure of
%   figure 1: DDS optimisation
if os==true
    saveas(figure(1),strcat(fout_folder,'DDS_opti'),'fig');
end
%   figure 2: DN frequency histogram
%   figure 3: classification
if is==false
    if cs==true
        saveas(figure(2),strcat(fout_folder,'DN_freq_hist'),'fig');
    end
    saveas(figure(3),strcat(fout_folder,'OverlayClass'),'fig');
elseif is==true
    if cs==false
        saveas(figure(3),strcat(fout_folder,fin_image{1,1}(1:fin_imageDot-1),...
            '_',num2str(thres_rgb(1,1)),'_',num2str(delta_rgb),'_manu'),'fig');
        clear fin_imageDot
    elseif cs==true
        saveas(figure(2),strcat(fout_folder,fin_image{1,1}(1:fin_imageDot-1),...
            '_DN_freq_hist'),'fig');
        saveas(figure(3),strcat(fout_folder,fin_image{1,1}(1:fin_imageDot-1),...
            '_',num2str(thres_b_image(1,1)),'_auto'),'fig');
        clear fin_imageDot
    end
end
toc
%%%%% loop for multiple photograph classification %%%%%
close all                                                                
end
%%%%%
% delete intermediate file for RMSE function
%   in case that GCP or DDS switch on
if os==true | gs==true
    delete('DDSintermediatedelete.mat')
end
%
disp('PRACTISE was executed successfully.')
clear all                                                                