function [ gaussian ] = gaussian( x, sigma, order, normalize )
        if isempty(normalize)
            normalize = false;
        end
        
        gaussian_base = exp(-(x.*x)/(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)
            gaussian_base = sigma*gaussian_base;
        end
            
        if order == 0            
            gaussian = 1/(sigma*sqrt(2*pi)) .* gaussian_base;
        elseif order == 1
            gaussian = -x./(sigma^2*sqrt(2*pi)) .* gaussian_base;
        elseif order == 2
            gaussian = - gaussian_base ./ (sigma^3*sqrt(2*pi)) + x.*x.*gaussian_base / (sigma^5*sqrt(2*pi));
        end
end