# 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:

• Please show what the difference is. The code looks to compute the central difference correctly, it does something special at the edges, but you don’t have enough data there anyway. Is that where the differences happen? Aug 11, 2019 at 14:13
• @CrisLuengo, I see that the edges are the same using both methods but the values in middle (where the object is) are different, I've added some illustration pictures of the results of both methods... Aug 14, 2019 at 7:55

## 3 Answers

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.

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.

• thank you sir for your answer but diff() uses different way of calculating the difference or the derivative (gradient() uses the Central Difference method) and fndir() calculate the Directional Derivative not like the gradient() function which calculate the Numerical Derivative. Aug 7, 2019 at 23:14

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);

• Is this not exactly what OP’s code does? Where is the difference? Aug 14, 2019 at 13:41