code: def rumba_deconv_global(
      data,
      kernel,
      mask,
      n_iter=600,
      recon_type="smf",
      n_coils=1,
      R=1,
      use_tv=True,
      verbose=False,
  ):
      r"""
      Fit fODF for all voxels simultaneously using RUMBA-SD.
  
      Deconvolves the kernel from the diffusion-weighted signal at each voxel by
      computing a maximum likelihood estimation of the fODF [1]_. Global fitting
      also permits the use of total variation regularization (RUMBA-SD + TV). The
      spatial dependence introduced by TV promotes smoother solutions (i.e.
      prevents oscillations), while still allowing for sharp discontinuities
      [2]_. This promotes smoothness and continuity along individual tracts while
      preventing smoothing of adjacent tracts.
  
      Generally, global_fit will proceed more quickly than the voxelwise fit
      provided that the computer has adequate RAM (>= 16 GB should be more than
      sufficient).
  
      Parameters
      ----------
      data : 4d ndarray (x, y, z, N)
          Signal values for entire brain. None of the volume dimensions x, y, z
          can be 1 if TV regularization is required.
      kernel : 2d ndarray (N, M)
          Deconvolution kernel mapping volume fractions of the M compartments to
          N-length signal. Last two columns should be for GM and CSF.
      mask : 3d ndarray(x, y, z)
          Binary mask specifying voxels of interest with 1; fODF will only be
          fit at these voxels (0 elsewhere).
      n_iter : int, optional
          Number of iterations for fODF estimation. Must be a positive int.
          Default: 600
      recon_type : {'smf', 'sos'}, optional
          MRI reconstruction method: spatial matched filter (SMF) or
          sum-of-squares (SoS). SMF reconstruction generates Rician noise while
          SoS reconstruction generates Noncentral Chi noise. Default: 'smf'
      n_coils : int, optional
          Number of coils in MRI scanner -- only relevant in SoS reconstruction.
          Must be a positive int. Default: 1
      use_tv : bool, optional
          If true, applies total variation regularization. This requires a brain
          volume with no singleton dimensions. Default: True
      verbose : bool, optional
          If true, logs updates on estimated signal-to-noise ratio after each
          iteration. Default: False
  
      Returns
      -------
      fit_array : 4d ndarray (x, y, z, M)
          fODF and GM/CSF volume fractions computed for each voxel. First M-2
          components are fODF, while last two are GM and CSf respectively.
  
      Notes
      -----
      TV modifies our cost function as follows:
  
      $J(\textbf{f}) = -\log{P(\textbf{S}|\textbf{H}, \textbf{f}, \sigma^2, n)})+
      \alpha_{TV}TV(\textbf{f})$
  
      where the first term is the negative log likelihood described in the notes
      of `rumba_deconv`, and the second term is the TV energy, or the sum of
      gradient absolute values for the fODF across the entire brain. This results
      in a new multiplicative factor in the iterative scheme, now becoming:
  
      $\textbf{f}^{k+1} = \textbf{f}^k \circ \frac{\textbf{H}^T\left[\textbf{S}
      \circ\frac{I_n(\textbf{S}\circ\textbf{Hf}^k/\sigma^2)} {I_{n-1}(\textbf{S}
      \circ\textbf{Hf}^k/\sigma^2)} \right ]} {\textbf{H}^T\textbf{Hf}^k}\circ
      \textbf{R}^k$
  
      where $\textbf{R}^k$ is computed voxelwise by:
  
      $(\textbf{R}^k)_j = \frac{1}{1 - \alpha_{TV}div\left(\frac{\triangledown[
      \textbf{f}^k_{3D}]_j}{\lvert\triangledown[\textbf{f}^k_{3D}]_j \rvert}
      \right)\biggr\rvert_{x, y, z}}$
  
      Here, $\triangledown$ is the symbol for the 3D gradient at any voxel.
  
      The regularization strength, $\alpha_{TV}$ is updated after each iteration
      by the discrepancy principle -- specifically, it is selected to match the
      estimated variance after each iteration [3]_.
  
      References
      ----------
      .. [1] Canales-Rodríguez, E. J., Daducci, A., Sotiropoulos, S. N., Caruyer,
             E., Aja-Fernández, S., Radua, J., Mendizabal, J. M. Y.,
             Iturria-Medina, Y., Melie-García, L., Alemán-Gómez, Y., Thiran,
             J.-P., Sarró, S., Pomarol-Clotet, E., & Salvador, R. (2015).
             Spherical Deconvolution of Multichannel Diffusion MRI Data with
             Non-Gaussian Noise Models and Spatial Regularization. PLOS ONE,
             10(10), e0138910. https://doi.org/10.1371/journal.pone.0138910
  
      .. [2] Rudin, L. I., Osher, S., & Fatemi, E. (1992). Nonlinear total
             variation based noise removal algorithms. Physica D: Nonlinear
             Phenomena, 60(1), 259–268.
             https://doi.org/10.1016/0167-2789(92)90242-F
  
      .. [3] Chambolle A. An algorithm for total variation minimization and
             applications. Journal of Mathematical Imaging and Vision. 2004;
             20:89–97.
      """
  
      # Crop data to reduce memory consumption
      dim_orig = data.shape
      ixmin, ixmax = bounding_box(mask)
      data = crop(data, ixmin, ixmax)
      mask = crop(mask, ixmin, ixmax)
  
      if np.any(np.array(data.shape[:3]) == 1) and use_tv:
          raise ValueError(
              "Cannot use TV regularization if any spatial"
              + "dimensions are 1; "
              + f"provided dimensions were {data.shape[:3]}"
          )
  
      epsilon = 1e-7
  
      n_grad = kernel.shape[0]  # gradient directions
      n_comp = kernel.shape[1]  # number of compartments
      dim = data.shape
      n_v_tot = np.prod(dim[:3])  # total number of voxels
  
      # Initial guess is iso-probable
      fodf0 = np.ones((n_comp, 1), dtype=np.float32)
      fodf0 = fodf0 / np.sum(fodf0, axis=0)
  
      if recon_type == "smf":
          n_order = (
              1  # Rician noise (same as Noncentral Chi with order 1)
          )
      elif recon_type == "sos":
          n_order = n_coils  # Noncentral Chi noise (order = # of coils)
      else:
          raise ValueError(
              "Invalid recon_type. Should be 'smf' or 'sos', "
              + f"received f{recon_type}"
          )
  
      mask_vec = np.ravel(mask)
      # Indices of target voxels
      index_mask = np.atleast_1d(np.squeeze(np.argwhere(mask_vec)))
      n_v_true = len(index_mask)  # number of target voxels
  
      data_2d = np.zeros((n_v_true, n_grad), dtype=np.float32)
      for i in range(n_grad):
          data_2d[:, i] = np.ravel(data[:, :, :, i])[
              index_mask
          ]  # only keep voxels of interest
  
      data_2d = data_2d.T
      fodf = np.tile(fodf0, (1, n_v_true))
      reblurred = np.matmul(kernel, fodf)
  
      # For use later
      kernel_t = kernel.T
      f_zero = 0
  
      # Initialize algorithm parameters
      sigma0 = 1 / 15
      sigma2 = sigma0**2
      tv_lambda = sigma2  # initial guess for TV regularization strength
  
      # Expand into matrix form for iterations
      sigma2 = sigma2 * np.ones(data_2d.shape, dtype=np.float32)
      tv_lambda_aux = np.zeros(n_v_tot, dtype=np.float32)
  
      reblurred_s = data_2d * reblurred / sigma2
  
      for i in range(n_iter):
          fodf_i = fodf
          ratio = mbessel_ratio(n_order, reblurred_s).astype(np.float32)
          rl_factor = np.matmul(kernel_t, data_2d * ratio) / (
              np.matmul(kernel_t, reblurred) + _EPS
          )
  
          if use_tv:  # apply TV regularization
              tv_factor = np.ones(fodf_i.shape, dtype=np.float32)
              fodf_4d = _reshape_2d_4d(fodf_i.T, mask)
              # Compute gradient, divergence
              gr = _grad(fodf_4d)
              d_inv = 1 / np.sqrt(
                  epsilon**2 + np.sum(gr**2, axis=3)
              )
              gr_norm = gr * d_inv[:, :, :, None, :]
              div_f = _divergence(gr_norm)
              g0 = np.abs(1 - tv_lambda * div_f)
              tv_factor_4d = 1 / (g0 + _EPS)
  
              for j in range(n_comp):
                  tv_factor_1d = np.ravel(tv_factor_4d[:, :, :, j])[
                      index_mask
                  ]
                  tv_factor[j, :] = tv_factor_1d
  
              # Apply TV regularization to iteration factor
              rl_factor = rl_factor * tv_factor
  
          fodf = fodf_i * rl_factor  # result of iteration
          fodf = np.maximum(f_zero, fodf)  # positivity constraint
  
          # Update other variables
          reblurred = np.matmul(kernel, fodf)
          reblurred_s = data_2d * reblurred / sigma2
  
          # Iterate variance
          sigma2_i = (1 / (n_grad * n_order)) * np.sum(
              (data_2d**2 + reblurred**2) / 2
              - (sigma2 * reblurred_s) * ratio,
              axis=0,
          )
          sigma2_i = np.minimum(
              (1 / 8) ** 2, np.maximum(sigma2_i, (1 / 80) ** 2)
          )
  
          if verbose:
              logger.info("Iteration %d of %d", i + 1, n_iter)
  
              snr_mean = np.mean(1 / np.sqrt(sigma2_i))
              snr_std = np.std(1 / np.sqrt(sigma2_i))
              logger.info(
                  "Mean SNR (S0/sigma) estimated to be %.3f +/- %.3f",
                  snr_mean,
                  snr_std,
              )
          # Expand into matrix
          sigma2 = np.tile(sigma2_i[None, :], (data_2d.shape[0], 1))
  
          # Update TV regularization strength using the discrepancy principle
          if use_tv:
              if R == 1:
                  tv_lambda = np.mean(sigma2_i)
  
                  if tv_lambda < (1 / 30) ** 2:
                      tv_lambda = (1 / 30) ** 2
              else:  # different factor for each voxel
                  tv_lambda_aux[index_mask] = sigma2_i
                  tv_lambda = np.reshape(tv_lambda_aux, (*dim[:3], 1))
  
      fodf = fodf.astype(np.float64)
      fodf = fodf / (
          np.sum(fodf, axis=0)[None, ...] + _EPS
      )  # normalize fODF
  
      # Extract compartments
      fit_array = np.zeros((*dim_orig[:3], n_comp))
      _reshape_2d_4d(
          fodf.T,
          mask,
          out=fit_array[
              ixmin[0] : ixmax[0],
              ixmin[1] : ixmax[1],
              ixmin[2] : ixmax[2],
          ],
      )
  
      return fit_array
  
imports: import logging
  import warnings
  import numpy
  from dipy.core.geometry import vec2vec_rotmat
  from dipy.core.gradients import gradient_table, unique_bvals_tolerance, get_bval_indices
  from dipy.core.onetime import auto_attr
  from dipy.core.sphere import Sphere
  from dipy.data import get_sphere
  from dipy.reconst.shm import lazy_index, normalize_data
  from dipy.reconst.odf import OdfModel, OdfFit
  from dipy.reconst.csdeconv import AxSymShResponse
  from dipy.segment.mask import bounding_box, crop
  from dipy.sims.voxel import single_tensor, all_tensor_evecs
label_desc: Use image node_id 4 for calling a global function with example usage: rumba_deconv_global(data, kernel, mask, n_iter, recon_type, n_coils, R, use_tv, verbose) and returns: fit_array