1 /* makesRGB.c -- build sRGB-to-linear and linear-to-sRGB conversion tables
2 *
3 * COPYRIGHT: Written by John Cunningham Bowler, 2013.
4 * To the extent possible under law, the author has waived all copyright and
5 * related or neighboring rights to this work. This work is published from:
6 * United States.
7 *
8 * Make a table to convert 8-bit sRGB encoding values into the closest 16-bit
9 * linear value.
10 *
11 * Make two tables to take a linear value scaled to 255*65535 and return an
12 * approximation to the 8-bit sRGB encoded value. Calculate the error in these
13 * tables and display it.
14 */
15
16 #define _C99_SOURCE 1
17 #include <stdio.h>
18 #include <math.h>
19 #include <stdlib.h>
20
21 /* pngpriv.h includes the definition of 'PNG_sRGB_FROM_LINEAR' which is required
22 * to verify the actual code.
23 */
24 #include "../../pngpriv.h"
25
26 #include "sRGB.h"
27
28 /* The tables are declared 'const' in pngpriv.h, so this redefines the tables to
29 * be used.
30 */
31 #define png_sRGB_table sRGB_table
32 #define png_sRGB_base sRGB_base
33 #define png_sRGB_delta sRGB_delta
34
35 static png_uint_16 png_sRGB_table[256];
36 static png_uint_16 png_sRGB_base[512];
37 static png_byte png_sRGB_delta[512];
38
39 static const unsigned int max_input = 255*65535;
40
41 double
42 fsRGB(double l)
43 {
44 return sRGB_from_linear(l/max_input);
45 }
46
47 double
48 sRGB(unsigned int i)
49 {
50 return fsRGB(i);
51 }
52
53 double
54 finvsRGB(unsigned int i)
55 {
56 return 65535 * linear_from_sRGB(i/255.);
57 }
58
59 png_uint_16
60 invsRGB(unsigned int i)
61 {
62 unsigned int x = nearbyint(finvsRGB(i));
63
64 if (x > 65535)
65 {
66 fprintf(stderr, "invsRGB(%u) overflows to %u\n", i, x);
67 exit(1);
68 }
69
70 return (png_uint_16)x;
71 }
72
73 int
74 main(int argc, char **argv)
75 {
76 unsigned int i, i16, ibase;
77 double min_error = 0;
78 double max_error = 0;
79 double min_error16 = 0;
80 double max_error16 = 0;
81 double adjust;
82 double adjust_lo = 0.4, adjust_hi = 0.6, adjust_mid = 0.5;
83 unsigned int ec_lo = 0, ec_hi = 0, ec_mid = 0;
84 unsigned int error_count = 0;
85 unsigned int error_count16 = 0;
86 int test_only = 0;
87
88 if (argc > 1)
89 test_only = strcmp("--test", argv[1]) == 0;
90
91 /* Initialize the encoding table first. */
92 for (i=0; i<256; ++i)
93 {
94 png_sRGB_table[i] = invsRGB(i);
95 }
96
97 /* Now work out the decoding tables (this is where the error comes in because
98 * there are 512 set points and 512 straight lines between them.)
99 */
100 for (;;)
101 {
102 if (ec_lo == 0)
103 adjust = adjust_lo;
104
105 else if (ec_hi == 0)
106 adjust = adjust_hi;
107
108 else if (ec_mid == 0)
109 adjust = adjust_mid;
110
111 else if (ec_mid < ec_hi)
112 adjust = (adjust_mid + adjust_hi)/2;
113
114 else if (ec_mid < ec_lo)
115 adjust = (adjust_mid + adjust_lo)/2;
116
117 else
118 {
119 fprintf(stderr, "not reached: %u .. %u .. %u\n", ec_lo, ec_mid, ec_hi);
120 exit(1);
121 }
122
123 /* Calculate the table using the current 'adjust' */
124 for (i=0; i<=511; ++i)
125 {
126 double lo = 255 * sRGB(i << 15);
127 double hi = 255 * sRGB((i+1) << 15);
128 unsigned int calc;
129
130 calc = nearbyint((lo+adjust) * 256);
131 if (calc > 65535)
132 {
133 fprintf(stderr, "table[%d][0]: overflow %08x (%d)\n", i, calc,
134 calc);
135 exit(1);
136 }
137 png_sRGB_base[i] = calc;
138
139 calc = nearbyint((hi-lo) * 32);
140 if (calc > 255)
141 {
142 fprintf(stderr, "table[%d][1]: overflow %08x (%d)\n", i, calc,
143 calc);
144 exit(1);
145 }
146 png_sRGB_delta[i] = calc;
147 }
148
149 /* Check the 16-bit linear values alone: */
150 error_count16 = 0;
151 for (i16=0; i16 <= 65535; ++i16)
152 {
153 unsigned int i = 255*i16;
154 unsigned int iexact = nearbyint(255*sRGB(i));
155 unsigned int icalc = PNG_sRGB_FROM_LINEAR(i);
156
157 if (icalc != iexact)
158 ++error_count16;
159 }
160
161 /* Now try changing the adjustment. */
162 if (ec_lo == 0)
163 ec_lo = error_count16;
164
165 else if (ec_hi == 0)
166 ec_hi = error_count16;
167
168 else if (ec_mid == 0)
169 {
170 ec_mid = error_count16;
171 printf("/* initial error counts: %u .. %u .. %u */\n", ec_lo, ec_mid,
172 ec_hi);
173 }
174
175 else if (error_count16 < ec_mid)
176 {
177 printf("/* adjust (mid ): %f: %u -> %u */\n", adjust, ec_mid,
178 error_count16);
179 ec_mid = error_count16;
180 adjust_mid = adjust;
181 }
182
183 else if (adjust < adjust_mid && error_count16 < ec_lo)
184 {
185 printf("/* adjust (low ): %f: %u -> %u */\n", adjust, ec_lo,
186 error_count16);
187 ec_lo = error_count16;
188 adjust_lo = adjust;
189 }
190
191 else if (adjust > adjust_mid && error_count16 < ec_hi)
192 {
193 printf("/* adjust (high): %f: %u -> %u */\n", adjust, ec_hi,
194 error_count16);
195 ec_hi = error_count16;
196 adjust_hi = adjust;
197 }
198
199 else
200 {
201 adjust = adjust_mid;
202 printf("/* adjust: %f: %u */\n", adjust, ec_mid);
203 break;
204 }
205 }
206
207 /* For each entry in the table try to adjust it to minimize the error count
208 * in that entry. Each entry corresponds to 128 input values.
209 */
210 for (ibase=0; ibase<65536; ibase+=128)
211 {
212 png_uint_16 base = png_sRGB_base[ibase >> 7], trybase = base, ob=base;
213 png_byte delta = png_sRGB_delta[ibase >> 7], trydelta = delta, od=delta;
214 unsigned int ecbase = 0, eco;
215
216 for (;;)
217 {
218 png_sRGB_base[ibase >> 7] = trybase;
219 png_sRGB_delta[ibase >> 7] = trydelta;
220
221 /* Check the 16-bit linear values alone: */
222 error_count16 = 0;
223 for (i16=ibase; i16 < ibase+128; ++i16)
224 {
225 unsigned int i = 255*i16;
226 unsigned int iexact = nearbyint(255*sRGB(i));
227 unsigned int icalc = PNG_sRGB_FROM_LINEAR(i);
228
229 if (icalc != iexact)
230 ++error_count16;
231 }
232
233 if (error_count16 == 0)
234 break;
235
236 if (ecbase == 0)
237 {
238 eco = ecbase = error_count16;
239 ++trybase; /* First test */
240 }
241
242 else if (error_count16 < ecbase)
243 {
244 if (trybase > base)
245 {
246 base = trybase;
247 ++trybase;
248 }
249 else if (trybase < base)
250 {
251 base = trybase;
252 --trybase;
253 }
254 else if (trydelta > delta)
255 {
256 delta = trydelta;
257 ++trydelta;
258 }
259 else if (trydelta < delta)
260 {
261 delta = trydelta;
262 --trydelta;
263 }
264 else
265 {
266 fprintf(stderr, "makesRGB: impossible\n");
267 exit(1);
268 }
269 ecbase = error_count16;
270 }
271
272 else
273 {
274 if (trybase > base)
275 trybase = base-1;
276 else if (trybase < base)
277 {
278 trybase = base;
279 ++trydelta;
280 }
281 else if (trydelta > delta)
282 trydelta = delta-1;
283 else if (trydelta < delta)
284 break; /* end of tests */
285 }
286 }
287
288 png_sRGB_base[ibase >> 7] = base;
289 png_sRGB_delta[ibase >> 7] = delta;
290 if (base != ob || delta != od)
291 {
292 printf("/* table[%u]={%u,%u} -> {%u,%u} %u -> %u errors */\n",
293 ibase>>7, ob, od, base, delta, eco, ecbase);
294 }
295 else if (0)
296 printf("/* table[%u]={%u,%u} %u errors */\n", ibase>>7, ob, od,
297 ecbase);
298 }
299
300 /* Only do the full (slow) test at the end: */
301 min_error = -.4999;
302 max_error = .4999;
303 error_count = 0;
304
305 for (i=0; i <= max_input; ++i)
306 {
307 unsigned int iexact = nearbyint(255*sRGB(i));
308 unsigned int icalc = PNG_sRGB_FROM_LINEAR(i);
309
310 if (icalc != iexact)
311 {
312 double err = 255*sRGB(i) - icalc;
313
314 if (err > (max_error+.001) || err < (min_error-.001))
315 {
316 printf(
317 "/* 0x%08x: exact: %3d, got: %3d [tables: %08x, %08x] (%f) */\n",
318 i, iexact, icalc, png_sRGB_base[i>>15],
319 png_sRGB_delta[i>>15], err);
320 }
321
322 ++error_count;
323 if (err > max_error)
324 max_error = err;
325 else if (err < min_error)
326 min_error = err;
327 }
328 }
329
330 /* Re-check the 16-bit cases too, including the warning if there is an error
331 * bigger than 1.
332 */
333 error_count16 = 0;
334 max_error16 = 0;
335 min_error16 = 0;
336 for (i16=0; i16 <= 65535; ++i16)
337 {
338 unsigned int i = 255*i16;
339 unsigned int iexact = nearbyint(255*sRGB(i));
340 unsigned int icalc = PNG_sRGB_FROM_LINEAR(i);
341
342 if (icalc != iexact)
343 {
344 double err = 255*sRGB(i) - icalc;
345
346 ++error_count16;
347 if (err > max_error16)
348 max_error16 = err;
349 else if (err < min_error16)
350 min_error16 = err;
351
352 if (abs(icalc - iexact) > 1)
353 printf(
354 "/* 0x%04x: exact: %3d, got: %3d [tables: %08x, %08x] (%f) */\n",
355 i16, iexact, icalc, png_sRGB_base[i>>15],
356 png_sRGB_delta[i>>15], err);
357 }
358 }
359
360 /* Check the round trip for each 8-bit sRGB value. */
361 for (i16=0; i16 <= 255; ++i16)
362 {
363 unsigned int i = 255 * png_sRGB_table[i16];
364 unsigned int iexact = nearbyint(255*sRGB(i));
365 unsigned int icalc = PNG_sRGB_FROM_LINEAR(i);
366
367 if (i16 != iexact)
368 {
369 fprintf(stderr, "8-bit rounding error: %d -> %d\n", i16, iexact);
370 exit(1);
371 }
372
373 if (icalc != i16)
374 {
375 double finv = finvsRGB(i16);
376
377 printf("/* 8-bit roundtrip error: %d -> %f -> %d(%f) */\n",
378 i16, finv, icalc, fsRGB(255*finv));
379 }
380 }
381
382
383 printf("/* error: %g - %g, %u (%g%%) of readings inexact */\n",
384 min_error, max_error, error_count, (100.*error_count)/max_input);
385 printf("/* 16-bit error: %g - %g, %u (%g%%) of readings inexact */\n",
386 min_error16, max_error16, error_count16, (100.*error_count16)/65535);
387
388 if (!test_only)
389 {
390 printf("const png_uint_16 png_sRGB_table[256] =\n{\n ");
391 for (i=0; i<255; )
392 {
393 do
394 {
395 printf("%d,", png_sRGB_table[i++]);
396 }
397 while ((i & 0x7) != 0 && i<255);
398 if (i<255) printf("\n ");
399 }
400 printf("%d\n};\n\n", png_sRGB_table[i]);
401
402
403 printf("const png_uint_16 png_sRGB_base[512] =\n{\n ");
404 for (i=0; i<511; )
405 {
406 do
407 {
408 printf("%d,", png_sRGB_base[i++]);
409 }
410 while ((i & 0x7) != 0 && i<511);
411 if (i<511) printf("\n ");
412 }
413 printf("%d\n};\n\n", png_sRGB_base[i]);
414
415 printf("const png_byte png_sRGB_delta[512] =\n{\n ");
416 for (i=0; i<511; )
417 {
418 do
419 {
420 printf("%d,", png_sRGB_delta[i++]);
421 }
422 while ((i & 0xf) != 0 && i<511);
423 if (i<511) printf("\n ");
424 }
425 printf("%d\n};\n\n", png_sRGB_delta[i]);
426 }
427
428 return 0;
429 }