% The following script produces a crevasse map of the image within this
% folder.
%
% B. Altena - 2018.

%% admin
clc;
fprintf('Crevasse detection in optical data by the use of the Radon-transform\n\n');
fprintf(['For application please cite:\n Gong et al. 2018.\n', ...
    'Simulating the roles of crevasse routing of surface water and\n', ...
    'basal friction on the surge evolution of Basin 3, Austfonna ice-cap.\n', ...
    'The Cryosphere.\n\n']);
fprintf('For the methodology please cite:\n Altena. 2018.\n');
fprintf('Observing change in glacier flow by using optical satellites\n');
fprintf('PhD thesis, University of Oslo.\n\n');

%% process images

% Radon - parameters
levels = 2;
alpha = 0.4;
blockSiz = 300; % meters

% Relaxation Labelling - parameters
options.iterMax = 5e2;
options.nomatch = false;
options.convFac = 1;        % convergence power of the iteration
options.mag = false;        % include magnitude in confidence function

cmapCirc = hsv(256);

imList = ['T35XMJ_20170926T114819_B02sub.tif'];

for i = 1:size(imList,1)
    imName = strtrim(imList(i,:));        
    fprintf(['reading imagery: ' imName '\n']);
    [ I, R, bbox ] = geotiffread(imName);
    gridSiz = round(blockSiz./10);

    % find no-data (a mask with water or clouds)
    if ~isempty(dir([imName(1:16) '*MSK*']))
        mskName = ls([imName(1:16) '*MSK*']);
        Imsk = imread(mskName);
    else
        figure();
        imagesc(I); axis equal; colormap(bone);
        title('draw polygon around water or clouds, if present');
        h = impoly();    M = createMask(h);
        imwrite(M, [imName(1:end-10) 'MSK.tif']);
        Imsk = M;
        clf; close all;
    end

    Inan = I==0;
    I = mat2gray(I);

    if exist('Imsk','var'), Inan = or(Imsk,Inan); end
    I(Inan) = NaN;

    clear Imsk mskName Inan mskName M R bbox
    %% Radon-transform
    fprintf('processing Radon Transform over imagery \n');
    fprintf('using slower but robust version \n');

    fprintf(['using square blocks of ' num2str(blockSiz) ' meter wide \n']); 

    warning off
    fun = @(block_struct) radonIce(block_struct.data,alpha,levels);
    Theta = blockproc(I, [gridSiz gridSiz], fun, 'UseParallel', true);
    warning on
    
    clear fun
    %% Relaxation labelling - get the most smooth transition
    primeDir = Theta(:,:,1);
    seconDir = Theta(:,:,3);

    dX = cat(3, -sind(2.*primeDir), -sind(2.*seconDir));
    dY = cat(3, cosd(2.*primeDir), cosd(2.*seconDir));
    dSNR = cat(3, Theta(:,:,2), Theta(:,:,4));
    
    fprintf('Applying Relaxation labelling,\n to get a smooth transition in estimate.\n');
    drho = relaxLabelling(dX,dY,dSNR, options);

    [val,idx] = max(drho,[],3);

    idx(idx==2) = 3;
    [M,N,B] = size(Theta);
    [Nm,Mm] = meshgrid(1:N,1:M);
    [f] = sub2ind([M,N,B],Mm(:),Nm(:),idx(:));
    primeDir = -reshape(Theta(f),M,N)+90;
    [f] = sub2ind([M,N,B],Mm(:),Nm(:),idx(:)+1);
    Score1 = reshape(Theta(f),M,N);
    idx(idx==3) = 2; idx(idx==1) = 3; idx(idx==2) = 3;
    [f] = sub2ind([M,N,B],Mm(:),Nm(:),idx(:));
    seconDir = -reshape(Theta(f),M,N)+90;
    [f] = sub2ind([M,N,B],Mm(:),Nm(:),idx(:)+1);
    Score2 = reshape(Theta(f),M,N);
    
    clear M N B f idx drho val dSNR dX dY Theta
    %% plot
    fprintf(['Estimated directions for ' imName(1:end-12) '\n']);

    figure();
    imagesc(I); axis equal; axis off; colormap(gray);
    title('Sentinel-2 image of Basin-3');
    figure();
    imagesc(primeDir); axis equal; axis off; colormap(cmapCirc);
    title('Dominant direction of crevasses [deg]');
    colorbar();
    figure();
    imagesc(Score1); axis equal; axis off; colormap(flip(pink,1));
    title('Strength of the crevasse signal');   
    
    clear imName
end
fprintf('done.\n');

clear cmapCirc levels alpha blockSiz options