Why does upsampling and interpolation by convolution introduce a shift compared to imresize?

Question

For purposes of understanding the process - not for any practical purpose for which I could use imresize - I wanted to show that 2x upsampling followed by convolution with an appropriate kernel (presumably the 2d equivalent of a "triangle") is equivalent to imresize(x,scalingFactor,'bilinear'), However I cannot do it perfectly. It "almost" works in the sense that qualitatively it looks very similar but it is as if it is slightly "shifted" cf. code below. Moreover for upsampling higher than ×2 I don't know how to build the correct kernel. I tried below, with zero "stuffing" but that is not correct. Any help and explanations/corrections are greatly appreciated.

N = 64;
%x = phantom(N)+0.25;
x = imresize(im2double(imread('cameraman.tif')),[N N]);
s = 2;

x_matlab_bilinear = imresize(x,s,'bilinear');

kernel1D = [0.5, 0.5*ones(1,s-2), 1, 0.5*ones(1,s-2), 0.5];
bilinearKernel =  kernel1D.'*kernel1D;

x_up = zeros(s*N,s*N);
x_up(1:s:end,1:s:end) = x;
x_custom_bilinear = conv2(x_up, bilinearKernel,'same');

figure;
imagesc(bilinearKernel);

figure;
subplot(141);
imagesc(x_up);colormap('gray');
subplot(142);
imagesc(x_matlab_bilinear);colormap('gray');title('x up matlab bilinear');
subplot(143);
imagesc(x_custom_bilinear);colormap('gray');title('x up custom bilinear');
subplot(144);
imagesc(x_matlab_bilinear-x_custom_bilinear);colorbar;title('diff');

For e.g. x8 upsampling I would use the following kernel (generated by my formula):

above is my kernel for interpolating a 2 times upsampled image.

For the delta-impulse image the results of 2 times upsampling are the following for the same above code but just replacing the first lines by

N=6; x=zeros(N); x(3,3)=1;

:

Answer 1

The linear interpolation kernel should be a triangle. You can construct it like this:

kernel1D = [(1:s)/s, (s-1:-1:1)/s];

The difference between MATLAB's imresize and your method for s=2 is a half-pixel shift. This is equivalent to shifting the kernel by half a pixel. So instead of a kernel [1/2, 1, 1/2], you'd use a kernel [1/4,3/4,3/4,1/4]. I was able to reproduce MATLAB's result like this (up to differences at the boundary, where MATLAB pads the input by repeating values, rather than adding zeros like conv2 does):

kernel1D = [1/4,3/4,3/4,1/4];  % kernel shifted by half a pixel!
bilinearKernel =  kernel1D.'*kernel1D;

x_up = zeros(s*N,s*N);
x_up(2:s:end,2:s:end) = x;   % note the shift here!
x_custom_bilinear = conv2(x_up, bilinearKernel,'same');

I'm guessing that the half-pixel shift makes the operation more symmetric: note that in your original result, you have a black column of pixels on the right. Shift by half a pixel, and now both the left and right columns of pixels contain the same fraction of "invented" (extrapolated) data.