function g  = gD(f, sigma, xorder, yorder)
X = [-2.5*sigma:2.5*sigma];
% Zeroth, first and second gaussian derivatives:
gD = [ Gauss1(sigma); ...
       (-X ./ (sigma^2)) .* Gauss1(sigma); ...
       (-1 / (sigma^2)) .* ((1 - (X.^2) / (sigma^2))) .* Gauss1(sigma) ];
%Apply both x and y derivatives to image 
%(note that gD of y has to be transposed)
g = imfilter(f, gD(xorder+1,:), 'conv', 'replicate');
g = imfilter(g, gD(yorder+1,:)', 'conv', 'replicate');