Line data Source code
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
13 6144 : float lut_interp_linear(double input_value, uint16_t *table, int length)
14 : {
15 : int upper, lower;
16 : float value;
17 6144 : input_value = input_value * (length - 1); // scale to length of the array
18 6144 : upper = ceil(input_value);
19 6144 : lower = floor(input_value);
20 : //XXX: can we be more performant here?
21 6144 : value = table[upper]*(1. - (upper - input_value)) + table[lower]*(upper - input_value);
22 : /* scale the value */
23 6144 : 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 */
27 133578 : 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 133578 : uint32_t value = (input_value * (length - 1));
32 133578 : uint32_t upper = (value + 65534) / 65535; /* equivalent to ceil(value/65535) */
33 133578 : uint32_t lower = value / 65535; /* equivalent to floor(value/65535) */
34 : /* interp is the distance from upper to value scaled to 0..65535 */
35 133578 : uint32_t interp = value % 65535;
36 :
37 133578 : value = (table[upper]*(interp) + table[lower]*(65535 - interp))/65535; // 0..65535*65535
38 :
39 133578 : 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
45 73728 : 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 73728 : uint32_t value = (input_value * (length - 1));
50 :
51 : /* equivalent to ceil(value/PRECACHE_OUTPUT_MAX) */
52 73728 : uint32_t upper = (value + PRECACHE_OUTPUT_MAX-1) / PRECACHE_OUTPUT_MAX;
53 : /* equivalent to floor(value/PRECACHE_OUTPUT_MAX) */
54 73728 : uint32_t lower = value / PRECACHE_OUTPUT_MAX;
55 : /* interp is the distance from upper to value scaled to 0..PRECACHE_OUTPUT_MAX */
56 73728 : uint32_t interp = value % PRECACHE_OUTPUT_MAX;
57 :
58 : /* the table values range from 0..65535 */
59 73728 : value = (table[upper]*(interp) + table[lower]*(PRECACHE_OUTPUT_MAX - interp)); // 0..(65535*PRECACHE_OUTPUT_MAX)
60 :
61 : /* round and scale */
62 73728 : value += (PRECACHE_OUTPUT_MAX*65535/255)/2;
63 73728 : value /= (PRECACHE_OUTPUT_MAX*65535/255); // scale to 0..255
64 73728 : return value;
65 : }
66 :
67 : /* value must be a value between 0 and 1 */
68 : //XXX: is the above a good restriction to have?
69 0 : float lut_interp_linear_float(float value, float *table, int length)
70 : {
71 : int upper, lower;
72 0 : value = value * (length - 1);
73 0 : upper = ceilf(value);
74 0 : lower = floorf(value);
75 : //XXX: can we be more performant here?
76 0 : value = table[upper]*(1. - (upper - value)) + table[lower]*(upper - value);
77 : /* scale the value */
78 0 : 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 :
98 0 : void compute_curve_gamma_table_type1(float gamma_table[256], uint16_t gamma)
99 : {
100 : unsigned int i;
101 0 : float gamma_float = u8Fixed8Number_to_float(gamma);
102 0 : for (i = 0; i < 256; i++) {
103 : // 0..1^(0..255 + 255/256) will always be between 0 and 1
104 0 : gamma_table[i] = pow(i/255., gamma_float);
105 : }
106 0 : }
107 :
108 24 : void compute_curve_gamma_table_type2(float gamma_table[256], uint16_t *table, int length)
109 : {
110 : unsigned int i;
111 6168 : for (i = 0; i < 256; i++) {
112 6144 : gamma_table[i] = lut_interp_linear(i/255., table, length);
113 : }
114 24 : }
115 :
116 0 : 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 0 : float y = parameter[0];
122 0 : if (count == 0) {
123 0 : a = 1;
124 0 : b = 0;
125 0 : c = 0;
126 0 : e = 0;
127 0 : f = 0;
128 0 : interval = -1;
129 0 : } else if(count == 1) {
130 0 : a = parameter[1];
131 0 : b = parameter[2];
132 0 : c = 0;
133 0 : e = 0;
134 0 : f = 0;
135 0 : interval = -1 * parameter[2] / parameter[1];
136 0 : } else if(count == 2) {
137 0 : a = parameter[1];
138 0 : b = parameter[2];
139 0 : c = 0;
140 0 : e = parameter[3];
141 0 : f = parameter[3];
142 0 : interval = -1 * parameter[2] / parameter[1];
143 0 : } else if(count == 3) {
144 0 : a = parameter[1];
145 0 : b = parameter[2];
146 0 : c = parameter[3];
147 0 : e = -c;
148 0 : f = 0;
149 0 : interval = parameter[4];
150 0 : } else if(count == 4) {
151 0 : a = parameter[1];
152 0 : b = parameter[2];
153 0 : c = parameter[3];
154 0 : e = parameter[5] - c;
155 0 : f = parameter[6];
156 0 : interval = parameter[4];
157 : } else {
158 0 : 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 0 : for (X = 0; X < 256; X++) {
167 0 : 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 0 : gamma_table[X] = clamp_float(pow(a * X / 255. + b, y) + c + e);
172 : } else {
173 0 : gamma_table[X] = clamp_float(c * X / 255. + f);
174 : }
175 : }
176 0 : }
177 :
178 0 : void compute_curve_gamma_table_type0(float gamma_table[256])
179 : {
180 : unsigned int i;
181 0 : for (i = 0; i < 256; i++) {
182 0 : gamma_table[i] = i/255.;
183 : }
184 0 : }
185 :
186 24 : float *build_input_gamma_table(struct curveType *TRC)
187 : {
188 : float *gamma_table;
189 :
190 24 : if (!TRC) return NULL;
191 24 : gamma_table = malloc(sizeof(float)*256);
192 24 : if (gamma_table) {
193 24 : if (TRC->type == PARAMETRIC_CURVE_TYPE) {
194 0 : compute_curve_gamma_table_type_parametric(gamma_table, TRC->parameter, TRC->count);
195 : } else {
196 24 : if (TRC->count == 0) {
197 0 : compute_curve_gamma_table_type0(gamma_table);
198 24 : } else if (TRC->count == 1) {
199 0 : compute_curve_gamma_table_type1(gamma_table, TRC->data[0]);
200 : } else {
201 24 : compute_curve_gamma_table_type2(gamma_table, TRC->data, TRC->count);
202 : }
203 : }
204 : }
205 24 : return gamma_table;
206 : }
207 :
208 16 : struct matrix build_colorant_matrix(qcms_profile *p)
209 : {
210 : struct matrix result;
211 16 : result.m[0][0] = s15Fixed16Number_to_float(p->redColorant.X);
212 16 : result.m[0][1] = s15Fixed16Number_to_float(p->greenColorant.X);
213 16 : result.m[0][2] = s15Fixed16Number_to_float(p->blueColorant.X);
214 16 : result.m[1][0] = s15Fixed16Number_to_float(p->redColorant.Y);
215 16 : result.m[1][1] = s15Fixed16Number_to_float(p->greenColorant.Y);
216 16 : result.m[1][2] = s15Fixed16Number_to_float(p->blueColorant.Y);
217 16 : result.m[2][0] = s15Fixed16Number_to_float(p->redColorant.Z);
218 16 : result.m[2][1] = s15Fixed16Number_to_float(p->greenColorant.Z);
219 16 : result.m[2][2] = s15Fixed16Number_to_float(p->blueColorant.Z);
220 16 : result.invalid = false;
221 16 : 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. */
228 9216 : uint16_fract_t lut_inverse_interp16(uint16_t Value, uint16_t LutTable[], int length)
229 : {
230 9216 : int l = 1;
231 9216 : int r = 0x10000;
232 9216 : 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 9216 : NumZeroes = 0;
245 27648 : while (LutTable[NumZeroes] == 0 && NumZeroes < length-1)
246 9216 : 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 9216 : if (NumZeroes == 0 && Value == 0)
252 0 : return 0;
253 :
254 9216 : NumPoles = 0;
255 27648 : while (LutTable[length-1- NumPoles] == 0xFFFF && NumPoles < length-1)
256 9216 : NumPoles++;
257 :
258 : // Does the curve belong to this case?
259 9216 : if (NumZeroes > 1 || NumPoles > 1)
260 : {
261 : int a, b;
262 :
263 : // Identify if value fall downto 0 or FFFF zone
264 0 : if (Value == 0) return 0;
265 : // if (Value == 0xFFFF) return 0xFFFF;
266 :
267 : // else restrict to valid zone
268 :
269 0 : if (NumZeroes > 1) {
270 0 : a = ((NumZeroes-1) * 0xFFFF) / (length-1);
271 0 : l = a - 1;
272 : }
273 0 : if (NumPoles > 1) {
274 0 : b = ((length-1 - NumPoles) * 0xFFFF) / (length-1);
275 0 : r = b + 1;
276 : }
277 : }
278 :
279 9216 : if (r <= l) {
280 : // If this happens LutTable is not invertible
281 0 : return 0;
282 : }
283 :
284 :
285 : // Seems not a degenerated case... apply binary search
286 148032 : while (r > l) {
287 :
288 133578 : x = (l + r) / 2;
289 :
290 133578 : res = (int) lut_interp_linear16((uint16_fract_t) (x-1), LutTable, length);
291 :
292 133578 : if (res == Value) {
293 :
294 : // Found exact match.
295 :
296 3978 : return (uint16_fract_t) (x - 1);
297 : }
298 :
299 129600 : if (res > Value) r = x - 1;
300 69948 : else l = x + 1;
301 : }
302 :
303 : // Not found, should we interpolate?
304 :
305 : // Get surrounding nodes
306 :
307 5238 : assert(x >= 1);
308 :
309 5238 : val2 = (length-1) * ((double) (x - 1) / 65535.0);
310 :
311 5238 : cell0 = (int) floor(val2);
312 5238 : cell1 = (int) ceil(val2);
313 :
314 5238 : if (cell0 == cell1) return (uint16_fract_t) x;
315 :
316 5238 : y0 = LutTable[cell0] ;
317 5238 : x0 = (65535.0 * cell0) / (length-1);
318 :
319 5238 : y1 = LutTable[cell1] ;
320 5238 : x1 = (65535.0 * cell1) / (length-1);
321 :
322 5238 : a = (y1 - y0) / (x1 - x0);
323 5238 : b = y0 - a * x0;
324 :
325 5238 : if (fabs(a) < 0.01) return (uint16_fract_t) x;
326 :
327 5238 : f = ((Value - b) / a);
328 :
329 5238 : if (f < 0.0) return (uint16_fract_t) 0;
330 5238 : if (f >= 65535.0) return (uint16_fract_t) 0xFFFF;
331 :
332 5229 : 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. */
349 9 : 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 9 : uint16_t *output = malloc(sizeof(uint16_t)*out_length);
355 9 : if (!output)
356 0 : return NULL;
357 :
358 9225 : for (i = 0; i < out_length; i++) {
359 9216 : double x = ((double) i * 65535.) / (double) (out_length - 1);
360 9216 : uint16_fract_t input = floor(x + .5);
361 9216 : output[i] = lut_inverse_interp16(input, table, length);
362 : }
363 9 : return output;
364 : }
365 :
366 0 : static void compute_precache_pow(uint8_t *output, float gamma)
367 : {
368 0 : uint32_t v = 0;
369 0 : for (v = 0; v < PRECACHE_OUTPUT_SIZE; v++) {
370 : //XXX: don't do integer/float conversion... and round?
371 0 : output[v] = 255. * pow(v/(double)PRECACHE_OUTPUT_MAX, gamma);
372 : }
373 0 : }
374 :
375 9 : void compute_precache_lut(uint8_t *output, uint16_t *table, int length)
376 : {
377 9 : uint32_t v = 0;
378 73737 : for (v = 0; v < PRECACHE_OUTPUT_SIZE; v++) {
379 73728 : output[v] = lut_interp_linear_precache_output(v, table, length);
380 : }
381 9 : }
382 :
383 0 : void compute_precache_linear(uint8_t *output)
384 : {
385 0 : uint32_t v = 0;
386 0 : for (v = 0; v < PRECACHE_OUTPUT_SIZE; v++) {
387 : //XXX: round?
388 0 : output[v] = v / (PRECACHE_OUTPUT_SIZE/256);
389 : }
390 0 : }
391 :
392 9 : qcms_bool compute_precache(struct curveType *trc, uint8_t *output)
393 : {
394 :
395 9 : 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 0 : int inverted_size = 256;
401 :
402 0 : compute_curve_gamma_table_type_parametric(gamma_table, trc->parameter, trc->count);
403 0 : for(i = 0; i < 256; i++) {
404 0 : 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 0 : if (inverted_size < 256)
412 0 : inverted_size = 256;
413 :
414 0 : inverted = invert_lut(gamma_table_uint, 256, inverted_size);
415 0 : if (!inverted)
416 0 : return false;
417 0 : compute_precache_lut(output, inverted, inverted_size);
418 0 : free(inverted);
419 : } else {
420 9 : if (trc->count == 0) {
421 0 : compute_precache_linear(output);
422 9 : } else if (trc->count == 1) {
423 0 : compute_precache_pow(output, 1./u8Fixed8Number_to_float(trc->data[0]));
424 : } else {
425 : uint16_t *inverted;
426 9 : 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 9 : if (inverted_size < 256)
432 0 : inverted_size = 256;
433 :
434 9 : inverted = invert_lut(trc->data, trc->count, inverted_size);
435 9 : if (!inverted)
436 0 : return false;
437 9 : compute_precache_lut(output, inverted, inverted_size);
438 9 : free(inverted);
439 : }
440 : }
441 9 : return true;
442 : }
443 :
444 :
445 0 : static uint16_t *build_linear_table(int length)
446 : {
447 : int i;
448 0 : uint16_t *output = malloc(sizeof(uint16_t)*length);
449 0 : if (!output)
450 0 : return NULL;
451 :
452 0 : for (i = 0; i < length; i++) {
453 0 : double x = ((double) i * 65535.) / (double) (length - 1);
454 0 : uint16_fract_t input = floor(x + .5);
455 0 : output[i] = input;
456 : }
457 0 : return output;
458 : }
459 :
460 0 : static uint16_t *build_pow_table(float gamma, int length)
461 : {
462 : int i;
463 0 : uint16_t *output = malloc(sizeof(uint16_t)*length);
464 0 : if (!output)
465 0 : return NULL;
466 :
467 0 : for (i = 0; i < length; i++) {
468 : uint16_fract_t result;
469 0 : double x = ((double) i) / (double) (length - 1);
470 0 : x = pow(x, gamma); //XXX turn this conversion into a function
471 0 : result = floor(x*65535. + .5);
472 0 : output[i] = result;
473 : }
474 0 : return output;
475 : }
476 :
477 0 : void build_output_lut(struct curveType *trc,
478 : uint16_t **output_gamma_lut, size_t *output_gamma_lut_length)
479 : {
480 0 : if (trc->type == PARAMETRIC_CURVE_TYPE) {
481 : float gamma_table[256];
482 : uint16_t i;
483 0 : uint16_t *output = malloc(sizeof(uint16_t)*256);
484 :
485 0 : if (!output) {
486 0 : *output_gamma_lut = NULL;
487 0 : return;
488 : }
489 :
490 0 : compute_curve_gamma_table_type_parametric(gamma_table, trc->parameter, trc->count);
491 0 : *output_gamma_lut_length = 256;
492 0 : for(i = 0; i < 256; i++) {
493 0 : output[i] = (uint16_t)(gamma_table[i] * 65535);
494 : }
495 0 : *output_gamma_lut = output;
496 : } else {
497 0 : if (trc->count == 0) {
498 0 : *output_gamma_lut = build_linear_table(4096);
499 0 : *output_gamma_lut_length = 4096;
500 0 : } else if (trc->count == 1) {
501 0 : float gamma = 1./u8Fixed8Number_to_float(trc->data[0]);
502 0 : *output_gamma_lut = build_pow_table(gamma, 4096);
503 0 : *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 0 : *output_gamma_lut_length = trc->count;
508 0 : if (*output_gamma_lut_length < 256)
509 0 : *output_gamma_lut_length = 256;
510 :
511 0 : *output_gamma_lut = invert_lut(trc->data, trc->count, *output_gamma_lut_length);
512 : }
513 : }
514 :
515 : }
516 :
|
