code: def radial_kurtosis(
      dki_params,
      min_kurtosis=-3.0 / 7,
      max_kurtosis=10,
      analytical=True,
  ):
      r"""Compute radial kurtosis (RK) of a diffusion kurtosis tensor [1]_, [2]_
  
      Parameters
      ----------
      dki_params : ndarray (x, y, z, 27) or (n, 27)
      All parameters estimated from the diffusion kurtosis model.
      Parameters are ordered as follows:
          1) Three diffusion tensor's eigenvalues
          2) Three lines of the eigenvector matrix each containing the first,
              second and third coordinates of the eigenvector
          3) Fifteen elements of the kurtosis tensor
      min_kurtosis : float (optional)
          To keep kurtosis values within a plausible biophysical range, radial
          kurtosis values that are smaller than `min_kurtosis` are replaced with
          `min_kurtosis`. Default = -3./7 (theoretical kurtosis limit for regions
          that consist of water confined to spherical pores [3]_)
      max_kurtosis : float (optional)
          To keep kurtosis values within a plausible biophysical range, radial
          kurtosis values that are larger than `max_kurtosis` are replaced with
          `max_kurtosis`. Default = 10
      analytical : bool (optional)
          If True, RK is calculated using its analytical solution, otherwise an
          exact numerical estimator is used (see Notes). Default is set to True.
  
      Returns
      -------
      rk : array
          Calculated RK.
  
      Notes
      -----
      RK is defined as the average of the directional kurtosis perpendicular
      to the fiber's main direction e1 [1]_, [2]_:
  
      .. math::
  
      RK \equiv \frac{1}{2\pi} \int d\Omega _\mathbf{\theta} K(\mathbf{\theta})
                \delta (\mathbf{\theta}\cdot \mathbf{e}_1)
  
      This equation can be numerically computed by averaging apparent
      directional kurtosis samples for directions perpendicular to e1.
  
      Otherwise, RK can be calculated from its analytical solution [2]_:
  
      .. math::
  
          K_{\bot} = G_1(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{2222} +
                     G_1(\lambda_1,\lambda_3,\lambda_2)\hat{W}_{3333} +
                     G_2(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{2233}
  
      where:
  
      .. math::
  
          G_1(\lambda_1,\lambda_2,\lambda_3)=
          \frac{(\lambda_1+\lambda_2+\lambda_3)^2}{18\lambda_2(\lambda_2-
          \lambda_3)} \left (2\lambda_2 +
          \frac{\lambda_3^2-3\lambda_2\lambda_3}{\sqrt{\lambda_2\lambda_3}}
          \right)
  
      and
  
      .. math::
  
          G_2(\lambda_1,\lambda_2,\lambda_3)=
          \frac{(\lambda_1+\lambda_2+\lambda_3)^2}{(\lambda_2-\lambda_3)^2}
          \left ( \frac{\lambda_2+\lambda_3}{\sqrt{\lambda_2\lambda_3}}-2\right )
  
      References
      ----------
      .. [1] Jensen, J.H., Helpern, J.A., 2010. MRI quantification of
             non-Gaussian water diffusion by kurtosis analysis. NMR in
             Biomedicine 23(7): 698-710
      .. [2] Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011.
             Estimation of tensors and tensor-derived measures in diffusional
             kurtosis imaging. Magn Reson Med. 65(3), 823-836
      .. [3] Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging:
             Robust estimation from DW-MRI using homogeneous polynomials.
             Proceedings of the 8th {IEEE} International Symposium on Biomedical
             Imaging: From Nano to Macro, ISBI 2011, 262-265.
             doi: 10.1109/ISBI.2011.5872402
  
      """
      outshape = dki_params.shape[:-1]
      dki_params = dki_params.reshape((-1, dki_params.shape[-1]))
  
      # Split the model parameters to three variable containing the evals,
      # evecs, and kurtosis elements
      evals, evecs, kt = split_dki_param(dki_params)
  
      if analytical:
          # Rotate the kurtosis tensor from the standard Cartesian coordinate
          # system to another coordinate system in which the 3 orthonormal
          # eigenvectors of DT are the base coordinate
          Wyyyy = Wrotate_element(kt, 1, 1, 1, 1, evecs)
          Wzzzz = Wrotate_element(kt, 2, 2, 2, 2, evecs)
          Wyyzz = Wrotate_element(kt, 1, 1, 2, 2, evecs)
  
          # Compute RK
          RK = (
              _G1m(evals[..., 0], evals[..., 1], evals[..., 2]) * Wyyyy
              + _G1m(evals[..., 0], evals[..., 2], evals[..., 1])
              * Wzzzz
              + _G2m(evals[..., 0], evals[..., 1], evals[..., 2])
              * Wyyzz
          )
  
      else:
          # Numerical Solution using 10 perpendicular directions samples
          npa = 10
  
          # Initialize RK
          RK = np.zeros(kt.shape[:-1])
  
          # select relevant voxels to process
          rel_i = _positive_evals(
              evals[..., 0], evals[..., 1], evals[..., 2]
          )
          dki_params = dki_params[rel_i]
          evecs = evecs[rel_i]
          RKi = RK[rel_i]
  
          # loop over all voxels
          KV = np.zeros((dki_params.shape[0], npa))
          for vox in range(len(dki_params)):
              V = perpendicular_directions(
                  np.array(evecs[vox, :, 0]), num=npa, half=True
              )
              sph = dps.Sphere(xyz=V)
              KV[vox, :] = apparent_kurtosis_coef(
                  dki_params[vox], sph, min_kurtosis=min_kurtosis
              )
          RKi = np.mean(KV, axis=-1)
  
          RK[rel_i] = RKi
  
      if min_kurtosis is not None:
          RK = RK.clip(min=min_kurtosis)
  
      if max_kurtosis is not None:
          RK = RK.clip(max=max_kurtosis)
  
      return RK.reshape(outshape)
  
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 13 for calling a global function with example usage: radial_kurtosis(dki_params, min_kurtosis, max_kurtosis, analytical) and returns: RK