Как свертка может быть выражена как матричное умножение (матричная форма)?

11

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

Если я правильно понял, фильтры Гаусса свернуты с изображением для уменьшения шума, так как они вычисляют средневзвешенное значение окрестности пикселя, и они очень полезны для обнаружения краев, поскольку вы можете применить размытие и получить изображение одновременно просто свертывается с производной гауссовской функции.

Но может ли кто-нибудь объяснить мне или дать некоторые рекомендации о том, как они рассчитываются?

Например , детектор краев Канни говорит о фильтре Гаусса 5x5, но как они получили эти конкретные числа? И как они прошли путь от непрерывной свертки до умножения матриц?

image-processing filtering gradient edge-detection Matteo
источник
dsp.stackexchange.com/questions/2969/…
Андрей Рубштейн
Я добавил ответ с полным кодом для генерации матрицы для Image Convolution.
Рой

Ответы:

3

Чтобы эта операция работала, вам нужно представить, что ваше изображение преобразовано в вектор. Затем этот вектор умножается слева на матрицу свертки, чтобы получить размытое изображение. Обратите внимание, что результатом также является вектор того же размера, что и входные данные, то есть изображение того же размера.

Каждая строка матрицы свертки соответствует одному пикселю во входном изображении. Он содержит вес вкладов всех других пикселей в изображении в размытый аналог рассматриваемого пикселя.

Давайте рассмотрим пример: размытие рамки размером пикселя на изображении размером пикселей. Измененное изображение представляет собой столбец из 36 элементов, а матрица размытия имеет размер .3×36×636×36

  • Давайте установим эту матрицу в 0 везде.
  • Теперь рассмотрим пиксель координат во входном изображении (для простоты не на его границе). Его размытый аналог получается путем приложения веса к себе и к каждому из его соседей в позициях .(i,j)1/9(i−1,j−1);(i−1,j),(i−1,j+1),…,(i+1,j+1)
  • В векторе столбца пиксель имеет позицию (при условии лексикографического упорядочения). мы сообщаем вес в -ой строке матрицы размытия.6 * я + J 1 / 9 ( 6 я + J )(i,j)6∗i+j1/9(6i+j)
  • Сделайте то же самое со всеми остальными пикселями.

Визуальную иллюстрацию тесно связанного процесса (свертка + вычитание) можно найти в этом посте (из моего личного блога).

sansuiso
источник
ссылка мертва
Гаут
2

Для приложений к изображениям или сверточным сетям, чтобы более эффективно использовать матричные умножители в современных графических процессорах, входные данные обычно преобразуются в столбцы матрицы активации, которые затем можно умножать на несколько фильтров / ядер одновременно.

Проверьте эту ссылку из CS231n Стэнфорда и прокрутите вниз до раздела «Реализация как матричное умножение» для получения подробной информации.

Процесс работает, принимая все локальные исправления на входном изображении или карте активации, которые будут умножены на ядро, и растягивая их в столбец новой матрицы X с помощью операции, обычно называемой im2col. Ядра также растягиваются, чтобы заполнить строки весовой матрицы W, чтобы при выполнении матричной операции W * X полученная матрица Y имела все результаты свертки. Наконец, матрица Y должна быть снова изменена путем преобразования столбцов обратно в изображения с помощью операции, обычно называемой cal2im.

Родриго
источник
1
Это очень хорошая ссылка, спасибо! Однако рекомендуется добавлять в ответ важные выдержки из ссылки, таким образом, ответ действителен, даже если ссылка не работает. Пожалуйста, рассмотрите возможность редактирования своего ответа, чтобы он был принят!
Маттео
1

Свертка во временной области равна матричному умножению в частотной области и наоборот.

Фильтрация эквивалентна свертке во временной области и, следовательно, умножению матриц в частотной области.

Что касается карт 5х5 или масок, то они исходят из дискретности операторов canny / sobel.

Naresh
источник
2
Я не согласен с тем, что фильтрация - это свертка в частотной области. Типы фильтров, о которых мы здесь говорим, - это свертки в пространственной области (то есть поэлементное умножение на отклик фильтра в частотной области).
pichenettes
Спасибо за исправление того, что я написал. Я сделал последующее редактирование. Я думаю, что я должен перепроверить мои ответы перед публикацией. Однако большая часть моего ответа остается в силе.
Нареш
Преобразование Фурье действительно превращает свертки в умножения (и наоборот). Тем не менее, они представляют собой умножения в пинте, тогда как вопрос заключается в умножении матрицы на вектор, которое получается путем изменения формы изображений.
Sansuiso
Я упоминал, что дискретизация операторов является причиной для матриц 5x5, полученных для операторов canny / sobel.
Нареш
1

Я написал функцию, которая решает эту проблему в моем репозитории StackOverflow Q2080835 GitHub (посмотрите CreateImageConvMtx()).
На самом деле эта функция может поддерживать любую форму свертки вы хотите - full, sameи valid.

Код выглядит следующим образом:

function [ mK ] = CreateImageConvMtx( mH, numRows, numCols, convShape )

CONVOLUTION_SHAPE_FULL  = 1;
CONVOLUTION_SHAPE_SAME  = 2;
CONVOLUTION_SHAPE_VALID = 3;

switch(convShape)
    case(CONVOLUTION_SHAPE_FULL)
        % Code for the 'full' case
        convShapeString = 'full';
    case(CONVOLUTION_SHAPE_SAME)
        % Code for the 'same' case
        convShapeString = 'same';
    case(CONVOLUTION_SHAPE_VALID)
        % Code for the 'valid' case
        convShapeString = 'valid';
end

mImpulse = zeros(numRows, numCols);

for ii = numel(mImpulse):-1:1
    mImpulse(ii)    = 1; %<! Create impulse image corresponding to i-th output matrix column
    mTmp            = sparse(conv2(mImpulse, mH, convShapeString)); %<! The impulse response
    cColumn{ii}     = mTmp(:);
    mImpulse(ii)    = 0;
end

mK = cell2mat(cColumn);


end

Я также создал функцию для создания матрицы для фильтрации изображений (идеи, аналогичные MATLAB imfilter()):

function [ mK ] = CreateImageFilterMtx( mH, numRows, numCols, operationMode, boundaryMode )
%UNTITLED6 Summary of this function goes here
%   Detailed explanation goes here

OPERATION_MODE_CONVOLUTION = 1;
OPERATION_MODE_CORRELATION = 2;

BOUNDARY_MODE_ZEROS         = 1;
BOUNDARY_MODE_SYMMETRIC     = 2;
BOUNDARY_MODE_REPLICATE     = 3;
BOUNDARY_MODE_CIRCULAR      = 4;

switch(operationMode)
    case(OPERATION_MODE_CONVOLUTION)
        mH = mH(end:-1:1, end:-1:1);
    case(OPERATION_MODE_CORRELATION)
        % mH = mH; %<! Default Code is correlation
end

switch(boundaryMode)
    case(BOUNDARY_MODE_ZEROS)
        mK = CreateConvMtxZeros(mH, numRows, numCols);
    case(BOUNDARY_MODE_SYMMETRIC)
        mK = CreateConvMtxSymmetric(mH, numRows, numCols);
    case(BOUNDARY_MODE_REPLICATE)
        mK = CreateConvMtxReplicate(mH, numRows, numCols);
    case(BOUNDARY_MODE_CIRCULAR)
        mK = CreateConvMtxCircular(mH, numRows, numCols);
end


end


function [ mK ] = CreateConvMtxZeros( mH, numRows, numCols )
%UNTITLED6 Summary of this function goes here
%   Detailed explanation goes here

numElementsImage    = numRows * numCols;
numRowsKernel       = size(mH, 1);
numColsKernel       = size(mH, 2);
numElementsKernel   = numRowsKernel * numColsKernel;

vRows = reshape(repmat(1:numElementsImage, numElementsKernel, 1), numElementsImage * numElementsKernel, 1);
vCols = zeros(numElementsImage * numElementsKernel, 1);
vVals = zeros(numElementsImage * numElementsKernel, 1);

kernelRadiusV = floor(numRowsKernel / 2);
kernelRadiusH = floor(numColsKernel / 2);

pxIdx       = 0;
elmntIdx    = 0;

for jj = 1:numCols
    for ii = 1:numRows
        pxIdx = pxIdx + 1;
        for ll = -kernelRadiusH:kernelRadiusH
            for kk = -kernelRadiusV:kernelRadiusV
                elmntIdx = elmntIdx + 1;

                pxShift = (ll * numCols) + kk;

                if((ii + kk <= numRows) && (ii + kk >= 1) && (jj + ll <= numCols) && (jj + ll >= 1))
                    vCols(elmntIdx) = pxIdx + pxShift;
                    vVals(elmntIdx) = mH(kk + kernelRadiusV + 1, ll + kernelRadiusH + 1);
                else
                    vCols(elmntIdx) = pxIdx;
                    vVals(elmntIdx) = 0; % See the accumulation property of 'sparse()'.
                end
            end
        end
    end
end

mK = sparse(vRows, vCols, vVals, numElementsImage, numElementsImage);


end


function [ mK ] = CreateConvMtxSymmetric( mH, numRows, numCols )
%UNTITLED6 Summary of this function goes here
%   Detailed explanation goes here

numElementsImage    = numRows * numCols;
numRowsKernel       = size(mH, 1);
numColsKernel       = size(mH, 2);
numElementsKernel   = numRowsKernel * numColsKernel;

vRows = reshape(repmat(1:numElementsImage, numElementsKernel, 1), numElementsImage * numElementsKernel, 1);
vCols = zeros(numElementsImage * numElementsKernel, 1);
vVals = zeros(numElementsImage * numElementsKernel, 1);

kernelRadiusV = floor(numRowsKernel / 2);
kernelRadiusH = floor(numColsKernel / 2);

pxIdx       = 0;
elmntIdx    = 0;

for jj = 1:numCols
    for ii = 1:numRows
        pxIdx = pxIdx + 1;
        for ll = -kernelRadiusH:kernelRadiusH
            for kk = -kernelRadiusV:kernelRadiusV
                elmntIdx = elmntIdx + 1;

                pxShift = (ll * numCols) + kk;

                if(ii + kk > numRows)
                    pxShift = pxShift - (2 * (ii + kk - numRows) - 1);
                end

                if(ii + kk < 1)
                    pxShift = pxShift + (2 * (1 -(ii + kk)) - 1);
                end

                if(jj + ll > numCols)
                    pxShift = pxShift - ((2 * (jj + ll - numCols) - 1) * numCols);
                end

                if(jj + ll < 1)
                    pxShift = pxShift + ((2 * (1 - (jj + ll)) - 1) * numCols);
                end

                vCols(elmntIdx) = pxIdx + pxShift;
                vVals(elmntIdx) = mH(kk + kernelRadiusV + 1, ll + kernelRadiusH + 1);

            end
        end
    end
end

mK = sparse(vRows, vCols, vVals, numElementsImage, numElementsImage);


end


function [ mK ] = CreateConvMtxReplicate( mH, numRows, numCols )
%UNTITLED6 Summary of this function goes here
%   Detailed explanation goes here

numElementsImage    = numRows * numCols;
numRowsKernel       = size(mH, 1);
numColsKernel       = size(mH, 2);
numElementsKernel   = numRowsKernel * numColsKernel;

vRows = reshape(repmat(1:numElementsImage, numElementsKernel, 1), numElementsImage * numElementsKernel, 1);
vCols = zeros(numElementsImage * numElementsKernel, 1);
vVals = zeros(numElementsImage * numElementsKernel, 1);

kernelRadiusV = floor(numRowsKernel / 2);
kernelRadiusH = floor(numColsKernel / 2);

pxIdx       = 0;
elmntIdx    = 0;

for jj = 1:numCols
    for ii = 1:numRows
        pxIdx = pxIdx + 1;
        for ll = -kernelRadiusH:kernelRadiusH
            for kk = -kernelRadiusV:kernelRadiusV
                elmntIdx = elmntIdx + 1;

                pxShift = (ll * numCols) + kk;

                if(ii + kk > numRows)
                    pxShift = pxShift - (ii + kk - numRows);
                end

                if(ii + kk < 1)
                    pxShift = pxShift + (1 -(ii + kk));
                end

                if(jj + ll > numCols)
                    pxShift = pxShift - ((jj + ll - numCols) * numCols);
                end

                if(jj + ll < 1)
                    pxShift = pxShift + ((1 - (jj + ll)) * numCols);
                end

                vCols(elmntIdx) = pxIdx + pxShift;
                vVals(elmntIdx) = mH(kk + kernelRadiusV + 1, ll + kernelRadiusH + 1);

            end
        end
    end
end

mK = sparse(vRows, vCols, vVals, numElementsImage, numElementsImage);


end


function [ mK ] = CreateConvMtxCircular( mH, numRows, numCols )
%UNTITLED6 Summary of this function goes here
%   Detailed explanation goes here

numElementsImage    = numRows * numCols;
numRowsKernel       = size(mH, 1);
numColsKernel       = size(mH, 2);
numElementsKernel   = numRowsKernel * numColsKernel;

vRows = reshape(repmat(1:numElementsImage, numElementsKernel, 1), numElementsImage * numElementsKernel, 1);
vCols = zeros(numElementsImage * numElementsKernel, 1);
vVals = zeros(numElementsImage * numElementsKernel, 1);

kernelRadiusV = floor(numRowsKernel / 2);
kernelRadiusH = floor(numColsKernel / 2);

pxIdx       = 0;
elmntIdx    = 0;

for jj = 1:numCols
    for ii = 1:numRows
        pxIdx = pxIdx + 1;
        for ll = -kernelRadiusH:kernelRadiusH
            for kk = -kernelRadiusV:kernelRadiusV
                elmntIdx = elmntIdx + 1;

                pxShift = (ll * numCols) + kk;

                if(ii + kk > numRows)
                    pxShift = pxShift - numRows;
                end

                if(ii + kk < 1)
                    pxShift = pxShift + numRows;
                end

                if(jj + ll > numCols)
                    pxShift = pxShift - (numCols * numCols);
                end

                if(jj + ll < 1)
                    pxShift = pxShift + (numCols * numCols);
                end

                vCols(elmntIdx) = pxIdx + pxShift;
                vVals(elmntIdx) = mH(kk + kernelRadiusV + 1, ll + kernelRadiusH + 1);

            end
        end
    end
end

mK = sparse(vRows, vCols, vVals, numElementsImage, numElementsImage);


end

Код был проверен на соответствие MATLAB imfilter().

Полный код доступен в моем хранилище StackOverflow Q2080835 GitHub .

Royi
источник