12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337 |
- /* Copyright (C) 2002 artofcode LLC. All rights reserved.
-
- This software is provided AS-IS with no warranty, either express or
- implied.
-
- This software is distributed under license and may not be copied,
- modified or distributed except as expressly authorized under the terms
- of the license contained in the file LICENSE in this distribution.
-
- For more information about licensing, please refer to
- http://www.ghostscript.com/licensing/. For information on
- commercial licensing, go to http://www.artifex.com/licensing/ or
- contact Artifex Software, Inc., 101 Lucas Valley Road #110,
- San Rafael, CA 94903, U.S.A., +1(415)492-9861.
- */
- /*$Id: gswts.c,v 1.5 2003/08/26 15:38:50 igor Exp $ */
- /* Screen generation for Well Tempered Screening. */
- #include "stdpre.h"
- #include <stdlib.h> /* for malloc */
- #include "gx.h"
- #include "gxstate.h"
- #include "gsht.h"
- #include "math_.h"
- #include "gserrors.h"
- #include "gxfrac.h"
- #include "gxwts.h"
- #include "gswts.h"
- #define VERBOSE
- #ifdef UNIT_TEST
- /**
- * This file can be compiled independently for unit testing purposes.
- * Try this invocation:
- *
- * gcc -I../obj -DUNIT_TEST gswts.c gxwts.c -o test_gswts -lm
- * ./test_gswts | sed '/P5/,$!d' | xv -
- **/
- #undef printf
- #undef stdout
- #undef dlprintf1
- #define dlprintf1 printf
- #undef dlprintf2
- #define dlprintf2 printf
- #undef dlprintf3
- #define dlprintf3 printf
- #undef dlprintf4
- #define dlprintf4 printf
- #undef dlprintf5
- #define dlprintf5 printf
- #undef dlprintf6
- #define dlprintf6 printf
- #undef dlprintf7
- #define dlprintf7 printf
- #endif
- /* A datatype used for representing the product of two 32 bit integers.
- If a 64 bit integer type is available, it may be a better choice. */
- typedef double wts_bigint;
- typedef struct gx_wts_cell_params_j_s gx_wts_cell_params_j_t;
- typedef struct gx_wts_cell_params_h_s gx_wts_cell_params_h_t;
- struct gx_wts_cell_params_j_s {
- gx_wts_cell_params_t base;
- int shift;
- double ufast_a;
- double vfast_a;
- double uslow_a;
- double vslow_a;
- /* Probabilities of "jumps". A and B jumps can happen when moving
- one pixel to the right. C and D can happen when moving one pixel
- down. */
- double pa;
- double pb;
- double pc;
- double pd;
- int xa;
- int ya;
- int xb;
- int yb;
- int xc;
- int yc;
- int xd;
- int yd;
- };
- struct gx_wts_cell_params_h_s {
- gx_wts_cell_params_t base;
- /* This is the exact value that x1 and (width-x1) approximates. */
- double px;
- /* Ditto y1 and (height-y1). */
- double py;
- int x1;
- int y1;
- };
- #define WTS_EXTRA_SORT_BITS 9
- #define WTS_SORTED_MAX (((frac_1) << (WTS_EXTRA_SORT_BITS)) - 1)
- typedef struct {
- int u;
- int v;
- int k;
- int l;
- } wts_vec_t;
- private void
- wts_vec_set(wts_vec_t *wv, int u, int v, int k, int l)
- {
- wv->u = u;
- wv->v = v;
- wv->k = k;
- wv->l = l;
- }
- private wts_bigint
- wts_vec_smag(const wts_vec_t *a)
- {
- wts_bigint u = a->u;
- wts_bigint v = a->v;
- return u * u + v * v;
- }
- private void
- wts_vec_add(wts_vec_t *a, const wts_vec_t *b, const wts_vec_t *c)
- {
- a->u = b->u + c->u;
- a->v = b->v + c->v;
- a->k = b->k + c->k;
- a->l = b->l + c->l;
- }
- private void
- wts_vec_sub(wts_vec_t *a, const wts_vec_t *b, const wts_vec_t *c)
- {
- a->u = b->u - c->u;
- a->v = b->v - c->v;
- a->k = b->k - c->k;
- a->l = b->l - c->l;
- }
- /**
- * wts_vec_gcd3: Compute 3-way "gcd" of three vectors.
- * @a, @b, @c: Vectors.
- *
- * Compute pair of vectors satisfying three constraints:
- * They are integer combinations of the source vectors.
- * The source vectors are integer combinations of the results.
- * The magnitude of the vectors is minimized.
- *
- * The algorithm for this computation is quite reminiscent of the
- * classic Euclid GCD algorithm for integers.
- *
- * On return, @a and @b contain the result. @c is destroyed.
- **/
- private void
- wts_vec_gcd3(wts_vec_t *a, wts_vec_t *b, wts_vec_t *c)
- {
- wts_vec_t d, e;
- double ra, rb, rc, rd, re;
- wts_vec_set(&d, 0, 0, 0, 0);
- wts_vec_set(&e, 0, 0, 0, 0);
- for (;;) {
- ra = wts_vec_smag(a);
- rb = wts_vec_smag(b);
- rc = wts_vec_smag(c);
- wts_vec_sub(&d, a, b);
- wts_vec_add(&e, a, b);
- rd = wts_vec_smag(&d);
- re = wts_vec_smag(&e);
- if (re && re < rd) {
- d = e;
- rd = re;
- }
- if (rd && rd < ra && ra >= rb) *a = d;
- else if (rd < rb && rb > ra) *b = d;
- else {
- wts_vec_sub(&d, a, c);
- wts_vec_add(&e, a, c);
- rd = wts_vec_smag(&d);
- re = wts_vec_smag(&e);
- if (re < rd) {
- d = e;
- rd = re;
- }
- if (rd && rd < ra && ra >= rc) *a = d;
- else if (rd < rc && rc > ra) *c = d;
- else {
- wts_vec_sub(&d, b, c);
- wts_vec_add(&e, b, c);
- rd = wts_vec_smag(&d);
- re = wts_vec_smag(&e);
- if (re < rd) {
- d = e;
- rd = re;
- }
- if (rd && rd < rb && rb >= rc) *b = d;
- else if (rd < rc && rc > rb) *c = d;
- else
- break;
- }
- }
- }
- }
- private wts_bigint
- wts_vec_cross(const wts_vec_t *a, const wts_vec_t *b)
- {
- wts_bigint au = a->u;
- wts_bigint av = a->v;
- wts_bigint bu = b->u;
- wts_bigint bv = b->v;
- return au * bv - av * bu;
- }
- private void
- wts_vec_neg(wts_vec_t *a)
- {
- a->u = -a->u;
- a->v = -a->v;
- a->k = -a->k;
- a->l = -a->l;
- }
- /* compute k mod m */
- private void
- wts_vec_modk(wts_vec_t *a, int m)
- {
- while (a->k < 0) a->k += m;
- while (a->k >= m) a->k -= m;
- }
- /* Compute modulo in rational cell. */
- private void
- wts_vec_modkls(wts_vec_t *a, int m, int i, int s)
- {
- while (a->l < 0) {
- a->l += i;
- a->k -= s;
- }
- while (a->l >= i) {
- a->l -= i;
- a->k += s;
- }
- while (a->k < 0) a->k += m;
- while (a->k >= m) a->k -= m;
- }
- private void
- wts_set_mat(gx_wts_cell_params_t *wcp, double sratiox, double sratioy,
- double sangledeg)
- {
- double sangle = sangledeg * M_PI / 180;
- wcp->ufast = cos(sangle) / sratiox;
- wcp->vfast = -sin(sangle) / sratiox;
- wcp->uslow = sin(sangle) / sratioy;
- wcp->vslow = cos(sangle) / sratioy;
- }
- /**
- * Calculate Screen H cell sizes.
- **/
- private void
- wts_cell_calc_h(double inc, int *px1, int *pxwidth, double *pp1, double memw)
- {
- double minrep = pow(2, memw) * 50 / pow(2, 7.5);
- int m1 = 0, m2 = 0;
- double e1, e2;
- int uacc;
- e1 = 1e5;
- e2 = 1e5;
- for (uacc = (int)ceil(minrep * inc); uacc <= floor(2 * minrep * inc); uacc++) {
- int mt;
- double et;
- mt = (int)floor(uacc / inc + 1e-5);
- et = uacc / inc - mt + mt * 0.001;
- if (et < e1) {
- e1 = et;
- m1 = mt;
- }
- mt = (int)ceil(uacc / inc - 1e-5);
- et = mt - uacc / inc + mt * 0.001;
- if (et < e2) {
- e2 = et;
- m2 = mt;
- }
- }
- if (m1 == m2) {
- *px1 = m1;
- *pxwidth = m1;
- *pp1 = 1.0;
- } else {
- *px1 = m1;
- *pxwidth = m1 + m2;
- e1 = fabs(m1 * inc - floor(0.5 + m1 * inc));
- e2 = fabs(m2 * inc - floor(0.5 + m2 * inc));
- *pp1 = e2 / (e1 + e2);
- }
- }
- /* Implementation for Screen H. This is optimized for 0 and 45 degree
- rotations. */
- private gx_wts_cell_params_t *
- wts_pick_cell_size_h(double sratiox, double sratioy, double sangledeg,
- double memw)
- {
- gx_wts_cell_params_h_t *wcph;
- double xinc, yinc;
- wcph = malloc(sizeof(gx_wts_cell_params_h_t));
- if (wcph == NULL)
- return NULL;
- wcph->base.t = WTS_SCREEN_H;
- wts_set_mat(&wcph->base, sratiox, sratioy, sangledeg);
- xinc = fabs(wcph->base.ufast);
- if (xinc == 0)
- xinc = fabs(wcph->base.vfast);
- wts_cell_calc_h(xinc, &wcph->x1, &wcph->base.width, &wcph->px, memw);
- yinc = fabs(wcph->base.uslow);
- if (yinc == 0)
- yinc = fabs(wcph->base.vslow);
- wts_cell_calc_h(yinc, &wcph->y1, &wcph->base.height, &wcph->py, memw);
- return &wcph->base;
- }
- private double
- wts_qart(double r, double rbase, double p, double pbase)
- {
- if (p < 1e-5) {
- return ((r + rbase) * p);
- } else {
- return ((r + rbase) * (p + pbase) - rbase * pbase);
- }
- }
- #ifdef VERBOSE
- private void
- wts_print_j_jump(const gx_wts_cell_params_j_t *wcpj, const char *name,
- double pa, int xa, int ya)
- {
- dlprintf6("jump %s: (%d, %d) %f, actual (%f, %f)\n",
- name, xa, ya, pa,
- wcpj->ufast_a * xa + wcpj->uslow_a * ya,
- wcpj->vfast_a * xa + wcpj->vslow_a * ya);
- }
- private void
- wts_j_add_jump(const gx_wts_cell_params_j_t *wcpj, double *pu, double *pv,
- double pa, int xa, int ya)
- {
- double jump_u = wcpj->ufast_a * xa + wcpj->uslow_a * ya;
- double jump_v = wcpj->vfast_a * xa + wcpj->vslow_a * ya;
- jump_u -= floor(jump_u + 0.5);
- jump_v -= floor(jump_v + 0.5);
- *pu += pa * jump_u;
- *pv += pa * jump_v;
- }
- private void
- wts_print_j(const gx_wts_cell_params_j_t *wcpj)
- {
- double uf, vf;
- double us, vs;
- dlprintf3("cell = %d x %d, shift = %d\n",
- wcpj->base.width, wcpj->base.height, wcpj->shift);
- wts_print_j_jump(wcpj, "a", wcpj->pa, wcpj->xa, wcpj->ya);
- wts_print_j_jump(wcpj, "b", wcpj->pb, wcpj->xb, wcpj->yb);
- wts_print_j_jump(wcpj, "c", wcpj->pc, wcpj->xc, wcpj->yc);
- wts_print_j_jump(wcpj, "d", wcpj->pd, wcpj->xd, wcpj->yd);
- uf = wcpj->ufast_a;
- vf = wcpj->vfast_a;
- us = wcpj->uslow_a;
- vs = wcpj->vslow_a;
- wts_j_add_jump(wcpj, &uf, &vf, wcpj->pa, wcpj->xa, wcpj->ya);
- wts_j_add_jump(wcpj, &uf, &vf, wcpj->pb, wcpj->xb, wcpj->yb);
- wts_j_add_jump(wcpj, &us, &vs, wcpj->pc, wcpj->xc, wcpj->yc);
- wts_j_add_jump(wcpj, &us, &vs, wcpj->pd, wcpj->xd, wcpj->yd);
- dlprintf6("d: %f, %f; a: %f, %f; err: %g, %g\n",
- wcpj->base.ufast, wcpj->base.vfast,
- wcpj->ufast_a, wcpj->vfast_a,
- wcpj->base.ufast - uf, wcpj->base.vfast - vf);
- dlprintf6("d: %f, %f; a: %f, %f; err: %g, %g\n",
- wcpj->base.uslow, wcpj->base.vslow,
- wcpj->uslow_a, wcpj->vslow_a,
- wcpj->base.uslow - us, wcpj->base.vslow - vs);
- }
- #endif
- /**
- * wts_set_scr_jxi_try: Try a width for setting Screen J parameters.
- * @wcpj: Screen J parameters.
- * @m: Width to try.
- * @qb: Best quality score seen so far.
- * @jmem: Weight given to memory usage.
- *
- * Evaluate the quality for width @i. If the quality is better than
- * @qb, then set the rest of the parameters in @wcpj.
- *
- * This routine corresponds to SetScrJXITrySP in the WTS source.
- *
- * Return value: Quality score for parameters chosen.
- **/
- private double
- wts_set_scr_jxi_try(gx_wts_cell_params_j_t *wcpj, int m, double qb,
- double jmem)
- {
- const double uf = wcpj->base.ufast;
- const double vf = wcpj->base.vfast;
- const double us = wcpj->base.uslow;
- const double vs = wcpj->base.vslow;
- wts_vec_t a, b, c;
- double ufj, vfj;
- const double rbase = 0.1;
- const double pbase = 0.001;
- double ug, vg;
- double pp, pa, pb;
- double rf;
- double xa, ya, ra, xb, yb, rb;
- double q, qx, qw, qxl, qbi;
- double pya, pyb;
- int i;
- bool jumpok;
- wts_vec_set(&a, (int)floor(uf * m + 0.5), (int)floor(vf * m + 0.5), 1, 0);
- if (a.u == 0 && a.v == 0)
- return qb + 1;
-
- ufj = a.u / (double)m;
- vfj = a.v / (double)m;
- /* (ufj, vfj) = movement in UV space from (0, 1) in XY space */
- wts_vec_set(&b, m, 0, 0, 0);
- wts_vec_set(&c, 0, m, 0, 0);
- wts_vec_gcd3(&a, &b, &c);
- ug = (uf - ufj) * m;
- vg = (vf - vfj) * m;
- pp = 1.0 / wts_vec_cross(&b, &a);
- pa = (b.u * vg - ug * b.v) * pp;
- pb = (ug * a.v - a.u * vg) * pp;
- if (pa < 0) {
- wts_vec_neg(&a);
- pa = -pa;
- }
- if (pb < 0) {
- wts_vec_neg(&b);
- pb = -pb;
- }
- wts_vec_modk(&a, m);
- wts_vec_modk(&b, m);
- /* Prefer nontrivial jumps. */
- jumpok = (wcpj->xa == 0 && wcpj->ya == 0) ||
- (wcpj->xb == 0 && wcpj->yb == 0) ||
- (a.k != 0 && b.k != 0);
- rf = (uf * uf + vf * vf) * m;
- xa = (a.u * uf + a.v * vf) / rf;
- ya = (a.v * uf - a.u * vf) / rf;
- xb = (b.u * uf + b.v * vf) / rf;
- yb = (b.v * uf - b.u * vf) / rf;
- ra = sqrt(xa * xa + 0.0625 * ya * ya);
- rb = sqrt(xb * xb + 0.0625 * yb * yb);
- qx = 0.5 * (wts_qart(ra, rbase, pa, pbase) +
- wts_qart(rb, rbase, pb, pbase));
- if (qx < 1.0 / 4000.0)
- qx *= 0.25;
- else
- qx -= 0.75 / 4000.0;
- if (m < 7500)
- qw = 0;
- else
- qw = 0.00025; /* cache penalty */
- qxl = qx + qw;
- if (qxl > qb)
- return qxl;
- /* width is ok, now try heights */
- pp = m / (double)wts_vec_cross(&b, &a);
- pya = (b.u * vs - us * b.v) * pp;
- pyb = (us * a.v - a.u * vs) * pp;
- qbi = qb;
- for (i = 1; qxl + m * i * jmem < qbi; i++) {
- int s = m * i;
- int ca, cb;
- double usj, vsj;
- double usg, vsg;
- wts_vec_t a1, b1, a2, b2;
- double pc, pd;
- int ck;
- double qy, qm;
- ca = (int)floor(i * pya + 0.5);
- cb = (int)floor(i * pyb + 0.5);
- wts_vec_set(&c, ca * a.u + cb * b.u, ca * a.v + cb * b.v, 0, 1);
- usj = c.u / (double)s;
- vsj = c.v / (double)s;
- usg = (us - usj);
- vsg = (vs - vsj);
- a1 = a;
- b1 = b;
- a1.u *= i;
- a1.v *= i;
- b1.u *= i;
- b1.v *= i;
- wts_vec_gcd3(&a1, &b1, &c);
- a2 = a1;
- b2 = b1;
- pp = s / (double)wts_vec_cross(&b1, &a1);
- pc = (b1.u * vsg - usg * b1.v) * pp;
- pd = (usg * a1.v - a1.u * vsg) * pp;
- if (pc < 0) {
- wts_vec_neg(&a1);
- pc = -pc;
- }
- if (pd < 0) {
- wts_vec_neg(&b1);
- pd = -pd;
- }
- ck = ca * a.k + cb * b.k;
- while (ck < 0) ck += m;
- while (ck >= m) ck -= m;
- wts_vec_modkls(&a1, m, i, ck);
- wts_vec_modkls(&b1, m, i, ck);
- rf = (us * us - vs * vs) * s;
- xa = (a1.u * us + a1.v * vs) / rf;
- ya = (a1.v * us - a1.u * vs) / rf;
- xb = (b1.u * us + b1.v * vs) / rf;
- yb = (b1.v * us - b1.u * vs) / rf;
- ra = sqrt(xa * xa + 0.0625 * ya * ya);
- rb = sqrt(xb * xb + 0.0625 * yb * yb);
- qy = 0.5 * (wts_qart(ra, rbase, pc, pbase) +
- wts_qart(rb, rbase, pd, pbase));
- if (qy < 1.0 / 100.0)
- qy *= 0.025;
- else
- qy -= 0.75 / 100.0;
- qm = s * jmem;
- /* first try a and b jumps within the scanline */
- q = qm + qw + qx + qy;
- if (q < qbi && jumpok) {
- #ifdef VERBOSE
- dlprintf7("m = %d, n = %d, q = %d, qx = %d, qy = %d, qm = %d, qw = %d\n",
- m, i, (int)(q * 1e6), (int)(qx * 1e6), (int)(qy * 1e6), (int)(qm * 1e6), (int)(qw * 1e6));
- #endif
- qbi = q;
- wcpj->base.width = m;
- wcpj->base.height = i;
- wcpj->shift = ck;
- wcpj->ufast_a = ufj;
- wcpj->vfast_a = vfj;
- wcpj->uslow_a = usj;
- wcpj->vslow_a = vsj;
- wcpj->xa = a.k;
- wcpj->ya = 0;
- wcpj->xb = b.k;
- wcpj->yb = 0;
- wcpj->xc = a1.k;
- wcpj->yc = a1.l;
- wcpj->xd = b1.k;
- wcpj->yd = b1.l;
- wcpj->pa = pa;
- wcpj->pb = pb;
- wcpj->pc = pc;
- wcpj->pd = pd;
- #ifdef VERBOSE
- wts_print_j(wcpj);
- #endif
- }
- /* then try unconstrained a and b jumps */
- if (i > 1) {
- double pa2, pb2, pp2;
- double qx2, qw2, q2;
- pp2 = pp;
- pa2 = (b2.u * vg - ug * b2.v) * pp2;
- pb2 = (ug * a2.v - a2.u * vg) * pp2;
- rf = (uf * uf + vf * vf) * s;
- xa = (a2.u * uf + a2.v * vf) / rf;
- ya = (a2.v * uf - a2.u * vf) / rf;
- xb = (b2.u * uf + b2.v * vf) / rf;
- yb = (b2.v * uf - b2.u * vf) / rf;
- ra = sqrt(xa * xa + 0.0625 * ya * ya);
- rb = sqrt(xb * xb + 0.0625 * yb * yb);
- if (pa2 < 0) {
- pa2 = -pa2;
- wts_vec_neg(&a2);
- }
- if (pb2 < 0) {
- pb2 = -pb2;
- wts_vec_neg(&b2);
- }
- wts_vec_modkls(&a2, m, i, ck);
- wts_vec_modkls(&b2, m, i, ck);
- qx2 = 0.5 * (wts_qart(ra, rbase, pa2, pbase) +
- wts_qart(rb, rbase, pb2, pbase));
- if (qx2 < 1.0 / 4000.0)
- qx2 *= 0.25;
- else
- qx2 -= 0.75 / 4000.0;
- if (s < 7500)
- qw2 = 0;
- else
- qw2 = 0.00025; /* cache penalty */
- q2 = qm + qw2 + qx2 + qy;
- if (q2 < qbi) {
- #ifdef VERBOSE
- dlprintf7("m = %d, n = %d, q = %d, qx2 = %d, qy = %d, qm = %d, qw2 = %d\n",
- m, i, (int)(q * 1e6), (int)(qx * 1e6), (int)(qy * 1e6), (int)(qm * 1e6), (int)(qw2 * 1e6));
- #endif
- if (qxl > qw2 + qx2)
- qxl = qw2 + qx2;
- qbi = q2;
- wcpj->base.width = m;
- wcpj->base.height = i;
- wcpj->shift = ck;
- wcpj->ufast_a = ufj;
- wcpj->vfast_a = vfj;
- wcpj->uslow_a = usj;
- wcpj->vslow_a = vsj;
- wcpj->xa = a2.k;
- wcpj->ya = a2.l;
- wcpj->xb = b2.k;
- wcpj->yb = a2.l;
- wcpj->xc = a1.k;
- wcpj->yc = a1.l;
- wcpj->xd = b1.k;
- wcpj->yd = b1.l;
- wcpj->pa = pa2;
- wcpj->pb = pb2;
- wcpj->pc = pc;
- wcpj->pd = pd;
- #ifdef VERBOSE
- wts_print_j(wcpj);
- #endif
- }
- } /* if (i > 1) */
- if (qx > 10 * qy)
- break;
- }
- return qbi;
- }
- private int
- wts_double_to_int_cap(double d)
- {
- if (d > 0x7fffffff)
- return 0x7fffffff;
- else return (int)d;
- }
- /**
- * wts_set_scr_jxi: Choose Screen J parameters.
- * @wcpj: Screen J parameters.
- * @jmem: Weight given to memory usage.
- *
- * Chooses a cell size based on the [uv]{fast,slow} parameters,
- * setting the rest of the parameters in @wcpj. Essentially, this
- * algorithm iterates through all plausible widths for the cell.
- *
- * This routine corresponds to SetScrJXISP in the WTS source.
- *
- * Return value: Quality score for parameters chosen.
- **/
- private double
- wts_set_scr_jxi(gx_wts_cell_params_j_t *wcpj, double jmem)
- {
- int i, imax;
- double q, qb;
- wcpj->xa = 0;
- wcpj->ya = 0;
- wcpj->xb = 0;
- wcpj->yb = 0;
- wcpj->xc = 0;
- wcpj->yc = 0;
- wcpj->xd = 0;
- wcpj->yd = 0;
- qb = 1.0;
- imax = wts_double_to_int_cap(qb / jmem);
- for (i = 1; i <= imax; i++) {
- if (i > 1) {
- q = wts_set_scr_jxi_try(wcpj, i, qb, jmem);
- if (q < qb) {
- qb = q;
- imax = wts_double_to_int_cap(q / jmem);
- if (imax >= 7500)
- imax = wts_double_to_int_cap((q - 0.00025) / jmem);
- if (imax < 7500) {
- imax = 7500;
- }
- }
- }
- }
- return qb;
- }
- /* Implementation for Screen J. This is optimized for general angles. */
- private gx_wts_cell_params_t *
- wts_pick_cell_size_j(double sratiox, double sratioy, double sangledeg,
- double memw)
- {
- gx_wts_cell_params_j_t *wcpj;
- double code;
- wcpj = malloc(sizeof(gx_wts_cell_params_j_t));
- if (wcpj == NULL)
- return NULL;
- wcpj->base.t = WTS_SCREEN_J;
- wts_set_mat(&wcpj->base, sratiox, sratioy, sangledeg);
- code = wts_set_scr_jxi(wcpj, pow(0.1, memw));
- if (code < 0) {
- free(wcpj);
- return NULL;
- }
- return &wcpj->base;
- }
- /**
- * wts_pick_cell_size: Choose cell size for WTS screen.
- * @ph: The halftone parameters.
- * @pmat: Transformation from 1/72" Y-up coords to device coords.
- *
- * Return value: The WTS cell parameters, or NULL on error.
- **/
- gx_wts_cell_params_t *
- wts_pick_cell_size(gs_screen_halftone *ph, const gs_matrix *pmat)
- {
- gx_wts_cell_params_t *result;
- /* todo: deal with landscape and mirrored coordinate systems */
- double sangledeg = ph->angle;
- double sratiox = 72.0 * fabs(pmat->xx) / ph->frequency;
- double sratioy = 72.0 * fabs(pmat->yy) / ph->frequency;
- double octangle;
- double memw = 8.0;
- octangle = sangledeg;
- while (octangle >= 45.0) octangle -= 45.0;
- while (octangle < -45.0) octangle += 45.0;
- if (fabs(octangle) < 1e-4)
- result = wts_pick_cell_size_h(sratiox, sratioy, sangledeg, memw);
- else
- result = wts_pick_cell_size_j(sratiox, sratioy, sangledeg, memw);
- if (result != NULL) {
- ph->actual_frequency = ph->frequency;
- ph->actual_angle = ph->angle;
- }
- return result;
- }
- struct gs_wts_screen_enum_s {
- wts_screen_type t;
- bits32 *cell;
- int width;
- int height;
- int size;
- int idx;
- };
- typedef struct {
- gs_wts_screen_enum_t base;
- const gx_wts_cell_params_j_t *wcpj;
- } gs_wts_screen_enum_j_t;
- typedef struct {
- gs_wts_screen_enum_t base;
- const gx_wts_cell_params_h_t *wcph;
- /* an argument can be made that these should be in the params */
- double ufast1, vfast1;
- double ufast2, vfast2;
- double uslow1, vslow1;
- double uslow2, vslow2;
- } gs_wts_screen_enum_h_t;
- private gs_wts_screen_enum_t *
- gs_wts_screen_enum_j_new(gx_wts_cell_params_t *wcp)
- {
- gs_wts_screen_enum_j_t *wsej;
- wsej = malloc(sizeof(gs_wts_screen_enum_j_t));
- wsej->base.t = WTS_SCREEN_J;
- wsej->wcpj = (const gx_wts_cell_params_j_t *)wcp;
- wsej->base.width = wcp->width;
- wsej->base.height = wcp->height;
- wsej->base.size = wcp->width * wcp->height;
- wsej->base.cell = malloc(wsej->base.size * sizeof(wsej->base.cell[0]));
- wsej->base.idx = 0;
- return (gs_wts_screen_enum_t *)wsej;
- }
- private int
- gs_wts_screen_enum_j_currentpoint(gs_wts_screen_enum_t *self,
- gs_point *ppt)
- {
- gs_wts_screen_enum_j_t *z = (gs_wts_screen_enum_j_t *)self;
- const gx_wts_cell_params_j_t *wcpj = z->wcpj;
-
- int x, y;
- double u, v;
- if (z->base.idx == z->base.size) {
- return 1;
- }
- x = z->base.idx % wcpj->base.width;
- y = z->base.idx / wcpj->base.width;
- u = wcpj->ufast_a * x + wcpj->uslow_a * y;
- v = wcpj->vfast_a * x + wcpj->vslow_a * y;
- u -= floor(u);
- v -= floor(v);
- ppt->x = 2 * u - 1;
- ppt->y = 2 * v - 1;
- return 0;
- }
- private gs_wts_screen_enum_t *
- gs_wts_screen_enum_h_new(gx_wts_cell_params_t *wcp)
- {
- gs_wts_screen_enum_h_t *wseh;
- const gx_wts_cell_params_h_t *wcph = (const gx_wts_cell_params_h_t *)wcp;
- int x1 = wcph->x1;
- int x2 = wcp->width - x1;
- int y1 = wcph->y1;
- int y2 = wcp->height - y1;
- wseh = malloc(sizeof(gs_wts_screen_enum_h_t));
- wseh->base.t = WTS_SCREEN_H;
- wseh->wcph = wcph;
- wseh->base.width = wcp->width;
- wseh->base.height = wcp->height;
- wseh->base.size = wcp->width * wcp->height;
- wseh->base.cell = malloc(wseh->base.size * sizeof(wseh->base.cell[0]));
- wseh->base.idx = 0;
- wseh->ufast1 = floor(0.5 + wcp->ufast * x1) / x1;
- wseh->vfast1 = floor(0.5 + wcp->vfast * x1) / x1;
- if (x2 > 0) {
- wseh->ufast2 = floor(0.5 + wcp->ufast * x2) / x2;
- wseh->vfast2 = floor(0.5 + wcp->vfast * x2) / x2;
- }
- wseh->uslow1 = floor(0.5 + wcp->uslow * y1) / y1;
- wseh->vslow1 = floor(0.5 + wcp->vslow * y1) / y1;
- if (y2 > 0) {
- wseh->uslow2 = floor(0.5 + wcp->uslow * y2) / y2;
- wseh->vslow2 = floor(0.5 + wcp->vslow * y2) / y2;
- }
- return &wseh->base;
- }
- private int
- gs_wts_screen_enum_h_currentpoint(gs_wts_screen_enum_t *self,
- gs_point *ppt)
- {
- gs_wts_screen_enum_h_t *z = (gs_wts_screen_enum_h_t *)self;
- const gx_wts_cell_params_h_t *wcph = z->wcph;
-
- int x, y;
- double u, v;
- if (self->idx == self->size) {
- return 1;
- }
- x = self->idx % wcph->base.width;
- y = self->idx / wcph->base.width;
- if (x < wcph->x1) {
- u = z->ufast1 * x;
- v = z->vfast1 * x;
- } else {
- u = z->ufast2 * (x - wcph->x1);
- v = z->vfast2 * (x - wcph->x1);
- }
- if (y < wcph->y1) {
- u += z->uslow1 * y;
- v += z->vslow1 * y;
- } else {
- u += z->uslow2 * (y - wcph->y1);
- v += z->vslow2 * (y - wcph->y1);
- }
- u -= floor(u);
- v -= floor(v);
- ppt->x = 2 * u - 1;
- ppt->y = 2 * v - 1;
- return 0;
- }
- gs_wts_screen_enum_t *
- gs_wts_screen_enum_new(gx_wts_cell_params_t *wcp)
- {
- if (wcp->t == WTS_SCREEN_J)
- return gs_wts_screen_enum_j_new(wcp);
- else if (wcp->t == WTS_SCREEN_H)
- return gs_wts_screen_enum_h_new(wcp);
- else
- return NULL;
- }
- int
- gs_wts_screen_enum_currentpoint(gs_wts_screen_enum_t *wse, gs_point *ppt)
- {
- if (wse->t == WTS_SCREEN_J)
- return gs_wts_screen_enum_j_currentpoint(wse, ppt);
- if (wse->t == WTS_SCREEN_H)
- return gs_wts_screen_enum_h_currentpoint(wse, ppt);
- else
- return -1;
- }
- int
- gs_wts_screen_enum_next(gs_wts_screen_enum_t *wse, floatp value)
- {
- bits32 sample;
- if (value < -1.0 || value > 1.0)
- return_error(gs_error_rangecheck);
- sample = (bits32) ((value + 1) * 0x7fffffff);
- wse->cell[wse->idx] = sample;
- wse->idx++;
- return 0;
- }
- /* Run the enum with a square dot. This is useful for testing. */
- #ifdef UNIT_TEST
- private void
- wts_run_enum_squaredot(gs_wts_screen_enum_t *wse)
- {
- int code;
- gs_point pt;
- floatp spot;
- for (;;) {
- code = gs_wts_screen_enum_currentpoint(wse, &pt);
- if (code)
- break;
- spot = 0.5 * (cos(pt.x * M_PI) + cos(pt.y * M_PI));
- gs_wts_screen_enum_next(wse, spot);
- }
- }
- #endif /* UNIT_TEST */
- private int
- wts_sample_cmp(const void *av, const void *bv) {
- const bits32 *const *a = (const bits32 *const *)av;
- const bits32 *const *b = (const bits32 *const *)bv;
- if (**a > **b) return 1;
- if (**a < **b) return -1;
- return 0;
- }
- /* This implementation simply sorts the threshold values (evening the
- distribution), without applying any moire reduction. */
- int
- wts_sort_cell(gs_wts_screen_enum_t *wse)
- {
- int size = wse->width * wse->height;
- bits32 *cell = wse->cell;
- bits32 **pcell;
- int i;
- pcell = malloc(size * sizeof(bits32 *));
- if (pcell == NULL)
- return -1;
- for (i = 0; i < size; i++)
- pcell[i] = &cell[i];
- qsort(pcell, size, sizeof(bits32 *), wts_sample_cmp);
- for (i = 0; i < size; i++)
- *pcell[i] = (bits32)floor(WTS_SORTED_MAX * (i + 0.5) / size);
- free(pcell);
- return 0;
- }
- /**
- * wts_blue_bump: Generate bump function for BlueDot.
- *
- * Return value: newly allocated bump.
- **/
- private bits32 *
- wts_blue_bump(gs_wts_screen_enum_t *wse)
- {
- const gx_wts_cell_params_t *wcp;
- int width = wse->width;
- int height = wse->height;
- int shift;
- int size = width * height;
- bits32 *bump;
- int i;
- double uf, vf;
- double am, eg;
- int z;
- int x0, y0;
- int x, y;
- if (wse->t == WTS_SCREEN_J) {
- gs_wts_screen_enum_j_t *wsej = (gs_wts_screen_enum_j_t *)wse;
- shift = wsej->wcpj->shift;
- wcp = &wsej->wcpj->base;
- } else if (wse->t == WTS_SCREEN_H) {
- gs_wts_screen_enum_h_t *wseh = (gs_wts_screen_enum_h_t *)wse;
- shift = 0;
- wcp = &wseh->wcph->base;
- } else
- return NULL;
- bump = (bits32 *)malloc(size * sizeof(bits32));
- if (bump == NULL)
- return NULL;
- for (i = 0; i < size; i++)
- bump[i] = 0;
- /* todo: more intelligence for anisotropic scaling */
- uf = wcp->ufast;
- vf = wcp->vfast;
- am = uf * uf + vf * vf;
- eg = (1 << 24) * 2.0 * sqrt (am);
- z = (int)(5 / sqrt (am));
- x0 = -(z / width) * shift - z;
- y0 = -(z % width);
- while (y0 < 0) {
- x0 -= shift;
- y0 += height;
- }
- while (x0 < 0) x0 += width;
- for (y = -z; y <= z; y++) {
- int x1 = x0;
- for (x = -z; x <= z; x++) {
- bump[y0 * width + x1] += (bits32)(eg * exp (-am * (x * x + y * y)));
- x1++;
- if (x1 == width)
- x1 = 0;
- }
- y0++;
- if (y0 == height) {
- x0 += shift;
- if (x0 >= width) x0 -= width;
- y0 = 0;
- }
- }
- return bump;
- }
- /**
- * wts_sort_blue: Sort cell using BlueDot.
- **/
- int
- wts_sort_blue(gs_wts_screen_enum_t *wse)
- {
- bits32 *cell = wse->cell;
- int width = wse->width;
- int height = wse->height;
- int shift;
- int size = width * height;
- bits32 *ref;
- bits32 **pcell;
- bits32 *bump;
- int i;
- if (wse->t == WTS_SCREEN_J) {
- gs_wts_screen_enum_j_t *wsej = (gs_wts_screen_enum_j_t *)wse;
- shift = wsej->wcpj->shift;
- } else
- shift = 0;
- ref = (bits32 *)malloc(size * sizeof(bits32));
- pcell = (bits32 **)malloc(size * sizeof(bits32 *));
- bump = wts_blue_bump(wse);
- if (ref == NULL || pcell == NULL || bump == NULL) {
- free(ref);
- free(pcell);
- free(bump);
- return -1;
- }
- for (i = 0; i < size; i++)
- pcell[i] = &cell[i];
- qsort(pcell, size, sizeof(bits32 *), wts_sample_cmp);
- /* set ref to sorted cell; pcell will now point to ref */
- for (i = 0; i < size; i++) {
- pcell[i] = (pcell[i] - cell) + ref;
- *pcell[i] = (bits32)floor((1 << 24) * (i + 0.5) / size);
- cell[i] = 0;
- }
- for (i = 0; i < size; i++) {
- bits32 gmin = *pcell[i];
- int j;
- int j_end = i + 5000;
- int jmin = i;
- int ix;
- int x0, y0;
- int x, y;
- int ref_ix, bump_ix;
- /* find minimum cell value, but prune search */
- if (j_end > size) j_end = size;
- for (j = i + 1; j < j_end; j++) {
- if (*pcell[j] < gmin) {
- gmin = *pcell[j];
- jmin = j;
- }
- }
- ix = pcell[jmin] - ref;
- pcell[jmin] = pcell[i];
- cell[ix] = (bits32)floor(WTS_SORTED_MAX * (i + 0.5) / size);
- x0 = ix % width;
- y0 = ix / width;
- /* Add in bump, centered at (x0, y0) */
- ref_ix = y0 * width;
- bump_ix = 0;
- for (y = 0; y < height; y++) {
- for (x = x0; x < width; x++)
- ref[ref_ix + x] += bump[bump_ix++];
- for (x = 0; x < x0; x++)
- ref[ref_ix + x] += bump[bump_ix++];
- ref_ix += width;
- y0++;
- if (y0 == height) {
- x0 += shift;
- if (x0 >= width) x0 -= width;
- y0 = 0;
- ref_ix = 0;
- }
- }
- /* Remove DC component to avoid integer overflow. */
- if ((i & 255) == 255 && i + 1 < size) {
- bits32 gmin = *pcell[i + 1];
- int j;
- for (j = i + 2; j < size; j++) {
- if (*pcell[j] < gmin) {
- gmin = *pcell[j];
- }
- }
- #ifdef VERBOSE
- if_debug1('h', "[h]gmin = %d\n", gmin);
- #endif
- for (j = i + 1; j < size; j++)
- *pcell[j] -= gmin;
-
- }
- }
- free(ref);
- free(pcell);
- free(bump);
- return 0;
- }
- private wts_screen_t *
- wts_screen_from_enum_j(const gs_wts_screen_enum_t *wse)
- {
- const gs_wts_screen_enum_j_t *wsej = (const gs_wts_screen_enum_j_t *)wse;
- wts_screen_j_t *wsj;
- wts_screen_sample_t *samples;
- int size;
- int i;
- wsj = malloc(sizeof(wts_screen_j_t));
- wsj->base.type = WTS_SCREEN_J;
- wsj->base.cell_width = wsej->base.width;
- wsj->base.cell_height = wsej->base.height;
- size = wsj->base.cell_width * wsj->base.cell_height;
- wsj->base.cell_shift = wsej->wcpj->shift;
- wsj->pa = (int)floor(wsej->wcpj->pa * (1 << 16) + 0.5);
- wsj->pb = (int)floor(wsej->wcpj->pb * (1 << 16) + 0.5);
- wsj->pc = (int)floor(wsej->wcpj->pc * (1 << 16) + 0.5);
- wsj->pd = (int)floor(wsej->wcpj->pd * (1 << 16) + 0.5);
- wsj->XA = wsej->wcpj->xa;
- wsj->YA = wsej->wcpj->ya;
- wsj->XB = wsej->wcpj->xb;
- wsj->YB = wsej->wcpj->yb;
- wsj->XC = wsej->wcpj->xc;
- wsj->YC = wsej->wcpj->yc;
- wsj->XD = wsej->wcpj->xd;
- wsj->YD = wsej->wcpj->yd;
- samples = malloc(sizeof(wts_screen_sample_t) * size);
- wsj->base.samples = samples;
- for (i = 0; i < size; i++) {
- samples[i] = wsej->base.cell[i] >> WTS_EXTRA_SORT_BITS;
- }
- return &wsj->base;
- }
- private wts_screen_t *
- wts_screen_from_enum_h(const gs_wts_screen_enum_t *wse)
- {
- const gs_wts_screen_enum_h_t *wseh = (const gs_wts_screen_enum_h_t *)wse;
- wts_screen_h_t *wsh;
- wts_screen_sample_t *samples;
- int size;
- int i;
- /* factor some of this out into a common init routine? */
- wsh = malloc(sizeof(wts_screen_h_t));
- wsh->base.type = WTS_SCREEN_H;
- wsh->base.cell_width = wseh->base.width;
- wsh->base.cell_height = wseh->base.height;
- size = wsh->base.cell_width * wsh->base.cell_height;
- wsh->base.cell_shift = 0;
- wsh->px = wseh->wcph->px;
- wsh->py = wseh->wcph->py;
- wsh->x1 = wseh->wcph->x1;
- wsh->y1 = wseh->wcph->y1;
- samples = malloc(sizeof(wts_screen_sample_t) * size);
- wsh->base.samples = samples;
- for (i = 0; i < size; i++) {
- samples[i] = wseh->base.cell[i] >> WTS_EXTRA_SORT_BITS;
- }
- return &wsh->base;
- }
- wts_screen_t *
- wts_screen_from_enum(const gs_wts_screen_enum_t *wse)
- {
- if (wse->t == WTS_SCREEN_J)
- return wts_screen_from_enum_j(wse);
- else if (wse->t == WTS_SCREEN_H)
- return wts_screen_from_enum_h(wse);
- else
- return NULL;
- }
- void
- gs_wts_free_enum(gs_wts_screen_enum_t *wse)
- {
- free(wse);
- }
- void
- gs_wts_free_screen(wts_screen_t * wts)
- {
- free(wts);
- }
- #ifdef UNIT_TEST
- private int
- dump_thresh(const wts_screen_t *ws, int width, int height)
- {
- int x, y;
- wts_screen_sample_t *s0;
- int dummy;
- wts_get_samples(ws, 0, 0, &s0, &dummy);
- printf("P5\n%d %d\n255\n", width, height);
- for (y = 0; y < height; y++) {
- for (x = 0; x < width;) {
- wts_screen_sample_t *samples;
- int n_samples, i;
- wts_get_samples(ws, x, y, &samples, &n_samples);
- #if 1
- for (i = 0; x + i < width && i < n_samples; i++)
- fputc(samples[i] >> 7, stdout);
- #else
- printf("(%d, %d): %d samples at %d\n",
- x, y, n_samples, samples - s0);
- #endif
- x += n_samples;
- }
- }
- return 0;
- }
- int
- main (int argc, char **argv)
- {
- gs_screen_halftone h;
- gs_matrix mat;
- double xres = 1200;
- double yres = 1200;
- gx_wts_cell_params_t *wcp;
- gs_wts_screen_enum_t *wse;
- wts_screen_t *ws;
- mat.xx = xres / 72.0;
- mat.xy = 0;
- mat.yx = 0;
- mat.yy = yres / 72.0;
- h.frequency = 121;
- h.angle = 45;
- wcp = wts_pick_cell_size(&h, &mat);
- dlprintf2("cell size = %d x %d\n", wcp->width, wcp->height);
- wse = gs_wts_screen_enum_new(wcp);
- wts_run_enum_squaredot(wse);
- #if 1
- wts_sort_blue(wse);
- #else
- wts_sort_cell(wse);
- #endif
- ws = wts_screen_from_enum(wse);
- dump_thresh(ws, 512, 512);
- return 0;
- }
- #endif
|