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 <mantiuk@gmail.com>
- % 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
|