function varargout = var_grad_comp(f)
% Slight modification of Matlab's gradient function
% to estimate the variance of the gradient components
% based on the gradient numerical scheme

perm = [2 1]; 
indx = size(f);
loc = cell(1, ndims(f));
for k = 1:ndims(f)
    loc(k) = {1:indx(k)};
end;

for k = 1:2
   [n,p] = size(f);
   h = loc{k}(:);   
   g  = zeros(size(f),class(f)); % case of singleton dimension

   if n > 1
      g(1,:) = (f(2,:) + f(1,:))/(h(2)-h(1)).^2;
      g(n,:) = (f(n,:) + f(n-1,:))/(h(end)-h(end-1)).^2;
   end

   if n > 2
      h = h(3:n) - h(1:n-2);
      g(2:n-1,:) = (f(3:n,:)+f(1:n-2,:))./h(:,ones(p,1)).^2;
   end

   varargout{k} = ipermute(g,[k:2 1:k-1]);

   f = permute(f,perm);
end 

  tmp = varargout{1};
  varargout{1} = varargout{2};
  varargout{2} = tmp;

