- start of external Imager API access:
[imager.git] / quant.c
CommitLineData
02d1d628
AMH
1/* quant.c - provides general image quantization
2 currently only used by gif.c, but maybe we'll support producing
3 8-bit (or bigger indexed) png files at some point
4*/
92bda632 5#include "imager.h"
02d1d628
AMH
6
7static void makemap_addi(i_quantize *, i_img **imgs, int count);
97c4effc 8static void makemap_mediancut(i_quantize *, i_img **imgs, int count);
02d1d628
AMH
9
10static
11void
12setcol(i_color *cl,unsigned char r,unsigned char g,unsigned char b,unsigned char a) {
13 cl->rgba.r=r;
14 cl->rgba.g=g;
15 cl->rgba.b=b;
16 cl->rgba.a=a;
17}
18
19
20
21/* make a colour map overwrites mc_existing/mc_count in quant Note
22 that i_makemap will be called once for each image if mc_perimage is
23 set and the format support multiple colour maps per image.
24
25 This means we don't need any special processing at this level to
26 handle multiple colour maps.
27*/
28
92bda632
TC
29/*
30=item i_quant_makemap(quant, imgs, count)
31
32=category Image quantization
33
34Analyzes the I<count> images in I<imgs> according to the rules in
35I<quant> to build a color map (optimal or not depending on
36quant->make_colors).
37
38=cut
39*/
40
02d1d628 41void
92bda632 42i_quant_makemap(i_quantize *quant, i_img **imgs, int count) {
97c4effc
TC
43
44 if (quant->translate == pt_giflib) {
45 /* giflib does it's own color table generation */
46 /* previously we used giflib's quantizer, but it didn't handle multiple
47 images, which made it hard to build a global color map
48 We've implemented our own median cut code so we can ignore
49 the giflib version */
50 makemap_mediancut(quant, imgs, count);
02d1d628 51 return;
97c4effc
TC
52 }
53
02d1d628
AMH
54 switch (quant->make_colors & mc_mask) {
55 case mc_none:
56 /* use user's specified map */
57 break;
58 case mc_web_map:
59 {
60 int r, g, b;
61 int i = 0;
62 for (r = 0; r < 256; r+=0x33)
63 for (g = 0; g < 256; g+=0x33)
64 for (b = 0; b < 256; b += 0x33)
b74d74ac 65 setcol(quant->mc_colors+i++, r, g, b, 255);
02d1d628
AMH
66 quant->mc_count = i;
67 }
68 break;
69
97c4effc
TC
70 case mc_median_cut:
71 makemap_mediancut(quant, imgs, count);
72 break;
73
02d1d628
AMH
74 case mc_addi:
75 default:
76 makemap_addi(quant, imgs, count);
77 break;
78 }
79}
80
02d1d628
AMH
81static void translate_closest(i_quantize *, i_img *, i_palidx *);
82static void translate_errdiff(i_quantize *, i_img *, i_palidx *);
83static void translate_addi(i_quantize *, i_img *, i_palidx *);
84
92bda632
TC
85/*
86=item i_quant_translate(quant, img)
87
88=category Image quantization
89
90Quantize the image given the palette in quant.
91
92On success returns a pointer to a memory block of img->xsize *
93img->ysize i_palidx entries.
94
95On failure returns NULL.
02d1d628 96
92bda632
TC
97You should call myfree() on the returned block when you're done with
98it.
99
100This function will fail if the supplied palette contains no colors.
101
102=cut
02d1d628 103*/
92bda632
TC
104i_palidx *
105i_quant_translate(i_quantize *quant, i_img *img) {
2ff8ed30 106 i_palidx *result;
f771d0ec
TC
107 int bytes;
108
a73aeb5f 109 mm_log((1, "quant_translate(quant %p, img %p)\n", quant, img));
2ff8ed30 110
1501d9b3
TC
111 /* there must be at least one color in the paletted (though even that
112 isn't very useful */
113 if (quant->mc_count == 0) {
114 i_push_error(0, "no colors available for translation");
115 return NULL;
116 }
117
f771d0ec
TC
118 bytes = img->xsize * img->ysize;
119 if (bytes / img->ysize != img->xsize) {
120 i_push_error(0, "integer overflow calculating memory allocation");
121 return NULL;
122 }
123 result = mymalloc(bytes);
2ff8ed30 124
02d1d628 125 switch (quant->translate) {
02d1d628 126 case pt_closest:
97c4effc 127 case pt_giflib:
02d1d628
AMH
128 translate_closest(quant, img, result);
129 break;
a73aeb5f 130
02d1d628
AMH
131 case pt_errdiff:
132 translate_errdiff(quant, img, result);
133 break;
a73aeb5f 134
02d1d628
AMH
135 case pt_perturb:
136 default:
137 translate_addi(quant, img, result);
138 break;
139 }
a73aeb5f 140
02d1d628
AMH
141 return result;
142}
143
02d1d628
AMH
144static void translate_closest(i_quantize *quant, i_img *img, i_palidx *out) {
145 quant->perturb = 0;
146 translate_addi(quant, img, out);
147}
148
149#define PWR2(x) ((x)*(x))
150
151typedef int (*cmpfunc)(const void*, const void*);
152
153typedef struct {
154 unsigned char r,g,b;
36e67d0b
TC
155 char fixed;
156 char used;
02d1d628
AMH
157 int dr,dg,db;
158 int cdist;
159 int mcount;
160} cvec;
161
162typedef struct {
163 int cnt;
164 int vec[256];
165} hashbox;
166
167typedef struct {
168 int boxnum;
169 int pixcnt;
170 int cand;
171 int pdc;
172} pbox;
173
18accb2a 174static void prescan(i_img **im,int count, int cnum, cvec *clr, i_sample_t *line);
02d1d628
AMH
175static void reorder(pbox prescan[512]);
176static int pboxcmp(const pbox *a,const pbox *b);
177static void boxcenter(int box,cvec *cv);
178static float frandn(void);
179static void boxrand(int box,cvec *cv);
180static void bbox(int box,int *r0,int *r1,int *g0,int *g1,int *b0,int *b1);
181static void cr_hashindex(cvec clr[256],int cnum,hashbox hb[512]);
182static int mindist(int boxnum,cvec *cv);
183static int maxdist(int boxnum,cvec *cv);
184
185/* Some of the simpler functions are kept here to aid the compiler -
186 maybe some of them will be inlined. */
187
188static int
189pixbox(i_color *ic) { return ((ic->channel[0] & 224)<<1)+ ((ic->channel[1]&224)>>2) + ((ic->channel[2] &224) >> 5); }
190
18accb2a
TC
191static int
192pixbox_ch(i_sample_t *chans) { return ((chans[0] & 224)<<1)+ ((chans[1]&224)>>2) + ((chans[2] &224) >> 5); }
193
02d1d628
AMH
194static unsigned char
195g_sat(int in) {
196 if (in>255) { return 255; }
197 else if (in>0) return in;
198 return 0;
199}
200
201static
202float
203frand(void) {
204 return rand()/(RAND_MAX+1.0);
205}
206
a659442a 207#ifdef NOTEF
02d1d628
AMH
208static
209int
210eucl_d(cvec* cv,i_color *cl) { return PWR2(cv->r-cl->channel[0])+PWR2(cv->g-cl->channel[1])+PWR2(cv->b-cl->channel[2]); }
a659442a 211#endif
02d1d628 212
18accb2a
TC
213static
214int
215eucl_d_ch(cvec* cv,i_sample_t *chans) {
216 return PWR2(cv->r - chans[0]) + PWR2(cv->g - chans[1])
217 + PWR2(cv->b - chans[2]);
218}
219
02d1d628
AMH
220static
221int
222ceucl_d(i_color *c1, i_color *c2) { return PWR2(c1->channel[0]-c2->channel[0])+PWR2(c1->channel[1]-c2->channel[1])+PWR2(c1->channel[2]-c2->channel[2]); }
223
18accb2a
TC
224static const int
225gray_samples[] = { 0, 0, 0 };
226
02d1d628
AMH
227/*
228
229This quantization algorithm and implementation routines are by Arnar
230M. Hrafnkelson. In case any new ideas are here they are mine since
231this was written from scratch.
232
233The algorithm uses local means in the following way:
234
235 For each point in the colormap we find which image points
236 have that point as it's closest point. We calculate the mean
237 of those points and in the next iteration it will be the new
238 entry in the colormap.
239
240In order to speed this process up (i.e. nearest neighbor problem) We
241divied the r,g,b space up in equally large 512 boxes. The boxes are
242numbered from 0 to 511. Their numbering is so that for a given vector
243it is known that it belongs to the box who is formed by concatenating the
2443 most significant bits from each component of the RGB triplet.
245
246For each box we find the list of points from the colormap who might be
247closest to any given point within the box. The exact solution
248involves finding the Voronoi map (or the dual the Delauny
249triangulation) and has many issues including numerical stability.
250
251So we use this approximation:
252
2531. Find which point has the shortest maximum distance to the box.
2542. Find all points that have a shorter minimum distance than that to the box
255
256This is a very simple task and is not computationally heavy if one
257takes into account that the minimum distances from a pixel to a box is
258always found by checking if it's inside the box or is closest to some
259side or a corner. Finding the maximum distance is also either a side
260or a corner.
261
262This approach results 2-3 times more than the actual points needed but
263is still a good gain over the complete space. Usually when one has a
264256 Colorcolor map a search over 30 is often obtained.
265
266A bit of an enhancement to this approach is to keep a seperate list
267for each side of the cube, but this will require even more memory.
268
269 Arnar M. Hrafnkelsson (addi@umich.edu);
270
271*/
272/*
273 Extracted from gifquant.c, removed dependencies on gif_lib,
274 and added support for multiple images.
275 starting from 1nov2000 by TonyC <tony@develop-help.com>.
276
277*/
278
279static void
280makemap_addi(i_quantize *quant, i_img **imgs, int count) {
281 cvec *clr;
36e67d0b 282 int cnum, i, x, y, bst_idx=0, ld, cd, iter, currhb, img_num;
18accb2a 283 i_sample_t *val;
02d1d628 284 float dlt, accerr;
9cfd5724 285 hashbox *hb;
18accb2a
TC
286 i_mempool mp;
287 int maxwidth = 0;
288 i_sample_t *line;
289 const int *sample_indices;
02d1d628 290
7ac6a2e9
TC
291 mm_log((1, "makemap_addi(quant %p { mc_count=%d, mc_colors=%p }, imgs %p, count %d)\n",
292 quant, quant->mc_count, quant->mc_colors, imgs, count));
293
18accb2a
TC
294 i_mempool_init(&mp);
295
296 clr = i_mempool_alloc(&mp, sizeof(cvec) * quant->mc_size);
297 hb = i_mempool_alloc(&mp, sizeof(hashbox) * 512);
02d1d628
AMH
298 for (i=0; i < quant->mc_count; ++i) {
299 clr[i].r = quant->mc_colors[i].rgb.r;
300 clr[i].g = quant->mc_colors[i].rgb.g;
301 clr[i].b = quant->mc_colors[i].rgb.b;
36e67d0b
TC
302 clr[i].fixed = 1;
303 clr[i].mcount = 0;
02d1d628
AMH
304 }
305 /* mymalloc doesn't clear memory, so I think we need this */
306 for (; i < quant->mc_size; ++i) {
7ac6a2e9
TC
307 /*clr[i].r = clr[i].g = clr[i].b = 0;*/
308 clr[i].dr = 0;
309 clr[i].dg = 0;
310 clr[i].db = 0;
36e67d0b
TC
311 clr[i].fixed = 0;
312 clr[i].mcount = 0;
02d1d628
AMH
313 }
314 cnum = quant->mc_size;
315 dlt = 1;
316
18accb2a
TC
317 for (img_num = 0; img_num < count; ++img_num) {
318 if (imgs[img_num]->xsize > maxwidth)
319 maxwidth = imgs[img_num]->xsize;
320 }
321 line = i_mempool_alloc(&mp, 3 * maxwidth * sizeof(*line));
322
323 prescan(imgs, count, cnum, clr, line);
02d1d628
AMH
324 cr_hashindex(clr, cnum, hb);
325
326 for(iter=0;iter<3;iter++) {
327 accerr=0.0;
328
36e67d0b
TC
329 for (img_num = 0; img_num < count; ++img_num) {
330 i_img *im = imgs[img_num];
18accb2a
TC
331 sample_indices = im->channels >= 3 ? NULL : gray_samples;
332 for(y=0;y<im->ysize;y++) {
333 i_gsamp(im, 0, im->xsize, y, line, sample_indices, 3);
334 val = line;
335 for(x=0;x<im->xsize;x++) {
336 ld=196608;
337 /*i_gpix(im,x,y,&val);*/
338 currhb=pixbox_ch(val);
339 /* printf("box = %d \n",currhb); */
340 for(i=0;i<hb[currhb].cnt;i++) {
341 /* printf("comparing: pix (%d,%d,%d) vec (%d,%d,%d)\n",val.channel[0],val.channel[1],val.channel[2],clr[hb[currhb].vec[i]].r,clr[hb[currhb].vec[i]].g,clr[hb[currhb].vec[i]].b); */
342
343 cd=eucl_d_ch(&clr[hb[currhb].vec[i]],val);
344 if (cd<ld) {
345 ld=cd; /* shortest distance yet */
346 bst_idx=hb[currhb].vec[i]; /* index of closest vector yet */
347 }
348 }
349
350 clr[bst_idx].mcount++;
351 accerr+=(ld);
352 clr[bst_idx].dr+=val[0];
353 clr[bst_idx].dg+=val[1];
354 clr[bst_idx].db+=val[2];
355
356 val += 3; /* next 3 samples (next pixel) */
357 }
02d1d628
AMH
358 }
359 }
18accb2a
TC
360
361 for(i=0;i<cnum;i++)
362 if (clr[i].mcount) {
363 clr[i].dr/=clr[i].mcount;
364 clr[i].dg/=clr[i].mcount;
365 clr[i].db/=clr[i].mcount;
366 }
367
02d1d628 368 /* for(i=0;i<cnum;i++) printf("vec(%d)=(%d,%d,%d) dest=(%d,%d,%d) matchcount=%d\n",
18accb2a
TC
369 i,clr[i].r,clr[i].g,clr[i].b,clr[i].dr,clr[i].dg,clr[i].db,clr[i].mcount); */
370
02d1d628 371 /* printf("total error: %.2f\n",sqrt(accerr)); */
18accb2a 372
02d1d628 373 for(i=0;i<cnum;i++) {
36e67d0b 374 if (clr[i].fixed) continue; /* skip reserved colors */
18accb2a 375
02d1d628 376 if (clr[i].mcount) {
18accb2a
TC
377 clr[i].used = 1;
378 clr[i].r=clr[i].r*(1-dlt)+dlt*clr[i].dr;
379 clr[i].g=clr[i].g*(1-dlt)+dlt*clr[i].dg;
380 clr[i].b=clr[i].b*(1-dlt)+dlt*clr[i].db;
02d1d628 381 } else {
18accb2a
TC
382 /* let's try something else */
383 clr[i].used = 0;
384 clr[i].r=rand();
385 clr[i].g=rand();
386 clr[i].b=rand();
02d1d628 387 }
18accb2a 388
02d1d628
AMH
389 clr[i].dr=0;
390 clr[i].dg=0;
391 clr[i].db=0;
392 clr[i].mcount=0;
393 }
394 cr_hashindex(clr,cnum,hb);
395 }
396
397
398#ifdef NOTEF
399 for(i=0;i<cnum;i++) {
400 cd=eucl_d(&clr[i],&val);
401 if (cd<ld) {
402 ld=cd;
403 bst_idx=i;
404 }
405 }
406#endif
407
36e67d0b
TC
408 /* if defined, we only include colours with an mcount or that were
409 supplied in the fixed palette, giving us a smaller output palette */
410#define ONLY_USE_USED
411#ifdef ONLY_USE_USED
412 /* transfer the colors back */
413 quant->mc_count = 0;
414 for (i = 0; i < cnum; ++i) {
415 if (clr[i].fixed || clr[i].used) {
416 /*printf("Adding %d (%d,%d,%d)\n", i, clr[i].r, clr[i].g, clr[i].b);*/
417 quant->mc_colors[quant->mc_count].rgb.r = clr[i].r;
418 quant->mc_colors[quant->mc_count].rgb.g = clr[i].g;
419 quant->mc_colors[quant->mc_count].rgb.b = clr[i].b;
420 ++quant->mc_count;
421 }
422 }
423#else
02d1d628
AMH
424 /* transfer the colors back */
425 for (i = 0; i < cnum; ++i) {
426 quant->mc_colors[i].rgb.r = clr[i].r;
427 quant->mc_colors[i].rgb.g = clr[i].g;
428 quant->mc_colors[i].rgb.b = clr[i].b;
429 }
430 quant->mc_count = cnum;
36e67d0b 431#endif
02d1d628 432
7ac6a2e9
TC
433#if 0
434 mm_log((1, "makemap_addi returns - quant.mc_count = %d\n", quant->mc_count));
435 for (i = 0; i < quant->mc_count; ++i)
436 mm_log((5, " map entry %d: (%d, %d, %d)\n", i, clr[i].r, clr[i].g, clr[i].b));
437#endif
438
18accb2a 439 i_mempool_destroy(&mp);
02d1d628
AMH
440}
441
97c4effc
TC
442typedef struct {
443 i_sample_t rgb[3];
444 int count;
445} quant_color_entry;
446
447#define MEDIAN_CUT_COLORS 32768
448
449#define MED_CUT_INDEX(c) ((((c).rgb.r & 0xF8) << 7) | \
450 (((c).rgb.g & 0xF8) << 2) | (((c).rgb.b & 0xF8) >> 3))
451
18accb2a
TC
452#define MED_CUT_GRAY_INDEX(c) ((((c).rgb.r & 0xF8) << 7) | \
453 (((c).rgb.r & 0xF8) << 2) | (((c).rgb.r & 0xF8) >> 3))
454
97c4effc
TC
455/* scale these to cover the whole range */
456#define MED_CUT_RED(index) ((((index) & 0x7C00) >> 10) * 255 / 31)
457#define MED_CUT_GREEN(index) ((((index) & 0x3E0) >> 5) * 255 / 31)
458#define MED_CUT_BLUE(index) (((index) & 0x1F) * 255 / 31)
459
460typedef struct {
461 i_sample_t min[3]; /* minimum for each channel */
462 i_sample_t max[3]; /* maximum for each channel */
463 i_sample_t width[3]; /* width for each channel */
464 int start, size; /* beginning and size of the partition */
465 int pixels; /* number of pixels represented by this partition */
466} medcut_partition;
467
468/*
469=item calc_part(part, colors)
470
471Calculates the new color limits for the given partition.
472
473Giflib assumes that the limits for the non-split channels stay the
474same, but this strikes me as incorrect, especially if the colors tend
475to be color ramps.
476
477Of course this could be optimized by not recalculating the channel we
478just sorted on, but it's not worth the effort right now.
479
480=cut
481*/
482static void calc_part(medcut_partition *part, quant_color_entry *colors) {
483 int i, ch;
484
485 for (ch = 0; ch < 3; ++ch) {
486 part->min[ch] = 255;
487 part->max[ch] = 0;
488 }
489 for (i = part->start; i < part->start + part->size; ++i) {
490 for (ch = 0; ch < 3; ++ch) {
491 if (part->min[ch] > colors[i].rgb[ch])
492 part->min[ch] = colors[i].rgb[ch];
493 if (part->max[ch] < colors[i].rgb[ch])
494 part->max[ch] = colors[i].rgb[ch];
495 }
496 }
497 for (ch = 0; ch < 3; ++ch) {
498 part->width[ch] = part->max[ch] - part->min[ch];
499 }
500}
501
502/* simple functions to sort by each channel - we could use a global, but
503 that would be bad */
504
505static int
506color_sort_red(void const *left, void const *right) {
507 return ((quant_color_entry *)left)->rgb[0] - ((quant_color_entry *)right)->rgb[0];
508}
509
510static int
511color_sort_green(void const *left, void const *right) {
512 return ((quant_color_entry *)left)->rgb[1] - ((quant_color_entry *)right)->rgb[1];
513}
514
515static int
516color_sort_blue(void const *left, void const *right) {
517 return ((quant_color_entry *)left)->rgb[2] - ((quant_color_entry *)right)->rgb[2];
518}
519
520static int (*sorters[])(void const *, void const *) =
521{
522 color_sort_red,
523 color_sort_green,
524 color_sort_blue,
525};
526
527static void
528makemap_mediancut(i_quantize *quant, i_img **imgs, int count) {
529 quant_color_entry *colors;
530 i_mempool mp;
531 int imgn, x, y, i, ch;
532 int max_width;
533 i_color *line;
534 int color_count;
535 int total_pixels;
536 medcut_partition *parts;
537 int part_num;
538 int in, out;
18accb2a
TC
539 /* number of channels we search for the best channel to partition
540 this isn't terribly efficient, but it should work */
541 int chan_count;
97c4effc
TC
542
543 /*printf("images %d pal size %d\n", count, quant->mc_size);*/
544
545 i_mempool_init(&mp);
546
547 colors = i_mempool_alloc(&mp, sizeof(*colors) * MEDIAN_CUT_COLORS);
548 for (i = 0; i < MEDIAN_CUT_COLORS; ++i) {
549 colors[i].rgb[0] = MED_CUT_RED(i);
550 colors[i].rgb[1] = MED_CUT_GREEN(i);
551 colors[i].rgb[2] = MED_CUT_BLUE(i);
552 colors[i].count = 0;
553 }
554
555 max_width = -1;
556 for (imgn = 0; imgn < count; ++imgn) {
557 if (imgs[imgn]->xsize > max_width)
558 max_width = imgs[imgn]->xsize;
559 }
560 line = i_mempool_alloc(&mp, sizeof(i_color) * max_width);
561
562 /* build the stats */
563 total_pixels = 0;
18accb2a 564 chan_count = 1; /* assume we just have grayscale */
97c4effc
TC
565 for (imgn = 0; imgn < count; ++imgn) {
566 total_pixels += imgs[imgn]->xsize * imgs[imgn]->ysize;
567 for (y = 0; y < imgs[imgn]->ysize; ++y) {
568 i_glin(imgs[imgn], 0, imgs[imgn]->xsize, y, line);
18accb2a
TC
569 if (imgs[imgn]->channels > 2) {
570 chan_count = 3;
571 for (x = 0; x < imgs[imgn]->xsize; ++x) {
572 ++colors[MED_CUT_INDEX(line[x])].count;
573 }
574 }
575 else {
576 /* a gray-scale image, just use the first channel */
577 for (x = 0; x < imgs[imgn]->xsize; ++x) {
578 ++colors[MED_CUT_GRAY_INDEX(line[x])].count;
579 }
97c4effc
TC
580 }
581 }
582 }
583
584 /* eliminate the empty colors */
585 out = 0;
586 for (in = 0; in < MEDIAN_CUT_COLORS; ++in) {
587 if (colors[in].count) {
588 colors[out++] = colors[in];
589 }
590 }
591 /*printf("out %d\n", out);
592
593 for (i = 0; i < out; ++i) {
594 if (colors[i].count) {
595 printf("%d: (%d,%d,%d) -> %d\n", i, colors[i].rgb[0], colors[i].rgb[1],
596 colors[i].rgb[2], colors[i].count);
597 }
598 }*/
599
600 if (out < quant->mc_size) {
601 /* just copy them into the color table */
602 for (i = 0; i < out; ++i) {
603 for (ch = 0; ch < 3; ++ch) {
604 quant->mc_colors[i].channel[ch] = colors[i].rgb[ch];
605 }
606 }
607 quant->mc_count = out;
608 }
609 else {
610 /* build the starting partition */
611 parts = i_mempool_alloc(&mp, sizeof(*parts) * quant->mc_size);
612 parts[0].start = 0;
613 parts[0].size = out;
614 parts[0].pixels = total_pixels;
615 calc_part(parts, colors);
616 color_count = 1;
617
618 while (color_count < quant->mc_size) {
619 int max_index, max_ch; /* index/channel with biggest spread */
620 int max_size;
621 medcut_partition *workpart;
622 int cum_total;
623 int half;
624
625 /* find the partition with the most biggest span with more than
626 one color */
627 max_size = -1;
628 for (i = 0; i < color_count; ++i) {
18accb2a 629 for (ch = 0; ch < chan_count; ++ch) {
97c4effc
TC
630 if (parts[i].width[ch] > max_size
631 && parts[i].size > 1) {
632 max_index = i;
633 max_ch = ch;
634 max_size = parts[i].width[ch];
635 }
636 }
637 }
638
639 /* nothing else we can split */
640 if (max_size == -1)
641 break;
642
643 workpart = parts+max_index;
644 /*printf("splitting partition %d (pixels %ld, start %d, size %d)\n", max_index, workpart->pixels, workpart->start, workpart->size);*/
645 qsort(colors + workpart->start, workpart->size, sizeof(*colors),
646 sorters[max_ch]);
647
648 /* find the median or something like it we need to make sure both
649 sides of the split have at least one color in them, so we don't
650 test at the first or last entry */
651 i = workpart->start;
652 cum_total = colors[i].count;
653 ++i;
654 half = workpart->pixels / 2;
655 while (i < workpart->start + workpart->size - 1
656 && cum_total < half) {
657 cum_total += colors[i++].count;
658 }
659 /*printf("Split at %d to make %d (half %ld, cumtotal %ld)\n", i, color_count, half, cum_total);*/
660
661 /* found the spot to split */
662 parts[color_count].start = i;
663 parts[color_count].size = workpart->start + workpart->size - i;
664 workpart->size = i - workpart->start;
665 parts[color_count].pixels = workpart->pixels - cum_total;
666 workpart->pixels = cum_total;
667
668 /* recalculate the limits */
669 calc_part(workpart, colors);
670 calc_part(parts+color_count, colors);
671 ++color_count;
672 }
673
674 /* fill in the color table - since we could still have partitions
675 that have more than one color, we need to average the colors */
676 for (part_num = 0; part_num < color_count; ++part_num) {
677 long sums[3];
678 medcut_partition *workpart;
679
680 workpart = parts+part_num;
681 for (ch = 0; ch < 3; ++ch)
682 sums[ch] = 0;
683
684 for (i = workpart->start; i < workpart->start + workpart->size; ++i) {
685 for (ch = 0; ch < 3; ++ch) {
686 sums[ch] += colors[i].rgb[ch] * colors[i].count;
687 }
688 }
689 for (ch = 0; ch < 3; ++ch) {
690 quant->mc_colors[part_num].channel[ch] = sums[ch] / workpart->pixels;
691 }
692 }
693 quant->mc_count = color_count;
694 }
695 /*printf("out %d colors\n", quant->mc_count);*/
696 i_mempool_destroy(&mp);
697}
698
02d1d628
AMH
699#define pboxjump 32
700
701/* Define one of the following 4 symbols to choose a colour search method
702 The idea is to try these out, including benchmarking, to see which
703 is fastest in a good spread of circumstances.
704 I'd expect IM_CFLINSEARCH to be fastest for very small palettes, and
705 IM_CFHASHBOX for large images with large palettes.
706
707 Some other possibilities include:
708 - search over entries sorted by luminance
709
710 Initially I was planning on testing using the macros and then
711 integrating the code directly into each function, but this means if
712 we find a bug at a late stage we will need to update N copies of
713 the same code. Also, keeping the code in the macros means that the
714 code in the translation functions is much more to the point,
715 there's no distracting colour search code to remove attention from
716 what makes _this_ translation function different. It may be
717 advisable to move the setup code into functions at some point, but
718 it should be possible to do this fairly transparently.
719
720 If IM_CF_COPTS is defined then CFLAGS must have an appropriate
721 definition.
722
723 Each option needs to define 4 macros:
724 CF_VARS - variables to define in the function
725 CF_SETUP - code to setup for the colour search, eg. allocating and
726 initializing lookup tables
727 CF_FIND - code that looks for the color in val and puts the best
728 matching index in bst_idx
729 CF_CLEANUP - code to clean up, eg. releasing memory
730*/
731#ifndef IM_CF_COPTS
732/*#define IM_CFLINSEARCH*/
733#define IM_CFHASHBOX
734/*#define IM_CFSORTCHAN*/
735/*#define IM_CFRAND2DIST*/
736#endif
737
738#ifdef IM_CFHASHBOX
739
740/* The original version I wrote for this used the sort.
741 If this is defined then we use a sort to extract the indices for
742 the hashbox */
743#define HB_SORT
744
745/* assume i is available */
9cfd5724 746#define CF_VARS hashbox *hb = mymalloc(sizeof(hashbox) * 512); \
02d1d628
AMH
747 int currhb; \
748 long ld, cd
749
750#ifdef HB_SORT
751
752static long *gdists; /* qsort is annoying */
753/* int might be smaller than long, so we need to do a real compare
754 rather than a subtraction*/
755static int distcomp(void const *a, void const *b) {
756 long ra = gdists[*(int const *)a];
757 long rb = gdists[*(int const *)b];
758 if (ra < rb)
759 return -1;
760 else if (ra > rb)
761 return 1;
762 else
763 return 0;
764}
765
766#endif
767
768/* for each hashbox build a list of colours that are in the hb or is closer
769 than other colours
770 This is pretty involved. The original gifquant generated the hashbox
771 as part of it's normal processing, but since the map generation is now
772 separated from the translation we need to do this on the spot.
773 Any optimizations, even if they don't produce perfect results would be
774 welcome.
775 */
776static void hbsetup(i_quantize *quant, hashbox *hb) {
a659442a 777 long *dists, mind, maxd;
02d1d628
AMH
778 int cr, cb, cg, hbnum, i;
779 i_color cenc;
780#ifdef HB_SORT
781 int *indices = mymalloc(quant->mc_count * sizeof(int));
782#endif
783
784 dists = mymalloc(quant->mc_count * sizeof(long));
785 for (cr = 0; cr < 8; ++cr) {
786 for (cg = 0; cg < 8; ++cg) {
787 for (cb = 0; cb < 8; ++cb) {
788 /* centre of the hashbox */
789 cenc.channel[0] = cr*pboxjump+pboxjump/2;
790 cenc.channel[1] = cg*pboxjump+pboxjump/2;
791 cenc.channel[2] = cb*pboxjump+pboxjump/2;
792 hbnum = pixbox(&cenc);
793 hb[hbnum].cnt = 0;
794 /* order indices in the order of distance from the hashbox */
795 for (i = 0; i < quant->mc_count; ++i) {
796#ifdef HB_SORT
797 indices[i] = i;
798#endif
799 dists[i] = ceucl_d(&cenc, quant->mc_colors+i);
800 }
801#ifdef HB_SORT
802 /* it should be possible to do this without a sort
803 but so far I'm too lazy */
804 gdists = dists;
805 qsort(indices, quant->mc_count, sizeof(int), distcomp);
806 /* any colors that can match are within mind+diagonal size of
807 a hashbox */
808 mind = dists[indices[0]];
809 i = 0;
810 maxd = (sqrt(mind)+pboxjump)*(sqrt(mind)+pboxjump);
811 while (i < quant->mc_count && dists[indices[i]] < maxd) {
812 hb[hbnum].vec[hb[hbnum].cnt++] = indices[i++];
813 }
814#else
815 /* work out the minimum */
816 mind = 256*256*3;
817 for (i = 0; i < quant->mc_count; ++i) {
818 if (dists[i] < mind) mind = dists[i];
819 }
820 /* transfer any colours that might be closest to a colour in
821 this hashbox */
822 maxd = (sqrt(mind)+pboxjump)*(sqrt(mind)+pboxjump);
823 for (i = 0; i < quant->mc_count; ++i) {
824 if (dists[i] < maxd)
825 hb[hbnum].vec[hb[hbnum].cnt++] = i;
826 }
827#endif
828 }
829 }
862b614c 830 }
02d1d628
AMH
831#ifdef HB_SORT
832 myfree(indices);
833#endif
834 myfree(dists) ;
835}
836#define CF_SETUP hbsetup(quant, hb)
837
838#define CF_FIND \
839 currhb = pixbox(&val); \
840 ld = 196608; \
841 for (i = 0; i < hb[currhb].cnt; ++i) { \
842 cd = ceucl_d(quant->mc_colors+hb[currhb].vec[i], &val); \
843 if (cd < ld) { ld = cd; bst_idx = hb[currhb].vec[i]; } \
844 }
845
9cfd5724 846#define CF_CLEANUP myfree(hb)
02d1d628
AMH
847
848#endif
849
850#ifdef IM_CFLINSEARCH
851/* as simple as it gets */
852#define CF_VARS long ld, cd
853#define CF_SETUP /* none needed */
854#define CF_FIND \
855 ld = 196608; \
856 for (i = 0; i < quant->mc_count; ++i) { \
857 cd = ceucl_d(quant->mc_colors+i, &val); \
858 if (cd < ld) { ld = cd; bst_idx = i; } \
859 }
860#define CF_CLEANUP
861#endif
862
863#ifdef IM_CFSORTCHAN
864static int gsortchan;
865static i_quantize *gquant;
866static int chansort(void const *a, void const *b) {
867 return gquant->mc_colors[*(int const *)a].channel[gsortchan] -
868 gquant->mc_colors[*(int const *)b].channel[gsortchan];
869}
870#define CF_VARS int *indices, sortchan, diff; \
871 long ld, cd; \
872 int vindex[256] /* where to find value i of chan */
873
874static void chansetup(i_img *img, i_quantize *quant, int *csortchan,
875 int *vindex, int **cindices) {
876 int *indices, sortchan, chan, i, chval;
877 int chanmins[MAXCHANNELS], chanmaxs[MAXCHANNELS], maxrange;
878
879 /* find the channel with the maximum range */
880 /* the maximum stddev would probably be better */
881 for (chan = 0; chan < img->channels; ++chan) {
882 chanmins[chan] = 256; chanmaxs[chan] = 0;
883 for (i = 0; i < quant->mc_count; ++i) {
884 if (quant->mc_colors[i].channel[chan] < chanmins[chan])
885 chanmins[chan] = quant->mc_colors[i].channel[chan];
886 if (quant->mc_colors[i].channel[chan] > chanmaxs[chan])
887 chanmaxs[chan] = quant->mc_colors[i].channel[chan];
888 }
889 }
890 maxrange = -1;
891 for (chan = 0; chan < img->channels; ++chan) {
892 if (chanmaxs[chan]-chanmins[chan] > maxrange) {
893 maxrange = chanmaxs[chan]-chanmins[chan];
894 sortchan = chan;
895 }
896 }
897 indices = mymalloc(quant->mc_count * sizeof(int)) ;
898 for (i = 0; i < quant->mc_count; ++i) {
899 indices[i] = i;
900 }
901 gsortchan = sortchan;
902 gquant = quant;
903 qsort(indices, quant->mc_count, sizeof(int), chansort) ;
904 /* now a lookup table to find entries faster */
905 for (chval=0, i=0; i < quant->mc_count; ++i) {
906 while (chval < 256 &&
907 chval < quant->mc_colors[indices[i]].channel[sortchan]) {
908 vindex[chval++] = i;
909 }
910 }
911 while (chval < 256) {
912 vindex[chval++] = quant->mc_count-1;
913 }
914 *csortchan = sortchan;
915 *cindices = indices;
916}
917
918#define CF_SETUP \
919 chansetup(img, quant, &sortchan, vindex, &indices)
920
921int chanfind(i_color val, i_quantize *quant, int *indices, int *vindex,
922 int sortchan) {
923 int i, bst_idx, diff, maxdiff;
924 long ld, cd;
925
926 i = vindex[val.channel[sortchan]];
927 bst_idx = indices[i];
928 ld = 196608;
929 diff = 0;
930 maxdiff = quant->mc_count;
931 while (diff < maxdiff) {
932 if (i+diff < quant->mc_count) {
933 cd = ceucl_d(&val, quant->mc_colors+indices[i+diff]);
934 if (cd < ld) {
935 bst_idx = indices[i+diff];
936 ld = cd;
937 maxdiff = sqrt(ld);
938 }
939 }
940 if (i-diff >= 0) {
941 cd = ceucl_d(&val, quant->mc_colors+indices[i-diff]);
942 if (cd < ld) {
943 bst_idx = indices[i-diff];
944 ld = cd;
945 maxdiff = sqrt(ld);
946 }
947 }
948 ++diff;
949 }
950
951 return bst_idx;
952}
953
954#define CF_FIND \
955 bst_idx = chanfind(val, quant, indices, vindex, sortchan)
956
957
958#define CF_CLEANUP myfree(indices)
959
960#endif
961
962#ifdef IM_CFRAND2DIST
963
964/* This is based on a method described by Addi in the #imager channel
965 on the 28/2/2001. I was about 1am Sydney time at the time, so I
966 wasn't at my most cogent. Well, that's my excuse :)
967
968<TonyC> what I have at the moment is: hashboxes, with optimum hash box
969filling; simple linear search; and a lookup in the widest channel
970(currently the channel with the maximum range)
971<Addi> There is one more way that might be simple to implement.
972<Addi> You want to hear?
973<TonyC> what's that?
974<purl> somebody said that was not true
975<Addi> For each of the colors in the palette start by creating a
976sorted list of the form:
977<Addi> [distance, color]
978<Addi> Where they are sorted by distance.
979<TonyC> distance to where?
980<Addi> Where the elements in the lists are the distances and colors of
981the other colors in the palette
982<TonyC> ok
983<Addi> So if you are at color 0
984<Addi> ok - now to search for the closest color when you are creating
985the final image is done like this:
986<Addi> a) pick a random color from the palette
987<Addi> b) calculate the distance to it
988<Addi> c) only check the vectors that are within double the distance
989in the list of the color you picked from the palette.
990<Addi> Does that seem logical?
991<Addi> Lets imagine that we only have grayscale to make an example:
992<Addi> Our palette has 1 4 10 20 as colors.
993<Addi> And we want to quantize the color 11
994<Addi> lets say we picked 10 randomly
995<Addi> the double distance is 2
996<Addi> since abs(10-11)*2 is 2
997<Addi> And the list at vector 10 is this:
998<Addi> [0, 10], [6 4], [9, 1], [10, 20]
999<Addi> so we look at the first one (but not the second one since 6 is
1000at a greater distance than 2.
1001<Addi> Any of that make sense?
1002<TonyC> yes, though are you suggesting another random jump to one of
1003the colours with the possible choices? or an exhaustive search?
1004<Addi> TonyC: It's possible to come up with a recursive/iterative
1005enhancement but this is the 'basic' version.
1006<Addi> Which would do an iterative search.
1007<Addi> You can come up with conditions where it pays to switch to a new one.
1008<Addi> And the 'random' start can be switched over to a small tree.
1009<Addi> So you would have a little index at the start.
1010<Addi> to get you into the general direction
1011<Addi> Perhaps just an 8 split.
1012<Addi> that is - split each dimension in half.
1013<TonyC> yep
1014<TonyC> I get the idea
1015<Addi> But this would seem to be a good approach in our case since we
1016usually have few codevectors.
1017<Addi> So we only need 256*256 entries in a table.
1018<Addi> We could even only index some of them that were deemed as good
1019candidates.
1020<TonyC> I was considering adding paletted output support for PNG and
1021TIFF at some point, which support 16-bit palettes
1022<Addi> ohh.
1023<Addi> 'darn' ;)
1024
1025
1026*/
1027
1028
1029typedef struct i_dists {
1030 int index;
1031 long dist;
1032} i_dists;
1033
1034#define CF_VARS \
1035 i_dists *dists;
1036
1037static int dists_sort(void const *a, void const *b) {
1038 return ((i_dists *)a)->dist - ((i_dists *)b)->dist;
1039}
1040
1041static void rand2dist_setup(i_quantize *quant, i_dists **cdists) {
1042 i_dists *dists =
1043 mymalloc(sizeof(i_dists)*quant->mc_count*quant->mc_count);
1044 int i, j;
1045 long cd;
1046 for (i = 0; i < quant->mc_count; ++i) {
1047 i_dists *ldists = dists + quant->mc_count * i;
1048 i_color val = quant->mc_colors[i];
1049 for (j = 0; j < quant->mc_count; ++j) {
1050 ldists[j].index = j;
1051 ldists[j].dist = ceucl_d(&val, quant->mc_colors+j);
1052 }
1053 qsort(ldists, quant->mc_count, sizeof(i_dists), dists_sort);
1054 }
1055 *cdists = dists;
1056}
1057
1058#define CF_SETUP \
1059 bst_idx = rand() % quant->mc_count; \
1060 rand2dist_setup(quant, &dists)
1061
1062static int rand2dist_find(i_color val, i_quantize *quant, i_dists *dists, int index) {
1063 i_dists *cdists;
1064 long cd, ld;
1065 long maxld;
1066 int i;
1067 int bst_idx;
1068
1069 cdists = dists + index * quant->mc_count;
1070 ld = 3 * 256 * 256;
1071 maxld = 8 * ceucl_d(&val, quant->mc_colors+index);
1072 for (i = 0; i < quant->mc_count && cdists[i].dist <= maxld; ++i) {
1073 cd = ceucl_d(&val, quant->mc_colors+cdists[i].index);
1074 if (cd < ld) {
1075 bst_idx = cdists[i].index;
1076 ld = cd;
1077 }
1078 }
1079 return bst_idx;
1080}
1081
1082#define CF_FIND bst_idx = rand2dist_find(val, quant, dists, bst_idx)
1083
1084#define CF_CLEANUP myfree(dists)
1085
1086
1087#endif
1088
1089static void translate_addi(i_quantize *quant, i_img *img, i_palidx *out) {
1090 int x, y, i, k, bst_idx;
1091 i_color val;
1092 int pixdev = quant->perturb;
1093 CF_VARS;
1094
1095 CF_SETUP;
1096
18accb2a
TC
1097 if (img->channels >= 3) {
1098 if (pixdev) {
1099 k=0;
1100 for(y=0;y<img->ysize;y++) for(x=0;x<img->xsize;x++) {
1101 i_gpix(img,x,y,&val);
1102 val.channel[0]=g_sat(val.channel[0]+(int)(pixdev*frandn()));
1103 val.channel[1]=g_sat(val.channel[1]+(int)(pixdev*frandn()));
1104 val.channel[2]=g_sat(val.channel[2]+(int)(pixdev*frandn()));
1105 CF_FIND;
1106 out[k++]=bst_idx;
1107 }
1108 } else {
1109 k=0;
1110 for(y=0;y<img->ysize;y++) for(x=0;x<img->xsize;x++) {
1111 i_gpix(img,x,y,&val);
1112 CF_FIND;
1113 out[k++]=bst_idx;
1114 }
02d1d628 1115 }
18accb2a
TC
1116 }
1117 else {
1118 if (pixdev) {
1119 k=0;
1120 for(y=0;y<img->ysize;y++) for(x=0;x<img->xsize;x++) {
1121 i_gpix(img,x,y,&val);
1122 val.channel[1] = val.channel[2] =
1123 val.channel[0]=g_sat(val.channel[0]+(int)(pixdev*frandn()));
1124 CF_FIND;
1125 out[k++]=bst_idx;
1126 }
1127 } else {
1128 k=0;
1129 for(y=0;y<img->ysize;y++) for(x=0;x<img->xsize;x++) {
1130 i_gpix(img,x,y,&val);
1131 val.channel[1] = val.channel[2] = val.channel[0];
1132 CF_FIND;
1133 out[k++]=bst_idx;
1134 }
02d1d628
AMH
1135 }
1136 }
1137 CF_CLEANUP;
1138}
1139
1140static int floyd_map[] =
1141{
1142 0, 0, 7,
1143 3, 5, 1
1144};
1145
1146static int jarvis_map[] =
1147{
1148 0, 0, 0, 7, 5,
1149 3, 5, 7, 5, 3,
1150 1, 3, 5, 3, 1
1151};
1152
1153static int stucki_map[] =
1154{
1155 0, 0, 0, 8, 4,
1156 2, 4, 8, 4, 2,
1157 1, 2, 4, 2, 1
1158};
1159
1160struct errdiff_map {
1161 int *map;
1162 int width, height, orig;
1163};
1164
1165static struct errdiff_map maps[] =
1166{
1167 { floyd_map, 3, 2, 1 },
1168 { jarvis_map, 5, 3, 2 },
1169 { stucki_map, 5, 3, 2 },
1170};
1171
1172typedef struct errdiff_tag {
1173 int r, g, b;
1174} errdiff_t;
1175
1176/* perform an error diffusion dither */
1177static
1178void
1179translate_errdiff(i_quantize *quant, i_img *img, i_palidx *out) {
1180 int *map;
1181 int mapw, maph, mapo;
1182 int i;
1183 errdiff_t *err;
1184 int errw;
1185 int difftotal;
1186 int x, y, dx, dy;
02d1d628
AMH
1187 int bst_idx;
1188 CF_VARS;
1189
1190 if ((quant->errdiff & ed_mask) == ed_custom) {
1191 map = quant->ed_map;
1192 mapw = quant->ed_width;
1193 maph = quant->ed_height;
1194 mapo = quant->ed_orig;
1195 }
1196 else {
1197 int index = quant->errdiff & ed_mask;
1198 if (index >= ed_custom) index = ed_floyd;
1199 map = maps[index].map;
1200 mapw = maps[index].width;
1201 maph = maps[index].height;
1202 mapo = maps[index].orig;
1203 }
1204
1205 errw = img->xsize+mapw;
1206 err = mymalloc(sizeof(*err) * maph * errw);
1207 /*errp = err+mapo;*/
1208 memset(err, 0, sizeof(*err) * maph * errw);
1209
1210 difftotal = 0;
1211 for (i = 0; i < maph * mapw; ++i)
1212 difftotal += map[i];
1213 /*printf("map:\n");
1214 for (dy = 0; dy < maph; ++dy) {
1215 for (dx = 0; dx < mapw; ++dx) {
1216 printf("%2d", map[dx+dy*mapw]);
1217 }
1218 putchar('\n');
1219 }*/
1220
1221 CF_SETUP;
1222
1223 for (y = 0; y < img->ysize; ++y) {
1224 for (x = 0; x < img->xsize; ++x) {
1225 i_color val;
1226 long ld, cd;
1227 errdiff_t perr;
1228 i_gpix(img, x, y, &val);
18accb2a
TC
1229 if (img->channels < 3) {
1230 val.channel[1] = val.channel[2] = val.channel[0];
1231 }
02d1d628
AMH
1232 perr = err[x+mapo];
1233 perr.r = perr.r < 0 ? -((-perr.r)/difftotal) : perr.r/difftotal;
1234 perr.g = perr.g < 0 ? -((-perr.g)/difftotal) : perr.g/difftotal;
1235 perr.b = perr.b < 0 ? -((-perr.b)/difftotal) : perr.b/difftotal;
1236 /*printf("x %3d y %3d in(%3d, %3d, %3d) di(%4d,%4d,%4d)\n", x, y, val.channel[0], val.channel[1], val.channel[2], perr.r, perr.g, perr.b);*/
1237 val.channel[0] = g_sat(val.channel[0]-perr.r);
1238 val.channel[1] = g_sat(val.channel[1]-perr.g);
1239 val.channel[2] = g_sat(val.channel[2]-perr.b);
1240 CF_FIND;
1241 /* save error */
1242 perr.r = quant->mc_colors[bst_idx].channel[0] - val.channel[0];
1243 perr.g = quant->mc_colors[bst_idx].channel[1] - val.channel[1];
1244 perr.b = quant->mc_colors[bst_idx].channel[2] - val.channel[2];
1245 /*printf(" out(%3d, %3d, %3d) er(%4d, %4d, %4d)\n", quant->mc_colors[bst_idx].channel[0], quant->mc_colors[bst_idx].channel[1], quant->mc_colors[bst_idx].channel[2], perr.r, perr.g, perr.b);*/
1246 for (dx = 0; dx < mapw; ++dx) {
1247 for (dy = 0; dy < maph; ++dy) {
1248 err[x+dx+dy*errw].r += perr.r * map[dx+mapw*dy];
1249 err[x+dx+dy*errw].g += perr.g * map[dx+mapw*dy];
1250 err[x+dx+dy*errw].b += perr.b * map[dx+mapw*dy];
1251 }
1252 }
1253 *out++ = bst_idx;
1254 }
1255 /* shift up the error matrix */
1256 for (dy = 0; dy < maph-1; ++dy) {
1257 memcpy(err+dy*errw, err+(dy+1)*errw, sizeof(*err)*errw);
1258 }
1259 memset(err+(maph-1)*errw, 0, sizeof(*err)*errw);
1260 }
1261 CF_CLEANUP;
7fd765fe 1262 myfree(err);
02d1d628
AMH
1263}
1264/* Prescan finds the boxes in the image that have the highest number of colors
1265 and that result is used as the initial value for the vectores */
1266
1267
18accb2a 1268static void prescan(i_img **imgs,int count, int cnum, cvec *clr, i_sample_t *line) {
02d1d628 1269 int i,k,j,x,y;
18accb2a
TC
1270 i_sample_t *val;
1271 const int *chans;
02d1d628
AMH
1272
1273 pbox prebox[512];
1274 for(i=0;i<512;i++) {
1275 prebox[i].boxnum=i;
1276 prebox[i].pixcnt=0;
1277 prebox[i].cand=1;
1278 }
1279
1280 /* process each image */
1281 for (i = 0; i < count; ++i) {
1282 i_img *im = imgs[i];
18accb2a
TC
1283 chans = im->channels >= 3 ? NULL : gray_samples;
1284 for(y=0;y<im->ysize;y++) {
1285 i_gsamp(im, 0, im->xsize, y, line, chans, 3);
1286 val = line;
1287 for(x=0;x<im->xsize;x++) {
1288 prebox[pixbox_ch(val)].pixcnt++;
1289 }
02d1d628
AMH
1290 }
1291 }
1292
1293 for(i=0;i<512;i++) prebox[i].pdc=prebox[i].pixcnt;
1294 qsort(prebox,512,sizeof(pbox),(cmpfunc)pboxcmp);
1295
1296 for(i=0;i<cnum;i++) {
1297 /* printf("Color %d\n",i);
1298 for(k=0;k<10;k++) { printf("box=%03d %04d %d %04d \n",prebox[k].boxnum,prebox[k].pixcnt,prebox[k].cand,prebox[k].pdc); }
1299 printf("\n\n"); */
1300 reorder(prebox);
1301 }
1302
1303 /* for(k=0;k<cnum;k++) { printf("box=%03d %04d %d %04d \n",prebox[k].boxnum,prebox[k].pixcnt,prebox[k].cand,prebox[k].pdc); } */
1304
1305 k=0;
1306 j=1;
1307 i=0;
1308 while(i<cnum) {
1309 /* printf("prebox[%d].cand=%d\n",k,prebox[k].cand); */
36e67d0b 1310 if (clr[i].fixed) { i++; continue; } /* reserved go to next */
02d1d628
AMH
1311 if (j>=prebox[k].cand) { k++; j=1; } else {
1312 if (prebox[k].cand == 2) boxcenter(prebox[k].boxnum,&(clr[i]));
1313 else boxrand(prebox[k].boxnum,&(clr[i]));
1314 /* printf("(%d,%d) %d %d -> (%d,%d,%d)\n",k,j,prebox[k].boxnum,prebox[k].pixcnt,clr[i].r,clr[i].g,clr[i].b); */
1315 j++;
1316 i++;
1317 }
1318 }
1319}
1320
1321
1322static void reorder(pbox prescan[512]) {
1323 int nidx;
1324 pbox c;
1325
1326 nidx=0;
1327 c=prescan[0];
1328
1329 c.cand++;
1330 c.pdc=c.pixcnt/(c.cand*c.cand);
1331 /* c.pdc=c.pixcnt/c.cand; */
1332 while(c.pdc < prescan[nidx+1].pdc && nidx < 511) {
1333 prescan[nidx]=prescan[nidx+1];
1334 nidx++;
1335 }
1336 prescan[nidx]=c;
1337}
1338
1339static int
1340pboxcmp(const pbox *a,const pbox *b) {
1341 if (a->pixcnt > b->pixcnt) return -1;
1342 if (a->pixcnt < b->pixcnt) return 1;
1343 return 0;
1344}
1345
1346static void
1347boxcenter(int box,cvec *cv) {
1348 cv->r=15+((box&448)>>1);
1349 cv->g=15+((box&56)<<2);
1350 cv->b=15+((box&7)<<5);
1351}
1352
1353static void
1354bbox(int box,int *r0,int *r1,int *g0,int *g1,int *b0,int *b1) {
1355 *r0=(box&448)>>1;
1356 *r1=(*r0)|31;
1357 *g0=(box&56)<<2;
1358 *g1=(*g0)|31;
1359 *b0=(box&7)<<5;
1360 *b1=(*b0)|31;
1361}
1362
1363static void
1364boxrand(int box,cvec *cv) {
1365 cv->r=6+(rand()%25)+((box&448)>>1);
1366 cv->g=6+(rand()%25)+((box&56)<<2);
1367 cv->b=6+(rand()%25)+((box&7)<<5);
1368}
1369
1370static float
1371frandn(void) {
1372
1373 float u1,u2,w;
1374
1375 w=1;
1376
1377 while (w >= 1 || w == 0) {
1378 u1 = 2 * frand() - 1;
1379 u2 = 2 * frand() - 1;
1380 w = u1*u1 + u2*u2;
1381 }
1382
1383 w = sqrt((-2*log(w))/w);
1384 return u1*w;
1385}
1386
1387/* Create hash index */
1388static
1389void
1390cr_hashindex(cvec clr[256],int cnum,hashbox hb[512]) {
1391
1392 int bx,mind,cd,cumcnt,bst_idx,i;
1393/* printf("indexing... \n");*/
1394
1395 cumcnt=0;
1396 for(bx=0; bx<512; bx++) {
1397 mind=196608;
1398 for(i=0; i<cnum; i++) {
1399 cd = maxdist(bx,&clr[i]);
1400 if (cd < mind) { mind=cd; bst_idx=i; }
1401 }
1402
1403 hb[bx].cnt=0;
1404 for(i=0;i<cnum;i++) if (mindist(bx,&clr[i])<mind) hb[bx].vec[hb[bx].cnt++]=i;
1405 /*printf("box %d -> approx -> %d\n",bx,hb[bx].cnt); */
1406 /* statbox(bx,cnum,clr); */
1407 cumcnt+=hb[bx].cnt;
1408 }
1409
1410/* printf("Average search space: %d\n",cumcnt/512); */
1411}
1412
1413static int
1414maxdist(int boxnum,cvec *cv) {
1415 int r0,r1,g0,g1,b0,b1;
1416 int r,g,b,mr,mg,mb;
1417
1418 r=cv->r;
1419 g=cv->g;
1420 b=cv->b;
1421
1422 bbox(boxnum,&r0,&r1,&g0,&g1,&b0,&b1);
1423
b33c08f8
TC
1424 mr=i_max(abs(b-b0),abs(b-b1));
1425 mg=i_max(abs(g-g0),abs(g-g1));
1426 mb=i_max(abs(r-r0),abs(r-r1));
02d1d628
AMH
1427
1428 return PWR2(mr)+PWR2(mg)+PWR2(mb);
1429}
1430
1431static int
1432mindist(int boxnum,cvec *cv) {
1433 int r0,r1,g0,g1,b0,b1;
1434 int r,g,b,mr,mg,mb;
1435
1436 r=cv->r;
1437 g=cv->g;
1438 b=cv->b;
1439
1440 bbox(boxnum,&r0,&r1,&g0,&g1,&b0,&b1);
1441
1442 /* printf("box %d, (%d,%d,%d)-(%d,%d,%d) vec (%d,%d,%d) ",boxnum,r0,g0,b0,r1,g1,b1,r,g,b); */
1443
1444 if (r0<=r && r<=r1 && g0<=g && g<=g1 && b0<=b && b<=b1) return 0;
1445
b33c08f8
TC
1446 mr=i_min(abs(b-b0),abs(b-b1));
1447 mg=i_min(abs(g-g0),abs(g-g1));
1448 mb=i_min(abs(r-r0),abs(r-r1));
02d1d628
AMH
1449
1450 mr=PWR2(mr);
1451 mg=PWR2(mg);
1452 mb=PWR2(mb);
1453
1454 if (r0<=r && r<=r1 && g0<=g && g<=g1) return mb;
1455 if (r0<=r && r<=r1 && b0<=b && b<=b1) return mg;
1456 if (b0<=b && b<=b1 && g0<=g && g<=g1) return mr;
1457
1458 if (r0<=r && r<=r1) return mg+mb;
1459 if (g0<=g && g<=g1) return mr+mb;
1460 if (b0<=b && b<=b1) return mg+mr;
1461
1462 return mr+mg+mb;
1463}
1464
1465static void transparent_threshold(i_quantize *, i_palidx *, i_img *, i_palidx);
1466static void transparent_errdiff(i_quantize *, i_palidx *, i_img *, i_palidx);
1467static void transparent_ordered(i_quantize *, i_palidx *, i_img *, i_palidx);
1468
92bda632
TC
1469/*
1470=item i_quant_transparent(quant, data, img, trans_index)
1471
1472=category Image quantization
1473
1474Dither the alpha channel on I<img> into the palette indexes in
1475I<data>. Pixels to be transparent are replaced with I<trans_pixel>.
1476
1477The method used depends on the tr_* members of quant.
1478
1479=cut
1480*/
1481
1482void
1483i_quant_transparent(i_quantize *quant, i_palidx *data, i_img *img,
02d1d628
AMH
1484 i_palidx trans_index)
1485{
1486 switch (quant->transp) {
1487 case tr_none:
1488 break;
1489
1490 default:
1491 quant->tr_threshold = 128;
1492 /* fall through */
1493 case tr_threshold:
1494 transparent_threshold(quant, data, img, trans_index);
1495 break;
1496
1497 case tr_errdiff:
1498 transparent_errdiff(quant, data, img, trans_index);
1499 break;
1500
1501 case tr_ordered:
1502 transparent_ordered(quant, data, img, trans_index);
1503 break;
1504 }
1505}
1506
1507static void
1508transparent_threshold(i_quantize *quant, i_palidx *data, i_img *img,
1509 i_palidx trans_index)
1510{
1511 int x, y;
18accb2a
TC
1512 i_sample_t *line = mymalloc(img->xsize * sizeof(i_sample_t));
1513 int trans_chan = img->channels > 2 ? 3 : 1;
02d1d628
AMH
1514
1515 for (y = 0; y < img->ysize; ++y) {
18accb2a 1516 i_gsamp(img, 0, img->xsize, y, line, &trans_chan, 1);
02d1d628 1517 for (x = 0; x < img->xsize; ++x) {
18accb2a 1518 if (line[x] < quant->tr_threshold)
02d1d628
AMH
1519 data[y*img->xsize+x] = trans_index;
1520 }
1521 }
18accb2a 1522 myfree(line);
02d1d628
AMH
1523}
1524
1525static void
1526transparent_errdiff(i_quantize *quant, i_palidx *data, i_img *img,
1527 i_palidx trans_index)
1528{
1529 int *map;
1530 int index;
1531 int mapw, maph, mapo;
1532 int errw, *err, *errp;
1533 int difftotal, out, error;
1534 int x, y, dx, dy, i;
18accb2a
TC
1535 i_sample_t *line;
1536 int trans_chan = img->channels > 2 ? 3 : 1;
02d1d628
AMH
1537
1538 /* no custom map for transparency (yet) */
1539 index = quant->tr_errdiff & ed_mask;
1540 if (index >= ed_custom) index = ed_floyd;
1541 map = maps[index].map;
1542 mapw = maps[index].width;
1543 maph = maps[index].height;
1544 mapo = maps[index].orig;
1545
1546 errw = img->xsize+mapw-1;
1547 err = mymalloc(sizeof(*err) * maph * errw);
1548 errp = err+mapo;
1549 memset(err, 0, sizeof(*err) * maph * errw);
1550
18accb2a 1551 line = mymalloc(img->xsize * sizeof(i_sample_t));
02d1d628
AMH
1552 difftotal = 0;
1553 for (i = 0; i < maph * mapw; ++i)
1554 difftotal += map[i];
1555 for (y = 0; y < img->ysize; ++y) {
18accb2a 1556 i_gsamp(img, 0, img->xsize, y, line, &trans_chan, 1);
02d1d628 1557 for (x = 0; x < img->xsize; ++x) {
18accb2a
TC
1558 line[x] = g_sat(line[x]-errp[x]/difftotal);
1559 if (line[x] < 128) {
02d1d628
AMH
1560 out = 0;
1561 data[y*img->xsize+x] = trans_index;
1562 }
1563 else {
1564 out = 255;
1565 }
18accb2a 1566 error = out - line[x];
02d1d628
AMH
1567 for (dx = 0; dx < mapw; ++dx) {
1568 for (dy = 0; dy < maph; ++dy) {
1569 errp[x+dx-mapo+dy*errw] += error * map[dx+mapw*dy];
1570 }
1571 }
1572 }
1573 /* shift up the error matrix */
1574 for (dy = 0; dy < maph-1; ++dy)
1575 memcpy(err+dy*errw, err+(dy+1)*errw, sizeof(*err)*errw);
1576 memset(err+(maph-1)*errw, 0, sizeof(*err)*errw);
1577 }
18accb2a
TC
1578 myfree(err);
1579 myfree(line);
02d1d628
AMH
1580}
1581
1582/* builtin ordered dither maps */
b33c08f8
TC
1583static unsigned char
1584orddith_maps[][64] =
02d1d628
AMH
1585{
1586 { /* random
1587 this is purely random - it's pretty awful
1588 */
1589 48, 72, 196, 252, 180, 92, 108, 52,
1590 228, 176, 64, 8, 236, 40, 20, 164,
1591 120, 128, 84, 116, 24, 28, 172, 220,
1592 68, 0, 188, 124, 184, 224, 192, 104,
1593 132, 100, 240, 200, 152, 160, 244, 44,
1594 96, 204, 144, 16, 140, 56, 232, 216,
1595 208, 4, 76, 212, 136, 248, 80, 168,
1596 156, 88, 32, 112, 148, 12, 36, 60,
1597 },
1598 {
1599 /* dot8
1600 perl spot.perl '($x-3.5)*($x-3.5)+($y-3.5)*($y-3.5)'
1601 */
1602 240, 232, 200, 136, 140, 192, 228, 248,
1603 220, 148, 100, 76, 80, 104, 152, 212,
1604 180, 116, 56, 32, 36, 60, 120, 176,
1605 156, 64, 28, 0, 8, 44, 88, 160,
1606 128, 92, 24, 12, 4, 40, 68, 132,
1607 184, 96, 48, 20, 16, 52, 108, 188,
1608 216, 144, 112, 72, 84, 124, 164, 224,
1609 244, 236, 196, 168, 172, 204, 208, 252,
1610 },
1611 { /* dot4
1612 perl spot.perl \
1613 'min(dist(1.5, 1.5),dist(5.5,1.5),dist(1.5,5.5),dist(5.5,5.5))'
1614 */
1615 196, 72, 104, 220, 200, 80, 112, 224,
1616 76, 4, 24, 136, 84, 8, 32, 144,
1617 108, 28, 52, 168, 116, 36, 56, 176,
1618 216, 140, 172, 244, 228, 148, 180, 248,
1619 204, 92, 124, 236, 192, 68, 96, 208,
1620 88, 12, 44, 156, 64, 0, 16, 128,
1621 120, 40, 60, 188, 100, 20, 48, 160,
1622 232, 152, 184, 252, 212, 132, 164, 240,
1623 },
1624 { /* hline
1625 perl spot.perl '$y-3'
1626 */
1627 160, 164, 168, 172, 176, 180, 184, 188,
1628 128, 132, 136, 140, 144, 148, 152, 156,
1629 32, 36, 40, 44, 48, 52, 56, 60,
1630 0, 4, 8, 12, 16, 20, 24, 28,
1631 64, 68, 72, 76, 80, 84, 88, 92,
1632 96, 100, 104, 108, 112, 116, 120, 124,
1633 192, 196, 200, 204, 208, 212, 216, 220,
1634 224, 228, 232, 236, 240, 244, 248, 252,
1635 },
1636 { /* vline
1637 perl spot.perl '$x-3'
1638 */
1639 180, 100, 40, 12, 44, 104, 184, 232,
1640 204, 148, 60, 16, 64, 128, 208, 224,
1641 212, 144, 76, 8, 80, 132, 216, 244,
1642 160, 112, 68, 20, 84, 108, 172, 236,
1643 176, 96, 72, 28, 88, 152, 188, 228,
1644 200, 124, 92, 0, 32, 116, 164, 240,
1645 168, 120, 36, 24, 48, 136, 192, 248,
1646 196, 140, 52, 4, 56, 156, 220, 252,
1647 },
1648 { /* slashline
1649 perl spot.perl '$y+$x-7'
1650 */
1651 248, 232, 224, 192, 140, 92, 52, 28,
1652 240, 220, 196, 144, 108, 60, 12, 64,
1653 216, 180, 148, 116, 76, 20, 80, 128,
1654 204, 152, 104, 44, 16, 72, 100, 160,
1655 164, 96, 68, 24, 56, 112, 168, 176,
1656 124, 40, 8, 36, 88, 136, 184, 212,
1657 84, 4, 32, 120, 156, 188, 228, 236,
1658 0, 48, 132, 172, 200, 208, 244, 252,
1659 },
1660 { /* backline
1661 perl spot.perl '$y-$x'
1662 */
1663 0, 32, 116, 172, 184, 216, 236, 252,
1664 56, 8, 72, 132, 136, 200, 228, 240,
1665 100, 36, 12, 40, 92, 144, 204, 220,
1666 168, 120, 60, 16, 44, 96, 156, 176,
1667 180, 164, 112, 48, 28, 52, 128, 148,
1668 208, 192, 152, 88, 84, 20, 64, 104,
1669 232, 224, 196, 140, 108, 68, 24, 76,
1670 248, 244, 212, 188, 160, 124, 80, 4,
1671 },
11e7329d
TC
1672 {
1673 /* tiny
1674 good for display, bad for print
1675 hand generated
1676 */
1677 0, 128, 32, 192, 8, 136, 40, 200,
1678 224, 64, 160, 112, 232, 72, 168, 120,
1679 48, 144, 16, 208, 56, 152, 24, 216,
1680 176, 96, 240, 80, 184, 104, 248, 88,
1681 12, 140, 44, 204, 4, 132, 36, 196,
1682 236, 76, 172, 124, 228, 68, 164, 116,
1683 60, 156, 28, 220, 52, 148, 20, 212,
1684 188, 108, 252, 92, 180, 100, 244, 84,
1685 },
02d1d628
AMH
1686};
1687
1688static void
1689transparent_ordered(i_quantize *quant, i_palidx *data, i_img *img,
1690 i_palidx trans_index)
1691{
1692 unsigned char *spot;
1693 int x, y;
18accb2a
TC
1694 i_sample_t *line;
1695 int trans_chan = img->channels > 2 ? 3 : 1;
02d1d628
AMH
1696 if (quant->tr_orddith == od_custom)
1697 spot = quant->tr_custom;
1698 else
1699 spot = orddith_maps[quant->tr_orddith];
18accb2a
TC
1700
1701 line = mymalloc(img->xsize * sizeof(i_sample_t));
02d1d628 1702 for (y = 0; y < img->ysize; ++y) {
18accb2a 1703 i_gsamp(img, 0, img->xsize, y, line, &trans_chan, 1);
02d1d628 1704 for (x = 0; x < img->xsize; ++x) {
18accb2a 1705 if (line[x] < spot[(x&7)+(y&7)*8])
02d1d628
AMH
1706 data[x+y*img->xsize] = trans_index;
1707 }
1708 }
18accb2a 1709 myfree(line);
02d1d628 1710}
18accb2a 1711