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