mix spectral colors

mix_colors.py

#this outputs 4 horizontal gradients in this order:
#spectral weighted geometric mean
#spectral arithmetic mean (normal)
#RGB weighted geometric mean
#RGB normal
#The illuminant is not weighted properly so these are incredibly bright
#also Meng does not constrain to <=1.0 so many RGB values will violate energy conservation
#Curiously the results are better without adapting samples to E before Meng recovery, and using D65 illuminant
#I think Meng would be better if we could intentionally bake the illuminant into the reflectance; choose D65, etc.

import colour
import numpy as np
import matplotlib
import matplotlib.pyplot as plt

#adjust these to mix new color combos
#avoid 0.0 which has infinite darkening power
color_a_sRGB = np.array([0.001,0.001,1])
color_b_sRGB = np.array([1,1,0.001])

color_a_XYZ = colour.sRGB_to_XYZ(color_a_sRGB)
color_b_XYZ = colour.sRGB_to_XYZ(color_b_sRGB)
Y_b = 20.0
L_A = 4.074366543152521
surround = colour.CAM16_VIEWING_CONDITIONS['Average']
#E reference illuminant
#illuminant_xy = colour.ILLUMINANTS['cie_2_1931']['E']
illuminant_spd = colour.ILLUMINANTS_SPDS['D65']
cmfs = colour.CMFS['CIE 1931 2 Degree Standard Observer']
XYZ_w = colour.xy_to_XYZ(colour.ILLUMINANTS['cie_2_1931']['D65'])
#XYZ_wr = colour.xy_to_XYZ(illuminant_xy)

#color_a_XYZ = colour.chromatic_adaptation_VonKries(color_a_XYZ, XYZ_w, XYZ_wr, transform="Bradford")
color_a_refl = colour.recovery.XYZ_to_spectral_Meng2015(color_a_XYZ)

#color_b_XYZ = colour.chromatic_adaptation_VonKries(color_b_XYZ, XYZ_w, XYZ_wr, transform="Bradford")
color_b_refl = colour.recovery.XYZ_to_spectral_Meng2015(color_b_XYZ)

iratio = 0
i = 0
srgb_colors = np.zeros((4,21,3))
while i <= 20:
    ratio = iratio / 100.0
    #weighted geometric mean and normal
    color_c_wgm = color_a_refl**ratio * color_b_refl**(1 - ratio)
    color_c_add = color_a_refl * ratio + color_b_refl * (1 - ratio)
    new_XYZ_wgm = colour.spectral_to_XYZ(color_c_wgm, cmfs, illuminant_spd)
    new_XYZ_add = colour.spectral_to_XYZ(color_c_add, cmfs, illuminant_spd)
    cam16 = colour.XYZ_to_CAM16(new_XYZ_wgm, XYZ_w, L_A, Y_b, surround)
    print("CAM16_WGM is  ", cam16.J, cam16.C, cam16.h)
    cam16_add = colour.XYZ_to_CAM16(new_XYZ_add, XYZ_w, L_A, Y_b, surround)
    print("CAM16_normal is  ", cam16_add.J, cam16_add.C, cam16_add.h)
    sRGB = colour.XYZ_to_sRGB(new_XYZ_wgm, apply_encoding_cctf=False)
    if np.any(sRGB < 0):
        w = - np.min(sRGB)
        sRGB += w
    if not np.all(sRGB==0):
        sRGB /= np.max(sRGB)
    sRGB = colour.oetf(sRGB, 'sRGB')
    srgb_colors[0][i] = [sRGB[0], sRGB[1], sRGB[2]]
    #do the same thing, normal mode
    sRGB = colour.XYZ_to_sRGB(new_XYZ_add, apply_encoding_cctf=False)
    if np.any(sRGB < 0):
        w = - np.min(sRGB)
        sRGB += w
    if not np.all(sRGB==0):
        sRGB /= np.max(sRGB)
    sRGB = colour.oetf(sRGB, 'sRGB')
    srgb_colors[1][i] = [sRGB[0], sRGB[1], sRGB[2]]
    #rgb WGM
    srgb_colors[2][i] = colour.oetf(colour.oetf_reverse(color_a_sRGB, function='sRGB')**ratio * colour.oetf_reverse(color_b_sRGB, function='sRGB')**(1 - ratio), function='sRGB')
    #rgb normal
    srgb_colors[3][i] = colour.oetf(colour.oetf_reverse(color_a_sRGB, function='sRGB') * ratio + colour.oetf_reverse(color_b_sRGB, function='sRGB') * (1 - ratio), function='sRGB')
    i += 1
    iratio += 5
    print(iratio, i)

plt.imshow(srgb_colors)
plt.show()