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

Question

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

Answer 1

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

Answer 2

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(:))))]);

[First answer]

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)