What is the algorithm used for interpolation in the matlab imresize function?

Question

I use the Matlab / Octave imresize() function, which oversambles the given 2D array. I want to understand how the specific interpolation algorithm used in imresize .

(I use an octave on the windows)

eg.

 A = 1 2 3 4 

is a two-dimensional array. Then i use the command

 b=imresize(a,2,'linear'); 

basically upsampling row and columns by 2.

Output signal

 1.0000 1.3333 1.6667 2.0000 1.6667 2.0000 2.3333 2.6667 2.3333 2.6667 3.0000 3.3333 3.0000 3.3333 3.6667 4.0000 

I do not understand how this linear interpolation works. They say that linear interpolation bi is used , but how does it fill the data at the borders and how does it get the result that it gets?

Second example: For

 A = 1 2 3 4 5 6 7 8 0 1 2 3 1 2 3 4 

how imresize(a,1.5,'linear') give the following output?

 1.00000 1.60000 2.20000 2.80000 3.40000 4.00000 3.40000 4.00000 4.60000 5.20000 5.80000 6.40000 4.00000 4.60000 5.20000 5.80000 6.40000 7.00000 1.00000 1.60000 2.20000 2.80000 3.40000 4.00000 0.40000 1.00000 1.60000 2.20000 2.80000 3.40000 1.00000 1.60000 2.20000 2.80000 3.40000 4.00000 

Answer 1

As you can see, in your example, each corner point is one of your original input values.

Intermediate values ​​are output via linear interpolation in each direction. So, for example, to calculate b(3,2) :

• b(1,2) is 1/3 of the way between b(1,1) and b(1,4) . So:

   b(1,2) = (1/3)*b(1,4) + (2/3)*b(1,1) 

• b(4,2) is 1/3 of the way between b(4,1) and b(4,4) . So:

   b(4,2) = (1/3)*b(4,4) + (2/3)*b(4,1) 

• b(3,2) is 2/3 of the path between b(1,2) and b(4,2) . So:

   b(3,2) = (2/3)*b(4,2) + (1/3)*b(1,2) 

Answer 2

The following code shows how to perform bilinear interpolation using INTERP2 :

 A = [1 2; 3 4]; SCALE = 2; xi = linspace(1,size(A,2),SCALE*size(A,2)); %# interpolated horizontal positions yi = linspace(1,size(A,1),SCALE*size(A,1)); %# interpolated vertical positions [XY] = meshgrid(1:size(A,2),1:size(A,1)); %# pixels X-/Y-coords [XI YI] = meshgrid(xi,yi); %# interpolated pixels X-/Y-coords B = interp2(X,Y,A, XI,YI, '*linear'); %# interp values at these positions 

the result is consistent with the output of the Octave code:

 B = 1 1.3333 1.6667 2 1.6667 2 2.3333 2.6667 2.3333 2.6667 3 3.3333 3 3.3333 3.6667 4 

---

I should mention that I have different results between MATLAB and Octave IMRESIZE output. For example, this is what I get when I do the following in MATLAB on the matrix A=[1 2; 3 4] A=[1 2; 3 4] :

 >> B = imresize([1 2; 3 4], 2, 'bilinear') B = 1 1.25 1.75 2 1.5 1.75 2.25 2.5 2.5 2.75 3.25 3.5 3 3.25 3.75 4 

which suggests that the MATLAB implementation does something extra ... Unfortunately, it is not easy to read the source code of IMRESIZE, especially since at some point it calls the function of the compiled MEX (without the form of the source code).

As a side note, there seems to be an older version of this function: IMRESIZE_OLD (purely implemented in m-code). From what I could understand, he performs some kind of affine transformation on the image. Perhaps someone more familiar with this technique may shed light on this subject ...

Answer 3

I adapted the MATLAB imresize function for Java:

 import java.util.ArrayList; import java.util.List; public class MatlabResize { private static final double TRIANGLE_KERNEL_WIDTH = 2; public static double[][] resizeMatlab(double[][] data, int out_y, int out_x) { double scale_x = ((double)out_x)/data[0].length; double scale_y = ((double)out_y)/data.length; double[][][] weights_indizes = contribution(data.length, out_y, scale_y, TRIANGLE_KERNEL_WIDTH); double[][] weights = weights_indizes[0]; double[][] indices = weights_indizes[1]; final double[][] result = new double[out_y][data[0].length]; double value = 0; for (int p=0; p<result[0].length; p++) { for (int i=0; i<weights.length; i++) { value = 0; for (int j=0; j<indices[0].length; j++) { value += weights[i][j] * data[(int)indices[i][j]][p]; } result[i][p] = value; } } weights_indizes = contribution(data[0].length, out_x, scale_x, TRIANGLE_KERNEL_WIDTH); weights = weights_indizes[0]; indices = weights_indizes[1]; final double[][] result2 = new double[result.length][out_x]; for (int p=0; p<result.length; p++) { for (int i=0; i<weights.length; i++) { value = 0; for (int j=0; j<indices[0].length; j++) { value += weights[i][j] * result[p][(int)indices[i][j]]; } result2[p][i] = value; } } return result2; } public static double[][] resizeMatlab(double[][] data, double scale) { int out_x = (int)Math.ceil(data[0].length * scale); int out_y = (int)Math.ceil(data.length * scale); double[][][] weights_indizes = contribution(data.length, out_y, scale, TRIANGLE_KERNEL_WIDTH); double[][] weights = weights_indizes[0]; double[][] indices = weights_indizes[1]; final double[][] result = new double[out_y][data[0].length]; double value = 0; for (int p=0; p<result[0].length; p++) { for (int i=0; i<weights.length; i++) { value = 0; for (int j=0; j<indices[0].length; j++) { value += weights[i][j] * data[(int)indices[i][j]][p]; } result[i][p] = value; } } weights_indizes = contribution(data[0].length, out_x, scale, TRIANGLE_KERNEL_WIDTH); weights = weights_indizes[0]; indices = weights_indizes[1]; final double[][] result2 = new double[result.length][out_x]; for (int p=0; p<result.length; p++) { for (int i=0; i<weights.length; i++) { value = 0; for (int j=0; j<indices[0].length; j++) { value += weights[i][j] * result[p][(int)indices[i][j]]; } result2[p][i] = value; } } return result2; } private static double[][][] contribution(int length, int output_size, double scale, double kernel_width) { if (scale < 1.0) { kernel_width = kernel_width/scale; } final double[] x = new double[output_size]; for (int i=0; i<x.length; i++) { x[i] = i+1; } final double[] u = new double[output_size]; for (int i=0; i<u.length; i++) { u[i] = x[i]/scale + 0.5*(1 - 1/scale); } final double[] left = new double[output_size]; for (int i=0; i<left.length; i++) { left[i] = Math.floor(u[i] - kernel_width/2); } int P = (int)Math.ceil(kernel_width) + 2; final double[][] indices = new double[left.length][P]; for (int i=0; i<left.length; i++) { for (int j=0; j<=P-1; j++) { indices[i][j] = left[i] + j; } } double[][] weights = new double[u.length][indices[0].length]; for (int i=0; i<u.length; i++) { for (int j=0; j<indices[i].length; j++) { weights[i][j] = u[i] - indices[i][j]; } } if (scale < 1.0) { weights = triangleAntiAliasing(weights, scale); } else { weights = triangle(weights); } double[] sum = Matlab.sum(weights, 2); for (int i=0; i<weights.length; i++) { for (int j=0; j<weights[i].length; j++) { weights[i][j] = weights[i][j] / sum[i]; } } for (int i=0; i<indices.length; i++) { for (int j=0; j<indices[i].length; j++) { indices[i][j] = Math.min(Math.max(indices[i][j], 1.0), length); } } sum = Matlab.sum(weights, 1); int a = 0; final List<Integer> list = new ArrayList<Integer>(); for (int i=0; i<sum.length; i++) { if (sum[i] != 0.0) { a++; list.add(i); } } final double[][][] result = new double[2][weights.length][a]; for (int i=0; i<weights.length; i++) { for (int j=0; j<list.size(); j++) { result[0][i][j] = weights[i][list.get(j)]; } } for (int i=0; i<indices.length; i++) { for (int j=0; j<list.size(); j++) { result[1][i][j] = indices[i][list.get(j)]-1; //java indices start by 0 and not by 1 } } return result; } private static double[][] triangle(final double[][] x) { for (int i=0; i<x.length; i++) { for (int j=0; j<x[i].length; j++) { if (-1.0 <= x[i][j] && x[i][j] < 0.0) { x[i][j] = x[i][j] + 1; } else if (0.0 <= x[i][j] && x[i][j] < 1.0) { x[i][j] = 1 - x[i][j]; } else { x[i][j] = 0; } } } return x; } private static double[][] triangleAntiAliasing(final double[][] x, final double scale) { for (int i=0; i<x.length; i++) { for (int j=0; j<x[i].length; j++) { x[i][j] = x[i][j] * scale; } } for (int i=0; i<x.length; i++) { for (int j=0; j<x[i].length; j++) { if (-1.0 <= x[i][j] && x[i][j] < 0.0) { x[i][j] = x[i][j] + 1; } else if (0.0 <= x[i][j] && x[i][j] < 1.0) { x[i][j] = 1 - x[i][j]; } else { x[i][j] = 0; } } } for (int i=0; i<x.length; i++) { for (int j=0; j<x[i].length; j++) { x[i][j] = x[i][j] * scale; } } return x; } }