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