123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475 ``````function Y = fast_gauss( X, sigma, do_norm ) % Low-pass filter image using the Gaussian filter % % Y = blur_gaussian( X, sigma, do_norm ) % % do_norm - normalize the results or not, 'true' by default % % Use FFT or spatial domain, whichever is faster %ksize2 = min(max(3, ceil(((sigma-0.8)/0.3 + 1)*2)), max(size(X))*2); if( ~exist( 'do_norm', 'var' ) ) do_norm = true; end if( do_norm ) norm_f = 1; else norm_f = sigma*sqrt(2*pi); end if( sigma >= 4.3 ) % Experimentally found threshold when FFT is faster ks = [size(X,1) size(X,2)]*2; NF = pi; % xx = [linspace( 0, NF, ks(2)/2 ) linspace(-NF, 0, ks(2)/2) ]; % yy = [linspace( 0, NF, ks(1)/2 ) linspace(-NF, 0, ks(1)/2) ]; xx = [linspace( 0, NF, ks(2)/2 ) linspace(-NF, -NF/(ks(2)/2), ks(2)/2) ]; yy = [linspace( 0, NF, ks(1)/2 ) linspace(-NF, -NF/(ks(1)/2), ks(1)/2) ]; [XX YY] = meshgrid( xx, yy ); K = exp( -0.5*(XX.^2 + YY.^2)*sigma^2 ) * norm_f; Y = zeros( size(X) ); for cc=1:size(X,3) Y(:,:,cc) = fast_conv_fft( X(:,:,cc), K, 'replicate' ); end else ksize = round(sigma*5); h = fspecial( 'gaussian', ksize, sigma ) * norm_f; Y = imfilter( X, h, 'replicate' ); end end function Y = fast_conv_fft( X, fH, pad_value ) % Convolve with a large support kernel in the Fourier domain. % % Y = fast_conv_fft( X, fH, pad_value ) % % X - image to be convolved (in spatial domain) % fH - filter to convolve with in the Fourier domain, idealy 2x size of X % pad_value - value to use for padding when expanding X to the size of fH % % (C) Rafal Mantiuk % This is an experimental code for internal use. Do not redistribute. pad_size = (size(fH)-size(X)); %mX = mean( X(:) ); fX = fft2( padarray( X, pad_size, pad_value, 'post' ) ); Yl = real(ifft2( fX.*fH, size(fX,1), size(fX,2), 'symmetric' )); Y = Yl(1:size(X,1),1:size(X,2)); end ``````