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