import numpy as np
import scipy as sp


def PoissonMixture(I, Mu=None, InitialTau=None, Tol=0.1):
    """Generates a Poisson mixture model to fit pixel intensities for
    foreground/background masking.

    Takes as input an array or intensity image 'I' and optimizes a
    two-component poisson model describing foreground and background intensity
    models. This model can be used to describe the probability that a pixel
    comes from foreground versus background. The poisson distribution assumes
    discrete values and so is suitable for integral valued intensity images.
    Assumes that foreground intensities are lower (darker) than background.

    Parameters
    ----------
    I : array_like
        A hematoxylin intensity image obtained from ColorDeconvolution.
    Mu : double
        Optional mean value of signal to optimize. Calculated from input if
        defined as 'None'. Default value = None.

    Returns
    -------
    Tau : double
        Optimal threshold for distinguishing foreground and background.
    Foreground : array_like
        An intensity image with values in the range [0, 1] representing
        foreground probabiities for each pixel.
    Foreground : array_like
        An intensity image with values in the range [0, 1] representing
        background probabiities for each pixel.

    See Also
    --------
    BinaryMixtureCut

    References
    ----------
    .. [1] Y. Al-Kofahi et al "Improved Automatic Detection and Segmentation
    of Cell Nuclei in Histopathology Images" in IEEE Transactions on Biomedical
    Engineering,vol.57,no.4,pp.847-52, 2010.
    """

    # check if intensity values in 'I' are integer type
    if not np.issubdtype(I.dtype, np.integer):
        raise TypeError('Inputs for Poisson mixture modeling should be integer'
                        ' type')

    # generate a small number for conditioning calculations
    Small = np.finfo(np.float).eps

    # generate histogram of inputs - assume range is 0, 255 (type uint8)
    H = np.histogram(np.ravel(I), bins=256, range=(0, 256))
    X = H[1]
    H = H[0].astype('float') / H[0].sum()

    if(Mu is None):
        Mu = np.dot(X[0: -1], H)

    # calculate cumulative sum along histogram counts
    Cumulative = np.cumsum(H)
    CumProd = np.cumsum(np.multiply(X[0: -1], H))

    # determine cost at each possible threshold
    P0 = Cumulative
    P0[P0 <= 0] = Small
    P1 = 1-Cumulative
    P1[P1 <= 0] = Small
    Mu0 = np.divide(CumProd, P0) + Small
    Mu1 = np.divide(CumProd[-1]-CumProd, P1) + Small
    Cost = Mu - np.multiply(P0, np.log(P0) + np.multiply(Mu0, np.log(Mu0))) -\
        np.multiply(P1, np.log(P1) + np.multiply(Mu1, np.log(Mu1)))

    # identify minimum cost threshold
    Tau = X[np.argmin(Cost)]
    Mu0 = Mu0[np.argmin(Cost)]

    # build probability distribution of foreground intensity values
    Poisson = sp.stats.poisson(Mu0)
    Model = Poisson.pmf(np.arange(0, 256))
    Foreground = Model[I]
    Poisson = sp.stats.poisson(Mu1)
    Model = Poisson.pmf(np.arange(0, 256))
    Background = Model[I]

    # return threshold and foreground probabilities
    return Tau, Foreground, Background

@@ -1175,28 +1175,28 @@
  pixel.%0A    
-Fore
+Back
 ground : arr

Update PoissonMixture.py