# Generate the Matrix Form of 2D Convolution Kernel

We know that a convolution can be replaced by a multiplication with a Toeplitz / Circulant Matrix. Meaning, assume I have convolution kernel $$h$$ and matrix $$I$$ (Of size $$m \times m$$ for example), then there is a matrix $$H$$ of size $$m^2 \times m^2$$ such that $$h \ast I$$ is the same as $$H I^{cs}$$ Where cs for column stacked image resized to $$m \times m$$ (reshape(H*I(:), [m, m])).

My question is - how to construct $$H$$? I know that every row of $$H$$ is responsible for single entry of the resulting multiplication (H(1,:)*I(:) gives the (1,1) entry of the convolution result), so I can try to work very hard to calculate each entry of $$H$$ (e.g. if the kernel $$h$$ is of size $$3 \times 3$$, than the convolution result in entry (i,j) is:

$$h_{11}\cdot I(i-1,j-1) + h_{12}\cdot I(i-1,j) + h_{13}\cdot I(i-1,j+1) + h_{21}\cdot I(i,j-1) + h_{22}\cdot I(i,j) + h_{23}\cdot I(i,j+1) + h_{31}\cdot I(i+1,j-1) + h_{32}\cdot I(i+1,j) + h_{33}\cdot I(i+1,j+1)$$

except for near the edges of $$I$$. so it's somewhat possible with a lot of hard work to construct $$H$$ such that when multiplying by $$I^{cs}$$ I will get that result). But is there a 'generic' way for this construction? Or perhaps a made MATLAB function?

I am aware of the convmtx2() function in MATLAB, but the resulting matrix is not in the proper dimensions (not even square).

• Whats the application? Depending on the matrix $m^2 \times m^2$ is going to be huge. e.g. 128^2 x 128^2 @ 32 bit is 1 GB. Why not use conv2? Jul 21, 2014 at 12:28
• @geometrikal I'm dealing with matrices of size $64 \times 64$ or less, so things arn't getting huge. I need it to build a wiener filter (and not the DFT version), and for that I need the full matrix representation of the kernel Jul 21, 2014 at 13:10
• The way I would code it is pad the image with half the kernel width, make a matrix of x and y shifts using ndgrid from +- kernel width, loop through the shift matrix and set the third dimension of a size [width I, height I, kernelWidth^2] temporary image with the shifted image, e.g. tempI(:,:,idx) = paddedI(xshift(idx):xshift(idx)+imageWidth-1,yshift(idx):yshift(idx)+imageWidth-1), then make into the column vector using permute to make the third dimension first, and (:) to get the column vector. Jul 21, 2014 at 21:15
• @Royi I honestly don't know... it was 5 years ago :/ Oct 16, 2019 at 14:24
• If imfilter() or / and conv2() does the work for you, this perfectly imitate it (See the unit test). So you can mark it as answered and the question will be helpful to others. Please show appreciation by marking it.
– Royi
May 14, 2021 at 9:43

I created a function to create a Matrix for Image Filtering (Similar ideas to MATLAB's imfilter()):

function [ mK ] = CreateImageFilterMtx( mH, numRows, numCols, operationMode, boundaryMode )
% ----------------------------------------------------------------------------------------------- %
% [ mK ] = CreateImageFilterMtx( mH, numRows, numCols, operationMode, boundaryMode )
% Generates an Image Filtering Matrix for the 2D Kernel (The Matrix mH)
% with support for different operations modes (Convolution / Correlation)
% and boundary conditions (Zeros, Symmetric, Replicate, Circular). The
% function should match the use of MATLAB's 'imfilter()' with the same
% parameters.
% Input:
%   - mH                -   Input 2D Convolution Kernel.
%                           Assumed to have odd dimensions.
%                           Structure: Matrix.
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - numRows           -   Number of Rows.
%                           Number of rows in the output convolution
%                           matrix.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3, ...}.
%   - numCols           -   Number of Columns.
%                           Number of columns in the output convolution
%                           matrix.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3, ...}.
%   - operationMode     -   Operation Mode.
%                           Sets whether to use Convolution or Correlation
%                           for the operation mode.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2}.
%   - boundaryMode      -   Boundary Condition Mode.
%                           Sets the boundary condition mode for the
%                           filtering. The options are - Zeros, Symmetric,
%                           Replicate and Circular.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3, 4}.
% Output:
%   - mK                -   Convolution Matrix.
%                           The output filtering matrix. Multiplying in
%                           the column stack form on an image should be
%                           equivalent to applying 'imfilter()' on the
%                           image.
%                           Structure: Matrix (Sparse).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References:
%   1.  MATLAB's 'convmtx2()' - https://www.mathworks.com/help/images/ref/convmtx2.html.
% Remarks:
%   1.  The height and width of 'mH' are assumed to be odd number. In case
%       either or both are even the user should pad the kernel with zeros
%       (Either a row, column or both). according to the anchor of the
%       kernel the user do the padding pre or post the kernel.
% TODO:
%   1.  Refactor the code to share the common operations of different
%       boundary modes.
% Release Notes:
%   -   1.0.001     30/12/2019  Royi Avital
%       *   Fixed some bugs related to using 'numCols' instead of 'numRows'
%           in the calculation of 'pixelShift' for the cases 'jj + ll >
%           numCols' and 'jj + ll < 1'.
%   -   1.0.000     16/01/2018  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

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);

pxIdx       = 0;
elmntIdx    = 0;

for jj = 1:numCols
for ii = 1:numRows
pxIdx = pxIdx + 1;
elmntIdx = elmntIdx + 1;

% Pixel Index Shift such that pxIdx + pxShift is the linear
% index of the pixel in the image
pxShift = (ll * numRows) + 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);

pxIdx       = 0;
elmntIdx    = 0;

for jj = 1:numCols
for ii = 1:numRows
pxIdx = pxIdx + 1;
elmntIdx = elmntIdx + 1;

% Pixel Index Shift such that pxIdx + pxShift is the linear
% index of the pixel in the image
pxShift = (ll * numRows) + 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) * numRows);
end

if(jj + ll < 1)
pxShift = pxShift + ((2 * (1 - (jj + ll)) - 1) * numRows);
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);

pxIdx       = 0;
elmntIdx    = 0;

for jj = 1:numCols
for ii = 1:numRows
pxIdx = pxIdx + 1;
elmntIdx = elmntIdx + 1;

% Pixel Index Shift such that pxIdx + pxShift is the linear
% index of the pixel in the image
pxShift = (ll * numRows) + 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) * numRows);
end

if(jj + ll < 1)
pxShift = pxShift + ((1 - (jj + ll)) * numRows);
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);

pxIdx       = 0;
elmntIdx    = 0;

for jj = 1:numCols
for ii = 1:numRows
pxIdx = pxIdx + 1;
elmntIdx = elmntIdx + 1;

% Pixel Index Shift such that pxIdx + pxShift is the linear
% index of the pixel in the image
pxShift = (ll * numRows) + kk;

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

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

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

if(jj + ll < 1)
pxShift = pxShift + (numCols * numRows);
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 StackExchnage Codes StackOverflow Q2080835 GitHub Repository (Look at the StackOverflow\Q2080835 folder).

• Related question - stackoverflow.com/questions/2080835.
– Royi
Jan 15, 2019 at 19:21
• Thanks for your answer. I have tried your code and it did generate the same result as imfilter in MATLAB. But I have another question, it seems that the code is to generate the matrix for linear convolution, instead of the circular convolution in discrete domain. Right? Dec 29, 2019 at 9:16
• @standerQiu, I'm happy to hear the code assisted you. Regarding your question, you have the boundaryMode parameter to chose the boundary condition. One of them is Circular Convolution.
– Royi
Dec 29, 2019 at 9:33
• Hi Royi, I find a bug in your code. In the sub-function CreateConvMtxCircular, the vVals(elmntIdx) = mH(kk + kernelRadiusV + 1, ll + kernelRadiusH + 1); would exceed the maximum index. Dec 30, 2019 at 2:46
• @standerQiu, In addition to what I wrote above (Supporting kernel with odd dimensions only and the way to get over it) I updated the code to fix few bugs. I also added Unit Test to verify results. It is indeed now replicate MATLAB's function perfectly.
– Royi
Dec 30, 2019 at 21:46