counting_without_bias.py

"""Counting objects without bias

http://blogs.mathworks.com/steve/2012/12/18/counting-objects-without-bias/
"""

from scipy.misc import imread, imshow
from scipy.ndimage import label
from scipy import asarray, ones, vstack, hstack
from sys import stdout

# ... counting objects

img = "images/rice-bw-2.png"
bw = imread(img, True)
labeled, n = label(bw, ones((3,3)))  # 8-connectivity

# >>> labeled
# array([[0, 0, 0, ..., 0, 0, 0],
#        [0, 0, 0, ..., 0, 0, 0],
#        [0, 0, 0, ..., 0, 0, 0],
#        ...,
#        [0, 0, 0, ..., 0, 0, 0],
#        [0, 0, 0, ..., 0, 0, 0],
#        [0, 0, 0, ..., 0, 0, 0]], dtype=int32)

# >>> n
# 95

# ... calculating biased density

objects_per_unit_area = n * 1.0 / labeled.size

# >>> objects_per_unit_area
# 0.0014495849609375

# ... to have better estimate of object density, we should not count
# objects touching two adjacent edges.

# Unlike Matlab, SciPy doesn't have a function to remove objects on a
# border. Instead we may pay two borders with white pixels, and
# all objects touching those borders will be counted as a single object.

height, width = bw.shape
pad_south = 255 * ones((1, width))
pad_east = 255 * ones((height + 1, 1))
bw_padded = hstack((vstack((bw, pad_south)), pad_east))

# ... counting all bottom-right boundary objects as one:

labeled2, n2 = label(bw_padded, ones((3,3)))
objects_per_unit_area2 = (n2-1) * 1.0 / bw.size

# >>> objects_per_unit_area2
# 0.0012969970703125

stdout.write("objects per unit area (naive):    %.4f\n" % objects_per_unit_area)
stdout.write("objects per unit area (unbiased): %.4f\n" % objects_per_unit_area2)