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