% Experiment 5: Study of sharpening filters
% a. Laplacian for image sharpening  b. Unsharp masking (and highboost filtering)
clear; close all; clc;

I = im2double(imread('moon.tif'));      % blurred image of the North Pole of the moon

%% Part A: Using inbuilt MATLAB functions
h4 = fspecial('laplacian', 0);          % [0 1 0; 1 -4 1; 0 1 0]
lapA = imfilter(I, h4, 'replicate');
sharpA = I - lapA;                      % c = -1 because the centre coefficient is negative
usmA = imsharpen(I, 'Radius', 2, 'Amount', 1);

figure('Name', 'A_inbuilt');
subplot(1, 4, 1); imshow(I);          title('Original (moon.tif)');
subplot(1, 4, 2); imshow(lapA, []);   title('Laplacian (scaled)');
subplot(1, 4, 3); imshow(sharpA);     title('Sharpened: f - \nabla^2 f');
subplot(1, 4, 4); imshow(usmA);       title('imsharpen (unsharp)');

%% Part B: Custom implementation
% a. Laplacian with the 4-neighbour and 8-neighbour masks
w4 = [0 1 0; 1 -4 1; 0 1 0];
w8 = [1 1 1; 1 -8 1; 1 1 1];
lap4 = myFilter2(I, w4);
lap8 = myFilter2(I, w8);
sharp4 = min(max(I - lap4, 0), 1);
sharp8 = min(max(I - lap8, 0), 1);
fprintf('Max |custom Laplacian - imfilter| = %.2e\n', max(abs(lap4(:) - lapA(:))));

figure('Name', 'B_laplacian');
subplot(2, 3, 1); imshow(I);        title('Original');
subplot(2, 3, 2); imshow(lap4, []); title('Laplacian, 4-neighbour');
subplot(2, 3, 3); imshow(sharp4);   title('Sharpened, 4-neighbour');
subplot(2, 3, 5); imshow(lap8, []); title('Laplacian, 8-neighbour');
subplot(2, 3, 6); imshow(sharp8);   title('Sharpened, 8-neighbour');

% b. Unsharp masking: blur, subtract to get the mask, add k * mask back
gk = myGaussKernel(5, 1.5);
blurred = myFilter2(I, gk);
gmask = I - blurred;
usm = min(max(I + 1 * gmask, 0), 1);        % k = 1: unsharp masking
hb  = min(max(I + 3 * gmask, 0), 1);        % k > 1: highboost filtering

figure('Name', 'B_unsharp');
subplot(2, 3, 1); imshow(I);         title('Original');
subplot(2, 3, 2); imshow(blurred);   title('Gaussian blurred (5x5, \sigma = 1.5)');
subplot(2, 3, 3); imshow(gmask, []); title('Unsharp mask g_{mask}');
subplot(2, 3, 4); imshow(usm);       title('Unsharp masking (k = 1)');
subplot(2, 3, 5); imshow(hb);        title('Highboost (k = 3)');
subplot(2, 3, 6); plot(I(270, :), 'k'); hold on; plot(hb(270, :), 'r'); hold off;
xlim([1 size(I, 2)]); legend('original', 'highboost'); title('Intensity profile, row 270');

fprintf('Mean gradient magnitude: original %.4f, Laplacian(4) %.4f, Laplacian(8) %.4f, unsharp %.4f, highboost %.4f\n', ...
    meanGrad(I), meanGrad(sharp4), meanGrad(sharp8), meanGrad(usm), meanGrad(hb));

%% Local functions
function out = myFilter2(img, w)
    % correlation of img with the mask w (replicate padding)
    [m, n] = size(w);
    a = (m - 1) / 2; b = (n - 1) / 2;
    P = padarray(img, [a b], 'replicate');
    out = zeros(size(img));
    for x = 1:size(img, 1)
        for y = 1:size(img, 2)
            region = P(x:x+m-1, y:y+n-1);
            out(x, y) = sum(region .* w, 'all');
        end
    end
end

function g = myGaussKernel(sz, sigma)
    r = (sz - 1) / 2;
    [X, Y] = meshgrid(-r:r, -r:r);
    g = exp(-(X.^2 + Y.^2) / (2 * sigma^2));
    g = g / sum(g(:));
end

function v = meanGrad(img)
    % average absolute first difference, used as a simple sharpness measure
    dx = diff(img, 1, 2); dy = diff(img, 1, 1);
    v = mean(abs(dx(:))) + mean(abs(dy(:)));
end