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