docs: def get_robust_background_threshold(image, 
                                      mask = None,
                                      lower_outlier_fraction = 0.05,
                                      upper_outlier_fraction = 0.05,
                                      deviations_above_average = 2.0,
                                      average_fn = np.mean,
                                      variance_fn = np.std):
      """Calculate threshold based on mean & standard deviation
         The threshold is calculated by trimming the top and bottom 5% of
         pixels off the image, then calculating the mean and standard deviation
         of the remaining image. The threshold is then set at 2 (empirical
         value) standard deviations above the mean.
  
         image - the image to threshold
         mask - mask of pixels to consider (default = all pixels)
         lower_outlier_fraction - after ordering the pixels by intensity, remove
             the pixels from 0 to len(image) * lower_outlier_fraction from
             the threshold calculation (default = .05).
          upper_outlier_fraction - remove the pixels from 
             len(image) * (1 - upper_outlier_fraction) to len(image) from 
             consideration (default = .05).
          deviations_above_average - calculate the standard deviation or MAD and
             multiply by this number and add to the average to get the final 
             threshold (default = 2)
          average_fn - function used to calculate the average intensity (e.g.
             np.mean, np.median or some sort of mode function). Default = np.mean
          variance_fn - function used to calculate the amount of variance.
                        Default = np.sd
      """ 
  
      cropped_image = np.array(image.flat) if mask is None else image[mask]
      n_pixels = np.product(cropped_image.shape)
      if n_pixels<3:
          return 0
      
      cropped_image.sort()
      if cropped_image[0] == cropped_image[-1]:
          return cropped_image[0]
      low_chop = int(round(n_pixels * lower_outlier_fraction))
      hi_chop = n_pixels - int(round(n_pixels * upper_outlier_fraction))
      im   = cropped_image if low_chop == 0 else cropped_image[low_chop:hi_chop]
      mean = average_fn(im)
      sd   = variance_fn(im)
      return mean+sd*deviations_above_average