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, Commented 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 — Roi Divon, Commented 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. — geometrikal, Commented Jul 21, 2014 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, Commented May 14, 2021 at 9:43

Royi · Accepted Answer · 2022-04-15 11:31:20Z

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

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

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

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, Commented 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, Commented 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. — stander Qiu, Commented 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, Commented Dec 30, 2019 at 21:46

Stack Exchange Network

Generate the Matrix Form of 2D Convolution Kernel

1 Answer 1

Your Answer

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

Linked

Hot Network Questions

Generate the Matrix Form of 2D Convolution Kernel

1 Answer 1

Your Answer

Sign up or log in

Post as a guest

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

Linked

Related

Hot Network Questions

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