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);
        }