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

40 ビュー (過去 30 日間)
Royi Avital
Royi Avital 2019 年 1 月 15 日
編集済み: Matt J 2019 年 1 月 22 日
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.
  3 件のコメント
Royi Avital
Royi Avital 2019 年 1 月 16 日
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).

サインインしてコメントする。

採用された回答

Matt J
Matt J 2019 年 1 月 15 日
編集済み: Matt J 2019 年 1 月 15 日
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
  12 件のコメント
Royi Avital
Royi Avital 2019 年 1 月 17 日
編集済み: Royi Avital 2019 年 1 月 17 日
By the way, instead of mK = cell2mat(cColumn); one could write mK = cat(2, cColumn{1, :});.

サインインしてコメントする。

その他の回答 (3 件)

Matt J
Matt J 2019 年 1 月 15 日
編集済み: Matt J 2019 年 1 月 15 日
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

Matt J
Matt J 2019 年 1 月 15 日
編集済み: Matt J 2019 年 1 月 15 日
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
  4 件のコメント
Royi Avital
Royi Avital 2019 年 1 月 16 日
If I understand it correctly what you do here is building the Impulse Response just in 3rd dimension. Is that correct?

サインインしてコメントする。


Royi Avital
Royi Avital 2019 年 1 月 19 日
編集済み: Royi Avital 2019 年 1 月 22 日
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.
  3 件のコメント
Matt J
Matt J 2019 年 1 月 22 日
編集済み: Matt J 2019 年 1 月 22 日
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.

サインインしてコメントする。

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by