# What is the most efficient way to perform Gaussian blurs for varying sizes of Gaussian kernels in an image?

I have detected corner points of images using Harris corner detector in GLSL and now, I intend on using Fast Retina Keypoint(FREAK) as keypoint descriptor to later match the keypoints between images.

FREAK uses circular sampling pattern, where each circles represent Gaussian kernels applied to the corresponding sampling points.

And, the size of Gaussian kernel increase exponentially for points farther away from central keypoint. I previously used nine hit Gaussian blur as mentioned here. So, what is the most efficient way to perform Gaussian blurs for a number of Gaussian kernels with varying sizes, which also overlap themselves, in an image?

• Is it important to have a perfect hexagonal sampling pattern? Are memory reads and writes fast or slow compared to calculations? Oct 30, 2015 at 14:31
• thank you very much for replying. I have implemented hexagonal sampling patterns as you mentioned. For varying size of kernel, I followed the given link above that uses 9 hit gaussian blur, but changed the radius form 9 to increasing powers of 2. I have yet to test the accuracy of output but it works well for matching between translated images only(not of differing orientations). Takes nearly 350ms from start of drawing procedure. Is that slow or fast? Nov 1, 2015 at 2:28
• With gaussians, one nice thing is that we can do separable convolution, First separate x and y dimensions, but also each dimension can be factored into a set of smaller filters which could speed up a lot. I think the speed would depend a lot on hardware. Does it take 350 ms for one frame? Nov 1, 2015 at 11:42
• yeah, 350 ms for a single frame. Well, I have to state that this time taken includes keypoint detection followed by description but I guess it's kind of slow. Nov 7, 2015 at 10:36

## 1 Answer

I'm not really used to thinking in texture fetches as I usually code high level non-hardware-specific signal processing, but low pass filtering also removes information you can couple it with sub-sampling which can give further gains. It is not clear to me if you have already done this.

In practice this means that after the first low-pass filtering, you can disregard $k$ out of $N$ of the samples in the subsequent low-pass filtering.

Simple example is the $2^n \times 2^n$ averaging-filter, first we can do 2x2 mean values and store. Then for the next pass, we only need to use the output of one of each $2 \times 2$ blocks. So we only need consider $1$ in $4^k$ data points at step $k$.

A 8x8 block this would mean:

$$I_1 = \left[\begin{array}{cccccccc} 1&1&1&1&1&1&1&1\\ 1&1&1&1&1&1&1&1\\ 1&1&1&1&1&1&1&1\\ 1&1&1&1&1&1&1&1\\ 1&1&1&1&1&1&1&1\\ 1&1&1&1&1&1&1&1\\ 1&1&1&1&1&1&1&1\\ 1&1&1&1&1&1&1&1 \end{array}\right]$$ First iteration, all are involved in 2x2 filtering. But they are only stored at the pixels which are ones here (so we only need to calculate 1/4 of them): $$I_2 = \left[\begin{array}{cccccccc} 1&0&1&0&1&0&1&0\\ 0&0&0&0&0&0&0&0\\ 1&0&1&0&1&0&1&0\\ 0&0&0&0&0&0&0&0\\ 1&0&1&0&1&0&1&0\\ 0&0&0&0&0&0&0&0\\ 1&0&1&0&1&0&1&0\\ 0&0&0&0&0&0&0&0 \end{array}\right]$$ Next step only involves the ones in the previous one, now our image is 4x4 instead of 8x8, so calculations are reduced a factor of 4. But we only need to store the results at pixels: $$I_3 = \left[\begin{array}{cccccccc} 1&0&0&0&1&0&0&0\\ 0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0\\ 1&0&0&0&1&0&0&0\\ 0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0 \end{array}\right]$$ We will end up with only one out-value, in this case upper left. And here's where the interpolations come into play. We can interpolate between these upper left results of neighboring "blocks". The interpolation kernel can be dependent on the low-pass filter, but I assume that this is where the hardware shines with it's strengths, having really fast ways to do this for us?

If you are interested in theory filter-banks, sub-band coding and maybe wavelets are interesting to read up on. The example I brought up here is the implementation low-pass filter of the basic but popular discrete Haar-wavelet. Gaussian pyramids in computer vision may be even more relevant. There's much to read about it so it's difficult to cram it all into just one answer.

• thanks for the answer. I used the method given in this link: rastergrid.com/blog/2010/09/… , so I can't actually reuse the values from previous samples as those samples aren't used. Nov 8, 2015 at 8:15