1; % dummy code to have local functions
% for MATLAB shift the code from bottom to here and add 'Antialiasing', false in imresize
% for explaination see also https://en.wikipedia.org/wiki/Bicubic_interpolation#Bicubic_convolution_algorithm

%%%%%%%%%%%%%%%%% helpers %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function w = cubic01(d, a)
    % for |d| <= 1
    absd = abs(d);
    w = (a+2) * absd.^3 - (a+3) * absd.^2 + 1;
end

function w = cubic12(d, a)
    % for 1 < |d| <= 2
    absd = abs(d);
    w = a * absd.^3 - 5*a * absd.^2 + 8*a * absd - 4*a;
end

function [idx_rows, idx_cols] = get_sample_grid(image, scale)
    out_size = ceil(size(image) * scale);
    offsets = 1 / scale / 2;
    idx_rows = 0.5 + offsets + (0:out_size(1)-1) / scale;
    idx_cols = 0.5 + offsets + (0:out_size(2)-1) / scale;
end

%%%%%%%%%%%%%%%%% scalar implementation %%%%%%%%%%%%%%%%%%%%%%%%%
function p = intpolcub_single(i, d, a)
    if nargin == 2
        a = -0.5;
    end
    p = i(1)*cubic12(-d-1, a) + i(2)*cubic01(-d, a) + i(3)*cubic01(-d+1, a) + i(4)*cubic12(-d+2, a);
end

function out = bicubic_conv_loop(img, scale, padding)
    if nargin == 2
        padding = 'symmetric';
    end
    
    % calculate the new pixel indices in terms of the old pixel indices
    [idx_rows, idx_cols] = get_sample_grid(img, scale);

    % padding
    padded_img = padarray(img, [2, 2], padding);
    padded_idx_rows = idx_rows + 2;
    padded_idx_cols = idx_cols + 2;
    
    % interpolate each point
    out = zeros(length(idx_rows), length(idx_cols));
    for out_idx_row = 1:numel(padded_idx_rows)
        ir = padded_idx_rows(out_idx_row);
        fir = floor(ir);
        dy = ir - fir;
        for out_idx_col = 1:numel(padded_idx_cols)
            ic = padded_idx_cols(out_idx_col);
            fic = floor(ic);
            dx = ic - fic;
            
            % 4 x 4 neighborhood
            im_part = padded_img(fir-1:fir+2, fic-1:fic+2);
            
            % interpolate 4x in y-direction and then the interpolated values in x-direction
            interped = [];
            for col = im_part
                interped(end+1) = intpolcub_single(col, dy);
            end
            out(out_idx_row, out_idx_col) = intpolcub_single(interped, dx);
        end
    end
end

%%%%%%%%%%%%%%%%% vector implementation using W(x) %%%%%%%%%%%%%%%%
function p = intpolcub_img(I1, I2, I3, I4, D, a)
    if nargin == 5
        a = -0.5;
    end
    p = I1 .* cubic12(-D-1, a) + I2 .* cubic01(-D, a) + I3 .* cubic01(-D+1, a) + I4 .* cubic12(-D+2, a);
end

function out = bicubic_conv_vec(img, scale, padding)
    if nargin == 2
        padding = 'symmetric';
    end
    
    % calculate the new pixel indices in terms of the old pixel indices
    [idx_rows, idx_cols] = get_sample_grid(img, scale);

    % padding
    P = padarray(img, [2, 2], padding);
    padded_idx_rows = idx_rows + 2;
    padded_idx_cols = idx_cols + 2;

    % interpolate in y-direction
    int_idx_rows = floor(padded_idx_rows);
    d_rows = padded_idx_rows - int_idx_rows;

    P_y = intpolcub_img(P(int_idx_rows-1, :, :), P(int_idx_rows, :, :), P(int_idx_rows+1, :, :), P(int_idx_rows+2, :, :), d_rows');

    % interpolate in x-direction
    int_idx_cols = floor(padded_idx_cols);
    d_cols = padded_idx_cols - int_idx_cols;

    out = intpolcub_img(P_y(:, int_idx_cols-1, :), P_y(:, int_idx_cols, :), P_y(:, int_idx_cols+1, :), P_y(:, int_idx_cols+2, :), d_cols);
end

%%%%%%%%%%%%%%%%%%%%%%%%% vector implementation 2
function out = bicubic_conv_matrix(img, scale, padding)
    if nargin == 2
        padding = 'symmetric';
    end

    % calculate the new pixel indices in terms of the old pixel indices
    [idx_rows, idx_cols] = get_sample_grid(img, scale);

    % padding
    P = padarray(img, [2, 2], padding);
    % padded_idx_rows = idx_rows' + 1; % for 'prepare all rows' (slower)
    padded_idx_rows = idx_rows' + 2;   % for 'prepare only indexed rows' (faster)
    padded_idx_cols = idx_cols + 1;    % for 'prepare all cols' (faster)
    % padded_idx_cols = idx_cols + 2;  % for 'prepare only indexed cols' (slower)

    % interpolate in y-direction
    int_idx_rows = floor(padded_idx_rows);
    d_rows = padded_idx_rows - int_idx_rows;

    % prepare all rows
    % I0 =                      2*P(2:end-2, :, :);
    % I1 =  -P(1:end-3, :, :)                      +   P(3:end-1, :, :);
    % I2 = 2*P(1:end-3, :, :) - 5*P(2:end-2, :, :) + 4*P(3:end-1, :, :) - P(4:end, :, :);
    % I3 =  -P(1:end-3, :, :) + 3*P(2:end-2, :, :) - 3*P(3:end-1, :, :) + P(4:end, :, :);

    % P_y = (I0(int_idx_rows, :, :) + d_rows .* I1(int_idx_rows, :, :) + d_rows.^2 .* I2(int_idx_rows, :, :) + d_rows.^3 .* I3(int_idx_rows, :, :)) / 2;

    % prepare only indexed rows
    P_1 = P(int_idx_rows-1, :, :);
    P0  = P(int_idx_rows, :, :);
    P1  = P(int_idx_rows+1, :, :);
    P2  = P(int_idx_rows+2, :, :);
    I0 =         2*P0;
    I1 =  -P_1        +   P1;
    I2 = 2*P_1 - 5*P0 + 4*P1 - P2;
    I3 =  -P_1 + 3*P0 - 3*P1 + P2;

    P_y = (I0 + d_rows .* I1 + d_rows.^2 .* I2 + d_rows.^3 .* I3) / 2;

    % interpolate in x-direction
    int_idx_cols = floor(padded_idx_cols);
    d_cols = padded_idx_cols - int_idx_cols;

    % prepare all cols
    I0 =                        2*P_y(:, 2:end-2, :);
    I1 =  -P_y(:, 1:end-3, :)                        +   P_y(:, 3:end-1, :);
    I2 = 2*P_y(:, 1:end-3, :) - 5*P_y(:, 2:end-2, :) + 4*P_y(:, 3:end-1, :) - P_y(:, 4:end, :);
    I3 =  -P_y(:, 1:end-3, :) + 3*P_y(:, 2:end-2, :) - 3*P_y(:, 3:end-1, :) + P_y(:, 4:end, :);

    out = (I0(:, int_idx_cols, :) + d_cols .* I1(:, int_idx_cols, :) + d_cols.^2 .* I2(:, int_idx_cols, :) + d_cols.^3 .* I3(:, int_idx_cols, :)) / 2;

    % prepare only indexed cols
    % P_y_1 = P_y(:, int_idx_cols-1, :);
    % P_y0  = P_y(:, int_idx_cols, :);
    % P_y1  = P_y(:, int_idx_cols+1, :);
    % P_y2  = P_y(:, int_idx_cols+2, :);
    % I0 =           2*P_y0;
    % I1 =  -P_y_1          +   P_y1;
    % I2 = 2*P_y_1 - 5*P_y0 + 4*P_y1 - P_y2;
    % I3 =  -P_y_1 + 3*P_y0 - 3*P_y1 + P_y2;

    % out = (I0 + d_cols .* I1 + d_cols.^2 .* I2 + d_cols.^3 .* I3) / 2;
end

% simple = repmat([4, 2, 1, 5], 4, 1);

% s3_bicubic_conv_loop = bicubic_conv_loop(simple, 3/4)
% s3_bicubic_conv_vec = bicubic_conv_vec(simple, 3/4)
% s3_bicubic_conv_matrix = bicubic_conv_matrix(simple, 3/4)
% s3_imresize = imresize(simple, 3/4, 'bicubic')

% s5_bicubic_conv_loop = bicubic_conv_loop(simple, 5/4)
% s5_bicubic_conv_vec = bicubic_conv_vec(simple, 5/4)
% s5_bicubic_conv_matrix = bicubic_conv_matrix(simple, 5/4)
% s5_imresize = imresize(simple, 5/4, 'bicubic')

m = magic(4);

m3_bicubic_conv_loop = bicubic_conv_loop(m, 3/4)
m3_bicubic_conv_vec = bicubic_conv_vec(m, 3/4)
m3_bicubic_conv_matrix = bicubic_conv_matrix(m, 3/4)
% m3_imresize = imresize(m, 3/4, 'bicubic')

m5_bicubic_conv_loop = bicubic_conv_loop(m, 5/4)
m5_bicubic_conv_vec = bicubic_conv_vec(m, 5/4)
m5_bicubic_conv_matrix = bicubic_conv_matrix(m, 5/4)
% m5_imresize = imresize(m, 5/4, 'bicubic')


r_big = randi(16000, 8000, 8000, 'int16');
% tic; bicubic_conv_loop(r_big, 100/8000); toc; % 4.40579 seconds.
tic; bicubic_conv_vec(r_big, 100/8000); toc;    % 0.144027 seconds.
tic; bicubic_conv_matrix(r_big, 100/8000); toc; % 4.79519 seconds. Now: 0.251252 seconds.



r_small = randi(16000, 100, 100, 'int16');
% tic; bicubic_conv_loop(r_small, 8000/100); toc; % not feasible
tic; bicubic_conv_vec(r_small, 8000/100); toc;    % 1.83731 seconds.
tic; bicubic_conv_matrix(r_small, 8000/100); toc; % 2.4178 seconds. Now: 2.44018 seconds.
