CIE color matching functions approximation

wavelength.py

import numpy as np
from scipy.optimize import curve_fit
from matplotlib import pyplot as plt


t = np.arange(380.0, 781.0, 5.0)
s = np.asarray([
    [0.0014, 0.0000, 0.0065], [0.0022, 0.0001, 0.0105], [0.0042, 0.0001, 0.0201],
    [0.0076, 0.0002, 0.0362], [0.0143, 0.0004, 0.0679], [0.0232, 0.0006, 0.1102],
    [0.0435, 0.0012, 0.2074], [0.0776, 0.0022, 0.3713], [0.1344, 0.0040, 0.6456],
    [0.2148, 0.0073, 1.0391], [0.2839, 0.0116, 1.3856], [0.3285, 0.0168, 1.6230],
    [0.3483, 0.0230, 1.7471], [0.3481, 0.0298, 1.7826], [0.3362, 0.0380, 1.7721],
    [0.3187, 0.0480, 1.7441], [0.2908, 0.0600, 1.6692], [0.2511, 0.0739, 1.5281],
    [0.1954, 0.0910, 1.2876], [0.1421, 0.1126, 1.0419], [0.0956, 0.1390, 0.8130],
    [0.0580, 0.1693, 0.6162], [0.0320, 0.2080, 0.4652], [0.0147, 0.2586, 0.3533],
    [0.0049, 0.3230, 0.2720], [0.0024, 0.4073, 0.2123], [0.0093, 0.5030, 0.1582],
    [0.0291, 0.6082, 0.1117], [0.0633, 0.7100, 0.0782], [0.1096, 0.7932, 0.0573],
    [0.1655, 0.8620, 0.0422], [0.2257, 0.9149, 0.0298], [0.2904, 0.9540, 0.0203],
    [0.3597, 0.9803, 0.0134], [0.4334, 0.9950, 0.0087], [0.5121, 1.0000, 0.0057],
    [0.5945, 0.9950, 0.0039], [0.6784, 0.9786, 0.0027], [0.7621, 0.9520, 0.0021],
    [0.8425, 0.9154, 0.0018], [0.9163, 0.8700, 0.0017], [0.9786, 0.8163, 0.0014],
    [1.0263, 0.7570, 0.0011], [1.0567, 0.6949, 0.0010], [1.0622, 0.6310, 0.0008],
    [1.0456, 0.5668, 0.0006], [1.0026, 0.5030, 0.0003], [0.9384, 0.4412, 0.0002],
    [0.8544, 0.3810, 0.0002], [0.7514, 0.3210, 0.0001], [0.6424, 0.2650, 0.0000],
    [0.5419, 0.2170, 0.0000], [0.4479, 0.1750, 0.0000], [0.3608, 0.1382, 0.0000],
    [0.2835, 0.1070, 0.0000], [0.2187, 0.0816, 0.0000], [0.1649, 0.0610, 0.0000],
    [0.1212, 0.0446, 0.0000], [0.0874, 0.0320, 0.0000], [0.0636, 0.0232, 0.0000],
    [0.0468, 0.0170, 0.0000], [0.0329, 0.0119, 0.0000], [0.0227, 0.0082, 0.0000],
    [0.0158, 0.0057, 0.0000], [0.0114, 0.0041, 0.0000], [0.0081, 0.0029, 0.0000],
    [0.0058, 0.0021, 0.0000], [0.0041, 0.0015, 0.0000], [0.0029, 0.0010, 0.0000],
    [0.0020, 0.0007, 0.0000], [0.0014, 0.0005, 0.0000], [0.0010, 0.0004, 0.0000],
    [0.0007, 0.0002, 0.0000], [0.0005, 0.0002, 0.0000], [0.0003, 0.0001, 0.0000],
    [0.0002, 0.0001, 0.0000], [0.0002, 0.0001, 0.0000], [0.0001, 0.0000, 0.0000],
    [0.0001, 0.0000, 0.0000], [0.0001, 0.0000, 0.0000], [0.0000, 0.0000, 0.0000]])

D65 = {
    380: 49.428921,
    385: 51.759850,
    390: 54.090779,
    395: 68.025953,
    400: 81.961127,
    405: 86.313504,
    410: 90.665880,
    415: 91.655559,
    420: 92.645238,
    425: 89.319349,
    430: 85.993460,
    435: 95.060474,
    440: 104.127488,
    445: 110.199807,
    450: 116.272126,
    455: 116.704928,
    460: 117.137730,
    465: 115.706664,
    470: 114.275597,
    475: 114.837711,
    480: 115.399824,
    485: 111.895876,
    490: 108.391928,
    495: 108.703479,
    500: 109.015030,
    505: 108.268521,
    510: 107.522013,
    515: 106.057800,
    520: 104.593588,
    525: 106.069007,
    530: 107.544427,
    535: 105.928756,
    540: 104.313086,
    545: 104.157277,
    550: 104.001468,
    555: 102.000734,
    560: 100.000000,
    565: 98.184891,
    570: 96.369783,
    575: 96.119049,
    580: 95.868315,
    585: 92.348821,
    590: 88.829328,
    595: 89.531304,
    600: 90.233281,
    605: 90.059665,
    610: 89.886050,
    615: 88.959839,
    620: 88.033629,
    625: 85.844284,
    630: 83.654939,
    635: 83.904045,
    640: 84.153151,
    645: 82.327113,
    650: 80.501074,
    655: 80.631971,
    660: 80.762868,
    665: 81.835533,
    670: 82.908199,
    675: 80.915868,
    680: 78.923538,
    685: 74.590608,
    690: 70.257678,
    695: 71.238575,
    700: 72.219472,
    705: 73.564167,
    710: 74.908862,
    715: 68.486759,
    720: 62.064656,
    725: 66.225988,
    730: 70.387321,
    735: 73.001133,
    740: 75.614945,
    745: 69.824749,
    750: 64.034552,
    755: 55.396835,
    760: 46.759117,
    765: 57.025699,
    770: 67.292280,
    775: 65.562176,
    780: 63.832072
}


def gauss(x, *p):
    A, mu, sigma = p
    return A * 1.0 / sigma / np.sqrt(2.0 * np.pi) * np.exp(-(x-mu)**2/(2.*sigma**2))


def gauss_mix(x, *p):
    A1, mu1, sigma1, A2, mu2, sigma2 = p
    return gauss(x, A1, mu1, sigma1) + gauss(x, A2, mu2, sigma2)


fig, ax = plt.subplots()


for i in range(3):
    p0 = [1., 450., 50., 1., 600., 50.] if i == 0 else [1., 450., 50.]
    d65s = [D65[x] * y for x, y in zip(t, s[:, i])]

    coeff, var_matrix = curve_fit(gauss_mix if i == 0 else gauss, t, d65s, p0=p0)
    print(*coeff)
    r = (gauss_mix if i == 0 else gauss)(t, *coeff)

    c = ['red', 'green', 'blue'][i]

    ax.plot(t, r, c=c)
    ax.scatter(t, d65s, c=c, label=c, alpha=0.9, edgecolors='none')


ax.set(xlabel='wave length[nm]', ylabel='Intensity',
       title='CIE Color matching functions')

ax.legend()
ax.grid(True)

fig.savefig("test.png")
plt.show()