Prima Homepage
ImaLab
Generating C++ Modules
User Manual
The Command Shell
Running Imalab
Plugin Process
Pixels and Images
Interactive selection
Graphics: plots, profiles
Image file I/O
Image display
Connectivity Analysis
Image Processing(1)
Gaussian operators
Technical Documentation
Creating New Modules
Tutorial
Tutorial Download

[PREV][SUIV]

Gaussian operators

This page describes image operators such as direct convolution of filter masks and images, recursive implementations of Gaussian derivatives, steering of Gaussian derivatives and scale normalisation.

We have two modules for Gaussian derivatives, based on two articles by VanVliet. The old implementation (module modFilter) uses 3 coefficients for the computation and the Abramowitz approximation. Direct formulas are given for derivatives up to order 3. There are significant boundary effects, especially in the upper left corner of the image. Important: you can not trust the values near the border of the image! The boundary effects are strong within a band of width 1.5*sigma pixels around the border of the image. The functions reduce the ROI of the image accordingly. Other major problems: the Laplacian is not sufficiently isotropic and the sigma of the recursive implementation is 1.177 times larger than the sigma of a true Gaussian.

The new implementation (module modTGFilter) uses 3,4 or 5 coefficients. This means that the derivatives are more precise but can take a little longer in computation. The better precision is noted in the higher derivatives. For example the Laplacian is much more round (isotropic) than with the old implementation and the sigma corresponds to the sigma of a true Gaussian. Better approximations are achieved for large sigmas. We have direct formulas for derivatives up to order 2. The 3rd order derivatives are approximated by an additional filtering of the 2nd order. This is ok, since the Gaussian derivatives are separable. It must be noted, that higher orders amplify the error between a direct convolution and the recursive approximation.

For both modules we have implemented automatic scale selection. The intrinsic scale of a feature is noted as a maximum in the Laplacian or Gradient profile. Such a profile is computed by measuring the energy of a Laplacian (or Gradient magnitude) placed at a image position (x,y) and varying over scales. The functions are loaded with the modules.

Convolution of filters and images

Author: V. Colin de Verdiere

In the module modFilter we have implemented functions to compute convolutions of filter masks with images. This kind of operation can be very costly, depending on the size of the mask and size of the image. For computing derivatives of images, you should use the recursive Gaussian derivative implementation described below.

Following functions are implemented: 


double ImgFilterPt(TBitmap& src, TBitmapFloat& filter, int x, int y);
void ImgFilter(TBitmap& src, TBitmapFloat& dst, TBitmapFloat& filter);
TBitmapFloat* ImgFilter (TBitmap& src, TBitmapFloat& filter)
These functions compute the response of a image src to a filter given by filter, either at one point (x,y) or an entire image. 

double ConvolutionPt(TBitmap& src, TBitmap& filter, int x, int y);
double ConvolutionPt(TBitmap& src, TBitmap& filter);
void Convolution(TBitmap& src, TBitmap& dst, TBitmap& filter);
TBitmapFloat* Convolution (TBitmap& src, TBitmap& filter)
This computes the convolution of src with filter. Same as ImgFilter, but filter is flipped before computation. 

void ConvolutionSepX(TBitmap& src, TBitmap& dst, TDoubleVector& filter);
void ConvolutionSepY(TBitmap& src, TBitmap& dst, TDoubleVector& filter);
TBitmapFloat* ConvolutionSepX (TBitmap& src, TDoubleVector& filter)
TBitmapFloat* ConvolutionSepY (TBitmap& src, TDoubleVector& filter)
These functions compute the convolution only along one axis (X or Y). These functions should always called in pairs. 

void ImgGaussX121(TBitmap& src, TBitmap& dst);
void ImgGaussY121(TBitmap& src, TBitmap& dst);
void ImgGradient (TBitmap& src, TBitmap& dst);
TBitmap* ImgGradient (TBitmap& src)
These functions apply a binomial 121 filter to the image. This is the coarsest approximation of a Gaussian you can find. The gradient is computed by a filter -101.

Additional functions: 


// Calculates the entropy for the window
double Entropy(TBitmapByte& window, int size);
// Calculates the energy for the window
double Energy(TBitmapByte& window, int size);
// Extracts the window whose center pixel is xo, yo
int ExtractCentered(TBitmapByte& src, int xo, int yo, TBitmapByte& r);

back to top

Recursive Gaussian derivatives

([VanVliet95]) Authors: V. Colin de Verdiere, O. Chomat

Currently there are functions implemented for derivatives up to order 3, for the gradient and for the Laplacian. These functions are loaded by require modFilter.

Declarations: 


TBitmapFloat* FGR_XY(TBitmap& fimg, double sigma, double = -1.0);
TBitmapFloat* FGR_XY2(TBitmapFloat& fimg, double sigma);
void FGR_XY2(TBitmapFloat& fimg, double sigma, TBitmapFloat *res);
TBitmapFloat* FGR_XY3(TBitmapFloat& fimg, double sigma);
compute the Gaussian of an TBitmapFloat image fimg with size sigma. The last arggument is not used. FGR_XY is the first version of the function. FGR_XY2 and FGR_XY3 are more recent implementations.

The procedure names give hints to the content. "d" means derivative. FGR_dX computes the derivative in direction X. There are 3 prototype functions. I would use when possible the FGR_dX2, because this is the currently optimised function. 


TBitmapFloat* FGR_dX(TBitmap&, double, double = -1.0);
void FGR_dX2(TBitmapFloat&, double,TBitmapFloat*);
TBitmapFloat* FGR_dX2(TBitmapFloat&, double);

The derivatives can be called by function 


TBitmapFloat* FGR(filterType f , TBitmapFloat& fimg, double sigma);
with filterType f noting the order of the derivative. 

typedef enum {kGaussian, kDx, kDy, kDxx, kDxy, kDyy, kDxxx, kDxxy, kDxyy, kDyyy} filterType;

Additional functions: 


void FGR_AngleFromDerivatives(TBitmapFloat&, TBitmapFloat&, TBitmapFloat&);
TBitmapFloat* FGR_AngleFromDerivatives(TBitmapFloat&, TBitmapFloat&);

void FGR_Angle(TBitmap&, double, TBitmapFloat&);
TBitmapFloat* FGR_Angle(TBitmap&, double);
These functions compute the dominant angle of all (pixelwise) local features in the image. The theory comes from [Freeman91]. 

void FGR_GradFromDerivatives(TBitmapFloat&, TBitmapFloat&, TBitmapFloat&);
TBitmapFloat* FGR_GradFromDerivatives(TBitmapFloat&, TBitmapFloat&);


void FGR_Grad(TBitmap&, double, TBitmapFloat&);
TBitmapFloat* FGR_Grad(TBitmap&, double);
This function computes the gradient magnitude of all local features. 


void FGR_LapFromDerivatives(TBitmapFloat&, TBitmapFloat&, TBitmapFloat&);
TBitmapFloat* FGR_LapFromDerivatives(TBitmapFloat&, TBitmapFloat&);

void FGR_Lap(TBitmap&, double, TBitmapFloat&);
TBitmapFloat* FGR_Lap(TBitmap&, double);
This function computes the Laplacian of all local features.

back to top

Recursive Gaussian derivatives

([VanVliet98]) Authors: S. Richetto, D. Hall

Currently there are functions implemented for derivatives up to order 3, for the gradient and for the Laplacian. These functions are loaded by require modTGFilter.

GaussianFilter3, GaussianFilter4, GaussianFilter5 are classes that inherit the functions from the parent class GaussianFilter.

The X in GaussianFilterX corresponds to the number of coefficients used for the computation. The GaussianFilter3 is slightly faster, but less precise than GaussianFilter5. 


GaussianFilter5( double Sigma);
The constructor sets the sigma of the current object. This sigma can be modified by 

void SetSigma(double Sigma)

The derivative functions are following. Important: these functions compute the derivatives directly without applying a smoothing before. This has the advantage that they can be concatenated to obtain higher order derivatives. 


  void FGR_Y  (TBitmapFloat&, TBitmapFloat&, int =1 ) ;
  void FGR_X  (TBitmapFloat&, TBitmapFloat&, int =1 ) ;
  void FGR_dX (TBitmapFloat&, TBitmapFloat&, int =1 ) ;
  void FGR_dY (TBitmapFloat&, TBitmapFloat&, int =1 ) ;
  void FGR_XdY (TBitmapFloat&, TBitmapFloat&, int =1 ) ;
  void FGR_YdX (TBitmapFloat&, TBitmapFloat&, int =1 ) ;
  void FGR_dYY(TBitmapFloat&, TBitmapFloat&, int =1 ) ;
  void FGR_dXX(TBitmapFloat&, TBitmapFloat&, int =1 ) ;

When you want to compute the derivatives of images up to order three use following functions. 


  void FGR_Lx(TBitmapFloat&, TBitmapFloat&, int=1);
  void FGR_Ly(TBitmapFloat&, TBitmapFloat&, int=1);
  void FGR_Lxx(TBitmapFloat&, TBitmapFloat&, int=1);
  void FGR_Lxy(TBitmapFloat&, TBitmapFloat&, int=1);
  void FGR_Lyy(TBitmapFloat&, TBitmapFloat&, int=1);
  void FGR_Lxxx(TBitmapFloat&, TBitmapFloat&, int=1);
  void FGR_Lxxy(TBitmapFloat&, TBitmapFloat&, int=1);
  void FGR_Lxyy(TBitmapFloat&, TBitmapFloat&, int=1);
  void FGR_Lyyy(TBitmapFloat&, TBitmapFloat&, int=1);

We have implemented an analogon to FGR() of recursive filter. This function calls the above functions according to filterType defined in recursiveFilter. 


void FGR_L(filterType, TBitmapFloat&, TBitmapFloat&, int=1);

Additional functions: 


  void FGR_Lap(TBitmapFloat&, TBitmapFloat&, int =1 ) ;
  void FGR_Grad(TBitmapFloat&, TBitmapFloat&, int =1 ) ;
  void FGR_XY  (TBitmapFloat&, TBitmapFloat&, int =1 ) ;
Computes Laplacian, Gradient and Gaussian.

back to top

Steering of the derivatives

Authors: V. Colin de Verdiere, D. Hall

Gaussian derivatives are steerable. This means that a derivative can be synthetized at arbitrary angle by a linear combination of the derivatives of the same order [Freeman91]. 


void FGR_d1d(TBitmapFloat& dx, TBitmapFloat& dy, TBitmapFloat& ang_img, float alpha, TBitmapFloat& dst);
TBitmapFloat* FGR_d1d(TBitmapFloat& dx , TBitmapFloat& dy, TBitmapFloat& ang_img, float = 0.0);
This function orients the first derivative according to the angles in ang_img with an offset of alpha. 


void FGR_d2d(TBitmapFloat&, TBitmapFloat&, TBitmapFloat&, TBitmapFloat&, float, TBitmapFloat&);
TBitmapFloat* FGR_d2d(TBitmapFloat&, TBitmapFloat&, TBitmapFloat&, TBitmapFloat&, float = 0.0);
void FGR_d2dxx(TBitmapFloat&, TBitmapFloat&, TBitmapFloat&, TBitmapFloat&, float, TBitmapFloat&);
void FGR_d2dxy(TBitmapFloat&, TBitmapFloat&, TBitmapFloat&, TBitmapFloat&, float, TBitmapFloat&);
void FGR_d2dyy(TBitmapFloat&, TBitmapFloat&, TBitmapFloat&, TBitmapFloat&, float, TBitmapFloat&);
TBitmapFloat* FGR_d2dxx(TBitmapFloat&, TBitmapFloat&, TBitmapFloat&, TBitmapFloat&, float = 0.0);
TBitmapFloat* FGR_d2dxy(TBitmapFloat&, TBitmapFloat&, TBitmapFloat&, TBitmapFloat&, float = 0.0);
TBitmapFloat* FGR_d2dyy(TBitmapFloat&, TBitmapFloat&, TBitmapFloat&, TBitmapFloat&, float = 0.0);
These functions do essentially the same for the 2nd order derivative. But we have observed that FGR_d2d orients only dxx. To steer the dxy and dyy use FGR_d2dxy, FGR_d2dyy and for completeness FGR_d2dxx. The correct used linear combinations are : 

dxx^theta = cos^2(theta)*dxx + 2*cos(theta)*sin(theta)*dxy + sin^2(theta)*dyy 
dxy^theta = (1-2*sin^2(theta))*dxy - sin(theta)*cos(theta)*dxx + sin(theta)*cos(theta)*dyy 
dyy^theta = -sin^2(theta)*dxx + 2*cos(theta)*sin(theta)*dxy + cos^2(theta)*dyy



void FGR_d3d(TBitmapFloat&, TBitmapFloat&, TBitmapFloat&, TBitmapFloat&, TBitmapFloat&, float, TBitmapFloat&);
TBitmapFloat* FGR_d3d(TBitmapFloat&, TBitmapFloat&, TBitmapFloat&, TBitmapFloat&, TBitmapFloat&, float = 0.0);
This function orients the 3rd derivative, but only dxxx. Steering of dxxy, dxyy, dyyy needs to be done!

The steering functions take as argument the filtered images. This means that they can be used with TGFilter. It is sufficient to hand the "tgfiltered" images as arguments to the function.

back to top

Automatic scale selection

Author: V. Colin de Verdiere, D. Hall

We have implemented functions that select automatically the smallest intrinsic scale of a Laplacian scale profile, normalise the feature to this scale([Lindeberg98], [Lindeberg98b]). This enables the matching of features invariant to the scale of observation. 


void GenerateSigmaImgTG5(TBitmapFloat* float_img, float sinit, float sfactor, long nb, long step, long from, long to, float GAUSSFAC, TBitmapByte& sigma_img);

This procedure computes the smallest intrinsic scale of all (pixelwise) image features according to the normalised Laplacian scale profile. The scale profile is computed for sigmas 


sigma = sinit*pow(sfactor,x)

within the range sigma_0 = sinit*pow(sfactor, from), sigma_1 = sinit*pow(sfactor, from+step), ..., sigma_n = sinit*pow(sfactor, min(from + k*step, to)). Important: nb >= to+1, else you risk a segmentation fault.

sigma_img contains the index x of the intrinsic scale. To obtain the true sigma, you need to compute sigma = sinit*pow(sfactor,x). We work with followin parameters: 


sinit = 0.5
sfactor = 1.02
nb = 141 // must be < 255 otherwise the index does not fit in a pixel of TBitmapByte 
from = 50
step = 10 // the smallest I have use was 3
to = 140

Once you have obtained the sigma_img, you can extract the normalised local features from any derivative you like. In the theses of Colin de Verdiere, Chomat, and Hall we used the intrinsic scale from the Laplacian profile. The intrinsic scale for edges is particularily unstable. In such cases the intrinsic scale should be computed from the Gardient magnitude ([Chomat00]). But for a first start you can continue to use the Laplacian profile.

For feature extraction call following functions in that order 


TBitmapFloatVector* MultiFilterTG5(filterType f, TBitmapFloat& bmp, float Sinit, float Sfactor, long nb, long freq, long init, long final);
void SelectOptSigmaTG5(TBitmapByte& sigma_img, TBitmapFloatVector& deriv_list, float sinit, float sfactor, long nb, float norm, TBitmapFloat& dst);

MultiFilterTG5 generates a (very large) TBitmapFloatVector that contains the fileterd images of all sigmas within the defined range for( int sig=init, sig < final, sig+=freq). SelectOptSigmaTG5 chooses for every pixel the value from deriv_list according to the index marked in sigma_img and places it in dst. Important: to achieve scale invariance, the derivatives must be normalised by sigma^(order of derivative). This means, norm must be set as follows:

  • Gaussian: norm = 0.0
  • 1st derivative: norm = 1.0
  • 2nd derivative: norm = 2.0
  • 3rd derivative: norm = 3.0
Otherwise, the features are not scale invariant.

Corresponding functions exist for the old implementation of the recursive filters, but due to the problem of the bad approximation of the Laplacian, I suggest to use the above describe functions.

back to top

Corresponding Literature: 


 [VanVliet95]
  author = 	 {I.T. Young and L.J. van Vliet},
  title = 	 {Recursive implementation of the Gaussian filter},
  journal = 	 {Signal Processing},
  year = 	 {1995},
  pages = 	 {139--151},

 [VanVliet98]
  author = 	 {L.J. van Vliet and I.T. Young and P.W. Verbeek},
  title = 	 {Recursive Gaussian Derivative Filters},
  booktitle = 	 {ICPR98},
  pages = 	 {509--514},
  year = 	 {1998},
  month = 	 {August},

 [Freeman91]
  author = 	 {W.T. Freeman and E.H. Adelson},
  title = 	 {The Design and Use of Steerable Filters},
  journal = 	 {Transactions on Pattern Analysis and Machine Intelligence},
  year = 	 {1991},
  volume = 	 {13},
  number = 	 {9},
  pages = 	 {891--906},
  month = 	 {September},

 [Lindeberg98]
  author = 	 {T. Lindeberg},
  title = 	 {Feature Detection with Automatic Scale Selection},
  journal = 	 {IJCV},
  year = 	 1998,
  volume = 	 30,
  number = 	 2,
  pages = 	 {79--116}

 [Lindeberg98b]
  author = 	 {T. Lindeberg},
  title = 	 {Edge detection and ridge detection with automatic scale selection},
  journal = 	 {IJCV},
  year = 	 {1998},
  volume = 	 {30},
  number = 	 {2},

 [Chomat00]
  author = 	 {O. Chomat and V. {Colin de Verdi\`ere} and D. Hall and J.L. Crowley},
  title = 	 {Local Scale Selection for Gaussian Based Description Techniques},
  booktitle = 	 {ECCV00},
  address =      {Dublin, Ireland},
  month =        {June},
  year =         2000,

back to top