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