Какой алгоритм используется для интерполяции в функции imresize Matlab?

Question

Я использую функцию Matlab / Octave imresize(), которая повторно дискретизирует данный 2D-массив.Я хочу понять, как работает конкретный алгоритм интерполяции, используемый в imresize.

(я использую октаву в окнах)

например,

A =  1 2 
     3 4

- это двумерный массив,Затем я использую команду

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

, в основном повышая дискретизацию строк и столбцов на 2.

Вывод

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

Я не понимаю, какова эта линейная интерполяцияза работой.Говорят, что он использует bi линейную интерполяцию, но как он дополняет данные на границах и как он получает вывод, который он получает?

Второй пример: Для

A = 

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

как imresize(a,1.5,'linear') дает следующий вывод?

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

В следующем коде показано, как выполнить билинейную интерполяцию с использованием 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
[X Y] = 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

результат согласуется с выводом кода Octave:

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

---

Я должен отметить, что я получаю разные результаты между MATLAB и Octave IMRESIZE.Например, это то, что я получаю, когда выполняю следующее в MATLAB для матрицы 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

, что говорит о том, что реализация MATLAB делает что-то дополнительное ... К сожалению, нелегко прочитатьIMRESIZE исходный код, особенно потому, что в какой-то момент он вызывает MEX-скомпилированную функцию (без формы исходного кода).

В качестве примечания, похоже, есть и более старая версия этой функции: IMRESIZE_OLD(чисто реализовано в м-коде).Из того, что я мог понять, он выполняет своего рода аффинное преобразование изображения.Возможно, кто-то, более знакомый с техникой, мог бы пролить свет на эту тему ...

Answer 2

Я адаптировал функцию MATLAB imresize для 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;
    }
}

Answer 3

Как вы можете видеть, в вашем примере каждая угловая точка является одним из ваших исходных входных значений.

Промежуточные значения получены с помощью линейной интерполяции в каждом направлении. Так, например, для расчета b(3,2):

• b(1,2) составляет 1/3 пути между b(1,1) и b(1,4). Итак:

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

• b(4,2) составляет 1/3 пути от b(4,1) до b(4,4). Итак:

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

• b(3,2) составляет 2/3 пути от b(1,2) до b(4,2). Итак:

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