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

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

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

e.g.

A =  1 2 
     3 4

представляет собой 2D-массив. Затем я использую команду

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

person goldenmean    schedule 22.06.2011    source источник


Ответы (4)


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

Промежуточные значения получаются с помощью линейной интерполяции в каждом направлении. Так, например, чтобы вычислить 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)
    
person Oliver Charlesworth    schedule 22.06.2011
comment
@Oli - Спасибо. Я понял. Я предполагаю, что расчет весов для другого коэффициента масштабирования примет соответствующие значения. т.е. для моего примера веса были 0,33,0,66, так как 2 заданных образца должны были быть интерполированы до 4 образцов всего, если это было 2 образца, которые должны быть интерполированы на общее количество образцов N, веса были бы 1 / (N-1), 2 / ( N-1) не так ли? - person goldenmean; 22.06.2011
comment
@OliCharlesworth - Отредактировал OP, чтобы добавить второй пример, который выполняет повторную выборку с использованием дробного коэффициента масштабирования. В этом случае неясно, как он рассчитывал свои результаты. Можете ли вы, пожалуйста. проверьте выше. - person goldenmean; 22.06.2011
comment
@Oli - Есть какие-нибудь отзывы о втором примере, который я добавил в OP? - person goldenmean; 22.06.2011
comment
@Oli - также в отредактированных вами точках ответа b (4,2) и b (3,2) вычисления кажутся мне некорректными (веса поменялись местами). Не так ли? - person goldenmean; 22.06.2011
comment
@goldenmean: Да, хорошее место, у меня был последний расчет задом наперед! - person Oliver Charlesworth; 22.06.2011
comment
@Oli - Есть какие-нибудь отзывы о втором примере, который я добавил в OP? - person goldenmean; 22.06.2011
comment
@Golden: это все еще линейная интерполяция для вашего второго примера - вы можете быть сбиты с толку, потому что все задействованные веса основаны на множителе (1/5), который (на единицу больше) желаемого размера вывода минус 1. - person Jonas Heidelberg; 03.02.2012

TL;DR

Результаты неверны, потому что в GNU Octave imresize была ошибка. Недавно это было исправлено для билинейной интерполяции (отчет об ошибке, совершить). Прокрутите вниз до Интерполяции пикселей для правильного объяснения того, как изображения интерполируются.

Интерполяция образцов

Начнем с линейной интерполяции отсчетов:

  0   1/3  2/3   1
  |    |    |    |
a=1----+----+----2=b

Вы переходите от a к b, используя

f(x) = (1 - x) * a + x * b.

Некоторые примеры:

  • 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

Это соответствует первой строке вашего первого примера. При билинейной интерполяции используется простая линейная интерполяция по осям x или y. Обычно простая линейная интерполяция не используется по диагонали или в любом другом направлении (ваш первый пример - вырожденный случай).

       0   1/3        1
       |    |    |    |
 0   a=1---f(x)--+----2=b
       |    |    |    |
      -+----+----+----+-
       |    |    |    |
2/3   -+---???---+----+-
       |    |    |    |
 1   c=3---g(x)--+----4=d
       |    |    |    |

Как тогда рассчитываются другие баллы? Мы используем простую линейную интерполяцию для верхней и нижней строки в направлении x, а затем интерполируем результаты в направлении y:

  • В направлении x, верхний ряд: f (x) = (1 - x) * a + x * b, например: f (1/3) = 4/3 = 1,3333
  • В направлении x нижний ряд: g (x) = (1 - x) * c + x * d, например: g (1/3) = 10/3 = 3,3333
  • В направлении y столбец интерполяции: h (y) = (1 - y) * f (x) + y * g (x), например: h (2/3) = 8/3 = 2,6667

Глядя на последнее уравнение, мы могли бы также подставить f (x) и g (x) и получить:

h(x, y) = (1 - x) * (1 - y) * a + x * (1 - y) * b + (1 - x) * y * c + x * y * d

Это то, что у вас есть.

Во втором примере точки немного отличаются, потому что вы преобразовываете из 4 точек в 6 в каждом направлении:

old: 0         1         2         3 (sample grid)
     |         |         |         |
     +-----+---+-+-----+-+---+-----+
     |     |     |     |     |     |
new: 0    3/5   6/5   9/5  12/5    3 (interpolation grid)

Это верно для направлений x и y во втором примере. Чтобы использовать приведенную выше формулу, вы должны сопоставить каждый квадрат с [0, 1] x [0, 1].

Это теория. Octave использует interp2 для внутренней билинейной интерполяции. Чтобы использовать interp2, вы указываете матрицу с выборками и сетку, которая определяет точки интерполяции:

A = [1, 2;
     3, 4];
xi = linspace(1, size(A, 2), 4);
yi = linspace(1, size(A, 1), 4)';
B = interp2(A, xi, yi)

Это дает результаты, которые вы получили, но они неверны!

Интерполяция пикселей

Основы билинейной интерполяции, как объяснено выше, по-прежнему действительны, но сетка интерполяции неверна. Это потому, что изображение состоит не из точек выборки, а из пикселей. Пиксели - это области, представленные ее средним значением. На самом деле пиксели изображения выглядят так:

    0.5        1        1.5        2        2.5
0.5  +-------------------+-------------------+
     |                   |                   |
     |                   |                   |
     |                   |                   |
 1   |         o         |         o         |
     |                   |                   |
     |                   |                   |
     |                   |                   |
1.5  +-------------------+-------------------+
     |                   |                   |
     |                   |                   |
     |                   |                   |
 2   |         o         |         o         |
     |                   |                   |
     |                   |                   |
     |                   |                   |
2.5  +-------------------+-------------------+

Таким образом, площадь верхнего левого пикселя составляет [0,5, 1,5] x [0,5, 1,5] с центром в (1, 1). И то, что вы хотите, увеличивая масштаб до 2, - это следующие новые пиксели (в координатном пространстве старой сетки, поскольку изображение покрывает все ту же область):

    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  +---------+---------+---------+---------+

Теперь вы берете новые центры x как сетку интерполяции, а старые центры o как сетку образца. Вы видите, что новые пиксели границы действительно нуждаются в экстраполяции. Мы предполагаем, что он экстраполирует константу, поэтому мы можем дополнить массив, чтобы снова выполнить интерполяцию или ограничить сетку интерполяции. В качестве кода с использованием interp2 вы можете:

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)

Вот более общее решение (действительно только для целочисленных размеров вывода), основанное на комментариях Амро ниже его ответ. Если вы разрешаете масштабные коэффициенты, которые приводят к выходному размеру с плавающей запятой.

При нецелочисленных выходных размерах количество новых пикселей будет таким, что последний пиксель будет перекрываться. Так, например, с коэффициентом масштабирования 5/4 = 1,25 размер пикселя будет 1 / (5/4) = 4/5 = 0,8. Следовательно, масштабирование изображения 2x2 на 1,25 дает изображение 3x3. Центры старых пикселей (сетка выборки) находятся в точках 1 и 2, а центры новых пикселей (сетка интерполяции) находятся в точках 0,9, 1,7 и 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

Вот код, чтобы показать это:

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

Это все, что касается изменения размера с помощью билинейной интерполяции. Для бикубической интерполяции Matlab использует симметричное заполнение и сверточный алгоритм, который взвешивает окрестности 4x4. Octave ведет себя иначе (патч поступает).

person John    schedule 11.04.2021
comment
Спасибо, что нашли время, чтобы поделиться этим, очень полезно (я буду ссылаться на него всякий раз, когда забуду все эти детали)! Чтение раздела интерполяции образцов и последней части напомнило мне предыдущий вопрос, в котором мы рассмотрели бикубическую интерполяцию и разницу в реализации. между MATLAB и OpenCV (в этом случае это была константа a=-0.5 vs a=-0.75, и, видя вышеуказанный отчет об ошибке, кажется, что Octave также выбрано другое постоянное значение) - person Amro; 12.04.2021
comment
@Amro Я видел этот вопрос, и есть даже похожие, где люди думают, что Matlab использует a = -1, но a=-0.5 правильно. Octave вообще не использует алгоритм сверточной бикубической интерполяции. Таким образом, нет a, чтобы определить, чтобы получить метод Octave. В настоящее время исправление предназначено только для imresize и imrotate. Может быть, изменения можно перенести на interp2, но, думаю, даже отчета об ошибке нет. - person John; 12.04.2021

В следующем коде показано, как выполнить билинейную интерполяцию с использованием 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 (полностью реализована в m-коде). Насколько я понял, он выполняет какое-то аффинное преобразование изображения. Возможно, кто-нибудь, более знакомый с этой техникой, сможет пролить свет на эту тему ...


ИЗМЕНИТЬ (2021 г.)

См. ответ Джона для более подробного объяснения.

person Amro    schedule 22.06.2011
comment
Я предполагаю, что Matlab запускает фильтр Гаусса для входных данных перед повышающей дискретизацией. Octave этого не делает (хотя octave упоминает это в справке imresize, но не делает этого, когда я проверял. - person goldenmean; 23.06.2011
comment
@goldenmean: это объясняет разницу. В любом случае, чтобы ответить на ваш исходный вопрос, статья в Википедии, на которую я ссылаюсь, хорошо объясняет процесс билинейной интерполяции на примере. - person Amro; 23.06.2011
comment
Я еще раз взглянул на функцию IMRESIZE, прочтите мой ответ на этот другой вопрос: stackoverflow.com/questions/7758078/ - person Amro; 14.10.2011
comment
Версия imresize в Octave делает это неправильно, а imresize в MATLAB делает это правильно. Ключ в том, что изображение 2x2 имеет верхний левый угол в (0,5, 0,5) и нижний правый угол в (2,5, 2,5). Эта область разделена на 2 на 2 пикселя. Затем вы хотите расширить его до 4 на 4 пикселя, что означает, что новые центры пикселей находятся в 0,75, 1,25, 1,75 и 2,25 с точки зрения старых координат. Это то, что вам нужно в xi и yi. См. Также этот отчет об ошибке: savannah.gnu.org/bugs/?51769 - person John; 01.04.2021
comment
@ Амро Да, именно так. Я немного упростил ваш скрипт и добавил отступы, а также imresize для сравнения, см. фрагмент здесь. Обратите внимание, что imresize в Octave online в настоящее время неверен, исправление исправлено отправлено только вчера. - person John; 02.04.2021
comment
@John nice, рассмотрите возможность публикации своего кода + объяснения в качестве нового ответа, у вас определенно есть мой голос :) - person Amro; 05.04.2021
comment
@Amro Хорошо, я добавил ответ с множеством объяснений и пиксельной графикой. Было довольно много работы. Спасибо за обсуждение здесь! - person John; 11.04.2021

Я адаптировал функцию 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;
    }
}
person lhlmgr    schedule 14.09.2011
comment
Так вы говорите, что ваш код эквивалентен коду MATLAB imresize? Вы каким-то образом перепроектировали скомпилированную MEX функцию, о которой упоминал Amro? Или вы просто реализовали аналогичную функцию, не проверяя, что она дает одинаковые результаты в разных угловых случаях? - person Jonas Heidelberg; 03.02.2012
comment
Хм .. эквивалент .. :) Я пытался понять функцию imresize в MATLAB и перенести ее на JAVA. (Исходный код этой функции доступен для чтения). С помощью этой функции я получил (для матриц 600x600, размер которых был изменен до 192x192) те же результаты, что и в MATLAB. Но я не тестировал его юнит-тестами. - person lhlmgr; 03.02.2012