Using Gibbs sampling to segment an image

Question

I have been reading some NumPy guides but can't seem to figure it out. My TA told me I should be able to speed up my code by using a NumPy array instead of a for loop in the following segment of code.

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)

Here is the entire code of the method it is in:

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

Answer 1

Since this is homework, I won't give you the exact answer, but here are some hints.

1. == is overloaded in numpy to return an array when you pass in an array. So you can do things like this:

   >>> numpy.arange(5) == 3
   array([False, False, False,  True, False], dtype=bool)
   >>> (numpy.arange(5) == 3) == False
   array([ True,  True,  True, False,  True], dtype=bool)
   

2. You can use boolean arrays to assign to specific locations in an array. For example:

   >>> mostly_true = (numpy.arange(5) == 3) == False
   >>> empty = numpy.zeros(5)
   >>> empty[mostly_true] = 5
   >>> empty
   array([ 5.,  5.,  5.,  0.,  5.])
   

3. You can also negate boolean arrays; together, these facts allow you to conditionally assign values to an array. (numpy.where can be used to do something similar.):

   >>> empty[~mostly_true] = 1
   >>> empty
   array([ 5.,  5.,  5.,  1.,  5.])
   

4. You can then multiply those values by other values:

   >>> empty * numpy.arange(5)
   array([  0.,   5.,  10.,   3.,  20.])
   

5. And many different numpy functions (really ufuncs) provide a reduce method that applies the function along the entire array:

   >>> results = empty * numpy.arange(5)
   >>> numpy.multiply.reduce(results)
   0.0
   

You should be able to completely eliminate that for loop using only the above techniques.