1; % dummy code to have local functions
% for MATLAB shift the code from bottom to here and add 'Antialiasing', false in imresize

% adapted from https://stackoverflow.com/a/39948966/2414411
function f = cubic(x, a)
    absx = abs(x);
    absx2 = absx.^2;
    absx3 = absx.^3;
    f = ((a+2)*absx3 - (a+3)*absx2 + 1) .* (absx <= 1) + ...
        (a*absx3 -5*a*absx2 + 8*a*absx - 4*a) .* ((1 < absx) & (absx <= 2));
end

% adapted from https://stackoverflow.com/a/39948966/2414411
function f = intpolcub(x, d, a)
    if nargin == 2
        a = -0.5;
    end
    f = x(1)*cubic(-d-1, a) + x(2)*cubic(-d, a) + x(3)*cubic(-d+1, a) + x(4)*cubic(-d+2, 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

% this method matches the results of MATLAB's imresize
function out = bicubic_conv(in, 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(in, scale);

    % apply padding
    padded = padarray(in, [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;
            
            im_part = padded(fir-1:fir+2, fic-1:fic+2);
            
            interped = [];
            for col = im_part
                interped(end+1) = intpolcub(col, dy);
            end
            out(out_idx_row, out_idx_col) = intpolcub(interped, dx);
        end
    end
end

% in MATLAB this gives the same results as bicubic_conv, in Octave not
function out = bicubic_interp(in, 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(in, scale);

    % apply padding
    pad_size = 2; % for Octave 1 is enough, for MATLAB 2 is required
    padded = padarray(in, [pad_size, pad_size], padding);
    padded_idx_rows = idx_rows + pad_size;
    padded_idx_cols = idx_cols + pad_size;

    % interpolate
    out = interp2(padded, padded_idx_rows', padded_idx_cols, 'cubic');
end

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

s3_bicubic_conv = bicubic_conv(simple, 3/4)
s3_bicubic_interp = bicubic_interp(simple, 3/4)
% s3_imresize = imresize(simple, 3/4, 'bicubic')

s5_bicubic_conv = bicubic_conv(simple, 5/4)
s5_bicubic_interp = bicubic_interp(simple, 5/4)
% s5_imresize = imresize(simple, 5/4, 'bicubic')

m = magic(4);

m3_bicubic_conv = bicubic_conv(m, 3/4)
m3_bicubic_interp = bicubic_interp(m, 3/4)
% m3_imresize = imresize(m, 3/4, 'bicubic')

m5_bicubic_conv = bicubic_conv(m, 5/4)
m5_bicubic_interp = bicubic_interp(m, 5/4)
% m5_imresize = imresize(m, 5/4, 'bicubic')

