TL;DR
The results are wrong, because there was a bug in GNU Octave's imresize. It has been fixed recently for bilinear interpolation (bug report, commit). Scroll down to "Interpolation of pixels" for the correct explanation how images are interpolated.
Interpolation of samples
Let's start with linear interpolation of samples:
0 1/3 2/3 1
| | | |
a=1----+----+----2=b
You blend over from a to b by using
f(x) = (1 - x) * a + x * b.
Some examples:
- f(0) = (1 - 0) * a + 0 * b = a = 1
- f(1/3) = (1 - 1/3) * a + 1/3 * b = 2/3 + 2/3 = 4/3 = 1.3333
- f(1/2) = (1 - 1/2) * a + 1/2 * b = (a + b) / 2 = 1.5
- f(2/3) = (1 - 2/3) * a + 2/3 * b = 1/3 + 4/3 = 5/3 = 1.6667
- f(1) = (1 - 1) * a + 1 * b = b = 2
This corresponds to the first line of your first example. In bilinear interpolation the simple linear interpolation is used in x or y direction. In general simple linear interpolation is not used in diagonal or any other direction (your first example is a degenerate case).
0 1/3 1
| | | |
0 a=1---f(x)--+----2=b
| | | |
-+----+----+----+-
| | | |
2/3 -+---???---+----+-
| | | |
1 c=3---g(x)--+----4=d
| | | |
How are the other points calculated then? We use the simple linear interpolation for top and bottom row in x direction and then interpolate the results in y direction:
- In x-direction, top row: f(x) = (1 - x) * a + x * b, example: f(1/3) = 4/3 = 1.3333
- In x-direction, bottom row: g(x) = (1 - x) * c + x * d, example: g(1/3) = 10/3 = 3.3333
- In y-direction, interpolation column: h(y) = (1 - y) * f(x) + y * g(x), example: h(2/3) = 8/3 = 2.6667
Looking at the last equation, we could also put in f(x) and g(x) and get:
h(x, y) = (1 - x) * (1 - y) * a + x * (1 - y) * b + (1 - x) * y * c + x * y * d
This is what you got.
In the second example the points are just a bit different, because you transform from 4 points to 6 in each direction:
old: 0 1 2 3 (sample grid)
| | | |
+-----+---+-+-----+-+---+-----+
| | | | | |
new: 0 3/5 6/5 9/5 12/5 3 (interpolation grid)
This is valid for x and y direction in your second example. To use the formula above, you must map each square to [0, 1] x [0, 1].
That's theory. Octave uses interp2 for bilinear interpolation internally. To use interp2 you specify the matrix with samples and the grid that defines the interpolation points:
A = [1, 2;
3, 4];
xi = linspace(1, size(A, 2), 4);
yi = linspace(1, size(A, 1), 4)';
B = interp2(A, xi, yi)
This gives the results you got, but they are wrong!
Interpolation of pixels
The basics for bilinear interpolation as explained above are still valid, but the interpolation grid is wrong. This is because an image does not consist of sample points, but of pixels. Pixels are areas represented by its mean value. So actually the pixels of the image look like this:
0.5 1 1.5 2 2.5
0.5 +-------------------+-------------------+
| | |
| | |
| | |
1 | o | o |
| | |
| | |
| | |
1.5 +-------------------+-------------------+
| | |
| | |
| | |
2 | o | o |
| | |
| | |
| | |
2.5 +-------------------+-------------------+
So the top left pixel's area is [0.5, 1.5] x [0.5, 1.5] with its center at (1, 1). And what you want by upscaling with 2 are the following new pixels (in coordinate space of the old grid, since the image covers still the same area):
0.5 0.75 1 1.25 1.5 1.75 2 2.25 2.5
0.5 +---------+---------+---------+---------+
| | | | |
0.75 | x | x | x | x |
| | | | |
1 +---------o---------+---------o---------+
| | | | |
1.25 | x | x | x | x |
| | | | |
1.5 +---------+---------+---------+---------+
| | | | |
1.75 | x | x | x | x |
| | | | |
2 +---------o---------+---------o---------+
| | | | |
2.25 | x | x | x | x |
| | | | |
2.5 +---------+---------+---------+---------+
Now you take the new centers x as interpolation grid and the old centers o as sample grid. You see, that the new border pixels actually need extrapolation. We assume that it extrapolates constant, so we could pad the array to make an interpolation again or limit the interpolation grid. As code using interp2 you can do:
A = [1, 2;
3, 4];
xi = linspace(0.75, 2.25, 4);
yi = linspace(0.75, 2.25, 4)';
xi(xi < 1) = 1; xi(xi > 2) = 2;
yi(yi < 1) = 1; yi(yi > 2) = 2;
B = interp2(A, xi, yi)
Here is a more general solution (only valid for integer output sizes), inspired by Amro's comments below his answer. If you allow scale factors that result in floating point output size.
On non-integer output sizes the number of new pixels will be such that the last pixel overlaps. So for example with a scale factor of 5/4 = 1.25 the pixel size will be 1 / (5/4) = 4/5 = 0.8. Hence, scaling a 2x2 image with 1.25 yields a 3x3 image. The old pixel centers (sample grid) are at 1 and 2 and the new pixel centers (interpolation grid) are at 0.9, 1.7 and 2.5.
0.5 1.5 2.5
| 1 | 2 |
old: +---------o---------+---------o---------+
new: +-------x-------+-------x-------+-------x-------+
| 0.9 | 1.7 | 2.5 |
0.5 1.3 2.1 2.9
Here is some code to show this:
img = [1, 2;
3, 4];
% make interpolation grid
scale = 1.25
pixel_size = 1 / scale
out_size = ceil(size(img) / pixel_size)
xi = 0.5 + pixel_size / 2 + (0:out_size(1)-1) / scale
yi = 0.5 + pixel_size / 2 + (0:out_size(2)-1) / scale
% limit interpolation grid to sample grid bounds
xi(xi < 1) = 1; xi(xi > size(img, 2)) = size(img, 2)
yi(yi < 1) = 1; yi(yi > size(img, 1)) = size(img, 1)
% interpolate
scaled_interp = interp2(img, xi, yi', 'linear')
% Octave's imresize does not support anti-aliasing yet
scaled_resize_octave = imresize(img, scale, 'bilinear')
% Matlab's imresize uses anti-aliasing for downscaling, switch off to keep is simple
scaled_resize_matlab = imresize(img, scale, 'bilinear', 'Antialiasing', false)
% yields:
% 1.0000 1.7000 2.0000
% 2.4000 3.1000 3.4000
% 3.0000 3.7000 4.0000
That's all for resizing with bilinear interpolation. For bicubic interpolation Matlab uses symmetric padding and a convolutional algorithm that weights a 4x4 neighborhood. Octave behaves differently (patch incoming).
