{"id":2493,"date":"2017-01-23T08:14:29","date_gmt":"2017-01-23T13:14:29","guid":{"rendered":"https:\/\/blogs.mathworks.com\/steve\/?p=2493"},"modified":"2019-11-01T17:04:39","modified_gmt":"2019-11-01T21:04:39","slug":"gaussian-filtering-with-imgaussfilt","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/steve\/2017\/01\/23\/gaussian-filtering-with-imgaussfilt\/","title":{"rendered":"Gaussian filtering with imgaussfilt"},"content":{"rendered":"

With the R2015a release a couple of years ago, the Image Processing Toolbox added the function imgaussfilt<\/tt><\/a>. This function performs 2-D Gaussian filtering on images. It features a heuristic that automatically switches between a spatial-domain implementation and a frequency-domain implementation.<\/p>

I wanted to check out the heuristic and see how well it works on my own computer (a 2015 MacBook Pro).<\/p>

You can use imgaussfilt<\/tt> this way:<\/p>

rgb = imread('peppers.png'<\/span>);\r\nsigma = 10;\r\nrgb_smoothed = imgaussfilt(rgb,sigma);\r\n\r\nimshow(rgb)\r\n<\/pre>\"\" 
imshow(rgb_smoothed)\r\ntitle('Gaussian smoothed, \\sigma = 10'<\/span>)\r\n<\/pre>\"\" 
With this call, imgaussfilt<\/tt> automatically chooses between a spatial-domain or a frequency-domain implementation. But you can also tell imgaussfilt<\/tt> which implementation to use.<\/p>
rgb_smoothed_s = imgaussfilt(rgb,sigma,'FilterDomain'<\/span>,'spatial'<\/span>);\r\nrgb_smoothed_f = imgaussfilt(rgb,sigma,'FilterDomain'<\/span>,'frequency'<\/span>);\r\n<\/pre>
I'll use this syntax, together with the timeit<\/tt> function, to get rough timing curves for the two implementation methods. The heuristic used by imgaussfilt<\/tt> uses a few different factors to decide, including image size, Gaussian kernel size, single or double precision, and the availability of processor-specific optimizations. For my computer, with a 2000-by-2000 image array, the cross-over point is at about $\\sigma = 50$. That is, imgaussfilt<\/tt> uses the spatial-domain method for $\\sigma < 50$, and it uses a frequency-domain method for $\\sigma > 50$.<\/p>
Let's check that out with a set of $\\sigma$ values ranging from 5 to 200.<\/p>
I = rand(2000,2000);\r\nww = 5:5:200;\r\nfor<\/span> k = 1:length(ww)\r\n    t_s(k) = timeit(@() imgaussfilt(I,ww(k),'FilterDomain'<\/span>,'spatial'<\/span>));\r\nend<\/span>\r\n\r\nfor<\/span> k = 1:length(ww)\r\n    t_f(k) = timeit(@() imgaussfilt(I,ww(k),'FilterDomain'<\/span>,'frequency'<\/span>));\r\nend<\/span>\r\n\r\nplot(ww,t_s)\r\nhold on<\/span>\r\nplot(ww,t_f)\r\nhold off<\/span>\r\nlegend('spatial'<\/span>,'frequency'<\/span>)\r\nxlabel('\\sigma'<\/span>)\r\nylabel('time (s)'<\/span>)\r\ntitle('imgaussfilt - 2000x2000 image'<\/span>)\r\n<\/pre>\"\" 
As you can see, the cross-over point of the timing curves is about the same as the threshold used by imgaussfilt<\/tt>. As computing hardware and optimized computation libraries evolve over time, though, the ideal cross-over point might change. The Image Processing Toolbox team will need to check this every few years and adjust the heuristic as needed.<\/p>