image derivatives using gaussian

Hi, this is the function I have for calculating the image derivatives. Please help me understand this code as I am new to this field. If anyone could give me some links to understand this concept, I'll be greatful. some doubts that i have -

1. Why are we using ndgrid here?
2. What are the directions 'x', 'y', 'xx', ('xy', 'yx'), 'yy' here?
3. And how and why does the formula for this gaussian change according to the directions?
4. Why are we using imfilter at the end?

 function D =calc_image_derivatives(I,sigma,direction)

[x,y]=ndgrid(floor(-3*sigma):ceil(3*sigma),floor(-3*sigma):ceil(3*sigma));

switch(direction)
  case 'x'
      DGauss=-(x./(2*pi*sigma^4)).*exp(-(x.^2+y.^2)/(2*sigma^2));
  case 'y'
      DGauss=-(y./(2*pi*sigma^4)).*exp(-(x.^2+y.^2)/(2*sigma^2));
  case 'xx'
      DGauss = 1/(2*pi*sigma^4) * (x.^2/sigma^2 - 1) .* exp(-(x.^2 + y.^2)/(2*sigma^2));
  case {'xy','yx'}
      DGauss = 1/(2*pi*sigma^6) * (x .* y)           .* exp(-(x.^2 + y.^2)/(2*sigma^2));
  case 'yy'
      DGauss = 1/(2*pi*sigma^4) * (y.^2/sigma^2 - 1) .* exp(-(x.^2 + y.^2)/(2*sigma^2));
end

 D = imfilter(I,DGauss,'conv','replicate');