GAUSSFFT convolves an image IMG with Gaussian kernels via FFT. Kernels with SIGMA smaller than 0.4 pixels will be neglected since the algorithm gets very slow and thus not suited for small scales. Kernels with SIGMA smaller than 10 pixels will be convolved without FFT transform. Input: img image to be convolved sigma kernel parameter of length 1 or 2 varargin first optional argument can be a look-up table for fft2 computational costs in order to predict optimal image sizes (vector format), all other arguments are kernel names, possible names are: 'G', 'Gx', 'Gy', 'Gxx', 'Gyy', 'Gxy', 'x2G', 'y2G', 'xyG' Output: varargout convolved images, corresponding to one kernel each Example: img = zeros(512); img(150, 100) = 1; fftResult = gaussFFT(img, 20, 'G'); imshow(fftResult, []); See also: buildScaleSpace, omega Licence: For internal use only. Warranty: No warranty for validity of this implementation. Author: Falko Schindler (falko.schindler@uni-bonn.de) Department of Photogrammetry Institute of Geodesy and Geoinformation University of Bonn Bonn, Germany Copyright 2009-2011
0001 function varargout = gaussFFT(img, sigma, varargin) 0002 % GAUSSFFT convolves an image IMG with Gaussian kernels via FFT. 0003 % 0004 % Kernels with SIGMA smaller than 0.4 pixels will be neglected since the 0005 % algorithm gets very slow and thus not suited for small scales. 0006 % Kernels with SIGMA smaller than 10 pixels will be convolved without FFT 0007 % transform. 0008 % 0009 % Input: 0010 % img image to be convolved 0011 % sigma kernel parameter of length 1 or 2 0012 % varargin first optional argument can be a look-up table for fft2 0013 % computational costs in order to predict optimal image 0014 % sizes (vector format), 0015 % all other arguments are kernel names, possible names are: 0016 % 'G', 0017 % 'Gx', 'Gy', 0018 % 'Gxx', 'Gyy', 'Gxy', 0019 % 'x2G', 'y2G', 'xyG' 0020 % 0021 % Output: 0022 % varargout convolved images, corresponding to one kernel each 0023 % 0024 % Example: 0025 % img = zeros(512); 0026 % img(150, 100) = 1; 0027 % fftResult = gaussFFT(img, 20, 'G'); 0028 % imshow(fftResult, []); 0029 % 0030 % See also: buildScaleSpace, omega 0031 % 0032 % Licence: 0033 % For internal use only. 0034 % 0035 % Warranty: 0036 % No warranty for validity of this implementation. 0037 % 0038 % Author: 0039 % Falko Schindler (falko.schindler@uni-bonn.de) 0040 % Department of Photogrammetry 0041 % Institute of Geodesy and Geoinformation 0042 % University of Bonn 0043 % Bonn, Germany 0044 % 0045 % Copyright 2009-2011 0046 0047 0048 %% get two values for sigma not smaller than 0.4 0049 % duplicate a single value 0050 if numel(sigma) == 1 0051 sigma = [sigma, sigma]; 0052 end 0053 % neglect small values 0054 sigma = sigma .* (sigma >= 0.4); 0055 0056 %% pad and transform image 0057 % image padding in spatial domain, depending on sigma 0058 padding = 4; 0059 % default size of padding 0060 pad = round(sigma * padding); 0061 % use LUT 0062 if isnumeric(varargin{1}) 0063 % use LUT for fft2 computational times 0064 time = varargin{1}; 0065 varargin = varargin(2:end); 0066 % minimum size of padding 0067 padMin = round(sigma * padding); 0068 % maximum size of padding 0069 padMax = (2.^ceil(log2(2 * padMin + size(img))) - size(img)) ./ 2; 0070 % optimal size of padding 0071 pad = zeros(1, 2); 0072 for dim = 1 : 2 0073 % possible image sizes 0074 sizes = size(img, dim) + 2 * (padMin(dim) : padMax(dim)); 0075 % fastest image size 0076 [t, idx] = min(time(sizes <= numel(time))); 0077 % optimal padding 0078 if ~isempty(idx), pad(dim) = (sizes(idx) - size(img, dim)) / 2; end 0079 end 0080 end 0081 % signal f: padded image 0082 f = img([ones(1, pad(1)), 1 : size(img, 1), end * ones(1, pad(1))], ... 0083 [ones(1, pad(2)), 1 : size(img, 2), end * ones(1, pad(2))]); 0084 % Fast Fourier Transformation of the image 0085 F = fft2(f); 0086 % size of padded image 0087 N = size(F); 0088 0089 %% gaussians Hc for rows and columns separately 0090 % transition function: gaussians 0091 Hc = cell(1, 2); 0092 % build them for rows and columns separately 0093 for dim = 1 : 2 0094 % circular repetitions in frequency domain 0095 if sigma(dim) > 0 0096 nShah = ceil(6 / sigma(dim)); 0097 else 0098 nShah = 0; 0099 end 0100 % multiple repetitions of each pixel 0101 [j, n] = ndgrid(-nShah : nShah, 0 : N(dim) - 1); 0102 % gaussian value at each position 0103 Hc{dim} = exp(-2 .* (sigma(dim) .* pi .* (j + n ./ N(dim))).^2); 0104 end 0105 % sum over all repetitions 0106 Hy_temp = sum(Hc{1}, 1); 0107 Hx_temp = sum(Hc{2}, 1); 0108 % normalization (divide by N, so it won't have to be done when using fft2 0109 % instead of ifft2) 0110 Hy = Hy_temp ./ Hy_temp(1, 1) ./ N(1); 0111 Hx = Hx_temp ./ Hx_temp(1, 1) ./ N(2); 0112 0113 %% construct kernels 0114 % see handwritten notes 0115 dy = - i * sin(2 * pi * (0 : N(1) - 1) ./ N(1)); 0116 dx = i * sin(2 * pi * (0 : N(2) - 1) ./ N(2)); 0117 dyy = - 2 * (1 - cos(2 * pi * (0 : N(1) - 1) ./ N(1))); 0118 dxx = - 2 * (1 - cos(2 * pi * (0 : N(2) - 1) ./ N(2))); 0119 for k = 1 : length(varargin) 0120 switch varargin{k} 0121 case 'G', H = Hy' * Hx; 0122 case 'Gy', H = (Hy .* dy)' * Hx; 0123 case 'Gx', H = Hy' * (Hx .* dx); 0124 case 'Gyy', H = (Hy .* dyy)' * Hx; 0125 case 'Gxx', H = Hy' * (Hx .* dxx); 0126 case 'Gxy', H = (Hy .* dy)' * (Hx .* dx); 0127 case 'y2G', H = (Hy .* dyy)' * Hx * sigma(1)^4 + ... 0128 Hy' * Hx * sigma(1)^2; 0129 case 'x2G', H = Hy' * (Hx .* dxx) * sigma(2)^4 + ... 0130 Hy' * Hx * sigma(2)^2; 0131 case 'xyG', H = (Hy .* dy)' * (Hx .* dx) * prod(sigma)^2; 0132 otherwise, error 'gaussFFT: unknown kernel name'; 0133 end 0134 % fft convolution (inverse fft transformation via ifft, normalization 0135 % is already done) 0136 G = real(fft2((F .* H)')'); 0137 % crop result 0138 varargout{k} = G(pad(1) + 1 : end - pad(1), pad(2) + 1 : end - pad(2)); 0139 end