% Experiment 2: Study of different image transformation techniques
% a. Image type conversion  b. Image resize  c. Image rotation  d. Image negative
clear; close all; clc;

rgbImg = imread('peppers.png');          % RGB image is read only once
[rows, cols, ~] = size(rgbImg);

%% Part A: Using inbuilt MATLAB functions (as given in the lab handout)
% a. Image type conversions
redCh   = rgbImg(:, :, 1);
greenCh = rgbImg(:, :, 2);
blueCh  = rgbImg(:, :, 3);
grayImg = rgb2gray(rgbImg);
binImg  = imbinarize(grayImg);                 % Otsu threshold by default
hsvImg  = rgb2hsv(rgbImg);
rgbFromHsv = hsv2rgb(hsvImg);
ycbcrImg = rgb2ycbcr(rgbImg);
[indImg, cmap] = rgb2ind(rgbImg, 256);
rgbFromInd = ind2rgb(indImg, cmap);
grayDouble = im2double(grayImg);               % uint8  -> double [0, 1]
grayUint8  = im2uint8(grayDouble);             % double -> uint8
grayUint16 = im2uint16(grayImg);               % uint8  -> uint16
fprintf('Class conversion: %s -> %s -> %s, uint16 max = %d\n', ...
    class(grayDouble), class(grayUint8), class(grayUint16), max(grayUint16(:)));

figure('Name', 'A_type_conversion');
subplot(3, 4, 1);  imshow(rgbImg);            title('Original RGB');
subplot(3, 4, 2);  imshow(redCh);             title('R channel');
subplot(3, 4, 3);  imshow(greenCh);           title('G channel');
subplot(3, 4, 4);  imshow(blueCh);            title('B channel');
subplot(3, 4, 5);  imshow(grayImg);           title('rgb2gray');
subplot(3, 4, 6);  imshow(binImg);            title('imbinarize');
subplot(3, 4, 7);  imshow(hsvImg(:, :, 1));   title('Hue (rgb2hsv)');
subplot(3, 4, 8);  imshow(rgbFromHsv);        title('hsv2rgb');
subplot(3, 4, 9);  imshow(ycbcrImg(:, :, 1)); title('Y (rgb2ycbcr)');
subplot(3, 4, 10); imshow(indImg, cmap);      title('Indexed (rgb2ind)');
subplot(3, 4, 11); imshow(rgbFromInd);        title('ind2rgb');
subplot(3, 4, 12); imshow(grayUint16);        title('im2uint16');

% b. Image resize and c. Image rotation
halfImg   = imresize(rgbImg, 0.5);
squareImg = imresize(rgbImg, [256 256]);
rot45     = imrotate(rgbImg, 45);              % counter-clockwise
rotCW30   = imrotate(rgbImg, -30);             % clockwise
fprintf('Resize: %dx%d -> %dx%d (0.5) and %dx%d ([256 256])\n', rows, cols, ...
    size(halfImg, 1), size(halfImg, 2), size(squareImg, 1), size(squareImg, 2));

figure('Name', 'A_resize_rotate');
subplot(2, 3, 1); imshow(rgbImg);    title(sprintf('Original %dx%d', rows, cols));
subplot(2, 3, 2); imshow(halfImg);   title(sprintf('imresize 0.5 (%dx%d)', size(halfImg, 1), size(halfImg, 2)));
subplot(2, 3, 3); imshow(squareImg); title('imresize [256 256]');
subplot(2, 3, 4); imshow(rot45);     title('imrotate 45 (CCW)');
subplot(2, 3, 5); imshow(rotCW30);   title('imrotate -30 (CW)');

% d. Image negative
negGray   = imcomplement(grayImg);
negRgb    = 255 - rgbImg;
negDouble = 1 - grayDouble;

figure('Name', 'A_negative');
subplot(2, 2, 1); imshow(grayImg);   title('Gray image');
subplot(2, 2, 2); imshow(negGray);   title('imcomplement');
subplot(2, 2, 3); imshow(negRgb);    title('255 - RGB');
subplot(2, 2, 4); imshow(negDouble); title('1 - double image');

%% Part B: Custom implementation (without the inbuilt functions)
myGray   = myRgb2Gray(rgbImg);
myBin    = myGray > mean(myGray(:));           % global threshold at the mean
myHalf   = myResize(rgbImg, round(rows / 2), round(cols / 2));
mySquare = myResize(rgbImg, 256, 256);
myRot45  = myRotate(rgbImg, 45);
myRotCW  = myRotate(rgbImg, -30);
myNeg    = myNegative(rgbImg);

fprintf('Max |myRgb2Gray - rgb2gray|  = %d gray levels\n', max(abs(double(myGray(:)) - double(grayImg(:)))));
fprintf('Max |myNegative - imcomplement| = %d\n', max(abs(double(myNeg(:)) - double(imcomplement(rgbImg(:))))));
fprintf('Custom resize size: %dx%d, custom rotation canvas (45 deg): %dx%d\n', ...
    size(myHalf, 1), size(myHalf, 2), size(myRot45, 1), size(myRot45, 2));

figure('Name', 'B_custom');
subplot(2, 4, 1); imshow(myGray);   title('Custom grayscale');
subplot(2, 4, 2); imshow(myBin);    title('Custom binary (T = mean)');
subplot(2, 4, 3); imshow(myHalf);   title('Custom resize 0.5');
subplot(2, 4, 4); imshow(mySquare); title('Custom resize 256x256');
subplot(2, 4, 5); imshow(myRot45);  title('Custom rotate 45 (CCW)');
subplot(2, 4, 6); imshow(myRotCW);  title('Custom rotate -30 (CW)');
subplot(2, 4, 7); imshow(myNeg);    title('Custom negative');
subplot(2, 4, 8); imshow(myNegative(myGray)); title('Custom gray negative');

%% Local functions
function g = myRgb2Gray(rgb)
    % luminance weights used by the ITU-R BT.601 standard
    rgb = double(rgb);
    g = uint8(0.2989 * rgb(:, :, 1) + 0.5870 * rgb(:, :, 2) + 0.1140 * rgb(:, :, 3));
end

function out = myResize(img, newR, newC)
    % nearest-neighbour interpolation: map each output pixel back to the input
    [r, c, ~] = size(img);
    rowIdx = min(r, max(1, round(((1:newR) - 0.5) * r / newR + 0.5)));
    colIdx = min(c, max(1, round(((1:newC) - 0.5) * c / newC + 0.5)));
    out = img(rowIdx, colIdx, :);
end

function out = myRotate(img, angleDeg)
    % rotation about the centre (positive = counter-clockwise) using inverse mapping
    [r, c, ch] = size(img);
    t = deg2rad(angleDeg);
    newR = ceil(abs(r * cos(t)) + abs(c * sin(t)));
    newC = ceil(abs(c * cos(t)) + abs(r * sin(t)));
    [X, Y] = meshgrid(1:newC, 1:newR);
    Xc = X - (newC + 1) / 2;
    Yc = Y - (newR + 1) / 2;
    xs = round( Xc * cos(t) - Yc * sin(t) + (c + 1) / 2);
    ys = round( Xc * sin(t) + Yc * cos(t) + (r + 1) / 2);
    valid = xs >= 1 & xs <= c & ys >= 1 & ys <= r;
    out = zeros(newR, newC, ch, 'like', img);
    for k = 1:ch
        src = img(:, :, k);
        dst = zeros(newR, newC, 'like', img);
        dst(valid) = src(sub2ind([r c], ys(valid), xs(valid)));
        out(:, :, k) = dst;
    end
end

function out = myNegative(img)
    % s = (L - 1) - r, with L = 256 for uint8 images
    L = double(intmax(class(img))) + 1;
    out = cast((L - 1) - double(img), class(img));
end