ImProc ass3: Added initial implementation of back projection.

improc/ass3/back_projection.py

+#!/usr/bin/env python
+from numpy import zeros
+from matplotlib.pyplot import imread, imshow, show
+from intersect import col2bin, domainIterator, colHist, histogramIntersect
+
+def D(image, radius, x, y):
+    """Get the average value of a circle around a pixel."""
+    disk_sum = 0.
+    pixels = 0
+    r_sq = radius ** 2
+    w, h = image.shape[:2]
+
+    # Loop to the square that surrounds the circle, and check if the pixel
+    # is inside the disk
+    for dx in xrange(-radius, radius + 1):
+        for dy in xrange(-radius, radius + 1):
+            cx = x + dx
+            cy = y + dy
+
+            if cx >= 0 and cy >= 0 and cx < w and cy < h \
+                    and dx ** 2 + dy ** 2 < r_sq:
+                disk_sum += image[cx, cy]
+                pixels += 1
+
+    return disk_sum / pixels if pixels else 0
+
+def convolution(image, radius):
+    """Calculate the convolution of an image with a specified disk radius."""
+    c = zeros(image.shape)
+
+    for x, y in domainIterator(image):
+        c[x, y] = D(image, radius, x, y)
+
+    return c
+
+def hbp(image, environment, bins, model, radius):
+    """Create the histogram back projection of two images."""
+    # Create image histograms
+    print 'Creating histograms...'
+    M = colHist(image, bins, 'rgb')
+    I = colHist(environment, bins, 'rgb')
+
+    # Create ratio histogram
+    R = zeros(I.shape)
+
+    for c in domainIterator(R, 3):
+        if (I[c] != 0).all():
+            R[c] = M[c] / I[c]
+
+    # Create back projection
+    print 'Creating back projection...'
+    b = zeros(environment.shape)
+    use = environment * bins
+
+    if model == 'rgb':
+        use /= 255
+    elif model == 'hsv':
+        # TODO: implement HSV color model
+        pass
+
+    for p in domainIterator(b):
+        b[p] = min(R[col2bin(use[p])], 1)
+
+    # Create convolution to create a peak value
+    print 'Creating convolution...'
+    return convolution(b, radius)
+
+def loc(x, y):
+    pass
+
+def find_image(image, environment, bins, model, radius):
+    """Find the location of the peak value of a back projection histogram."""
+    b = hbp(image, environment, bins, model, radius)
+
+    return loc(environment, b.max())
+
+if __name__ == '__main__':
+    print 'Reading images...'
+    waldo = imread('waldo.tiff')
+    env = imread('waldo_env.tiff')
+
+    p = hbp(waldo, env, [8] * 3, 'rgb', 5)
+    print 'done'
+    imshow(p, cmap='gray')
+
+    # Draw a rectangle around the found center pixel
+    #x, y = find_image(waldo, env, [8] * 3, 'rgb')
+    #w, h = waldo.shape[:2]
+    #plot([x - w / 2, x + w / 2], [y - h / 2, y + h / 2], 'r-')
+    #imshow(env)
+
+    show()

improc/ass3/intersect.py

@@ -6,11 +6,17 @@ def col2bin(color):
     """Get the histogram bin coordinates of a color."""
     return tuple(map(lambda x: round(x - 1), color))
 
-def domainIterator(image):
-    """Pixel iterator for colHist."""
-    for x in xrange(image.shape[0]):
-        for y in xrange(image.shape[1]):
-            yield x, y
+def domainIterator(image, dim=2):
+    """Pixel iterator for arrays of with 2 or 3 dimensions."""
+    if dim == 2:
+        for x in xrange(image.shape[0]):
+            for y in xrange(image.shape[1]):
+                yield x, y
+    elif dim == 3:
+        for x in xrange(image.shape[0]):
+            for y in xrange(image.shape[1]):
+                for z in xrange(image.shape[2]):
+                    yield x, y, z
 
 def colHist(image, bins, model):
     """Create the color histogram of an image."""
@@ -20,7 +26,10 @@ def colHist(image, bins, model):
     if model == 'rgb':
         use /= 255
     elif model == 'hsv':
+        # TODO: implement HSV color model
         pass
+    else:
+        raise ValueError('Color model "%s" is not supported.' % model)
 
     for p in domainIterator(image):
         h[col2bin(use[p])] += 1