# 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.