OLD | NEW |
| (Empty) |
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 Softwar
e | |
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? | |
39 float lut_interp_linear(double value, uint16_t *table, int 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 - valu
e); | |
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 fro
m 0..1 */ | |
52 uint16_t lut_interp_linear16(uint16_t input_value, uint16_t *table, int 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 uint32_t value = (input_value * (length - 1)); | |
57 uint32_t upper = (value + 65534) / 65535; /* equivalent to ceil(value/65
535) */ | |
58 uint32_t lower = value / 65535; /* equivalent to floor(value/6
5535) */ | |
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 | |
70 uint8_t lut_interp_linear_precache_output(uint32_t input_value, uint16_t *table,
int 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 uint32_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_OUTP
UT_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 - int
erp)); // 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? | |
94 float lut_interp_linear_float(float value, float *table, int 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 - valu
e); | |
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 c
an 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 fro
m 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/40
96) */ | |
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 | |
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 | |
131 void compute_curve_gamma_table_type2(float gamma_table[256], uint16_t *table, in
t 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 | |
139 void compute_curve_gamma_table_type_parametric(float gamma_table[256], float par
ameter[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 expressi
on. | |
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 | |
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 | |
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 | |
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 | |
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 float *build_input_gamma_table(struct curveType *TRC) | |
239 { | |
240 float *gamma_table; | |
241 | |
242 if (!TRC) return NULL; | |
243 gamma_table = malloc(sizeof(float)*256); | |
244 if (gamma_table) { | |
245 if (TRC->type == PARAMETRIC_CURVE_TYPE) { | |
246 compute_curve_gamma_table_type_parametric(gamma_table, T
RC->parameter, TRC->count); | |
247 } else { | |
248 if (TRC->count == 0) { | |
249 compute_curve_gamma_table_type0(gamma_table); | |
250 } else if (TRC->count == 1) { | |
251 compute_curve_gamma_table_type1(gamma_table, u8F
ixed8Number_to_float(TRC->data[0])); | |
252 } else { | |
253 compute_curve_gamma_table_type2(gamma_table, TRC
->data, TRC->count); | |
254 } | |
255 } | |
256 } | |
257 return gamma_table; | |
258 } | |
259 | |
260 struct matrix build_colorant_matrix(qcms_profile *p) | |
261 { | |
262 struct matrix result; | |
263 result.m[0][0] = s15Fixed16Number_to_float(p->redColorant.X); | |
264 result.m[0][1] = s15Fixed16Number_to_float(p->greenColorant.X); | |
265 result.m[0][2] = s15Fixed16Number_to_float(p->blueColorant.X); | |
266 result.m[1][0] = s15Fixed16Number_to_float(p->redColorant.Y); | |
267 result.m[1][1] = s15Fixed16Number_to_float(p->greenColorant.Y); | |
268 result.m[1][2] = s15Fixed16Number_to_float(p->blueColorant.Y); | |
269 result.m[2][0] = s15Fixed16Number_to_float(p->redColorant.Z); | |
270 result.m[2][1] = s15Fixed16Number_to_float(p->greenColorant.Z); | |
271 result.m[2][2] = s15Fixed16Number_to_float(p->blueColorant.Z); | |
272 result.invalid = false; | |
273 return result; | |
274 } | |
275 | |
276 /* The following code is copied nearly directly from lcms. | |
277 * I think it could be much better. For example, Argyll seems to have better cod
e in | |
278 * icmTable_lookup_bwd and icmTable_setup_bwd. However, for now this is a quick
way | |
279 * to a working solution and allows for easy comparing with lcms. */ | |
280 uint16_fract_t lut_inverse_interp16(uint16_t Value, uint16_t LutTable[], int len
gth) | |
281 { | |
282 int l = 1; | |
283 int r = 0x10000; | |
284 int x = 0, res; // 'int' Give spacing for negative values | |
285 int NumZeroes, NumPoles; | |
286 int cell0, cell1; | |
287 double val2; | |
288 double y0, y1, x0, x1; | |
289 double a, b, f; | |
290 | |
291 // July/27 2001 - Expanded to handle degenerated curves with an arbitrar
y | |
292 // number of elements containing 0 at the begining of the table (Zeroes) | |
293 // and another arbitrary number of poles (FFFFh) at the end. | |
294 // First the zero and pole extents are computed, then value is compared. | |
295 | |
296 NumZeroes = 0; | |
297 while (LutTable[NumZeroes] == 0 && NumZeroes < length-1) | |
298 NumZeroes++; | |
299 | |
300 // There are no zeros at the beginning and we are trying to find a zero,
so | |
301 // return anything. It seems zero would be the less destructive choice | |
302 /* I'm not sure that this makes sense, but oh well... */ | |
303 if (NumZeroes == 0 && Value == 0) | |
304 return 0; | |
305 | |
306 NumPoles = 0; | |
307 while (LutTable[length-1- NumPoles] == 0xFFFF && NumPoles < length-1) | |
308 NumPoles++; | |
309 | |
310 // Does the curve belong to this case? | |
311 if (NumZeroes > 1 || NumPoles > 1) | |
312 { | |
313 int a, b; | |
314 | |
315 // Identify if value fall downto 0 or FFFF zone | |
316 if (Value == 0) return 0; | |
317 // if (Value == 0xFFFF) return 0xFFFF; | |
318 | |
319 // else restrict to valid zone | |
320 | |
321 a = ((NumZeroes-1) * 0xFFFF) / (length-1); | |
322 b = ((length-1 - NumPoles) * 0xFFFF) / (length-1); | |
323 | |
324 l = a - 1; | |
325 r = b + 1; | |
326 } | |
327 | |
328 | |
329 // Seems not a degenerated case... apply binary search | |
330 | |
331 while (r > l) { | |
332 | |
333 x = (l + r) / 2; | |
334 | |
335 res = (int) lut_interp_linear16((uint16_fract_t) (x-1), LutTable
, length); | |
336 | |
337 if (res == Value) { | |
338 | |
339 // Found exact match. | |
340 | |
341 return (uint16_fract_t) (x - 1); | |
342 } | |
343 | |
344 if (res > Value) r = x - 1; | |
345 else l = x + 1; | |
346 } | |
347 | |
348 // Not found, should we interpolate? | |
349 | |
350 | |
351 // Get surrounding nodes | |
352 | |
353 val2 = (length-1) * ((double) (x - 1) / 65535.0); | |
354 | |
355 cell0 = (int) floor(val2); | |
356 cell1 = (int) ceil(val2); | |
357 | |
358 if (cell0 == cell1) return (uint16_fract_t) x; | |
359 | |
360 y0 = LutTable[cell0] ; | |
361 x0 = (65535.0 * cell0) / (length-1); | |
362 | |
363 y1 = LutTable[cell1] ; | |
364 x1 = (65535.0 * cell1) / (length-1); | |
365 | |
366 a = (y1 - y0) / (x1 - x0); | |
367 b = y0 - a * x0; | |
368 | |
369 if (fabs(a) < 0.01) return (uint16_fract_t) x; | |
370 | |
371 f = ((Value - b) / a); | |
372 | |
373 if (f < 0.0) return (uint16_fract_t) 0; | |
374 if (f >= 65535.0) return (uint16_fract_t) 0xFFFF; | |
375 | |
376 return (uint16_fract_t) floor(f + 0.5); | |
377 | |
378 } | |
379 | |
380 /* | |
381 The number of entries needed to invert a lookup table should not | |
382 necessarily be the same as the original number of entries. This is | |
383 especially true of lookup tables that have a small number of entries. | |
384 | |
385 For example: | |
386 Using a table like: | |
387 {0, 3104, 14263, 34802, 65535} | |
388 invert_lut will produce an inverse of: | |
389 {3, 34459, 47529, 56801, 65535} | |
390 which has an maximum error of about 9855 (pixel difference of ~38.346) | |
391 | |
392 For now, we punt the decision of output size to the caller. */ | |
393 static uint16_t *invert_lut(uint16_t *table, int length, int out_length) | |
394 { | |
395 int i; | |
396 /* for now we invert the lut by creating a lut of size out_length | |
397 * and attempting to lookup a value for each entry using lut_inverse_int
erp16 */ | |
398 uint16_t *output = malloc(sizeof(uint16_t)*out_length); | |
399 if (!output) | |
400 return NULL; | |
401 | |
402 for (i = 0; i < out_length; i++) { | |
403 double x = ((double) i * 65535.) / (double) (out_length - 1); | |
404 uint16_fract_t input = floor(x + .5); | |
405 output[i] = lut_inverse_interp16(input, table, length); | |
406 } | |
407 return output; | |
408 } | |
409 | |
410 static void compute_precache_pow(uint8_t *output, float gamma) | |
411 { | |
412 uint32_t v = 0; | |
413 for (v = 0; v < PRECACHE_OUTPUT_SIZE; v++) { | |
414 //XXX: don't do integer/float conversion... and round? | |
415 output[v] = 255. * pow(v/(double)PRECACHE_OUTPUT_MAX, gamma); | |
416 } | |
417 } | |
418 | |
419 void compute_precache_lut(uint8_t *output, uint16_t *table, int length) | |
420 { | |
421 uint32_t v = 0; | |
422 for (v = 0; v < PRECACHE_OUTPUT_SIZE; v++) { | |
423 output[v] = lut_interp_linear_precache_output(v, table, length); | |
424 } | |
425 } | |
426 | |
427 void compute_precache_linear(uint8_t *output) | |
428 { | |
429 uint32_t v = 0; | |
430 for (v = 0; v < PRECACHE_OUTPUT_SIZE; v++) { | |
431 //XXX: round? | |
432 output[v] = v / (PRECACHE_OUTPUT_SIZE/256); | |
433 } | |
434 } | |
435 | |
436 qcms_bool compute_precache(struct curveType *trc, uint8_t *output) | |
437 { | |
438 | |
439 if (trc->type == PARAMETRIC_CURVE_TYPE) { | |
440 float gamma_table[256]; | |
441 uint16_t gamma_table_uint[256]; | |
442 uint16_t i; | |
443 uint16_t *inverted; | |
444 int inverted_size = 256; | |
445 | |
446 compute_curve_gamma_table_type_parametric(gamma_table, t
rc->parameter, trc->count); | |
447 for(i = 0; i < 256; i++) { | |
448 gamma_table_uint[i] = (uint16_t)(gamma_table[i]
* 65535); | |
449 } | |
450 | |
451 //XXX: the choice of a minimum of 256 here is not backed
by any theory, | |
452 // measurement or data, howeve r it is what lcms use
s. | |
453 // the maximum number we would need is 65535 because
that's the | |
454 // accuracy used for computing the pre cache table | |
455 if (inverted_size < 256) | |
456 inverted_size = 256; | |
457 | |
458 inverted = invert_lut(gamma_table_uint, 256, inverted_si
ze); | |
459 if (!inverted) | |
460 return false; | |
461 compute_precache_lut(output, inverted, inverted_size); | |
462 free(inverted); | |
463 } else { | |
464 if (trc->count == 0) { | |
465 compute_precache_linear(output); | |
466 } else if (trc->count == 1) { | |
467 compute_precache_pow(output, 1./u8Fixed8Number_to_float(
trc->data[0])); | |
468 } else { | |
469 uint16_t *inverted; | |
470 int inverted_size = trc->count; | |
471 //XXX: the choice of a minimum of 256 here is not backed
by any theory, | |
472 // measurement or data, howeve r it is what lcms use
s. | |
473 // the maximum number we would need is 65535 because
that's the | |
474 // accuracy used for computing the pre cache table | |
475 if (inverted_size < 256) | |
476 inverted_size = 256; | |
477 | |
478 inverted = invert_lut(trc->data, trc->count, inverted_si
ze); | |
479 if (!inverted) | |
480 return false; | |
481 compute_precache_lut(output, inverted, inverted_size); | |
482 free(inverted); | |
483 } | |
484 } | |
485 return true; | |
486 } | |
487 | |
488 | |
489 static uint16_t *build_linear_table(int length) | |
490 { | |
491 int i; | |
492 uint16_t *output = malloc(sizeof(uint16_t)*length); | |
493 if (!output) | |
494 return NULL; | |
495 | |
496 for (i = 0; i < length; i++) { | |
497 double x = ((double) i * 65535.) / (double) (length - 1); | |
498 uint16_fract_t input = floor(x + .5); | |
499 output[i] = input; | |
500 } | |
501 return output; | |
502 } | |
503 | |
504 static uint16_t *build_pow_table(float gamma, int length) | |
505 { | |
506 int i; | |
507 uint16_t *output = malloc(sizeof(uint16_t)*length); | |
508 if (!output) | |
509 return NULL; | |
510 | |
511 for (i = 0; i < length; i++) { | |
512 uint16_fract_t result; | |
513 double x = ((double) i) / (double) (length - 1); | |
514 x = pow(x, gamma); //XXX turn this conversion int
o a function | |
515 result = floor(x*65535. + .5); | |
516 output[i] = result; | |
517 } | |
518 return output; | |
519 } | |
520 | |
521 void build_output_lut(struct curveType *trc, | |
522 uint16_t **output_gamma_lut, size_t *output_gamma_lut_length) | |
523 { | |
524 if (trc->type == PARAMETRIC_CURVE_TYPE) { | |
525 float gamma_table[256]; | |
526 uint16_t i; | |
527 uint16_t *output = malloc(sizeof(uint16_t)*256); | |
528 | |
529 if (!output) { | |
530 *output_gamma_lut = NULL; | |
531 return; | |
532 } | |
533 | |
534 compute_curve_gamma_table_type_parametric(gamma_table, trc->para
meter, trc->count); | |
535 *output_gamma_lut_length = 256; | |
536 for(i = 0; i < 256; i++) { | |
537 output[i] = (uint16_t)(gamma_table[i] * 65535); | |
538 } | |
539 *output_gamma_lut = output; | |
540 } else { | |
541 if (trc->count == 0) { | |
542 *output_gamma_lut = build_linear_table(4096); | |
543 *output_gamma_lut_length = 4096; | |
544 } else if (trc->count == 1) { | |
545 float gamma = 1./u8Fixed8Number_to_float(trc->data[0]); | |
546 *output_gamma_lut = build_pow_table(gamma, 4096); | |
547 *output_gamma_lut_length = 4096; | |
548 } else { | |
549 //XXX: the choice of a minimum of 256 here is not backed
by any theory, | |
550 // measurement or data, however it is what lcms uses
. | |
551 *output_gamma_lut_length = trc->count; | |
552 if (*output_gamma_lut_length < 256) | |
553 *output_gamma_lut_length = 256; | |
554 | |
555 *output_gamma_lut = invert_lut(trc->data, trc->count, *o
utput_gamma_lut_length); | |
556 } | |
557 } | |
558 | |
559 } | |
OLD | NEW |