# PixelMath Python code to implement Histogram Equalization.
# Although Python isn't covered until later in the book than
# Chapter 3, this can be used as a "canned" demonstration
# of the before-and-after effects of the method.
# Updated version of April 17, 2013.
# S. Tanimoto
# This one runs the whole process.
# The old version assumed two copies of the image have been set up,
# one as Source1, the other as Destination.
# First create a monochrome version of the original
sourceWin = pmOpenTwo("Herons256.jpg")
destWin = sourceWin + 1
pmSetFormula("(Red1(x,y)+Green1(x,y)+Blue1(x,y))/3")
pmCompute()
# Compute a histogram of the monochrome.
# We assume that this image is in window number 2.
hist = pmHistogramRed(destWin)
# Now compute a cumulative histogram. We are relying on
# long integers working correctly here, which they seem to in
# this version of the Jython implementation.
cum_hist = hist[:]
sum = 0
for i in range(256):
sum += hist[i]
cum_hist[i] = sum
# Store the total number of pixels in its own scalar variable:
total = cum_hist[-1] # It's the last element of cum. hist.
# Create an array of 256 elements to hold the mapping table.
table = hist[:]
j = 0 # Represents the old cumulative histogram index
# Scan through the range of pixel values, filling in the
# mapping table in a way that will try to keep the cumulative
# histogram of the mapped image as close as possible to the
# ideal (linear) cumulative histogram.
for i in range(256):
ideal_cum_val = total * (i+1)/ 256.0
if j < 255:
while cum_hist[j] <= ideal_cum_val:
table[j] = i
# For debugging:
print "Writing "+str(i)+" into table position "+str(j)
j += 1
if j > 254: break
# Complete the table, even if no pixels in the original image
# have these (high) values:
while j < 256:
table[j] = 255
# For debugging:
print "Writing "+str(i)+" into table position "+str(j)
j += 1
# Create a "mapping table" image as a PixelMath image:
# The width is arbitrary, but it should be bigger than 1, because
# PixelMath can't handle images that narrow, and besides, it is nice
# to have something wide enough to see.
width = 32
newWin = pmNewImage(0, "Mapping Table", width, 256, 0,0,0)
# The following code fills in the pixel values from the mapping table.
for i in range(256):
p = table[i]
for x in range(width):
pmSetPixel(newWin, x, i, p, p, p)
# To use the mapping table, set up the source1 to be the
# monochrome image to be equalized,
# Source2 to be the mapping table image,
# and Destination to be an image the same size as the original.
# Then apply the formula
# Source2(0, Source1(x,y))
pmSetSource2(newWin)
pmSetFormula("S2(0, S1(x,y))")
finalWindow = pmCloneImage(sourceWin)
pmSetDestination(finalWindow)
pmCompute()
# To see how the histogram looks after equalization, open a Histogram
# view for the resulting image.