Creating Convolution Matrix of 2D Kernel for Different Shapes of Convolution

MATLAB, for thos who have access to Image Processing Toolbox offers the function convmtx2().

Yet there are 2 issues:

• It is only available to those who purchased Image PRocessing Toolbox.
• It creates the matrix for full convolution shape only.

I implemented the matrix form for imfiter() in Generate the Matrix Form of 2D Convolution Kernel. It was written in simple form (No vectorization tricks) for clarity and simplicity for thos who want to learn.

What I'm after is doing somthing similar for Convolution Matrices for the different shapes: full, same, valid.

Namely a function with the following form:

function [ mK ] = CreateImageConvMtx( mH, numRows, numCols, convShape )
%UNTITLED6 Summary of this function goes here
%   Detailed explanation goes here

CONVOLUTION_SHAPE_FULL         = 1;
CONVOLUTION_SHAPE_SAME         = 2;
CONVOLUTION_SHAPE_VALID        = 3;

switch(convShape)
    case(CONVOLUTION_SHAPE_FULL)
        % Code for the 'full' case
    case(CONVOLUTION_SHAPE_SAME)
        % Code for the 'same' case
    case(CONVOLUTION_SHAPE_VALID)
        % Code for the 'valid' case
end


end

I would be happy of someone could assist with that.

Again, prefer clarity over performance.

Thank You.

[Royi Avital]

Added Cody - https://www.mathworks.com/matlabcentral/cody/problems/44832.

[Image Analyst]

I guess I don't understand. If one of the disadvantages is that it requires the Image Processing Toolbox, then how/why are you using imfilter(), which also requires the same toolbox? So if you have the toolbox, why do you care about that?

I have never used convmtx2(). Why would I need to? I don't see any advantage to using that over conv2()m, especially if you say it doesn't have 'same', 'full', and 'valid' options.

If you're trying to come up with a kernel to filter your 2-D array with, why don't you just create it yourself, either with fspecial() or just manually hard-code in all the weights?

[Royi Avital]

It's not for practical use as applying convolution using Matrix isn't efficient (Unless done for multiple images at once as done in Deep Learning).

It is for educational use hence it should stand on its own.

My other code doesn't use imfilter() it just build imfilter in matrix form (Without using it).

Answer 1

[Matt J]

Readability and clarity over performance.

OK, can't get much simpler and more readable then the following:

function [ mK ] = CreateImageConvMtx( mH, nRows, nCols, convShape )
%convShape is 'full', 'same', or 'valid'

   impulse=zeros(nRows,nCols);
   
   for i=numel(impulse):-1:1
  
         impulse(i)=1;  %Create impulse image corresponding to i-th output matrix column
         
         tmp=sparse( conv2(impulse,mH,convShape) );  %impulse response
   
         Column{i}=tmp(:);
         
         impulse(i)=0;
   end

   mK=cell2mat(Column);
   
end
   

[Royi Avital]

I think you had a typo with Z.

[Matt J]

Fixed.

[Royi Avital]

I voted this up (Validated it on https://www.mathworks.com/matlabcentral/cody/problems/44832-generate-convolution-matrix-of-2d-kernel-with-different-convolution-shapes-full-same-valid/solutions/1707450).

Is there a solution which is more C style like?

Without using the Cell2Mat? Just building the Sparse Matrix directly?

Thank You.

[Matt J]

Is there a solution which is more C style like?...Just building the Sparse Matrix directly?

I think you mean something like below, but it's a really unrealistic approach, whether in C or anything else. It is horribly inefficient to build a sparse matrix by incrementally loading in the non-zero entries.

function [ mK ] = CreateImageConvMtx( mH, nRows, nCols, convShape )
%convShape is 'full', 'same', or 'valid'
   
   impulse=zeros(nRows,nCols);
   tmp=conv2(impulse,mH,convShape);
   
   Kcol=numel(impulse);
   Krow=numel(tmp);
   
   
   mK=spalloc(Krow,Kcol,numel(mH)*Kcol);
   
   for k=Kcol:-1:1
  
         impulse(k)=1;  %Create impulse image corresponding to i-th output matrix column
         
         
         tmp=conv2(impulse,mH,convShape); %impulse response
         
         [I,~,S]=find(tmp(:));
         
         mK(I,k)=S;
     
         impulse(k)=0; %reset
   end   
   
   
end

[Royi Avital]

What about building the Rows / Coloumns index vecotrs and the Vals vector?

[Matt J]

That would look as follows,

function [ mK ] = CreateImageConvMtx( mH, nRows, nCols, convShape )
%convShape is 'full', 'same', or 'valid'
   
   impulse=zeros(nRows,nCols); 
   
   Kcols=numel(impulse);
   
   for k=Kcols:-1:1
  
         impulse(k)=1;  %Create impulse image corresponding to i-th output matrix column
         
         
         tmp=conv2(impulse,mH,convShape); %impulse response
         
         [I{k},~,Vals{k}]=find(tmp(:));
         
         J{k}=I{k}; J{k}(:)=k;
     
         impulse(k)=0; %reset
   end   
   
   I=vertcat(I{:}); 
   J=vertcat(J{:});
   Vals=vertcat(Vals{:});
   Krows=numel(tmp);
   
      
   mK=sparse(  I, J, Vals, Krows, Kcols );
   
end

[Royi Avital]

Matt, Would you build the matrix in that way in C?

I thought there would be a way to generate the Rows / Cols and Vals in a direct way (Without doing the nice trick of impulse response). Something that will scale for large images as well.

[Matt J]

That's what my proposed solution with ndSparse does, but IMO it defeats the educational priorities you say you are pursuing here. It is much clearer what the computations are doing, and the code is much more compact, when expressed as a series of convolutions with impulses.

[Royi Avital]

I do agree the impulse response is very clear.

Yet it is very inefficient. Is there a way to vectorize the build of the I, J and Vals without using too exotic operations?

The idea is the Impulse Solution is Step I. Now in Step II we can do something more efficient. The ndSparse is something like Step V.

[Matt J]

Yes, there are ways to generate the (i,j,vals) data efficiently in a stand-alone piece of code. However, since we're now shifting priority to efficiency, the code would not look very C-like.

Could I ask who the audience of this code is going to be? Is it for a programming course that you are teaching? If the idea is to show how this would look in C, I'm having trouble seeing why you are pursuing this in Matlab at all.

[Royi Avital]

Hi Matt, First, thank you fo the effort.

Let me try explain. I would like to be able to explain in asimple way how to build the convolution matrix. Indeed the Impulse Response is probably the simplest way out there.

The full case is usually handled by doubly block toeplitz matrix. Yet I couldn't find a code which shows it clearly (I guess this is what convmtx2() does but the code isn't clear, at leats in my opinion).

Next step would be being more efficinent. I think it can diverge into 2 ways:

• C Style - The code is element wise to construct the I, J, Vals vectors. Preferebly the loop will show the inner structure of the large matrix.
• MATLAB Style - The code is more vectorized yet it is still clear what is being done in order to construct the matrix.

What do you think? How would you show the steps so others could be able to fully understand the code, the structure, what's beiong done and later even implement it, in C let's say.

[Royi Avital]

By the way, instead of mK = cell2mat(cColumn); one could write mK = cat(2, cColumn{1, :});.

Answer 2

[Matt J]

You can use my func2mat (Download) utility. It will find the matrix form of any linear function and doesn't require any toolboxes. However, there are much better options if your convolution kernels are separable.

function [ mK ] = CreateImageConvMtx( mH, nRows, nCols, convShape )
%convShape is 'full', 'same', or 'valid'

  Xtypical=zeros([nRows,nCols]);
  
  fun=@(x) conv2(x,mH,convShape);
  
  mK=func2mat(fun,Xtypical);
  
end

[Royi Avital]

@Matt, any thoughts on the code I posted at - https://www.mathworks.com/matlabcentral/answers/439928-creating-convolution-matrix-of-2d-kernel-for-different-shapes-of-convolution#answer_357161?

Answer 3

[Matt J]

You could also use the convn method of my ndSparse class (Download).

function [ mK ] = CreateImageConvMtx( mH, nRows, nCols, convShape )
%convShape is 'full', 'same', or 'valid'

  E=ndSparse(speye(nRows*nCols), [nRows,nCols,nRows,nCols]);
  
  mK=convn(E,mH, convShape);

  mK=sparse2d( reshape(mK,   nRows*nCols, []  ) ) ;

end

[Royi Avital]

Hi Matt,

I'm looking for a function which is stand alone (MATLAB Base only).

Also Id appericiate code, as in my code for `imfilter()` which is simple to grasp (For others to look at). Readability and clarity over performance.

Thank You.

[Matt J]

Yes, both of my solutions run in Matlab Base (no toolboxes required).

[Royi Avital]

But they require additional functions. I want it to be a stand alone solution within a single function.

Also if someone wants to learn how build this matrix, it will be really hard form those functions.

[Royi Avital]

If I understand it correctly what you do here is building the Impulse Response just in 3rd dimension. Is that correct?

Answer 4

[Royi Avital]

Here is my solution which build Doubly Block Toeplitz Matrix:

function [ mK ] = Create2DKernelConvMtxSparse( mH, numRows, numCols, convShape )

CONVOLUTION_SHAPE_FULL  = 1;
CONVOLUTION_SHAPE_SAME  = 2;
CONVOLUTION_SHAPE_VALID = 3;

numColsKernel   = size(mH, 2);
numBlockMtx     = numColsKernel;

cBlockMtx = cell(numBlockMtx, 1);

for ii = 1:numBlockMtx
    cBlockMtx{ii} = CreateConvMtxSparse(mH(:, ii), numRows, convShape);
end

switch(convShape)
    case(CONVOLUTION_SHAPE_FULL)
        diagIdx     = 0;
        numRowsKron = numCols + numColsKernel - 1;
    case(CONVOLUTION_SHAPE_SAME)
        diagIdx     = floor(numColsKernel / 2);
        numRowsKron = numCols;
    case(CONVOLUTION_SHAPE_VALID)
        diagIdx     = numColsKernel - 1;
        numRowsKron = numCols - numColsKernel + 1;
end

vI = ones(min(numRowsKron, numCols), 1);

mK = kron(spdiags(vI, diagIdx, numRowsKron, numCols), cBlockMtx{1});
for ii = 2:numBlockMtx
    diagIdx = diagIdx - 1;
    mK = mK + kron(spdiags(vI, diagIdx, numRowsKron, numCols), cBlockMtx{ii});
end


end

It relies on a function called CreateConvMtxSparse which creates a convolution matrix for 1D Kernel which goes:

function [ mK ] = CreateConvMtxSparse( vK, numElements, convShape )

CONVOLUTION_SHAPE_FULL         = 1;
CONVOLUTION_SHAPE_SAME         = 2;
CONVOLUTION_SHAPE_VALID        = 3;

kernelLength    = length(vK);

switch(convShape)
    case(CONVOLUTION_SHAPE_FULL)
        rowIdxFirst = 1;
        rowIdxLast  = numElements + kernelLength - 1;
        outputSize  = numElements + kernelLength - 1;
    case(CONVOLUTION_SHAPE_SAME)
        rowIdxFirst = 1 + floor(kernelLength / 2);
        rowIdxLast  = rowIdxFirst + numElements - 1;
        outputSize  = numElements;
    case(CONVOLUTION_SHAPE_VALID)
        rowIdxFirst = kernelLength;
        rowIdxLast  = (numElements + kernelLength - 1) - kernelLength + 1;
        outputSize  = numElements - kernelLength + 1;
end

mtxIdx = 0;

% The sparse matrix constructor ignores valus of zero yet the Row / Column
% indices must be valid indices (Positive integers). Hence 'vI' and 'vJ'
% are initialized to 1 yet for invalid indices 'vV' will be 0 hence it has
% no effect.
vI = ones(numElements * kernelLength, 1);
vJ = ones(numElements * kernelLength, 1);
vV = zeros(numElements * kernelLength, 1);

for jj = 1:numElements
    for ii = 1:kernelLength
        if((ii + jj - 1 >= rowIdxFirst) && (ii + jj - 1 <= rowIdxLast))
            % Valid otuput matrix row index
            mtxIdx = mtxIdx + 1;
            vI(mtxIdx) = ii + jj - rowIdxFirst;
            vJ(mtxIdx) = jj;
            vV(mtxIdx) = vK(ii);
        end
    end
end

mK = sparse(vI, vJ, vV, outputSize, numElements);


end

The above pass the test (As the code by Matt) in Cody.

I wonder if there is more efficient way to create the matrices while keeping the main idea of creating the Doubly Block Toeplitz Matrix. Maybe even doing vI, vJ and vV directly to the 2D convolution matrix.

[Matt J]

@Matt, any thoughts on the code I posted at -

Overall, I think it's a fine blend of simplicity and efficiency (so +1!), but as I still don't know who this code is for, it's hard for me to say whether it's the best you can do.

I wonder if there is more efficient way to create the matrices while keeping the main idea of creating the Doubly Block Toeplitz Matrix.

I would consider this:

function [ mK ] = CreateConvMtxSparse( vK, numElements, convShape )
%conShape is 'full','same', or'valid'

  mK = sparse( conv2(eye(numElements),vk,convShape) );

[Royi Avital]

Hi Matt,

The problem with mK = sparse( conv2(eye(numElements),vk,convShape) ); that in case numElements is a large number this becomes inefficient both computationaly and memory wise.

You can see issues I'm having with Improving Speed and Reducing Memory Consumption with Creation of 2D Sparse Convolution Matrix.

Regaridng the code audience, What I'm looking for is a code which makes sense Math wise (Namely you can see what it does matches the math of the case) while being efficient, not tricky and if someone want to reprodcue it in C or any other low level language he can do that easily.

[Matt J]

that in case numElements is a large number this becomes inefficient both computationaly and memory wise.

Less efficient, yes, but clearer (are you still prioritizing clarity over performance?), and still very fast even for absurdly large image sizes.

N=5000; %image length
K=7; %kernel length

vK=rand(K,1);
E=eye(N);

tic
mK1 = CreateConvMtxSparse( vK, N, 3 );
toc; %Elapsed time is 0.002846 seconds.

tic; 
mK2=sparse(conv2(E,vK,'valid')); 
toc %Elapsed time is 0.275314 seconds.