third_party/qcms/src/transform_util.c - Issue 1438093002: Remove //third_party/qcms.

Side by Side Diff: third_party/qcms/src/transform_util.c

Issue 1438093002: Remove //third_party/qcms. (Closed) Base URL: https://github.com/domokit/mojo.git@master

Patch Set: Created 5 years, 1 month ago

Use n/p to move between diff chunks; N/P to move between comments. Draft comments are only viewable by you.

Jump to:

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, 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 - 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, size_t lengt h)

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/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, 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_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, 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 - 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 /* 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, T RC->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, u8F ixed8Number_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 cod e 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 len gth)

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 arbitrar y

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 The number of entries needed to invert a lookup table should not

399 necessarily be the same as the original number of entries. This is

400 especially true of lookup tables that have a small number of entries.

401

402 For example:

403 Using a table like:

404 {0, 3104, 14263, 34802, 65535}

405 invert_lut will produce an inverse of:

406 {3, 34459, 47529, 56801, 65535}

407 which has an maximum error of about 9855 (pixel difference of ~38.346)

408

409 For now, we punt the decision of output size to the caller. */

410 static uint16_t invert_lut(uint16_t table, int length, size_t out_length)

411 {

412 int i;

413 /* for now we invert the lut by creating a lut of size out_length

414 * and attempting to lookup a value for each entry using lut_inverse_int erp16 */

415 uint16_t output = malloc(sizeof(uint16_t)out_length);

416 if (!output)

417 return NULL;

418

419 for (i = 0; i < out_length; i++) {

420 double x = ((double) i * 65535.) / (double) (out_length - 1);

421 uint16_fract_t input = floor(x + .5);

422 output[i] = lut_inverse_interp16(input, table, length);

423 }

424 return output;

425 }

426

427 static void compute_precache_pow(uint8_t *output, float gamma)

428 {

429 uint32_t v = 0;

430 for (v = 0; v < PRECACHE_OUTPUT_SIZE; v++) {

431 //XXX: don't do integer/float conversion... and round?

432 output[v] = 255. * pow(v/(double)PRECACHE_OUTPUT_MAX, gamma);

433 }

434 }

435

436 void compute_precache_lut(uint8_t output, uint16_t table, int length)

437 {

438 uint32_t v = 0;

439 for (v = 0; v < PRECACHE_OUTPUT_SIZE; v++) {

440 output[v] = lut_interp_linear_precache_output(v, table, length);

441 }

442 }

443

444 void compute_precache_linear(uint8_t *output)

445 {

446 uint32_t v = 0;

447 for (v = 0; v < PRECACHE_OUTPUT_SIZE; v++) {

448 //XXX: round?

449 output[v] = v / (PRECACHE_OUTPUT_SIZE/256);

450 }

451 }

452

453 qcms_bool compute_precache(struct curveType trc, uint8_t output)

454 {

455

456 if (trc->type == PARAMETRIC_CURVE_TYPE) {

457 float gamma_table[256];

458 uint16_t gamma_table_uint[256];

459 uint16_t i;

460 uint16_t *inverted;

461 int inverted_size = 256;

462

463 compute_curve_gamma_table_type_parametric(gamma_table, t rc->parameter, trc->count);

464 for(i = 0; i < 256; i++) {

465 gamma_table_uint[i] = (uint16_t)(gamma_table[i] * 65535);

466 }

467

468 //XXX: the choice of a minimum of 256 here is not backed by any theory,

469 // measurement or data, howeve r it is what lcms use s.

470 // the maximum number we would need is 65535 because that's the

471 // accuracy used for computing the pre cache table

472 if (inverted_size < 256)

473 inverted_size = 256;

474

475 inverted = invert_lut(gamma_table_uint, 256, inverted_si ze);

476 if (!inverted)

477 return false;

478 compute_precache_lut(output, inverted, inverted_size);

479 free(inverted);

480 } else {

481 if (trc->count == 0) {

482 compute_precache_linear(output);

483 } else if (trc->count == 1) {

484 compute_precache_pow(output, 1./u8Fixed8Number_to_float( trc->data[0]));

485 } else {

486 uint16_t *inverted;

487 int inverted_size = trc->count;

488 //XXX: the choice of a minimum of 256 here is not backed by any theory,

489 // measurement or data, howeve r it is what lcms use s.

490 // the maximum number we would need is 65535 because that's the

491 // accuracy used for computing the pre cache table

492 if (inverted_size < 256)

493 inverted_size = 256;

494

495 inverted = invert_lut(trc->data, trc->count, inverted_si ze);

496 if (!inverted)

497 return false;

498 compute_precache_lut(output, inverted, inverted_size);

499 free(inverted);

500 }

501 }

502 return true;

503 }

504

505

506 static uint16_t *build_linear_table(int length)

507 {

508 int i;

509 uint16_t output = malloc(sizeof(uint16_t)length);

510 if (!output)

511 return NULL;

512

513 for (i = 0; i < length; i++) {

514 double x = ((double) i * 65535.) / (double) (length - 1);

515 uint16_fract_t input = floor(x + .5);

516 output[i] = input;

517 }

518 return output;

519 }

520

521 static uint16_t *build_pow_table(float gamma, int length)

522 {

523 int i;

524 uint16_t output = malloc(sizeof(uint16_t)length);

525 if (!output)

526 return NULL;

527

528 for (i = 0; i < length; i++) {

529 uint16_fract_t result;

530 double x = ((double) i) / (double) (length - 1);

531 x = pow(x, gamma); //XXX turn this conversion int o a function

532 result = floor(x*65535. + .5);

533 output[i] = result;

534 }

535 return output;

536 }

537

538 void build_output_lut(struct curveType *trc,

539 uint16_t *output_gamma_lut, size_t output_gamma_lut_length)

540 {

541 if (trc->type == PARAMETRIC_CURVE_TYPE) {

542 float gamma_table[256];

543 uint16_t i;

544 uint16_t output = malloc(sizeof(uint16_t)256);

545

546 if (!output) {

547 *output_gamma_lut = NULL;

548 return;

549 }

550

551 compute_curve_gamma_table_type_parametric(gamma_table, trc->para meter, trc->count);

552 *output_gamma_lut_length = 256;

553 for(i = 0; i < 256; i++) {

554 output[i] = (uint16_t)(gamma_table[i] * 65535);

555 }

556 *output_gamma_lut = output;

557 } else {

558 if (trc->count == 0) {

559 *output_gamma_lut = build_linear_table(4096);

560 *output_gamma_lut_length = 4096;

561 } else if (trc->count == 1) {

562 float gamma = 1./u8Fixed8Number_to_float(trc->data[0]);

563 *output_gamma_lut = build_pow_table(gamma, 4096);

564 *output_gamma_lut_length = 4096;

565 } else {

566 //XXX: the choice of a minimum of 256 here is not backed by any theory,

567 // measurement or data, however it is what lcms uses .

568 *output_gamma_lut_length = trc->count;

569 if (*output_gamma_lut_length < 256)

570 *output_gamma_lut_length = 256;

571

572 output_gamma_lut = invert_lut(trc->data, trc->count, o utput_gamma_lut_length);

573 }

574 }

575

576 }

OLD	NEW

« no previous file with comments | « third_party/qcms/src/transform_util.h ('k') | no next file » | no next file with comments »