## Author: Hakkel Tamás
## Created: 2016-10-16

% Parameters:
%   image_matrix: uint8 type matrix of image
%   iteration_number: number of iterations during the process of anisotropic diffusion
%   K: constant for noise level
%   lamba: parameter for the strength of anisotropic diffusion

% This function uses 'exponential' form of g function that privileges sharp edges

%Create a function that:
%  ● realizes anisotropic diffusion,
%  ● parameters of this function should be as follows:
%    ○ input1: image matrix,
%    ○ input2: iteration number,
%    ○ input3: K,
%    ○ input4: lambda,
%    ○ output: processed image,
%  ● it should work with grayscale images (check if this holds),
%  ● although formulas (7), (8), (10) and the two different versions of function g (at the
%  bottom right of page 633) are enough to the implementation, I would like to
%  recommend to read the first 5 sections of the article.

% Reference: The original article (P. Perona, J Malik (July 1990). "Scale-space and edge detection using anisotropic diffusion".
% IEEE Tr. PAMI, 12 (7): 629–639.) can be found here: http://image.diku.dk/imagecanon/material/PeronaMalik1990.pdf


function [out] = anisotropic_exp (image_matrix, iteration_number, K, lambda)

if size(image_matrix,3) ~= 1
  error 'This picture is not grayscale!'
end

% Későbbi nearest-neighbour differences-beli circshift miatt, 
% hogy a széleken ne legyen gond. Ezt a végén úgyis levágom.
image_matrix = double(padarray(image_matrix, [1 1], "replicate"));
for i = 1:iteration_number

  % Olvasmány (8)-as képlet (nearest neighbour differences)
  imDiff_n = circshift(image_matrix, [-1 0]) - image_matrix;
  imDiff_w = circshift(image_matrix, [0 +1]) - image_matrix;
  imDiff_s = circshift(image_matrix, [+1 0]) - image_matrix;
  imDiff_e = circshift(image_matrix, [0 -1]) - image_matrix;
  
  % Olvasmány (10)-es képlet és az exponenciális g függvény
  [imGradMagn_n imGradDir ] = imgradient(uint8(imDiff_n));
  [imGradMagn_w imGradDir ] = imgradient(uint8(imDiff_w));
  [imGradMagn_s imGradDir ] = imgradient(uint8(imDiff_s));
  [imGradMagn_e imGradDir ] = imgradient(uint8(imDiff_e));
  c_n = exp(-((imGradMagn_n / K) .^ 2));
  c_w = exp(-((imGradMagn_w / K) .^ 2));
  c_s = exp(-((imGradMagn_s / K) .^ 2));
  c_e = exp(-((imGradMagn_e / K) .^ 2));
  
  % Olvasmány (7)-es képlet
  image_matrix = image_matrix + lambda .* ...
    (c_n .* imDiff_n + c_w .* imDiff_w + c_s .* imDiff_s + c_e .* imDiff_e);
    
end

out = uint8(image_matrix(2:end-2,2:end-2));

endfunction