code: def mask_brain(
      volume,
      template_name=None,
      mask_threshold=None,
      mask_self=True,
  ):
      """Mask the simulated volume
      This creates a mask specifying the approximate likelihood that a voxel is
      part of the brain. All values are bounded to the range of 0 to 1. An
      appropriate threshold to isolate brain voxels is >0.2. Critically,
      the data that should be used to create a template shouldn't already be
      masked/skull stripped. If it is then it will give in accurate estimates
      of non-brain noise and corrupt estimations of SNR.
  
      Parameters
      ----------
  
      volume : multidimensional array
          Either numpy array of a volume or a tuple describing the dimensions
          of the mask to be created
  
      template_name : str
          What is the path to the template to be loaded? If empty then it
          defaults to an MNI152 grey matter mask. This is ignored if mask_self
          is True.
  
      mask_threshold : float
          What is the threshold (0 -> 1) for including a voxel in the mask? If
          None then the program will try and identify the last wide peak in a
          histogram of the template (assumed to be the brain voxels) and takes
          the minima before that peak as the threshold. Won't work when the
          data is not bimodal.
  
      mask_self : bool or None
          If set to true then it makes a mask from the volume supplied (by
          averaging across time points and changing the range). If it is set
          to false then it will use the template_name as an input.
  
      Returns
      ----------
  
      mask : 3 dimensional array, binary
          The masked brain, thresholded to distinguish brain and non-brain
  
      template : 3 dimensional array, float
          A continuous (0 -> 1) volume describing the likelihood a voxel is in
          the brain. This can be used to contrast the brain and non brain.
  
      """
  
      # If the volume supplied is a 1d array then output a volume of the
      # supplied dimensions
      if len(volume.shape) == 1:
          volume = np.ones(volume)
  
      # Load in the mask
      if mask_self is True:
          mask_raw = volume
      elif template_name is None:
          mfn = resource_stream(
              __name__, "sim_parameters/grey_matter_mask.npy"
          )
          mask_raw = np.load(mfn)
      else:
          mask_raw = np.load(template_name)
  
      # Make the masks 3dremove_baseline
      if len(mask_raw.shape) == 3:
          mask_raw = np.array(mask_raw)
      elif len(mask_raw.shape) == 4 and mask_raw.shape[3] == 1:
          mask_raw = np.array(mask_raw[:, :, :, 0])
      else:
          mask_raw = np.mean(mask_raw, 3)
  
      # Find the max value (so you can calulate these as proportions)
      mask_max = mask_raw.max()
  
      # Make sure the mask values range from 0 to 1 (make out of max of volume
      #  so that this is invertible later)
      mask_raw = mask_raw / mask_max
  
      # If there is only one brain volume then make this a forth dimension
      if len(volume.shape) == 3:
          temp = np.zeros(
              [volume.shape[0], volume.shape[1], volume.shape[2], 1]
          )
          temp[:, :, :, 0] = volume
          volume = temp
  
      # Reshape the mask to be the size as the brain
      brain_dim = volume.shape
      mask_dim = mask_raw.shape
  
      zoom_factor = (
          brain_dim[0] / mask_dim[0],
          brain_dim[1] / mask_dim[1],
          brain_dim[2] / mask_dim[2],
      )
  
      # Scale the mask according to the input brain
      # You might get a warning if the zoom_factor is not an integer but you
      # can safely ignore that.
      template = ndimage.zoom(mask_raw, zoom_factor, order=2)
      template[template < 0] = 0
  
      # If the mask threshold is not supplied then guess it is a minima
      # between the two peaks of the bimodal distribution of voxel activity
      if mask_threshold is None:
          # How many bins on either side of a peak will be compared
          order = 5
  
          # Make the histogram
          template_vector = template.reshape(
              brain_dim[0] * brain_dim[1] * brain_dim[2]
          )
          template_hist = np.histogram(template_vector, 100)
  
          # Zero pad the values
          binval = np.concatenate(
              [np.zeros((order,)), template_hist[0]]
          )
          bins = np.concatenate([np.zeros((order,)), template_hist[1]])
  
          # Identify the first two peaks
          peaks = signal.argrelmax(binval, order=order)[0][0:2]
  
          # What is the minima between peaks
          minima = binval[peaks[0] : peaks[1]].min()
  
          # What is the index of the last idx with this min value (since if
          # zero, there may be many)
          minima_idx = (
              np.where(binval[peaks[0] : peaks[1]] == minima) + peaks[0]
          )[-1]
  
          # Convert the minima into a threshold
          mask_threshold = bins[minima_idx][0]
  
      # Mask the template based on the threshold
      mask = np.zeros(template.shape)
      mask[template > mask_threshold] = 1
  
      return mask, template
  
imports: import logging
  from itertools import product
  from statsmodels.tsa.arima.model import ARIMA
  import math
  import numpy
  from numpy.linalg import LinAlgError
  from pkg_resources import resource_stream
  from scipy import stats
  from scipy import signal
  import scipy.ndimage
  import copy
  from scipy import optimize
label_desc: Use image node_id 22 for calling a global function with example usage: mask_brain(volume, template_name, mask_threshold, mask_) and returns: mask,  template