% clear workspace and command line
clear;clc;

% 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(max(myImage));
minValue = min(min(myImage));
myRange = maxValue - minValue;
normImage = ((myImage - minValue)./(maxValue-minValue)).*1;


% PART2: Perona-Malik diffusion process
% first add noice from 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 = 2;       %exponent in conduction coefficient
K = 0.02;        %denominator in conduction coefficient
lambda = 0.025;   %time step

% n = size(A,1);
n=numel(noisyImage);
% u0 = double(A);
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);


u_exp = u0;
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_exp) ./K).^(1+alpha)), 0, n, n);
    Cbx = spdiags(1 ./ (1 + (abs(Dbx*u_exp) ./K).^(1+alpha)), 0, n, n);
    Cux = spdiags(1 ./ (1 + (abs(Dux*u_exp) ./K).^(1+alpha)), 0, n, n);
    Cdx = spdiags(1 ./ (1 + (abs(Ddx*u_exp) ./K).^(1+alpha)), 0, n, n);
    
    % iterate using explicit formulation
    %
    A = I + lambda*(Cfx*Dfx + Cbx*Dbx + Cux*Dux + Cdx+Ddx);
    u_exp = A*u_exp;
    
    % iterate using implicit formulation (derivatives applied to u(t+1))
    %
    B = I - lambda*(Cfx*Dfx + Cbx*Dbx + Cux*Dux + Cdx*Cdx);
    u_imp = B\u_imp;
    
        %plot results
        
        figure;
        subplot(2,3,1);
        image(myImage);
        title('Actual Image');
        
        subplot(2,3,2);
        imshow(normImage);
        title('Normalized Image');
        
        subplot(2,3,3);
        imshow(noisyImage);
        title('Noisy Image');
        
        subplot(2,3,4);
        imshow(reshape(u_exp,[181 217]));
        %     axis([0, n, -0.5, 1.5]);
        title(['PM exp. -> iter:'  num2str(i) ' alpha:' num2str(alpha) ' K:' num2str(K) ' lambda:' num2str(lambda)]);
        
        subplot(2,3,5);
        imshow(reshape(u_imp,[181 217]));
        %     axis([0, n, -0.5, 1.5]);
        title(['PM imp. -> iter:'  num2str(i) ' alpha:' num2str(alpha) ' K:' num2str(K) ' lambda:' num2str(lambda)]);
        
    
    %     drawnow;
end