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

3 comments:

  1. Εγώ ξέρω την τούτη ρε:

    mvnpdf()

    έχει και παράδειγμα στο help της MATLAB.

    ReplyDelete
  2. Ωραία , δεν την ήξερα.
    Το help search της matlab μερικές φορές δε βοηθάει.

    ReplyDelete
  3. Εν θυμούμαι που την ήβρα ρε. Εχρειάστηκα την στην διπλωματική σε κάποια φάση αλλά τελικά χρησιμοποιήσα κάτι άλλο.

    ReplyDelete

Blog Directory Hostgator promo codes
Premium Trick