code: def mask_for_response_msmt(
      gtab,
      data,
      roi_center=None,
      roi_radii=10,
      wm_fa_thr=0.7,
      gm_fa_thr=0.2,
      csf_fa_thr=0.1,
      gm_md_thr=0.0007,
      csf_md_thr=0.002,
  ):
      """Computation of masks for multi-shell multi-tissue (msmt) response
          function using FA and MD.
  
      Parameters
      ----------
      gtab : GradientTable
      data : ndarray
          diffusion data (4D)
      roi_center : array-like, (3,)
          Center of ROI in data. If center is None, it is assumed that it is
          the center of the volume with shape `data.shape[:3]`.
      roi_radii : int or array-like, (3,)
          radii of cuboid ROI
      wm_fa_thr : float
          FA threshold for WM.
      gm_fa_thr : float
          FA threshold for GM.
      csf_fa_thr : float
          FA threshold for CSF.
      gm_md_thr : float
          MD threshold for GM.
      csf_md_thr : float
          MD threshold for CSF.
  
      Returns
      -------
      mask_wm : ndarray
          Mask of voxels within the ROI and with FA above the FA threshold
          for WM.
      mask_gm : ndarray
          Mask of voxels within the ROI and with FA below the FA threshold
          for GM and with MD below the MD threshold for GM.
      mask_csf : ndarray
          Mask of voxels within the ROI and with FA below the FA threshold
          for CSF and with MD below the MD threshold for CSF.
  
      Notes
      -----
      In msmt-CSD there is an important pre-processing step: the estimation of
      every tissue's response function. In order to do this, we look for voxels
      corresponding to WM, GM and CSF. This function aims to accomplish that by
      returning a mask of voxels within a ROI and who respect some threshold
      constraints, for each tissue. More precisely, the WM mask must have a FA
      value above a given threshold. The GM mask and CSF mask must have a FA
      below given thresholds and a MD below other thresholds. To get the FA and
      MD, we need to fit a Tensor model to the datasets.
      """
  
      if len(data.shape) < 4:
          msg = """Data must be 4D (3D image + directions). To use a 2D image,
          please reshape it into a (N, N, 1, ndirs) array."""
          raise ValueError(msg)
  
      if isinstance(roi_radii, numbers.Number):
          roi_radii = (roi_radii, roi_radii, roi_radii)
  
      if roi_center is None:
          roi_center = np.array(data.shape[:3]) // 2
  
      roi_radii = _roi_in_volume(
          data.shape, np.asarray(roi_center), np.asarray(roi_radii)
      )
  
      roi_mask = _mask_from_roi(data.shape[:3], roi_center, roi_radii)
  
      list_bvals = unique_bvals_tolerance(gtab.bvals)
      if not np.all(list_bvals <= 1200):
          msg_bvals = """Some b-values are higher than 1200.
          The DTI fit might be affected."""
          warnings.warn(msg_bvals, UserWarning)
  
      ten = TensorModel(gtab)
      tenfit = ten.fit(data, mask=roi_mask)
      fa = fractional_anisotropy(tenfit.evals)
      fa[np.isnan(fa)] = 0
      md = mean_diffusivity(tenfit.evals)
      md[np.isnan(md)] = 0
  
      mask_wm = np.zeros(fa.shape, dtype=np.int64)
      mask_wm[fa > wm_fa_thr] = 1
      mask_wm *= roi_mask
  
      md_mask_gm = np.zeros(md.shape, dtype=np.int64)
      md_mask_gm[(md < gm_md_thr)] = 1
  
      fa_mask_gm = np.zeros(fa.shape, dtype=np.int64)
      fa_mask_gm[(fa < gm_fa_thr) & (fa > 0)] = 1
  
      mask_gm = md_mask_gm * fa_mask_gm
      mask_gm *= roi_mask
  
      md_mask_csf = np.zeros(md.shape, dtype=np.int64)
      md_mask_csf[(md < csf_md_thr) & (md > 0)] = 1
  
      fa_mask_csf = np.zeros(fa.shape, dtype=np.int64)
      fa_mask_csf[(fa < csf_fa_thr) & (fa > 0)] = 1
  
      mask_csf = md_mask_csf * fa_mask_csf
      mask_csf *= roi_mask
  
      msg = """No voxel with a {0} than {1} were found.
      Try a larger roi or a {2} threshold for {3}."""
  
      if np.sum(mask_wm) == 0:
          msg_fa = msg.format(
              "FA higher", str(wm_fa_thr), "lower FA", "WM"
          )
          warnings.warn(msg_fa, UserWarning)
  
      if np.sum(mask_gm) == 0:
          msg_fa = msg.format(
              "FA lower", str(gm_fa_thr), "higher FA", "GM"
          )
          msg_md = msg.format(
              "MD lower", str(gm_md_thr), "higher MD", "GM"
          )
          warnings.warn(msg_fa, UserWarning)
          warnings.warn(msg_md, UserWarning)
  
      if np.sum(mask_csf) == 0:
          msg_fa = msg.format(
              "FA lower", str(csf_fa_thr), "higher FA", "CSF"
          )
          msg_md = msg.format(
              "MD lower", str(csf_md_thr), "higher MD", "CSF"
          )
          warnings.warn(msg_fa, UserWarning)
          warnings.warn(msg_md, UserWarning)
  
      return mask_wm, mask_gm, mask_csf
  
imports: from packaging.version import Version
  import numbers
  import warnings
  import numpy
  from dipy.core import geometry
  from dipy.core.gradients import GradientTable, gradient_table, unique_bvals_tolerance, get_bval_indices
  from dipy.data import default_sphere
  from dipy.reconst import shm
  from dipy.reconst.csdeconv import response_from_mask_ssst
  from dipy.reconst.dti import TensorModel, fractional_anisotropy, mean_diffusivity
  from dipy.reconst.multi_voxel import multi_voxel_fit
  from dipy.reconst.utils import _roi_in_volume, _mask_from_roi
  from dipy.sims.voxel import single_tensor
  from dipy.utils.optpkg import optional_package
label_desc: Use image node_id 6 for calling a global function with example usage: mask_for_response_msmt(gtab, data, roi_center, roi_radii, wm_fa_thr, gm_fa_thr, csf_fa_thr, gm_md_thr, csf_md_thr) and returns: mask_wm,  mask_gm,  mask_csf