1 files changed, 90 insertions, 0 deletions
diff --git a/linpoly.py b/linpoly.py
new file mode 100755
index 0000000..3419a07
--- /dev/null
+++ b/linpoly.py
@@ -0,0 +1,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()