% 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));


% 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.2;        %denominator in conduction coefficient
lambda = 1.95;   %time step

[rows,cols] = size(myImage);        %number of samples
diff = double(noisyImage);

for i = 1:max_it
    i % print iteration number
    
    % Construct diffl which is the same as diff but
    % has an extra padding of zeros around it.
    diffl = zeros(rows+2, cols+2);
    diffl(2:rows+1, 2:cols+1) = diff;
    
    % North, South, East and West differences
    deltaN = diffl(1:rows,2:cols+1)   - diff;
    deltaS = diffl(3:rows+2,2:cols+1) - diff;
    deltaE = diffl(2:rows+1,3:cols+2) - diff;
    deltaW = diffl(2:rows+1,1:cols)   - diff;
    
    cN = 1./(1 + (deltaN/K).^2);
    cS = 1./(1 + (deltaS/K).^2);
    cE = 1./(1 + (deltaE/K).^2);
    cW = 1./(1 + (deltaW/K).^2);
    
    diff = diff + lambda*(cN.*deltaN + cS.*deltaS + cE.*deltaE + cW.*deltaW);

    if mod(i,10) == 0
        %plot results
        
        figure;
        subplot(2,3,1);
        imshow(myImage);
        title('Actual Image');
        
        subplot(2,3,2);
        image(normImage);
        title('Normalized Image');
        
        subplot(2,3,3);
        image(noisyImage);
        title('Noisy Image');
        
        subplot(2,3,4);
        imshow(diff);
        %     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);
%         image(u_imp);
%         %     axis([0, n, -0.5, 1.5]);
%         title(['PM imp. -> iter:'  num2str(i) ' alpha:' num2str(alpha) ' K:' num2str(K) ' lambda:' num2str(lambda)]);
        
    end
    
    %     drawnow;
end


