%============================================================================
% Copyright (C) 2014, Heikki Hyyti
%
% Modified (06/2014) by Jakke Kulovesi to implement the method from:
% http://www.reliefshading.com/analytical/shading_methods.html
%
% Permission is hereby granted, free of charge, to any person obtaining a
% copy of this software and associated documentation files (the "Software"),
% to deal in the Software without restriction, including without limitation
% the rights to use, copy, modify, merge, publish, distribute, sublicense,
% and/or sell copies of the Software, and to permit persons to whom the
% Software is furnished to do so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
% FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
% DEALINGS IN THE SOFTWARE.
%============================================================================
%
% generateHillShadeMap(inputMapFilename, inputXYZMatFilename, outputFilename)
%
% Description:
%    Converts map and depthmap and generates hillshaded output map image.
%    Needs wgs2utm by Alexandre Schimel, which can be downloaded from there:
%    http://www.mathworks.com/matlabcentral/fileexchange/14804-wgs2utm--version-2-
%
% Input:
%    inputMapFilename: Filename of input map image (got using fetch_map python script).
%    inputXYZMatFilename: Filename of XYZ.mat file generated using
%       loadXYZtoMat.m function.
%    outputFilename (optional): Filename to save resulting map image.
%    slopeMinMax (optional): To enable aspect shading of deep slopes. The value 
%       is a vector of slope minimum value in radians (default: 0) and Slope 
%       maximum value in radians (default: pi/4). If only aspect shading 
%       (without diffuse shading is needed, use minimum and maximum value of 0.
%
function generateHillShadeMap(inputMapFilename, inputXYZMatFilename, outputFilename, slopeMinMax)
    % compute name for map file from image filename (change ending)
    [pieces, delims] = regexp(inputMapFilename, '\.', 'split', 'match');
    j = [pieces(1:end-1); [delims(1:end-1) {''}]];
    calibrationFilename = [[j{:}] '.map'];
    
    % test input filenames
    if (~exist(inputMapFilename,'file'))
        disp(['Input Map image (' inputMapFilename ') not found!']);
        return;
    elseif (~exist(calibrationFilename,'file'))
        disp(['Input Map calibration file (' calibrationFilename ') not found!']);
        return;
    elseif (~exist(inputXYZMatFilename,'file'))
        disp(['Input depth data file (' inputXYZMatFilename ') not found!']);
        return;
    end
    
    % Load map from image file (got using fetch_map python script)
    map = imread(inputMapFilename);
    N_N = size(map,1); 
    N_E = size(map,2); 
    aspectShading = false;
    
    if (exist('slopeMinMax','var') && length(slopeMinMax) > 1)
        slopeMin = slopeMinMax(1);
        slopeMax = slopeMinMax(2);
        aspectShading = true;
    end
    
    % values from calibration map-file (got using fetch_map python script)
    calibration = loadMAPFile(calibrationFilename);

    Lon = calibration.MMPLL(:,1);
    Lat = calibration.MMPLL(:,2);
    [MMP_E,MMP_N] = wgs2utm(Lat,Lon,calibration.utmzone,'N'); %todo: N / S

    maxE = max(MMP_E);
    minE = min(MMP_E);
    maxN = max(MMP_N);
    minN = min(MMP_N);

    % compute depth map from mat-file generated using loadXYZtoMat.m
    data = load(inputXYZMatFilename);
    dxy = 10; %native resolution of xyz data

    N_x = ceil( (maxE-minE) / dxy );
    N_y = ceil( (maxN-minN) / dxy );

    depthMap = nan(N_y,N_x);

    %filter out z=-999 and x and y indexes which are out of area
    accepted = data.z > -999 & data.x >= minE & data.x <= maxE & data.y >= minN & data.y <= maxN;

    x = data.x(accepted);
    y = (minN+maxN) - data.y(accepted); %reverse y (same as map image)
    z = data.z(accepted);

    x_idx = 1 + floor((x - minE) / dxy);
    y_idx = 1 + floor((y - minN) / dxy);

    for i = 1:size(x,1)
        depthMap(y_idx(i), x_idx(i)) = z(i);
    end

    %remove some variables to reduce memory usage
    clear x y z x_idx y_idx accepted;

    % compute hillshade (exact diffuseShade implemented by Jakke Kulovesi)
    delta_x = [depthMap(:,2:end) - depthMap(:,1:end-1) zeros(size(depthMap,1),1)] ./ dxy;
    delta_y = [depthMap(2:end,:) - depthMap(1:end-1,:); zeros(1, size(depthMap,2))] ./ dxy;
    %hillshade = round(128 + 128*(delta_y+delta_x)); %(dummy approximation using gradients)
    diffuseShade = (1/sqrt(3)).*(delta_x + delta_y + 1)./sqrt(delta_x.^2 + delta_y.^2 + 1);
    diffuseShade = diffuseShade - (1/sqrt(3)); % reduce mean from diffuseShade
    
    if (aspectShading)
        % compute hillshade (aspectShade implemented by Jakke Kulovesi)
        lightingAzimuth = 1/4*pi;	
        steepness = abs(atan(delta_x + delta_y));    
        aspectShade = 0.5 .* (cos(atan2(delta_x,delta_y) - lightingAzimuth) + 1);
        aspectShade = 0.5*(aspectShade - 0.5); %reduce mean from aspectShade and reduce gain

        mask = (steepness-slopeMin) / slopeMax;
        mask(mask < 0) = 0;
        mask(mask > 1) = 1;

        % combine the two different shading methods
        diffuseShade = diffuseShade .* (1-mask) + aspectShade .* mask;
    end
    
    %compute hillshade image from shading(s)
    hillshade = round(128 + 256*diffuseShade); 
     
    nanidx = isnan(hillshade);
    hillshade(nanidx) = 128;
    
    %remove some variables to reduce memory usage
    clear diffuseShade aspectShade steepness mask nanidx;
    
    % resize hillshade to same resolution as map
    hillshade = imresize(uint8(hillshade), [N_N N_E]);

    % integrate hillshade to map image
    gain = 1.2;
    for i = 1:3
        map(:,:,i) = uint8(1.5*((double(map(:,:,i)).*(0.3342+double(hillshade)/383)).^gain)./(255^(gain-1)));
    end

    % save image to file
    if (exist('outputFilename','var') && ~isempty(outputFilename))
        imwrite(map, outputFilename);
    end

    % plot result
    figure(1); clf;
    imshow(map);
    axis equal;
end