1 /* Copyright (C) 2001-2006 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, modified
8    or distributed except as expressly authorized under the terms of that
9    license.  Refer to licensing information at http://www.artifex.com/
10    or contact Artifex Software, Inc.,  7 Mt. Lassen Drive - Suite A-134,
11    San Rafael, CA  94903, U.S.A., +1(415)492-9861, for further information.
12 */
13 /*$Id: gswts.c 10104 2009-10-01 01:22:52Z masaki $ */
14 /* Screen generation for Well Tempered Screening. */
15 #include "stdpre.h"
16 #include <stdlib.h> /* for malloc */
17 #include "gx.h"
18 #include "gp.h" /* for persistent cache */
19 #include "gxstate.h"
20 #include "gsht.h"
21 #include "math_.h"
22 #include "string_.h"
23 #include "gserrors.h"
24 #include "gxfrac.h"
25 #include "gxwts.h"
26 #include "gswts.h"
27 
28 #define noDUMP_WTS
29 #ifdef DUMP_WTS
30 #include "fcntl_.h"
31 #endif
32 
33 #define noVERBOSE
34 
35 #ifdef UNIT_TEST
36 /**
37  * This file can be compiled independently for unit testing purposes.
38  * Try this invocation:
39  *
40  * gcc -I../obj -DUNIT_TEST gswts.c gxwts.c -o test_gswts -lm
41  * ./test_gswts | sed '/P5/,$!d' | xv -
42  **/
43 #undef printf
44 #undef stdout
45 #undef dlprintf1
46 #define dlprintf1 printf
47 #undef dlprintf2
48 #define dlprintf2 printf
49 #undef dlprintf3
50 #define dlprintf3 printf
51 #undef dlprintf4
52 #define dlprintf4 printf
53 #undef dlprintf5
54 #define dlprintf5 printf
55 #undef dlprintf6
56 #define dlprintf6 printf
57 #undef dlprintf7
58 #define dlprintf7 printf
59 
60 #endif
61 
62 /* A datatype used for representing the product of two 32 bit integers.
63    If a 64 bit integer type is available, it may be a better choice. */
64 typedef double wts_bigint;
65 
66 typedef struct gx_wts_cell_params_j_s gx_wts_cell_params_j_t;
67 typedef struct gx_wts_cell_params_h_s gx_wts_cell_params_h_t;
68 
69 struct gx_wts_cell_params_j_s {
70     gx_wts_cell_params_t base;
71     int shift;
72 
73     double ufast_a;
74     double vfast_a;
75     double uslow_a;
76     double vslow_a;
77 
78     /* Probabilities of "jumps". A and B jumps can happen when moving
79        one pixel to the right. C and D can happen when moving one pixel
80        down. */
81     double pa;
82     double pb;
83     double pc;
84     double pd;
85 
86     int xa;
87     int ya;
88     int xb;
89     int yb;
90     int xc;
91     int yc;
92     int xd;
93     int yd;
94 };
95 
96 struct gx_wts_cell_params_h_s {
97     gx_wts_cell_params_t base;
98 
99     /* This is the exact value that x1 and (width-x1) approximates. */
100     double px;
101     /* Ditto y1 and (height-y1). */
102     double py;
103 
104     int x1;
105     int y1;
106 };
107 
108 #define WTS_EXTRA_SORT_BITS 9
109 #define WTS_SORTED_MAX (((frac_1) << (WTS_EXTRA_SORT_BITS)) - 1)
110 
111 typedef struct {
112     int u;
113     int v;
114     int k;
115     int l;
116 } wts_vec_t;
117 
118 static void
wts_vec_set(wts_vec_t * wv,int u,int v,int k,int l)119 wts_vec_set(wts_vec_t *wv, int u, int v, int k, int l)
120 {
121     wv->u = u;
122     wv->v = v;
123     wv->k = k;
124     wv->l = l;
125 }
126 
127 static wts_bigint
wts_vec_smag(const wts_vec_t * a)128 wts_vec_smag(const wts_vec_t *a)
129 {
130     wts_bigint u = a->u;
131     wts_bigint v = a->v;
132     return u * u + v * v;
133 }
134 
135 static void
wts_vec_add(wts_vec_t * a,const wts_vec_t * b,const wts_vec_t * c)136 wts_vec_add(wts_vec_t *a, const wts_vec_t *b, const wts_vec_t *c)
137 {
138     a->u = b->u + c->u;
139     a->v = b->v + c->v;
140     a->k = b->k + c->k;
141     a->l = b->l + c->l;
142 }
143 
144 static void
wts_vec_sub(wts_vec_t * a,const wts_vec_t * b,const wts_vec_t * c)145 wts_vec_sub(wts_vec_t *a, const wts_vec_t *b, const wts_vec_t *c)
146 {
147     a->u = b->u - c->u;
148     a->v = b->v - c->v;
149     a->k = b->k - c->k;
150     a->l = b->l - c->l;
151 }
152 
153 /**
154  * wts_vec_gcd3: Compute 3-way "gcd" of three vectors.
155  * @a, @b, @c: Vectors.
156  *
157  * Compute pair of vectors satisfying three constraints:
158  * They are integer combinations of the source vectors.
159  * The source vectors are integer combinations of the results.
160  * The magnitude of the vectors is minimized.
161  *
162  * The algorithm for this computation is quite reminiscent of the
163  * classic Euclid GCD algorithm for integers.
164  *
165  * On return, @a and @b contain the result. @c is destroyed.
166  **/
167 static void
wts_vec_gcd3(wts_vec_t * a,wts_vec_t * b,wts_vec_t * c)168 wts_vec_gcd3(wts_vec_t *a, wts_vec_t *b, wts_vec_t *c)
169 {
170     wts_vec_t d, e;
171     double ra, rb, rc, rd, re;
172 
173     wts_vec_set(&d, 0, 0, 0, 0);
174     wts_vec_set(&e, 0, 0, 0, 0);
175     for (;;) {
176 	ra = wts_vec_smag(a);
177 	rb = wts_vec_smag(b);
178 	rc = wts_vec_smag(c);
179 	wts_vec_sub(&d, a, b);
180 	wts_vec_add(&e, a, b);
181 	rd = wts_vec_smag(&d);
182 	re = wts_vec_smag(&e);
183 	if (re && re < rd) {
184 	    d = e;
185 	    rd = re;
186 	}
187 	if (rd && rd < ra && ra >= rb) *a = d;
188 	else if (rd < rb && rb > ra) *b = d;
189 	else {
190 	    wts_vec_sub(&d, a, c);
191 	    wts_vec_add(&e, a, c);
192 	    rd = wts_vec_smag(&d);
193 	    re = wts_vec_smag(&e);
194 	    if (re < rd) {
195 		d = e;
196 		rd = re;
197 	    }
198 	    if (rd && rd < ra && ra >= rc) *a = d;
199 	    else if (rd < rc && rc > ra) *c = d;
200 	    else {
201 		wts_vec_sub(&d, b, c);
202 		wts_vec_add(&e, b, c);
203 		rd = wts_vec_smag(&d);
204 		re = wts_vec_smag(&e);
205 		if (re < rd) {
206 		    d = e;
207 		    rd = re;
208 		}
209 		if (rd && rd < rb && rb >= rc) *b = d;
210 		else if (rd < rc && rc > rb) *c = d;
211 		else
212 		    break;
213 	    }
214 	}
215     }
216 }
217 
218 static wts_bigint
wts_vec_cross(const wts_vec_t * a,const wts_vec_t * b)219 wts_vec_cross(const wts_vec_t *a, const wts_vec_t *b)
220 {
221     wts_bigint au = a->u;
222     wts_bigint av = a->v;
223     wts_bigint bu = b->u;
224     wts_bigint bv = b->v;
225     return au * bv - av * bu;
226 }
227 
228 static void
wts_vec_neg(wts_vec_t * a)229 wts_vec_neg(wts_vec_t *a)
230 {
231     a->u = -a->u;
232     a->v = -a->v;
233     a->k = -a->k;
234     a->l = -a->l;
235 }
236 
237 /* compute k mod m */
238 static void
wts_vec_modk(wts_vec_t * a,int m)239 wts_vec_modk(wts_vec_t *a, int m)
240 {
241     while (a->k < 0) a->k += m;
242     while (a->k >= m) a->k -= m;
243 }
244 
245 /* Compute modulo in rational cell. */
246 static void
wts_vec_modkls(wts_vec_t * a,int m,int i,int s)247 wts_vec_modkls(wts_vec_t *a, int m, int i, int s)
248 {
249     while (a->l < 0) {
250 	a->l += i;
251 	a->k -= s;
252     }
253     while (a->l >= i) {
254 	a->l -= i;
255 	a->k += s;
256     }
257     while (a->k < 0) a->k += m;
258     while (a->k >= m) a->k -= m;
259 }
260 
261 static void
wts_set_mat(gx_wts_cell_params_t * wcp,double sratiox,double sratioy,double sangledeg)262 wts_set_mat(gx_wts_cell_params_t *wcp, double sratiox, double sratioy,
263 	    double sangledeg)
264 {
265     double sangle = sangledeg * M_PI / 180;
266 
267     wcp->ufast = cos(sangle) / sratiox;
268     wcp->vfast = -sin(sangle) / sratiox;
269     wcp->uslow = sin(sangle) / sratioy;
270     wcp->vslow = cos(sangle) / sratioy;
271 }
272 
273 
274 /**
275  * Calculate Screen H cell sizes.
276  **/
277 static void
wts_cell_calc_h(double inc,int * px1,int * pxwidth,double * pp1,double memw)278 wts_cell_calc_h(double inc, int *px1, int *pxwidth, double *pp1, double memw)
279 {
280     double minrep = pow(2, memw) * 50 / pow(2, 7.5);
281     int m1 = 0, m2 = 0;
282     double e1, e2;
283 
284     int uacc;
285 
286     e1 = 1e5;
287     e2 = 1e5;
288     for (uacc = (int)ceil(minrep * inc); uacc <= floor(2 * minrep * inc); uacc++) {
289 	int mt;
290 	double et;
291 
292 	mt = (int)floor(uacc / inc + 1e-5);
293 	et = uacc / inc - mt + mt * 0.001;
294 	if (et < e1) {
295 	    e1 = et;
296 	    m1 = mt;
297 	}
298 	mt = (int)ceil(uacc / inc - 1e-5);
299 	et = mt - uacc / inc + mt * 0.001;
300 	if (et < e2) {
301 	    e2 = et;
302 	    m2 = mt;
303 	}
304     }
305     if (m1 == m2) {
306 	*px1 = m1;
307 	*pxwidth = m1;
308 	*pp1 = 1.0;
309     } else {
310 	*px1 = m1;
311 	*pxwidth = m1 + m2;
312 	e1 = fabs(m1 * inc - floor(0.5 + m1 * inc));
313 	e2 = fabs(m2 * inc - floor(0.5 + m2 * inc));
314 	*pp1 = e2 / (e1 + e2);
315     }
316 }
317 
318 /* Implementation for Screen H. This is optimized for 0 and 45 degree
319    rotations. */
320 static gx_wts_cell_params_t *
wts_pick_cell_size_h(double sratiox,double sratioy,double sangledeg,double memw)321 wts_pick_cell_size_h(double sratiox, double sratioy, double sangledeg,
322 		     double memw)
323 {
324     gx_wts_cell_params_h_t *wcph;
325     double xinc, yinc;
326 
327     wcph = malloc(sizeof(gx_wts_cell_params_h_t));
328     if (wcph == NULL)
329 	return NULL;
330 
331     wcph->base.t = WTS_SCREEN_H;
332     wts_set_mat(&wcph->base, sratiox, sratioy, sangledeg);
333 
334     xinc = fabs(wcph->base.ufast);
335     if (xinc == 0)
336 	xinc = fabs(wcph->base.vfast);
337 
338     wts_cell_calc_h(xinc, &wcph->x1, &wcph->base.width, &wcph->px, memw);
339 
340     yinc = fabs(wcph->base.uslow);
341     if (yinc == 0)
342 	yinc = fabs(wcph->base.vslow);
343     wts_cell_calc_h(yinc, &wcph->y1, &wcph->base.height, &wcph->py, memw);
344 
345     return &wcph->base;
346 }
347 
348 static double
wts_qart(double r,double rbase,double p,double pbase)349 wts_qart(double r, double rbase, double p, double pbase)
350 {
351    if (p < 1e-5) {
352       return ((r + rbase) * p);
353    } else {
354       return ((r + rbase) * (p + pbase) - rbase * pbase);
355    }
356 }
357 
358 #ifdef VERBOSE
359 static void
wts_print_j_jump(const gx_wts_cell_params_j_t * wcpj,const char * name,double pa,int xa,int ya)360 wts_print_j_jump(const gx_wts_cell_params_j_t *wcpj, const char *name,
361 		 double pa, int xa, int ya)
362 {
363     dlprintf6("jump %s: (%d, %d) %f, actual (%f, %f)\n",
364 	      name, xa, ya, pa,
365 	      wcpj->ufast_a * xa + wcpj->uslow_a * ya,
366 	      wcpj->vfast_a * xa + wcpj->vslow_a * ya);
367 }
368 
369 static void
wts_j_add_jump(const gx_wts_cell_params_j_t * wcpj,double * pu,double * pv,double pa,int xa,int ya)370 wts_j_add_jump(const gx_wts_cell_params_j_t *wcpj, double *pu, double *pv,
371 	       double pa, int xa, int ya)
372 {
373     double jump_u = wcpj->ufast_a * xa + wcpj->uslow_a * ya;
374     double jump_v = wcpj->vfast_a * xa + wcpj->vslow_a * ya;
375     jump_u -= floor(jump_u + 0.5);
376     jump_v -= floor(jump_v + 0.5);
377     *pu += pa * jump_u;
378     *pv += pa * jump_v;
379 }
380 
381 static void
wts_print_j(const gx_wts_cell_params_j_t * wcpj)382 wts_print_j(const gx_wts_cell_params_j_t *wcpj)
383 {
384     double uf, vf;
385     double us, vs;
386 
387     dlprintf3("cell = %d x %d, shift = %d\n",
388 	      wcpj->base.width, wcpj->base.height, wcpj->shift);
389     wts_print_j_jump(wcpj, "a", wcpj->pa, wcpj->xa, wcpj->ya);
390     wts_print_j_jump(wcpj, "b", wcpj->pb, wcpj->xb, wcpj->yb);
391     wts_print_j_jump(wcpj, "c", wcpj->pc, wcpj->xc, wcpj->yc);
392     wts_print_j_jump(wcpj, "d", wcpj->pd, wcpj->xd, wcpj->yd);
393     uf = wcpj->ufast_a;
394     vf = wcpj->vfast_a;
395     us = wcpj->uslow_a;
396     vs = wcpj->vslow_a;
397     wts_j_add_jump(wcpj, &uf, &vf, wcpj->pa, wcpj->xa, wcpj->ya);
398     wts_j_add_jump(wcpj, &uf, &vf, wcpj->pb, wcpj->xb, wcpj->yb);
399     wts_j_add_jump(wcpj, &us, &vs, wcpj->pc, wcpj->xc, wcpj->yc);
400     wts_j_add_jump(wcpj, &us, &vs, wcpj->pd, wcpj->xd, wcpj->yd);
401     dlprintf6("d: %f, %f; a: %f, %f; err: %g, %g\n",
402 	    wcpj->base.ufast, wcpj->base.vfast,
403 	    wcpj->ufast_a, wcpj->vfast_a,
404 	    wcpj->base.ufast - uf, wcpj->base.vfast - vf);
405     dlprintf6("d: %f, %f; a: %f, %f; err: %g, %g\n",
406 	    wcpj->base.uslow, wcpj->base.vslow,
407 	    wcpj->uslow_a, wcpj->vslow_a,
408 	    wcpj->base.uslow - us, wcpj->base.vslow - vs);
409 }
410 #endif
411 
412 /**
413  * wts_set_scr_jxi_try: Try a width for setting Screen J parameters.
414  * @wcpj: Screen J parameters.
415  * @m: Width to try.
416  * @qb: Best quality score seen so far.
417  * @jmem: Weight given to memory usage.
418  *
419  * Evaluate the quality for width @i. If the quality is better than
420  * @qb, then set the rest of the parameters in @wcpj.
421  *
422  * This routine corresponds to SetScrJXITrySP in the WTS source.
423  *
424  * Return value: Quality score for parameters chosen.
425  **/
426 static double
wts_set_scr_jxi_try(gx_wts_cell_params_j_t * wcpj,int m,double qb,double jmem)427 wts_set_scr_jxi_try(gx_wts_cell_params_j_t *wcpj, int m, double qb,
428 		    double jmem)
429 {
430     const double uf = wcpj->base.ufast;
431     const double vf = wcpj->base.vfast;
432     const double us = wcpj->base.uslow;
433     const double vs = wcpj->base.vslow;
434     wts_vec_t a, b, c;
435     double ufj, vfj;
436     const double rbase = 0.1;
437     const double pbase = 0.001;
438     double ug, vg;
439     double pp, pa, pb;
440     double rf;
441     double xa, ya, ra, xb, yb, rb;
442     double q, qx, qw, qxl, qbi;
443     double pya, pyb;
444     int i;
445     bool jumpok;
446 
447     wts_vec_set(&a, (int)floor(uf * m + 0.5), (int)floor(vf * m + 0.5), 1, 0);
448     if (a.u == 0 && a.v == 0)
449 	return qb + 1;
450 
451     ufj = a.u / (double)m;
452     vfj = a.v / (double)m;
453     /* (ufj, vfj) = movement in UV space from (0, 1) in XY space */
454 
455     wts_vec_set(&b, m, 0, 0, 0);
456     wts_vec_set(&c, 0, m, 0, 0);
457     wts_vec_gcd3(&a, &b, &c);
458     ug = (uf - ufj) * m;
459     vg = (vf - vfj) * m;
460     pp = 1.0 / wts_vec_cross(&b, &a);
461     pa = (b.u * vg - ug * b.v) * pp;
462     pb = (ug * a.v - a.u * vg) * pp;
463     if (pa < 0) {
464 	wts_vec_neg(&a);
465 	pa = -pa;
466     }
467     if (pb < 0) {
468 	wts_vec_neg(&b);
469 	pb = -pb;
470     }
471     wts_vec_modk(&a, m);
472     wts_vec_modk(&b, m);
473 
474     /* Prefer nontrivial jumps. */
475     jumpok = (wcpj->xa == 0 && wcpj->ya == 0) ||
476       (wcpj->xb == 0 && wcpj->yb == 0) ||
477       (a.k != 0 && b.k != 0);
478 
479     rf = (uf * uf + vf * vf) * m;
480     xa = (a.u * uf + a.v * vf) / rf;
481     ya = (a.v * uf - a.u * vf) / rf;
482     xb = (b.u * uf + b.v * vf) / rf;
483     yb = (b.v * uf - b.u * vf) / rf;
484     ra = sqrt(xa * xa + 0.0625 * ya * ya);
485     rb = sqrt(xb * xb + 0.0625 * yb * yb);
486     qx = 0.5 * (wts_qart(ra, rbase, pa, pbase) +
487 		wts_qart(rb, rbase, pb, pbase));
488     if (qx < 1.0 / 4000.0)
489 	qx *= 0.25;
490     else
491 	qx -= 0.75 / 4000.0;
492     if (m < 7500)
493 	qw = 0;
494     else
495 	qw = 0.00025; /* cache penalty */
496     qxl = qx + qw;
497     if (qxl > qb)
498 	return qxl;
499 
500     /* width is ok, now try heights */
501 
502     pp = m / (double)wts_vec_cross(&b, &a);
503     pya = (b.u * vs - us * b.v) * pp;
504     pyb = (us * a.v - a.u * vs) * pp;
505 
506     qbi = qb;
507     for (i = 1; qxl + m * i * jmem < qbi; i++) {
508 	int s = m * i;
509 	int ca, cb;
510 	double usj, vsj;
511 	double usg, vsg;
512 	wts_vec_t a1, b1, a2, b2;
513 	double pc, pd;
514 	int ck;
515 	double qy, qm;
516 
517 	ca = (int)floor(i * pya + 0.5);
518 	cb = (int)floor(i * pyb + 0.5);
519 	wts_vec_set(&c, ca * a.u + cb * b.u, ca * a.v + cb * b.v, 0, 1);
520 	usj = c.u / (double)s;
521 	vsj = c.v / (double)s;
522 	usg = (us - usj);
523 	vsg = (vs - vsj);
524 
525 	a1 = a;
526 	b1 = b;
527 	a1.u *= i;
528 	a1.v *= i;
529 	b1.u *= i;
530 	b1.v *= i;
531 	wts_vec_gcd3(&a1, &b1, &c);
532 	a2 = a1;
533 	b2 = b1;
534 	pp = s / (double)wts_vec_cross(&b1, &a1);
535 	pc = (b1.u * vsg - usg * b1.v) * pp;
536 	pd = (usg * a1.v - a1.u * vsg) * pp;
537 	if (pc < 0) {
538 	    wts_vec_neg(&a1);
539 	    pc = -pc;
540 	}
541 	if (pd < 0) {
542 	    wts_vec_neg(&b1);
543 	    pd = -pd;
544 	}
545 	ck = ca * a.k + cb * b.k;
546 	while (ck < 0) ck += m;
547 	while (ck >= m) ck -= m;
548 	wts_vec_modkls(&a1, m, i, ck);
549 	wts_vec_modkls(&b1, m, i, ck);
550 	rf = (us * us - vs * vs) * s;
551 	xa = (a1.u * us + a1.v * vs) / rf;
552 	ya = (a1.v * us - a1.u * vs) / rf;
553 	xb = (b1.u * us + b1.v * vs) / rf;
554 	yb = (b1.v * us - b1.u * vs) / rf;
555 	ra = sqrt(xa * xa + 0.0625 * ya * ya);
556 	rb = sqrt(xb * xb + 0.0625 * yb * yb);
557 	qy = 0.5 * (wts_qart(ra, rbase, pc, pbase) +
558 		    wts_qart(rb, rbase, pd, pbase));
559 	if (qy < 1.0 / 100.0)
560 	    qy *= 0.025;
561 	else
562 	    qy -= 0.75 / 100.0;
563 	qm = s * jmem;
564 
565 	/* first try a and b jumps within the scanline */
566 	q = qm + qw + qx + qy;
567 	if (q < qbi && jumpok) {
568 #ifdef VERBOSE
569 	    dlprintf7("m = %d, n = %d, q = %d, qx = %d, qy = %d, qm = %d, qw = %d\n",
570 		      m, i, (int)(q * 1e6), (int)(qx * 1e6), (int)(qy * 1e6), (int)(qm * 1e6), (int)(qw * 1e6));
571 #endif
572 	    qbi = q;
573 	    wcpj->base.width = m;
574 	    wcpj->base.height = i;
575 	    wcpj->shift = ck;
576 	    wcpj->ufast_a = ufj;
577 	    wcpj->vfast_a = vfj;
578 	    wcpj->uslow_a = usj;
579 	    wcpj->vslow_a = vsj;
580 	    wcpj->xa = a.k;
581 	    wcpj->ya = 0;
582 	    wcpj->xb = b.k;
583 	    wcpj->yb = 0;
584 	    wcpj->xc = a1.k;
585 	    wcpj->yc = a1.l;
586 	    wcpj->xd = b1.k;
587 	    wcpj->yd = b1.l;
588 	    wcpj->pa = pa;
589 	    wcpj->pb = pb;
590 	    wcpj->pc = pc;
591 	    wcpj->pd = pd;
592 #ifdef VERBOSE
593 	    wts_print_j(wcpj);
594 #endif
595 	}
596 
597 	/* then try unconstrained a and b jumps */
598 	if (i > 1) {
599 	    double pa2, pb2, pp2;
600 	    double qx2, qw2, q2;
601 
602 	    pp2 = pp;
603 	    pa2 = (b2.u * vg - ug * b2.v) * pp2;
604 	    pb2 = (ug * a2.v - a2.u * vg) * pp2;
605 	    rf = (uf * uf + vf * vf) * s;
606 	    xa = (a2.u * uf + a2.v * vf) / rf;
607 	    ya = (a2.v * uf - a2.u * vf) / rf;
608 	    xb = (b2.u * uf + b2.v * vf) / rf;
609 	    yb = (b2.v * uf - b2.u * vf) / rf;
610 	    ra = sqrt(xa * xa + 0.0625 * ya * ya);
611 	    rb = sqrt(xb * xb + 0.0625 * yb * yb);
612 	    if (pa2 < 0) {
613 		pa2 = -pa2;
614 		wts_vec_neg(&a2);
615 	    }
616 	    if (pb2 < 0) {
617 		pb2 = -pb2;
618 		wts_vec_neg(&b2);
619 	    }
620 	    wts_vec_modkls(&a2, m, i, ck);
621 	    wts_vec_modkls(&b2, m, i, ck);
622 	    qx2 = 0.5 * (wts_qart(ra, rbase, pa2, pbase) +
623 			 wts_qart(rb, rbase, pb2, pbase));
624 	    if (qx2 < 1.0 / 4000.0)
625 		qx2 *= 0.25;
626 	    else
627 		qx2 -= 0.75 / 4000.0;
628 	    if (s < 7500)
629 		qw2 = 0;
630 	    else
631 		qw2 = 0.00025; /* cache penalty */
632 	    q2 = qm + qw2 + qx2 + qy;
633 	    if (q2 < qbi) {
634 #ifdef VERBOSE
635 		dlprintf7("m = %d, n = %d, q = %d, qx2 = %d, qy = %d, qm = %d, qw2 = %d\n",
636 			  m, i, (int)(q * 1e6), (int)(qx * 1e6), (int)(qy * 1e6), (int)(qm * 1e6), (int)(qw2 * 1e6));
637 #endif
638 		if (qxl > qw2 + qx2)
639 		    qxl = qw2 + qx2;
640 		qbi = q2;
641 		wcpj->base.width = m;
642 		wcpj->base.height = i;
643 		wcpj->shift = ck;
644 		wcpj->ufast_a = ufj;
645 		wcpj->vfast_a = vfj;
646 		wcpj->uslow_a = usj;
647 		wcpj->vslow_a = vsj;
648 		wcpj->xa = a2.k;
649 		wcpj->ya = a2.l;
650 		wcpj->xb = b2.k;
651 		wcpj->yb = a2.l;
652 		wcpj->xc = a1.k;
653 		wcpj->yc = a1.l;
654 		wcpj->xd = b1.k;
655 		wcpj->yd = b1.l;
656 		wcpj->pa = pa2;
657 		wcpj->pb = pb2;
658 		wcpj->pc = pc;
659 		wcpj->pd = pd;
660 #ifdef VERBOSE
661 		wts_print_j(wcpj);
662 #endif
663 	    }
664 	} /* if (i > 1) */
665 	if (qx > 10 * qy)
666 	    break;
667     }
668     return qbi;
669 }
670 
671 static int
wts_double_to_int_cap(double d)672 wts_double_to_int_cap(double d)
673 {
674     if (d > 0x7fffffff)
675 	return 0x7fffffff;
676     else return (int)d;
677 }
678 
679 /**
680  * wts_set_scr_jxi: Choose Screen J parameters.
681  * @wcpj: Screen J parameters.
682  * @jmem: Weight given to memory usage.
683  *
684  * Chooses a cell size based on the [uv]{fast,slow} parameters,
685  * setting the rest of the parameters in @wcpj. Essentially, this
686  * algorithm iterates through all plausible widths for the cell.
687  *
688  * This routine corresponds to SetScrJXISP in the WTS source.
689  *
690  * Return value: Quality score for parameters chosen.
691  **/
692 static double
wts_set_scr_jxi(gx_wts_cell_params_j_t * wcpj,double jmem)693 wts_set_scr_jxi(gx_wts_cell_params_j_t *wcpj, double jmem)
694 {
695     int i, imax;
696     double q, qb;
697 
698     wcpj->xa = 0;
699     wcpj->ya = 0;
700     wcpj->xb = 0;
701     wcpj->yb = 0;
702     wcpj->xc = 0;
703     wcpj->yc = 0;
704     wcpj->xd = 0;
705     wcpj->yd = 0;
706 
707     qb = 1.0;
708     imax = wts_double_to_int_cap(qb / jmem);
709     for (i = 1; i <= imax; i++) {
710 	if (i > 1) {
711 	    q = wts_set_scr_jxi_try(wcpj, i, qb, jmem);
712 	    if (q < qb) {
713 		qb = q;
714 		imax = wts_double_to_int_cap(q / jmem);
715 		if (imax >= 7500)
716 		    imax = wts_double_to_int_cap((q - 0.00025) / jmem);
717 		if (imax < 7500) {
718 		    imax = 7500;
719 		}
720 	    }
721 	}
722     }
723     return qb;
724 }
725 
726 /* Implementation for Screen J. This is optimized for general angles. */
727 static gx_wts_cell_params_t *
wts_pick_cell_size_j(double sratiox,double sratioy,double sangledeg,double memw)728 wts_pick_cell_size_j(double sratiox, double sratioy, double sangledeg,
729 		     double memw)
730 {
731     gx_wts_cell_params_j_t *wcpj;
732     double code;
733 
734     wcpj = malloc(sizeof(gx_wts_cell_params_j_t));
735     if (wcpj == NULL)
736 	return NULL;
737 
738     wcpj->base.t = WTS_SCREEN_J;
739     wts_set_mat(&wcpj->base, sratiox, sratioy, sangledeg);
740 
741     code = wts_set_scr_jxi(wcpj, pow(0.1, memw));
742     if (code < 0) {
743 	free(wcpj);
744 	return NULL;
745     }
746 
747     return &wcpj->base;
748 }
749 
750 typedef struct {
751     double sratiox;
752     double sratioy;
753     double sangledeg;
754     double memw;
755 } wts_cell_size_key;
756 
757 static void *
wts_cache_alloc_callback(void * data,int bytes)758 wts_cache_alloc_callback(void *data, int bytes)
759 {
760     return malloc(bytes);
761 }
762 
763 static void
store_be32(byte * ptr,int x)764 store_be32(byte *ptr, int x)
765 {
766     ptr[0] = (x >> 24) & 0xff;
767     ptr[1] = (x >> 16) & 0xff;
768     ptr[2] = (x >> 8) & 0xff;
769     ptr[3] = x & 0xff;
770 }
771 
772 /**
773  * wts_pick_cell_size: Choose cell size for WTS screen.
774  * @ph: The halftone parameters.
775  * @pmat: Transformation from 1/72" Y-up coords to device coords.
776  *
777  * Return value: The WTS cell parameters, or NULL on error.
778  **/
779 gx_wts_cell_params_t *
wts_pick_cell_size(gs_screen_halftone * ph,const gs_matrix * pmat)780 wts_pick_cell_size(gs_screen_halftone *ph, const gs_matrix *pmat)
781 {
782     gx_wts_cell_params_t *result;
783 
784     /* todo: deal with landscape and mirrored coordinate systems */
785     double sangledeg = ph->angle;
786     double sratiox = 72.0 * fabs(pmat->xx) / ph->frequency;
787     double sratioy = 72.0 * fabs(pmat->yy) / ph->frequency;
788     double octangle;
789     double memw = 8.0;
790     wts_cell_size_key key;
791     int len;
792     byte *keyptr = NULL;
793     byte intkey[12];
794     int keysize;
795 
796     octangle = sangledeg;
797     while (octangle >= 45.0) octangle -= 45.0;
798     while (octangle < -45.0) octangle += 45.0;
799 
800     /* try persistent cache */
801     if (sangledeg == 45.0) {
802 	int srxi, sryi;
803 
804 	srxi = (int)floor(sratiox / sqrt(2) + 0.5);
805 	sryi = (int)floor(sratioy / sqrt(2) + 0.5);
806 	if (fabs(srxi * sqrt(2) - sratiox) < 2e-6 &&
807 	    fabs(sryi * sqrt(2) - sratioy) < 2e-6) {
808 	    store_be32(intkey, (int)sangledeg);
809 	    store_be32(intkey + 4, srxi);
810 	    store_be32(intkey + 8, sryi);
811 	    keyptr = intkey;
812 	    keysize = sizeof(intkey);
813 	}
814     }
815     if (keyptr == NULL) {
816 	key.sratiox = sratiox;
817 	key.sratioy = sratioy;
818 	key.sangledeg = sangledeg;
819 	key.memw = memw;
820 	keyptr = (byte *)&key;
821 	keysize = sizeof(key);
822     }
823     len = gp_cache_query(GP_CACHE_TYPE_WTS_SIZE, keyptr, keysize,
824 			 (void **)&result, wts_cache_alloc_callback, NULL);
825     if (len >= 0)
826 	return result;
827 
828     if (fabs(octangle) < 1e-4)
829 	result = wts_pick_cell_size_h(sratiox, sratioy, sangledeg, memw);
830     else
831 	result = wts_pick_cell_size_j(sratiox, sratioy, sangledeg, memw);
832 
833     if (result != NULL) {
834 	int resultsize = 0;
835 
836 	/* insert computed cell size into cache */
837 	if (result->t == WTS_SCREEN_H)
838 	    resultsize = sizeof(gx_wts_cell_params_h_t);
839 	else if (result->t == WTS_SCREEN_J)
840 	    resultsize = sizeof(gx_wts_cell_params_j_t);
841 	if (resultsize)
842 	    gp_cache_insert(GP_CACHE_TYPE_WTS_SIZE, (byte *)&key, sizeof(key),
843 			    (void *)result, resultsize);
844 
845 	ph->actual_frequency = ph->frequency;
846 	ph->actual_angle = ph->angle;
847     }
848     return result;
849 }
850 
851 struct gs_wts_screen_enum_s {
852     wts_screen_type t;
853     bits32 *cell;
854     int width;
855     int height;
856     int size;
857 
858     int idx;
859 };
860 
861 typedef struct {
862     gs_wts_screen_enum_t base;
863     const gx_wts_cell_params_j_t *wcpj;
864 } gs_wts_screen_enum_j_t;
865 
866 typedef struct {
867     gs_wts_screen_enum_t base;
868     const gx_wts_cell_params_h_t *wcph;
869 
870     /* an argument can be made that these should be in the params */
871     double ufast1, vfast1;
872     double ufast2, vfast2;
873     double uslow1, vslow1;
874     double uslow2, vslow2;
875 } gs_wts_screen_enum_h_t;
876 
877 static gs_wts_screen_enum_t *
gs_wts_screen_enum_j_new(gx_wts_cell_params_t * wcp)878 gs_wts_screen_enum_j_new(gx_wts_cell_params_t *wcp)
879 {
880     gs_wts_screen_enum_j_t *wsej;
881 
882     wsej = malloc(sizeof(gs_wts_screen_enum_j_t));
883     wsej->base.t = WTS_SCREEN_J;
884     wsej->wcpj = (const gx_wts_cell_params_j_t *)wcp;
885     wsej->base.width = wcp->width;
886     wsej->base.height = wcp->height;
887     wsej->base.size = wcp->width * wcp->height;
888     wsej->base.cell = malloc(wsej->base.size * sizeof(wsej->base.cell[0]));
889     wsej->base.idx = 0;
890 
891     return (gs_wts_screen_enum_t *)wsej;
892 }
893 
894 static int
gs_wts_screen_enum_j_currentpoint(gs_wts_screen_enum_t * self,gs_point * ppt)895 gs_wts_screen_enum_j_currentpoint(gs_wts_screen_enum_t *self,
896 				  gs_point *ppt)
897 {
898     gs_wts_screen_enum_j_t *z = (gs_wts_screen_enum_j_t *)self;
899     const gx_wts_cell_params_j_t *wcpj = z->wcpj;
900 
901     int x, y;
902     double u, v;
903 
904     if (z->base.idx == z->base.size) {
905 	return 1;
906     }
907     x = z->base.idx % wcpj->base.width;
908     y = z->base.idx / wcpj->base.width;
909     u = wcpj->ufast_a * x + wcpj->uslow_a * y;
910     v = wcpj->vfast_a * x + wcpj->vslow_a * y;
911     u -= floor(u);
912     v -= floor(v);
913     ppt->x = 2 * u - 1;
914     ppt->y = 2 * v - 1;
915     return 0;
916 }
917 
918 static gs_wts_screen_enum_t *
gs_wts_screen_enum_h_new(gx_wts_cell_params_t * wcp)919 gs_wts_screen_enum_h_new(gx_wts_cell_params_t *wcp)
920 {
921     gs_wts_screen_enum_h_t *wseh;
922     const gx_wts_cell_params_h_t *wcph = (const gx_wts_cell_params_h_t *)wcp;
923     int x1 = wcph->x1;
924     int x2 = wcp->width - x1;
925     int y1 = wcph->y1;
926     int y2 = wcp->height - y1;
927 
928     wseh = malloc(sizeof(gs_wts_screen_enum_h_t));
929     wseh->base.t = WTS_SCREEN_H;
930     wseh->wcph = wcph;
931     wseh->base.width = wcp->width;
932     wseh->base.height = wcp->height;
933     wseh->base.size = wcp->width * wcp->height;
934     wseh->base.cell = malloc(wseh->base.size * sizeof(wseh->base.cell[0]));
935     wseh->base.idx = 0;
936 
937     wseh->ufast1 = floor(0.5 + wcp->ufast * x1) / x1;
938     wseh->vfast1 = floor(0.5 + wcp->vfast * x1) / x1;
939     if (x2 > 0) {
940 	wseh->ufast2 = floor(0.5 + wcp->ufast * x2) / x2;
941 	wseh->vfast2 = floor(0.5 + wcp->vfast * x2) / x2;
942     }
943     wseh->uslow1 = floor(0.5 + wcp->uslow * y1) / y1;
944     wseh->vslow1 = floor(0.5 + wcp->vslow * y1) / y1;
945     if (y2 > 0) {
946 	wseh->uslow2 = floor(0.5 + wcp->uslow * y2) / y2;
947 	wseh->vslow2 = floor(0.5 + wcp->vslow * y2) / y2;
948     }
949 
950     return &wseh->base;
951 }
952 
953 static int
gs_wts_screen_enum_h_currentpoint(gs_wts_screen_enum_t * self,gs_point * ppt)954 gs_wts_screen_enum_h_currentpoint(gs_wts_screen_enum_t *self,
955 				  gs_point *ppt)
956 {
957     gs_wts_screen_enum_h_t *z = (gs_wts_screen_enum_h_t *)self;
958     const gx_wts_cell_params_h_t *wcph = z->wcph;
959 
960     int x, y;
961     double u, v;
962 
963     if (self->idx == self->size) {
964 	return 1;
965     }
966     x = self->idx % wcph->base.width;
967     y = self->idx / wcph->base.width;
968     if (x < wcph->x1) {
969 	u = z->ufast1 * x;
970 	v = z->vfast1 * x;
971     } else {
972 	u = z->ufast2 * (x - wcph->x1);
973 	v = z->vfast2 * (x - wcph->x1);
974     }
975     if (y < wcph->y1) {
976 	u += z->uslow1 * y;
977 	v += z->vslow1 * y;
978     } else {
979 	u += z->uslow2 * (y - wcph->y1);
980 	v += z->vslow2 * (y - wcph->y1);
981     }
982     u -= floor(u);
983     v -= floor(v);
984     ppt->x = 2 * u - 1;
985     ppt->y = 2 * v - 1;
986     return 0;
987 }
988 
989 gs_wts_screen_enum_t *
gs_wts_screen_enum_new(gx_wts_cell_params_t * wcp)990 gs_wts_screen_enum_new(gx_wts_cell_params_t *wcp)
991 {
992     if (wcp->t == WTS_SCREEN_J)
993 	return gs_wts_screen_enum_j_new(wcp);
994     else if (wcp->t == WTS_SCREEN_H)
995 	return gs_wts_screen_enum_h_new(wcp);
996     else
997 	return NULL;
998 }
999 
1000 int
gs_wts_screen_enum_currentpoint(gs_wts_screen_enum_t * wse,gs_point * ppt)1001 gs_wts_screen_enum_currentpoint(gs_wts_screen_enum_t *wse, gs_point *ppt)
1002 {
1003     if (wse->t == WTS_SCREEN_J)
1004 	return gs_wts_screen_enum_j_currentpoint(wse, ppt);
1005     if (wse->t == WTS_SCREEN_H)
1006 	return gs_wts_screen_enum_h_currentpoint(wse, ppt);
1007     else
1008 	return -1;
1009 }
1010 
1011 int
gs_wts_screen_enum_next(gs_wts_screen_enum_t * wse,floatp value)1012 gs_wts_screen_enum_next(gs_wts_screen_enum_t *wse, floatp value)
1013 {
1014     bits32 sample;
1015 
1016     if (value < -1.0 || value > 1.0)
1017 	return_error(gs_error_rangecheck);
1018     sample = (bits32) ((value + 1) * 0x7fffffff);
1019     wse->cell[wse->idx] = sample;
1020     wse->idx++;
1021     return 0;
1022 }
1023 
1024 /* Run the enum with a square dot. This is useful for testing. */
1025 #ifdef UNIT_TEST
1026 static void
wts_run_enum_squaredot(gs_wts_screen_enum_t * wse)1027 wts_run_enum_squaredot(gs_wts_screen_enum_t *wse)
1028 {
1029     int code;
1030     gs_point pt;
1031     floatp spot;
1032 
1033     for (;;) {
1034 	code = gs_wts_screen_enum_currentpoint(wse, &pt);
1035 	if (code)
1036 	    break;
1037 	spot = 0.5 * (cos(pt.x * M_PI) + cos(pt.y * M_PI));
1038 	gs_wts_screen_enum_next(wse, spot);
1039     }
1040 }
1041 #endif /* UNIT_TEST */
1042 
1043 static int
wts_sample_cmp(const void * av,const void * bv)1044 wts_sample_cmp(const void *av, const void *bv) {
1045     const bits32 *const *a = (const bits32 *const *)av;
1046     const bits32 *const *b = (const bits32 *const *)bv;
1047 
1048     if (**a > **b) return 1;
1049     if (**a < **b) return -1;
1050     return 0;
1051 }
1052 
1053 /* This implementation simply sorts the threshold values (evening the
1054    distribution), without applying any moire reduction. */
1055 int
wts_sort_cell(gs_wts_screen_enum_t * wse)1056 wts_sort_cell(gs_wts_screen_enum_t *wse)
1057 {
1058     int size = wse->width * wse->height;
1059     bits32 *cell = wse->cell;
1060     bits32 **pcell;
1061     int i;
1062 
1063     pcell = malloc(size * sizeof(bits32 *));
1064     if (pcell == NULL)
1065 	return -1;
1066     for (i = 0; i < size; i++)
1067 	pcell[i] = &cell[i];
1068     qsort(pcell, size, sizeof(bits32 *), wts_sample_cmp);
1069     for (i = 0; i < size; i++)
1070 	*pcell[i] = (bits32)floor(WTS_SORTED_MAX * (i + 0.5) / size);
1071     free(pcell);
1072     return 0;
1073 }
1074 
1075 /**
1076  * wts_blue_bump: Generate bump function for BlueDot.
1077  *
1078  * Return value: newly allocated bump.
1079  **/
1080 static bits32 *
wts_blue_bump(const gs_wts_screen_enum_t * wse)1081 wts_blue_bump(const gs_wts_screen_enum_t *wse)
1082 {
1083     const gx_wts_cell_params_t *wcp;
1084     int width = wse->width;
1085     int height = wse->height;
1086     int shift;
1087     int size = width * height;
1088     bits32 *bump;
1089     int i;
1090     double uf, vf;
1091     double am, eg;
1092     int z;
1093     int x0, y0;
1094     int x, y;
1095 
1096     if (wse->t == WTS_SCREEN_J) {
1097 	gs_wts_screen_enum_j_t *wsej = (gs_wts_screen_enum_j_t *)wse;
1098 	shift = wsej->wcpj->shift;
1099 	wcp = &wsej->wcpj->base;
1100     } else if (wse->t == WTS_SCREEN_H) {
1101 	gs_wts_screen_enum_h_t *wseh = (gs_wts_screen_enum_h_t *)wse;
1102 	shift = 0;
1103 	wcp = &wseh->wcph->base;
1104     } else
1105 	return NULL;
1106 
1107     bump = (bits32 *)malloc(size * sizeof(bits32));
1108     if (bump == NULL)
1109 	return NULL;
1110 
1111     for (i = 0; i < size; i++)
1112 	bump[i] = 0;
1113     /* todo: more intelligence for anisotropic scaling */
1114     uf = wcp->ufast;
1115     vf = wcp->vfast;
1116 
1117     am = uf * uf + vf * vf;
1118     eg = (1 << 24) * 2.0 * sqrt (am);
1119     z = (int)(5 / sqrt (am));
1120 
1121     x0 = -(z / width) * shift - z;
1122     y0 = -(z % width);
1123     while (y0 < 0) {
1124 	x0 -= shift;
1125 	y0 += height;
1126     }
1127     while (x0 < 0) x0 += width;
1128     for (y = -z; y <= z; y++) {
1129 	int x1 = x0;
1130 	for (x = -z; x <= z; x++) {
1131 	    bump[y0 * width + x1] += (bits32)(eg * exp (-am * (x * x + y * y)));
1132 	    x1++;
1133 	    if (x1 == width)
1134 		x1 = 0;
1135 	}
1136 	y0++;
1137 	if (y0 == height) {
1138 	    x0 += shift;
1139 	    if (x0 >= width) x0 -= width;
1140 	    y0 = 0;
1141 	}
1142     }
1143     return bump;
1144 }
1145 
1146 /**
1147  * wts_sort_blue: Sort cell using BlueDot.
1148  **/
1149 static int
wts_sort_blue(const gs_wts_screen_enum_t * wse)1150 wts_sort_blue(const gs_wts_screen_enum_t *wse)
1151 {
1152     bits32 *cell = wse->cell;
1153     int width = wse->width;
1154     int height = wse->height;
1155     int shift;
1156     int size = width * height;
1157     bits32 *ref;
1158     bits32 **pcell;
1159     bits32 *bump;
1160     int i;
1161 
1162     if (wse->t == WTS_SCREEN_J) {
1163 	gs_wts_screen_enum_j_t *wsej = (gs_wts_screen_enum_j_t *)wse;
1164 	shift = wsej->wcpj->shift;
1165     } else
1166 	shift = 0;
1167 
1168     ref = (bits32 *)malloc(size * sizeof(bits32));
1169     pcell = (bits32 **)malloc(size * sizeof(bits32 *));
1170     bump = wts_blue_bump(wse);
1171     if (ref == NULL || pcell == NULL || bump == NULL) {
1172 	free(ref);
1173 	free(pcell);
1174 	free(bump);
1175 	return -1;
1176     }
1177     for (i = 0; i < size; i++)
1178 	pcell[i] = &cell[i];
1179     qsort(pcell, size, sizeof(bits32 *), wts_sample_cmp);
1180     /* set ref to sorted cell; pcell will now point to ref */
1181     for (i = 0; i < size; i++) {
1182 	pcell[i] = (pcell[i] - cell) + ref;
1183 	*pcell[i] = (bits32)floor((1 << 24) * (i + 0.5) / size);
1184 	cell[i] = 0;
1185     }
1186 
1187     for (i = 0; i < size; i++) {
1188 	bits32 gmin = *pcell[i];
1189 	int j;
1190 	int j_end = i + 5000;
1191 	int jmin = i;
1192 	int ix;
1193 	int x0, y0;
1194 	int x, y;
1195 	int ref_ix, bump_ix;
1196 
1197 	/* find minimum cell value, but prune search */
1198 	if (j_end > size) j_end = size;
1199 	for (j = i + 1; j < j_end; j++) {
1200 	    if (*pcell[j] < gmin) {
1201 		gmin = *pcell[j];
1202 		jmin = j;
1203 	    }
1204 	}
1205 	ix = pcell[jmin] - ref;
1206 	pcell[jmin] = pcell[i];
1207 	cell[ix] = (bits32)floor(WTS_SORTED_MAX * (i + 0.5) / size);
1208 
1209 	x0 = ix % width;
1210 	y0 = ix / width;
1211 
1212 	/* Add in bump, centered at (x0, y0) */
1213 	ref_ix = y0 * width;
1214 	bump_ix = 0;
1215 	for (y = 0; y < height; y++) {
1216 	    for (x = x0; x < width; x++)
1217 		ref[ref_ix + x] += bump[bump_ix++];
1218 	    for (x = 0; x < x0; x++)
1219 		ref[ref_ix + x] += bump[bump_ix++];
1220 	    ref_ix += width;
1221 	    y0++;
1222 	    if (y0 == height) {
1223 		x0 += shift;
1224 		if (x0 >= width) x0 -= width;
1225 		y0 = 0;
1226 		ref_ix = 0;
1227 	    }
1228 	}
1229 
1230 	/* Remove DC component to avoid integer overflow. */
1231 	if ((i & 255) == 255 && i + 1 < size) {
1232 	    bits32 gmin = *pcell[i + 1];
1233 	    int j;
1234 
1235 	    for (j = i + 2; j < size; j++) {
1236 		if (*pcell[j] < gmin) {
1237 		    gmin = *pcell[j];
1238 		}
1239 	    }
1240 #ifdef VERBOSE
1241 	    if_debug1('h', "[h]gmin = %d\n", gmin);
1242 #endif
1243 	    for (j = i + 1; j < size; j++)
1244 		*pcell[j] -= gmin;
1245 
1246 	}
1247     }
1248 
1249     free(ref);
1250     free(pcell);
1251     free(bump);
1252     return 0;
1253 }
1254 
1255 static wts_screen_t *
wts_screen_from_enum_j(const gs_wts_screen_enum_t * wse)1256 wts_screen_from_enum_j(const gs_wts_screen_enum_t *wse)
1257 {
1258     const gs_wts_screen_enum_j_t *wsej = (const gs_wts_screen_enum_j_t *)wse;
1259     wts_screen_j_t *wsj;
1260     wts_screen_sample_t *samples;
1261     int size;
1262     int i;
1263 
1264     wsj = malloc(sizeof(wts_screen_j_t));
1265     wsj->base.type = WTS_SCREEN_J;
1266     wsj->base.cell_width = wsej->base.width;
1267     wsj->base.cell_height = wsej->base.height;
1268     size = wsj->base.cell_width * wsj->base.cell_height;
1269     wsj->base.cell_shift = wsej->wcpj->shift;
1270     wsj->pa = (int)floor(wsej->wcpj->pa * (1 << 16) + 0.5);
1271     wsj->pb = (int)floor(wsej->wcpj->pb * (1 << 16) + 0.5);
1272     wsj->pc = (int)floor(wsej->wcpj->pc * (1 << 16) + 0.5);
1273     wsj->pd = (int)floor(wsej->wcpj->pd * (1 << 16) + 0.5);
1274     wsj->XA = wsej->wcpj->xa;
1275     wsj->YA = wsej->wcpj->ya;
1276     wsj->XB = wsej->wcpj->xb;
1277     wsj->YB = wsej->wcpj->yb;
1278     wsj->XC = wsej->wcpj->xc;
1279     wsj->YC = wsej->wcpj->yc;
1280     wsj->XD = wsej->wcpj->xd;
1281     wsj->YD = wsej->wcpj->yd;
1282 
1283     samples = malloc(sizeof(wts_screen_sample_t) * size);
1284     wsj->base.samples = samples;
1285     for (i = 0; i < size; i++) {
1286 	samples[i] = wsej->base.cell[i] >> WTS_EXTRA_SORT_BITS;
1287     }
1288 
1289 #ifdef WTS_CACHE_SIZE_X
1290     for (i = 0; i < WTS_CACHE_SIZE_X; i++)
1291 	wsj->xcache[i].tag = -1;
1292     for (i = 0; i < WTS_CACHE_SIZE_Y; i++)
1293 	wsj->ycache[i].tag = -1;
1294 #endif
1295 
1296     return &wsj->base;
1297 }
1298 
1299 static wts_screen_t *
wts_screen_from_enum_h(const gs_wts_screen_enum_t * wse)1300 wts_screen_from_enum_h(const gs_wts_screen_enum_t *wse)
1301 {
1302     const gs_wts_screen_enum_h_t *wseh = (const gs_wts_screen_enum_h_t *)wse;
1303     wts_screen_h_t *wsh;
1304     wts_screen_sample_t *samples;
1305     int size;
1306     int i;
1307 
1308     /* factor some of this out into a common init routine? */
1309     wsh = malloc(sizeof(wts_screen_h_t));
1310     wsh->base.type = WTS_SCREEN_H;
1311     wsh->base.cell_width = wseh->base.width;
1312     wsh->base.cell_height = wseh->base.height;
1313     size = wsh->base.cell_width * wsh->base.cell_height;
1314     wsh->base.cell_shift = 0;
1315     wsh->px = wseh->wcph->px;
1316     wsh->py = wseh->wcph->py;
1317     wsh->x1 = wseh->wcph->x1;
1318     wsh->y1 = wseh->wcph->y1;
1319 
1320     samples = malloc(sizeof(wts_screen_sample_t) * size);
1321     wsh->base.samples = samples;
1322     for (i = 0; i < size; i++) {
1323 	samples[i] = wseh->base.cell[i] >> WTS_EXTRA_SORT_BITS;
1324     }
1325 
1326     return &wsh->base;
1327 }
1328 
1329 typedef struct {
1330     wts_screen_type t;
1331     int width;
1332     int height;
1333 } wts_key_j;
1334 
1335 typedef struct {
1336     wts_screen_type t;
1337     int width;
1338     int height;
1339 } wts_key_h;
1340 
1341 wts_screen_t *
wts_screen_from_enum(const gs_wts_screen_enum_t * wse)1342 wts_screen_from_enum(const gs_wts_screen_enum_t *wse)
1343 {
1344     wts_screen_t *result = NULL;
1345     byte *key = NULL;
1346     int key_size = 0; /* A stub. Was uninitialized when wse->t != WTS_SCREEN_J */
1347     int cell_off;
1348     int cell_len = -1;
1349     byte *cell_result;
1350 
1351     if (wse->t == WTS_SCREEN_J) {
1352 	wts_key_j k;
1353 	k.t = wse->t;
1354 	k.width = wse->width;
1355 	k.height = wse->height;
1356 	cell_off = sizeof(k);
1357 
1358 	key_size = cell_off + wse->size * sizeof(bits32);
1359 	key = (byte *)malloc(key_size);
1360 	/* todo: more params */
1361 	memcpy(key, &k, cell_off);
1362 	memcpy(key + cell_off, wse->cell, wse->size * sizeof(bits32));
1363     } else if (wse->t == WTS_SCREEN_H) {
1364 	wts_key_h k;
1365 	k.t = wse->t;
1366 	k.width = wse->width;
1367 	k.height = wse->height;
1368 	key_size = sizeof(k);
1369 	key = (byte *)malloc(key_size);
1370 	/* todo: more params */
1371 	memcpy(key, &k, key_size);
1372     }
1373 
1374     if (key != NULL)
1375 	cell_len = gp_cache_query(GP_CACHE_TYPE_WTS_CELL, key, key_size,
1376 				  (void **)&cell_result,
1377 				  wts_cache_alloc_callback, NULL);
1378     if (cell_len >= 0) {
1379 	memcpy(wse->cell, cell_result, cell_len);
1380 	free(cell_result);
1381     } else {
1382 	wts_sort_blue(wse);
1383 	cell_len = wse->size * sizeof(bits32);
1384 	gp_cache_insert(GP_CACHE_TYPE_WTS_CELL, key, key_size,
1385 			(void *)wse->cell, cell_len);
1386     }
1387     free(key);
1388 
1389     if (wse->t == WTS_SCREEN_J)
1390 	result = wts_screen_from_enum_j(wse);
1391     else if (wse->t == WTS_SCREEN_H)
1392 	result = wts_screen_from_enum_h(wse);
1393 
1394 #ifdef DUMP_WTS
1395     {
1396 	static int dump_idx = 0;
1397 	char dump_fn[128];
1398 	int dump_fd;
1399 	byte *buf;
1400 	int size;
1401 
1402 	size = gs_wts_to_buf(result, &buf);
1403 	sprintf(dump_fn, "wts_dump_%d", dump_idx++);
1404 	dump_fd = open(dump_fn, O_WRONLY | O_CREAT | O_TRUNC, 0666);
1405 	write(dump_fd, buf, size);
1406 	close(dump_fd);
1407 	free(buf);
1408     }
1409 #endif
1410 
1411     return result;
1412 }
1413 
1414 void
gs_wts_free_enum(gs_wts_screen_enum_t * wse)1415 gs_wts_free_enum(gs_wts_screen_enum_t *wse)
1416 {
1417     free(wse);
1418 }
1419 
1420 void
gs_wts_free_screen(wts_screen_t * wts)1421 gs_wts_free_screen(wts_screen_t * wts)
1422 {
1423     /* todo: free cell */
1424     free(wts);
1425 }
1426 
1427 /* 20090929:mu: dirty fix to deal with bug #690710.
1428    wts plane files (wts_plane_{0,1,2,3}) are not maintained long time
1429    and contains littile endian 32 bit machine specific data.
1430    this quick fix only adresses reading this file on 64 bit and big endian machines.
1431    writing of those files are not supported this time. */
1432 /* short le_s16(byte *p), ushort le_u16(byte *p) */
1433 #define le_u16(p) (ushort)(*(byte *)(p) | *((byte *)(p) + 1) << 8)
1434 #define le_s16(p) (short)le_u16(p)
1435 /* int le_s32(byte *p), uint le_u32(byte *p) */
1436 #define le_u32(p) (uint)(*(byte *)(p) | *((byte *)(p) + 1) << 8 | *((byte *)(p) + 2) << 16 | *((byte *)(p) + 3) << 24)
1437 #define le_s32(p) (int)le_u32(p)
1438 
1439 int
wts_size(const wts_screen_t * ws)1440 wts_size(const wts_screen_t *ws)
1441 {
1442     int size = 0; /* A stub. Was uninitialized when none of 3 cases below. */
1443     int type = le_s32(&ws->type);
1444 
1445     if (type == WTS_SCREEN_RAT) {
1446 	size = sizeof(wts_screen_t);
1447     } else if (type == WTS_SCREEN_J) {
1448 	size = sizeof(wts_screen_j_t);
1449     } else if (type == WTS_SCREEN_H) {
1450 	size = sizeof(wts_screen_h_t);
1451     }
1452     return size;
1453 }
1454 
1455 wts_screen_t *
gs_wts_from_buf(const byte * buf,int bufsize)1456 gs_wts_from_buf(const byte *buf, int bufsize)
1457 {
1458     const wts_screen_t *ws = (const wts_screen_t *)buf;
1459     wts_screen_t *result;
1460     int size = wts_size(ws);
1461     int hdr_size;
1462     int cell_size; /* size of cell in bytes */
1463 
1464     result = (wts_screen_t *)malloc(size);
1465     if (result == NULL)
1466 	return NULL;
1467 
1468     hdr_size = offset_of(wts_screen_t, samples) + sizeof(int);
1469     if (bufsize < hdr_size ) {
1470 	free(result);
1471 	return NULL;
1472     }
1473     result->type        = le_s32(&ws->type);
1474     result->cell_width  = le_s32(&ws->cell_width);
1475     result->cell_height = le_s32(&ws->cell_height);
1476     result->cell_shift  = le_s32(&ws->cell_shift);
1477     result->samples     = NULL;
1478     if (result->type == WTS_SCREEN_J) {
1479 	wts_screen_j_t *wsj = (wts_screen_j_t *)result;
1480 	const int *wsj_params = (const int *)((const byte *)ws + hdr_size);
1481 
1482 	hdr_size += sizeof(int) * 12;
1483 	if (bufsize < hdr_size ) {
1484 	    free(result);
1485 	    return NULL;
1486 	}
1487 	wsj->pa         = le_s32(&wsj_params[0]);
1488 	wsj->pb         = le_s32(&wsj_params[1]);
1489 	wsj->pc         = le_s32(&wsj_params[2]);
1490 	wsj->pd         = le_s32(&wsj_params[3]);
1491 	wsj->XA         = le_s32(&wsj_params[4]);
1492 	wsj->YA         = le_s32(&wsj_params[5]);
1493 	wsj->XB         = le_s32(&wsj_params[6]);
1494 	wsj->YB         = le_s32(&wsj_params[7]);
1495 	wsj->XC         = le_s32(&wsj_params[8]);
1496 	wsj->YC         = le_s32(&wsj_params[9]);
1497 	wsj->XD         = le_s32(&wsj_params[10]);
1498 	wsj->YD         = le_s32(&wsj_params[11]);
1499 	/* 20090929:mu: In last version, we didn't copied these entries unless
1500 	   WTS_SCREEN_J_SIZE_NOCACHE has been defined.
1501 	   But, judging from gs_wts_to_buf() code below,
1502 	   those are written into file even WTS_SCREEN_J_SIZE_NOCACHE haven't defined
1503 	   and thereofore I think these should be read. */
1504     }
1505     /* 20090929:mu: This code doesn't care about WTS_SCREEN_H type file.
1506        If you attenpt to read such files, probably end up with errors. */
1507     cell_size = result->cell_width * result->cell_height * sizeof(wts_screen_sample_t);
1508 
1509     if (bufsize < (cell_size + hdr_size) ||
1510 	(result->samples = (wts_screen_sample_t *)malloc(cell_size)) == NULL) {
1511 	free(result);
1512 	return NULL;
1513     }
1514 #ifdef WTS_SCREEN_J_SIZE_NOCACHE
1515     if (ws->type == WTS_SCREEN_J) {
1516 	wts_screen_j_t *wsj = (wts_screen_j_t *)result;
1517 	int i;
1518 
1519 	for (i = 0; i < WTS_CACHE_SIZE_X; i++)
1520 	    wsj->xcache[i].tag = -1;
1521 	for (i = 0; i < WTS_CACHE_SIZE_Y; i++)
1522 	    wsj->ycache[i].tag = -1;
1523     }
1524 #endif
1525     { /* 20090929:mu: memcpy(result->samples, buf + hdr_size, cell_size); */
1526 	/* assuming wts_screen_sample_t == ushort */
1527 	wts_screen_sample_t *p = result->samples, *q = (wts_screen_sample_t *)(buf + hdr_size);
1528 	int i = result->cell_width * result->cell_height;
1529 	for ( ; i > 0; i-- ) {
1530 		*p++ = le_u16(q);
1531 		q++;
1532 	}
1533     }
1534 
1535     return result;
1536 }
1537 
1538 #if 0 /* Never called. */
1539 /* Return value is size of buf in bytes */
1540 static int
1541 gs_wts_to_buf(const wts_screen_t *ws, byte **pbuf)
1542 {
1543     int size = wts_size(ws);
1544     int cell_size = ws->cell_width * ws->cell_height * sizeof(wts_screen_sample_t);
1545     byte *buf;
1546 
1547     buf = (byte *)malloc(size + cell_size);
1548     if (buf == NULL)
1549 	return -1;
1550     memcpy(buf, ws, size);
1551     ((wts_screen_t *)buf)->samples = NULL;
1552     memcpy(buf + size, ws->samples, cell_size);
1553     *pbuf = buf;
1554 
1555     return size + cell_size;
1556 }
1557 #endif
1558 
1559 #ifdef UNIT_TEST
1560 static int
dump_thresh(const wts_screen_t * ws,int width,int height)1561 dump_thresh(const wts_screen_t *ws, int width, int height)
1562 {
1563     int x, y;
1564     wts_screen_sample_t *s0;
1565     int cx, cy
1566     int dummy;
1567 
1568     wts_get_samples(ws, 0, 0, &cx, &cy, &dummy);
1569     s0 = ws->samples + cy * ws->cell_width + cx;
1570 
1571     printf("P5\n%d %d\n255\n", width, height);
1572     for (y = 0; y < height; y++) {
1573 	for (x = 0; x < width;) {
1574 	    wts_screen_sample_t *samples;
1575 	    int n_samples, i;
1576 
1577 	    wts_get_samples(ws, x, y, &cx, &cy, &n_samples);
1578 	    samples = ws->samples + cy * ws->cell_width + cx;
1579 #if 1
1580 	    for (i = 0; x + i < width && i < n_samples; i++)
1581 		fputc(samples[i] >> 7, stdout);
1582 #else
1583 	    printf("(%d, %d): %d samples at %d\n",
1584 		   x, y, n_samples, samples - s0);
1585 #endif
1586 	    x += n_samples;
1587 	}
1588     }
1589     return 0;
1590 }
1591 
1592 int
main(int argc,char ** argv)1593 main (int argc, char **argv)
1594 {
1595     gs_screen_halftone h;
1596     gs_matrix mat;
1597     double xres = 1200;
1598     double yres = 1200;
1599     gx_wts_cell_params_t *wcp;
1600     gs_wts_screen_enum_t *wse;
1601     wts_screen_t *ws;
1602 
1603     mat.xx = xres / 72.0;
1604     mat.xy = 0;
1605     mat.yx = 0;
1606     mat.yy = yres / 72.0;
1607 
1608     h.frequency = 121;
1609     h.angle = 45;
1610 
1611     wcp = wts_pick_cell_size(&h, &mat);
1612     dlprintf2("cell size = %d x %d\n", wcp->width, wcp->height);
1613     wse = gs_wts_screen_enum_new(wcp);
1614     wts_run_enum_squaredot(wse);
1615 #if 1
1616     wts_sort_blue(wse);
1617 #else
1618     wts_sort_cell(wse);
1619 #endif
1620     ws = wts_screen_from_enum(wse);
1621 
1622     dump_thresh(ws, 512, 512);
1623     return 0;
1624 }
1625 #endif
1626