Calculating the numerical gradient of a volume in one direction at a time

Question

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:

Answer 1

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.

Answer 2

The gradient by definition is a vector formed by three directional derivatives, ∂F/∂x, ∂F/∂y and ∂F/∂z, so if you only need the z component, at least we have two perspectives, which are, in order of recommendation:

1.- You can calculate it directly with diff(F,z). Check the docs https://www.mathworks.com/help/symbolic/differentiation.html

2.- Or you can calculate the directional derivative in the z direction using fndir(F,Vz), where Vz is a vector in the desired direction, z in this case. Check the docs https://www.mathworks.com/help/curvefit/fndir.html

This perspectives have the advantage of being based on built in, optimized functions.

Answer 3

From MATLAB documentation, subsection Algorithms, the gradient G along the second axis of A is calculated as the central difference:

For j in 2:N-1, where N = size(A,2):

G(:,j) = 0.5*(A(:,j+1) - A(:,j-1));

At the edges where central difference cannot be used, the gradient is calculated by:

G(:,1) = A(:,2) - A(:,1);
G(:,N) = A(:,N) - A(:,N-1);