% ENSC 180 - Introduction to Engineering Analysis Tools
% Solution to Assignment #5 - Basic Image Processing

% ----- Task 1 ----- 
X = imread('baboon.tiff');
X = double(X);
Y = wextend('2D', 'sym', X, [1,1]);
[r,c,p] = size(Y);
for k = 1:p
    for i = 2:r-1
        for j = 2:c-1
            Y(i,j,k) = (Y(i-1,j-1,k) + Y(i-1,j,k) + Y(i-1,j+1,k) + ...
                        Y(i,j-1,k)   + Y(i,j,k)   + Y(i,j+1,k)   + ...
                        Y(i+1,j-1,k) + Y(i+1,j,k) + Y(i+1,j+1,k))/9;
        end
    end
end

X1 = Y(2:r-1, 2:c-1, :);
figure(1); imshow(uint8(X));
figure(2); imshow(uint8(X1));


% ----- Task 2 ----- 
Y = wextend('2D', 'sym', X, [2,2]); 
[r,c,p] = size(Y);
for k = 1:p
    for i = 3:r-2
        for j = 3:c-2
            Y(i,j,k) = (Y(i-2,j-2,k) + Y(i-2,j-1,k) + Y(i-2,j,k) + Y(i-2,j+1,k) + Y(i-2,j+2,k) + ...
                        Y(i-1,j-2,k) + Y(i-1,j-1,k) + Y(i-1,j,k) + Y(i-1,j+1,k) + Y(i-1,j+2,k) + ...
                        Y(i,  j-2,k) + Y(i,  j-1,k) + Y(i,  j,k) + Y(i,  j+1,k) + Y(i,  j+2,k) + ...
                        Y(i+1,j-2,k) + Y(i+1,j-1,k) + Y(i+1,j,k) + Y(i+1,j+1,k) + Y(i+1,j+2,k) + ...
                        Y(i+2,j-2,k) + Y(i+2,j-1,k) + Y(i+2,j,k) + Y(i+2,j+1,k) + Y(i+2,j+2,k))/25;
        end
    end
end

X2 = Y(3:r-2, 3:c-2, :);
figure(3); imshow(uint8(X2));


% ----- Task 3 ----- 
% 3x3 averaging kernel
H = fspecial('average',[3 3]);
X3_1 = imfilter(X,H,'symmetric','same');
figure(4); imshow(uint8(X3_1));
%max(max(max(abs(X3_1-X1)))) % print the maximum difference
% 5x5 averaging kernel
H = fspecial('average',[5 5]);
X3_2 = imfilter(X,H,'symmetric','same');
figure(5); imshow(uint8(X3_2));
%max(max(max(abs(X3_2-X2)))) % print the maximum difference
% 5x5 Gaussian kernel with sigma = 1
H = fspecial('gaussian',[5 5], 1);
X3_3 = imfilter(X,H,'symmetric','same');
figure(6); imshow(uint8(X3_3));
%max(max(max(abs(X3_3-X2)))) % print the maximum difference


% ----- Task 4 -----
c = 0.57; sigma = 2;
H = fspecial('gaussian',[5 5], sigma);
Xg = imfilter(X2,H,'symmetric','same');
X4 = X2*c/(2*c-1) - Xg*(1-c)/(2*c-1);
figure(7); imshow(uint8(X4));
