docs: def get_mog_threshold(image, mask=None, object_fraction = 0.2):
      """Compute a background using a mixture of gaussians
      
      This function finds a suitable
      threshold for the input image Block. It assumes that the pixels in the
      image belong to either a background class or an object class. 'pObject'
      is an initial guess of the prior probability of an object pixel, or
      equivalently, the fraction of the image that is covered by objects.
      Essentially, there are two steps. First, a number of Gaussian
      distributions are estimated to match the distribution of pixel
      intensities in OrigImage. Currently 3 Gaussian distributions are
      fitted, one corresponding to a background class, one corresponding to
      an object class, and one distribution for an intermediate class. The
      distributions are fitted using the Expectation-Maximization (EM)
      algorithm, a procedure referred to as Mixture of Gaussians modeling.
      When the 3 Gaussian distributions have been fitted, it's decided
      whether the intermediate class models background pixels or object
      pixels based on the probability of an object pixel 'pObject' given by
      the user.        
      """
      cropped_image = np.array(image.flat) if mask is None else image[mask]
      pixel_count = np.product(cropped_image.shape)
      max_count   = 512**2 # maximum # of pixels analyzed
      #
      # We need at least 3 pixels to keep from crashing because the highest 
      # and lowest are chopped out below.
      #
      object_fraction = float(object_fraction)
      background_fraction = 1.0-object_fraction
      if pixel_count < 3/min(object_fraction,background_fraction):
          return 1
      if np.max(cropped_image)==np.min(cropped_image):
          return cropped_image[0]
      number_of_classes = 3
      if pixel_count > max_count:
          np.random.seed(0)
          pixel_indices = np.random.permutation(pixel_count)[:max_count]
          cropped_image = cropped_image[pixel_indices]
      # Initialize mean and standard deviations of the three Gaussian
      # distributions by looking at the pixel intensities in the original
      # image and by considering the percentage of the image that is
      # covered by object pixels. Class 1 is the background class and Class
      # 3 is the object class. Class 2 is an intermediate class and we will
      # decide later if it encodes background or object pixels. Also, for
      # robustness the we remove 1% of the smallest and highest intensities
      # in case there are any quantization effects that have resulted in
      # unnaturally many 0:s or 1:s in the image.
      cropped_image.sort()
      one_percent = (np.product(cropped_image.shape) + 99)//100
      cropped_image=cropped_image[one_percent:-one_percent]
      pixel_count = np.product(cropped_image.shape)
      # Guess at the class means for the 3 classes: background,
      # in-between and object
      bg_pixel = cropped_image[int(round(pixel_count * background_fraction/2.0))]
      fg_pixel = cropped_image[int(round(pixel_count * (1-object_fraction/2)))]
      class_mean = np.array([bg_pixel, (bg_pixel+fg_pixel)/2,fg_pixel])
      class_std = np.ones((3,)) * 0.15
      # Initialize prior probabilities of a pixel belonging to each class.
      # The intermediate class steals some probability from the background
      # and object classes.
      class_prob = np.array([3.0/4.0 * background_fraction ,
                                1.0/4.0,
                                3.0/4.0 * object_fraction])
      # Expectation-Maximization algorithm for fitting the three Gaussian
      # distributions/classes to the data. Note, the code below is general
      # and works for any number of classes. Iterate until parameters don't
      # change anymore.
      class_count = np.prod(class_mean.shape)
      #
      # Do a coarse iteration on subsampled data and a fine iteration on the real
      # data
      #
      r = np.random.RandomState()
      r.seed(np.frombuffer(cropped_image[:100].data, np.uint8).tolist())
      for data in (
          r.permutation(cropped_image)[0:(len(cropped_image) // 10)],
          cropped_image):
          delta = 1
          pixel_count = len(data)
          while delta > 0.001:
              old_class_mean = class_mean.copy()
              # Update probabilities of a pixel belonging to the background or
              # object1 or object2
              pixel_class_prob = np.ndarray((pixel_count,class_count))
              for k in range(class_count):
                  norm = scipy.stats.norm(class_mean[k],class_std[k])
                  pixel_class_prob[:,k] = class_prob[k] * norm.pdf(data)
              pixel_class_normalizer = np.sum(pixel_class_prob,1)+.000000000001
              for k in range(class_count):
                  pixel_class_prob[:,k] = pixel_class_prob[:,k] / pixel_class_normalizer
                  # Update parameters in Gaussian distributions
                  class_prob[k] = np.mean(pixel_class_prob[:,k])
                  class_mean[k] = (np.sum(pixel_class_prob[:,k] * data) /
                                   (class_prob[k] * pixel_count))
                  class_std[k] = \
                      math.sqrt(np.sum(pixel_class_prob[:,k] * 
                                          (data-class_mean[k])**2)/
                                (pixel_count * class_prob[k])) + .000001
              delta = np.sum(np.abs(old_class_mean - class_mean))
      # Now the Gaussian distributions are fitted and we can describe the
      # histogram of the pixel intensities as the sum of these Gaussian
      # distributions. To find a threshold we first have to decide if the
      # intermediate class 2 encodes background or object pixels. This is
      # done by choosing the combination of class probabilities "class_prob"
      # that best matches the user input "object_fraction".
      
      # Construct an equally spaced array of values between the background
      # and object mean
      ndivisions = 10000
      level = (np.arange(ndivisions) *
               ((class_mean[2]-class_mean[0]) / ndivisions)
               + class_mean[0])
      class_gaussian = np.ndarray((ndivisions,class_count))
      for k in range(class_count):
          norm = scipy.stats.norm(class_mean[k],class_std[k])
          class_gaussian[:,k] = class_prob[k] * norm.pdf(level)
      if (abs(class_prob[1]+class_prob[2]-object_fraction) <
          abs(class_prob[2]-object_fraction)):
          # classifying the intermediate as object more closely models
          # the user's desired object fraction
          background_distribution = class_gaussian[:,0]
          object_distribution = class_gaussian[:,1]+class_gaussian[:,2]
      else:
          background_distribution = class_gaussian[:,0]+class_gaussian[:,1]
          object_distribution = class_gaussian[:,2]
      # Now, find the threshold at the intersection of the background
      # distribution and the object distribution.
      index = np.argmin(np.abs(background_distribution - object_distribution))
      return level[index]