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