Otsu's method of thresholding. Partially based on the implementation shown at: http://www.labbookpages.co.uk/software/imgProc/otsuThreshold.html

otsu.py

#!/usr/local/bin/python

#######################################################
# Otsu's Method
# Author: Samuel Jackson ([email protected])
# Date: 21/07/2013
# Description: Performs Otsu's method of thresholding
# using the between class variance.
#######################################################

def otsu(hist, totalPixels):
	currentMax = 0
	threshold = 0
	sumTotal, sumForeground, sumBackground = 0, 0, 0
	weightBackground, weightForeground = 0, 0

	# Calculate the total of the data
	for i,t in enumerate(hist): sumTotal += i * hist[i] 

	for i,t in enumerate(hist):

		# Calculate the weight of the background
		weightBackground += hist[i]
		if( weightBackground == 0 ): continue

		# Calculate the weight of the foreground
		weightForeground = totalPixels - weightBackground
		if ( weightForeground == 0 ): break

		sumBackground += i*hist[i]

		# Calculate the mean of the classes
		meanB = sumBackground / weightBackground
		meanF = (sumTotal - sumBackground) / weightForeground

		# Calculate variance between classes
		varBetween = weightBackground*weightForeground
		varBetween *= (meanB-meanF)*(meanB-meanF)

		# Check if the variance between classes is greater than
		# the current best
		if(varBetween > currentMax):
			currentMax = varBetween
			threshold = i

	return threshold

# Otsu demonstration code	
if __name__ == "__main__":
	import matplotlib.pyplot as plt
	import matplotlib.image as mpimg
	import numpy as np
	import sys

	#Load image and create histogram
	image = mpimg.imread(sys.argv[1])
	img = image.flatten()
	hist, bins = np.histogram(img, bins=256, range=(0,255))

	#Calculate thresholding using otsu
	threshold = otsu(hist,len(img))

	#show histogram of image data
	xs = np.arange(len(hist))
	plt.bar(xs, hist, color="black")
	plt.plot([threshold, threshold],[0, plt.ylim()[1]],linewidth=2.0, color='r')
	plt.show()

	#show converted image using threshold
	plt.imshow(image > threshold)
	plt.show()