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.