# Separable convolution 26

Posted by **Steve Eddins**,

The sci.image.processing newsgroup had a discussion this week on separable filters, which reminded me that separability has been on my blog topic ideas list for a while now.

### Contents

### What is a separable filter?

A two-dimensional filter kernel is *separable* if it can be expressed as the outer product of two vectors. For example, let's look at a Sobel kernel.

s = [-1 0 1; -2 0 2; -1 0 1]

s = -1 0 1 -2 0 2 -1 0 1

This kernel can be written as a matrix product of a column and a row vector.

v = [1; 2; 1]

v = 1 2 1

h = [-1 0 1]

h = -1 0 1

v * h

ans = -1 0 1 -2 0 2 -1 0 1

### Associativity of convolution

As it turns out, the matrix product of a column vector and a row vector is equivalent to the two-dimensional convolution of the two vectors.

conv2(v, h)

ans = -1 0 1 -2 0 2 -1 0 1

This matters because convolution is associative. That is,

(I've used the asterix here to mean convolution.) So if a filter *s* is separable:

then you can filter with *s* by filtering first with *v*, and then filtering the result with *h*.

### Computational advantage of separable convolution

Why would you want to filter this way? Because you can do it faster. Filtering an M-by-N image with a P-by-Q filter kernel
requires roughly *MNPQ* multiplies and adds (assuming we aren't using an implementation based on the FFT). If the kernel is separable, you can filter
in two steps. The first step requires about *MNP* multiplies and adds. The second requires about *MNQ* multiplies and adds, for a total of *MN(P + Q)*.

The computational advantage of separable convolution versus nonseparable convolution is therefore:

For a 9-by-9 filter kernel, that's a theoretical speed-up of 4.5.

That's enough for now. Next time, I'll write about how to determine whether a filter kernel is separable, and what MATLAB and toolbox functions test automatically for separability.

Get
the MATLAB code

Published with MATLAB® 7.3

## 26 CommentsOldest to Newest

Hi Steve,

I would like to know more about 2D seperable filters, have u got any reference sites to recommend?

In addition, if i am using remez () in matlab to design 1D FIR filter and then usetrans() to transform it to a 2D Filter, is this a seperable or non-seperable filter? In general, how do I create a 2D separable filter?

Thanks alot

Foo

Foo – I have reference books, not reference sites. :-) For example,

Multidimensional Digital Signal Processing, by Dudgeon and Mersereau. The bookDigital Image Processing Using MATLABtalks about separability, but not in any more detail than I’ve given here. You should also look at my separability – part 2 post. It discusses how to detect separability.I’m not sure which function you are referring to for creating a 2-D filter. It looks like there’s a typo in your comment. Can you clarify?

A filter is separable if it has rank 1, so you can use the MATLAB function

rank>to see whether a given filter is separable.You can create a separable filter by forming the outer product of two 1-D filters.

Hi Steve:

I wanna construct a blur matrix by using ‘conv2′ command,however, my matrix is sparse matrix that ‘conv2′ doesn’t support, So, is there other similar functions I can use?

Woo – I’m not aware of MATLAB functions for performing convolution on sparse matrices. You could possibly implement something yourself, based the notion of convolution as a sum of weighted and shifted versions of the convolution filter.

Dear sir how to perform Polygonal approximation of an image using matlab?

Shan—Can you be more specific? What are your inputs; what specific problem do you want to solve; how do you want your outputs expressed?

Steve,

You talked about advantage of being separable without using fft. How about the comparison of separable and nonseparble 2D convolution with the aid of fft? Thanks.

Dan—Filter separability makes no practical difference when using the FFT to perform convolution. The 2-D FFT is separable anyway, regardless of the properties of the filter. In fact, MATLAB performs all multidimensional FFTs by exploiting the separability of the transform.

Hi Steve

When you convolve two matrices f and g using conv2 the dimension of f*g is more than f and g. so if i wanna deconvolve what do i need to do? g=iff2(fft2(f*g)./fft2(f)) does not agree in dimension. If I do zero padding then it does not give the right result. Thanks

Kaushal—Even if you get the padding right, you’ll have problems if fft2(f) has zero values or very small values. You might want to take a look at Stan Reeves’ deblurring series in this blog.

sir

how can I convolve an image with blurring filter [0.25 .5 0.25]

Reji—Use

imfilterin the Image Processing Toolbox, or useconv2orfilter2in MATLAB.sir

i used imfilter to convolve the image in each dimension.i’m trying to implement an algorithm for my research work.in the algorithm I want to blur and scale down an image.in the algorithm they used the filter[0.25 0.5 0.25].sir what type filter is this?

I = imread(‘peppers.png’)

figure;imshow(I);

h=[0.25 0.5 0.25];

blur1_m=imfilter(I, h,’conv’);

figure;imshow(blur1_m);

I used this program to convolve the image with [0.25 0.5 0.25] in each direction.the result i’m getting had no much difference with the original one

reji

Reji—If you filter with [0.25 0.5 0.25], each pixel is replaced 0.5 times itself plus 0.25 times its left neighbor plus 0.25 times its right neighbor. The filter will have a small smoothing effect.

Sir, I am trying to convolve an image with a gaussian derivative kernel for edge detection. Is there any fast way to do this? Thanks

Ray—If you search for “fast Gaussian filtering algorithms” you should see several papers describing techniques.

Can you help me to perform 2d convolution of a large image (10000×10000) with conv2 matlab function.

Hi Steve,

I am just entering the image processing realm and have what is probably a very basic question.

Let G_MxM be a separable filter given by

Gx_MxM = d_Mx1 * s_1xM,

where “*” denotes matrix multiplication.

The result of filtering [1] an input matrix X at row n and column m by

this Gx matrix is then

Yx_M(n,m) = X_M(n,m) * Gx

= X_M(n,m) * d * s

= (X_M(n,m) * d) * s,

where X_M(n,m) and Yx_M(n,x) are the MxM matrices at n, m.

Further expanding this gives

Yx(n,m) = s(m) * \sum_{k = 0}{M – 1} X(n, m + k) d(k),

where X(n,m) and Yx(n,m) are the elements of the matrices X and Yx at row

n and column m, respectively.

Now note that the product of X_M and d is a 1D convolution of the filter

d along the nth row of X_M, and the value s(m) simply scales that

convolution.

Thus if d is the simple 1D difference filter [1 0 -1] (as it is for a

Sobel filter), then this operation produces the derivative of X_M along

the nth row of X_M, or what I would call the gradient in the x

direction.

Fine. Great.

But then, if Gy = Gx’, as it is here for a Sobel filter,

http://en.wikipedia.org/wiki/Sobel_operator

then the gradient in the y direction is

Yy_M(n,m) = X_M(n,m) * Gy

and therefore

Yy(n,m) = d(m) * \sum_{k = 0}{M – 1} X(n, m+k) s(k).

By the same mathematics as above, we see that the product

of X_M and s is a 1D convolution of the filter s along

the nth row of X_M.

But this is NOT the gradient of X in the y direction, since such a

gradient would have to be performed on the COLUMNS of X.

So from this reasoning, it seems that to determine the x-gradient

of X we should compute

X * Gx

but to determine the y-gradient of X we should compute

Gy * X.

Is this correct? If so, then isn’t is also true that to perform these

gradients correctly one must have a side-specific image filter

operation? Specifically, my question is, if this is true, how can the

Matlab/Octave image processing library’s imfilter(X,Gx) and

imfilter(X,Gy) compute both gradients correctly, since it presumably

multiplies on only one “side” and thus one would be incorrect?

–Randy

[1] We ignore various details such as edge effects and the sign-reversal

that standard convolution requires here to simplify the notation and

presentation.

Steve, ignore that previous post – i was getting confused between the multiplication operator and the convolution operator.

Hey Steve, thanks for the insight into separable kernels. I would like to know how what you have written in the initial posts applies to 3D.

For instance I have a 3D mask and intend to convolve with a guassian. I know that the gaussian is separable, how do I go about separating the kernel for convolution in 3D ?

regards

GT

GT—Just filter three times in succession. Use a one-dimensional Gaussian filter oriented along the 1st dimension, then along the 2nd dimension, then along the 3rd.

Hi Steve,

Concerning the built-in n-D convolution methods:

MATLAB (v2008b) keeps failing due to out-of-memory errors when convolving large 3D arrays (e.g. (256)^3) with a ~16^3 (general, non-separable) kernel.

I’ve tried imfilter, convn or convnfft from the filexchange, same result …

Would it make sense to add an functionality that transparently breaks down the task into filtering smaller sub-arrays and stitches together the results, depending on the free contiguous memory blocks, to prevent memory errors? Or can you think of another solution ?

Thanks

Vasek

Vasek—I believe convn does not use extra memory besides that required for the input and output arrays. So if you are running out of memory with convn, a block-based approach will not help you. You might want to look at the “Avoiding Out of Memory Errors” support page for more information and tips.

Hi Steve,

You describe how to determine if a filter is separable in the time domain. What if I have an filter in the frequency domain only, and want to recognize if it is separable, and possibly write the h1 and h2 components of the filter in the time domain, whithout using ifft2?

Hi Steve, I was wondering whether or not this content of this article relates at all to a Separability filter – read it in a paper and I have not been able to get a clear definition.

Thanks!

Jonathan—If you mean “separable filter,” then yes, that is the same thing as separable convolution.