dual image script

combo.py

#!/usr/bin/env python2
# dual image script
# requires numpy and PIL
# tested in Python 2.7.13

# Copyright (c) 2017 contributors.
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.

import sys
import numpy as np
from PIL import Image, PngImagePlugin

THUMBSIZE = 250
EPS = np.finfo(np.float32).eps

def bisect(f, target, args=(), prec=12):
    decreasing = (f(0, *args) > f(1, *args))
    shape = np.array(target).shape
    a = np.zeros(shape, np.float32)
    b = np.ones(shape, np.float32)
    for i in xrange(prec):
        x = (a + b)/2
        fx = np.array(f(x, *args))
        comp = (target < fx) ^ decreasing
        a = np.where(comp, a, x)
        b = np.where(comp, x, b)
    return (a + b)/2

def sRGB_to_linear(c):
    return np.where(c <= .04045, c/12.92, ((c+.055)/1.055)**2.4)

if not 4 <= len(sys.argv) <= 9:
    print "combo.py thumb.png hidden.png output.png [resize=100] [gamma=.023] [critdensity=.25] [maxleak1=.001] [maxleak2=.001]"
    exit()
im1f = sys.argv[1]
im2f = sys.argv[2]
im3f = sys.argv[3]
resize = float(sys.argv[4]) if len(sys.argv) > 4 else 100.
resize /= 100
gamma = float(sys.argv[5]) if len(sys.argv) > 5 else .023
ig = 1/gamma
critdensity = float(sys.argv[6]) if len(sys.argv) > 6 else .25
if gamma >= 1:
    critdensity = 1 - critdensity
maxleak1 = float(sys.argv[7]) if len(sys.argv) > 7 else .001
maxleak2 = float(sys.argv[8]) if len(sys.argv) > 8 else .001

im1 = Image.open(im1f).convert("RGB")
im2 = Image.open(im2f).convert("RGB")

size2 = tuple(int(round(min(x)*resize)) for x in zip(im1.size, im2.size))
if im1.size != size2:
    im1 = im1.resize(size2, Image.ANTIALIAS)
if im2.size != size2:
    im2 = im2.resize(size2, Image.ANTIALIAS)

im1d = np.array(im1.getdata(), np.float32)/255
im2d = np.array(im2.getdata(), np.float32)/255

def adjust(scale, fixpoint, imd):
    if gamma >= 1:
        fixpoint = 1 - fixpoint
    return (imd - fixpoint)*scale + fixpoint

def adjust2(scale, fixpoint, imd):
    return [adjust(scale, 1 - fixpoint, imd[0]), adjust(scale, fixpoint, imd[1])]

def leakage1(scale, imd):
    im1a, im2a = adjust2(scale, 0, imd)
    return np.sum(np.maximum(np.sign(1-gamma)*(sRGB_to_linear(im2a) - im1a), 0))/im1a.size

def leakage2(scale, imd):
    im1a, im2a = adjust2(scale, .5, imd)
    return np.sum(np.maximum(np.sign(1-gamma)*(im1a - im2a**gamma), 0))/im1a.size

thumbs = [im.copy() for im in (im1, im2)]
for thumb in thumbs:
    thumb.thumbnail((THUMBSIZE, THUMBSIZE), Image.ANTIALIAS)
thumbs = [np.array(thumb, np.float32)/255 for thumb in thumbs]
scale1 = bisect(leakage1, maxleak1, (thumbs,))
thumbs = adjust2(scale1, 0, thumbs)
thumbs[1] = sRGB_to_linear(thumbs[1])
scale2 = bisect(leakage2, maxleak2, (thumbs,))

im1d, im2d = adjust2(scale1, 0, [im1d, im2d])
im2d = sRGB_to_linear(im2d)
im1d, im2d = adjust2(scale2, .5, [im1d, im2d])

if gamma < 1:
    im1d = im1d.clip(im2d, im2d**gamma)
else:
    im1d = im1d.clip(im2d**gamma, im2d)

inputvals = np.arange(0, 256, dtype=np.float32)/255
def invert(x, outputvals):
    if outputvals[0] > outputvals[-1]:
        return np.interp(x, outputvals[::-1], inputvals[::-1])
    else:
        return np.interp(x, outputvals, inputvals)

imratios_fixwhite = (1-inputvals**ig)/np.maximum(1-inputvals, EPS)
imratios_fixwhite[-1] = ig
vdark_fixwhite = invert((1-im2d)/(1-im1d+EPS), imratios_fixwhite)
if gamma < 1:
    density_fixwhite = (im2d - vdark_fixwhite**ig)/np.maximum(1 - vdark_fixwhite**ig, EPS)
    vdark_fixwhite = (im1d - density_fixwhite)/np.maximum(1-density_fixwhite, EPS)
else:
    density_fixwhite = (im1d - vdark_fixwhite)/np.maximum(1 - vdark_fixwhite, EPS)
    vdark_fixwhite = ((im2d - density_fixwhite)/np.maximum(1-density_fixwhite, EPS))**gamma

vlight_fixblack = (im2d/np.maximum(im1d, EPS))**(gamma/(1-gamma))
density_fixblack = im2d/np.maximum(vlight_fixblack**ig, EPS)

vlightratios_fixdensity = (1 - (1-critdensity)*inputvals)/max(critdensity, EPS)
imratios_fixdensity = (1-critdensity)*np.maximum(inputvals, EPS)**ig + critdensity*np.maximum(vlightratios_fixdensity, EPS)**ig
vdark_fixdensity = invert(im2d/np.maximum(im1d**ig, EPS), imratios_fixdensity) * im1d
if gamma < 1:
    vlight_fixdensity = ((im2d - (1-critdensity)*vdark_fixdensity**ig)/max(critdensity, EPS))**gamma
    vdark_fixdensity = (im1d - critdensity*vlight_fixdensity)/max(1-critdensity, EPS)
else:
    vlight_fixdensity = (im1d - (1-critdensity)*vdark_fixdensity)/max(critdensity, EPS)
    vdark_fixdensity = ((im2d - critdensity*vlight_fixdensity**ig)/max(1-critdensity, EPS))**gamma

density = np.full(im1d.shape, critdensity, np.float32).clip(density_fixwhite, density_fixblack).clip(0, 1)
vdark = np.select([density_fixwhite >= critdensity, density_fixblack <= critdensity, 1], [vdark_fixwhite, 0, vdark_fixdensity])
if gamma < 1:
    vlight = ((im2d - (1-density)*vdark**ig)/np.maximum(density, EPS))**gamma
else:
    vlight = (im1d - (1-density)*vdark)/np.maximum(density, EPS)

dith = Image.new("F", im1.size)
def dither(indata):
    outdata = np.zeros(indata.shape, np.uint8)
    for i in xrange(3):
        dith.putdata(255*indata[:, i])
        outdata[:, i] = np.array(dith.convert("1").getdata())
    return outdata

def ditherround(data):
    data *= 255
    intp = np.floor(data)
    data = intp + (dither(data - intp) >> 7)
    return data/255

if gamma < 1:
    vlight = ditherround(vlight.clip(0, 1))
    density = ((im2d - vdark**ig)/np.maximum(vlight**ig - vdark**ig, EPS)).clip(0, 1)
    vdark = (im1d - density*vlight)/np.maximum(1 - density, EPS)
    vlight = ditherround(vlight.clip(0, 1))
else:
    vdark = ditherround(vdark.clip(0, 1))
    density = ((im2d - vdark**ig)/np.maximum(vlight**ig - vdark**ig, EPS)).clip(0, 1)
    vlight = (im1d - (1-density)*vdark)/np.maximum(density, EPS)
    vlight = ditherround(vlight.clip(0, 1))

im3d = np.where(dither(density), vlight, vdark)

im3c = [Image.new("L", im1.size) for i in xrange(3)]
for i in xrange(3):
    im3c[i].putdata(np.round(im3d[:, i]*255))
im3 = Image.merge("RGB", im3c)

info = PngImagePlugin.PngInfo()
info.add("gAMA", PngImagePlugin.o32(int(round(gamma*1e5))))
im3.save(im3f, "PNG", pnginfo=info)