1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
|
#!/usr/bin/env -S python3
import numpy
from numpy.polynomial import Polynomial, Chebyshev
from matplotlib.pyplot import *
from scipy.signal import windows
from scipy import interpolate
measured = numpy.asarray([
[1.0000, 0.9241, 0.9978, 0.9802],
[0.9666, 0.9209, 0.9966, 0.9774],
[0.9330, 0.9174, 0.9952, 0.9746],
[0.9000, 0.9135, 0.9936, 0.9715],
[0.8666, 0.9098, 0.9916, 0.9680],
[0.8333, 0.9057, 0.9893, 0.9646],
[0.8000, 0.9011, 0.9867, 0.9609],
[0.7666, 0.8968, 0.9834, 0.9569],
[0.7333, 0.8920, 0.9801, 0.9528],
[0.7000, 0.8869, 0.9764, 0.9483],
[0.6666, 0.8813, 0.9724, 0.9433],
[0.6333, 0.8753, 0.9681, 0.9384],
[0.6000, 0.8691, 0.9632, 0.9325],
[0.5666, 0.8622, 0.9580, 0.9264],
[0.5333, 0.8551, 0.9522, 0.9199],
[0.5000, 0.8471, 0.9464, 0.9129],
[0.4666, 0.8370, 0.9393, 0.9051],
[0.4333, 0.8234, 0.9321, 0.8964],
[0.4000, 0.8053, 0.9239, 0.8871],
[0.3666, 0.7829, 0.9143, 0.8763],
[0.3333, 0.7557, 0.9038, 0.8641],
[0.3000, 0.7223, 0.8917, 0.8500],
[0.2666, 0.6789, 0.8777, 0.8295],
[0.2333, 0.6224, 0.8610, 0.7973],
[0.2000, 0.5434, 0.8399, 0.7504],
[0.1666, 0.4629, 0.7992, 0.6823],
[0.1333, 0.3816, 0.7318, 0.5687],
[0.1000, 0.3009, 0.6064, 0.4406],
[0.0666, 0.2207, 0.4530, 0.3140],
[0.0333, 0.1400, 0.2403, 0.1872],
[0.0000, 0.0625, 0.0649, 0.0629],
])
xs = measured[:,0]
xs_ = numpy.linspace(-0.10, 1.15, 100)
tukey = windows.tukey(measured.shape[0]);
int_r = interpolate.PchipInterpolator(numpy.flip(xs), numpy.flip(measured[:,1]), extrapolate=True)
int_g = interpolate.PchipInterpolator(numpy.flip(xs), numpy.flip(measured[:,2]), extrapolate=True)
int_b = interpolate.PchipInterpolator(numpy.flip(xs), numpy.flip(measured[:,3]), extrapolate=True)
off_r = int_r(0)
off_g = int_g(0)
off_b = int_b(0)
k = 1./max( int_r(1)-off_r, int_g(1)-off_g, int_b(1)-off_b )
fit_r = Polynomial.fit(xs_, k*(int_r(xs_)-off_r), deg=16, domain=[0,1], window=[0,1], w=windows.tukey(xs_.shape[0]))
fit_g = Polynomial.fit(xs_, k*(int_g(xs_)-off_g), deg=16, domain=[0,1], window=[0,1], w=windows.tukey(xs_.shape[0]))
fit_b = Polynomial.fit(xs_, k*(int_b(xs_)-off_b), deg=16, domain=[0,1], window=[0,1], w=windows.tukey(xs_.shape[0]))
def c_arr_print(label, arr):
print( label+"[] = {" + ', '.join([f'\n\t{v:e}' if 0 == i%4 else f'{v:e}' for i,v in enumerate(arr)]) + "};\n" )
c_arr_print('GLfloat const coef_r', fit_r.convert().coef[1:])
c_arr_print('GLfloat const coef_g', fit_g.convert().coef[1:])
c_arr_print('GLfloat const coef_b', fit_b.convert().coef[1:])
def polyeval(c, xs):
C_ = numpy.flip(c)
ys = numpy.ones(xs.size)*C_[0]
for c in C_[1:]:
ys = ys*xs + c
return ys*xs
if True:
plot(xs, measured[:,1], '+')
plot(xs_, fit_r(xs_))
plot(xs_, polyeval(fit_r.convert().coef[1:], xs_), '*')
plot(xs, measured[:,2], '+')
plot(xs_, fit_g(xs_))
plot(xs_, polyeval(fit_g.convert().coef[1:], xs_), '*')
plot(xs, measured[:,3], '+')
plot(xs_, fit_b(xs_))
plot(xs_, polyeval(fit_b.convert().coef[1:], xs_), '*')
xlim(0,1)
ylim(0,1)
show()
|
