Calculating the numerical gradient of a volume in one direction at a time
The only difference I can see between your gradient3 function and MATLAB's gradient is that the latter returns the horizontal derivative as the first output, and your code returns as "x" derivative the vertical derivative.

Note that MATLAB arrays are stored such that the first index is vertical, and the second index is horizontal. Therefore, the code under case 'x', D(2:k-1,:,:) = (F(3:k,:,:)-F(1:k-2,:,:))/2 computes the vertical difference, which intuitively we'd call y derivative.

Comparing:



F = randn(20,20,20);
[FX,FY,FZ] = gradient(F);
isequal(FX,gradient3(F,'y')) % returns true




By the way, using convn you can speed up the computation of the finite difference derivative by about a factor 2x:



D = convn(F,[0.5,0,-0.5],'same');


This produces a different result at the edges of the image, which might or might not be a problem for you but is typically unimportant.