Skip to content

Instantly share code, notes, and snippets.

@jart
Last active Mar 13, 2021
Embed
What would you like to do?
Rational Illumination for PC and TV w/ IBM Scratchpad
#!/bin/sh
fricas -nosman <<EOF
E := [_
-- Standard Illumination Model for Computers_
--_
-- Is defined as a system of linear equations, where negative_
-- colors don't exist and is solved by computing the point at_
-- which they all intersect the one which needs to be defined_
-- as the Planckian locus of the illuminant._
Yw = 1,_
_
-- sRGB chromaticity magnums_
-- See ITU REC BT.701, W3C, Wikipedia, etc._
xr=64/100, yr=33/100,_
xg=30/100, yg=60/100,_
xb=15/100, yb= 6/100,_
_
-- D65 illumination magnums_
-- sRGB and MPEG were designed to assume an ideal correlated color_
-- temperature 11,247°F averaged w/ the average white western male_
-- ability to see radiation for all wavelengths compared with that_
-- of a black body hotter than the sun using 100 y/o measurements._
xw=31271/100000, yw=32902/100000,_
_
-- XYZ/rgb model_
zr = 1 - xr - yr,_
zg = 1 - xg - yg,_
zb = 1 - xb - yb,_
zw = 1 - xw - yw,_
xr = Xr/(Xr+Yr+Zr),_
yr = Yr/(Xr+Yr+Zr),_
zr = Zr/(Xr+Yr+Zr),_
xg = Xg/(Xg+Yg+Zg),_
yg = Yg/(Xg+Yg+Zg),_
zg = Zg/(Xg+Yg+Zg),_
xb = Xb/(Xb+Yb+Zb),_
yb = Yb/(Xb+Yb+Zb),_
zb = Zb/(Xb+Yb+Zb),_
xw = Xw/(Xw+Yw+Zw),_
yw = Yw/(Xw+Yw+Zw),_
zw = Zw/(Xw+Yw+Zw),_
_
-- Derive XYZ to rgb coefficients_
Yw = C00*Xw + C01*Yw + C02*Zw,_
Yw = C00*Xw + C01*Yw + C02*Zw,_
Yw = C00*Xr + C01*Yr + C02*Zr,_
0 = C00*Xg + C01*Yg + C02*Zg,_
0 = C00*Xb + C01*Yb + C02*Zb,_
Yw = C10*Xw + C11*Yw + C12*Zw,_
0 = C10*Xr + C11*Yr + C12*Zr,_
Yw = C10*Xg + C11*Yg + C12*Zg,_
0 = C10*Xb + C11*Yb + C12*Zb,_
Yw = C20*Xw + C21*Yw + C22*Zw,_
0 = C20*Xr + C21*Yr + C22*Zr,_
0 = C20*Xg + C21*Yg + C22*Zg,_
Yw = C20*Xb + C21*Yb + C22*Zb,_
_
-- Derive rgb to XYZ coefficients as C⁻¹_
Cdet = C00 * C11 * C22 + C10 * C21 * C02 + C20 * C01 * C12-_
C20 * C11 * C02 - C10 * C01 * C22 - C00 * C21 * C12,_
Ci00 = (C11 * C22 - C21 * C12) / Cdet,_
Ci01 = (C21 * C02 - C01 * C22) / Cdet,_
Ci02 = (C01 * C12 - C11 * C02) / Cdet,_
Ci10 = (C20 * C12 - C10 * C22) / Cdet,_
Ci11 = (C00 * C22 - C20 * C02) / Cdet,_
Ci12 = (C10 * C02 - C00 * C12) / Cdet,_
Ci20 = (C10 * C21 - C20 * C11) / Cdet,_
Ci21 = (C20 * C01 - C00 * C21) / Cdet,_
Ci22 = (C00 * C11 - C10 * C01) / Cdet_
]
)set output fortran on
)set fortran optlevel 0
)set fortran startindex 0
)set fortran ints2floats off
concat!(solve(E))
EOF
# output of above
Yw=1 # 1.000000000000000
xr=16/25 # .640000000000000
yr=33/100 # .330000000000000
xg=3/10 # .300000000000000
yg=3/5 # .600000000000000
xb=3/20 # .150000000000000
yb=3/50 # .060000000000000
xw=31271/100000 # .312710000000000
yw=16451/50000 # .329020000000000
zr=3/100 # .030000000000000
zg=1/10 # .100000000000000
zb=79/100 # .790000000000000
zw=35827/100000 # .358270000000000
Zr=79184/4096299 # .019330620152484
Yr=871024/4096299 # .212636821677324
Xr=5067776/12288897 # .412386563252992
Zg=4394405/36866691 # .119197163640208
Yg=8788810/12288897 # .715182981841251
Xg=4394405/12288897 # .357591490920625
Zb=70074185/73733382 # .950372587005435
Yb=887015/12288897 # .072180196481425
Xb=4435075/24577794 # .180450491203564
Zw=35827/32902 # 1.088900370798128
Xw=31271/32902 # .950428545377181
C02=-49353/98980 # -.498615881996363
C01=-608687/395920 # -1.537398969488786
C00=641589/197960 # 3.241003232976359
C12=1826061/43944050 # .041554226340085
C11=82435961/43944050 # 1.875929983695176
C10=-42591639/43944050 # -.969224252202517
C22=49353/46685 # 1.057148977187534
C21=-180961/887015 # -.204011206123910
C20=49353/887015 # .055639419851975
Ci00=5067776/12288897 # .412386563252992
Ci01=4394405/12288897 # .357591490920625
Ci02=4435075/24577794 # .180450491203564
Ci10=871024/4096299 # .212636821677324
Ci11=8788810/12288897 # .715182981841251
Ci12=887015/12288897 # .072180196481425
Ci20=79184/4096299 # .019330620152484
Ci21=4394405/36866691 # .119197163640208
Ci22=70074185/73733382 # .950372587005435
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment