function [ gaussian2d ] = gaussian2d( x, y, sigma, order_x, order_y, normalize )
	if (nargin<6)
        normalize = false;
    end        
        
    gaussian2d_base = exp( -(x.*x + y.*y) / (2*sigma^2) );
    
    % we will do the same as ITK. this is just done to be able to
    % compare gaussian results with those coming from ITK.
	if (normalize)
    	gaussian2d_base = sigma^2*gaussian2d_base;
	end    
    
    if (order_x == 0 && order_y ==0)
        scale = 1 / (2*pi*sigma^2);
        gaussian2d = scale .* gaussian2d_base;
    elseif (order_x == 1 && order_y == 0)
        scale = - x / (sigma^4*2*pi);
        gaussian2d = scale .* gaussian2d_base;
    elseif (order_x == 0 && order_y == 1)
        scale = - y / (sigma^4*2*pi);
        gaussian2d = scale .* gaussian2d_base;    
    elseif (order_x == 2 && order_y == 0)
        t1 = - gaussian2d_base / (sigma^4*2*pi);
        t2 = x .* x .* gaussian2d_base / (sigma^6*2*pi);
        gaussian2d = t1 + t2;
    elseif (order_x == 0 && order_y == 2)
        t1 = - gaussian2d_base / (sigma^4*2*pi);
        t2 = y .* y .* gaussian2d_base / (sigma^6*2*pi);
        gaussian2d = t1 + t2;
    elseif (order_x == 1 && order_y == 1)
       gaussian2d = x .* y .* gaussian2d_base / (2*pi*sigma^6);
    end
end
