code: def carlson_rd(x, y, z, errtol=1e-4):
      r"""Compute the Carlson's incomplete elliptic integral of the second kind
  
      Carlson's incomplete elliptic integral of the second kind is defined as:
  
      .. math::
  
          R_D = \frac{3}{2} \int_{0}^{\infty} (t+x)^{-\frac{1}{2}}
          (t+y)^{-\frac{1}{2}}(t+z)  ^{-\frac{3}{2}}
  
      Parameters
      ----------
      x : ndarray
          First independent variable of the integral.
      y : ndarray
          Second independent variable of the integral.
      z : ndarray
          Third independent variable of the integral.
      errtol : float
          Error tolerance. Integral is computed with relative error less in
          magnitude than the defined value
  
      Returns
      -------
      RD : ndarray
          Value of the incomplete second order elliptic integral
  
      Notes
      -----
      x, y, and z have to be nonnegative and at most x or y is zero.
  
      """
      xn = x.copy()
      yn = y.copy()
      zn = z.copy()
      A0 = (xn + yn + 3.0 * zn) / 5.0
      An = A0.copy()
      Q = (errtol / 4.0) ** (-1 / 6.0) * np.max(
          np.abs([An - xn, An - yn, An - zn]), axis=0
      )
      sum_term = np.zeros(x.shape, dtype=x.dtype)
      n = np.zeros(x.shape)
  
      # Convergence has to be done voxel by voxel
      index = ndindex(x.shape)
      for v in index:
          # Convergence condition
          while 4.0 ** (-n[v]) * Q[v] > abs(An[v]):
              xnroot = np.sqrt(xn[v])
              ynroot = np.sqrt(yn[v])
              znroot = np.sqrt(zn[v])
              lamda = xnroot * (ynroot + znroot) + ynroot * znroot
              sum_term[v] = sum_term[v] + 4.0 ** (-n[v]) / (
                  znroot * (zn[v] + lamda)
              )
              n[v] = n[v] + 1
              xn[v] = (xn[v] + lamda) * 0.250
              yn[v] = (yn[v] + lamda) * 0.250
              zn[v] = (zn[v] + lamda) * 0.250
              An[v] = (An[v] + lamda) * 0.250
  
      # post convergence calculation
      X = (A0 - x) / (4.0**n * An)
      Y = (A0 - y) / (4.0**n * An)
      Z = -(X + Y) / 3.0
      E2 = X * Y - 6.0 * Z * Z
      E3 = (3.0 * X * Y - 8.0 * Z * Z) * Z
      E4 = 3.0 * (X * Y - Z * Z) * Z**2.0
      E5 = X * Y * Z**3.0
      RD = (
          4 ** (-n)
          * An ** (-3 / 2.0)
          * (
              1
              - 3 / 14.0 * E2
              + 1 / 6.0 * E3
              + 9 / 88.0 * (E2**2)
              - 3 / 22.0 * E4
              - 9 / 52.0 * E2 * E3
              + 3 / 26.0 * E5
          )
          + 3 * sum_term
      )
  
      return RD
  
imports: import numpy
  import warnings
  import scipy.optimize
  import dipy.core.sphere
  from dipy.reconst.multi_voxel import multi_voxel_fit
  from dipy.reconst.dti import TensorFit, mean_diffusivity, from_lower_triangular, lower_triangular, decompose_tensor, MIN_POSITIVE_SIGNAL, nlls_fit_tensor, restore_fit_tensor
  from dipy.reconst.utils import dki_design_matrix
  from dipy.reconst.recspeed import local_maxima
  from dipy.reconst.base import ReconstModel
  from dipy.core.ndindex import ndindex
  from dipy.core.geometry import sphere2cart, cart2sphere, perpendicular_directions
  from dipy.core.optimize import PositiveDefiniteLeastSquares
  from dipy.data import get_sphere, get_fnames, load_sdp_constraints
  from dipy.reconst.vec_val_sum import vec_val_vect
  from dipy.core.gradients import check_multi_b
label_desc: Use image node_id 3 for calling a global function with example usage: carlson_rd(x, y, z, errtol) and returns: RD