How to compute 2D convolution using 1D convolution over rows and columns?

Question

I'm looking for a method to compute the same result that conv2 will give, by using conv in MATLAB. (I'm doing the convolution in C code and I need to compare the result between MATLAB and my code.)

I have heard that it's possible to use 1D convolutional vector multiplication to achieve the same result as 2D convolutional matrix multiplication, only if the kernel K is separable e.g gaussian kernel. Also, by doing this way, 1D-way is much faster than computing the 2D way.

Here I have collected some Matlab code:

% Create data
A = randn(10, 10);

% Create gaussian kernel
sigma = 1;
kernel_size = round(6 * sigma);

% Create mesh grid
[x, y] = meshgrid(-kernel_size:kernel_size, -kernel_size:kernel_size);

% Create gaussian 2D kernel
K = 1/(2*pi*sigma^2)*exp(-(x.^2 + y.^2)/(2*sigma^2));

I wonder how I can achieve the same result as conv2(A, K, 'same') if I'm using conv only. Should I first do conv(A(i, :); K(i, :), 'same') for each row and then conv(A(:, i); K(:, i), 'same') for each column?

Answer 1

If the 2-D kernel K is separable into a row vector kr and a column vector kc such that K = kc*kr, the result of conv2(A, K, 'same') can be obtained more efficiently as conv2(kr, kc, A, 'same') (or equivalently as conv2(kc, kr, A, 'same')).

If you need to use only conv, you have to convolve each row of A with kr, and then each column of the result with kc (or equivalently in the opposite order).

So, in your example:

A = randn(10, 10);
sigma = 1;
kernel_size = round(6 * sigma);
[x, y] = meshgrid(-kernel_size:kernel_size, -kernel_size:kernel_size);
K = 1/(2*pi*sigma^2)*exp(-(x.^2 + y.^2)/(2*sigma^2));

result = conv2(A, K, 'same');

you can obtain the same result with

x = -kernel_size:kernel_size;
y = x.';
kr = 1/sqrt(2*pi)/sigma * exp(-x.^2/(2*sigma^2));
kc = 1/sqrt(2*pi)/sigma * exp(-y.^2/(2*sigma^2));

result_sep = NaN(size(A));
for r = 1:size(A,1)
    result_sep(r,:) = conv(A(r,:), kr, 'same');
end
for c = 1:size(A,1)
    result_sep(:,c) = conv(result_sep(:,c), kc, 'same');
end

Check:

max(max(abs(result-result_sep)))

gives a value of the order of eps, so the two results are the same up to floating-point inaccuracies.