OfflineValidation/python/DivergingColor.py

0001 import os
0002 import math
0003 import numpy as np
0004
0005 # reference: https://www.kennethmoreland.com/color-maps/ColorMapsExpanded.pdf
0006
0007 # RGB to XYZ matrix
0008 # 0.4124 | 0.2126 | 0.0193
0009 # 0.3576 | 0.7152 | 0.1192
0010 # 0.1805 | 0.0722 | 0.9505
0011
0012 # Inverse
0013 #  3.24063  | -0.968931  |  0.0557101
0014 # -1.53721  |  1.87576   | -0.204021
0015 # -0.498629 |  0.0415175 |  1.057
0016
0017
0018
0019
0020 def fRGB(x):
0021     if x > 0.00313080495356037152: val = math.pow(x, 1./2.4)*1.055 - 0.055
0022     else: val = x*12.92
0023     return val*255.
0024
0025 def fRGBinv(x):
0026     if x > 0.04045: val = math.pow((x + 0.055)/1.055, 2.4)
0027     else: val = (x+0.)/12.92
0028     return val*100
0029
0030 def lRGB2sRGB(r, g, b):
0031     #def rgb_lin(x):
0032     #    if x > 0.00313080495356037152: val = math.pow(x, 1./2.4)*1.055 - 0.055
0033     #    else: val = x*12.92
0034     #    return val*255.
0035     #return rgb_lin((r+0.)/100.), rgb_lin((g+0.)/100.), rgb_lin((b+0.)/100.)
0036     return fRGB((r+0.)/100.), fRGB((g+0.)/100.), fRGB((b+0.)/100.)
0037
0038 def sRGB2lRGB(r, g, b):
0039     #def srgb_lrgb(x):
0040     #    if x > 0.04045: val = math.pow((x + 0.055)/1.055, 2.4)
0041     #    else: val = (x+0.)/12.92
0042     #    return val*100
0043     #return srgb_lrgb((r+0.)/255.), srgb_lrgb((g+0.)/255.), srgb_lrgb((b+0.)/255.)
0044     return fRGBinv((r+0.)/255.), fRGBinv((g+0.)/255.), fRGBinv((b+0.)/255.)
0045
0046
0047 def rgb2xyz(r, g, b):
0048     sr, sg, sb = sRGB2lRGB(r, g, b)
0049     x = 0.4124*sr + 0.3576*sg + 0.1805*sb
0050     y = 0.2126*sr + 0.7152*sg + 0.0722*sb
0051     z = 0.0193*sr + 0.1192*sg + 0.9505*sb
0052     return x, y, z
0053
0054 def xyz2rgb(x, y, z):
0055     r =    3.24063*x -  1.53721*y -  0.498629*z
0056     g =  -0.968931*x +  1.87576*y + 0.0415175*z
0057     b =  0.0557101*x - 0.204021*y +     1.057*z
0058     #m = max(max(r, g), b)
0059     #if m > 1.:
0060     #    r = (r+0.)/(m+0.)
0061     #    g = (g+0.)/(m+0.)
0062     #    b = (b+0.)/(m+0.)
0063     return lRGB2sRGB(r, g, b)
0064
0065
0066
0067 def F(v):
0068     if v > 0.008856: return math.pow(v, 1./3.)
0069     else: return 7.787*v + 16./116.
0070
0071 def Finv(v):
0072     if v > 0.20689270648: return math.pow(v, 3)
0073     else: return (v - 16./116.)/7.787
0074
0075
0076 def xyz2Lab(x, y, z, refW):
0077     xn = refW[0]
0078     yn = refW[1]
0079     zn = refW[2]
0080     #xn, yn, zn = [95.047, 100.0, 108.883]
0081     #def F(v):
0082     #    if v > 0.008856: return math.pow(v, 1./3.)
0083     #    else: return 7.787*v + 16./116.
0084     L = 116*(F((y+0.)/(yn +0.)) - 16./116.)
0085     a = 500*(F((x+0.)/(xn+0.)) - F((y+0.)/(yn+0.)))
0086     b = 200*(F((y+0.)/(yn+0.)) - F((z+0.)/(zn+0.)))
0087     return L, a, b
0088
0089 def Lab2xyz(L, a, b, refW):
0090     xn = refW[0]
0091     yn = refW[1]
0092     zn = refW[2]
0093     #xn, yn, zn = [95.047, 100.0, 108.883]
0094     #def Finv(v):
0095     #    if v > 0.20689270648: return math.pow(v, 3)
0096     #    else: return (v - 16./116.)/7.787
0097     x = Finv((a+0.)/500. + (L + 16.)/116.)*xn
0098     y = Finv((L + 16.)/116.)*yn
0099     z = Finv((L + 16.)/116. - (b+0.)/200.)*zn
0100     return x, y, z
0101
0102 def Lab2Msh(L, a, b):
0103     M = math.sqrt(math.pow(L,2) + math.pow(a,2) + math.pow(b,2))
0104     s = math.acos((L+0.)/(M+0.))
0105     h = math.atan2(b,a)
0106     return M, s, h
0107
0108 def Msh2Lab(M, s, h):
0109     L = M*math.cos(s)
0110     a = M*math.sin(s)*math.cos(h)
0111     b = M*math.sin(s)*math.sin(h)
0112     return L, a, b
0113
0114 def rgb2Msh(r, g, b, refW):
0115     x, y, z = rgb2xyz(r, g, b)
0116     xr, yr, zr = rgb2xyz(refW[0], refW[1], refW[2])
0117     L, a, b = xyz2Lab(x, y, z, [xr, yr, zr])
0118     return Lab2Msh(L, a, b)
0119
0120 def Msh2rgb(M, s, h, refW):
0121     xr, yr, zr = rgb2xyz(refW[0], refW[1], refW[2])
0122     L, a, b = Msh2Lab(M, s, h)
0123     x, y, z = Lab2xyz(L, a, b, [xr, yr, zr])
0124     return xyz2rgb(x, y, z)
0125
0126 def AdjustHue(Ms, ss, hs, Mu):
0127     #print('Adjusting Hue')
0128     if Ms >= Mu: return hs
0129     h = ss*math.sqrt(math.pow(Mu, 2.) - math.pow(Ms, 2.))/(Ms*math.sin(ss)+0.)
0130     if hs > -math.pi/3.: return hs + h
0131     else: return hs - h
0132
0133 def radDiff(a1, a2):
0134     diff = abs(a1 - a2)
0135     if diff > math.pi: return abs(diff - 2*math.pi)
0136     else: return diff
0137
0138 #('red: ', (117.34353643868656, 1.099939672641069, 0.698178814103516))
0139 #('blue: ', (137.64998152940237, 1.333915268336423, -0.9374394027523394))
0140
0141 def DivergingColor(col1, col2, white, frac):
0142     M1, s1, h1 = rgb2Msh(col1[0], col1[1], col1[2], white)
0143     M2, s2, h2 = rgb2Msh(col2[0], col2[1], col2[2], white)
0144
0145     #if s1 > 0.05 and s2 > 0.05 and radDiff(h1,h2) > math.pi/3.:
0146     if s1 > 0.05 and s2 > 0.05 and abs(h1 - h2) > math.pi/3.:
0147         Mmid = max(max(M1,M2),88.)
0148         if frac < .5:
0149             M2 = Mmid
0150             s2 = 0.
0151             h2 = 0.
0152             frac = 2*frac
0153         else:
0154             M1 = Mmid
0155             s1 = 0.
0156             h1 = 0.
0157             frac = 2*frac - 1
0158     if s1 < 0.05 and s2 > 0.05: h1 = AdjustHue(M2, s2, h2, M1)
0159     elif s2 < 0.05 and s1 > 0.05: h2 = AdjustHue(M1, s1, h1, M2)
0160
0161     M = (1 - frac)*M1 + frac*M2
0162     s = (1 - frac)*s1 + frac*s2
0163     h = (1 - frac)*h1 + frac*h2
0164     #print('temp', M, s, h, h1, h2, frac)
0165     return Msh2rgb(M, s, h, white)
0166
0167
0168 if __name__ == '__main__':
0169     Msh = [83.9912098 ,  0.54009147, -0.18776355]
0170     Msh_np = np.array(Msh)
0171     #red = [243.59789395465015, 146.5213165050506, 192.51678151291404]
0172     red   = [59, 76, 192]
0173     blue  = [180, 4, 38]
0174     white = [1, 1, 1]
0175     frac = 0.75
0176     print('my val: ', DivergingColor(blue, red, white, frac))
0177     print(xyz2rgb(95.047, 100.0, 108.883))
0178