code: def gradient(f):
      """
      Return the gradient of an N-dimensional array.
  
      The gradient is computed using central differences in the interior
      and first differences at the boundaries. The returned gradient hence has
      the same shape as the input array.
  
      Parameters
      ----------
      f : array_like
        An N-dimensional array containing samples of a scalar function.
  
      Returns
      -------
      gradient : ndarray
        N arrays of the same shape as `f` giving the derivative of `f` with
        respect to each dimension.
  
      Examples
      --------
      >>> x = np.array([1, 2, 4, 7, 11, 16], dtype=float)
      >>> gradient(x)
      array([ 1. ,  1.5,  2.5,  3.5,  4.5,  5. ])
  
      >>> gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=float))
      [array([[ 2.,  2., -1.],
             [ 2.,  2., -1.]]), array([[ 1. ,  2.5,  4. ],
             [ 1. ,  1. ,  1. ]])]
  
      Notes
      -----
      This is a simplified implementation of gradient that is part of numpy
      1.8. In order to mitigate the effects of changes added to this
      implementation in version 1.9 of numpy, we include this implementation
      here.
      """
      f = np.asanyarray(f)
      N = len(f.shape)  # number of dimensions
      dx = [1.0] * N
  
      # use central differences on interior and first differences on endpoints
      outvals = []
  
      # create slice objects --- initially all are [:, :, ..., :]
      slice1 = [slice(None)] * N
      slice2 = [slice(None)] * N
      slice3 = [slice(None)] * N
  
      for axis in range(N):
          # select out appropriate parts for this dimension
          out = np.empty_like(f)
          slice1[axis] = slice(1, -1)
          slice2[axis] = slice(2, None)
          slice3[axis] = slice(None, -2)
  
          # 1D equivalent -- out[1:-1] = (f[2:] - f[:-2])/2.0
          out[tuple(slice1)] = (
              f[tuple(slice2)] - f[tuple(slice3)]
          ) / 2.0
          slice1[axis] = 0
          slice2[axis] = 1
          slice3[axis] = 0
          # 1D equivalent -- out[0] = (f[1] - f[0])
          out[tuple(slice1)] = f[tuple(slice2)] - f[tuple(slice3)]
          slice1[axis] = -1
          slice2[axis] = -1
          slice3[axis] = -2
          # 1D equivalent -- out[-1] = (f[-1] - f[-2])
          out[tuple(slice1)] = f[tuple(slice2)] - f[tuple(slice3)]
  
          # divide by step size
          outvals.append(out / dx[axis])
          # reset the slice object in this dimension to ":"
          slice1[axis] = slice(None)
          slice2[axis] = slice(None)
          slice3[axis] = slice(None)
  
      if N == 1:
          return outvals[0]
      else:
          return outvals
  
imports: import numpy
  import scipy.sparse
  import scipy.linalg
  from dipy.reconst.base import ReconstModel, ReconstFit
  from dipy.tracking.utils import unique_rows
  from dipy.tracking.streamline import transform_streamlines
  from dipy.tracking.vox2track import _voxel2streamline
  import dipy.data
  import dipy.core.optimize
label_desc: Use image node_id 1 for calling a global function with example usage: gradient(f) and returns: outvals, outvals