code: def _calc_sfnr(volume,
                 mask,
                 ):
      
  
      # Make a matrix of brain voxels by time
      brain_voxels = volume[mask > 0]
  
      # Take the means of each voxel over time
      mean_voxels = np.nanmean(brain_voxels, 1)
  
      # Detrend (second order polynomial) the voxels over time and then
      # calculate the standard deviation.
      order = 2
      seq = np.linspace(1, brain_voxels.shape[1], brain_voxels.shape[1])
      detrend_poly = np.polyfit(seq, brain_voxels.transpose(), order)
  
      # Detrend for each voxel
      detrend_voxels = np.zeros(brain_voxels.shape)
      for voxel in range(brain_voxels.shape[0]):
          trend = detrend_poly[0, voxel] * seq ** 2 + detrend_poly[1, voxel] * \
                                                     seq + detrend_poly[2, voxel]
          detrend_voxels[voxel, :] = brain_voxels[voxel, :] - trend
  
      std_voxels = np.nanstd(detrend_voxels, 1)
  
      # Calculate the sfnr of all voxels across the brain
      sfnr_voxels = mean_voxels / std_voxels
  
      # Return the average sfnr
      return np.mean(sfnr_voxels)
docs: Calculate the the SFNR of a volume
      Calculates the Signal to Fluctuation Noise Ratio, the mean divided
      by the detrended standard deviation of each brain voxel. Based on
      Friedman and Glover, 2006
  
      Parameters
      ----------
  
      volume : 4d array, float
          Take a volume time series
  
      mask : 3d array, binary
          A binary mask the same size as the volume
  
      Returns
      -------
  
      snr : float 
          The SFNR of the volume