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]
def get_mog_threshold ( image , mask = None , object_fraction = 0.2 ) : cropped_image = np . array ( image . flat ) if mask is None else image [ mask ] pixel_count = np . product ( cropped_image . shape ) max_count = 512 ** 2 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 ] 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 ) 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 class_prob = np . array ( [ 3.0 / 4.0 * background_fraction , 1.0 / 4.0 , 3.0 / 4.0 * object_fraction ] ) class_count = np . prod ( class_mean . shape ) 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 ( ) 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 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 ) ) 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 ) ) : 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 ] index = np . argmin ( np . abs ( background_distribution - object_distribution ) ) return level [ index ]
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.