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? 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 */ 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 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? 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 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, 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 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 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 /* 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. */ 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 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 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. */ 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. */ 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 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 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 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 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 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 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 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
