1 /* Copyright (C) 2001-2019 Artifex Software, Inc.
2 All Rights Reserved.
3
4 This software is provided AS-IS with no warranty, either express or
5 implied.
6
7 This software is distributed under license and may not be copied,
8 modified or distributed except as expressly authorized under the terms
9 of the license contained in the file LICENSE in this distribution.
10
11 Refer to licensing information at http://www.artifex.com or contact
12 Artifex Software, Inc., 1305 Grant Avenue - Suite 200, Novato,
13 CA 94945, U.S.A., +1(415)492-9861, for further information.
14 */
15
16
17 /* CIE color rendering cache management */
18 #include "math_.h"
19 #include "memory_.h"
20 #include "gx.h"
21 #include "gserrors.h"
22 #include "gsstruct.h"
23 #include "gsmatrix.h" /* for gscolor2.h */
24 #include "gxcspace.h" /* for gxcie.c */
25 #include "gscolor2.h" /* for gs_set/currentcolorrendering */
26 #include "gxarith.h"
27 #include "gxcie.h"
28 #include "gxdevice.h" /* for gxcmap.h */
29 #include "gxcmap.h"
30 #include "gzstate.h"
31 #include "gsicc.h"
32
33 /*
34 * Define whether to optimize the CIE mapping process by combining steps.
35 * This should only be disabled (commented out) for debugging.
36 */
37 #define OPTIMIZE_CIE_MAPPING
38
39 /* Forward references */
40 static int cie_joint_caches_init(gx_cie_joint_caches *,
41 const gs_cie_common *,
42 gs_cie_render *);
43 static void cie_joint_caches_complete(gx_cie_joint_caches *,
44 const gs_cie_common *,
45 const gs_cie_abc *,
46 const gs_cie_render *);
47 static void cie_cache_restrict(cie_cache_floats *, const gs_range *);
48 static void cie_invert3(const gs_matrix3 *, gs_matrix3 *);
49 static void cie_matrix_init(gs_matrix3 *);
50
51 /* Allocator structure types */
52 private_st_joint_caches();
53 extern_st(st_gs_gstate);
54
55 #define RESTRICTED_INDEX(v, n, itemp)\
56 ((uint)(itemp = (int)(v)) >= (n) ?\
57 (itemp < 0 ? 0 : (n) - 1) : itemp)
58
59 /* Define cache interpolation threshold values. */
60 #ifdef CIE_CACHE_INTERPOLATE
61 # ifdef CIE_INTERPOLATE_THRESHOLD
62 # define CACHE_THRESHOLD CIE_INTERPOLATE_THRESHOLD
63 # else
64 # define CACHE_THRESHOLD 0 /* always interpolate */
65 # endif
66 #else
67 # define CACHE_THRESHOLD 1.0e6 /* never interpolate */
68 #endif
69 #ifdef CIE_RENDER_TABLE_INTERPOLATE
70 # define RENDER_TABLE_THRESHOLD 0
71 #else
72 # define RENDER_TABLE_THRESHOLD 1.0e6
73 #endif
74
75 /*
76 * Determine whether a function is a linear transformation of the form
77 * f(x) = scale * x + origin.
78 */
79 static bool
cache_is_linear(cie_linear_params_t * params,const cie_cache_floats * pcf)80 cache_is_linear(cie_linear_params_t *params, const cie_cache_floats *pcf)
81 {
82 double origin = pcf->values[0];
83 double diff = pcf->values[countof(pcf->values) - 1] - origin;
84 double scale = diff / (countof(pcf->values) - 1);
85 int i;
86 double test = origin + scale;
87
88 for (i = 1; i < countof(pcf->values) - 1; ++i, test += scale)
89 if (fabs(pcf->values[i] - test) >= 0.5 / countof(pcf->values))
90 return (params->is_linear = false);
91 params->origin = origin - pcf->params.base;
92 params->scale =
93 diff * pcf->params.factor / (countof(pcf->values) - 1);
94 return (params->is_linear = true);
95 }
96
97 static void
cache_set_linear(cie_cache_floats * pcf)98 cache_set_linear(cie_cache_floats *pcf)
99 {
100 if (pcf->params.is_identity) {
101 if_debug1('c', "[c]is_linear(0x%lx) = true (is_identity)\n",
102 (ulong)pcf);
103 pcf->params.linear.is_linear = true;
104 pcf->params.linear.origin = 0;
105 pcf->params.linear.scale = 1;
106 } else if (cache_is_linear(&pcf->params.linear, pcf)) {
107 if (pcf->params.linear.origin == 0 &&
108 fabs(pcf->params.linear.scale - 1) < 0.00001)
109 pcf->params.is_identity = true;
110 if_debug4('c',
111 "[c]is_linear(0x%lx) = true, origin = %g, scale = %g%s\n",
112 (ulong)pcf, pcf->params.linear.origin,
113 pcf->params.linear.scale,
114 (pcf->params.is_identity ? " (=> is_identity)" : ""));
115 }
116 #ifdef DEBUG
117 else
118 if_debug1('c', "[c]linear(0x%lx) = false\n", (ulong)pcf);
119 #endif
120 }
121 static void
cache3_set_linear(gx_cie_vector_cache3_t * pvc)122 cache3_set_linear(gx_cie_vector_cache3_t *pvc)
123 {
124 cache_set_linear(&pvc->caches[0].floats);
125 cache_set_linear(&pvc->caches[1].floats);
126 cache_set_linear(&pvc->caches[2].floats);
127 }
128
129 #ifdef DEBUG
130 static void
if_debug_vector3(const char * str,const gs_vector3 * vec)131 if_debug_vector3(const char *str, const gs_vector3 *vec)
132 {
133 if_debug4('c', "%s[%g %g %g]\n", str, vec->u, vec->v, vec->w);
134 }
135 static void
if_debug_matrix3(const char * str,const gs_matrix3 * mat)136 if_debug_matrix3(const char *str, const gs_matrix3 *mat)
137 {
138 if_debug10('c', "%s [%g %g %g] [%g %g %g] [%g %g %g]\n", str,
139 mat->cu.u, mat->cu.v, mat->cu.w,
140 mat->cv.u, mat->cv.v, mat->cv.w,
141 mat->cw.u, mat->cw.v, mat->cw.w);
142 }
143 #else
144 # define if_debug_vector3(str, vec) DO_NOTHING
145 # define if_debug_matrix3(str, mat) DO_NOTHING
146 #endif
147
148 /* ------ Default values for CIE dictionary elements ------ */
149
150 /* Default transformation procedures. */
151
152 float
a_identity(double in,const gs_cie_a * pcie)153 a_identity(double in, const gs_cie_a * pcie)
154 {
155 return in;
156 }
157 static float
a_from_cache(double in,const gs_cie_a * pcie)158 a_from_cache(double in, const gs_cie_a * pcie)
159 {
160 return gs_cie_cached_value(in, &pcie->caches.DecodeA.floats);
161 }
162
163 float
abc_identity(double in,const gs_cie_abc * pcie)164 abc_identity(double in, const gs_cie_abc * pcie)
165 {
166 return in;
167 }
168 static float
abc_from_cache_0(double in,const gs_cie_abc * pcie)169 abc_from_cache_0(double in, const gs_cie_abc * pcie)
170 {
171 return gs_cie_cached_value(in, &pcie->caches.DecodeABC.caches[0].floats);
172 }
173 static float
abc_from_cache_1(double in,const gs_cie_abc * pcie)174 abc_from_cache_1(double in, const gs_cie_abc * pcie)
175 {
176 return gs_cie_cached_value(in, &pcie->caches.DecodeABC.caches[1].floats);
177 }
178 static float
abc_from_cache_2(double in,const gs_cie_abc * pcie)179 abc_from_cache_2(double in, const gs_cie_abc * pcie)
180 {
181 return gs_cie_cached_value(in, &pcie->caches.DecodeABC.caches[2].floats);
182 }
183
184 static float
def_identity(double in,const gs_cie_def * pcie)185 def_identity(double in, const gs_cie_def * pcie)
186 {
187 return in;
188 }
189 static float
def_from_cache_0(double in,const gs_cie_def * pcie)190 def_from_cache_0(double in, const gs_cie_def * pcie)
191 {
192 return gs_cie_cached_value(in, &pcie->caches_def.DecodeDEF[0].floats);
193 }
194 static float
def_from_cache_1(double in,const gs_cie_def * pcie)195 def_from_cache_1(double in, const gs_cie_def * pcie)
196 {
197 return gs_cie_cached_value(in, &pcie->caches_def.DecodeDEF[1].floats);
198 }
199 static float
def_from_cache_2(double in,const gs_cie_def * pcie)200 def_from_cache_2(double in, const gs_cie_def * pcie)
201 {
202 return gs_cie_cached_value(in, &pcie->caches_def.DecodeDEF[2].floats);
203 }
204
205 static float
defg_identity(double in,const gs_cie_defg * pcie)206 defg_identity(double in, const gs_cie_defg * pcie)
207 {
208 return in;
209 }
210 static float
defg_from_cache_0(double in,const gs_cie_defg * pcie)211 defg_from_cache_0(double in, const gs_cie_defg * pcie)
212 {
213 return gs_cie_cached_value(in, &pcie->caches_defg.DecodeDEFG[0].floats);
214 }
215 static float
defg_from_cache_1(double in,const gs_cie_defg * pcie)216 defg_from_cache_1(double in, const gs_cie_defg * pcie)
217 {
218 return gs_cie_cached_value(in, &pcie->caches_defg.DecodeDEFG[1].floats);
219 }
220 static float
defg_from_cache_2(double in,const gs_cie_defg * pcie)221 defg_from_cache_2(double in, const gs_cie_defg * pcie)
222 {
223 return gs_cie_cached_value(in, &pcie->caches_defg.DecodeDEFG[2].floats);
224 }
225 static float
defg_from_cache_3(double in,const gs_cie_defg * pcie)226 defg_from_cache_3(double in, const gs_cie_defg * pcie)
227 {
228 return gs_cie_cached_value(in, &pcie->caches_defg.DecodeDEFG[3].floats);
229 }
230
231 float
common_identity(double in,const gs_cie_common * pcie)232 common_identity(double in, const gs_cie_common * pcie)
233 {
234 return in;
235 }
236 static float
lmn_from_cache_0(double in,const gs_cie_common * pcie)237 lmn_from_cache_0(double in, const gs_cie_common * pcie)
238 {
239 return gs_cie_cached_value(in, &pcie->caches.DecodeLMN[0].floats);
240 }
241 static float
lmn_from_cache_1(double in,const gs_cie_common * pcie)242 lmn_from_cache_1(double in, const gs_cie_common * pcie)
243 {
244 return gs_cie_cached_value(in, &pcie->caches.DecodeLMN[1].floats);
245 }
246 static float
lmn_from_cache_2(double in,const gs_cie_common * pcie)247 lmn_from_cache_2(double in, const gs_cie_common * pcie)
248 {
249 return gs_cie_cached_value(in, &pcie->caches.DecodeLMN[2].floats);
250 }
251
252 /* Transformation procedures for accessing an already-loaded cache. */
253
254 float
gs_cie_cached_value(double in,const cie_cache_floats * pcache)255 gs_cie_cached_value(double in, const cie_cache_floats *pcache)
256 {
257 /*
258 * We need to get the same results when we sample an already-loaded
259 * cache, so we need to round the index just a tiny bit.
260 */
261 int index =
262 (int)((in - pcache->params.base) * pcache->params.factor + 0.0001);
263
264 CIE_CLAMP_INDEX(index);
265 return pcache->values[index];
266 }
267
268 /* Default vectors and matrices. */
269
270 const gs_range3 Range3_default = {
271 { {0, 1}, {0, 1}, {0, 1} }
272 };
273 const gs_range4 Range4_default = {
274 { {0, 1}, {0, 1}, {0, 1}, {0, 1} }
275 };
276 const gs_cie_defg_proc4 DecodeDEFG_default = {
277 {defg_identity, defg_identity, defg_identity, defg_identity}
278 };
279 const gs_cie_defg_proc4 DecodeDEFG_from_cache = {
280 {defg_from_cache_0, defg_from_cache_1, defg_from_cache_2, defg_from_cache_3}
281 };
282 const gs_cie_def_proc3 DecodeDEF_default = {
283 {def_identity, def_identity, def_identity}
284 };
285 const gs_cie_def_proc3 DecodeDEF_from_cache = {
286 {def_from_cache_0, def_from_cache_1, def_from_cache_2}
287 };
288 const gs_cie_abc_proc3 DecodeABC_default = {
289 {abc_identity, abc_identity, abc_identity}
290 };
291 const gs_cie_abc_proc3 DecodeABC_from_cache = {
292 {abc_from_cache_0, abc_from_cache_1, abc_from_cache_2}
293 };
294 const gs_cie_common_proc3 DecodeLMN_default = {
295 {common_identity, common_identity, common_identity}
296 };
297 const gs_cie_common_proc3 DecodeLMN_from_cache = {
298 {lmn_from_cache_0, lmn_from_cache_1, lmn_from_cache_2}
299 };
300 const gs_matrix3 Matrix3_default = {
301 {1, 0, 0},
302 {0, 1, 0},
303 {0, 0, 1},
304 1 /*true */
305 };
306 const gs_range RangeA_default = {0, 1};
307 const gs_cie_a_proc DecodeA_default = a_identity;
308 const gs_cie_a_proc DecodeA_from_cache = a_from_cache;
309 const gs_vector3 MatrixA_default = {1, 1, 1};
310 const gs_vector3 BlackPoint_default = {0, 0, 0};
311
312 /* Initialize a CIE color. */
313 /* This only happens on setcolorspace. */
314 void
gx_init_CIE(gs_client_color * pcc,const gs_color_space * pcs)315 gx_init_CIE(gs_client_color * pcc, const gs_color_space * pcs)
316 {
317 gx_init_paint_4(pcc, pcs);
318 /* (0...) may not be within the range of allowable values. */
319 (*pcs->type->restrict_color)(pcc, pcs);
320 }
321
322 /* Restrict CIE colors. */
323
324 static inline void
cie_restrict(float * pv,const gs_range * range)325 cie_restrict(float *pv, const gs_range *range)
326 {
327 if (*pv <= range->rmin)
328 *pv = range->rmin;
329 else if (*pv >= range->rmax)
330 *pv = range->rmax;
331 }
332
333 void
gx_restrict_CIEDEFG(gs_client_color * pcc,const gs_color_space * pcs)334 gx_restrict_CIEDEFG(gs_client_color * pcc, const gs_color_space * pcs)
335 {
336 const gs_cie_defg *pcie = pcs->params.defg;
337
338 cie_restrict(&pcc->paint.values[0], &pcie->RangeDEFG.ranges[0]);
339 cie_restrict(&pcc->paint.values[1], &pcie->RangeDEFG.ranges[1]);
340 cie_restrict(&pcc->paint.values[2], &pcie->RangeDEFG.ranges[2]);
341 cie_restrict(&pcc->paint.values[3], &pcie->RangeDEFG.ranges[3]);
342 }
343 void
gx_restrict_CIEDEF(gs_client_color * pcc,const gs_color_space * pcs)344 gx_restrict_CIEDEF(gs_client_color * pcc, const gs_color_space * pcs)
345 {
346 const gs_cie_def *pcie = pcs->params.def;
347
348 cie_restrict(&pcc->paint.values[0], &pcie->RangeDEF.ranges[0]);
349 cie_restrict(&pcc->paint.values[1], &pcie->RangeDEF.ranges[1]);
350 cie_restrict(&pcc->paint.values[2], &pcie->RangeDEF.ranges[2]);
351 }
352 void
gx_restrict_CIEABC(gs_client_color * pcc,const gs_color_space * pcs)353 gx_restrict_CIEABC(gs_client_color * pcc, const gs_color_space * pcs)
354 {
355 const gs_cie_abc *pcie = pcs->params.abc;
356
357 cie_restrict(&pcc->paint.values[0], &pcie->RangeABC.ranges[0]);
358 cie_restrict(&pcc->paint.values[1], &pcie->RangeABC.ranges[1]);
359 cie_restrict(&pcc->paint.values[2], &pcie->RangeABC.ranges[2]);
360 }
361 void
gx_restrict_CIEA(gs_client_color * pcc,const gs_color_space * pcs)362 gx_restrict_CIEA(gs_client_color * pcc, const gs_color_space * pcs)
363 {
364 const gs_cie_a *pcie = pcs->params.a;
365
366 cie_restrict(&pcc->paint.values[0], &pcie->RangeA);
367 }
368
369 /* ================ Table setup ================ */
370
371 /* ------ Install a CIE color space ------ */
372
373 static void cie_cache_mult(gx_cie_vector_cache *, const gs_vector3 *,
374 const cie_cache_floats *, double);
375 static bool cie_cache_mult3(gx_cie_vector_cache3_t *,
376 const gs_matrix3 *, double);
377
378 int
gx_install_cie_abc(gs_cie_abc * pcie,gs_gstate * pgs)379 gx_install_cie_abc(gs_cie_abc *pcie, gs_gstate * pgs)
380 {
381 if_debug_matrix3("[c]CIE MatrixABC =", &pcie->MatrixABC);
382 cie_matrix_init(&pcie->MatrixABC);
383 CIE_LOAD_CACHE_BODY(pcie->caches.DecodeABC.caches, pcie->RangeABC.ranges,
384 &pcie->DecodeABC, DecodeABC_default, pcie,
385 "DecodeABC");
386 gx_cie_load_common_cache(&pcie->common, pgs);
387 gs_cie_abc_complete(pcie);
388 return gs_cie_cs_complete(pgs, true);
389 }
390
391 int
gx_install_CIEDEFG(gs_color_space * pcs,gs_gstate * pgs)392 gx_install_CIEDEFG(gs_color_space * pcs, gs_gstate * pgs)
393 {
394 gs_cie_defg *pcie = pcs->params.defg;
395
396 CIE_LOAD_CACHE_BODY(pcie->caches_defg.DecodeDEFG, pcie->RangeDEFG.ranges,
397 &pcie->DecodeDEFG, DecodeDEFG_default, pcie,
398 "DecodeDEFG");
399 return gx_install_cie_abc((gs_cie_abc *)pcie, pgs);
400 }
401
402 int
gx_install_CIEDEF(gs_color_space * pcs,gs_gstate * pgs)403 gx_install_CIEDEF(gs_color_space * pcs, gs_gstate * pgs)
404 {
405 gs_cie_def *pcie = pcs->params.def;
406
407 CIE_LOAD_CACHE_BODY(pcie->caches_def.DecodeDEF, pcie->RangeDEF.ranges,
408 &pcie->DecodeDEF, DecodeDEF_default, pcie,
409 "DecodeDEF");
410 return gx_install_cie_abc((gs_cie_abc *)pcie, pgs);
411 }
412
413 int
gx_install_CIEABC(gs_color_space * pcs,gs_gstate * pgs)414 gx_install_CIEABC(gs_color_space * pcs, gs_gstate * pgs)
415 {
416 return gx_install_cie_abc(pcs->params.abc, pgs);
417 }
418
419 int
gx_install_CIEA(gs_color_space * pcs,gs_gstate * pgs)420 gx_install_CIEA(gs_color_space * pcs, gs_gstate * pgs)
421 {
422 gs_cie_a *pcie = pcs->params.a;
423 gs_sample_loop_params_t lp;
424 int i;
425
426 gs_cie_cache_init(&pcie->caches.DecodeA.floats.params, &lp,
427 &pcie->RangeA, "DecodeA");
428 for (i = 0; i <= lp.N; ++i) {
429 float in = SAMPLE_LOOP_VALUE(i, lp);
430
431 pcie->caches.DecodeA.floats.values[i] = (*pcie->DecodeA)(in, pcie);
432 if_debug3m('C', pgs->memory, "[C]DecodeA[%d] = %g => %g\n",
433 i, in, pcie->caches.DecodeA.floats.values[i]);
434 }
435 gx_cie_load_common_cache(&pcie->common, pgs);
436 gs_cie_a_complete(pcie);
437 return gs_cie_cs_complete(pgs, true);
438 }
439
440 /* Load the common caches when installing the color space. */
441 /* This routine is exported for the benefit of gsicc.c */
442 void
gx_cie_load_common_cache(gs_cie_common * pcie,gs_gstate * pgs)443 gx_cie_load_common_cache(gs_cie_common * pcie, gs_gstate * pgs)
444 {
445 if_debug_matrix3("[c]CIE MatrixLMN =", &pcie->MatrixLMN);
446 cie_matrix_init(&pcie->MatrixLMN);
447 CIE_LOAD_CACHE_BODY(pcie->caches.DecodeLMN, pcie->RangeLMN.ranges,
448 &pcie->DecodeLMN, DecodeLMN_default, pcie,
449 "DecodeLMN");
450 }
451
452 /* Complete loading the common caches. */
453 /* This routine is exported for the benefit of gsicc.c */
454 void
gx_cie_common_complete(gs_cie_common * pcie)455 gx_cie_common_complete(gs_cie_common *pcie)
456 {
457 int i;
458
459 for (i = 0; i < 3; ++i)
460 cache_set_linear(&pcie->caches.DecodeLMN[i].floats);
461 }
462
463 /*
464 * Restrict the DecodeDEF[G] cache according to RangeHIJ[K], and scale to
465 * the dimensions of Table.
466 */
467 static void
gs_cie_defx_scale(float * values,const gs_range * range,int dim)468 gs_cie_defx_scale(float *values, const gs_range *range, int dim)
469 {
470 double scale = (dim - 1.0) / (range->rmax - range->rmin);
471 int i;
472
473 for (i = 0; i < gx_cie_cache_size; ++i) {
474 float value = values[i];
475
476 values[i] =
477 (value <= range->rmin ? 0 :
478 value >= range->rmax ? dim - 1 :
479 (value - range->rmin) * scale);
480 }
481 }
482
483 /* Complete loading a CIEBasedDEFG color space. */
484 /* This routine is NOT idempotent. */
485 void
gs_cie_defg_complete(gs_cie_defg * pcie)486 gs_cie_defg_complete(gs_cie_defg * pcie)
487 {
488 int j;
489
490 for (j = 0; j < 4; ++j)
491 gs_cie_defx_scale(pcie->caches_defg.DecodeDEFG[j].floats.values,
492 &pcie->RangeHIJK.ranges[j], pcie->Table.dims[j]);
493 gs_cie_abc_complete((gs_cie_abc *)pcie);
494 }
495
496 /* Complete loading a CIEBasedDEF color space. */
497 /* This routine is NOT idempotent. */
498 void
gs_cie_def_complete(gs_cie_def * pcie)499 gs_cie_def_complete(gs_cie_def * pcie)
500 {
501 int j;
502
503 for (j = 0; j < 3; ++j)
504 gs_cie_defx_scale(pcie->caches_def.DecodeDEF[j].floats.values,
505 &pcie->RangeHIJ.ranges[j], pcie->Table.dims[j]);
506 gs_cie_abc_complete((gs_cie_abc *)pcie);
507 }
508
509 /* Complete loading a CIEBasedABC color space. */
510 /* This routine is idempotent. */
511 void
gs_cie_abc_complete(gs_cie_abc * pcie)512 gs_cie_abc_complete(gs_cie_abc * pcie)
513 {
514 cache3_set_linear(&pcie->caches.DecodeABC);
515 pcie->caches.skipABC =
516 cie_cache_mult3(&pcie->caches.DecodeABC, &pcie->MatrixABC,
517 CACHE_THRESHOLD);
518 gx_cie_common_complete((gs_cie_common *)pcie);
519 }
520
521 /* Complete loading a CIEBasedA color space. */
522 /* This routine is idempotent. */
523 void
gs_cie_a_complete(gs_cie_a * pcie)524 gs_cie_a_complete(gs_cie_a * pcie)
525 {
526 cie_cache_mult(&pcie->caches.DecodeA, &pcie->MatrixA,
527 &pcie->caches.DecodeA.floats,
528 CACHE_THRESHOLD);
529 cache_set_linear(&pcie->caches.DecodeA.floats);
530 gx_cie_common_complete((gs_cie_common *)pcie);
531 }
532
533 /*
534 * Set the ranges where interpolation is required in a vector cache.
535 * This procedure is idempotent.
536 */
537 typedef struct cie_cache_range_temp_s {
538 cie_cached_value prev;
539 int imin, imax;
540 } cie_cache_range_temp_t;
541 static inline void
check_interpolation_required(cie_cache_range_temp_t * pccr,cie_cached_value cur,int i,double threshold)542 check_interpolation_required(cie_cache_range_temp_t *pccr,
543 cie_cached_value cur, int i, double threshold)
544 {
545 cie_cached_value prev = pccr->prev;
546
547 if (cie_cached_abs(cur - prev) > threshold * min(cie_cached_abs(prev), cie_cached_abs(cur))) {
548 if (i - 1 < pccr->imin)
549 pccr->imin = i - 1;
550 if (i > pccr->imax)
551 pccr->imax = i;
552 }
553 pccr->prev = cur;
554 }
555 static void
cie_cache_set_interpolation(gx_cie_vector_cache * pcache,double threshold)556 cie_cache_set_interpolation(gx_cie_vector_cache *pcache, double threshold)
557 {
558 cie_cached_value base = pcache->vecs.params.base;
559 cie_cached_value factor = pcache->vecs.params.factor;
560 cie_cache_range_temp_t temp[3];
561 int i, j;
562
563 for (j = 0; j < 3; ++j)
564 temp[j].imin = gx_cie_cache_size, temp[j].imax = -1;
565 temp[0].prev = pcache->vecs.values[0].u;
566 temp[1].prev = pcache->vecs.values[0].v;
567 temp[2].prev = pcache->vecs.values[0].w;
568
569 for (i = 0; i < gx_cie_cache_size; ++i) {
570 check_interpolation_required(&temp[0], pcache->vecs.values[i].u, i,
571 threshold);
572 check_interpolation_required(&temp[1], pcache->vecs.values[i].v, i,
573 threshold);
574 check_interpolation_required(&temp[2], pcache->vecs.values[i].w, i,
575 threshold);
576 }
577
578 for (j = 0; j < 3; ++j) {
579 pcache->vecs.params.interpolation_ranges[j].rmin =
580 base + (cie_cached_value)((double)temp[j].imin / factor);
581 pcache->vecs.params.interpolation_ranges[j].rmax =
582 base + (cie_cached_value)((double)temp[j].imax / factor);
583 if_debug3('c', "[c]interpolation_ranges[%d] = %g, %g\n", j,
584 cie_cached2float(pcache->vecs.params.interpolation_ranges[j].rmin),
585 cie_cached2float(pcache->vecs.params.interpolation_ranges[j].rmax));
586 }
587
588 }
589
590 /*
591 * Convert a scalar cache to a vector cache by multiplying the scalar
592 * values by a vector. Also set the range where interpolation is needed.
593 * This procedure is idempotent.
594 */
595 static void
cie_cache_mult(gx_cie_vector_cache * pcache,const gs_vector3 * pvec,const cie_cache_floats * pcf,double threshold)596 cie_cache_mult(gx_cie_vector_cache * pcache, const gs_vector3 * pvec,
597 const cie_cache_floats * pcf, double threshold)
598 {
599 float u = pvec->u, v = pvec->v, w = pvec->w;
600 int i;
601
602 pcache->vecs.params.base = float2cie_cached(pcf->params.base);
603 pcache->vecs.params.factor = float2cie_cached(pcf->params.factor);
604 pcache->vecs.params.limit =
605 float2cie_cached((gx_cie_cache_size - 1) / pcf->params.factor +
606 pcf->params.base);
607 for (i = 0; i < gx_cie_cache_size; ++i) {
608 float f = pcf->values[i];
609
610 pcache->vecs.values[i].u = float2cie_cached(f * u);
611 pcache->vecs.values[i].v = float2cie_cached(f * v);
612 pcache->vecs.values[i].w = float2cie_cached(f * w);
613 }
614 cie_cache_set_interpolation(pcache, threshold);
615 }
616
617 /*
618 * Set the interpolation ranges in a 3-vector cache, based on the ranges in
619 * the individual vector caches. This procedure is idempotent.
620 */
621 static void
cie_cache3_set_interpolation(gx_cie_vector_cache3_t * pvc)622 cie_cache3_set_interpolation(gx_cie_vector_cache3_t * pvc)
623 {
624 int j, k;
625
626 /* Iterate over output components. */
627 for (j = 0; j < 3; ++j) {
628 /* Iterate over sub-caches. */
629 cie_interpolation_range_t *p =
630 &pvc->caches[0].vecs.params.interpolation_ranges[j];
631 cie_cached_value rmin = p->rmin, rmax = p->rmax;
632
633 for (k = 1; k < 3; ++k) {
634 p = &pvc->caches[k].vecs.params.interpolation_ranges[j];
635 rmin = min(rmin, p->rmin), rmax = max(rmax, p->rmax);
636 }
637 pvc->interpolation_ranges[j].rmin = rmin;
638 pvc->interpolation_ranges[j].rmax = rmax;
639 if_debug3('c', "[c]Merged interpolation_ranges[%d] = %g, %g\n",
640 j, rmin, rmax);
641 }
642 }
643
644 /*
645 * Convert 3 scalar caches to vector caches by multiplying by a matrix.
646 * Return true iff the resulting cache is an identity transformation.
647 * This procedure is idempotent.
648 */
649 static bool
cie_cache_mult3(gx_cie_vector_cache3_t * pvc,const gs_matrix3 * pmat,double threshold)650 cie_cache_mult3(gx_cie_vector_cache3_t * pvc, const gs_matrix3 * pmat,
651 double threshold)
652 {
653 cie_cache_mult(&pvc->caches[0], &pmat->cu, &pvc->caches[0].floats, threshold);
654 cie_cache_mult(&pvc->caches[1], &pmat->cv, &pvc->caches[1].floats, threshold);
655 cie_cache_mult(&pvc->caches[2], &pmat->cw, &pvc->caches[2].floats, threshold);
656 cie_cache3_set_interpolation(pvc);
657 return pmat->is_identity & pvc->caches[0].floats.params.is_identity &
658 pvc->caches[1].floats.params.is_identity &
659 pvc->caches[2].floats.params.is_identity;
660 }
661
662 /* ------ Install a rendering dictionary ------ */
663
664 bool
vector_equal(const gs_vector3 * p1,const gs_vector3 * p2)665 vector_equal(const gs_vector3 *p1, const gs_vector3 *p2)
666 {
667 if (p1->u != p2->u)
668 return false;
669 if (p1->v != p2->v)
670 return false;
671 if (p1->w != p2->w)
672 return false;
673 return true;
674 }
675
676 bool
matrix_equal(const gs_matrix3 * p1,const gs_matrix3 * p2)677 matrix_equal(const gs_matrix3 *p1, const gs_matrix3 *p2)
678 {
679 if (p1->is_identity != p2->is_identity)
680 return false;
681 if (!vector_equal(&(p1->cu), &(p2->cu)))
682 return false;
683 if (!vector_equal(&(p1->cv), &(p2->cv)))
684 return false;
685 if (!vector_equal(&(p1->cw), &(p2->cw)))
686 return false;
687 return true;
688 }
689
690 static bool
transform_equal(const gs_cie_transform_proc3 * p1,const gs_cie_transform_proc3 * p2)691 transform_equal(const gs_cie_transform_proc3 *p1, const gs_cie_transform_proc3 *p2)
692 {
693 if (p1->proc != p2->proc)
694 return false;
695 if (p1->proc_data.size != p2->proc_data.size)
696 return false;
697 if (memcmp(p1->proc_data.data, p2->proc_data.data, p1->proc_data.size) != 0)
698 return false;
699 if (p1->driver_name != p2->driver_name)
700 return false;
701 if (p1->proc_name != p2->proc_name)
702 return false;
703 return true;
704 }
705
706 bool
range_equal(const gs_range3 * p1,const gs_range3 * p2)707 range_equal(const gs_range3 *p1, const gs_range3 *p2)
708 {
709 int k;
710
711 for (k = 0; k < 3; k++) {
712 if (p1->ranges[k].rmax != p2->ranges[k].rmax)
713 return false;
714 if (p1->ranges[k].rmin != p2->ranges[k].rmin)
715 return false;
716 }
717 return true;
718 }
719
720 /* setcolorrendering */
721 int
gs_setcolorrendering(gs_gstate * pgs,gs_cie_render * pcrd)722 gs_setcolorrendering(gs_gstate * pgs, gs_cie_render * pcrd)
723 {
724 int code = gs_cie_render_complete(pcrd);
725 const gs_cie_render *pcrd_old = pgs->cie_render;
726 bool joint_ok;
727
728 if (code < 0)
729 return code;
730 if (pcrd_old != 0 && pcrd->id == pcrd_old->id)
731 return 0; /* detect needless reselecting */
732 joint_ok =
733 pcrd_old != 0 &&
734 vector_equal(&pcrd->points.WhitePoint, &pcrd_old->points.WhitePoint) &&
735 vector_equal(&pcrd->points.BlackPoint, &pcrd_old->points.BlackPoint) &&
736 matrix_equal(&pcrd->MatrixPQR, &pcrd_old->MatrixPQR) &&
737 range_equal(&pcrd->RangePQR, &pcrd_old->RangePQR) &&
738 transform_equal(&pcrd->TransformPQR, &pcrd_old->TransformPQR);
739 rc_assign(pgs->cie_render, pcrd, "gs_setcolorrendering");
740 /* Initialize the joint caches if needed. */
741 if (!joint_ok)
742 code = gs_cie_cs_complete(pgs, true);
743 gx_unset_dev_color(pgs);
744 return code;
745 }
746
747 /* currentcolorrendering */
748 const gs_cie_render *
gs_currentcolorrendering(const gs_gstate * pgs)749 gs_currentcolorrendering(const gs_gstate * pgs)
750 {
751 return pgs->cie_render;
752 }
753
754 /* Unshare (allocating if necessary) the joint caches. */
755 gx_cie_joint_caches *
gx_unshare_cie_caches(gs_gstate * pgs)756 gx_unshare_cie_caches(gs_gstate * pgs)
757 {
758 gx_cie_joint_caches *pjc = pgs->cie_joint_caches;
759
760 rc_unshare_struct(pgs->cie_joint_caches, gx_cie_joint_caches,
761 &st_joint_caches, pgs->memory,
762 return 0, "gx_unshare_cie_caches");
763 if (pgs->cie_joint_caches != pjc) {
764 pjc = pgs->cie_joint_caches;
765 pjc->cspace_id = pjc->render_id = gs_no_id;
766 pjc->id_status = pjc->status = CIE_JC_STATUS_BUILT;
767 }
768 return pjc;
769 }
770
771 gx_cie_joint_caches *
gx_get_cie_caches_ref(gs_gstate * pgs,gs_memory_t * mem)772 gx_get_cie_caches_ref(gs_gstate * pgs, gs_memory_t * mem)
773 {
774 gx_cie_joint_caches *pjc = pgs->cie_joint_caches;
775
776 /* Take a reference here, to allow for the one that
777 * rc_unshare_struct might drop if it has to copy it.
778 * Whatever happens we will have taken 1 net new
779 * reference which we return to the caller. */
780 rc_increment(pgs->cie_joint_caches);
781 rc_unshare_struct(pjc, gx_cie_joint_caches,
782 &st_joint_caches, mem,
783 return NULL, "gx_unshare_cie_caches");
784 return pjc;
785 }
786
787 /* Compute the parameters for loading a cache, setting base and factor. */
788 /* This procedure is idempotent. */
789 void
gs_cie_cache_init(cie_cache_params * pcache,gs_sample_loop_params_t * pslp,const gs_range * domain,client_name_t cname)790 gs_cie_cache_init(cie_cache_params * pcache, gs_sample_loop_params_t * pslp,
791 const gs_range * domain, client_name_t cname)
792 {
793 /*
794 We need to map the values in the range [domain->rmin..domain->rmax].
795 However, if rmin < 0 < rmax and the function is non-linear, this can
796 lead to anomalies at zero, which is the default value for CIE colors.
797 The "correct" way to approach this is to run the mapping functions on
798 demand, but we don't want to deal with the complexities of the
799 callbacks this would involve (especially in the middle of rendering
800 images); instead, we adjust the range so that zero maps precisely to a
801 cache slot. Define:
802
803 A = domain->rmin;
804 B = domain->rmax;
805 N = gx_cie_cache_size - 1;
806
807 R = B - A;
808 h(v) = N * (v - A) / R; // the index of v in the cache
809 X = h(0).
810
811 If X is not an integer, we can decrease A and/increase B to make it
812 one. Let A' and B' be the adjusted values of A and B respectively,
813 and let K be the integer derived from X (either floor(X) or ceil(X)).
814 Define
815
816 f(K) = (K * B' + (N - K) * A') / N).
817
818 We want f(K) = 0. This occurs precisely when, for any real number
819 C != 0,
820
821 A' = -K * C;
822 B' = (N - K) * C.
823
824 In order to ensure A' <= A and B' >= B, we require
825
826 C >= -A / K;
827 C >= B / (N - K).
828
829 Since A' and B' must be exactly representable as floats, we round C
830 upward to ensure that it has no more than M mantissa bits, where
831
832 M = ARCH_FLOAT_MANTISSA_BITS - ceil(log2(N)).
833 */
834 float A = domain->rmin, B = domain->rmax;
835 double R = B - A, delta;
836 #define NN (gx_cie_cache_size - 1) /* 'N' is a member name, see end of proc */
837 #define N NN
838 #define CEIL_LOG2_N CIE_LOG2_CACHE_SIZE
839
840 /* Adjust the range if necessary. */
841 if (A < 0 && B >= 0) {
842 const double X = -N * A / R; /* know X > 0 */
843 /* Choose K to minimize range expansion. */
844 const int K = (int)(A + B < 0 ? floor(X) : ceil(X)); /* know 0 < K < N */
845 const double Ca = -A / K, Cb = B / (N - K); /* know Ca, Cb > 0 */
846 double C = max(Ca, Cb); /* know C > 0 */
847 const int M = ARCH_FLOAT_MANTISSA_BITS - CEIL_LOG2_N;
848 int cexp;
849 const double cfrac = frexp(C, &cexp);
850
851 if_debug4('c', "[c]adjusting cache_init(%8g, %8g), X = %8g, K = %d:\n",
852 A, B, X, K);
853 /* Round C to no more than M significant bits. See above. */
854 C = ldexp(ceil(ldexp(cfrac, M)), cexp - M);
855 /* Finally, compute A' and B'. */
856 A = -K * C;
857 B = (N - K) * C;
858 if_debug2('c', "[c] => %8g, %8g\n", A, B);
859 R = B - A;
860 }
861 delta = R / N;
862 #ifdef CIE_CACHE_INTERPOLATE
863 pcache->base = A; /* no rounding */
864 #else
865 pcache->base = A - delta / 2; /* so lookup will round */
866 #endif
867 /*
868 * If size of the domain is zero, then use 1.0 as the scaling
869 * factor. This prevents divide by zero errors in later calculations.
870 * This should only occurs with zero matrices. It does occur with
871 * Genoa test file 050-01.ps.
872 */
873 pcache->factor = (any_abs(delta) < 1e-30 ? 1.0 : N / R);
874 if_debug4('c', "[c]cache %s 0x%lx base=%g, factor=%g\n",
875 (const char *)cname, (ulong) pcache,
876 pcache->base, pcache->factor);
877 pslp->A = A;
878 pslp->B = B;
879 #undef N
880 pslp->N = NN;
881 #undef NN
882 }
883
884 /* ------ Complete a rendering structure ------ */
885
886 /*
887 * Compute the derived values in a CRD that don't involve the cached
888 * procedure values. This procedure is idempotent.
889 */
890 static void cie_transform_range3(const gs_range3 *, const gs_matrix3 *,
891 gs_range3 *);
892 int
gs_cie_render_init(gs_cie_render * pcrd)893 gs_cie_render_init(gs_cie_render * pcrd)
894 {
895 gs_matrix3 PQR_inverse;
896
897 if (pcrd->status >= CIE_RENDER_STATUS_INITED)
898 return 0; /* init already done */
899 if_debug_matrix3("[c]CRD MatrixLMN =", &pcrd->MatrixLMN);
900 cie_matrix_init(&pcrd->MatrixLMN);
901 if_debug_matrix3("[c]CRD MatrixABC =", &pcrd->MatrixABC);
902 cie_matrix_init(&pcrd->MatrixABC);
903 if_debug_matrix3("[c]CRD MatrixPQR =", &pcrd->MatrixPQR);
904 cie_matrix_init(&pcrd->MatrixPQR);
905 cie_invert3(&pcrd->MatrixPQR, &PQR_inverse);
906 cie_matrix_mult3(&pcrd->MatrixLMN, &PQR_inverse,
907 &pcrd->MatrixPQR_inverse_LMN);
908 cie_transform_range3(&pcrd->RangePQR, &pcrd->MatrixPQR_inverse_LMN,
909 &pcrd->DomainLMN);
910 cie_transform_range3(&pcrd->RangeLMN, &pcrd->MatrixABC,
911 &pcrd->DomainABC);
912 cie_mult3(&pcrd->points.WhitePoint, &pcrd->MatrixPQR, &pcrd->wdpqr);
913 cie_mult3(&pcrd->points.BlackPoint, &pcrd->MatrixPQR, &pcrd->bdpqr);
914 pcrd->status = CIE_RENDER_STATUS_INITED;
915 return 0;
916 }
917
918 /*
919 * Sample the EncodeLMN, EncodeABC, and RenderTableT CRD procedures, and
920 * load the caches. This procedure is idempotent.
921 */
922 int
gs_cie_render_sample(gs_cie_render * pcrd)923 gs_cie_render_sample(gs_cie_render * pcrd)
924 {
925 int code;
926
927 if (pcrd->status >= CIE_RENDER_STATUS_SAMPLED)
928 return 0; /* sampling already done */
929 code = gs_cie_render_init(pcrd);
930 if (code < 0)
931 return code;
932 CIE_LOAD_CACHE_BODY(pcrd->caches.EncodeLMN.caches, pcrd->DomainLMN.ranges,
933 &pcrd->EncodeLMN, Encode_default, pcrd, "EncodeLMN");
934 cache3_set_linear(&pcrd->caches.EncodeLMN);
935 CIE_LOAD_CACHE_BODY(pcrd->caches.EncodeABC, pcrd->DomainABC.ranges,
936 &pcrd->EncodeABC, Encode_default, pcrd, "EncodeABC");
937 if (pcrd->RenderTable.lookup.table != 0) {
938 int i, j, m = pcrd->RenderTable.lookup.m;
939 gs_sample_loop_params_t lp;
940 bool is_identity = true;
941
942 for (j = 0; j < m; j++) {
943 gs_cie_cache_init(&pcrd->caches.RenderTableT[j].fracs.params,
944 &lp, &Range3_default.ranges[0],
945 "RenderTableT");
946 is_identity &= pcrd->RenderTable.T.procs[j] ==
947 RenderTableT_default.procs[j];
948 }
949 pcrd->caches.RenderTableT_is_identity = is_identity;
950 /*
951 * Unfortunately, we defined the first argument of the RenderTable
952 * T procedures as being a byte, limiting the number of distinct
953 * cache entries to 256 rather than gx_cie_cache_size.
954 * We confine this decision to this loop, rather than propagating
955 * it to the procedures that use the cached data, so that we can
956 * change it more easily at some future time.
957 */
958 for (i = 0; i < gx_cie_cache_size; i++) {
959 #if gx_cie_log2_cache_size >= 8
960 byte value = i >> (gx_cie_log2_cache_size - 8);
961 #else
962 byte value = (i << (8 - gx_cie_log2_cache_size)) +
963 (i >> (gx_cie_log2_cache_size * 2 - 8));
964 #endif
965 for (j = 0; j < m; j++) {
966 pcrd->caches.RenderTableT[j].fracs.values[i] =
967 (*pcrd->RenderTable.T.procs[j])(value, pcrd);
968 if_debug3('C', "[C]RenderTableT[%d,%d] = %g\n",
969 i, j,
970 frac2float(pcrd->caches.RenderTableT[j].fracs.values[i]));
971 }
972 }
973 }
974 pcrd->status = CIE_RENDER_STATUS_SAMPLED;
975 return 0;
976 }
977
978 /* Transform a set of ranges. */
979 static void
cie_transform_range(const gs_range3 * in,double mu,double mv,double mw,gs_range * out)980 cie_transform_range(const gs_range3 * in, double mu, double mv, double mw,
981 gs_range * out)
982 {
983 float umin = mu * in->ranges[0].rmin, umax = mu * in->ranges[0].rmax;
984 float vmin = mv * in->ranges[1].rmin, vmax = mv * in->ranges[1].rmax;
985 float wmin = mw * in->ranges[2].rmin, wmax = mw * in->ranges[2].rmax;
986 float temp;
987
988 if (umin > umax)
989 temp = umin, umin = umax, umax = temp;
990 if (vmin > vmax)
991 temp = vmin, vmin = vmax, vmax = temp;
992 if (wmin > wmax)
993 temp = wmin, wmin = wmax, wmax = temp;
994 out->rmin = umin + vmin + wmin;
995 out->rmax = umax + vmax + wmax;
996 }
997 static void
cie_transform_range3(const gs_range3 * in,const gs_matrix3 * mat,gs_range3 * out)998 cie_transform_range3(const gs_range3 * in, const gs_matrix3 * mat,
999 gs_range3 * out)
1000 {
1001 cie_transform_range(in, mat->cu.u, mat->cv.u, mat->cw.u,
1002 &out->ranges[0]);
1003 cie_transform_range(in, mat->cu.v, mat->cv.v, mat->cw.v,
1004 &out->ranges[1]);
1005 cie_transform_range(in, mat->cu.w, mat->cv.w, mat->cw.w,
1006 &out->ranges[2]);
1007 }
1008
1009 /*
1010 * Finish preparing a CRD for installation, by restricting and/or
1011 * transforming the cached procedure values.
1012 * This procedure is idempotent.
1013 */
1014 int
gs_cie_render_complete(gs_cie_render * pcrd)1015 gs_cie_render_complete(gs_cie_render * pcrd)
1016 {
1017 int code;
1018
1019 if (pcrd->status >= CIE_RENDER_STATUS_COMPLETED)
1020 return 0; /* completion already done */
1021 code = gs_cie_render_sample(pcrd);
1022 if (code < 0)
1023 return code;
1024 /*
1025 * Since range restriction happens immediately after
1026 * the cache lookup, we can save a step by restricting
1027 * the values in the cache entries.
1028 *
1029 * If there is no lookup table, we want the final ABC values
1030 * to be fracs; if there is a table, we want them to be
1031 * appropriately scaled ints.
1032 */
1033 pcrd->MatrixABCEncode = pcrd->MatrixABC;
1034 {
1035 int c;
1036 double f;
1037
1038 for (c = 0; c < 3; c++) {
1039 gx_cie_float_fixed_cache *pcache = &pcrd->caches.EncodeABC[c];
1040
1041 cie_cache_restrict(&pcrd->caches.EncodeLMN.caches[c].floats,
1042 &pcrd->RangeLMN.ranges[c]);
1043 cie_cache_restrict(&pcrd->caches.EncodeABC[c].floats,
1044 &pcrd->RangeABC.ranges[c]);
1045 if (pcrd->RenderTable.lookup.table == 0) {
1046 cie_cache_restrict(&pcache->floats,
1047 &Range3_default.ranges[0]);
1048 gs_cie_cache_to_fracs(&pcache->floats, &pcache->fixeds.fracs);
1049 pcache->fixeds.fracs.params.is_identity = false;
1050 } else {
1051 int i;
1052 int n = pcrd->RenderTable.lookup.dims[c];
1053
1054 #ifdef CIE_RENDER_TABLE_INTERPOLATE
1055 # define SCALED_INDEX(f, n, itemp)\
1056 RESTRICTED_INDEX(f * (1 << _cie_interpolate_bits),\
1057 (n) << _cie_interpolate_bits, itemp)
1058 #else
1059 int m = pcrd->RenderTable.lookup.m;
1060 int k =
1061 (c == 0 ? 1 : c == 1 ?
1062 m * pcrd->RenderTable.lookup.dims[2] : m);
1063 # define SCALED_INDEX(f, n, itemp)\
1064 (RESTRICTED_INDEX(f, n, itemp) * k)
1065 #endif
1066 const gs_range *prange = pcrd->RangeABC.ranges + c;
1067 double scale = (n - 1) / (prange->rmax - prange->rmin);
1068
1069 for (i = 0; i < gx_cie_cache_size; ++i) {
1070 float v =
1071 (pcache->floats.values[i] - prange->rmin) * scale
1072 #ifndef CIE_RENDER_TABLE_INTERPOLATE
1073 + 0.5
1074 #endif
1075 ;
1076 int itemp;
1077
1078 if_debug5('c',
1079 "[c]cache[%d][%d] = %g => %g => %d\n",
1080 c, i, pcache->floats.values[i], v,
1081 SCALED_INDEX(v, n, itemp));
1082 pcache->fixeds.ints.values[i] =
1083 SCALED_INDEX(v, n, itemp);
1084 }
1085 pcache->fixeds.ints.params = pcache->floats.params;
1086 pcache->fixeds.ints.params.is_identity = false;
1087 #undef SCALED_INDEX
1088 }
1089 }
1090 /* Fold the scaling of the EncodeABC cache index */
1091 /* into MatrixABC. */
1092 #define MABC(i, t)\
1093 f = pcrd->caches.EncodeABC[i].floats.params.factor;\
1094 pcrd->MatrixABCEncode.cu.t *= f;\
1095 pcrd->MatrixABCEncode.cv.t *= f;\
1096 pcrd->MatrixABCEncode.cw.t *= f;\
1097 pcrd->EncodeABC_base[i] =\
1098 float2cie_cached(pcrd->caches.EncodeABC[i].floats.params.base * f)
1099 MABC(0, u);
1100 MABC(1, v);
1101 MABC(2, w);
1102 #undef MABC
1103 pcrd->MatrixABCEncode.is_identity = 0;
1104 }
1105 cie_cache_mult3(&pcrd->caches.EncodeLMN, &pcrd->MatrixABCEncode,
1106 CACHE_THRESHOLD);
1107 pcrd->status = CIE_RENDER_STATUS_COMPLETED;
1108 return 0;
1109 }
1110
1111 /* Apply a range restriction to a cache. */
1112 static void
cie_cache_restrict(cie_cache_floats * pcache,const gs_range * prange)1113 cie_cache_restrict(cie_cache_floats * pcache, const gs_range * prange)
1114 {
1115 int i;
1116
1117 for (i = 0; i < gx_cie_cache_size; i++) {
1118 float v = pcache->values[i];
1119
1120 if (v < prange->rmin)
1121 pcache->values[i] = prange->rmin;
1122 else if (v > prange->rmax)
1123 pcache->values[i] = prange->rmax;
1124 }
1125 }
1126
1127 /* Convert a cache from floats to fracs. */
1128 /* Note that the two may be aliased. */
1129 void
gs_cie_cache_to_fracs(const cie_cache_floats * pfloats,cie_cache_fracs * pfracs)1130 gs_cie_cache_to_fracs(const cie_cache_floats *pfloats, cie_cache_fracs *pfracs)
1131 {
1132 int i;
1133
1134 /* Loop from bottom to top so that we don't */
1135 /* overwrite elements before they're used. */
1136 for (i = 0; i < gx_cie_cache_size; ++i)
1137 pfracs->values[i] = float2frac(pfloats->values[i]);
1138 pfracs->params = pfloats->params;
1139 }
1140
1141 /* ------ Fill in the joint cache ------ */
1142
1143 /* If the current color space is a CIE space, or has a CIE base space, */
1144 /* return a pointer to the common part of the space; otherwise return 0. */
1145 static const gs_cie_common *
cie_cs_common_abc(const gs_color_space * pcs_orig,const gs_cie_abc ** ppabc)1146 cie_cs_common_abc(const gs_color_space *pcs_orig, const gs_cie_abc **ppabc)
1147 {
1148 const gs_color_space *pcs = pcs_orig;
1149
1150 *ppabc = 0;
1151 do {
1152 switch (pcs->type->index) {
1153 case gs_color_space_index_CIEDEF:
1154 *ppabc = (const gs_cie_abc *)pcs->params.def;
1155 return &pcs->params.def->common;
1156 case gs_color_space_index_CIEDEFG:
1157 *ppabc = (const gs_cie_abc *)pcs->params.defg;
1158 return &pcs->params.defg->common;
1159 case gs_color_space_index_CIEABC:
1160 *ppabc = pcs->params.abc;
1161 return &pcs->params.abc->common;
1162 case gs_color_space_index_CIEA:
1163 return &pcs->params.a->common;
1164 default:
1165 pcs = gs_cspace_base_space(pcs);
1166 break;
1167 }
1168 } while (pcs != 0);
1169
1170 return 0;
1171 }
1172 const gs_cie_common *
gs_cie_cs_common(const gs_gstate * pgs)1173 gs_cie_cs_common(const gs_gstate * pgs)
1174 {
1175 const gs_cie_abc *ignore_pabc;
1176
1177 return cie_cs_common_abc(gs_currentcolorspace_inline(pgs), &ignore_pabc);
1178 }
1179
1180 /*
1181 * Mark the joint caches as needing completion. This is done lazily,
1182 * when a color is being mapped. However, make sure the joint caches
1183 * exist now.
1184 */
1185 int
gs_cie_cs_complete(gs_gstate * pgs,bool init)1186 gs_cie_cs_complete(gs_gstate * pgs, bool init)
1187 {
1188 gx_cie_joint_caches *pjc = gx_unshare_cie_caches(pgs);
1189
1190 if (pjc == 0)
1191 return_error(gs_error_VMerror);
1192 pjc->status = (init ? CIE_JC_STATUS_BUILT : CIE_JC_STATUS_INITED);
1193 return 0;
1194 }
1195 /* Actually complete the joint caches. */
1196 int
gs_cie_jc_complete(const gs_gstate * pgs,const gs_color_space * pcs)1197 gs_cie_jc_complete(const gs_gstate *pgs, const gs_color_space *pcs)
1198 {
1199 const gs_cie_abc *pabc;
1200 const gs_cie_common *common = cie_cs_common_abc(pcs, &pabc);
1201 gs_cie_render *pcrd = pgs->cie_render;
1202 gx_cie_joint_caches *pjc = pgs->cie_joint_caches;
1203
1204 if (pjc->cspace_id == pcs->id &&
1205 pjc->render_id == pcrd->id)
1206 pjc->status = pjc->id_status;
1207 switch (pjc->status) {
1208 case CIE_JC_STATUS_BUILT: {
1209 int code = cie_joint_caches_init(pjc, common, pcrd);
1210
1211 if (code < 0)
1212 return code;
1213 }
1214 /* falls through */
1215 case CIE_JC_STATUS_INITED:
1216 cie_joint_caches_complete(pjc, common, pabc, pcrd);
1217 pjc->cspace_id = pcs->id;
1218 pjc->render_id = pcrd->id;
1219 pjc->id_status = pjc->status = CIE_JC_STATUS_COMPLETED;
1220 /* falls through */
1221 case CIE_JC_STATUS_COMPLETED:
1222 break;
1223 }
1224 return 0;
1225 }
1226
1227 /*
1228 * Compute the source and destination WhitePoint and BlackPoint for
1229 * the TransformPQR procedure.
1230 */
1231 int
gs_cie_compute_points_sd(gx_cie_joint_caches * pjc,const gs_cie_common * pcie,const gs_cie_render * pcrd)1232 gs_cie_compute_points_sd(gx_cie_joint_caches *pjc,
1233 const gs_cie_common * pcie,
1234 const gs_cie_render * pcrd)
1235 {
1236 gs_cie_wbsd *pwbsd = &pjc->points_sd;
1237
1238 pwbsd->ws.xyz = pcie->points.WhitePoint;
1239 cie_mult3(&pwbsd->ws.xyz, &pcrd->MatrixPQR, &pwbsd->ws.pqr);
1240 pwbsd->bs.xyz = pcie->points.BlackPoint;
1241 cie_mult3(&pwbsd->bs.xyz, &pcrd->MatrixPQR, &pwbsd->bs.pqr);
1242 pwbsd->wd.xyz = pcrd->points.WhitePoint;
1243 pwbsd->wd.pqr = pcrd->wdpqr;
1244 pwbsd->bd.xyz = pcrd->points.BlackPoint;
1245 pwbsd->bd.pqr = pcrd->bdpqr;
1246 return 0;
1247 }
1248
1249 /*
1250 * Sample the TransformPQR procedure for the joint caches.
1251 * This routine is idempotent.
1252 */
1253 static int
cie_joint_caches_init(gx_cie_joint_caches * pjc,const gs_cie_common * pcie,gs_cie_render * pcrd)1254 cie_joint_caches_init(gx_cie_joint_caches * pjc,
1255 const gs_cie_common * pcie,
1256 gs_cie_render * pcrd)
1257 {
1258 bool is_identity;
1259 int j;
1260
1261 gs_cie_compute_points_sd(pjc, pcie, pcrd);
1262 /*
1263 * If a client pre-loaded the cache, we can't adjust the range.
1264 * ****** WRONG ******
1265 */
1266 if (pcrd->TransformPQR.proc == TransformPQR_from_cache.proc)
1267 return 0;
1268 is_identity = pcrd->TransformPQR.proc == TransformPQR_default.proc;
1269 for (j = 0; j < 3; j++) {
1270 int i;
1271 gs_sample_loop_params_t lp;
1272
1273 gs_cie_cache_init(&pjc->TransformPQR.caches[j].floats.params, &lp,
1274 &pcrd->RangePQR.ranges[j], "TransformPQR");
1275 for (i = 0; i <= lp.N; ++i) {
1276 float in = SAMPLE_LOOP_VALUE(i, lp);
1277 float out;
1278 int code = (*pcrd->TransformPQR.proc)(j, in, &pjc->points_sd,
1279 pcrd, &out);
1280
1281 if (code < 0)
1282 return code;
1283 pjc->TransformPQR.caches[j].floats.values[i] = out;
1284 if_debug4('C', "[C]TransformPQR[%d,%d] = %g => %g\n",
1285 j, i, in, out);
1286 }
1287 pjc->TransformPQR.caches[j].floats.params.is_identity = is_identity;
1288 }
1289 return 0;
1290 }
1291
1292 /*
1293 * Complete the loading of the joint caches.
1294 * This routine is idempotent.
1295 */
1296 static void
cie_joint_caches_complete(gx_cie_joint_caches * pjc,const gs_cie_common * pcie,const gs_cie_abc * pabc,const gs_cie_render * pcrd)1297 cie_joint_caches_complete(gx_cie_joint_caches * pjc,
1298 const gs_cie_common * pcie,
1299 const gs_cie_abc * pabc /* NULL if CIEA */,
1300 const gs_cie_render * pcrd)
1301 {
1302 gs_matrix3 mat3, mat2;
1303 gs_matrix3 MatrixLMN_PQR;
1304 int j;
1305
1306 pjc->remap_finish = gx_cie_real_remap_finish;
1307
1308 /*
1309 * We number the pipeline steps as follows:
1310 * 1 - DecodeABC/MatrixABC
1311 * 2 - DecodeLMN/MatrixLMN/MatrixPQR
1312 * 3 - TransformPQR/MatrixPQR'/MatrixLMN
1313 * 4 - EncodeLMN/MatrixABC
1314 * 5 - EncodeABC, RenderTable (we don't do anything with this here)
1315 * We work from back to front, combining steps where possible.
1316 * Currently we only combine steps if a procedure is the identity
1317 * transform, but we could do it whenever the procedure is linear.
1318 * A project for another day....
1319 */
1320
1321 /* Step 4 */
1322
1323 #ifdef OPTIMIZE_CIE_MAPPING
1324 if (pcrd->caches.EncodeLMN.caches[0].floats.params.is_identity &&
1325 pcrd->caches.EncodeLMN.caches[1].floats.params.is_identity &&
1326 pcrd->caches.EncodeLMN.caches[2].floats.params.is_identity
1327 ) {
1328 /* Fold step 4 into step 3. */
1329 if_debug0('c', "[c]EncodeLMN is identity, folding MatrixABC(Encode) into MatrixPQR'+LMN.\n");
1330 cie_matrix_mult3(&pcrd->MatrixABCEncode, &pcrd->MatrixPQR_inverse_LMN,
1331 &mat3);
1332 pjc->skipEncodeLMN = true;
1333 } else
1334 #endif /* OPTIMIZE_CIE_MAPPING */
1335 {
1336 if_debug0('c', "[c]EncodeLMN is not identity.\n");
1337 mat3 = pcrd->MatrixPQR_inverse_LMN;
1338 pjc->skipEncodeLMN = false;
1339 }
1340
1341 /* Step 3 */
1342
1343 cache3_set_linear(&pjc->TransformPQR);
1344 cie_matrix_mult3(&pcrd->MatrixPQR, &pcie->MatrixLMN,
1345 &MatrixLMN_PQR);
1346
1347 #ifdef OPTIMIZE_CIE_MAPPING
1348 if (pjc->TransformPQR.caches[0].floats.params.is_identity &
1349 pjc->TransformPQR.caches[1].floats.params.is_identity &
1350 pjc->TransformPQR.caches[2].floats.params.is_identity
1351 ) {
1352 /* Fold step 3 into step 2. */
1353 if_debug0('c', "[c]TransformPQR is identity, folding MatrixPQR'+LMN into MatrixLMN+PQR.\n");
1354 cie_matrix_mult3(&mat3, &MatrixLMN_PQR, &mat2);
1355 pjc->skipPQR = true;
1356 } else
1357 #endif /* OPTIMIZE_CIE_MAPPING */
1358 {
1359 if_debug0('c', "[c]TransformPQR is not identity.\n");
1360 mat2 = MatrixLMN_PQR;
1361 for (j = 0; j < 3; j++) {
1362 cie_cache_restrict(&pjc->TransformPQR.caches[j].floats,
1363 &pcrd->RangePQR.ranges[j]);
1364 }
1365 cie_cache_mult3(&pjc->TransformPQR, &mat3, CACHE_THRESHOLD);
1366 pjc->skipPQR = false;
1367 }
1368
1369 /* Steps 2 & 1 */
1370
1371 #ifdef OPTIMIZE_CIE_MAPPING
1372 if (pcie->caches.DecodeLMN[0].floats.params.is_identity &
1373 pcie->caches.DecodeLMN[1].floats.params.is_identity &
1374 pcie->caches.DecodeLMN[2].floats.params.is_identity
1375 ) {
1376 if_debug0('c', "[c]DecodeLMN is identity, folding MatrixLMN+PQR into MatrixABC.\n");
1377 if (!pabc) {
1378 pjc->skipDecodeLMN = mat2.is_identity;
1379 pjc->skipDecodeABC = false;
1380 if (!pjc->skipDecodeLMN) {
1381 for (j = 0; j < 3; j++) {
1382 cie_cache_mult(&pjc->DecodeLMN.caches[j], &mat2.cu + j,
1383 &pcie->caches.DecodeLMN[j].floats,
1384 CACHE_THRESHOLD);
1385 }
1386 cie_cache3_set_interpolation(&pjc->DecodeLMN);
1387 }
1388 } else {
1389 /*
1390 * Fold step 2 into step 1. This is a little different because
1391 * the data for step 1 are in the color space structure.
1392 */
1393 gs_matrix3 mat1;
1394
1395 cie_matrix_mult3(&mat2, &pabc->MatrixABC, &mat1);
1396 for (j = 0; j < 3; j++) {
1397 cie_cache_mult(&pjc->DecodeLMN.caches[j], &mat1.cu + j,
1398 &pabc->caches.DecodeABC.caches[j].floats,
1399 CACHE_THRESHOLD);
1400 }
1401 cie_cache3_set_interpolation(&pjc->DecodeLMN);
1402 pjc->skipDecodeLMN = false;
1403 pjc->skipDecodeABC = true;
1404 }
1405 } else
1406 #endif /* OPTIMIZE_CIE_MAPPING */
1407 {
1408 if_debug0('c', "[c]DecodeLMN is not identity.\n");
1409 for (j = 0; j < 3; j++) {
1410 cie_cache_mult(&pjc->DecodeLMN.caches[j], &mat2.cu + j,
1411 &pcie->caches.DecodeLMN[j].floats,
1412 CACHE_THRESHOLD);
1413 }
1414 cie_cache3_set_interpolation(&pjc->DecodeLMN);
1415 pjc->skipDecodeLMN = false;
1416 pjc->skipDecodeABC = pabc != 0 && pabc->caches.skipABC;
1417 }
1418
1419 }
1420
1421 /*
1422 * Initialize (just enough of) an gs_gstate so that "concretizing" colors
1423 * using this gs_gstate will do only the CIE->XYZ mapping. This is a
1424 * semi-hack for the PDF writer.
1425 */
1426 int
gx_cie_to_xyz_alloc(gs_gstate ** ppgs,const gs_color_space * pcs,gs_memory_t * mem)1427 gx_cie_to_xyz_alloc(gs_gstate **ppgs, const gs_color_space *pcs,
1428 gs_memory_t *mem)
1429 {
1430 /*
1431 * In addition to the gs_gstate itself, we need the joint caches.
1432 */
1433 gs_gstate *pgs =
1434 gs_alloc_struct(mem, gs_gstate, &st_gs_gstate,
1435 "gx_cie_to_xyz_alloc(gs_gstate)");
1436 gx_cie_joint_caches *pjc;
1437 const gs_cie_abc *pabc;
1438 const gs_cie_common *pcie = cie_cs_common_abc(pcs, &pabc);
1439 int j;
1440
1441 if (pgs == 0)
1442 return_error(gs_error_VMerror);
1443 memset(pgs, 0, sizeof(*pgs)); /* mostly paranoia */
1444 pgs->memory = mem;
1445 GS_STATE_INIT_VALUES(pgs, 1.0);
1446 gs_gstate_initialize(pgs, mem);
1447
1448 pjc = gs_alloc_struct(mem, gx_cie_joint_caches, &st_joint_caches,
1449 "gx_cie_to_xyz_free(joint caches)");
1450 if (pjc == 0) {
1451 gs_free_object(mem, pgs, "gx_cie_to_xyz_alloc(gs_gstate)");
1452 return_error(gs_error_VMerror);
1453 }
1454 rc_init(pjc, mem, 1);
1455
1456 /*
1457 * Perform an abbreviated version of cie_joint_caches_complete.
1458 * Don't bother with any optimizations.
1459 */
1460 for (j = 0; j < 3; j++) {
1461 cie_cache_mult(&pjc->DecodeLMN.caches[j], &pcie->MatrixLMN.cu + j,
1462 &pcie->caches.DecodeLMN[j].floats,
1463 CACHE_THRESHOLD);
1464 }
1465 cie_cache3_set_interpolation(&pjc->DecodeLMN);
1466 pjc->skipDecodeLMN = false;
1467 pjc->skipDecodeABC = pabc != 0 && pabc->caches.skipABC;
1468 /* Mark the joint caches as completed. */
1469 pjc->remap_finish = gx_cie_xyz_remap_finish;
1470 pjc->cspace_id = pcs->id;
1471 pjc->status = CIE_JC_STATUS_COMPLETED;
1472 pgs->cie_joint_caches = pjc;
1473 pgs->cie_to_xyz = true;
1474 *ppgs = pgs;
1475 return 0;
1476 }
1477 void
gx_cie_to_xyz_free(gs_gstate * pgs)1478 gx_cie_to_xyz_free(gs_gstate *pgs)
1479 {
1480 gs_memory_t *mem = pgs->memory;
1481
1482 rc_decrement(pgs->cie_joint_caches,"gx_cie_to_xyz_free");
1483
1484 /* Free up the ICC objects if created */ /* FIXME: does this need to be thread safe */
1485 if (pgs->icc_link_cache != NULL) {
1486 rc_decrement(pgs->icc_link_cache,"gx_cie_to_xyz_free");
1487 }
1488 if (pgs->icc_manager != NULL) {
1489 rc_decrement(pgs->icc_manager,"gx_cie_to_xyz_free");
1490 }
1491 if (pgs->icc_profile_cache != NULL) {
1492 rc_decrement(pgs->icc_profile_cache,"gx_cie_to_xyz_free");
1493 }
1494 gs_free_object(mem, pgs, "gx_cie_to_xyz_free(gs_gstate)");
1495 }
1496
1497 /* ================ Utilities ================ */
1498
1499 /* Multiply a vector by a matrix. */
1500 /* Note that we are computing M * V where v is a column vector. */
1501 void
cie_mult3(const gs_vector3 * in,register const gs_matrix3 * mat,gs_vector3 * out)1502 cie_mult3(const gs_vector3 * in, register const gs_matrix3 * mat,
1503 gs_vector3 * out)
1504 {
1505 if_debug_vector3("[c]mult", in);
1506 if_debug_matrix3(" *", mat);
1507 {
1508 float u = in->u, v = in->v, w = in->w;
1509
1510 out->u = (u * mat->cu.u) + (v * mat->cv.u) + (w * mat->cw.u);
1511 out->v = (u * mat->cu.v) + (v * mat->cv.v) + (w * mat->cw.v);
1512 out->w = (u * mat->cu.w) + (v * mat->cv.w) + (w * mat->cw.w);
1513 }
1514 if_debug_vector3(" =", out);
1515 }
1516
1517 /*
1518 * Multiply two matrices. Note that the composition of the transformations
1519 * M1 followed by M2 is M2 * M1, not M1 * M2. (See gscie.h for details.)
1520 */
1521 void
cie_matrix_mult3(const gs_matrix3 * ma,const gs_matrix3 * mb,gs_matrix3 * mc)1522 cie_matrix_mult3(const gs_matrix3 *ma, const gs_matrix3 *mb, gs_matrix3 *mc)
1523 {
1524 gs_matrix3 mprod;
1525 gs_matrix3 *mp = (mc == ma || mc == mb ? &mprod : mc);
1526
1527 if_debug_matrix3("[c]matrix_mult", ma);
1528 if_debug_matrix3(" *", mb);
1529 cie_mult3(&mb->cu, ma, &mp->cu);
1530 cie_mult3(&mb->cv, ma, &mp->cv);
1531 cie_mult3(&mb->cw, ma, &mp->cw);
1532 cie_matrix_init(mp);
1533 if_debug_matrix3(" =", mp);
1534 if (mp != mc)
1535 *mc = *mp;
1536 }
1537
1538 /*
1539 * Transpose a 3x3 matrix. In and out can not be the same
1540 */
1541 void
cie_matrix_transpose3(const gs_matrix3 * in,gs_matrix3 * out)1542 cie_matrix_transpose3(const gs_matrix3 *in, gs_matrix3 *out)
1543 {
1544 out->cu.u = in->cu.u;
1545 out->cu.v = in->cv.u;
1546 out->cu.w = in->cw.u;
1547
1548 out->cv.u = in->cu.v;
1549 out->cv.v = in->cv.v;
1550 out->cv.w = in->cw.v;
1551
1552 out->cw.u = in->cu.w;
1553 out->cw.v = in->cv.w;
1554 out->cw.w = in->cw.w;
1555 }
1556
1557 /* Invert a matrix. */
1558 /* The output must not be an alias for the input. */
1559 static void
cie_invert3(const gs_matrix3 * in,gs_matrix3 * out)1560 cie_invert3(const gs_matrix3 *in, gs_matrix3 *out)
1561 { /* This is a brute force algorithm; maybe there are better. */
1562 /* We label the array elements */
1563 /* [ A B C ] */
1564 /* [ D E F ] */
1565 /* [ G H I ] */
1566 #define A cu.u
1567 #define B cv.u
1568 #define C cw.u
1569 #define D cu.v
1570 #define E cv.v
1571 #define F cw.v
1572 #define G cu.w
1573 #define H cv.w
1574 #define I cw.w
1575 double coA = in->E * in->I - in->F * in->H;
1576 double coB = in->F * in->G - in->D * in->I;
1577 double coC = in->D * in->H - in->E * in->G;
1578 double det = in->A * coA + in->B * coB + in->C * coC;
1579
1580 if_debug_matrix3("[c]invert", in);
1581 out->A = coA / det;
1582 out->D = coB / det;
1583 out->G = coC / det;
1584 out->B = (in->C * in->H - in->B * in->I) / det;
1585 out->E = (in->A * in->I - in->C * in->G) / det;
1586 out->H = (in->B * in->G - in->A * in->H) / det;
1587 out->C = (in->B * in->F - in->C * in->E) / det;
1588 out->F = (in->C * in->D - in->A * in->F) / det;
1589 out->I = (in->A * in->E - in->B * in->D) / det;
1590 if_debug_matrix3(" =", out);
1591 #undef A
1592 #undef B
1593 #undef C
1594 #undef D
1595 #undef E
1596 #undef F
1597 #undef G
1598 #undef H
1599 #undef I
1600 out->is_identity = in->is_identity;
1601 }
1602
1603 /* Set the is_identity flag that accelerates multiplication. */
1604 static void
cie_matrix_init(register gs_matrix3 * mat)1605 cie_matrix_init(register gs_matrix3 * mat)
1606 {
1607 mat->is_identity =
1608 mat->cu.u == 1.0 && is_fzero2(mat->cu.v, mat->cu.w) &&
1609 mat->cv.v == 1.0 && is_fzero2(mat->cv.u, mat->cv.w) &&
1610 mat->cw.w == 1.0 && is_fzero2(mat->cw.u, mat->cw.v);
1611 }
1612
1613 bool
gx_color_space_needs_cie_caches(const gs_color_space * pcs)1614 gx_color_space_needs_cie_caches(const gs_color_space * pcs)
1615 {
1616 switch (pcs->type->index) {
1617 case gs_color_space_index_CIEDEFG:
1618 case gs_color_space_index_CIEDEF:
1619 case gs_color_space_index_CIEABC:
1620 case gs_color_space_index_CIEA:
1621 return true;
1622 case gs_color_space_index_ICC:
1623 return false;
1624 case gs_color_space_index_DevicePixel:
1625 case gs_color_space_index_DeviceN:
1626 case gs_color_space_index_Separation:
1627 case gs_color_space_index_Indexed:
1628 case gs_color_space_index_Pattern:
1629 return gx_color_space_needs_cie_caches(pcs->base_space);
1630 default:
1631 return false;
1632 }
1633 }
1634