#!/usr/bin/env python
# threshold.py <input file> <output file basename>
import sys, Image

im1 = Image.open(sys.argv[1]).convert("L")
hist = im1.histogram()
sum = 0
for i in hist:
    sum += i

# compute median
median = 0
temp = 0
while temp < sum / 2:
    temp += hist[median]
    median += 1

# compute "extremes"
significant = 0.05
lower = 0
temp = 0
while temp < significant * sum:
    temp += hist[lower]
    lower += 1
while lower > 0 and hist[lower] > 0:
    lower -= 1
higher = 255
temp = 0
while temp < significant * sum:
    temp += hist[higher]
    higher -= 1
while higher < 255 and hist[higher] > 0:
    higher += 1

# determine background
if median - lower < (higher - median) / 2:
    bgmin, fgmin = lower, higher
else:
    bgmin, fgmin = higher, lower

# compute fg/bg "edge"
edge = median + (median - bgmin) / 2

thres1 = edge + (fgmin - edge) * 3 / 16
thres2 = edge + (fgmin - edge) * 3 / 8
thres3 = edge + (fgmin - edge) * 3 / 4
if thres1 > thres3: thres1, thres3 = thres3, thres1

# build new image
im2 = Image.new("P", im1.size)
im2.putpalette([0]*3 + [85]*3 + [170]*3 + [255]*3 + [0]*244*3)
for i in range(0, im1.size[0]):
    for j in range(0, im1.size[1]):
        if im1.getpixel((i,j)) < thres1:
            c = 0
        elif im1.getpixel((i,j)) < thres2:
            c = 1
        elif im1.getpixel((i,j)) < thres3:
            c = 2
        else:
            c = 3
        im2.putpixel((i,j), c)

im2.save(sys.argv[2]+".gif")
im2.convert("RGB").save(sys.argv[2]+".jpg", quality = 100)
