Separable convolutions
The convolution is an important tool in image processing. All linear filters are convolutions (like, for example, the Gaussian filter). One of the reasons linear filters are so prevalent is that they imitate physical systems. For example, an analog electric circuit (containing resistors, capacitors and inductors) can be described by a convolution. The projection of a cell on a slide, through the microscope’s optical systems, onto the CCD sensor, can be described by a convolution. Even the sampling and averaging that occur in the CCD can be described by a convolution.
Two properties of the convolution are quite interesting when looking for an efficient implementation:

The convolution is a multiplication in the Fourier domain,
This means that you can compute the convolution by applying the Fourier transform to both the image and the convolution kernel, multiplying the two results, then inverse transforming the result.

The convolution is associative:
This post is about the repercussions of this property.
Given this associativity property, if we have a filtering kernel that can be decomposed into two simpler kernels such that , then we can compute the convolution as a sequence of two cheaper convolutions:
The most trivial example of this is a square, uniform filter:
a = double(deltaim(40,40));
dipshow(a,'lin')
h1 = ones(15,1)/15;
h2 = ones(1,15)/15;
b = conv2(a,h1,'same');
dipshow(b,'lin')
c = conv2(b,h2,'same');
dipshow(c,'lin')
Convolving with two onedimensional uniform filters in sequence is equivalent to convolving with a square, 2D uniform filter. But instead of doing 15x15=225 multiplications and additions per pixel, we’re only doing 15+15=30. Whenever we can decompose an ndimensional filter into a set of n onedimensional filters, we say that the filter is separable. Some of the most useful linear filters are separable.
Decomposing a separable filter kernel
Given that is a column vector (a filter for the yaxis), and a row vector, their product is a matrix and
identical to the convolution of the one with the other.
That is, .
We can thus see how each column of a separable kernel is a scaled version of any other column (and conversely,
each row is a scaled version of any other row). In MATLAB terms, h(:,i)=h1*h2(i)
.
To decompose the kernel into two vectors, one thus can take any column of h
as h1
, and determine the values
of h2
by simple division:
h = [2,1,0,1,2]'*[3,2,1];
h1 = h(:,1);
h2 = h(1,:)/h1(1);
isequal(h1*h2,h)
ans =
1
Note, however, that the h1
and h2
computed in the code above are not identical to the two vectors that we multiplied
together to obtain h
. We can multiply h1
by any constant, as long as we divide h2
by that same constant.
The problem with this method is that if he happened to pick a column for h1
that is all zeros, h2
will be all
infinity, and we wouldn’t get anywhere. Also, if some columns of h1
have very small values, the decomposition result
would be less accurate if we used that column compared to any other column. By adding a little bit of logic to the code
above, we can easily get the most accurate decomposition:
[~,I] = max(sum(abs(h))); % Pick the column with largest values
h1 = h(:,I);
[~,I] = max(sum(abs(h),2)); % Pick the row with largest values
h2 = h(I,:)/h1(I);
isequal(h1*h2,h)
ans =
1
Let’s look at a slightly more complex example for a separable filter: the Gaussian derivative. We generate the kernel simply by applying the filter to a delta image, an image with a single pixel set:
h = deltaim(31,31);
h = dx(h,5);
h = double(h);
[~,I] = max(sum(abs(h)));
h1 = h(:,I);
[~,I] = max(sum(abs(h),2));
h2 = h(I,:)/h1(I);
isequal(h1*h2,h)
ans =
0
The filter, though it is separable, failed the isequal
test. We shouldn’t expect an exact equality, though, because we
are dealing with floatingpoint values. The function dx
returns singleprecision floating point values, this is the
maximum accuracy we should expect in the end result:
q = h1*h2h;
q = abs(q(:))/max(abs(h(:)));
max(q) < eps('single')
ans =
1
A test for separability
Consider the matrix h
as a system of linear equations. Because each row is a scaled version of every other row, this
is actually a set of identical equations. That is, the rank of the system is 1:
h = [2,1,0,1,2]'*[3,2,1];
rank(h)
ans =
1
But, again, in general we cannot expect infinite accuracy, we should check the rank taking a certain tolerance into account:
h = double(dx(deltaim(31,31),5));
s = svd(h);
sum(s > eps('single'))
ans =
1
If we decide to check the separability of the kernel before computing its components, we might as well use svd
to its
full power. svd
computes not only the singular values, but also two matrices such that h=u*s*v'
. These two
matrices, u
and v
, contain the h1
and h2
that we are looking for:
[u,s,v] = svd(h);
s = diag(s);
if sum(s > eps('single'))==1
h1 = u(:,1)*s(1);
h2 = v(:,1)';
else
error('Kernel not separable!')
end
q = h1*h2h;
q = abs(q(:))/max(abs(h(:)));
max(q) < eps('single')
ans =
1
This last bit of code, though more elegant, will produce slightly less accurate results than the first method in those cases where the first method can compute an exact result. But if there are roundoff errors, the SVD method will, in general, produce a more accurate result.
To round off this post, I wanted to point you to two entries on the MATLAB File Exchange. One post, that I submitted, computes the kernel decomposition using the SVD method explained here. Interesting in that file is that it extends the method to ndimensional kernels by iterating over the last dimension. The other post, by DirkJan Kroon, introduces a whole different approach to kernel decomposition. I think his method is more complex and not any better than the two methods discussed here, but it is interesting because it can produce, using a least squares approach, an approximation for nonseparable kernels.