## Thursday, November 17, 2011

### Gaussian 2d

I needed to compute a 2-dimensional Gaussian distribution which is very common when using Gabor filters.
First google result Custom 2D Gauss  provided a quick solution but upon first look the implementation didn't take advantage of any of matlab's features (i.e matrix manipulation) or included functions so it is a bit slow.

I ended rigging it up a bit making it very fast.

One problem I didn't address is adding a scale for the pixels but that can be fixed inside applications by scaling the standard deviation.

The code in Matlab:
function mat = Gaussian2d(...
gsize, center, sigmax, sigmay, theta, offset, factor)
% GAUSSIAN2d - Produces a 2-dimensional gaussian matrix
%
%   Input
%   gsize   :  Size of the output 'gauss', should be a 1x2 vector
%   center  :  The center position of the gaussian [row col]
%   sigmax  :  Std. dev. in the X direction
%   sigmay  :  Std. dev. in the Y direction
%   theta   :  Rotation in degrees
%   offset  :  Minimum value in output
%   factor  :  Related to maximum value of output, should be
%                    different from zero
%
%   Output
%   mat     : probability vector for each angle bin
%-------------------------------------------------------------------------
%%  Validate input arguments
if ~exist('gsize','var')
error('Error : gsize argument not set');
end

if ~exist('center','var')
center = gsize / 2;
end

if ~exist('sigmax','var')
sigmax = 1;
end

if ~exist('sigmay','var')
sigmay = 1;
end

if ~exist('theta','var')
theta = 0;
end

if ~exist('offset','var')
offset = 0;
end

if ~exist('factor','var')
factor = 1;
end
theta  = (theta/180)*pi;
%-------------------------------------------------------------------------
%%  Compute
[X,Y] = meshgrid(1:gsize(1),1:gsize(2));
Xc  = X - center(1);
Yc  = Y - center(2);
xm  = (Xc)*cos(theta) - (Yc)*sin(theta);
ym  = (Xc)*sin(theta) + (Yc)*cos(theta);
u   = (xm / sigmax).^2 + (ym / sigmay).^2;
mat = offset + factor * exp(-u/2);
end


The code in C# in the OpenCv framework :
        /// Create a rotated 2d Gaussian distribution
/// </summary>
/// <param name="_width">Width of the result image in pixels</param>
/// <param name="_height">Height of the result image in pixels</param>
/// <param name="_centerX">X-coordinate of the center of the distribution in pixels</param>
/// <param name="_centerY">Y-coordinate of the center of the distribution in pixels</param>
/// <param name="_sigmaX">Standard deviation of X-axis in pixels</param>
/// <param name="_sigmaY">Standard deviation of Y-axis in pixels</param>
/// <param name="_theta">Rotation of the distribution in degrees</param>
/// <param name="_offset">Resulting value offset</param>
/// <param name="_factor">Resulting value multiplier</param>
/// <returns></returns>
public static Image<Gray, float> Gaussian2D(
int   _width    = 10,
int   _height   = 10,
float _centerX  = 5,
float _centerY  = 5,
float _sigmaX   = 1.0f,
float _sigmaY   = 1.0f,
float _theta    = 0.0f,
float _offset   = 0.0f,
float _factor   = 1.0f)
{
float[, ,] data = new float[_height, _width, 1];
float xc, yc, xm, ym , tmpX, tmpY, u,
cosTheta = (float)Math.Cos(_theta / 180.0f * Math.PI),
sinTheta = (float)Math.Sin(_theta / 180.0f * Math.PI);
_sigmaX = 1.0f / _sigmaX;
_sigmaY = 1.0f / _sigmaY;
for (int i = 0; i < _width; i++)
{
xc = i - _centerX;
for (int j = 0; j < _height; j++)
{
yc = j - _centerY;
xm = xc * cosTheta - yc * sinTheta;
ym = xc * sinTheta + yc * cosTheta;
tmpX = xm * _sigmaX;
tmpY = ym * _sigmaY;
u = tmpX * tmpX + tmpY * tmpY;
data[j, i, 0] = _offset + _factor * (float)Math.Exp(-u / 2);
}
}
return new Image<Gray,float>(data);
}