3D Numerical differentiation with spline approximation

I have three 3D matrices X, Y, and Z that define a matrix V of the same size over some region. The matrices are regularly spaced. I'm trying to compute the gradient of V. I have read that interpolating and computing derivatives with splines leads to better results than using central differences. For instance, I have worked before with splinefit and ppdiff (http://www.mathworks.com/matlabcentral/fileexchange/13812-splinefit). The problem is that I can't find code to do this in 3D.

Assuming I only want the derivatives at the sampled locations define by the X, Y, and Z matrices, could I do 1D spline approximations for each dimension and compute the partial derivatives that way? I've tried it with a simple 2D Matlab example, and this idea works; however, I wanted to see if it made mathematical sense