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

• Sam Mussmann: What is neighbor? What does observation_model(obs,True) return?
• Coolguy123: neighbor is a boolean corresponding to whether the pixel is currently believed to be in the foreground or background. observation_model(obs,True) returns a number corresponding to the probability of observing True given that you are at obs (where True corresponds to the pixel being in the foreground).

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.

• Coolguy123: If it is True then I want to multiply by ALPHA and if it is false I multiply by 1-ALPHA but how can I do this in one line?
• Coolguy123: Nevermind I looked at np.where