Comment une convolution peut-elle être exprimée comme une multiplication matricielle (forme matricielle)?

Je sais que cette question peut ne pas être très pertinente pour la programmation, mais si je ne comprends pas la théorie derrière le traitement d'image, je ne pourrai jamais mettre en œuvre quelque chose dans la pratique.

Si je comprends bien, les filtres gaussiens sont convolués avec une image pour la réduction du bruit car ils calculent une moyenne pondérée du voisinage d'un pixel et ils sont très utiles pour la détection des contours, car vous pouvez appliquer un flou et dériver l'image en même temps en convoluant simplement avec la dérivée d'une fonction gaussienne.

Mais quelqu'un peut-il m'expliquer ou me donner des références sur la façon dont ils sont calculés?

Par exemple, le détecteur de bord de Canny parle d'un filtre gaussien 5x5, mais comment ont-ils obtenu ces nombres particuliers? Et comment sont-ils passés d'une convolution continue à une multiplication matricielle?

Pour que cette opération fonctionne, vous devez imaginer que votre image est remodelée en tant que vecteur. Ensuite, ce vecteur est multiplié à sa gauche par la matrice de convolution afin d'obtenir l'image floue. Notez que le résultat est également un vecteur de la même taille que l'entrée, c'est-à-dire une image de la même taille.

Chaque ligne de la matrice de convolution correspond à un pixel dans l'image d'entrée. Il contient le poids des contributions de tous les autres pixels de l'image à la contrepartie floue du pixel considéré.

Prenons un exemple: boîte flou de taille pixels sur une image de taille 6 × 6 pixels. L'image remodelée est une colonne de 36 choix, tandis que la matrice de flou a une taille 36 × 36 .$3 \times 3$$6 \times 6$$36 \times 36$

• Initions cette matrice à 0 partout.
• Maintenant, considérez le pixel de coordonnées dans l'image d'entrée (pas sur sa bordure pour plus de simplicité). Son homologue floue est obtenue en appliquant un poids de 1 / neuf à lui - même et chacun de ses voisins au niveau des positions ( i - 1 , j - 1 ) ; ( i - 1 , j ) , ( i - 1 , j + 1 ) , … , ( i + 1 $(i,j)$$1/9$ .$(i-1,j-1); (i-1,j), (i-1,j+1),\ldots,(i+1,j+1)$
• Dans le vecteur colonne, le pixel a la position 6 ∗ i + j (en supposant un ordre lexicographique). Nous rapportons le poids une / 9 dans le ( 6 i + j ) -ième ligne de la matrice de flou.$(i,j)$$6*i + j$$1/9$$(6i+j)$
• Faites de même avec tous les autres pixels.

Une illustration visuelle d'un processus étroitement lié (convolution + soustraction) peut être trouvée sur ce billet de blog (à partir de mon blog personnel).

Pour les applications aux images ou aux réseaux de convolution, pour utiliser plus efficacement les multiplicateurs matriciels dans les GPU modernes, les entrées sont généralement remodelées en colonnes d'une matrice d'activation qui peut ensuite être multipliée avec plusieurs filtres / noyaux à la fois.

Consultez ce lien depuis le CS231n de Stanford et faites défiler jusqu'à la section «Implémentation en tant que multiplication matricielle» pour plus de détails.

Le processus fonctionne en prenant tous les correctifs locaux sur une image d'entrée ou une carte d'activation, ceux qui seraient multipliés avec le noyau, et en les étirant dans une colonne d'une nouvelle matrice X via une opération communément appelée im2col. Les noyaux sont également étirés pour remplir les rangées d'une matrice de poids W de sorte que lors de l'exécution de l'opération de matrice W * X, la matrice résultante Y ait tous les résultats de la convolution. Enfin, la matrice Y doit être remodelée à nouveau en reconvertissant les colonnes en images par une opération généralement appelée cal2im.

La convolution dans le domaine temporel est égale à la multiplication matricielle dans le domaine fréquentiel et vice versa.

Le filtrage est équivalent à la convolution dans le domaine temporel et donc à la multiplication matricielle dans le domaine fréquentiel.

Quant aux cartes ou masques 5x5, ils proviennent de la discrétisation des opérateurs canny / sobel.

J'ai écrit une fonction qui résout cela dans mon référentiel GitHub StackOverflow Q2080835 (Jetez un oeil à CreateImageConvMtx()).
En fait, la fonction peut prendre en charge n'importe quelle forme de convolution que vous souhaitez - full, sameet valid.

Le code est le suivant:

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

J'ai également créé une fonction pour créer une matrice pour le filtrage d'images (idées similaires à 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

Le code a été validé par rapport à MATLAB imfilter().

Le code complet est disponible dans mon référentiel GitHub StackOverflow Q2080835 .