1 #define _ISOC99_SOURCE /* for INFINITY */
5 #include <string.h> //memcpy
7 #include "transform_util.h"
10 #if !defined(INFINITY)
11 #define INFINITY HUGE_VAL
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
)
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
);
29 return value
* (1.f
/65535.f
);
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
)
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
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 */
51 uint8_t lut_interp_linear_precache_output(uint32_t input_value
, uint16_t *table
, int length
)
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)
68 value
+= (PRECACHE_OUTPUT_MAX
*65535/255)/2;
69 value
/= (PRECACHE_OUTPUT_MAX
*65535/255); // scale to 0..255
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
)
78 value
= value
* (length
- 1);
80 lower
= floorf(value
);
81 //XXX: can we be more performant here?
82 value
= table
[upper
]*(1. - (upper
- value
)) + table
[lower
]*(upper
- value
);
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
)
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
104 void compute_curve_gamma_table_type1(float gamma_table
[256], uint16_t gamma
)
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
);
114 void compute_curve_gamma_table_type2(float gamma_table
[256], uint16_t *table
, int length
)
117 for (i
= 0; i
< 256; i
++) {
118 gamma_table
[i
] = lut_interp_linear(i
/255., table
, length
);
122 void compute_curve_gamma_table_type_parametric(float gamma_table
[256], float parameter
[7], int count
)
127 float y
= parameter
[0];
134 interval
= -INFINITY
;
135 } else if(count
== 1) {
141 interval
= -1 * parameter
[2] / parameter
[1];
142 } else if(count
== 2) {
148 interval
= -1 * parameter
[2] / parameter
[1];
149 } else if(count
== 3) {
155 interval
= parameter
[4];
156 } else if(count
== 4) {
160 e
= parameter
[5] - c
;
162 interval
= parameter
[4];
164 assert(0 && "invalid parametric function type.");
170 interval
= -INFINITY
;
172 for (X
= 0; X
< 256; X
++) {
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
);
179 gamma_table
[X
] = clamp_float(c
* X
/ 255. + f
);
184 void compute_curve_gamma_table_type0(float gamma_table
[256])
187 for (i
= 0; i
< 256; i
++) {
188 gamma_table
[i
] = i
/255.;
192 float *build_input_gamma_table(struct curveType
*TRC
)
196 if (!TRC
) return NULL
;
197 gamma_table
= malloc(sizeof(float)*256);
199 if (TRC
->type
== PARAMETRIC_CURVE_TYPE
) {
200 compute_curve_gamma_table_type_parametric(gamma_table
, TRC
->parameter
, TRC
->count
);
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]);
207 compute_curve_gamma_table_type2(gamma_table
, TRC
->data
, TRC
->count
);
214 struct matrix
build_colorant_matrix(qcms_profile
*p
)
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;
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
)
238 int x
= 0, res
; // 'int' Give spacing for negative values
239 int NumZeroes
, NumPoles
;
242 double y0
, y1
, x0
, x1
;
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.
251 while (LutTable
[NumZeroes
] == 0 && NumZeroes
< length
-1)
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)
261 while (LutTable
[length
-1- NumPoles
] == 0xFFFF && NumPoles
< length
-1)
264 // Does the curve belong to this case?
265 if (NumZeroes
> 1 || NumPoles
> 1)
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);
283 // Seems not a degenerated case... apply binary search
289 res
= (int) lut_interp_linear16((uint16_fract_t
) (x
-1), LutTable
, length
);
293 // Found exact match.
295 return (uint16_fract_t
) (x
- 1);
298 if (res
> Value
) r
= x
- 1;
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
);
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);
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.
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
)
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
);
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
);
364 static void compute_precache_pow(uint8_t *output
, float gamma
)
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
);
373 void compute_precache_lut(uint8_t *output
, uint16_t *table
, int length
)
376 for (v
= 0; v
< PRECACHE_OUTPUT_SIZE
; v
++) {
377 output
[v
] = lut_interp_linear_precache_output(v
, table
, length
);
381 void compute_precache_linear(uint8_t *output
)
384 for (v
= 0; v
< PRECACHE_OUTPUT_SIZE
; v
++) {
386 output
[v
] = v
/ (PRECACHE_OUTPUT_SIZE
/256);
390 qcms_bool
compute_precache(struct curveType
*trc
, uint8_t *output
)
393 if (trc
->type
== PARAMETRIC_CURVE_TYPE
) {
394 float gamma_table
[256];
395 uint16_t gamma_table_uint
[256];
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);
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)
412 inverted
= invert_lut(gamma_table_uint
, 256, inverted_size
);
415 compute_precache_lut(output
, inverted
, inverted_size
);
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]));
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)
432 inverted
= invert_lut(trc
->data
, trc
->count
, inverted_size
);
435 compute_precache_lut(output
, inverted
, inverted_size
);
443 static uint16_t *build_linear_table(int length
)
446 uint16_t *output
= malloc(sizeof(uint16_t)*length
);
450 for (i
= 0; i
< length
; i
++) {
451 double x
= ((double) i
* 65535.) / (double) (length
- 1);
452 uint16_fract_t input
= floor(x
+ .5);
458 static uint16_t *build_pow_table(float gamma
, int length
)
461 uint16_t *output
= malloc(sizeof(uint16_t)*length
);
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);
475 void build_output_lut(struct curveType
*trc
,
476 uint16_t **output_gamma_lut
, size_t *output_gamma_lut_length
)
478 if (trc
->type
== PARAMETRIC_CURVE_TYPE
) {
479 float gamma_table
[256];
481 uint16_t *output
= malloc(sizeof(uint16_t)*256);
484 *output_gamma_lut
= NULL
;
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);
493 *output_gamma_lut
= output
;
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;
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
);