clear, clc

%% initialize

% load and distort sample image (scale to 0..1 range)
I = double(imread('cameraman.tif'))/255;
[height, width] = size(I);

% parameters
sigma = 0.05;           % variance of noise
lambda = 1.0;           % 'smoothness' parameter of optimization
numit = 120;            % number of interations in the optimization
display_every_nth = 5;  % displaying intermediate results for debugging
init_zero = 0;          % initialize with zero or with the noisy image

% distort with gaussian noise
f = I + randn(size(I))*sigma;

% initialize with the distorted image or with zero (result should be equal)
if init_zero
    % initialize with black image
    u = zeros(height, width);
else
    u = f;
end

% % symmetric neighbourhood
N = [1  0; -1  0; 0  1; 0 -1];

%% run Gauss-Seidel
 
for i = 1:numit
    
    u_prev = u; % u_prev = u^(k), u = u^(k+1)
    
    for y = 1:height
        for x = 1:width
            val = f(y,x);
            count = 0;
            for j = 1:size(N,1)
                xp = x + N(j,1);
                yp = y + N(j,2);
                % ignore values outside the image bounds
                if (xp >= 1) && (xp <= width) && (yp >= 1) && (yp <= height)
                    val = val + lambda * u(yp, xp);
                    count = count +1;
                end
            end
            u(y,x) = val / (1 + lambda*count);
        end
    end

    disp(['restored (iteration ',num2str(i),...
        ', change: ',num2str(norm(u(:)-u_prev(:)))])
    
    % display images
    if mod(i, display_every_nth) == 0
        
        subplot(1,3,1)
        imshow(I)
        title("original")

        subplot(1,3,2)
        imshow(f)
        title("noisy")

        subplot(1,3,3)
        imshow(u)
        title("restored")
        
        drawnow
    end
end