Nearest-neighbor interpolation algorithm in MATLAB

Question

I am trying to write my own function for scaling up an input image by using the Nearest-neighbor interpolation algorithm. The bad part is I am able to see how it works but cannot find the algorithm itself. I will be grateful for any help.

Here's what I tried for scaling up the input image by a factor of 2:

function output = nearest(input)
[x,y]=size(input);
output = repmat(uint8(0),x*2,y*2);
[newwidth,newheight]=size(output);
for i=1:y
    for j=1:x
        xloc = round ((j * (newwidth+1)) / (x+1));
        yloc = round ((i * (newheight+1)) / (y+1));
        output(xloc,yloc) = input(j,i);
    end
end

Here is the output after Mark's suggestion

Answer 1

This answer is more explanatory than trying to be concise and efficient. I think gnovice's solution is best in that regard. In case you are trying to understand how it works, keep reading...

Now the problem with your code is that you are mapping locations from the input image to the output image, which is why you are getting the spotty output. Consider an example where input image is all white and output initialized to black, we get the following:

What you should be doing is the opposite (from output to input). To illustrate, consider the following notation:

1           c         1                 scaleC*c
+-----------+ 1       +----------------------+ 1
|    |      |         |        |             |
|----o      |   <===  |        |             |
|  (ii,jj)  |         |--------o             |
+-----------+ r       |      (i,j)           |
  inputImage          |                      |
                      |                      |
                      +----------------------+ scaleR*r
                            ouputImage

Note: I am using matrix notation (row/col), so:
  i ranges on [1,scaleR*r] , and j on [1,scaleC*c]
  and ii on [1,r], jj on [1,c]

The idea is that for each location (i,j) in the output image, we want to map it to the "nearest" location in the input image coordinates. Since this is a simple mapping we use the formula that maps a given x to y (given all the other params):

 x-minX      y-minY
--------- = ---------
maxX-minX   maxY-minY

in our case, x is the i/j coordinate and y is the ii/jj coordinate. Therefore substituting for each gives us:

jj = (j-1)*(c-1)/(scaleC*c-1) + 1
ii = (i-1)*(r-1)/(scaleR*r-1) + 1

Putting pieces together, we get the following code:

% read a sample image
inputI = imread('coins.png');
[r,c] = size(inputI);
scale = [2 2];        % you could scale each dimension differently

outputI = zeros(scale(1)*r,scale(2)*c, class(inputI));

for i=1:scale(1)*r
    for j=1:scale(2)*c
        % map from output image location to input image location
        ii = round( (i-1)*(r-1)/(scale(1)*r-1)+1 );
        jj = round( (j-1)*(c-1)/(scale(2)*c-1)+1 );

        % assign value
        outputI(i,j) = inputI(ii,jj);
    end
end

figure(1), imshow(inputI)
figure(2), imshow(outputI)

Answer 2

A while back I went through the code of the imresize function in the MATLAB Image Processing Toolbox to create a simplified version for just nearest neighbor interpolation of images. Here's how it would be applied to your problem:

%# Initializations:

scale = [2 2];              %# The resolution scale factors: [rows columns]
oldSize = size(inputImage);                   %# Get the size of your image
newSize = max(floor(scale.*oldSize(1:2)),1);  %# Compute the new image size

%# Compute an upsampled set of indices:

rowIndex = min(round(((1:newSize(1))-0.5)./scale(1)+0.5),oldSize(1));
colIndex = min(round(((1:newSize(2))-0.5)./scale(2)+0.5),oldSize(2));

%# Index old image to get new image:

outputImage = inputImage(rowIndex,colIndex,:);

Another option would be to use the built-in interp2 function, although you mentioned not wanting to use built-in functions in one of your comments.

EDIT: EXPLANATION

In case anyone is interested, I thought I'd explain how the solution above works...

newSize = max(floor(scale.*oldSize(1:2)),1);

First, to get the new row and column sizes the old row and column sizes are multiplied by the scale factor. This result is rounded down to the nearest integer with floor. If the scale factor is less than 1 you could end up with a weird case of one of the size values being 0, which is why the call to max is there to replace anything less than 1 with 1.

rowIndex = min(round(((1:newSize(1))-0.5)./scale(1)+0.5),oldSize(1));
colIndex = min(round(((1:newSize(2))-0.5)./scale(2)+0.5),oldSize(2));

Next, a new set of indices is computed for both the rows and columns. First, a set of indices for the upsampled image is computed: 1:newSize(...). Each image pixel is considered as having a given width, such that pixel 1 spans from 0 to 1, pixel 2 spans from 1 to 2, etc.. The "coordinate" of the pixel is thus treated as the center, which is why 0.5 is subtracted from the indices. These coordinates are then divided by the scale factor to give a set of pixel-center coordinates for the original image, which then have 0.5 added to them and are rounded off to get a set of integer indices for the original image. The call to min ensures that none of these indices are larger than the original image size oldSize(...).

outputImage = inputImage(rowIndex,colIndex,:);

Finally, the new upsampled image is created by simply indexing into the original image.