Generate the Matrix Form of 2D Convolution Kernel

Question

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? — geometrikal, Jul 21 '14 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 — Roi Divon, Jul 21 '14 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. — geometrikal, Jul 21 '14 at 21:15
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 '21 at 09:43
@RoiDivon, The answer got 6 up votes. It seems to do what you asked for (Which is very explicit). It will be appreciated if you marked it. — Royi, Mar 14 '22 at 20:12

Royi · Answer 1 · 2022-04-15T11:31:20.493

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

            if((ii + kk &lt;= numRows) &amp;&amp; (ii + kk &gt;= 1) &amp;&amp; (jj + ll &lt;= numCols) &amp;&amp; (jj + ll &gt;= 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;
            % Pixel Index Shift such that pxIdx + pxShift is the linear
            % index of the pixel in the image
            pxShift = (ll * numRows) + kk;

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

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

            if(jj + ll &gt; numCols)
                pxShift = pxShift - ((2 * (jj + ll - numCols) - 1) * numRows);
            end

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

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

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

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

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

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

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

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

            if(jj + ll &lt; 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 - https://stackoverflow.com/questions/2080835. — Royi, Jan 15 '19 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? — stander Qiu, Dec 29 '19 at 09: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 '19 at 09: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. — stander Qiu, Dec 30 '19 at 02:46
@standerQiu, Could you give me an example code to reproduce the issue? — Royi, Dec 30 '19 at 05:01
Just use h = randn(20,20) to produce the kernel, and output is other parameters are 50,50,1,4. — stander Qiu, Dec 30 '19 at 05:15
I can see my code assumes the kernel to have odd dimensions. It doesn't support even Height / Width for the kernel. Simple solution it to add row and columns of zeros to your kernel to the anchor will be in the middle. — Royi, Dec 30 '19 at 20:42
@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 '19 at 21:46

Generate the Matrix Form of 2D Convolution Kernel

1 Answers1

Linked