问利用Gibbs采样分割图像

EN

for neighbor in get_neighbors(estimates,i,j):
            pXgivenX_ *= edge_model(True,neighbor)*observation_model(obs,True)
            pX_givenX_  *= edge_model(False,neighbor)*observation_model(obs,False)

def gibbs_segmentation(image, burnin, collect_frequency, n_samples):
"""
Uses Gibbs sampling to segment an image into foreground and background.

Inputs
------
image : a numpy array with the image. Should be Nx x Ny x 3
burnin : Number of iterations to run as 'burn-in' before collecting data
collect_frequency : How many samples in between collected samples
n_samples : how many samples to collect in total

Returns
-------
A distribution of the collected samples: a numpy array with a value between
0 and 1 (inclusive) at every pixel.
"""
(Nx, Ny, _) = image.shape
total_iterations = burnin + (collect_frequency * (n_samples - 1))
pixel_indices = list(itertools.product(xrange(Nx),xrange(Ny)))
# The distribution that you will return
distribution = np.zeros( (Nx, Ny) )

# Initialize binary estimates at every pixel randomly. Your code should
# update this array pixel by pixel at each iteration.
estimates = np.random.random( (Nx, Ny) ) > .5

# PreProcessing
preProObs = {}
for (i,j) in pixel_indices:
        preProObs[(i,j)] = []
        preProObs[(i,j)].append(observation_model(image[i][j],False))
        preProObs[(i,j)].append(observation_model(image[i][j],True))

for iteration in xrange(total_iterations):

    # Loop over entire grid, using a random order for faster convergence
    random.shuffle(pixel_indices)

    for (i,j) in pixel_indices:

        pXgivenX_ = 1
        pX_givenX_ = 1

        for neighbor in get_neighbors(estimates,i,j):
            pXgivenX_ *= edge_model(True,neighbor)*preProObs[(i,j)][1]
            pX_givenX_  *= edge_model(False,neighbor)*preProObs[(i,j)][0]
        estimates[i][j] = np.random.random() > pXgivenX_/(pXgivenX_+pX_givenX_)

    if iteration > burnin and (iteration-burnin)%collect_frequency == 0:
        distribution += estimates

return distribution / n_samples

def edge_model(label1, label2):
    """
    Given the values at two pixels, returns the edge potential between
    those two pixels.

    Hint: there might be a more efficient way to compute this for an array
    of values using numpy!
    """
    if label1 == label2:
        return ALPHA
    else:
        return 1-ALPHA