% INTEGGAUSSFILT - Approximate Gaussian filtering using integral filters
%
% This function approximates Gaussian filtering by repeatedly applying
% averaging filters. The averaging is performed via integral images which
% results in a fixed and very low computational cost that is independent of
% the Gaussian size.
%
% Usage: [fim, sigmaActual] = integgaussfilt(im, sigma, nFilt)
%
% Arguments:
% im - Image to be Gaussian smoothed
% sigma - Desired standard deviation of Gaussian filter
% nFilt - The number of average filterings to be used to
% approximate the Gaussian. This should be a minimum of
% 3, using 4 is better. If the smoothed image is to be
% differentiated an additional averaging should be applied
% for each derivative. Eg if a second derivative is to be
% taken at least 5 averagings should be applied. If omitted
% this parameter defaults to 5.
%
% Returns:
% fim - Smoothed image
% sigmaActual - Actual standard deviation of approximate Gaussian filter
% that was used
%
% Notes:
% 1. The desired standard deviation will not be achieved exactly. A combination
% of different sized averaging filters are applied to approximate it as closely
% as possible. If nFilt is 5 the deviation from the desired standard deviation
% will be at most about 0.15 pixels.
%
% 2. Values of sigma less than about 1.8 cannot be well approximated by
% repeated averagings. For sigma < 1.8 the smoothing is performed using
% conventional Gaussian convolution.
%
% See also: INTEGAVERAGE, SOLVEINTEG, INTEGRALIMAGE, GAUSSFILT
% Copyright (c) 2009-2010 Peter Kovesi
% www.peterkovesi.com/matlabfns/
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.
% September 2009 - Original version
% April 2010 - Added return of actual standard deviation of effective
% filter used. Use standard convolution for small sigma
function [im, sigmaActual] = integgaussfilt(im, sigma, nFilt)
if ~exist('nFilt', 'var')
nFilt = 5;
end
% First check if sigma is too small to be well represented by repeated
% averagings. 5 averagings with a width 3 filter produces an equivalent
% sigma of ~1.8 This represents the minimum threshold. For sigma less
% than this we use conventional convolution
if sigma < 1.8
im = gaussfilt(im, sigma);
sigmaActual = sigma;
else % Use repeated averagings via integral images
% Solve for the combination of averaging filter sizes that will result
% in the closest approximation of sigma given nFilt.
[wl, wu, m, sigmaActual] = solveinteg(sigma, nFilt);
radl = (wl-1)/2;
radu = (wu-1)/2;
% Apply the averaging filters via integral images.
for i = 1:m
im = integaverage(im,radl);
% im = runningaverage(im,radl);
end
for n = 1:(nFilt-m)
im = integaverage(im,radu);
% im = runningaverage(im,radu);
end
end
%-------------------------------------------------------------------
% GAUSSFILT - Small wrapper function for convenient Gaussian filtering
%
% Usage: smim = gaussfilt(im, sigma)
%
function smim = gaussfilt(im, sigma)
sze = ceil(6*sigma);
if ~mod(sze,2) % Ensure filter size is odd
sze = sze+1;
end
sze = max(sze,1); % and make sure it is at least 1
h = fspecial('gaussian', [sze sze], sigma);
smim = filter2(h, im);