How does nearest neighbor interpolation work in MATLAB?

Question

I want to find out how nearest neighbor interpolation works in MATLAB. I have input data :

A = [1 4 7 4 3 6]     % 6 digit vector

I use the following MATLAB code :

B = imresize(A,[1 9],'nearest');

I get the following result :

[1,4,4,7,4,4,3,6,6]

Solving by hand, I get this result :

[1 4 4 7 4 4 3 3 6]

Can you please guide me? Am I going wrong somewhere?

Solution

If you apply regular interpolation using interp1, it will give you the result you computed by hand:

>> N = 9;
>> B = interp1(linspace(0,1,numel(A)), A, linspace(0,1,N), 'nearest')
B =
     1     4     4     7     4     4     3     3     6

Some time ago, I went through the source code of IMRESIZE trying to understand how it works. See this post for a summary. At some point the code will call an private MEX-function (no corresponding source code available), but the comments are enough to understand the implementation.

For what it's worth, there an also a function imresize_old which provides an older implementation of imresize (used in version R2006b and earlier). It gave yet another different result:

>> B = imresize(A, [1 N], 'nearest')
B =
     1     4     4     7     4     4     3     6     6

>> B = imresize_old(A, [1 N], 'nearest')
B =
     1     4     4     7     7     4     3     6     6

What's more it was previously observed that the implementation between MATLAB and Octave also differed in some cases.

---

EDIT:

As you noted, in some cases you have to be mindful about floating-point limitations when working with interp1. So we could do the interpolation by choosing x-numbers to be between [0,1] range, or a more stable range like [1,numel(A)]. Because of rounding errors in edge cases, this might give different results.

For example compare the two codes below:

% interpolation in [0,1]
N = 11;
y = [1 4 7 4 3 6];
x = linspace(0,1,numel(y));
xi = linspace(0,1,N);
yi = interp1(x, y, xi, 'nearest');

% print numbers with extended precision
fprintf('%.17f %g\n',[x;y])
fprintf('%.17f %g\n',[xi;yi])

against:

% interpolation in [1,k]
N = 11;
y = [1 4 7 4 3 6];
x = 1:numel(y);
xi = linspace(1,numel(y),N);
yi = interp1(x, y, xi, 'nearest');

% print numbers with extended precision
fprintf('%.17f %g\n',[x;y])
fprintf('%.17f %g\n',[xi;yi])

Here is the output nicely formatted:

--------------------------------------------------------
       [0,1] RANGE        |         [1,k] RANGE
--------------------------------------------------------     
        xi           yi   |            xi            yi
--------------------------------------------------------
0.00000000000000000  1    |     1.00000000000000000  1 |
0.20000000000000001  4    |     2.00000000000000000  4 |
0.40000000000000002  7    |     3.00000000000000000  7 |
0.59999999999999998  4    |     4.00000000000000000  4 |  INPUT
0.80000000000000004  3    |     5.00000000000000000  3 |
1.00000000000000000  6    |     6.00000000000000000  6 |
--------------------------------------------------------
0.00000000000000000  1    |     1.00000000000000000  1 |
0.10000000000000001  4    |     1.50000000000000000  4 |
0.20000000000000001  4    |     2.00000000000000000  4 |
0.29999999999999999  4    |     2.50000000000000000  7 |
0.40000000000000002  7    |     3.00000000000000000  7 |
0.50000000000000000  4    |     3.50000000000000000  4 | OUTPUT
0.59999999999999998  4    |     4.00000000000000000  4 |
0.69999999999999996  4    |     4.50000000000000000  3 |
0.80000000000000004  3    |     5.00000000000000000  3 |
0.90000000000000002  6    |     5.50000000000000000  6 |
1.00000000000000000  6    |     6.00000000000000000  6 |
--------------------------------------------------------

So you can see that some numbers are not exactly representable in double-precision when working in the [0,1] range. So 0.3 which is supposed to be in the middle [0.2, 0.4], turns to be closer to the lower end 0.2 than 0.4 because of roundoff error. While on the other side, 2.5 is exactly in the middle of [2,3] (all numbers exactly represented), and is assigned to the upper end 3 using nearest neighbor.

Also be aware that colon and linspace can produce different outputs sometimes:

>> (0:0.1:1)' - linspace(0,1,11)'
ans =
            0
            0
            0
   5.5511e-17
            0
            0
            0
            0
            0
            0
            0