Calculating the numerical gradient of a volume in one direction at a time I want to calculate the gradient of an image volume in one direction at a time.

Using the built-in function of Matlab gradient() I can get the ∂F/∂x for an F volume (the differences in the x direction) and I can get all the differences along all other directions let say (∂F/∂y and ∂F/∂z) for 3D volume.



FX = gradient(F)
[FX,FY] = gradient(F)
[FX,FY,FZ,...,FN] = gradient(F)


The point is that I need the gradient in just one direction lets say (along the z direction) and I can't do that using this function (to my knowledge and according to the documentation) plus calculating the gradient in all direction is CPU and memory intensive for large volumes (like my case).

I come across the snippet of code  below that calculate the gradient along one direction at a time, but when I compared the results of this with the built-in function of Matlab I didn't get the same results!



function D = gradient3(F,option)
% Example:
% Fx = gradient3(F,'x');

[k,l,m] = size(F);
D  = zeros(size(F),class(F)); 

switch lower(option)
case 'x'
    % Take forward differences on left and right edges
    D(1,:,:) = (F(2,:,:) - F(1,:,:));
    D(k,:,:) = (F(k,:,:) - F(k-1,:,:));
    % Take centered differences on interior points
    D(2:k-1,:,:) = (F(3:k,:,:)-F(1:k-2,:,:))/2;
case 'y'
    D(:,1,:) = (F(:,2,:) - F(:,1,:));
    D(:,l,:) = (F(:,l,:) - F(:,l-1,:));
    D(:,2:l-1,:) = (F(:,3:l,:)-F(:,1:l-2,:))/2;
case 'z'
    D(:,:,1) = (F(:,:,2) - F(:,:,1));
    D(:,:,m) = (F(:,:,m) - F(:,:,m-1));
    D(:,:,2:m-1) = (F(:,:,3:m)-F(:,:,1:m-2))/2;
otherwise
    disp('Unknown option')
end


Does this code calculate the gradient properly and the one that comes with Matlab calculate it differently? or is there is a better way to do it?

Illustrations of both methods results on the x direction on a volume and displaying a slice from that volume and inspecting some of the same pixels from both images, they both obviously calculate the gradient but the central values aren't the same!

Using Matlab's gradient function gradient(): 


Using the above code:
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.