2D convolution of image with filter as successive 1D convolutions

I want to prove (or more precisely experiment with) the idea that a 2D convoltion as produced by the Matlab conv2() function between an image I (2D matrix) and a kernel (smaller 2D matrix) can be implemented as some 1D conv i.e. the Matlab conv() function and NOT conv2(). Of course possibly some reshapes and matrix multiply might be needed but no conv2().

And to make it clear, I am NOT refering to that kind if thing:

s1=[1,0,-1]'
s2=[1 2 1]
diff=conv2(x,y)-conv2(conv2(x,s1),s2)

diff is = 0 everywhere

Rather, I want to do something like

conv(conv(x(:), filter1)filter2) ...

• please? Maybe @Fat32 ? Nov 23 '19 at 15:26
• Anyone please?. Nov 23 '19 at 17:25
• Your code was not implicit. Care should be taken on borders Nov 23 '19 at 17:48
• Math stackexchange may be better. In mathematical terms it is sums of kronecker products of matrices. Nov 23 '19 at 23:00

When a 2D filter $$h[n,m]$$ is separable; i.e., $$h[n,m] = f[n]g[m]$$, then the 2D convolution of an image $$I[n,m]$$ with that filter can be decomposed into 1D convolutions between rows and columns of the image and the 1D filters $$f[n]$$ and $$g[m]$$ respectively.

Let me give you the MATLAB / OCTAVE code, I hope this is what you wanted to show ?

clc; clear all; close all;

N1 = 8;        % input x[n1,n2] row-count
N2 = 5;        % input x[n1,n2] clm-count
M1 = 4;        % impulse response h[n1,n2] row-count
M2 = 3;        % impulse response h[n1,n2] clm-count
L1 = N1+M1-1;  % output row-count
L2 = N2+M2-1;  % output clm-count

x = rand(N1,N2);  % input signal
f = rand(1,M2);   % f[n1] = row vector
g = rand(M1,1);   % g[n1] = column vector
h = g*f;          % h[n1,n2] = f[n1]*g[n2]
y = zeros(L1,L2); % output signal

% S1 - Implement Separable Convolution
% ------------------------------------
for k1 = 1:N2       % I - Convolve COLUMNS of x[:,k] with g[k]
y(:,k1) = conv(x(:,k1),g);   % intermediate output
end

for k2 = 1:L1   % II- Convolve ROWS of yi[k,:] with f[k]
y(k2,:) = conv(y(k2,1:N2),f);
end

% S2 - Matlab conv2() :
% ---------------------
y2 = conv2(x,h);   % check for matlab conv2()

% S3 - Display the Results
% ------------------------
title('The Difference y[n,m] - y2[n,m]');
• thanks for taking time to answer, great! Nov 23 '19 at 18:44
• @Machupicchu ok. But it seems to late :-) Nov 23 '19 at 18:45
• no no its always good to have multiple answers Nov 23 '19 at 18:48
• Yes indeed by Ray Charles Nov 23 '19 at 18:55
• I added a comment here since there is no way to directly send messages to members (which I regret), Dear Fat32, I would be very happy I you had a look at this new question of mine: datascience.stackexchange.com/questions/81923/… Sep 18 '20 at 19:01

If a 2D $$K_2$$ filter kernel is of rank $$0$$ or $$1$$, it can be written as a separable product of $$2$$ 1D kernels $$K_1^r$$ and $$K_1^c$$ on rows and columns. As such, it can implemented by 1D convolutions, as long as one properly reshape the 2D matrices into 1D ones, and take care about "out-of-range" values, to avoid wrap-around. For instance, you can pad in every direction by the size of the filter, and make sure the convolution does not add unwanted information.

Assuming that you know that you have a separable 2D filter, the following code does the job. A one-liner would be:

xRowFull = reshape(conv(reshape(reshape( conv(x(:),s1,'same'),nRow,nCol)',nRow*nCol,1),s2,'same'),nRow,nCol)';

And the code is:

% https://dsp.stackexchange.com/questions/62115/2d-convolution-of-image-with-filter-as-successive-1d-convolutions
%% Initialization
clear all
nRow = 16;
nCol = 16;
HalfSizeCentralImageKernel = 1;
x = zeros(nRow,nCol);
x(nRow/2-HalfSizeCentralImageKernel:nRow/2+HalfSizeCentralImageKernel,nCol/2-HalfSizeCentralImageKernel:nCol/2+HalfSizeCentralImageKernel)=rand(2*HalfSizeCentralImageKernel+1);

%% Original 2D version
s1=[1,0,-1]';
s2=[1 2 1];
y = s1*s2;

%% Step by step 2x1D version
xRowFlat1 = x(:);
xRowFlat1FiltCol = conv(xRowFlat1,s1,'same');
xRowFlat2 = (reshape(xRowFlat1FiltCol,nRow,nCol))';
xRowFlat2 = xRowFlat2(:);
xRowFlat2FiltRowFlat = conv(xRowFlat2,s2,'same');
xRowFlatFilt2Row = reshape(xRowFlat2FiltRowFlat,nRow,nCol)';

%% Compact vectorized 1D version
xRowFull = reshape(conv(reshape(reshape( conv(x(:),s1,'same'),nRow,nCol)',nRow*nCol,1),s2,'same'),nRow,nCol)';

%% Display
figure(1);
imagesc(x);

figure(2);
subplot(1,3,1)
imagesc([conv2(x,y,'same')]); xlabel('Original')
subplot(1,3,2)
imagesc(xRowFlatFilt2Row); xlabel('Separable, step by step')
subplot(1,3,3)
imagesc(xRowFull); xlabel('Separable, one-liner')

diff1=conv2(x,y,'same')-conv2(conv2(x,s1,'same'),s2,'same');
disp(['Max error 1: ',num2str(max(abs(diff1(:))))]);

diff2=conv2(x,y,'same')-xRowFlatFilt2Row;
disp(['Max error 2: ',num2str(max(abs(diff2(:))))]);

Here is a crude Matlab code. Can you test it, and if OK, i'll send a one-liner (if I can).

nRow = 8;
nCol = 8;
HalfSizeCentralKernel = 1;
x = zeros(nRow,nCol);
x(nRow/2-HalfSizeCentralKernel:nRow/2+HalfSizeCentralKernel,nCol/2-HalfSizeCentralKernel:nCol/2+HalfSizeCentralKernel)=rand(2*HalfSizeCentralKernel+1);
figure(1);
imagesc(x);

% 2D version
s1=[1,0,-1]';
s2=[1 2 1];
y = s1*s2;
diff1=conv2(x,y,'same')-conv2(conv2(x,s1,'same'),s2,'same');
disp(['Max error 1: ',num2str(max(abs(diff1(:))))]);

% 1D version
xRowFlat1 = x(:);
xRowFlat1FiltCol = conv(xRowFlat1,s1,'same');
xRowFlat2 = (reshape(xRowFlat1FiltCol,nRow,nCol))';
xRowFlat2 = xRowFlat2(:);
xRowFlat2FiltRow = conv(xRowFlat2,s2,'same');
xRowFlatFilt2Row = reshape(xRowFlat2FiltRow,nRow,nCol)';

figure(2);
subplot(1,2,1)
imagesc([conv2(x,y,'same')])
subplot(1,2,2)
imagesc(xRowFlatFilt2Row)
• thanks it seems to work, fantastic ! Im now trying to understand why you are using the 'same' conv and how you pad etc Nov 23 '19 at 18:42
• fantastic! great Nov 23 '19 at 18:43
• Updated with more details Nov 23 '19 at 20:22
• thanks for the great answers! If you can, maybe have a look to a new question, related to Normalized Cross Correlation (NCC). -> dsp.stackexchange.com/questions/62259/… Nov 29 '19 at 10:03