1 /*
2 * Testbed implementation of Even Better Screening.
3 * Updated to function properly by Robin Watts and Michael Vrhel
4 *
5 * Code in this module is covered by US Patents 5055942 and
6 * 5917614 and corresponding international patents. Contact
7 * Artifex for more information, or see:
8 *
9 * http://www.artifex.com/page/imaging-science.html
10 */
11
12 #include <stdio.h>
13 #include <stdlib.h>
14 #include <string.h>
15 #include <math.h>
16 #include "ets.h"
17
18 /* source for threshold matrix - need to improve build process */
19 #include "tm.h"
20
21 #define ETS_VERSION 150
22
23 #define ETS_SHIFT 16
24 #define IMO_SHIFT 14
25
26 #define FANCY_COUPLING
27
28 typedef struct _ETS_PlaneCtx ETS_PlaneCtx;
29 typedef unsigned int uint32;
30
31 struct _ETS_Ctx {
32 int width;
33 int n_planes;
34 int levels; /* Number of levels on output, <= 256 */
35 ETS_PlaneCtx **plane_ctx;
36 int aspect_x;
37 int aspect_y;
38 int *strengths;
39 int elo;
40 int ehi;
41 int *c_line;
42
43 int ets_style;
44 int r_style;
45
46 uint32 seed1;
47 uint32 seed2;
48
49 FILE *dump_file;
50 ETS_DumpLevel dump_level;
51
52 /* Threshold modulation array */
53 int y;
54 int tmwid;
55 int tmheight;
56 const signed char *tmmat;
57 };
58
59 struct _ETS_PlaneCtx {
60 int width; /* Width of plane in pixels */
61 int *dst_line; /* Output pointer */
62 int *err_line; /* Total error carried out of pixel in the line above */
63 int *r_line; /* expected distance table (see paper for details) */
64 int *a_line; /* expected distance intermediate table (see paper) */
65 int *b_line; /* expected distance intermediate table (see paper) */
66 int *lut; /* Table to map from input source value to internal
67 * intensity level. Internal intensity level is 0 to
68 * 1<<ETS_SHIFT. */
69 int *dist_lut; /* A table of "expected distance between set pixels"
70 * values, stored in fixed point format with (ETS_SHIFT-c1)
71 * fractional bits. Values outside of the 'level 0-1' band
72 * will be set to 0 to avoid ETS weighting being used. */
73 char *rs_lut; /* Random noise table; values between 0 and 24. x meaning
74 * use 32-x bits of random noise, */
75 int c1; /* Shift adjustment for the dist_lut. */
76 int tm_offset;/* Plane offset within tm data */
77 };
78
79 void *
ets_malloc_aligned(int size,int align)80 ets_malloc_aligned(int size, int align)
81 {
82 void *result;
83 void *alloced = malloc(size + align);
84 int pad;
85
86 if (alloced == 0)
87 return 0;
88 pad = ((12 - (int)alloced) & 15) + 4;
89 result = (void *)(pad + (char *)alloced);
90 ((int *)result)[-1] = pad;
91 return result;
92 }
93
94 void
ets_free_aligned(void * p)95 ets_free_aligned(void *p)
96 {
97 int pad = ((int *)p)[-1];
98 free((char*)p - pad);
99 }
100
101 static double
compute_distscale(const ETS_Params * params)102 compute_distscale(const ETS_Params *params)
103 {
104 double distscale = params->distscale;
105
106 if (distscale == 0.0)
107 {
108 distscale = -1;
109 switch(params->aspect_x)
110 {
111 case 1:
112 switch(params->aspect_y)
113 {
114 case 1:
115 distscale = 0.95;
116 break;
117 case 2:
118 distscale = 1.8;
119 break;
120 case 3:
121 distscale = 2.4; /* FIXME */
122 break;
123 case 4:
124 distscale = 3.6;
125 break;
126 }
127 break;
128 case 2:
129 switch(params->aspect_y)
130 {
131 case 1:
132 distscale = 1.8;
133 break;
134 case 2:
135 break;
136 case 3:
137 distscale = 1.35; /* FIXME */
138 break;
139 case 4:
140 break;
141 }
142 break;
143 case 3:
144 switch(params->aspect_y)
145 {
146 case 1:
147 distscale = 2.4; /* FIXME */
148 break;
149 case 2:
150 distscale = 1.35; /* FIXME */
151 break;
152 case 3:
153 break;
154 case 4:
155 distscale = 0.675; /* FIXME */
156 break;
157 }
158 break;
159 case 4:
160 switch(params->aspect_y)
161 {
162 case 1:
163 distscale = 3.6;
164 break;
165 case 2:
166 break;
167 case 3:
168 distscale = 0.675; /* FIXME */
169 break;
170 case 4:
171 break;
172 }
173 break;
174 }
175 if (distscale == -1)
176 {
177 fprintf(stderr, "aspect ratio of %d:%d not supported\n",
178 params->aspect_x, params->aspect_y);
179 exit(1);
180 }
181 }
182 return distscale;
183 }
184
185 static int
ets_log2(int x)186 ets_log2(int x)
187 {
188 int y = 0;
189 int z;
190
191 for (z = x; z > 1; z = z >> 1)
192 y++;
193 return y;
194 }
195
196 static int
ets_log2up(int x)197 ets_log2up(int x)
198 {
199 return ets_log2(x-1)+1;
200 }
201
202 static int
compute_randshift(int nl,int rs_base,int levels)203 compute_randshift(int nl, int rs_base, int levels)
204 {
205 int rs = rs_base;
206
207 if ((nl > (90 << (ETS_SHIFT - 10)) &&
208 nl < (129 << (ETS_SHIFT - 10))) ||
209 (nl > (162 << (ETS_SHIFT - 10)) &&
210 nl < (180 << (ETS_SHIFT - 10))))
211 rs--;
212 else if (nl > (321 << (ETS_SHIFT - 10)) &&
213 nl < (361 << (ETS_SHIFT - 10)))
214 {
215 rs--;
216 if (nl > (331 << (ETS_SHIFT - 10)) &&
217 nl < (351 << (ETS_SHIFT - 10)))
218 rs--;
219 }
220 else if ((nl == (levels - 1) << ETS_SHIFT) &&
221 nl > (((levels - 1) << ETS_SHIFT) -
222 (1 << (ETS_SHIFT - 2))))
223 {
224 /* don't add randomness in extreme shadows */
225 }
226 else if ((nl > (3 << (ETS_SHIFT - 2))))
227 {
228 nl -= (nl + (1 << (ETS_SHIFT - 2))) & -(1 << (ETS_SHIFT - 1));
229 if (nl < 0) nl = -nl;
230 if (nl < (1 << (ETS_SHIFT - 4))) rs--;
231 if (nl < (1 << (ETS_SHIFT - 5))) rs--;
232 if (nl < (1 << (ETS_SHIFT - 6))) rs--;
233 }
234 else
235 {
236 if (nl < (3 << (ETS_SHIFT - 3))) nl += 1 << (ETS_SHIFT - 2);
237 nl = nl - (1 << (ETS_SHIFT - 1));
238 if (nl < 0) nl = -nl;
239 if (nl < (1 << (ETS_SHIFT - 4))) rs--;
240 if (nl < (1 << (ETS_SHIFT - 5))) rs--;
241 if (nl < (1 << (ETS_SHIFT - 6))) rs--;
242 }
243 return rs;
244 }
245
246 /* Maximum number of planes, but actually we want to dynamically
247 allocate all scratch buffers that depend on this. */
248 #define M 16
249
250 /**
251 * ets_line: Screen a line using Even TonenFS screeing.
252 * @ctx: An #EBPlaneCtx context.
253 * @dest: Array of destination buffers, 8 bpp pixels each.
254 * @src: Array of source buffer, ET_SrcPixel pixels each.
255 *
256 * Screens a single line using Even ToneFS screening.
257 **/
258 void
ets_line(ETS_Ctx * etc,uchar ** dest,const ETS_SrcPixel * const * src)259 ets_line(ETS_Ctx *etc, uchar **dest, const ETS_SrcPixel *const *src)
260 {
261 int a[M], b[M];
262 int e_1_0[M], e_m1_1[M], e_0_1[M], e_1_1[M];
263 int i;
264 int im;
265 int *pa, *pb, *piir, *pr;
266 int r[M], rg;
267 int xd;
268 uint32 seed1 = etc->seed1;
269 uint32 seed2 = etc->seed2;
270 uint32 sum;
271 int plane_idx;
272 int n_planes = etc->n_planes;
273 int levels = etc->levels;
274 #ifdef OLD_QUANT
275 int dith_mul = levels << 8;
276 #else
277 int dith_mul = (levels - 1) << 8;
278 #endif
279 int imo_mul = (1 << (ETS_SHIFT + IMO_SHIFT)) / (levels - 1);
280 int aspect_x2 = etc->aspect_x * etc->aspect_x;
281 int aspect_y2 = etc->aspect_y * etc->aspect_y;
282 int *strengths = etc->strengths;
283 int elo = etc->elo;
284 int ehi = etc->ehi;
285 int coupling;
286 int *c_line = etc->c_line;
287 int rand_shift;
288 const signed char *tmline = etc->tmmat + (etc->y % etc->tmheight) * etc->tmwid;
289 int tmmask = etc->tmwid - 1;
290
291 xd = etc->width;
292
293 /* Setup initial conditions for walking across the scanline. Because we
294 * are dealing with multiple planes, we have arrays of each variable,
295 * indexed by p = plane number.
296 * a[p] = 2x+1 (where x is the horizontal distance to the nearest set pixel)
297 * b[p] = 2y+1 (where y is the vertical distance to the nearest set pixel)
298 * r[p] = distance^2 to the nearest set pixel in this plane.
299 * e_0_1[p] = error from pixel above
300 * e_1_0[p] = error from pixel to the left
301 * e_m1_1[p]= error from pixel above right
302 * e_1_1[p] = error from pixel above left
303 * Because the source line is run length compressed, we need to run
304 * length decode it. We do this using the following:
305 * count[p] = number of pixels left in this run
306 * src_idx[p] = which run length entry we are currently on.
307 */
308
309 /* A potted recap of the distance calculations in the paper for easy
310 * reference.
311 * distance to last dot = SQR( (aspect_y * x)^2 + (aspect_x * y)^2 )
312 * r = distance^2 = (aspect_y * x)^2 + (aspect_x * y)^2
313 * = aspect_y^2 * x^2 + aspect_x^2 * y^2
314 * r_below - r = (aspect_x^2 * (y+1)^2) - (aspect_x^2 * y^2)
315 * = aspect_x^2 * ( (y+1)^2 - y^2 )
316 * = aspect_x^2 * ( 2y + 1 )
317 * r_under - r = (aspect_y^2 * (x+1)^2) - (aspect_y^2 * x^2)
318 * = aspect_y^2 * ( (x+1)^2 - x^2 )
319 * = aspect_y^2 * ( 2x + 1 )
320 * So, we keep:
321 * a = aspect_y^2 * (2x+1)
322 * b = aspect_x^2 * (2y+1)
323 * And we can then update r by adding either a or b at each stage.
324 */
325 for (plane_idx = 0; plane_idx < n_planes; plane_idx++)
326 {
327 ETS_PlaneCtx *ctx = etc->plane_ctx[plane_idx];
328
329 a[plane_idx] = aspect_y2; /* aspect_y^2 * (2x + 1) where x = 0 */
330 b[plane_idx] = aspect_x2; /* aspect_x^2 * (2y + 1) where y = 0 */
331 r[plane_idx] = 0;
332 e_0_1[plane_idx] = 0;
333 e_1_0[plane_idx] = 0;
334 e_m1_1[plane_idx] = ctx->err_line[0];
335 }
336
337 coupling = 0;
338
339 for (i = 0; i < xd; i++)
340 {
341 #ifdef FANCY_COUPLING
342 coupling += c_line[i];
343 #else
344 coupling = 0;
345 #endif
346 /* Lookup image data and compute R for all planes. */
347 for (plane_idx = 0; plane_idx < n_planes; plane_idx++)
348 {
349 ETS_PlaneCtx *ctx = etc->plane_ctx[plane_idx];
350 ETS_SrcPixel src_pixel = src[plane_idx][i];
351 int new_r;
352 int c1 = ctx->c1;
353 int rlimit = 1 << (30 - ETS_SHIFT + c1);
354 uchar *dst_ptr = dest[plane_idx];
355 int new_e_1_0;
356 int achieved_error;
357 int err;
358 int imo;
359 int expected_r;
360
361 pr = ctx->r_line;
362 pa = ctx->a_line;
363 pb = ctx->b_line;
364 piir = ctx->err_line;
365
366 im = ctx->lut[src_pixel]; /* image pixel (ink level) */
367 expected_r = ctx->dist_lut[src_pixel]; /* expected distance */
368 rand_shift = ctx->rs_lut[src_pixel]; /* random noise shift */
369
370 /* Forward pass distance computation; equation 2 from paper */
371 if (r[plane_idx] + a[plane_idx] < pr[i])
372 {
373 r[plane_idx] += a[plane_idx];
374 a[plane_idx] += 2*aspect_y2;
375 }
376 else
377 {
378 a[plane_idx] = pa[i];
379 b[plane_idx] = pb[i];
380 r[plane_idx] = pr[i];
381 }
382
383 /* Shuffle all the errors and read the next one. */
384 e_1_1[plane_idx] = e_0_1[plane_idx];
385 e_0_1[plane_idx] = e_m1_1[plane_idx];
386 e_m1_1[plane_idx] = i == xd - 1 ? 0 : piir[i + 1];
387 /* Reuse of variables here; new_e_1_0 is the total error passed
388 * into this pixel, with the traditional fs weights. */
389 new_e_1_0 = ((e_1_0[plane_idx] * 7 + e_m1_1[plane_idx] * 3 +
390 e_0_1[plane_idx] * 5 + e_1_1[plane_idx] * 1) >> 4);
391
392 /* White pixels stay white */
393 if (im == 0)
394 {
395 dst_ptr[i] = 0;
396 /* If we are forcing white pixels to stay white, we should
397 * not propagate errors through them. Or at the very least
398 * we should attenuate such errors. */
399 new_e_1_0 = 0;
400 }
401 else
402 {
403 /* The guts of ets (Equation 5) */
404 int ets_bias;
405
406 if (expected_r == 0)
407 {
408 ets_bias = 0;
409 }
410 else
411 {
412 /* Read the current distance, and clamp to avoid overflow
413 * in subsequent calculations. */
414 new_r = r[plane_idx];
415 if (new_r > rlimit)
416 new_r = rlimit;
417 /* Should we store back with the limit? */
418
419 /* Change the units on the distance to match our lut
420 * and subtract our actual distance (rg) from the expected
421 * distance (expected_r). */
422 rg = new_r << (ETS_SHIFT - c1);
423 ets_bias = rg - expected_r;
424
425 /* So ets_bias is the difference that we want to base our
426 * threshold modulation on (section 2.1 of the paper).
427 * Exactly how do we do that? We present various options
428 * here.
429 * 0 no modulation
430 * 1 what the code did when it came to me. No reference
431 * to this in the paper.
432 * 2 use it unchanged.
433 * 3 like 1, but same shift either side of 0.
434 * 4+ scale the modulation down.
435 */
436 switch (etc->ets_style)
437 {
438 case 0:
439 ets_bias = 0;
440 break;
441 case 1:
442 if (ets_bias > 0) ets_bias >>= 3;
443 break;
444 case 2:
445 break;
446 case 3:
447 ets_bias >>= 3;
448 break;
449 default:
450 ets_bias /= etc->ets_style-3;
451 }
452 }
453
454 /* Non white pixels get biased, and have the error
455 * applied. The error starts from the total error passed
456 * in. */
457 err = new_e_1_0;
458
459 /* Plus any ETS bias (calculated above) */
460 err += ets_bias;
461
462 /* Plus any random noise. Again various options here:
463 * 0 No random noise
464 * 1 The code as it came to me, using lookup table
465 * 2 commented out when it came to me; using pseudo
466 * random numbers generated from seed.
467 */
468 switch(etc->r_style)
469 {
470 case 0:
471 break;
472 case 2:
473 /* Add the two seeds together */
474 sum = seed1 + seed2;
475
476 /* If the add generated a carry, increment
477 * the result of the addition.
478 */
479 if (sum < seed1 || sum < seed2) sum++;
480
481 /* Seed2 becomes old seed1, seed1 becomes result */
482 seed2 = seed1;
483 seed1 = sum;
484
485 err -= (sum >> rand_shift) - (0x80000000 >> rand_shift);
486 break;
487 case 1:
488 err += tmline[(i+ctx->tm_offset) & tmmask] << (24 - rand_shift);
489 break;
490 }
491
492 /* Clamp the error; this is explained in the paper in
493 * section 6 just after equation 7. */
494 /* FIXME: Understand this better */
495 if (err < elo)
496 err = elo;
497 else if (err > ehi)
498 err = ehi;
499
500 /* Add the coupling to our combined 'error + bias' value */
501 /* FIXME: Are we sure this shouldn't be clamped? */
502 err += coupling;
503
504 /* Calculate imo = the quantised image value (Equation 7) */
505 #ifdef OLD_QUANT
506 imo = ((err + im) * dith_mul) >> (ETS_SHIFT + 8);
507 #else
508 imo = ((err + im) * dith_mul + (1 << (ETS_SHIFT + 7))) >> (ETS_SHIFT + 8);
509 #endif
510 /* Clamp to allow for over/underflow due to large errors */
511 if (imo < 0) imo = 0;
512 else if (imo > levels - 1) imo = levels - 1;
513
514 /* Store final output pixel */
515 dst_ptr[i] = imo;
516
517 /* Calculate the error between the desired and the obtained
518 * pixel values. */
519 achieved_error = im - ((imo * imo_mul) >> IMO_SHIFT);
520
521 /* And the error passed in is updated with the error for
522 * this pixel. */
523 new_e_1_0 += achieved_error;
524
525 /* Do the magic coupling here; strengths is 0 when
526 * multiplane optimisation is turned off, hence coupling
527 * remains 0 always. Equation 6. */
528 coupling += (achieved_error * strengths[plane_idx]) >> 8;
529
530 /* If we output a set pixel, then reset our distances. */
531 if (imo != 0)
532 {
533 a[plane_idx] = aspect_y2;
534 b[plane_idx] = aspect_x2;
535 r[plane_idx] = 0;
536 }
537 }
538
539 /* Store the values back for the next pass (Equation 3) */
540 pa[i] = a[plane_idx];
541 pb[i] = b[plane_idx];
542 pr[i] = r[plane_idx];
543 piir[i] = new_e_1_0;
544 e_1_0[plane_idx] = new_e_1_0;
545 }
546 #ifdef FANCY_COUPLING
547 coupling = coupling >> 1;
548 c_line[i] = coupling;
549 #endif
550 }
551
552 /* Note: this isn't white optimized, but the payoff is probably not
553 that important. */
554 #ifdef FANCY_COUPLING
555 coupling = 0;
556 for (i = xd - 1; i >= 0; i--)
557 {
558 coupling = (coupling + c_line[i]) >> 1;
559 c_line[i] = (coupling - (coupling >> 4));
560 }
561 #endif
562
563 /* Update distances. Reverse scanline pass. */
564 for (plane_idx = 0; plane_idx < n_planes; plane_idx++)
565 {
566 ETS_PlaneCtx *ctx = etc->plane_ctx[plane_idx];
567 int av, bv, rv;
568 int c1 = ctx->c1;
569 int rlimit = 1 << (30 - ETS_SHIFT + c1);
570
571 pr = ctx->r_line;
572 pa = ctx->a_line;
573 pb = ctx->b_line;
574
575 av = aspect_y2;
576 bv = aspect_x2;
577 rv = 0;
578 for (i = xd - 1; i >= 0; i--)
579 {
580 /* Equation 4 from the paper */
581 if (rv + bv + av < pr[i] + pb[i])
582 {
583 rv += av;
584 av += (aspect_y2<<1);
585 }
586 else
587 {
588 rv = pr[i];
589 av = pa[i];
590 bv = pb[i];
591 }
592 if (rv > rlimit) rv = rlimit;
593 pa[i] = av;
594 pb[i] = bv + (aspect_x2 << 1);
595 pr[i] = rv + bv;
596 }
597 }
598
599 etc->seed1 = seed1;
600 etc->seed2 = seed2;
601 etc->y++;
602 }
603
604 /**
605 * ets_plane_free: Free an #EBPlaneCtx context.
606 * @ctx: The #EBPlaneCtx context to free.
607 *
608 * Frees @ctx.
609 **/
610 static void
ets_plane_free(ETS_PlaneCtx * ctx)611 ets_plane_free(ETS_PlaneCtx *ctx)
612 {
613 free(ctx->err_line);
614 free(ctx->r_line);
615 free(ctx->a_line);
616 free(ctx->b_line);
617 free(ctx->lut);
618 free(ctx->dist_lut);
619 free(ctx->rs_lut);
620 free(ctx);
621 }
622
623 /**
624 * ets_new: Create new Even ToneFS screening context.
625 * @source_width: Width of source buffer.
626 * @dest_width: Width of destination buffer, in pixels.
627 * @lut: Lookup table for gray values.
628 *
629 * Creates a new context for Even ToneFS screening.
630 *
631 * If @dest_width is larger than @source_width, then input lines will
632 * be expanded using nearest-neighbor sampling.
633 *
634 * @lut should be an array of 256 values, one for each possible input
635 * gray value. @lut is a lookup table for gray values. Output is from
636 * 0 for white (no ink) to ....
637 *
638 *
639 * Return value: The new #EBPlaneCtx context.
640 **/
641 static ETS_PlaneCtx *
ets_plane_new(const ETS_Params * params,ETS_Ctx * etc,int plane_idx)642 ets_plane_new(const ETS_Params *params, ETS_Ctx *etc, int plane_idx)
643 {
644 int width = params->width;
645 int *lut = params->luts[plane_idx];
646 ETS_PlaneCtx *result;
647 int i;
648 int *new_lut;
649 int *dist_lut;
650 char *rs_lut;
651 double distscale = compute_distscale(params);
652 int c1;
653 int rlimit;
654 int log2_levels, log2_aspect;
655 int rs_base;
656
657 result = (ETS_PlaneCtx *)malloc(sizeof(ETS_PlaneCtx));
658
659 result->width = width;
660
661 log2_levels = ets_log2(params->levels);
662 log2_aspect = ets_log2(params->aspect_x) + ets_log2(params->aspect_y); /* FIXME */
663 c1 = 6 + log2_aspect + log2_levels;
664 if (params->c1_scale)
665 c1 -= params->c1_scale[plane_idx];
666 result->c1 = c1;
667 rlimit = 1 << (30 - ETS_SHIFT + c1);
668 result->tm_offset = TM_WIDTH/ets_log2up(params->n_planes);
669
670 /* Set up a lut to map input values from the source domain to the
671 * amount of ink. Callers can provide a lut of their own, which can be
672 * used for gamma correction etc. In the absence of this, a linear
673 * distribution is assumed. The user supplied lut should map from
674 * 'amount of light' to 'gamma adjusted amount of light', as the code
675 * subtracts the final value from (1<<ETS_SHIFT) (typically 65536) to
676 * get 'amount of ink'. */
677 new_lut = (int *)malloc((ETS_SRC_MAX + 1) * sizeof(int));
678 for (i = 0; i < ETS_SRC_MAX + 1; i++)
679 {
680 int nli;
681
682 if (lut == NULL)
683 {
684 #if ETS_SRC_MAX == 255
685 nli = (i * 65793 + (i >> 7)) >> (24 - ETS_SHIFT);
686 #else
687 nli = (i * ((double) (1 << ETS_SHIFT)) / ETS_SRC_MAX) + 0.5;
688 #endif
689 }
690 else
691 nli = lut[i] >> (24 - ETS_SHIFT);
692 if (params->polarity == ETS_BLACK_IS_ZERO)
693 new_lut[i] = (1 << ETS_SHIFT) - nli;
694 else
695 new_lut[i] = nli;
696 }
697
698 /* Here we calculate 2 more lookup tables. These could be separated out
699 * into 2 different loops, but are done in 1 to avoid a small amount of
700 * recalculation.
701 * dist_lut[i] = expected distance between dots for a greyscale of level i
702 * rs_lut[i] = whacky random noise scale factor.
703 */
704 dist_lut = (int *)malloc((ETS_SRC_MAX + 1) * sizeof(int));
705 rs_lut = (char *)malloc((ETS_SRC_MAX + 1) * sizeof(int));
706
707 rs_base = 35 - ETS_SHIFT + log2_levels - params->rand_scale;
708
709 /* The paper says that the expected 'value' for a grayshade g is:
710 * d_avg = 0.95 / 0.95/(g^2)
711 * This seems wrong to me. Let's consider some common cases; for a given
712 * greyscale, lay out the 'ideal' dithering, then consider removing each
713 * set pixel in turn and measuring the distance between that pixel and
714 * the closest set pixel.
715 *
716 * g = 1/2 #.#.#.#. visibly, expected distance = SQR(2)
717 * .#.#.#.#
718 * #.#.#.#.
719 * .#.#.#.#
720 *
721 * g = 1/4 #.#.#.#. expected distance = 2
722 * ........
723 * #.#.#.#.
724 * ........
725 *
726 * g = 1/16 #...#... expected distance = 4
727 * ........
728 * ........
729 * ........
730 * #...#...
731 * ........
732 * ........
733 * ........
734 *
735 * This rough approach leads us to suspect that we should be finding
736 * values roughly proportional to 1/SQR(g). Given the algorithm works in
737 * terms of square distance, this means 1/g. This is at odds with the
738 * value given in the paper. Being charitable and assuming that the paper
739 * means 'squared distance' when it says 'value', we are still a square
740 * off.
741 *
742 * Nonetheless, the code as supplied uses 0.95/g for the squared distance
743 * (i.e. it appears to agree with our logic here).
744 */
745 for (i = 0; i <= ETS_SRC_MAX; i++)
746 {
747 double dist;
748 int nl = new_lut[i] * (params->levels - 1);
749 int rs;
750
751 /* This is (or is supposed to be) equation 5 from the paper. If nl
752 * is g, why aren't we dividing by nl*nl ? */
753 if (nl == 0)
754 {
755 /* The expected distance for an ink level of 0 is infinite. Just
756 * put 0! */
757 dist = 0;
758 }
759 else if (nl >= ((1<<ETS_SHIFT)/(params->levels-1)))
760 {
761 /* New from RJW: Our distance measurements are only meaningful
762 * within the bottom 'level band' of the output. Do not apply
763 * ETS to higher ink levels. */
764 dist = 0;
765 }
766 else
767 {
768 dist = (distscale * (1 << (2 * ETS_SHIFT - c1))) / nl;
769 if (dist > rlimit << (ETS_SHIFT - c1))
770 dist = rlimit << (ETS_SHIFT - c1);
771 }
772
773 if (params->rand_scale_luts == NULL)
774 {
775 rs = compute_randshift(nl, rs_base, params->levels);
776 rs_lut[i] = rs;
777 }
778 else
779 {
780 int val = params->rand_scale_luts[plane_idx][i];
781
782 rs_lut[i] = rs_base + 16 - ets_log2(val + (val >> 1));
783 }
784 dist_lut[i] = (int)dist;
785 }
786
787 result->lut = new_lut;
788 result->dist_lut = dist_lut;
789 result->rs_lut = rs_lut;
790
791 result->err_line = (int *)calloc(width, sizeof(int));
792 result->r_line = (int *)calloc(width, sizeof(int));
793 result->a_line = (int *)calloc(width, sizeof(int));
794 result->b_line = (int *)calloc(width, sizeof(int));
795 for (i = 0; i < width; i++)
796 {
797 result->a_line[i] = 1;
798 result->b_line[i] = 1;
799 /* Initialize error with a non zero random value to ensure dots don't
800 land on dots when we have same planes with same gray level and
801 the plane interaction option is turned off. Ideally the level
802 of this error should be based upon the values of the first line
803 to ensure that things get primed properly */
804 result->err_line[i] = -((rand () & 0x7fff) << 6) >> (24 - ETS_SHIFT);
805 }
806
807 return result;
808 }
809
810 ETS_Ctx *
ets_new(const ETS_Params * params)811 ets_new(const ETS_Params *params)
812 {
813 ETS_Ctx *result = (ETS_Ctx *)malloc(sizeof(ETS_Ctx));
814 int n_planes = params->n_planes;
815 int i;
816 int using_vectors = 0;
817
818 if (params->dump_file)
819 {
820 int header[5];
821
822 header[0] = 0x70644245;
823 header[1] = 'M' * 0x1010000 + 'I' * 0x101;
824 header[2] = ETS_VERSION;
825 header[3] = ETS_SRC_MAX;
826 header[4] = sizeof(ETS_SrcPixel);
827 fwrite(header, sizeof(int), sizeof(header) / sizeof(header[0]),
828 params->dump_file);
829 if (params->dump_level >= ETS_DUMP_PARAMS)
830 {
831 fwrite(params, 1, sizeof(ETS_Params), params->dump_file);
832 }
833 if (params->dump_level >= ETS_DUMP_LUTS)
834 {
835 for (i = 0; i < params->n_planes; i++)
836 fwrite(params->luts[i], sizeof(int), ETS_SRC_MAX + 1,
837 params->dump_file);
838 }
839 }
840
841 result->width = params->width;
842 result->n_planes = n_planes;
843 result->levels = params->levels;
844
845 result->aspect_x = params->aspect_x;
846 result->aspect_y = params->aspect_y;
847
848 result->ehi = (int)(0.6 * (1 << ETS_SHIFT) / (params->levels - 1));
849 result->elo = -result->ehi;
850
851 result->strengths = (int *)malloc(sizeof(int) * n_planes);
852 memcpy(result->strengths, params->strengths,
853 sizeof(int) * n_planes);
854
855 result->ets_style = params->ets_style;
856 result->r_style = params->r_style;
857
858 result->c_line = (int *)calloc(params->width, sizeof(int));
859
860 result->seed1 = 0x5324879f;
861 result->seed2 = 0xb78d0945;
862
863 result->dump_file = params->dump_file;
864 result->dump_level = params->dump_level;
865
866 result->plane_ctx = (ETS_PlaneCtx **)malloc(sizeof(ETS_PlaneCtx *) * n_planes);
867 for (i = 0; i < n_planes; i++)
868 result->plane_ctx[i] = ets_plane_new(params, result, i);
869 result->y = 0;
870 result->tmmat = tmmat;
871 result->tmwid = TM_WIDTH;
872 result->tmheight = TM_HEIGHT;
873 return result;
874 }
875
876 /**
877 * even_better_free: Free an #EvenBetterCtx context.
878 * @ctx: The #EvenBetterCtx context to free.
879 *
880 * Frees @ctx.
881 **/
882 void
ets_free(ETS_Ctx * ctx)883 ets_free(ETS_Ctx *ctx)
884 {
885 int i;
886 int n_planes = ctx->n_planes;
887
888 if (ctx->dump_file)
889 fclose(ctx->dump_file);
890
891 for (i = 0; i < n_planes; i++)
892 ets_plane_free(ctx->plane_ctx[i]);
893 free(ctx->plane_ctx);
894 free(ctx->strengths);
895 free(ctx->c_line);
896
897 free(ctx);
898 }
899