Home > Matching_SYM_LSM > src > Functions > LSM_Functions > gaussFFT.m

gaussFFT

PURPOSE ^

GAUSSFFT convolves an image IMG with Gaussian kernels via FFT.

SYNOPSIS ^

function varargout = gaussFFT(img, sigma, varargin)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Sun 19-Jul-2020 23:00:25 by m2html © 2005