/pliant/graphic/color/spectrum.pli
1	abstract
2	[Color spectrum handling functions]
3
4	module "/pliant/language/unsafe.pli"
5	module "/pliant/language/stream.pli"
6	module "/pliant/math/functions.pli"
7
8	constant X_spec "from 380 to 775 step 5 spectrum 0.0014 0.0022 0.0042 0.0076 0.0143 0.0232 0.0435 0.0776 0.1344 0.2148 0.2839 0.3285 0.3483 0.3481 0.3362 0.3187 0.2908 0.2511 0.1954 0.1421 0.0956 0.058 0.032 0.0147 0.0049 0.0024 0.0093 0.0291 0.0633 0.1096 0.1655 0.2257 0.2904 0.3597 0.4334 0.5121 0.5945 0.6784 0.7621 0.8425 0.9163 0.9786 1.0263 1.0567 1.0622 1.0456 1.0026 0.9384 0.8544 0.7514 0.6424 0.5419 0.4479 0.3608 0.2835 0.2187 0.1649 0.1212 0.0874 0.0636 0.0468 0.0329 0.0227 0.0158 0.0114 0.0081 0.0058 0.0041 0.0029 0.002 0.0014 0.001 0.0007 0.0005 0.0003 0.0002 0.0002 0.0001 0.0001 0.0001"
9	constant Y_spec "from 385 to 760 step 5 spectrum 0.0001 0.0001 0.0002 0.0004 0.0006 0.0012 0.0022 0.004 0.0073 0.0116 0.0168 0.023 0.0298 0.038 0.048 0.06 0.0739 0.091 0.1126 0.139 0.1693 0.208 0.2586 0.323 0.4073 0.503 0.6082 0.71 0.7932 0.862 0.9149 0.954 0.9808 0.995 1 0.995 0.9786 0.952 0.9154 0.87 0.8163 0.757 0.6949 0.631 0.5668 0.503 0.4412 0.381 0.321 0.265 0.217 0.175 0.1382 0.107 0.0816 0.061 0.0446 0.032 0.0232 0.017 0.0119 0.0082 0.0057 0.0041 0.0029 0.0021 0.0015 0.001 0.0007 0.0005 0.0004 0.0002 0.0002 0.0001 0.0001 0.0001"
10	constant Z_spec "from 380 to 625 step 5 spectrum 0.0065 0.0105 0.0201 0.0362 0.0679 0.1102 0.2074 0.3713 0.6456 1.0391 1.3856 1.623 1.7471 1.7826 1.7721 1.7441 1.6692 1.5281 1.2876 1.0419 0.813 0.6162 0.4652 0.3533 0.272 0.2123 0.1582 0.1117 0.0782 0.0573 0.0422 0.0298 0.0203 0.0134 0.0087 0.0057 0.0039 0.0027 0.0021 0.0018 0.0017 0.0014 0.0011 0.001 0.0008 0.0006 0.0003 0.0002 0.0002 0.0001"
11	constant D50_spec "from 300 to 830 step 5 spectrum 0.02 1.03 2.05 4.91 7.78 11.26 14.75 16.35 17.95 19.48 21.01 22.48 23.94 25.45 26.96 25.72 24.49 27.18 29.87 39.59 49.31 52.91 56.51 58.27 60.03 58.93 57.82 66.32 74.82 81.04 87.25 88.93 90.61 90.99 91.37 93.24 95.11 93.54 91.96 93.84 95.72 96.17 96.61 96.87 97.13 99.61 102.1 101.43 100.75 101.54 102.32 101.16 100 98.87 97.74 98.33 98.92 96.21 93.5 95.59 97.69 98.48 99.27 99.16 99.04 97.38 95.72 97.29 98.86 97.26 95.67 96.93 98.19 100.6 103 101.07 99.13 93.26 87.38 89.49 91.6 92.25 92.89 84.87 76.85 81.68 86.51 89.55 92.58 85.4 78.23 67.96 57.69 70.31 82.92 80.6 78.27 78.91 79.55 76.48 73.4 68.66 63.92 67.35 70.78 72.61 74.44"
12	# constant D55_spec "from 400 to 700 step 10 spectrum 60.95 68.55 71.58 67.91 85.61 97.99 100.46 99.91 102.74 98.08 100.68 100.70 99.99 104.21 102.10 102.97 100.00 97.22 97.75 91.43 94.42 95.14 94.22 90.45 92.33 88.85 90.32 93.95 89.96 79.68 82.84"
13	constant D55_spec "from 360 to 830 step 10 spectrum 30.5985 34.2841 32.5540 38.0560 60.9021 68.5057 71.5311 67.8734 85.5608 97.9481 100.4224 99.8776 102.7072 98.0524 100.6593 100.6783 99.9744 104.2004 102.0963 102.9651 100.0000 97.2186 97.7535 91.4415 94.4348 95.1597 94.2431 90.4725 92.3609 88.8867 90.3550 93.9930 90.0006 79.7135 82.8818 84.8819 70.2666 79.3356 85.0292 71.9107 52.8168 75.9602 71.8484 72.9666 67.3762 58.7513 65.0203 68.3353"
14	constant D65_spec "from 300 to 830 step 5 spectrum  300:0.03 1.66 3.29 11.77 20.24  325:28.64 37.05 38.50  39.95 42.43  350:44.91 45.78 46.64 49.36 52.09  375:51.03 49.98 52.31 54.65 68.70  400:82.75 87.12 91.49 92.46 93.43  425:90.06 86.68 95.77 104.86 110.94  450:117.01 117.41 117.81 116.34 114.86  475:115.39 115.92 112.37 108.81 109.08  500:109.35 108.58 107.80 106.30 104.79  525:106.24 107.69 106.05 104.41 104.23  550:104.05 102.02 100.00 98.17 96.33  575:96.06 95.79 92.24 88.69 89.35  600:90.01 89.80 89.60 88.65 87.70  625:85.49 83.29 83.49 83.70 81.86  650:80.03 80.12 80.21 81.25 82.28  675:80.28 78.28 74.00 69.72 70.67  700:71.61 72.98 74.35 67.98 61.60  725:65.74 69.89 72.49 75.09 69.34  750:63.59 55.01 46.42 56.61 66.81  775:65.09 63.38 63.84 64.30 61.88  800:59.45 55.71 51.96 54.70 57.44  825:58.88 60.31"
15
16	constant X10_spec "from 380 to 770 step 5 spectrum 0.0002 0.0007 0.0024 0.0072 0.0191 0.043 0.0847 0.1406 0.2045 0.2647 0.3147 0.3577 0.3837 0.3867 0.3707 0.3430 0.3023 0.2541 0.1956 0.1323 0.0805 0.0411 0.0162 0.0051 0.0038 0.0154 0.0375 0.0714 0.1177 0.1730 0.2365 0.3042 0.3768 0.4516 0.5298 0.6161 0.7052 0.7938 0.8787 0.9512 1.0142 1.0743 1.1185 1.1343 1.1240 1.0891 1.0305 0.9507 0.8563 0.7549 0.6475 0.5351 0.4316 0.3437 0.2683 0.2043 0.1526 0.1122 0.0813 0.0579 0.0409 0.0286 0.0199 0.0138 0.0096 0.0066 0.0046 0.0031 0.0022 0.0015 0.0010 0.0007 0.0005 0.0004 0.0003 0.0002 0.0001 0.0001 0.0001"
17	constant Y10_spec "from 385 to 755 step 5 spectrum 0.0001 0.0003 0.0008 0.0020 0.0045 0.0088 0.0145 0.0214 0.0295 0.0387 0.0496 0.0621 0.0747 0.0895 0.1063 0.1282 0.1528 0.1852 0.2199 0.2536 0.2977 0.3391 0.3954 0.4608 0.5314 0.6067 0.6857 0.7618 0.8233 0.8752 0.9238 0.9620 0.9822 0.9918 0.9991 0.9973 0.9824 0.9556 0.9152 0.8689 0.8256 0.7774 0.7204 0.6583 0.5939 0.5280 0.4618 0.3981 0.3396 0.2835 0.2283 0.1798 0.1402 0.1076 0.0812 0.0603 0.0441 0.0318 0.0226 0.0159 0.0111 0.0077 0.0054 0.0037 0.0026 0.0018 0.0012 0.0008 0.0006 0.0004 0.0003 0.0002 0.0001 0.0001 0.0001"
18	constant Z10_spec "from 380 to 555 step 5 spectrum 0.0007 0.0029 0.0105 0.0323 0.0860 0.1971 0.3894 0.6568 0.9725 1.2825 1.5535 1.7985 1.9673 2.0273 1.9948 1.9007 1.7454 1.5549 1.3176 1.0302 0.7721 0.5701 0.4153 0.3024 0.2185 0.1592 0.1120 0.0822 0.0607 0.0431 0.0305 0.0206 0.0137 0.0079 0.0040 0.0011"
19
20
21
22	doc
23	fixed [ColorSectrum ] ; [is storing a color as a sprectal curve, with arbitrary step.]
24
25
26	type ColorSpectrumMeasure
27	field Float measure
28
29	function build  m
30	arg_w ColorSpectrumMeasure m
31	m measure := 0
32
33	type ColorSpectrum
34	field Array:ColorSpectrumMeasure spectrum
35	field Float step <- undefined
36
37	function build  cs
38	arg_w ColorSpectrum cs
39	cs step := undefined
40
41
42	function 'cast Status' cs -> status
43	arg ColorSpectrum cs ; arg Status status
44	extension
45	status := cast cs:step Status
46
47
48	method cs set_step nm
49	arg_rw ColorSpectrum cs ; arg Float nm
50	cs:spectrum size := 0
51	cs step := nm
52
53	method cs set_measure nm m
54	arg_rw ColorSpectrum cs ; arg Float nm ; arg Float m
55	check cs:step=defined
56	var Int i := cast nm/cs:step Int
57	check i>0
58	check (abs i-nm/cs:step)<1e-6
59	if i>=cs:spectrum:size
60	cs:spectrum size := i+1
61	cs:spectrum:i measure := m
62
63
64	method cs measure i -> m
65	arg ColorSpectrum cs ; arg Int i ; arg Float m
66	if i>=0 and i<cs:spectrum:size
67	m := cs:spectrum:i:measure
68	else
69	m := 0
70
71	method cs '' nm -> intensity
72	arg ColorSpectrum cs ; arg Float nm intensity
73	check cs:step=defined
74	var Int i := cast nm/cs:step-0.5 Int ; var Float r := nm/cs:step-i
75	return (1.0-r)*(cs measure i)+r*(cs measure i+1)
76
77
78	export ColorSpectrum '' 'cast Status'
79	export '. set_step' '. set_measure'
80	export '. step' ''
81
82
83
84
85
86	method cs 'to string' options -> s
87	arg ColorSpectrum cs ; arg Str options s
88	if cs=undefined
89	return "?"
90	var Int i := 0
91	while i<cs:spectrum:size and cs:spectrum:i:measure=0
92	i += 1
93	s := ""
94	s += "from "+(string i*cs:step)
95	s += " to "+(string (cs:spectrum:size-1)*cs:step)
96	s += " step "+(string cs:step)
97	s += " spectrum"
98	while i<cs:spectrum:size
99	s += " "+(string cs:spectrum:i:measure)
100	i += 1
101
102	method cs 'from string' string options may_skip skiped offset -> status
103	arg_w ColorSpectrum cs ; arg Str string options ; arg CBool may_skip ; arg_w Int skiped offset ; arg Status status
104	var Float mini := string option "from" Float
105	var Float maxi := string option "to" Float
106	var Float step := string option "step" Float
107	var Int p := string option_position "spectrum" undefined
108	if mini=undefined or maxi=undefined or step=undefined or p=undefined
109	return failure
110	cs set_step step
111	(string p string:len) parse "spectrum" any:(var Str measures)
112	var Float nm := mini
113	while (measures parse (var Int check_nm) ":" (var Float m1) any:(var Str remain)) or { check_nm := undefined ; measures parse (var Float m1) any:(var Str remain) }
114	if check_nm=defined and nm<>check_nm
115	return failure
116	cs set_measure nm m1
117	nm += step
118	measures := remain
119	if (abs nm-step-maxi)>1e-6
120	return failure
121	skiped := 0 ; offset := string len ; status := success
122
123
124
125
126	doc
127	fixed [ColorSectrum32 ] ; [is a more efficient implementation, using 10 nm step, ranging from 400 to 700.] ; eol
128	[All computations (+ - * / ^ integral and modulus) are defined for] ; fixed [ ColorSectrum32 ] ; [data type.]
129
130
131	constant nm_min 400
132	constant nm_max 700
133	constant nm_step 10
134	constant nm_count (nm_max-nm_min)\nm_step+1
135
136
137	type ColorSpectrum32
138	field (Array Float32 nm_count) measure
139
140	function build  s32
141	arg_w ColorSpectrum32 s32
142	s32:measure 0 := undefined
143
144	function 'cast Status' s32 -> status
145	arg ColorSpectrum32 s32 ; arg Status status
146	status := cast s32:measure:0 Status
147
148
149	function 'cast ColorSpectrum32' c -> s32
150	arg Float c ; arg ColorSpectrum32 s32
151	explicit
152	for (var Int i) 0 nm_count-1
153	s32:measure i := c
154
155	function 'cast ColorSpectrum32' s -> s32
156	arg ColorSpectrum s ; arg ColorSpectrum32 s32
157	extension
158	if s=failure
159	s32:measure 0 := undefined
160	return
161	for (var Int i) 0 nm_count-1
162	s32:measure i := s nm_min+i*nm_step
163
164	function 'cast ColorSpectrum' s32 -> s
165	arg ColorSpectrum32 s32 ; arg ColorSpectrum s
166	explicit
167	if s32=failure
168	s := var ColorSpectrum no_spectrum
169	return
170	s set_step nm_step
171	for (var Int i) 0 nm_count-1
172	s set_measure nm_min+i*nm_step s32:measure:i
173
174	function color_spectrum32 ascii -> s32
175	arg Str ascii ; arg ColorSpectrum32 s32
176	if (ascii parse (var ColorSpectrum s))
177	s32 := s
178	else
179	s32:measure 0 := undefined
180
181
182	function '+' s1 s2 -> s
183	arg ColorSpectrum32 s1 s2 s
184	for (var Int i) 0 nm_count-1
185	s:measure i := s1:measure:i + s2:measure:i
186
187	function '-' s1 s2 -> s
188	arg ColorSpectrum32 s1 s2 s
189	for (var Int i) 0 nm_count-1
190	s:measure i := s1:measure:i - s2:measure:i
191
192	function '*' c s1 -> s
193	arg Float c ; arg ColorSpectrum32 s1 s
194	for (var Int i) 0 nm_count-1
195	s:measure i := c * s1:measure:i
196
197	function '*' s1 s2 -> s
198	arg ColorSpectrum32 s1 s2 s
199	for (var Int i) 0 nm_count-1
200	s:measure i := s1:measure:i * s2:measure:i
201
202	function '/' s1 s2 -> s
203	arg ColorSpectrum32 s1 s2 s
204	for (var Int i) 0 nm_count-1
205	var Float f := s1:measure:i / s2:measure:i
206	if f=undefined
207	f := 1
208	s:measure i := f
209
210	function '^' s1 c -> s
211	arg ColorSpectrum32 s1 ; arg Float c ; arg ColorSpectrum32 s
212	for (var Int i) 0 nm_count-1
213	var Float f := s1:measure:i ^ c
214	if f=undefined
215	# console s1:measure:i " ^ " c " -> " f eol
216	if c=0
217	f := 1
218	else
219	f := s1:measure:i
220	s:measure i := f
221
222	method s integral -> f
223	arg ColorSpectrum32 s ; arg Float f
224	f := 0
225	for (var Int i) 0 nm_count-1
226	f += s:measure i
227	f *= nm_step
228
229	method s modulus -> f
230	arg ColorSpectrum32 s ; arg Float f
231	f := 0
232	for (var Int i) 0 nm_count-1
233	f += abs s:measure:i
234	f *= nm_step
235
236	function min s1 s2 -> s
237	arg ColorSpectrum32 s1 s2 s
238	for (var Int i) 0 nm_count-1
239	s:measure i := min s1:measure:i s2:measure:i
240
241	function max s1 s2 -> s
242	arg ColorSpectrum32 s1 s2 s
243	for (var Int i) 0 nm_count-1
244	s:measure i := max s1:measure:i s2:measure:i
245
246	export ColorSpectrum32 'cast Status' 'cast ColorSpectrum32' 'cast ColorSpectrum' color_spectrum32 '+' '-' '*' '/' '^' '. integral' '. modulus' min max
247
248
249	function log s1 -> s
250	arg ColorSpectrum32 s1 s
251	for (var Int i) 0 nm_count-1
252	s:measure i := log s1:measure:i
253
254	function exp s1 -> s
255	arg ColorSpectrum32 s1 s
256	for (var Int i) 0 nm_count-1
257	s:measure i := exp s1:measure:i
258
259	export log exp
260
261
262
263
264	doc
265	[Correction formulae]
266
267
268	function exposure x gg -> y
269	arg Float x gg y
270	var Float g := -8*gg
271	y := ((exp x*g)-1)/(exp:g-1)
272	if y=undefined
273	y := x
274
275	function unexposure y gg -> x
276	arg Float y gg x
277	var Float g := -8*gg
278	# console "exp:g-1 = " exp:g-1 eol
279	# console "pre log = " y*(exp:g-1)+1 eol
280	x := (log y*(exp:g-1)+1)/g
281	if x=undefined or abs:gg<0.01 and (abs x-y)>abs:gg
282	x := y
283
284	function exposure s1 e -> s
285	arg ColorSpectrum32 s1 ; arg Float e ; arg ColorSpectrum32 s
286	for (var Int i) 0 nm_count-1
287	s:measure i := exposure s1:measure:i e
288
289	function unexposure s1 u -> s
290	arg ColorSpectrum32 s1 ; arg Float u ; arg ColorSpectrum32 s
291	for (var Int i) 0 nm_count-1
292	s:measure i := unexposure s1:measure:i u
293
294	export exposure unexposure
295
296
297
298
299	doc
300	[Now we define the CIE X, Y and Z value associated with a color spectrum, as defined by the CIE]
301
302
303	public
304	gvar ColorSpectrum32 X_spectrum Y_spectrum Z_spectrum illuminant_spectrum
305
306	X_spectrum := color_spectrum32 X_spec
307	Y_spectrum := color_spectrum32 Y_spec
308	Z_spectrum := color_spectrum32 Z_spec
309	# X_spectrum := color_spectrum32 X10_spec
310	# Y_spectrum := color_spectrum32 Y10_spec
311	# Z_spectrum := color_spectrum32 Z10_spec
312	illuminant_spectrum := 1/(color_spectrum32:D50_spec*Y_spectrum integral) * color_spectrum32:D50_spec
313	# illuminant_spectrum := 1/(color_spectrum32:D55_spec*Y_spectrum integral) * color_spectrum32:D55_spec
314	# illuminant_spectrum := 1/(color_spectrum32:D65_spec*Y_spectrum integral) * color_spectrum32:D65_spec
315
316
317	method s X -> f
318	arg ColorSpectrum32 s ; arg Float f
319	f := 0
320	for (var Int i) 0 nm_count-1
321	f += s:measure:i * X_spectrum:measure:i
322	f *= nm_step
323
324	method s Y -> f
325	arg ColorSpectrum32 s ; arg Float f
326	f := 0
327	for (var Int i) 0 nm_count-1
328	f += s:measure:i * Y_spectrum:measure:i
329	f *= nm_step
330
331	method s Z -> f
332	arg ColorSpectrum32 s ; arg Float f
333	f := 0
334	for (var Int i) 0 nm_count-1
335	f += s:measure:i * Z_spectrum:measure:i
336	f *= nm_step
337
338
339	export '. X' '. Y' '. Z'
340
341	if false
342	console "D50 illuminant" eol
343	console "x = " illuminant_spectrum:X/(illuminant_spectrum:X+illuminant_spectrum:Y+illuminant_spectrum:Z) eol
344	console "y = " illuminant_spectrum:Y/(illuminant_spectrum:X+illuminant_spectrum:Y+illuminant_spectrum:Z) eol
345	console "Y = " illuminant_spectrum:Y eol
346