% clear workspace and command line
clear;clc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% START OF PART3: Perona-Malik diffusion process (implicit form)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% read file
fid = fopen('myDataset.rawb','r');
A = fread(fid);

% reshape the dataset
B = reshape(A,[181 217 181]);

% extract the slice at z=90
myImage = B(:,:,90);


% normalize image to 0-1 interval
maxValue = max(myImage(:));
minValue = min(myImage(:));
myRange = maxValue - minValue;
normImage = ((myImage - minValue)./(maxValue-minValue)).*1;

% then add noise from a normal distribution of N(0,0.005)
noisyImage = normImage;
for i = 1:size(noisyImage,1)
    for j = 1:size(noisyImage,2)
        noisyImage(i,j) = normImage(i,j) + normrnd(0,0.005);
    end
end

%set algorithm parameters
max_it = 100;
alpha = 10;       %exponent in conduction coefficient
K = 3.5;        %denominator in conduction coefficient
lambda = 0.25;   %time step

n=numel(noisyImage);
% u0 = noisyImage(:);
u0 = noisyImage(:);
I = speye(n,n);
e = ones(n,1);

% Dfx computes forward differences in the x-direction
Dfx = spdiags([-e e],[0:1], n, n);
% Dbx computes backward differences in the x-direction
Dbx = spdiags([-e e],[-1:0], n, n);
% Dux computes upward differences in the y-direction
Dux = spdiags([-e e],[0,size(noisyImage,1)], n, n);
% Ddx computes downward differences in the y-direction
Ddx = spdiags([-e e],[-size(noisyImage,1),0], n, n);

% initialize the initial images of explicit version of
% Perona-Malik diffusion process
u_imp = u0;

for i = 1:max_it
    i % print iteration number
    
    % Compute conduction coefficients and place on the main diagonal of
    % a sparse matrix. These matrices can be used to scale a vector of
    % derivatives
    % f, b, u, d are for foward, backward, up and down
    Cfx = spdiags(1 ./ (1 + (abs(Dfx*u_) ./K).^(alpha)), 0, n, n);
    Cbx = spdiags(1 ./ (1 + (abs(Dbx*u_exp) ./K).^(alpha)), 0, n, n);
    Cux = spdiags(1 ./ (1 + (abs(Dux*u_exp) ./K).^(alpha)), 0, n, n);
    Cdx = spdiags(1 ./ (1 + (abs(Ddx*u_exp) ./K).^(alpha)), 0, n, n);
    
    % iterate using implicit formulation
    B = I - lambda*(Cfx*Dfx + Cbx*Dbx + Cux*Dux + Cdx*Cdx);
    u_imp = B\u_imp;
    if mod(i,10) == 0
        %plot all images and save every 10th image
        imagesc(reshape(u_imp,[181 217]));colormap(gray);
        title(['PeronaMalik Implicit -> Iteration:'  num2str(i) ' alpha:' num2str(alpha) ' K:' num2str(K) ' lambda:' num2str(lambda)]);
        drawnow;
    end
    
end