1 // qcms
2 // Copyright (C) 2009 Mozilla Foundation
3 //
4 // Permission is hereby granted, free of charge, to any person obtaining
5 // a copy of this software and associated documentation files (the "Software"),
6 // to deal in the Software without restriction, including without limitation
7 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
8 // and/or sell copies of the Software, and to permit persons to whom the Software
9 // is furnished to do so, subject to the following conditions:
10 //
11 // The above copyright notice and this permission notice shall be included in
12 // all copies or substantial portions of the Software.
13 //
14 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
15 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO
16 // THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
17 // NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
18 // LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
19 // OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
20 // WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
21
22 #define _ISOC99_SOURCE /* for INFINITY */
23
24 #include <math.h>
25 #include <assert.h>
26 #include <string.h> //memcpy
27 #include "qcmsint.h"
28 #include "transform_util.h"
29 #include "matrix.h"
30
31 #if !defined(INFINITY)
32 #define INFINITY HUGE_VAL
33 #endif
34
35 #define PARAMETRIC_CURVE_TYPE 0x70617261 //'para'
36
37 /* value must be a value between 0 and 1 */
38 //XXX: is the above a good restriction to have?
lut_interp_linear(double value,uint16_t * table,size_t length)39 float lut_interp_linear(double value, uint16_t *table, size_t length)
40 {
41 int upper, lower;
42 value = value * (length - 1); // scale to length of the array
43 upper = ceil(value);
44 lower = floor(value);
45 //XXX: can we be more performant here?
46 value = table[upper]*(1. - (upper - value)) + table[lower]*(upper - value);
47 /* scale the value */
48 return value * (1./65535.);
49 }
50
51 /* same as above but takes and returns a uint16_t value representing a range from 0..1 */
lut_interp_linear16(uint16_t input_value,uint16_t * table,size_t length)52 uint16_t lut_interp_linear16(uint16_t input_value, uint16_t *table, size_t length)
53 {
54 /* Start scaling input_value to the length of the array: 65535*(length-1).
55 * We'll divide out the 65535 next */
56 uintptr_t value = (input_value * (length - 1));
57 uint32_t upper = (value + 65534) / 65535; /* equivalent to ceil(value/65535) */
58 uint32_t lower = value / 65535; /* equivalent to floor(value/65535) */
59 /* interp is the distance from upper to value scaled to 0..65535 */
60 uint32_t interp = value % 65535;
61
62 value = (table[upper]*(interp) + table[lower]*(65535 - interp))/65535; // 0..65535*65535
63
64 return value;
65 }
66
67 /* same as above but takes an input_value from 0..PRECACHE_OUTPUT_MAX
68 * and returns a uint8_t value representing a range from 0..1 */
69 static
lut_interp_linear_precache_output(uint32_t input_value,uint16_t * table,size_t length)70 uint8_t lut_interp_linear_precache_output(uint32_t input_value, uint16_t *table, size_t length)
71 {
72 /* Start scaling input_value to the length of the array: PRECACHE_OUTPUT_MAX*(length-1).
73 * We'll divide out the PRECACHE_OUTPUT_MAX next */
74 uintptr_t value = (input_value * (length - 1));
75
76 /* equivalent to ceil(value/PRECACHE_OUTPUT_MAX) */
77 uint32_t upper = (value + PRECACHE_OUTPUT_MAX-1) / PRECACHE_OUTPUT_MAX;
78 /* equivalent to floor(value/PRECACHE_OUTPUT_MAX) */
79 uint32_t lower = value / PRECACHE_OUTPUT_MAX;
80 /* interp is the distance from upper to value scaled to 0..PRECACHE_OUTPUT_MAX */
81 uint32_t interp = value % PRECACHE_OUTPUT_MAX;
82
83 /* the table values range from 0..65535 */
84 value = (table[upper]*(interp) + table[lower]*(PRECACHE_OUTPUT_MAX - interp)); // 0..(65535*PRECACHE_OUTPUT_MAX)
85
86 /* round and scale */
87 value += (PRECACHE_OUTPUT_MAX*65535/255)/2;
88 value /= (PRECACHE_OUTPUT_MAX*65535/255); // scale to 0..255
89 return value;
90 }
91
92 /* value must be a value between 0 and 1 */
93 //XXX: is the above a good restriction to have?
lut_interp_linear_float(float value,float * table,size_t length)94 float lut_interp_linear_float(float value, float *table, size_t length)
95 {
96 int upper, lower;
97 value = value * (length - 1);
98 upper = ceil(value);
99 lower = floor(value);
100 //XXX: can we be more performant here?
101 value = table[upper]*(1. - (upper - value)) + table[lower]*(upper - value);
102 /* scale the value */
103 return value;
104 }
105
106 #if 0
107 /* if we use a different representation i.e. one that goes from 0 to 0x1000 we can be more efficient
108 * because we can avoid the divisions and use a shifting instead */
109 /* same as above but takes and returns a uint16_t value representing a range from 0..1 */
110 uint16_t lut_interp_linear16(uint16_t input_value, uint16_t *table, int length)
111 {
112 uint32_t value = (input_value * (length - 1));
113 uint32_t upper = (value + 4095) / 4096; /* equivalent to ceil(value/4096) */
114 uint32_t lower = value / 4096; /* equivalent to floor(value/4096) */
115 uint32_t interp = value % 4096;
116
117 value = (table[upper]*(interp) + table[lower]*(4096 - interp))/4096; // 0..4096*4096
118
119 return value;
120 }
121 #endif
122
compute_curve_gamma_table_type1(float gamma_table[256],double gamma)123 void compute_curve_gamma_table_type1(float gamma_table[256], double gamma)
124 {
125 unsigned int i;
126 for (i = 0; i < 256; i++) {
127 gamma_table[i] = pow(i/255., gamma);
128 }
129 }
130
compute_curve_gamma_table_type2(float gamma_table[256],uint16_t * table,int length)131 void compute_curve_gamma_table_type2(float gamma_table[256], uint16_t *table, int length)
132 {
133 unsigned int i;
134 for (i = 0; i < 256; i++) {
135 gamma_table[i] = lut_interp_linear(i/255., table, length);
136 }
137 }
138
compute_curve_gamma_table_type_parametric(float gamma_table[256],float parameter[7],int count)139 void compute_curve_gamma_table_type_parametric(float gamma_table[256], float parameter[7], int count)
140 {
141 size_t X;
142 float interval;
143 float a, b, c, e, f;
144 float y = parameter[0];
145 if (count == 0) {
146 a = 1;
147 b = 0;
148 c = 0;
149 e = 0;
150 f = 0;
151 interval = -INFINITY;
152 } else if(count == 1) {
153 a = parameter[1];
154 b = parameter[2];
155 c = 0;
156 e = 0;
157 f = 0;
158 interval = -1 * parameter[2] / parameter[1];
159 } else if(count == 2) {
160 a = parameter[1];
161 b = parameter[2];
162 c = 0;
163 e = parameter[3];
164 f = parameter[3];
165 interval = -1 * parameter[2] / parameter[1];
166 } else if(count == 3) {
167 a = parameter[1];
168 b = parameter[2];
169 c = parameter[3];
170 e = -c;
171 f = 0;
172 interval = parameter[4];
173 } else if(count == 4) {
174 a = parameter[1];
175 b = parameter[2];
176 c = parameter[3];
177 e = parameter[5] - c;
178 f = parameter[6];
179 interval = parameter[4];
180 } else {
181 assert(0 && "invalid parametric function type.");
182 a = 1;
183 b = 0;
184 c = 0;
185 e = 0;
186 f = 0;
187 interval = -INFINITY;
188 }
189 for (X = 0; X < 256; X++) {
190 if (X >= interval) {
191 // XXX The equations are not exactly as definied in the spec but are
192 // algebraic equivilent.
193 // TODO Should division by 255 be for the whole expression.
194 gamma_table[X] = pow(a * X / 255. + b, y) + c + e;
195 } else {
196 gamma_table[X] = c * X / 255. + f;
197 }
198 }
199 }
200
compute_curve_gamma_table_type0(float gamma_table[256])201 void compute_curve_gamma_table_type0(float gamma_table[256])
202 {
203 unsigned int i;
204 for (i = 0; i < 256; i++) {
205 gamma_table[i] = i/255.;
206 }
207 }
208
209
clamp_float(float a)210 float clamp_float(float a)
211 {
212 if (a > 1.)
213 return 1.;
214 else if (a < 0)
215 return 0;
216 else
217 return a;
218 }
219
clamp_u8(float v)220 unsigned char clamp_u8(float v)
221 {
222 if (v > 255.)
223 return 255;
224 else if (v < 0)
225 return 0;
226 else
227 return floor(v+.5);
228 }
229
u8Fixed8Number_to_float(uint16_t x)230 float u8Fixed8Number_to_float(uint16_t x)
231 {
232 // 0x0000 = 0.
233 // 0x0100 = 1.
234 // 0xffff = 255 + 255/256
235 return x/256.;
236 }
237
238 /* The SSE2 code uses min & max which let NaNs pass through.
239 We want to try to prevent that here by ensuring that
240 gamma table is within expected values. */
validate_gamma_table(float gamma_table[256])241 void validate_gamma_table(float gamma_table[256])
242 {
243 int i;
244 for (i = 0; i < 256; i++) {
245 // Note: we check that the gamma is not in range
246 // instead of out of range so that we catch NaNs
247 if (!(gamma_table[i] >= 0.f && gamma_table[i] <= 1.f)) {
248 gamma_table[i] = 0.f;
249 }
250 }
251 }
252
build_input_gamma_table(struct curveType * TRC)253 float *build_input_gamma_table(struct curveType *TRC)
254 {
255 float *gamma_table;
256
257 if (!TRC) return NULL;
258 gamma_table = malloc(sizeof(float)*256);
259 if (gamma_table) {
260 if (TRC->type == PARAMETRIC_CURVE_TYPE) {
261 compute_curve_gamma_table_type_parametric(gamma_table, TRC->parameter, TRC->count);
262 } else {
263 if (TRC->count == 0) {
264 compute_curve_gamma_table_type0(gamma_table);
265 } else if (TRC->count == 1) {
266 compute_curve_gamma_table_type1(gamma_table, u8Fixed8Number_to_float(TRC->data[0]));
267 } else {
268 compute_curve_gamma_table_type2(gamma_table, TRC->data, TRC->count);
269 }
270 }
271 }
272
273 validate_gamma_table(gamma_table);
274
275 return gamma_table;
276 }
277
build_colorant_matrix(qcms_profile * p)278 struct matrix build_colorant_matrix(qcms_profile *p)
279 {
280 struct matrix result;
281 result.m[0][0] = s15Fixed16Number_to_float(p->redColorant.X);
282 result.m[0][1] = s15Fixed16Number_to_float(p->greenColorant.X);
283 result.m[0][2] = s15Fixed16Number_to_float(p->blueColorant.X);
284 result.m[1][0] = s15Fixed16Number_to_float(p->redColorant.Y);
285 result.m[1][1] = s15Fixed16Number_to_float(p->greenColorant.Y);
286 result.m[1][2] = s15Fixed16Number_to_float(p->blueColorant.Y);
287 result.m[2][0] = s15Fixed16Number_to_float(p->redColorant.Z);
288 result.m[2][1] = s15Fixed16Number_to_float(p->greenColorant.Z);
289 result.m[2][2] = s15Fixed16Number_to_float(p->blueColorant.Z);
290 result.invalid = false;
291 return result;
292 }
293
294 /* The following code is copied nearly directly from lcms.
295 * I think it could be much better. For example, Argyll seems to have better code in
296 * icmTable_lookup_bwd and icmTable_setup_bwd. However, for now this is a quick way
297 * to a working solution and allows for easy comparing with lcms. */
lut_inverse_interp16(uint16_t Value,uint16_t LutTable[],int length)298 uint16_fract_t lut_inverse_interp16(uint16_t Value, uint16_t LutTable[], int length)
299 {
300 int l = 1;
301 int r = 0x10000;
302 int x = 0, res; // 'int' Give spacing for negative values
303 int NumZeroes, NumPoles;
304 int cell0, cell1;
305 double val2;
306 double y0, y1, x0, x1;
307 double a, b, f;
308
309 // July/27 2001 - Expanded to handle degenerated curves with an arbitrary
310 // number of elements containing 0 at the begining of the table (Zeroes)
311 // and another arbitrary number of poles (FFFFh) at the end.
312 // First the zero and pole extents are computed, then value is compared.
313
314 NumZeroes = 0;
315 while (LutTable[NumZeroes] == 0 && NumZeroes < length-1)
316 NumZeroes++;
317
318 // There are no zeros at the beginning and we are trying to find a zero, so
319 // return anything. It seems zero would be the less destructive choice
320 /* I'm not sure that this makes sense, but oh well... */
321 if (NumZeroes == 0 && Value == 0)
322 return 0;
323
324 NumPoles = 0;
325 while (LutTable[length-1- NumPoles] == 0xFFFF && NumPoles < length-1)
326 NumPoles++;
327
328 // Does the curve belong to this case?
329 if (NumZeroes > 1 || NumPoles > 1)
330 {
331 int a, b;
332
333 // Identify if value fall downto 0 or FFFF zone
334 if (Value == 0) return 0;
335 // if (Value == 0xFFFF) return 0xFFFF;
336
337 // else restrict to valid zone
338
339 a = ((NumZeroes-1) * 0xFFFF) / (length-1);
340 b = ((length-1 - NumPoles) * 0xFFFF) / (length-1);
341
342 l = a - 1;
343 r = b + 1;
344 }
345
346
347 // Seems not a degenerated case... apply binary search
348
349 while (r > l) {
350
351 x = (l + r) / 2;
352
353 res = (int) lut_interp_linear16((uint16_fract_t) (x-1), LutTable, length);
354
355 if (res == Value) {
356
357 // Found exact match.
358
359 return (uint16_fract_t) (x - 1);
360 }
361
362 if (res > Value) r = x - 1;
363 else l = x + 1;
364 }
365
366 // Not found, should we interpolate?
367
368
369 // Get surrounding nodes
370
371 val2 = (length-1) * ((double) (x - 1) / 65535.0);
372
373 cell0 = (int) floor(val2);
374 cell1 = (int) ceil(val2);
375
376 if (cell0 == cell1) return (uint16_fract_t) x;
377
378 y0 = LutTable[cell0] ;
379 x0 = (65535.0 * cell0) / (length-1);
380
381 y1 = LutTable[cell1] ;
382 x1 = (65535.0 * cell1) / (length-1);
383
384 a = (y1 - y0) / (x1 - x0);
385 b = y0 - a * x0;
386
387 if (fabs(a) < 0.01) return (uint16_fract_t) x;
388
389 f = ((Value - b) / a);
390
391 if (f < 0.0) return (uint16_fract_t) 0;
392 if (f >= 65535.0) return (uint16_fract_t) 0xFFFF;
393
394 return (uint16_fract_t) floor(f + 0.5);
395
396 }
397
398 /*
399 The number of entries needed to invert a lookup table should not
400 necessarily be the same as the original number of entries. This is
401 especially true of lookup tables that have a small number of entries.
402
403 For example:
404 Using a table like:
405 {0, 3104, 14263, 34802, 65535}
406 invert_lut will produce an inverse of:
407 {3, 34459, 47529, 56801, 65535}
408 which has an maximum error of about 9855 (pixel difference of ~38.346)
409
410 For now, we punt the decision of output size to the caller. */
invert_lut(uint16_t * table,int length,size_t out_length)411 static uint16_t *invert_lut(uint16_t *table, int length, size_t out_length)
412 {
413 int i;
414 /* for now we invert the lut by creating a lut of size out_length
415 * and attempting to lookup a value for each entry using lut_inverse_interp16 */
416 uint16_t *output = malloc(sizeof(uint16_t)*out_length);
417 if (!output)
418 return NULL;
419
420 for (i = 0; i < out_length; i++) {
421 double x = ((double) i * 65535.) / (double) (out_length - 1);
422 uint16_fract_t input = floor(x + .5);
423 output[i] = lut_inverse_interp16(input, table, length);
424 }
425 return output;
426 }
427
compute_precache_pow(uint8_t * output,float gamma)428 static void compute_precache_pow(uint8_t *output, float gamma)
429 {
430 uint32_t v = 0;
431 for (v = 0; v < PRECACHE_OUTPUT_SIZE; v++) {
432 //XXX: don't do integer/float conversion... and round?
433 output[v] = 255. * pow(v/(double)PRECACHE_OUTPUT_MAX, gamma);
434 }
435 }
436
compute_precache_lut(uint8_t * output,uint16_t * table,int length)437 void compute_precache_lut(uint8_t *output, uint16_t *table, int length)
438 {
439 uint32_t v = 0;
440 for (v = 0; v < PRECACHE_OUTPUT_SIZE; v++) {
441 output[v] = lut_interp_linear_precache_output(v, table, length);
442 }
443 }
444
compute_precache_linear(uint8_t * output)445 void compute_precache_linear(uint8_t *output)
446 {
447 uint32_t v = 0;
448 for (v = 0; v < PRECACHE_OUTPUT_SIZE; v++) {
449 //XXX: round?
450 output[v] = v / (PRECACHE_OUTPUT_SIZE/256);
451 }
452 }
453
compute_precache(struct curveType * trc,uint8_t * output)454 qcms_bool compute_precache(struct curveType *trc, uint8_t *output)
455 {
456
457 if (trc->type == PARAMETRIC_CURVE_TYPE) {
458 float gamma_table[256];
459 uint16_t gamma_table_uint[256];
460 uint16_t i;
461 uint16_t *inverted;
462 int inverted_size = 256;
463
464 compute_curve_gamma_table_type_parametric(gamma_table, trc->parameter, trc->count);
465 for(i = 0; i < 256; i++) {
466 gamma_table_uint[i] = (uint16_t)(gamma_table[i] * 65535);
467 }
468
469 //XXX: the choice of a minimum of 256 here is not backed by any theory,
470 // measurement or data, howeve r it is what lcms uses.
471 // the maximum number we would need is 65535 because that's the
472 // accuracy used for computing the pre cache table
473 if (inverted_size < 256)
474 inverted_size = 256;
475
476 inverted = invert_lut(gamma_table_uint, 256, inverted_size);
477 if (!inverted)
478 return false;
479 compute_precache_lut(output, inverted, inverted_size);
480 free(inverted);
481 } else {
482 if (trc->count == 0) {
483 compute_precache_linear(output);
484 } else if (trc->count == 1) {
485 compute_precache_pow(output, 1./u8Fixed8Number_to_float(trc->data[0]));
486 } else {
487 uint16_t *inverted;
488 int inverted_size = trc->count;
489 //XXX: the choice of a minimum of 256 here is not backed by any theory,
490 // measurement or data, howeve r it is what lcms uses.
491 // the maximum number we would need is 65535 because that's the
492 // accuracy used for computing the pre cache table
493 if (inverted_size < 256)
494 inverted_size = 256;
495
496 inverted = invert_lut(trc->data, trc->count, inverted_size);
497 if (!inverted)
498 return false;
499 compute_precache_lut(output, inverted, inverted_size);
500 free(inverted);
501 }
502 }
503 return true;
504 }
505
506
build_linear_table(int length)507 static uint16_t *build_linear_table(int length)
508 {
509 int i;
510 uint16_t *output = malloc(sizeof(uint16_t)*length);
511 if (!output)
512 return NULL;
513
514 for (i = 0; i < length; i++) {
515 double x = ((double) i * 65535.) / (double) (length - 1);
516 uint16_fract_t input = floor(x + .5);
517 output[i] = input;
518 }
519 return output;
520 }
521
build_pow_table(float gamma,int length)522 static uint16_t *build_pow_table(float gamma, int length)
523 {
524 int i;
525 uint16_t *output = malloc(sizeof(uint16_t)*length);
526 if (!output)
527 return NULL;
528
529 for (i = 0; i < length; i++) {
530 uint16_fract_t result;
531 double x = ((double) i) / (double) (length - 1);
532 x = pow(x, gamma); //XXX turn this conversion into a function
533 result = floor(x*65535. + .5);
534 output[i] = result;
535 }
536 return output;
537 }
538
build_output_lut(struct curveType * trc,uint16_t ** output_gamma_lut,size_t * output_gamma_lut_length)539 void build_output_lut(struct curveType *trc,
540 uint16_t **output_gamma_lut, size_t *output_gamma_lut_length)
541 {
542 if (trc->type == PARAMETRIC_CURVE_TYPE) {
543 float gamma_table[256];
544 uint16_t i;
545 uint16_t *output = malloc(sizeof(uint16_t)*256);
546
547 if (!output) {
548 *output_gamma_lut = NULL;
549 return;
550 }
551
552 compute_curve_gamma_table_type_parametric(gamma_table, trc->parameter, trc->count);
553 *output_gamma_lut_length = 256;
554 for(i = 0; i < 256; i++) {
555 output[i] = (uint16_t)(gamma_table[i] * 65535);
556 }
557 *output_gamma_lut = output;
558 } else {
559 if (trc->count == 0) {
560 *output_gamma_lut = build_linear_table(4096);
561 *output_gamma_lut_length = 4096;
562 } else if (trc->count == 1) {
563 float gamma = 1./u8Fixed8Number_to_float(trc->data[0]);
564 *output_gamma_lut = build_pow_table(gamma, 4096);
565 *output_gamma_lut_length = 4096;
566 } else {
567 //XXX: the choice of a minimum of 256 here is not backed by any theory,
568 // measurement or data, however it is what lcms uses.
569 *output_gamma_lut_length = trc->count;
570 if (*output_gamma_lut_length < 256)
571 *output_gamma_lut_length = 256;
572
573 *output_gamma_lut = invert_lut(trc->data, trc->count, *output_gamma_lut_length);
574 }
575 }
576
577 }
578