Skip to content

Instantly share code, notes, and snippets.

@podgorskiy

podgorskiy/wavelength.py

Created October 27, 2020 05:09
Show Gist options
  • Save podgorskiy/99c283773f7cee8e71386aa8ef622fdf to your computer and use it in GitHub Desktop.
Save podgorskiy/99c283773f7cee8e71386aa8ef622fdf to your computer and use it in GitHub Desktop.
Download ZIP
CIE color matching functions approximation
Raw
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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment