Help with convolution in matlab

There are many convolutions in my code and I am looking for a way to reduce the computation of the convolution. Since I am using conv2(a,b,'same'), I was wondering if it is possible to just calculate the center part that I want? Also, 'a' is a 2D matrix while 'b' is a 1D vector so will it be faster if I repmat 'b' to the same length as 'a'?

Looking at the convolution theorem, I do not understand how the 2D convolution works. I guess I understand the idea of sliding and summing them but from the convolution theorem, shouldn't the size of both matrix be the same?

Answer 1

[Image Analyst]

Convolution is sliding a matrix ("template" or "kernel") of any size over another matrix of any size. Usually the "sliding" matrix is the smaller one and it can be a 1-D vector or a 2-D matrix. If you want a 1-D vector then just use it because if you replicate it in the other dimension the result will be different because you'd be including different elements in the computation. For example a 1x3 array give the average of 3 elements only in that row while a 3x3 kernel/template gives the average of 9 elements - both in the same row and the row above and below the target row, so they're not the same.

If you're using images, you'll need to cast from integer to double. And to keep the image in the same range you'll want the sum of the sliding elements (template) to be 1.

grayImage = imread('cameraman.tif');
template = ones(9,9); % Box blur / average in a 9x9 window.
template = template / sum(template(:))

template = 9×9

0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123 0.0123

blurredImage = conv2(double(grayImage), template, 'same'); % Blur the image.
subplot(1, 2, 1);
imshow(grayImage)
subplot(1, 2, 2);
imshow(blurredImage, [])

If you want to filter only a part of the image/matrix you can crop it out, then convolve it, then assign it back in. See attached copy and paste demo.

The conv2 function, or imfilter function, are highly optimized, expecially for common situations (like uniform or symmetric templates). I doubt you would be able to improve on their speed, unless, like you said, you wanted to just filter a very small part of a larger matrix.

If you don't want edge effects, where the sliding template goes off the edge of the image you can use the 'valid' option to use only those parts where the sliding template is completely within the main image. If you use imfilter there are a lot more options for how to handle edge effects than conv2 has.

[Jiayun Liu]

Thanks for the detailed explanation. Yeah I guessed these base functions are highly optimized and it will be difficult to come up with a better one. But do you know how 2D convolution with 1D vector is calculated in matlab?

Currently I am using symmetric extension to solve the edge effects and it seems to be working for my application.

[Image Analyst]

For example, to blur in just one direction (horizontally), try this.

grayImage = imread('cameraman.tif');
template = ones(1,29); % Box blur / average in a 1x29 window.
template = template / sum(template(:))

template = 1×29

0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345

blurredImage = conv2(double(grayImage), template, 'same'); % Blur the image.
subplot(1, 2, 1);
imshow(grayImage)
subplot(1, 2, 2);
imshow(blurredImage, [])

If you want to know how convolution works, you basically get a window of the image and multiply it and add it to your sliding template. When the template moves over by one pixel, you can just add in the new pixels and subtract off the pixels no longer in the window. That way you don't have to process all of the pixels in the window, saving a lot of time.

[Jiayun Liu]

Ok so matlab doesn't using convolution theorem (FT signal and kernel multiply them and IFT) to do the calculation?

[Image Analyst]

I'm not sure what they do internally. Theoretically that is mathematically the same, and may be a faster way for larger arrays, but I'm not sure how the accuracy compares, considering truncation and roundoff errors.

[Jiayun Liu]

I see. Thanks for your help

[Walter Roberson]

If I recall correctly, historically conv2(A,B) worked by convolving the columns of A with the columns of B, and then transposing the results and transposing B, and doing the same process.

I do not recall whether historically the 2D convolution with a vector was optimized. It might potentially have been, as that could be implemented in a more vectorized way... I just do not recall now.

Also, under at least some circumstances that were not clear to me, it would find seperable kernels and had more efficient ways of dealing with those. Details were buried in a mex file. There were hints that fft's might be involved.

[Image Analyst]

I also don't remember the mathematical details but if I remember correctly (and may not since last time I delved into it was literally 20 years ago) there were certain types of kernels that lent themselves to "the method of separable kernels", like the kernel matrix was symmetric or something like that, and if the kernel is right you can use the method which is like two passes over the whole image with a 1-D kernel : one in the horizontal direction and the second pass in the vertical direction. If the method could not be used you had to do the way where you added in the new right edge and/or bottom edge and subtract off the old left edge and/or top edge. Or you could do it via Fourier Transform. There is some size where if it's bigger, you're faster doing the Fourier way and if it's smaller it's faster to do it in the spatial domain. If they do switch domains for speed I'm sure they will check the matrix sizes to ensure that they use the fastest way. If the arrays are integers, they may do it only in the spatial domain instead of the Fourier domain, for accuracy, but I'm not sure what they do or what decisions are made internally in the mex file.

[Jiayun Liu]

How did you delve into it before? As in how can I have access to the internal calculation?

So matlab might use the FT method to calculate convolution? If that is the case, can it be done with kernal size that is different from the signal size? I mean it should be IFT[G(kx)*H(kx)] where G(kx)&H(kx) are the FT of g(x) & h(x). How to do the multiplication if the sizes are different?

[Image Analyst]

You can just Google convolution methods. If the matrixes are different sizes you have to zero pad one.

[Walter Roberson]

In recent versions, conv2 is built-in, completely compiled.

But there was a time when conv2 was a .m file that could be read.

My memory is that it used to have a couple of code branches, one of which was to call conv(), transpose, call conv() again.

... but now when I look, I see that conv() these days is implemented in terms of calling built-in conv2.

[Jiayun Liu]

I see. In that case I guess I will just have to try writing and compare the speed.

Thanks

[Walter Roberson]

I just checked back to R2013a, and conv2() was already built-in then...