How Does a Convolution of an Image Can Be Expressed as a Matrix Multiplication (Matrix Form)?

Question

I know this question may not be very relevant to programming, but if I don't understand the theory behind image processing I'll never be able to implement something in practice.

If I got it right Gaussian filters are convolved with an image for noise reduction since they compute a weighed average of a pixel's neighborhood and they are very useful in edge-detection, since you can apply a blur and derive the image at the same time by simply convolving with the derivative of a Gaussian function.

But can anyone explain me, or give me some references on how are they computed?

E.g. Canny's edge detector talks about a 5x5 Gaussian filter, but how did they get those particular numbers? And how did they go from a continuous convolution to a Matrix multiplication?

I added an answer with full code for generating a matrix for Image Convolution. — Royi, Commented Nov 28, 2019 at 6:35
See How can the convolution operation be implemented as a matrix-vector multiplication?. — user40095, Commented Jun 14, 2020 at 14:00
Related: dsp.stackexchange.com/questions/81949, dsp.stackexchange.com/questions/17418, dsp.stackexchange.com/questions/1774. — Royi, Commented Dec 3, 2022 at 15:53

lennon310 · Accepted Answer · 2021-02-20 14:24:44Z

For this operation to work, you need to imagine that your image is reshaped as a vector. Then, this vector is multiplied on its left by the convolution matrix in order to obtain the blurred image. Note that the result is also a vector the same size as the input, i.e., an image of the same size.

Each row of the convolution matrix corresponds to one pixel in the input image. It contains the weight of the contributions of all the other pixels in the image to the blurred counterpart of the considered pixel.

Let's take an example: box blur of size $3 \times 3$ pixels on an image of size $6 \times 6$ pixels. The reshaped image is a column of 36 elects, while the blur matrix has size $36 \times 36$.

Let's initialize this matrix to 0 everywhere.
Now, consider the pixel of coordinates $(i,j)$ in the input image (not on its border for simplicity). Its blurred counterpart is obtained by applying a weight of $1/9$ to itself and each of its neighbours at the positions $(i-1,j-1); (i-1,j), (i-1,j+1),\ldots,(i+1,j+1)$.
In the column vector, the pixel $(i,j)$ has the position $6*i + j$ (assuming lexicographic ordering). we report the weight $1/9$ in the $(6i+j)$-th line of the blur matrix.
Do the same with all other pixels.

A visual illustration of a closely related process (convolution + subtraction) can be found on my personal blog Computers Don't See (Yet) - Linearized model of LBPs.

Really nice blog. I think you should translate the posts from French to English. — Mark, Commented Feb 20, 2021 at 6:45

Naresh · Accepted Answer · 2013-03-17 04:14:42Z

2

Convolution in Time domain equals matrix multiplication in the frequency domain and vice versa.

Filtering is equivalent to convolution in the time domain and hence matrix multiplication in the frequency domain.

As for the 5x5 maps or masks, they come from discretizing the canny/sobel operators.

edited Mar 17, 2013 at 4:14

answered Mar 16, 2013 at 16:30

Naresh

7584 silver badges20 bronze badges

2

$\begingroup$ I don't agree with the fact that filtering is a convolution in the frequency domain. The kind of filters we are talking about here are convolutions in the spatial domain (that is to say, element-wise multiplication by the filter response in the frequency domain). $\endgroup$
– pichenettes
Commented Mar 16, 2013 at 16:51
$\begingroup$ Thanks for correcting what I wrote. I made a subsequent edit. I guess I should double-check my answers before posting. However, the majority of my answer still stands. $\endgroup$
– Naresh
Commented Mar 17, 2013 at 4:15
$\begingroup$ The Fourier transform does indeed turn convolutions into multiplications (and vice versa). However, they are pint wise multiplications, while the question is about matrix-vector multiplications that are obtained by reshaping the images. $\endgroup$
– sansuiso
Commented Mar 17, 2013 at 16:34
$\begingroup$ I did mention how discretizing the operators is the reason for the 5x5 matrices obtained for the canny/sobel operators. $\endgroup$
– Naresh
Commented Mar 17, 2013 at 17:49

Add a comment |

Rodrigo · Accepted Answer · 2017-02-24 06:07:33Z

2

For applications to images or convolution networks, to more efficiently use the matrix multipliers in modern GPUs, the inputs are typically reshaped into columns of an activation matrix that can then be multiplied with multiple filters/kernels at once.

Check out this link from Stanford's CS231n, and scroll down to the section on "Implementation as Matrix Multiplication" for details.

The process works by taking all the local patches on an input image or activation map, the ones that would be multiplied with the kernel, and stretching them into a column of a new matrix X through an operation commonly called im2col. The kernels are also stretched to populate the rows of a weight matrix W so that when performing the matrix operation W*X, the resulting matrix Y has all the results of the convolution. Finally, the Y matrix must be reshaped again by converting the columns back into images by an operation typically called cal2im.

edited Feb 24, 2017 at 6:07

answered Feb 22, 2017 at 22:58

Rodrigo

292 bronze badges

1

$\begingroup$ This is a very good link, thanks! However it is good practice to add the important extracts from the link into the answer, this way the answer is valid even if the link breaks. Please consider editing your answer to get it accepted! $\endgroup$
– Matteo
Commented Feb 23, 2017 at 18:27
$\begingroup$ Using im2ncol() means a lot of memory goes to waste. There are much more efficient methods using Matrix Multiplication. $\endgroup$
– Mark
Commented Feb 20, 2021 at 6:40
$\begingroup$ @Mark. Can you share pointers or descriptions of those more efficient methods and their applicability to GPUs or matrix multipliers? $\endgroup$
– Rodrigo
Commented May 7, 2021 at 8:43
$\begingroup$ Look for Convolution Matrix in this site. There are some answers about this. $\endgroup$
– Mark
Commented May 9, 2021 at 13:45

Add a comment |

lennon310 · Accepted Answer · 2021-02-19 19:03:01Z

1

After looking for the answer to this post for weeks, this is by far the best detailed answer I have found about how to convert a convolution into a matrix multiplication (GEMM) operator.

edited Feb 19, 2021 at 19:03

lennon310

3,59019 gold badges25 silver badges27 bronze badges

answered Feb 19, 2021 at 18:32

Ginés Hidalgo

1112 bronze badges

Add a comment |

Royi · Accepted Answer · 2022-12-02 06:10:57Z

I wrote a function which solves this in my StackOverflow Q2080835 GitHub Repository (Have a look at CreateImageConvMtx()).
Actually the function can support any convolution shape you'd like - full, same and valid.

The idea is simple:

Model the image as a single long vector. In the MATLAB / Fortran / Julia context we assume the memory is contiguous along columns. In Python / C we assume it is along the rows.
Each row of the matrix output a single pixel in the output.
Each row of the matrix basically combines the pixels of interest in the input image.

Now it is just a game of indices to catch the correct pixels represented by the kernel. The code is well documented and clear to follow.

The code is as following:

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

I also created a function to create a Matrix for Image Filtering (Similar ideas to MATLAB's 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

The code was validated against MATLAB imfilter().

Full code is available in my StackOverflow Q2080835 GitHub Repository.

Stack Exchange Network

How Does a Convolution of an Image Can Be Expressed as a Matrix Multiplication (Matrix Form)?

5 Answers 5

Your Answer

Not the answer you're looking for? Browse other questions tagged
image-processing
convolution
filtering
edge-detection
linear-algebra
or ask your own question.

Linked

Hot Network Questions

How Does a Convolution of an Image Can Be Expressed as a Matrix Multiplication (Matrix Form)?

5 Answers 5

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged image-processingconvolutionfilteringedge-detectionlinear-algebra or ask your own question.

Linked

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
image-processing
convolution
filtering
edge-detection
linear-algebra
or ask your own question.