gfx/qcms/transform_util.c@97036ab72558
gfx/qcms/transform_util.c
Tue, 06 Jan 2015 21:39:09 +0100
- author
- Michael Schloh von Bennewitz <michael@schloh.com>
- date
- Tue, 06 Jan 2015 21:39:09 +0100
- branch
- TOR_BUG_9701
- changeset 8
- 97036ab72558
- permissions
- -rw-r--r--
Conditionally force memory storage according to privacy.thirdparty.isolate;
This solves Tor bug #9701, complying with disk avoidance documented in
https://www.torproject.org/projects/torbrowser/design/#disk-avoidance.
1 #define _ISOC99_SOURCE /* for INFINITY */
3 #include <math.h>
4 #include <assert.h>
5 #include <string.h> //memcpy
6 #include "qcmsint.h"
7 #include "transform_util.h"
8 #include "matrix.h"
10 #if !defined(INFINITY)
11 #define INFINITY HUGE_VAL
12 #endif
14 #define PARAMETRIC_CURVE_TYPE 0x70617261 //'para'
16 /* value must be a value between 0 and 1 */
17 //XXX: is the above a good restriction to have?
18 // the output range of this functions is 0..1
19 float lut_interp_linear(double input_value, uint16_t *table, int length)
20 {
21 int upper, lower;
22 float value;
23 input_value = input_value * (length - 1); // scale to length of the array
24 upper = ceil(input_value);
25 lower = floor(input_value);
26 //XXX: can we be more performant here?
27 value = table[upper]*(1. - (upper - input_value)) + table[lower]*(upper - input_value);
28 /* scale the value */
29 return value * (1.f/65535.f);
30 }
32 /* same as above but takes and returns a uint16_t value representing a range from 0..1 */
33 uint16_t lut_interp_linear16(uint16_t input_value, uint16_t *table, int length)
34 {
35 /* Start scaling input_value to the length of the array: 65535*(length-1).
36 * We'll divide out the 65535 next */
37 uint32_t value = (input_value * (length - 1));
38 uint32_t upper = (value + 65534) / 65535; /* equivalent to ceil(value/65535) */
39 uint32_t lower = value / 65535; /* equivalent to floor(value/65535) */
40 /* interp is the distance from upper to value scaled to 0..65535 */
41 uint32_t interp = value % 65535;
43 value = (table[upper]*(interp) + table[lower]*(65535 - interp))/65535; // 0..65535*65535
45 return value;
46 }
48 /* same as above but takes an input_value from 0..PRECACHE_OUTPUT_MAX
49 * and returns a uint8_t value representing a range from 0..1 */
50 static
51 uint8_t lut_interp_linear_precache_output(uint32_t input_value, uint16_t *table, int length)
52 {
53 /* Start scaling input_value to the length of the array: PRECACHE_OUTPUT_MAX*(length-1).
54 * We'll divide out the PRECACHE_OUTPUT_MAX next */
55 uint32_t value = (input_value * (length - 1));
57 /* equivalent to ceil(value/PRECACHE_OUTPUT_MAX) */
58 uint32_t upper = (value + PRECACHE_OUTPUT_MAX-1) / PRECACHE_OUTPUT_MAX;
59 /* equivalent to floor(value/PRECACHE_OUTPUT_MAX) */
60 uint32_t lower = value / PRECACHE_OUTPUT_MAX;
61 /* interp is the distance from upper to value scaled to 0..PRECACHE_OUTPUT_MAX */
62 uint32_t interp = value % PRECACHE_OUTPUT_MAX;
64 /* the table values range from 0..65535 */
65 value = (table[upper]*(interp) + table[lower]*(PRECACHE_OUTPUT_MAX - interp)); // 0..(65535*PRECACHE_OUTPUT_MAX)
67 /* round and scale */
68 value += (PRECACHE_OUTPUT_MAX*65535/255)/2;
69 value /= (PRECACHE_OUTPUT_MAX*65535/255); // scale to 0..255
70 return value;
71 }
73 /* value must be a value between 0 and 1 */
74 //XXX: is the above a good restriction to have?
75 float lut_interp_linear_float(float value, float *table, int length)
76 {
77 int upper, lower;
78 value = value * (length - 1);
79 upper = ceilf(value);
80 lower = floorf(value);
81 //XXX: can we be more performant here?
82 value = table[upper]*(1. - (upper - value)) + table[lower]*(upper - value);
83 /* scale the value */
84 return value;
85 }
87 #if 0
88 /* if we use a different representation i.e. one that goes from 0 to 0x1000 we can be more efficient
89 * because we can avoid the divisions and use a shifting instead */
90 /* same as above but takes and returns a uint16_t value representing a range from 0..1 */
91 uint16_t lut_interp_linear16(uint16_t input_value, uint16_t *table, int length)
92 {
93 uint32_t value = (input_value * (length - 1));
94 uint32_t upper = (value + 4095) / 4096; /* equivalent to ceil(value/4096) */
95 uint32_t lower = value / 4096; /* equivalent to floor(value/4096) */
96 uint32_t interp = value % 4096;
98 value = (table[upper]*(interp) + table[lower]*(4096 - interp))/4096; // 0..4096*4096
100 return value;
101 }
102 #endif
104 void compute_curve_gamma_table_type1(float gamma_table[256], uint16_t gamma)
105 {
106 unsigned int i;
107 float gamma_float = u8Fixed8Number_to_float(gamma);
108 for (i = 0; i < 256; i++) {
109 // 0..1^(0..255 + 255/256) will always be between 0 and 1
110 gamma_table[i] = pow(i/255., gamma_float);
111 }
112 }
114 void compute_curve_gamma_table_type2(float gamma_table[256], uint16_t *table, int length)
115 {
116 unsigned int i;
117 for (i = 0; i < 256; i++) {
118 gamma_table[i] = lut_interp_linear(i/255., table, length);
119 }
120 }
122 void compute_curve_gamma_table_type_parametric(float gamma_table[256], float parameter[7], int count)
123 {
124 size_t X;
125 float interval;
126 float a, b, c, e, f;
127 float y = parameter[0];
128 if (count == 0) {
129 a = 1;
130 b = 0;
131 c = 0;
132 e = 0;
133 f = 0;
134 interval = -INFINITY;
135 } else if(count == 1) {
136 a = parameter[1];
137 b = parameter[2];
138 c = 0;
139 e = 0;
140 f = 0;
141 interval = -1 * parameter[2] / parameter[1];
142 } else if(count == 2) {
143 a = parameter[1];
144 b = parameter[2];
145 c = 0;
146 e = parameter[3];
147 f = parameter[3];
148 interval = -1 * parameter[2] / parameter[1];
149 } else if(count == 3) {
150 a = parameter[1];
151 b = parameter[2];
152 c = parameter[3];
153 e = -c;
154 f = 0;
155 interval = parameter[4];
156 } else if(count == 4) {
157 a = parameter[1];
158 b = parameter[2];
159 c = parameter[3];
160 e = parameter[5] - c;
161 f = parameter[6];
162 interval = parameter[4];
163 } else {
164 assert(0 && "invalid parametric function type.");
165 a = 1;
166 b = 0;
167 c = 0;
168 e = 0;
169 f = 0;
170 interval = -INFINITY;
171 }
172 for (X = 0; X < 256; X++) {
173 if (X >= interval) {
174 // XXX The equations are not exactly as definied in the spec but are
175 // algebraic equivilent.
176 // TODO Should division by 255 be for the whole expression.
177 gamma_table[X] = clamp_float(pow(a * X / 255. + b, y) + c + e);
178 } else {
179 gamma_table[X] = clamp_float(c * X / 255. + f);
180 }
181 }
182 }
184 void compute_curve_gamma_table_type0(float gamma_table[256])
185 {
186 unsigned int i;
187 for (i = 0; i < 256; i++) {
188 gamma_table[i] = i/255.;
189 }
190 }
192 float *build_input_gamma_table(struct curveType *TRC)
193 {
194 float *gamma_table;
196 if (!TRC) return NULL;
197 gamma_table = malloc(sizeof(float)*256);
198 if (gamma_table) {
199 if (TRC->type == PARAMETRIC_CURVE_TYPE) {
200 compute_curve_gamma_table_type_parametric(gamma_table, TRC->parameter, TRC->count);
201 } else {
202 if (TRC->count == 0) {
203 compute_curve_gamma_table_type0(gamma_table);
204 } else if (TRC->count == 1) {
205 compute_curve_gamma_table_type1(gamma_table, TRC->data[0]);
206 } else {
207 compute_curve_gamma_table_type2(gamma_table, TRC->data, TRC->count);
208 }
209 }
210 }
211 return gamma_table;
212 }
214 struct matrix build_colorant_matrix(qcms_profile *p)
215 {
216 struct matrix result;
217 result.m[0][0] = s15Fixed16Number_to_float(p->redColorant.X);
218 result.m[0][1] = s15Fixed16Number_to_float(p->greenColorant.X);
219 result.m[0][2] = s15Fixed16Number_to_float(p->blueColorant.X);
220 result.m[1][0] = s15Fixed16Number_to_float(p->redColorant.Y);
221 result.m[1][1] = s15Fixed16Number_to_float(p->greenColorant.Y);
222 result.m[1][2] = s15Fixed16Number_to_float(p->blueColorant.Y);
223 result.m[2][0] = s15Fixed16Number_to_float(p->redColorant.Z);
224 result.m[2][1] = s15Fixed16Number_to_float(p->greenColorant.Z);
225 result.m[2][2] = s15Fixed16Number_to_float(p->blueColorant.Z);
226 result.invalid = false;
227 return result;
228 }
230 /* The following code is copied nearly directly from lcms.
231 * I think it could be much better. For example, Argyll seems to have better code in
232 * icmTable_lookup_bwd and icmTable_setup_bwd. However, for now this is a quick way
233 * to a working solution and allows for easy comparing with lcms. */
234 uint16_fract_t lut_inverse_interp16(uint16_t Value, uint16_t LutTable[], int length)
235 {
236 int l = 1;
237 int r = 0x10000;
238 int x = 0, res; // 'int' Give spacing for negative values
239 int NumZeroes, NumPoles;
240 int cell0, cell1;
241 double val2;
242 double y0, y1, x0, x1;
243 double a, b, f;
245 // July/27 2001 - Expanded to handle degenerated curves with an arbitrary
246 // number of elements containing 0 at the begining of the table (Zeroes)
247 // and another arbitrary number of poles (FFFFh) at the end.
248 // First the zero and pole extents are computed, then value is compared.
250 NumZeroes = 0;
251 while (LutTable[NumZeroes] == 0 && NumZeroes < length-1)
252 NumZeroes++;
254 // There are no zeros at the beginning and we are trying to find a zero, so
255 // return anything. It seems zero would be the less destructive choice
256 /* I'm not sure that this makes sense, but oh well... */
257 if (NumZeroes == 0 && Value == 0)
258 return 0;
260 NumPoles = 0;
261 while (LutTable[length-1- NumPoles] == 0xFFFF && NumPoles < length-1)
262 NumPoles++;
264 // Does the curve belong to this case?
265 if (NumZeroes > 1 || NumPoles > 1)
266 {
267 int a, b;
269 // Identify if value fall downto 0 or FFFF zone
270 if (Value == 0) return 0;
271 // if (Value == 0xFFFF) return 0xFFFF;
273 // else restrict to valid zone
275 a = ((NumZeroes-1) * 0xFFFF) / (length-1);
276 b = ((length-1 - NumPoles) * 0xFFFF) / (length-1);
278 l = a - 1;
279 r = b + 1;
280 }
283 // Seems not a degenerated case... apply binary search
285 while (r > l) {
287 x = (l + r) / 2;
289 res = (int) lut_interp_linear16((uint16_fract_t) (x-1), LutTable, length);
291 if (res == Value) {
293 // Found exact match.
295 return (uint16_fract_t) (x - 1);
296 }
298 if (res > Value) r = x - 1;
299 else l = x + 1;
300 }
302 // Not found, should we interpolate?
305 // Get surrounding nodes
307 val2 = (length-1) * ((double) (x - 1) / 65535.0);
309 cell0 = (int) floor(val2);
310 cell1 = (int) ceil(val2);
312 if (cell0 == cell1) return (uint16_fract_t) x;
314 y0 = LutTable[cell0] ;
315 x0 = (65535.0 * cell0) / (length-1);
317 y1 = LutTable[cell1] ;
318 x1 = (65535.0 * cell1) / (length-1);
320 a = (y1 - y0) / (x1 - x0);
321 b = y0 - a * x0;
323 if (fabs(a) < 0.01) return (uint16_fract_t) x;
325 f = ((Value - b) / a);
327 if (f < 0.0) return (uint16_fract_t) 0;
328 if (f >= 65535.0) return (uint16_fract_t) 0xFFFF;
330 return (uint16_fract_t) floor(f + 0.5);
332 }
334 /*
335 The number of entries needed to invert a lookup table should not
336 necessarily be the same as the original number of entries. This is
337 especially true of lookup tables that have a small number of entries.
339 For example:
340 Using a table like:
341 {0, 3104, 14263, 34802, 65535}
342 invert_lut will produce an inverse of:
343 {3, 34459, 47529, 56801, 65535}
344 which has an maximum error of about 9855 (pixel difference of ~38.346)
346 For now, we punt the decision of output size to the caller. */
347 static uint16_t *invert_lut(uint16_t *table, int length, int out_length)
348 {
349 int i;
350 /* for now we invert the lut by creating a lut of size out_length
351 * and attempting to lookup a value for each entry using lut_inverse_interp16 */
352 uint16_t *output = malloc(sizeof(uint16_t)*out_length);
353 if (!output)
354 return NULL;
356 for (i = 0; i < out_length; i++) {
357 double x = ((double) i * 65535.) / (double) (out_length - 1);
358 uint16_fract_t input = floor(x + .5);
359 output[i] = lut_inverse_interp16(input, table, length);
360 }
361 return output;
362 }
364 static void compute_precache_pow(uint8_t *output, float gamma)
365 {
366 uint32_t v = 0;
367 for (v = 0; v < PRECACHE_OUTPUT_SIZE; v++) {
368 //XXX: don't do integer/float conversion... and round?
369 output[v] = 255. * pow(v/(double)PRECACHE_OUTPUT_MAX, gamma);
370 }
371 }
373 void compute_precache_lut(uint8_t *output, uint16_t *table, int length)
374 {
375 uint32_t v = 0;
376 for (v = 0; v < PRECACHE_OUTPUT_SIZE; v++) {
377 output[v] = lut_interp_linear_precache_output(v, table, length);
378 }
379 }
381 void compute_precache_linear(uint8_t *output)
382 {
383 uint32_t v = 0;
384 for (v = 0; v < PRECACHE_OUTPUT_SIZE; v++) {
385 //XXX: round?
386 output[v] = v / (PRECACHE_OUTPUT_SIZE/256);
387 }
388 }
390 qcms_bool compute_precache(struct curveType *trc, uint8_t *output)
391 {
393 if (trc->type == PARAMETRIC_CURVE_TYPE) {
394 float gamma_table[256];
395 uint16_t gamma_table_uint[256];
396 uint16_t i;
397 uint16_t *inverted;
398 int inverted_size = 256;
400 compute_curve_gamma_table_type_parametric(gamma_table, trc->parameter, trc->count);
401 for(i = 0; i < 256; i++) {
402 gamma_table_uint[i] = (uint16_t)(gamma_table[i] * 65535);
403 }
405 //XXX: the choice of a minimum of 256 here is not backed by any theory,
406 // measurement or data, howeve r it is what lcms uses.
407 // the maximum number we would need is 65535 because that's the
408 // accuracy used for computing the pre cache table
409 if (inverted_size < 256)
410 inverted_size = 256;
412 inverted = invert_lut(gamma_table_uint, 256, inverted_size);
413 if (!inverted)
414 return false;
415 compute_precache_lut(output, inverted, inverted_size);
416 free(inverted);
417 } else {
418 if (trc->count == 0) {
419 compute_precache_linear(output);
420 } else if (trc->count == 1) {
421 compute_precache_pow(output, 1./u8Fixed8Number_to_float(trc->data[0]));
422 } else {
423 uint16_t *inverted;
424 int inverted_size = trc->count;
425 //XXX: the choice of a minimum of 256 here is not backed by any theory,
426 // measurement or data, howeve r it is what lcms uses.
427 // the maximum number we would need is 65535 because that's the
428 // accuracy used for computing the pre cache table
429 if (inverted_size < 256)
430 inverted_size = 256;
432 inverted = invert_lut(trc->data, trc->count, inverted_size);
433 if (!inverted)
434 return false;
435 compute_precache_lut(output, inverted, inverted_size);
436 free(inverted);
437 }
438 }
439 return true;
440 }
443 static uint16_t *build_linear_table(int length)
444 {
445 int i;
446 uint16_t *output = malloc(sizeof(uint16_t)*length);
447 if (!output)
448 return NULL;
450 for (i = 0; i < length; i++) {
451 double x = ((double) i * 65535.) / (double) (length - 1);
452 uint16_fract_t input = floor(x + .5);
453 output[i] = input;
454 }
455 return output;
456 }
458 static uint16_t *build_pow_table(float gamma, int length)
459 {
460 int i;
461 uint16_t *output = malloc(sizeof(uint16_t)*length);
462 if (!output)
463 return NULL;
465 for (i = 0; i < length; i++) {
466 uint16_fract_t result;
467 double x = ((double) i) / (double) (length - 1);
468 x = pow(x, gamma); //XXX turn this conversion into a function
469 result = floor(x*65535. + .5);
470 output[i] = result;
471 }
472 return output;
473 }
475 void build_output_lut(struct curveType *trc,
476 uint16_t **output_gamma_lut, size_t *output_gamma_lut_length)
477 {
478 if (trc->type == PARAMETRIC_CURVE_TYPE) {
479 float gamma_table[256];
480 uint16_t i;
481 uint16_t *output = malloc(sizeof(uint16_t)*256);
483 if (!output) {
484 *output_gamma_lut = NULL;
485 return;
486 }
488 compute_curve_gamma_table_type_parametric(gamma_table, trc->parameter, trc->count);
489 *output_gamma_lut_length = 256;
490 for(i = 0; i < 256; i++) {
491 output[i] = (uint16_t)(gamma_table[i] * 65535);
492 }
493 *output_gamma_lut = output;
494 } else {
495 if (trc->count == 0) {
496 *output_gamma_lut = build_linear_table(4096);
497 *output_gamma_lut_length = 4096;
498 } else if (trc->count == 1) {
499 float gamma = 1./u8Fixed8Number_to_float(trc->data[0]);
500 *output_gamma_lut = build_pow_table(gamma, 4096);
501 *output_gamma_lut_length = 4096;
502 } else {
503 //XXX: the choice of a minimum of 256 here is not backed by any theory,
504 // measurement or data, however it is what lcms uses.
505 *output_gamma_lut_length = trc->count;
506 if (*output_gamma_lut_length < 256)
507 *output_gamma_lut_length = 256;
509 *output_gamma_lut = invert_lut(trc->data, trc->count, *output_gamma_lut_length);
510 }
511 }
513 }
