from __future__ import division
from math import sin,cos,atan2,sqrt
from random import random

# WARNING: Untested code is never correct.  This code is untested.


def gray_rgb(gray,r,g,b):
    """
    Return a similar RGB with specified intensity.

    Given a gray and rgb values (all [0-1]),
    return a perceptually similar rgb,
    which has an intensity near gray.

    KNOWN BUGS:
    - Intensities are currently sometimes off by > 10%.
      Cause unknown.
    - gray_rgb(gray,0,0,0) throws a divide by zero exception.

    """
    return gamut_sRGB_from_gray_sRGB(gray,r,g,b)    


# Used http://www.colour.org/tc8-05/Docs/colorspace/61966-2-1.pdf
# http://www.w3.org/Graphics/Color/sRGB is slightly different

def linear_RGB_from_sRGB(r,g,b):
    def op(v):
        if v <= 0.04045:
            V = v / 12.92
        else:
            V = ((v + 0.055)/1.055)**2.4
        return V
    return op(r),op(g),op(b)

def sRGB_from_linear_RGB(r,g,b):
    def op(v):
        if v <= 0.0031308:
            V = 12.92 * v
        else:
            V = 1.055 * v**(1.0/2.4) - 0.055
        return V
    return op(r),op(g),op(b)

def XYZ_from_linear_RGB(R,G,B):
    X = 0.4124*R + 0.3576*G + 0.1805*B
    Y = 0.2126*R + 0.7152*G + 0.0722*B
    Z = 0.0193*R + 0.1192*G + 0.9505*B
    return [e*100 for e in (X,Y,Z)]

def linear_RGB_from_XYZ(X,Y,Z):
    X,Y,Z = [e/100. for e in (X,Y,Z)]
    R =  3.2406*X + -1.5372*Y + -0.4986*Z
    G = -0.9689*X +  1.8758*Y +  0.0415*Z
    B =  0.0557*X + -0.2040*Y +  1.0570*Z
    return R,G,B

#

def UV_from_XYZ(X,Y,Z):
    d = (X + (15*Y) + (3*Z))
    U = (4*X) / d
    V = (9*Y) / d
    return U,V

# D65
wx,wy,wz = 0.3127, 0.3290, 0.3583
wX,wY,wZ = 95.05, 100.0, 108.9
wU,wV = UV_from_XYZ(wX,wY,wZ)

def Luv_from_XYZ(X,Y,Z):
    U,V = UV_from_XYZ(X,Y,Z)
    Y_wY = Y / wY
    if Y_wY > 0.008856:
        L = 116 * Y_wY**(1/3) - 16
    else:
        L = 903.3 * Y_wY
    u = 13*L * (U - wU)
    v = 13*L * (V - wV)
    return L,u,v

def XYZ_from_Luv(L,u,v):
    U = u/(13*L) + wU
    V = v/(13*L) + wV
    y = (L + 16) / 116
    if L <= 8:
        Y = wY * L * 0.001107
    else:
        Y = wY * ((L+16)/116)**3
    X = -9*Y*U / ((U - 4)*V - U*V)
    Z = (9*Y - 15*V*Y - V*X) / (3*V)
    return X,Y,Z

def LCHuv_from_Luv(L,u,v):
    C = sqrt(u**2 + v**2)
    H = atan2(v,u)
    return L,C,H

def Luv_from_LCHuv(L,C,H):
    u = C*cos(H)
    v = C*sin(H)
    return L,u,v

#

def XYZ_from_sRGB(r,g,b):
    return XYZ_from_linear_RGB(*linear_RGB_from_sRGB(r,g,b))

def sRGB_from_XYZ(X,Y,Z):
    return sRGB_from_linear_RGB(*linear_RGB_from_XYZ(X,Y,Z))

def LCHuv_from_sRGB(r,g,b):
    return LCHuv_from_Luv(*Luv_from_XYZ(*XYZ_from_sRGB(r,g,b)))

def sRGB_from_LCHuv(L,C,H):
    return sRGB_from_XYZ(*XYZ_from_Luv(*Luv_from_LCHuv(L,C,H)))

#

def gamut_contains_sRGB(r,g,b):
    return (0 <= r and r <= 1 and
            0 <= g and g <= 1 and
            0 <= b and b <= 1)

def gamut_contains_LCHuv(L,C,H):
    return gamut_contains_sRGB(*sRGB_from_LCHuv(L,C,H))

def sRGB_clip_epsilon(r,g,b):
    if gamut_contains_sRGB(r,g,b):
        return r,g,b
    e = 0.0001
    def clip(v):
        if v < 0 and v > -e:
            v = 0
        if v > 1 and v < 1+e:
            v = 1
        return v
    r = clip(r)
    g = clip(g)
    b = clip(b)
    return r,g,b

#

def gamut_sRGB_from_LCHuv(L,C,H):
    rgb = sRGB_from_LCHuv(L,C,H)
    if gamut_contains_sRGB(*rgb):
        return rgb
    rgb_trimmed = sRGB_clip_epsilon(*rgb)
    if gamut_contains_sRGB(*rgb_trimmed):
        return rgb_trimmed
    for C1 in range(int(round(C)),0,-5)+[0]:
        rgb1 = sRGB_from_LCHuv(L,C1,H)
        rgb1_trimmed = sRGB_clip_epsilon(*rgb1)
        if gamut_contains_sRGB(*rgb1_trimmed):
            return rgb1_trimmed
    assert(False)

def gamut_sRGB_from_gray_sRGB(gray,r,g,b):
    assert(0 <= gray and gray <= 1)
    L1 = gray * 100.
    L,C,H = LCHuv_from_sRGB(r,g,b)
    return gamut_sRGB_from_LCHuv(L1,C,H)

# Testing

def test_roundtrip():
    def neq(expected,actual):
        return abs(expected-actual) > 0.001
    def test(r,g,b):
        r1,g1,b1 = sRGB_from_LCHuv(*LCHuv_from_sRGB(r,g,b))
        if neq(r,r1) or neq(g,g1) or neq(b,b1):
            print ""
            print r,g,b
            print r1,g1,b1
    test(0,0,0)
    test(1,1,1)
    test(1,0,0)
    test(1,1,0)
    test(0,1,1)
    for i in range(0,100000):
        r,g,b = random(),random(),random()
        test(r,g,b)
def test_gamut_view_colors():
    # python thisfile.py | sort > tmp.html
    import colorsys
    def hexrgb(r,g,b):
        return "#%02x%02x%02x" % (r*255,g*255,b*255)
    def td(r,g,b):
        h = hexrgb(r,g,b)
        return "<td bgcolor=\"%s\">%s</td>" % (h,h)
    print "<table>"
    print "<tr><th>gray</th><th>result</th><th>given</th>",
    print "<th>alt</th>",
    print "<th>got gray</th><th>gray</th><th>alt gray</th>",
    print "</tr>"
    for i in range(0,1000):
        gray,r,g,b = random(),random(),random(),random()
        r1,g1,b1 = gamut_sRGB_from_gray_sRGB(gray,r,g,b)
        y1 = 0.3*r1 + 0.6*g1 + 0.1*b1

        # Check a simple alternative to all this
        h,l,s = colorsys.rgb_to_hls(r,g,b)
        rx,gx,bx = colorsys.hls_to_rgb(h,gray,s)
        yx = 0.3*rx + 0.6*gx + 0.1*bx
        
        print ("<tr>"+td(gray,gray,gray)+td(r1,g1,b1)+td(r,g,b)+
               td(rx,gx,bx)+
               td(y1,y1,y1)+
               td(gray,gray,gray)+
               td(yx,yx,yx)+
               "</tr>")
    print "</table>"
def test_gamut_gray_match():
    for i in range(0,1000):
        gray,r,g,b = random(),random(),random(),random()
        r1,g1,b1 = gamut_sRGB_from_gray_sRGB(gray,r,g,b)
        Lg,Cg,Hg = LCHuv_from_sRGB(gray,gray,gray)
        L1,C1,H1 = LCHuv_from_sRGB(r1,g1,b1)
        y1 = 0.3*r1 + 0.6*g1 + 0.1*b1
        e = 0.1  #XXX why isn't this better?
        if abs(gray - y1) > e:
            print ""
            print Lg,L1
            print gray, r,g,b
            print y1, r1,g1,b1
        

#test_roundtrip()
#test_gamut_view_colors()
#test_gamut_gray_match()