#!/usr/bin/python

# http://www.imatest.com/docs/lightfall.html

# Copyright (c) 2009 by Reid Priedhorsky, <reid@reidster.net>.
#
# This script is distributed under the terms of the GNU General Public
# License; see http://www.gnu.org/licenses/gpl.txt for more information.


import math
import re
import sys

from PIL import Image

# Camera response gamma.
#
# WARNING: this is a "typical" value taken from the web page above, not
# anything principled.
GAMMA = 0.5

# Luminosity contour step size
LUM_STEP = 0.03333333333

# F-stop contour step size (1/3 stops)
F_STEP = 0.333333333333

# How big is the patch to figure center/corner values?
PATCH_SIZE = 5

# Color map (from colorbrewer.org)
COLORS = ((153, 153, 153),
          (247, 129, 191),
          (166, 86, 40),
          (255, 255, 51),
          (255, 127, 0),
          (152, 78, 163),
          (77, 175, 74),
          (55, 126, 184),
          (228, 26, 28))


def main():
   filename = sys.argv[1]
   im_in = Image.open(filename)
   pix_in = im_in.load()
   (width, height) = im_in.size

   # compute variation between patches in center and corners

   lum_center = lum(avg(pix_in,
                        width/2 - PATCH_SIZE/2, height/2 - PATCH_SIZE/2,
                        PATCH_SIZE, PATCH_SIZE))
   lum_center_r = 1.0

   lum_ul = lum(avg(pix_in, 0, 0,
                    PATCH_SIZE, PATCH_SIZE))
   lum_ul_r = lum_ul / lum_center
   ev_ul = ev(lum_ul, lum_center)

   lum_ur = lum(avg(pix_in, width - PATCH_SIZE, 0,
                    PATCH_SIZE, PATCH_SIZE))
   lum_ur_r = lum_ur / lum_center
   ev_ur = ev(lum_ur, lum_center)

   lum_ll = lum(avg(pix_in, 0, height - PATCH_SIZE,
                    PATCH_SIZE, PATCH_SIZE))
   lum_ll_r = lum_ll / lum_center
   ev_ll = ev(lum_ll, lum_center)

   lum_lr = lum(avg(pix_in, width - PATCH_SIZE, height - PATCH_SIZE,
                    PATCH_SIZE, PATCH_SIZE))
   lum_lr_r = lum_lr / lum_center
   ev_lr = ev(lum_lr, lum_center)

   lum_vary = (max(lum_ul, lum_ur, lum_ll, lum_lr)
               / min(lum_ul, lum_ur, lum_ll, lum_lr))
   ev_vary = (max(ev_ul, ev_ur, ev_ll, ev_lr)
              - min(ev_ul, ev_ur, ev_ll, ev_lr))

   # create luminosity variation visualization

   lum_out = Image.new('RGB', im_in.size)
   pix_out = lum_out.load()
   for x in range(width):
      for y in range(height):
         pix_out[x, y] = color_map(lum(pix_in[x, y]) / lum_center)
   lum_out.save(re.sub(r'\.(png|tiff?)$', '.lum.png', filename), optimize=True)
   

   # create ev variation visualization

   # print

   print '''\
%(filename)s:\
 LUM %(lum_ul_r).2f %(lum_ur_r).2f %(lum_ll_r).2f %(lum_lr_r).2f\
 vary %(lum_vary).2f\
 EV %(ev_ul).2f %(ev_ur).2f %(ev_ll).2f %(ev_lr).2f\
 vary %(ev_vary).2f''' % locals()


def avg(pix, x_min, y_min, width, height):
   r_tot = 0.0
   g_tot = 0.0
   b_tot = 0.0
   for x in range(x_min, x_min + width):
      for y in range(y_min, y_min + width):
         (r, g, b) = pix[x, y]
         r_tot += r
         g_tot += g
         b_tot += b
   pix_ct = width * height
   return (r_tot / pix_ct, g_tot / pix_ct, b_tot / pix_ct)


def color_map(x):
   return COLORS[max(0, min(len(COLORS) - 1, int((1-x)/LUM_STEP)))]


def ev(x, ref):
   'Return the EV loss of x relative to ref. Formula from Imatest.'
   return math.log(x/ref, 1/GAMMA)


def lum((r, g, b)):
   'Return the luminosity of RGB tuple. Formula from Imatest.'
   return (0.30*r + 0.59*g + 0.11*b) / 255





if (__name__ == '__main__'):
   main()