avoid a possible sign-extension for offsets/sizes in SGI
[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     }
637     quant->mc_count = out;
638   }
639   else {
640     /* build the starting partition */
641     parts = i_mempool_alloc(&mp, sizeof(*parts) * quant->mc_size);
642     parts[0].start = 0;
643     parts[0].size = out;
644     parts[0].pixels = total_pixels;
645     calc_part(parts, colors);
646     color_count = 1;
647     
648     while (color_count < quant->mc_size) {
649       /* initialized to avoid compiler warnings */
650       int max_index = 0, max_ch = 0; /* index/channel with biggest spread */
651       int max_size;
652       medcut_partition *workpart;
653       i_img_dim cum_total;
654       i_img_dim half;
655       
656       /* find the partition with the most biggest span with more than 
657          one color */
658       max_size = -1;
659       for (i = 0; i < color_count; ++i) {
660         for (ch = 0; ch < chan_count; ++ch) {
661           if (parts[i].width[ch] > max_size 
662               && parts[i].size > 1) {
663             max_index = i;
664             max_ch = ch;
665             max_size = parts[i].width[ch];
666           }
667         }
668       }
669       
670       /* nothing else we can split */
671       if (max_size == -1)
672         break;
673       
674       workpart = parts+max_index;
675       /*printf("splitting partition %d (pixels %ld, start %d, size %d)\n", max_index, workpart->pixels, workpart->start, workpart->size);*/
676       qsort(colors + workpart->start, workpart->size, sizeof(*colors),
677             sorters[max_ch]);
678       
679       /* find the median or something like it we need to make sure both
680          sides of the split have at least one color in them, so we don't
681          test at the first or last entry */
682       i = workpart->start;
683       cum_total = colors[i].count;
684       ++i;
685       half = workpart->pixels / 2;
686       while (i < workpart->start + workpart->size - 1
687              && cum_total < half) {
688         cum_total += colors[i++].count;
689       }
690       /*printf("Split at %d to make %d (half %ld, cumtotal %ld)\n", i, color_count, half, cum_total);*/
691       
692       /* found the spot to split */
693       parts[color_count].start = i;
694       parts[color_count].size = workpart->start + workpart->size - i;
695       workpart->size = i - workpart->start;
696       parts[color_count].pixels = workpart->pixels - cum_total;
697       workpart->pixels = cum_total;
698       
699       /* recalculate the limits */
700       calc_part(workpart, colors);
701       calc_part(parts+color_count, colors);
702       ++color_count;
703     }
704     
705     /* fill in the color table - since we could still have partitions
706        that have more than one color, we need to average the colors */
707     for (part_num = 0; part_num < color_count; ++part_num) {
708       double sums[3];
709       medcut_partition *workpart;
710       
711       workpart = parts+part_num;
712       for (ch = 0; ch < 3; ++ch)
713         sums[ch] = 0;
714       
715       for (i = workpart->start; i < workpart->start + workpart->size; ++i) {
716         for (ch = 0; ch < 3; ++ch) {
717           sums[ch] += (int)(colors[i].rgb[ch]) * colors[i].count;
718         }
719       }
720       for (ch = 0; ch < 3; ++ch) {
721         quant->mc_colors[part_num].channel[ch] = sums[ch] / workpart->pixels;
722       }
723     }
724     quant->mc_count = color_count;
725   }
726   /*printf("out %d colors\n", quant->mc_count);*/
727   i_mempool_destroy(&mp);
728
729   mm_log((1, "makemap_mediancut() - %d colors\n", quant->mc_count));
730 }
731
732 static void
733 makemap_mono(i_quantize *quant) {
734   quant->mc_colors[0].rgba.r = 0;
735   quant->mc_colors[0].rgba.g = 0;
736   quant->mc_colors[0].rgba.b = 0;
737   quant->mc_colors[0].rgba.a = 255;
738   quant->mc_colors[1].rgba.r = 255;
739   quant->mc_colors[1].rgba.g = 255;
740   quant->mc_colors[1].rgba.b = 255;
741   quant->mc_colors[1].rgba.a = 255;
742   quant->mc_count = 2;
743 }
744
745 static void
746 makemap_gray(i_quantize *quant, int step) {
747   int gray = 0;
748   int i = 0;
749
750   while (gray < 256) {
751     setcol(quant->mc_colors+i, gray, gray, gray, 255);
752     ++i;
753     gray += step;
754   }
755   quant->mc_count = i;
756 }
757
758 static void
759 makemap_webmap(i_quantize *quant) {
760   int r, g, b;
761
762   int i = 0;
763   for (r = 0; r < 256; r+=0x33)
764     for (g = 0; g < 256; g+=0x33)
765       for (b = 0; b < 256; b += 0x33)
766         setcol(quant->mc_colors+i++, r, g, b, 255);
767   quant->mc_count = i;
768 }
769
770 static int 
771 in_palette(i_color *c, i_quantize *quant, int size) {
772   int i;
773
774   for (i = 0; i < size; ++i) {
775     if (c->channel[0] == quant->mc_colors[i].channel[0]
776         && c->channel[1] == quant->mc_colors[i].channel[1]
777         && c->channel[2] == quant->mc_colors[i].channel[2]) {
778       return i;
779     }
780   }
781
782   return -1;
783 }
784
785 /*
786 =item makemap_palette(quant, imgs, count)
787
788 Tests if all the given images are paletted and have a common palette,
789 if they do it builds that palette.
790
791 A possible improvement might be to eliminate unused colors in the
792 images palettes.
793
794 =cut
795 */
796
797 static int
798 makemap_palette(i_quantize *quant, i_img **imgs, int count) {
799   int size = quant->mc_count;
800   int i;
801   int imgn;
802   char used[256];
803   int col_count;
804
805   mm_log((1, "makemap_palette(quant %p { mc_count=%d, mc_colors=%p }, imgs %p, count %d)\n", 
806           quant, quant->mc_count, quant->mc_colors, imgs, count));
807   /* we try to build a common palette here, if we can manage that, then
808      that's the palette we use */
809   for (imgn = 0; imgn < count; ++imgn) {
810     int eliminate_unused;
811     if (imgs[imgn]->type != i_palette_type) {
812       mm_log((1, "makemap_palette() -> 0 (non-palette image)\n"));
813       return 0;
814     }
815
816     if (!i_tags_get_int(&imgs[imgn]->tags, "gif_eliminate_unused", 0, 
817                         &eliminate_unused)) {
818       eliminate_unused = 1;
819     }
820
821     if (eliminate_unused) {
822       i_palidx *line = mymalloc(sizeof(i_palidx) * imgs[imgn]->xsize);
823       i_img_dim x, y;
824       memset(used, 0, sizeof(used));
825
826       for (y = 0; y < imgs[imgn]->ysize; ++y) {
827         i_gpal(imgs[imgn], 0, imgs[imgn]->xsize, y, line);
828         for (x = 0; x < imgs[imgn]->xsize; ++x)
829           used[line[x]] = 1;
830       }
831
832       myfree(line);
833     }
834     else {
835       /* assume all are in use */
836       memset(used, 1, sizeof(used));
837     }
838
839     col_count = i_colorcount(imgs[imgn]);
840     for (i = 0; i < col_count; ++i) {
841       i_color c;
842       
843       i_getcolors(imgs[imgn], i, &c, 1);
844       if (used[i]) {
845         if (in_palette(&c, quant, size) < 0) {
846           if (size < quant->mc_size) {
847             quant->mc_colors[size++] = c;
848           }
849           else {
850             mm_log((1, "makemap_palette() -> 0 (too many colors)\n"));
851             return 0;
852           }
853         }
854       }
855     }
856   }
857
858   mm_log((1, "makemap_palette() -> 1 (%d total colors)\n", size));
859   quant->mc_count = size;
860
861   return 1;
862 }
863
864 #define pboxjump 32
865
866 /* Define one of the following 4 symbols to choose a colour search method
867    The idea is to try these out, including benchmarking, to see which
868    is fastest in a good spread of circumstances.
869    I'd expect IM_CFLINSEARCH to be fastest for very small palettes, and
870    IM_CFHASHBOX for large images with large palettes.
871
872    Some other possibilities include:
873     - search over entries sorted by luminance
874
875    Initially I was planning on testing using the macros and then
876    integrating the code directly into each function, but this means if
877    we find a bug at a late stage we will need to update N copies of
878    the same code.  Also, keeping the code in the macros means that the
879    code in the translation functions is much more to the point,
880    there's no distracting colour search code to remove attention from
881    what makes _this_ translation function different.  It may be
882    advisable to move the setup code into functions at some point, but
883    it should be possible to do this fairly transparently.
884
885    If IM_CF_COPTS is defined then CFLAGS must have an appropriate 
886    definition.
887
888    Each option needs to define 4 macros:
889     CF_VARS - variables to define in the function
890     CF_SETUP - code to setup for the colour search, eg. allocating and
891       initializing lookup tables
892     CF_FIND - code that looks for the color in val and puts the best 
893       matching index in bst_idx
894     CF_CLEANUP - code to clean up, eg. releasing memory
895 */
896 #ifndef IM_CF_COPTS
897 /*#define IM_CFLINSEARCH*/
898 #define IM_CFHASHBOX
899 /*#define IM_CFSORTCHAN*/
900 /*#define IM_CFRAND2DIST*/
901 #endif
902
903 /* return true if the color map contains only grays */
904 static int
905 is_gray_map(const i_quantize *quant) {
906   int i;
907
908   for (i = 0; i < quant->mc_count; ++i) {
909     if (quant->mc_colors[i].rgb.r != quant->mc_colors[i].rgb.g
910         || quant->mc_colors[i].rgb.r != quant->mc_colors[i].rgb.b) {
911       mm_log((1, "  not a gray map\n"));
912       return 0;
913     }
914   }
915
916   mm_log((1, "  is a gray map\n"));
917   return 1;
918 }
919
920 #ifdef IM_CFHASHBOX
921
922 /* The original version I wrote for this used the sort.
923    If this is defined then we use a sort to extract the indices for 
924    the hashbox */
925 #define HB_SORT
926
927 /* assume i is available */
928 #define CF_VARS hashbox *hb = mymalloc(sizeof(hashbox) * 512); \
929                int currhb;  \
930                long ld, cd
931
932 #ifdef HB_SORT
933
934 static long *gdists; /* qsort is annoying */
935 /* int might be smaller than long, so we need to do a real compare 
936    rather than a subtraction*/
937 static int distcomp(void const *a, void const *b) {
938   long ra = gdists[*(int const *)a];
939   long rb = gdists[*(int const *)b];
940   if (ra < rb)
941     return -1;
942   else if (ra > rb)
943     return 1;
944   else
945     return 0;
946 }
947
948 #endif
949
950 /* for each hashbox build a list of colours that are in the hb or is closer
951    than other colours
952    This is pretty involved.  The original gifquant generated the hashbox
953    as part of it's normal processing, but since the map generation is now 
954    separated from the translation we need to do this on the spot.
955    Any optimizations, even if they don't produce perfect results would be
956    welcome.
957  */
958 static void hbsetup(i_quantize *quant, hashbox *hb) {
959   long *dists, mind, maxd;
960   int cr, cb, cg, hbnum, i;
961   i_color cenc;
962 #ifdef HB_SORT
963   int *indices = mymalloc(quant->mc_count * sizeof(int)); 
964 #endif
965
966   dists = mymalloc(quant->mc_count * sizeof(long)); 
967   for (cr = 0; cr < 8; ++cr) { 
968     for (cg = 0; cg < 8; ++cg) { 
969       for (cb = 0; cb < 8; ++cb) { 
970         /* centre of the hashbox */ 
971         cenc.channel[0] = cr*pboxjump+pboxjump/2; 
972         cenc.channel[1] = cg*pboxjump+pboxjump/2; 
973         cenc.channel[2] = cb*pboxjump+pboxjump/2; 
974         hbnum = pixbox(&cenc); 
975         hb[hbnum].cnt = 0; 
976         /* order indices in the order of distance from the hashbox */ 
977         for (i = 0; i < quant->mc_count; ++i) { 
978 #ifdef HB_SORT
979           indices[i] = i; 
980 #endif
981           dists[i] = ceucl_d(&cenc, quant->mc_colors+i); 
982         } 
983 #ifdef HB_SORT
984         /* it should be possible to do this without a sort 
985            but so far I'm too lazy */
986         gdists = dists; 
987         qsort(indices, quant->mc_count, sizeof(int), distcomp); 
988         /* any colors that can match are within mind+diagonal size of 
989            a hashbox */ 
990         mind = dists[indices[0]]; 
991         i = 0; 
992         maxd = (sqrt(mind)+pboxjump)*(sqrt(mind)+pboxjump);
993         while (i < quant->mc_count && dists[indices[i]] < maxd) { 
994           hb[hbnum].vec[hb[hbnum].cnt++] = indices[i++]; 
995         } 
996 #else
997         /* work out the minimum */
998         mind = 256*256*3;
999         for (i = 0; i < quant->mc_count; ++i) {
1000           if (dists[i] < mind) mind = dists[i];
1001         }
1002         /* transfer any colours that might be closest to a colour in 
1003            this hashbox */
1004         maxd = (sqrt(mind)+pboxjump)*(sqrt(mind)+pboxjump);
1005         for (i = 0; i < quant->mc_count; ++i) {
1006           if (dists[i] < maxd)
1007             hb[hbnum].vec[hb[hbnum].cnt++] = i;
1008         }
1009 #endif
1010       } 
1011     } 
1012   }
1013 #ifdef HB_SORT
1014   myfree(indices); 
1015 #endif
1016   myfree(dists) ;
1017 }
1018 #define CF_SETUP hbsetup(quant, hb)
1019
1020 #define CF_FIND \
1021   currhb = pixbox(&val); \
1022   ld = 196608; \
1023   for (i = 0; i < hb[currhb].cnt; ++i) { \
1024     cd = ceucl_d(quant->mc_colors+hb[currhb].vec[i], &val); \
1025     if (cd < ld) { ld = cd; bst_idx = hb[currhb].vec[i]; } \
1026   }
1027
1028 #define CF_CLEANUP myfree(hb)
1029   
1030 #endif
1031
1032 #ifdef IM_CFLINSEARCH
1033 /* as simple as it gets */
1034 #define CF_VARS long ld, cd
1035 #define CF_SETUP /* none needed */
1036 #define CF_FIND \
1037    ld = 196608; \
1038    for (i = 0; i < quant->mc_count; ++i) { \
1039      cd = ceucl_d(quant->mc_colors+i, &val); \
1040      if (cd < ld) { ld = cd; bst_idx = i; } \
1041    }
1042 #define CF_CLEANUP
1043 #endif
1044
1045 #ifdef IM_CFSORTCHAN
1046 static int gsortchan;
1047 static i_quantize *gquant;
1048 static int chansort(void const *a, void const *b) {
1049   return gquant->mc_colors[*(int const *)a].channel[gsortchan] -
1050     gquant->mc_colors[*(int const *)b].channel[gsortchan];
1051 }
1052 #define CF_VARS int *indices, sortchan, diff; \
1053                 long ld, cd; \
1054                 int vindex[256] /* where to find value i of chan */
1055
1056 static void chansetup(i_img *img, i_quantize *quant, int *csortchan, 
1057                       int *vindex, int **cindices) {
1058   int *indices, sortchan, chan, i, chval;
1059   int chanmins[MAXCHANNELS], chanmaxs[MAXCHANNELS], maxrange;
1060
1061   /* find the channel with the maximum range */ 
1062   /* the maximum stddev would probably be better */
1063   for (chan = 0; chan < img->channels; ++chan) { 
1064     chanmins[chan] = 256; chanmaxs[chan] = 0; 
1065     for (i = 0; i < quant->mc_count; ++i) { 
1066       if (quant->mc_colors[i].channel[chan] < chanmins[chan]) 
1067         chanmins[chan] = quant->mc_colors[i].channel[chan]; 
1068       if (quant->mc_colors[i].channel[chan] > chanmaxs[chan]) 
1069         chanmaxs[chan] = quant->mc_colors[i].channel[chan]; 
1070     } 
1071   } 
1072   maxrange = -1; 
1073   for (chan = 0; chan < img->channels; ++chan) { 
1074     if (chanmaxs[chan]-chanmins[chan] > maxrange) { 
1075       maxrange = chanmaxs[chan]-chanmins[chan]; 
1076       sortchan = chan; 
1077     } 
1078   } 
1079   indices = mymalloc(quant->mc_count * sizeof(int)) ;
1080   for (i = 0; i < quant->mc_count; ++i) { 
1081     indices[i] = i; 
1082   } 
1083   gsortchan = sortchan; 
1084   gquant = quant; 
1085   qsort(indices, quant->mc_count, sizeof(int), chansort) ;
1086   /* now a lookup table to find entries faster */ 
1087   for (chval=0, i=0; i < quant->mc_count; ++i) { 
1088     while (chval < 256 && 
1089            chval < quant->mc_colors[indices[i]].channel[sortchan]) { 
1090       vindex[chval++] = i; 
1091     } 
1092   } 
1093   while (chval < 256) { 
1094     vindex[chval++] = quant->mc_count-1; 
1095   }
1096   *csortchan = sortchan;
1097   *cindices = indices;
1098 }
1099
1100 #define CF_SETUP \
1101   chansetup(img, quant, &sortchan, vindex, &indices)
1102
1103 int chanfind(i_color val, i_quantize *quant, int *indices, int *vindex, 
1104              int sortchan) {
1105   int i, bst_idx, diff, maxdiff;
1106   long ld, cd;
1107
1108   i = vindex[val.channel[sortchan]];
1109   bst_idx = indices[i];
1110   ld = 196608;
1111   diff = 0;
1112   maxdiff = quant->mc_count;
1113   while (diff < maxdiff) {
1114     if (i+diff < quant->mc_count) {
1115       cd = ceucl_d(&val, quant->mc_colors+indices[i+diff]); 
1116       if (cd < ld) {
1117         bst_idx = indices[i+diff];
1118         ld = cd;
1119         maxdiff = sqrt(ld);
1120       }
1121     }
1122     if (i-diff >= 0) {
1123       cd = ceucl_d(&val, quant->mc_colors+indices[i-diff]); 
1124       if (cd < ld) {
1125         bst_idx = indices[i-diff];
1126         ld = cd;
1127         maxdiff = sqrt(ld);
1128       }
1129     }
1130     ++diff;
1131   }
1132
1133   return bst_idx;
1134 }
1135
1136 #define CF_FIND \
1137   bst_idx = chanfind(val, quant, indices, vindex, sortchan)
1138   
1139
1140 #define CF_CLEANUP myfree(indices)
1141
1142 #endif
1143
1144 #ifdef IM_CFRAND2DIST
1145
1146 /* This is based on a method described by Addi in the #imager channel 
1147    on the 28/2/2001.  I was about 1am Sydney time at the time, so I 
1148    wasn't at my most cogent.  Well, that's my excuse :)
1149
1150 <TonyC> what I have at the moment is: hashboxes, with optimum hash box
1151 filling; simple linear search; and a lookup in the widest channel
1152 (currently the channel with the maximum range)
1153 <Addi> There is one more way that might be simple to implement.
1154 <Addi> You want to hear?
1155 <TonyC> what's that?
1156 <purl> somebody said that was not true
1157 <Addi> For each of the colors in the palette start by creating a
1158 sorted list of the form:
1159 <Addi> [distance, color]
1160 <Addi> Where they are sorted by distance.
1161 <TonyC> distance to where?
1162 <Addi> Where the elements in the lists are the distances and colors of
1163 the other colors in the palette
1164 <TonyC> ok
1165 <Addi> So if you are at color 0
1166 <Addi> ok - now to search for the closest color when you are creating
1167 the final image is done like this:
1168 <Addi> a) pick a random color from the palette
1169 <Addi> b) calculate the distance to it
1170 <Addi> c) only check the vectors that are within double the distance
1171 in the list of the color you picked from the palette.
1172 <Addi> Does that seem logical?
1173 <Addi> Lets imagine that we only have grayscale to make an example:
1174 <Addi> Our palette has 1 4 10 20 as colors.
1175 <Addi> And we want to quantize the color 11
1176 <Addi> lets say we picked 10 randomly
1177 <Addi> the double distance is 2
1178 <Addi> since abs(10-11)*2 is 2
1179 <Addi> And the list at vector 10 is this:
1180 <Addi> [0, 10], [6 4], [9, 1], [10, 20]
1181 <Addi> so we look at the first one (but not the second one since 6 is
1182 at a greater distance than 2.
1183 <Addi> Any of that make sense?
1184 <TonyC> yes, though are you suggesting another random jump to one of
1185 the colours with the possible choices? or an exhaustive search?
1186 <Addi> TonyC: It's possible to come up with a recursive/iterative 
1187 enhancement but this is the 'basic' version.
1188 <Addi> Which would do an iterative search.
1189 <Addi> You can come up with conditions where it pays to switch to a new one.
1190 <Addi> And the 'random' start can be switched over to a small tree.
1191 <Addi> So you would have a little index at the start.
1192 <Addi> to get you into the general direction
1193 <Addi> Perhaps just an 8 split.
1194 <Addi> that is - split each dimension in half.
1195 <TonyC> yep
1196 <TonyC> I get the idea
1197 <Addi> But this would seem to be a good approach in our case since we 
1198 usually have few codevectors.
1199 <Addi> So we only need 256*256 entries in a table.
1200 <Addi> We could even only index some of them that were deemed as good 
1201 candidates.
1202 <TonyC> I was considering adding paletted output support for PNG and 
1203 TIFF at some point, which support 16-bit palettes
1204 <Addi> ohh.
1205 <Addi> 'darn' ;)
1206
1207
1208 */
1209
1210
1211 typedef struct i_dists {
1212   int index;
1213   long dist;
1214 } i_dists;
1215
1216 #define CF_VARS \
1217     i_dists *dists;
1218
1219 static int dists_sort(void const *a, void const *b) {
1220   return ((i_dists *)a)->dist - ((i_dists *)b)->dist;
1221 }
1222
1223 static void rand2dist_setup(i_quantize *quant, i_dists **cdists) {
1224   i_dists *dists = 
1225     mymalloc(sizeof(i_dists)*quant->mc_count*quant->mc_count);
1226   int i, j;
1227   long cd;
1228   for (i = 0; i < quant->mc_count; ++i) {
1229     i_dists *ldists = dists + quant->mc_count * i;
1230     i_color val = quant->mc_colors[i];
1231     for (j = 0; j < quant->mc_count; ++j) {
1232       ldists[j].index = j;
1233       ldists[j].dist = ceucl_d(&val, quant->mc_colors+j);
1234     }
1235     qsort(ldists, quant->mc_count, sizeof(i_dists), dists_sort);
1236   }
1237   *cdists = dists;
1238 }
1239
1240 #define CF_SETUP \
1241                 bst_idx = rand() % quant->mc_count; \
1242                 rand2dist_setup(quant, &dists)
1243
1244 static int rand2dist_find(i_color val, i_quantize *quant, i_dists *dists, int index) {
1245   i_dists *cdists;
1246   long cd, ld;
1247   long maxld;
1248   int i;
1249   int bst_idx;
1250
1251   cdists = dists + index * quant->mc_count;
1252   ld = 3 * 256 * 256;
1253   maxld = 8 * ceucl_d(&val, quant->mc_colors+index);
1254   for (i = 0; i < quant->mc_count && cdists[i].dist <= maxld; ++i) {
1255     cd = ceucl_d(&val, quant->mc_colors+cdists[i].index);
1256     if (cd < ld) {
1257       bst_idx = cdists[i].index;
1258       ld = cd;
1259     }
1260   }
1261   return bst_idx;
1262 }
1263
1264 #define CF_FIND bst_idx = rand2dist_find(val, quant, dists, bst_idx)
1265
1266 #define CF_CLEANUP myfree(dists)
1267
1268
1269 #endif
1270
1271 static void translate_addi(i_quantize *quant, i_img *img, i_palidx *out) {
1272   i_img_dim x, y, k;
1273   int i, bst_idx = 0;
1274   i_color val;
1275   int pixdev = quant->perturb;
1276   CF_VARS;
1277
1278   CF_SETUP;
1279
1280   if (img->channels >= 3) {
1281     if (pixdev) {
1282       k=0;
1283       for(y=0;y<img->ysize;y++) for(x=0;x<img->xsize;x++) {
1284         i_gpix(img,x,y,&val);
1285         val.channel[0]=g_sat(val.channel[0]+(int)(pixdev*frandn()));
1286         val.channel[1]=g_sat(val.channel[1]+(int)(pixdev*frandn()));
1287         val.channel[2]=g_sat(val.channel[2]+(int)(pixdev*frandn()));
1288         CF_FIND;
1289         out[k++]=bst_idx;
1290       }
1291     } else {
1292       k=0;
1293       for(y=0;y<img->ysize;y++) for(x=0;x<img->xsize;x++) {
1294         i_gpix(img,x,y,&val);
1295         CF_FIND;
1296         out[k++]=bst_idx;
1297       }
1298     }
1299   }
1300   else {
1301     if (pixdev) {
1302       k=0;
1303       for(y=0;y<img->ysize;y++) for(x=0;x<img->xsize;x++) {
1304         i_gpix(img,x,y,&val);
1305         val.channel[1] = val.channel[2] =
1306           val.channel[0]=g_sat(val.channel[0]+(int)(pixdev*frandn()));
1307         CF_FIND;
1308         out[k++]=bst_idx;
1309       }
1310     } else {
1311       k=0;
1312       for(y=0;y<img->ysize;y++) for(x=0;x<img->xsize;x++) {
1313         i_gpix(img,x,y,&val);
1314         val.channel[1] = val.channel[2] = val.channel[0];
1315         CF_FIND;
1316         out[k++]=bst_idx;
1317       }
1318     }
1319   }
1320   CF_CLEANUP;
1321 }
1322
1323 static int floyd_map[] =
1324 {
1325   0, 0, 7,
1326   3, 5, 1
1327 };
1328
1329 static int jarvis_map[] =
1330 {
1331   0, 0, 0, 7, 5,
1332   3, 5, 7, 5, 3,
1333   1, 3, 5, 3, 1
1334 };
1335
1336 static int stucki_map[] =
1337 {
1338   0, 0, 0, 8, 4,
1339   2, 4, 8, 4, 2,
1340   1, 2, 4, 2, 1
1341 };
1342
1343 struct errdiff_map {
1344   int *map;
1345   int width, height, orig;
1346 };
1347
1348 static struct errdiff_map maps[] =
1349 {
1350   { floyd_map, 3, 2, 1 },
1351   { jarvis_map, 5, 3, 2 },
1352   { stucki_map, 5, 3, 2 },
1353 };
1354
1355 typedef struct errdiff_tag {
1356   int r, g, b;
1357 } errdiff_t;
1358
1359 /* perform an error diffusion dither */
1360 static int
1361 translate_errdiff(i_quantize *quant, i_img *img, i_palidx *out) {
1362   int *map;
1363   int mapw, maph, mapo;
1364   int i;
1365   errdiff_t *err;
1366   i_img_dim errw;
1367   int difftotal;
1368   i_img_dim x, y, dx, dy;
1369   int bst_idx = 0;
1370   int is_gray = is_gray_map(quant);
1371   CF_VARS;
1372
1373   if ((quant->errdiff & ed_mask) == ed_custom) {
1374     map = quant->ed_map;
1375     mapw = quant->ed_width;
1376     maph = quant->ed_height;
1377     mapo = quant->ed_orig;
1378   }
1379   else {
1380     int index = quant->errdiff & ed_mask;
1381     if (index >= ed_custom) index = ed_floyd;
1382     map = maps[index].map;
1383     mapw = maps[index].width;
1384     maph = maps[index].height;
1385     mapo = maps[index].orig;
1386   }
1387   
1388   difftotal = 0;
1389   for (i = 0; i < maph * mapw; ++i) {
1390     if (map[i] < 0) {
1391       i_push_errorf(0, "errdiff_map values must be non-negative, errdiff[%d] is negative", i);
1392       return 0;
1393     }
1394     difftotal += map[i];
1395   }
1396
1397   if (!difftotal) {
1398     i_push_error(0, "error diffusion map must contain some non-zero values");
1399     return 0;
1400   }
1401
1402   errw = img->xsize+mapw;
1403   err = mymalloc(sizeof(*err) * maph * errw);
1404   /*errp = err+mapo;*/
1405   memset(err, 0, sizeof(*err) * maph * errw);
1406   
1407   /*printf("map:\n");
1408  for (dy = 0; dy < maph; ++dy) {
1409    for (dx = 0; dx < mapw; ++dx) {
1410      printf("%2d", map[dx+dy*mapw]);
1411    }
1412    putchar('\n');
1413    }*/
1414
1415   CF_SETUP;
1416
1417   for (y = 0; y < img->ysize; ++y) {
1418     for (x = 0; x < img->xsize; ++x) {
1419       i_color val;
1420       errdiff_t perr;
1421       i_gpix(img, x, y, &val);
1422       if (img->channels < 3) {
1423         val.channel[1] = val.channel[2] = val.channel[0];
1424       }
1425       else if (is_gray) {
1426         int gray = 0.5 + color_to_grey(&val);
1427         val.channel[0] = val.channel[1] = val.channel[2] = gray;
1428       }
1429       perr = err[x+mapo];
1430       perr.r = perr.r < 0 ? -((-perr.r)/difftotal) : perr.r/difftotal;
1431       perr.g = perr.g < 0 ? -((-perr.g)/difftotal) : perr.g/difftotal;
1432       perr.b = perr.b < 0 ? -((-perr.b)/difftotal) : perr.b/difftotal;
1433       /*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);*/
1434       val.channel[0] = g_sat(val.channel[0]-perr.r);
1435       val.channel[1] = g_sat(val.channel[1]-perr.g);
1436       val.channel[2] = g_sat(val.channel[2]-perr.b);
1437       CF_FIND;
1438       /* save error */
1439       perr.r = quant->mc_colors[bst_idx].channel[0] - val.channel[0];
1440       perr.g = quant->mc_colors[bst_idx].channel[1] - val.channel[1];
1441       perr.b = quant->mc_colors[bst_idx].channel[2] - val.channel[2];
1442       /*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);*/
1443       for (dx = 0; dx < mapw; ++dx) {
1444         for (dy = 0; dy < maph; ++dy) {
1445           err[x+dx+dy*errw].r += perr.r * map[dx+mapw*dy];
1446           err[x+dx+dy*errw].g += perr.g * map[dx+mapw*dy];
1447           err[x+dx+dy*errw].b += perr.b * map[dx+mapw*dy];
1448         }
1449       }
1450       *out++ = bst_idx;
1451     }
1452     /* shift up the error matrix */
1453     for (dy = 0; dy < maph-1; ++dy) {
1454       memcpy(err+dy*errw, err+(dy+1)*errw, sizeof(*err)*errw);
1455     }
1456     memset(err+(maph-1)*errw, 0, sizeof(*err)*errw);
1457   }
1458   CF_CLEANUP;
1459   myfree(err);
1460
1461   return 1;
1462 }
1463 /* Prescan finds the boxes in the image that have the highest number of colors 
1464    and that result is used as the initial value for the vectores */
1465
1466
1467 static void prescan(i_img **imgs,int count, int cnum, cvec *clr, i_sample_t *line) {
1468   int i,k,j;
1469   i_img_dim x,y;
1470   i_sample_t *val;
1471   const int *chans;
1472
1473   pbox prebox[512];
1474   for(i=0;i<512;i++) {
1475     prebox[i].boxnum=i;
1476     prebox[i].pixcnt=0;
1477     prebox[i].cand=1;
1478   }
1479
1480   /* process each image */
1481   for (i = 0; i < count; ++i) {
1482     i_img *im = imgs[i];
1483     chans = im->channels >= 3 ? NULL : gray_samples;
1484     for(y=0;y<im->ysize;y++) {
1485       i_gsamp(im, 0, im->xsize, y, line, chans, 3);
1486       val = line;
1487       for(x=0;x<im->xsize;x++) {
1488         prebox[pixbox_ch(val)].pixcnt++;
1489       }
1490     }
1491   }
1492
1493   for(i=0;i<512;i++) prebox[i].pdc=prebox[i].pixcnt;
1494   qsort(prebox,512,sizeof(pbox),(cmpfunc)pboxcmp);
1495
1496   for(i=0;i<cnum;i++) {
1497     /*      printf("Color %d\n",i); 
1498             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); } 
1499             printf("\n\n"); */
1500     reorder(prebox);
1501   }
1502   
1503   /*    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); } */
1504   
1505   k=0;
1506   j=1;
1507   i=0;
1508   while(i<cnum) {
1509     /*    printf("prebox[%d].cand=%d\n",k,prebox[k].cand); */
1510     if (clr[i].fixed) { i++; continue; } /* reserved go to next */
1511     if (j>=prebox[k].cand) { k++; j=1; } else {
1512       if (prebox[k].cand == 2) boxcenter(prebox[k].boxnum,&(clr[i]));
1513       else boxrand(prebox[k].boxnum,&(clr[i]));
1514       /*      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); */
1515       j++;
1516       i++;
1517     }
1518   }
1519 }
1520   
1521
1522 static void reorder(pbox prescan[512]) {
1523   int nidx;
1524   pbox c;
1525
1526   nidx=0;
1527   c=prescan[0];
1528   
1529   c.cand++;
1530   c.pdc=c.pixcnt/(c.cand*c.cand); 
1531   /*  c.pdc=c.pixcnt/c.cand; */
1532   while(nidx < 511 && c.pdc < prescan[nidx+1].pdc) {
1533     prescan[nidx]=prescan[nidx+1];
1534     nidx++;
1535   }
1536   prescan[nidx]=c;
1537 }
1538
1539 static int
1540 pboxcmp(const pbox *a,const pbox *b) {
1541   if (a->pixcnt > b->pixcnt) return -1;
1542   if (a->pixcnt < b->pixcnt) return 1;
1543   return 0;
1544 }
1545
1546 static void
1547 boxcenter(int box,cvec *cv) {
1548   cv->r=15+((box&448)>>1);
1549   cv->g=15+((box&56)<<2);
1550   cv->b=15+((box&7)<<5);
1551 }
1552
1553 static void
1554 bbox(int box,int *r0,int *r1,int *g0,int *g1,int *b0,int *b1) {
1555   *r0=(box&448)>>1;
1556   *r1=(*r0)|31;
1557   *g0=(box&56)<<2;
1558   *g1=(*g0)|31;
1559   *b0=(box&7)<<5;
1560   *b1=(*b0)|31;
1561 }
1562
1563 static void
1564 boxrand(int box,cvec *cv) {
1565   cv->r=6+(rand()%25)+((box&448)>>1);
1566   cv->g=6+(rand()%25)+((box&56)<<2);
1567   cv->b=6+(rand()%25)+((box&7)<<5);
1568 }
1569
1570 static float
1571 frandn(void) {
1572
1573   float u1,u2,w;
1574   
1575   w=1;
1576   
1577   while (w >= 1 || w == 0) {
1578     u1 = 2 * frand() - 1;
1579     u2 = 2 * frand() - 1;
1580     w = u1*u1 + u2*u2;
1581   }
1582   
1583   w = sqrt((-2*log(w))/w);
1584   return u1*w;
1585 }
1586
1587 /* Create hash index */
1588 static
1589 void
1590 cr_hashindex(cvec clr[256],int cnum,hashbox hb[512]) {
1591   
1592   int bx,mind,cd,cumcnt,i;
1593 /*  printf("indexing... \n");*/
1594   
1595   cumcnt=0;
1596   for(bx=0; bx<512; bx++) {
1597     mind=196608;
1598     for(i=0; i<cnum; i++) { 
1599       cd = maxdist(bx,&clr[i]);
1600       if (cd < mind) { mind=cd; } 
1601     }
1602     
1603     hb[bx].cnt=0;
1604     for(i=0;i<cnum;i++) if (mindist(bx,&clr[i])<mind) hb[bx].vec[hb[bx].cnt++]=i;
1605     /*printf("box %d -> approx -> %d\n",bx,hb[bx].cnt); */
1606     /*  statbox(bx,cnum,clr); */
1607     cumcnt+=hb[bx].cnt;
1608   }
1609   
1610 /*  printf("Average search space: %d\n",cumcnt/512); */
1611 }
1612
1613 static int
1614 maxdist(int boxnum,cvec *cv) {
1615   int r0,r1,g0,g1,b0,b1;
1616   int r,g,b,mr,mg,mb;
1617
1618   r=cv->r;
1619   g=cv->g;
1620   b=cv->b;
1621   
1622   bbox(boxnum,&r0,&r1,&g0,&g1,&b0,&b1);
1623
1624   mr=i_max(abs(b-b0),abs(b-b1));
1625   mg=i_max(abs(g-g0),abs(g-g1));
1626   mb=i_max(abs(r-r0),abs(r-r1));
1627   
1628   return PWR2(mr)+PWR2(mg)+PWR2(mb);
1629 }
1630
1631 static int
1632 mindist(int boxnum,cvec *cv) {
1633   int r0,r1,g0,g1,b0,b1;
1634   int r,g,b,mr,mg,mb;
1635
1636   r=cv->r;
1637   g=cv->g;
1638   b=cv->b;
1639   
1640   bbox(boxnum,&r0,&r1,&g0,&g1,&b0,&b1);
1641
1642   /*  printf("box %d, (%d,%d,%d)-(%d,%d,%d) vec (%d,%d,%d) ",boxnum,r0,g0,b0,r1,g1,b1,r,g,b); */
1643
1644   if (r0<=r && r<=r1 && g0<=g && g<=g1 && b0<=b && b<=b1) return 0;
1645
1646   mr=i_min(abs(b-b0),abs(b-b1));
1647   mg=i_min(abs(g-g0),abs(g-g1));
1648   mb=i_min(abs(r-r0),abs(r-r1));
1649   
1650   mr=PWR2(mr);
1651   mg=PWR2(mg);
1652   mb=PWR2(mb);
1653
1654   if (r0<=r && r<=r1 && g0<=g && g<=g1) return mb;
1655   if (r0<=r && r<=r1 && b0<=b && b<=b1) return mg;
1656   if (b0<=b && b<=b1 && g0<=g && g<=g1) return mr;
1657
1658   if (r0<=r && r<=r1) return mg+mb;
1659   if (g0<=g && g<=g1) return mr+mb;
1660   if (b0<=b && b<=b1) return mg+mr;
1661
1662   return mr+mg+mb;
1663 }
1664
1665 static void transparent_threshold(i_quantize *, i_palidx *, i_img *, i_palidx);
1666 static void transparent_errdiff(i_quantize *, i_palidx *, i_img *, i_palidx);
1667 static void transparent_ordered(i_quantize *, i_palidx *, i_img *, i_palidx);
1668
1669 /*
1670 =item i_quant_transparent(C<quant>, C<data>, C<img>, C<trans_index>)
1671
1672 =category Image quantization
1673
1674 Dither the alpha channel on C<img> into the palette indexes in
1675 C<data>.  Pixels to be transparent are replaced with C<trans_pixel>.
1676
1677 The method used depends on the tr_* members of C<quant>.
1678
1679 =cut
1680 */
1681
1682 void 
1683 i_quant_transparent(i_quantize *quant, i_palidx *data, i_img *img,
1684                        i_palidx trans_index)
1685 {
1686   switch (quant->transp) {
1687   case tr_none:
1688     break;
1689     
1690   default:
1691     quant->tr_threshold = 128;
1692     /* fall through */
1693   case tr_threshold:
1694     transparent_threshold(quant, data, img, trans_index);
1695     break;
1696     
1697   case tr_errdiff:
1698     transparent_errdiff(quant, data, img, trans_index);
1699     break;
1700
1701   case tr_ordered:
1702     transparent_ordered(quant, data, img, trans_index);
1703     break;
1704   }
1705 }
1706
1707 static void
1708 transparent_threshold(i_quantize *quant, i_palidx *data, i_img *img,
1709                       i_palidx trans_index)
1710 {
1711   i_img_dim x, y;
1712   i_sample_t *line = mymalloc(img->xsize * sizeof(i_sample_t));
1713   int trans_chan = img->channels > 2 ? 3 : 1;
1714   
1715   for (y = 0; y < img->ysize; ++y) {
1716     i_gsamp(img, 0, img->xsize, y, line, &trans_chan, 1);
1717     for (x = 0; x < img->xsize; ++x) {
1718       if (line[x] < quant->tr_threshold)
1719         data[y*img->xsize+x] = trans_index;
1720     }
1721   }
1722   myfree(line);
1723 }
1724
1725 static void
1726 transparent_errdiff(i_quantize *quant, i_palidx *data, i_img *img,
1727                     i_palidx trans_index)
1728 {
1729   int *map;
1730   int index;
1731   int mapw, maph, mapo;
1732   int errw, *err, *errp;
1733   int difftotal, out, error;
1734   i_img_dim x, y, dx, dy;
1735   int i;
1736   i_sample_t *line;
1737   int trans_chan = img->channels > 2 ? 3 : 1;
1738
1739   /* no custom map for transparency (yet) */
1740   index = quant->tr_errdiff & ed_mask;
1741   if (index >= ed_custom) index = ed_floyd;
1742   map = maps[index].map;
1743   mapw = maps[index].width;
1744   maph = maps[index].height;
1745   mapo = maps[index].orig;
1746
1747   errw = img->xsize+mapw-1;
1748   err = mymalloc(sizeof(*err) * maph * errw);
1749   errp = err+mapo;
1750   memset(err, 0, sizeof(*err) * maph * errw);
1751
1752   line = mymalloc(img->xsize * sizeof(i_sample_t));
1753   difftotal = 0;
1754   for (i = 0; i < maph * mapw; ++i)
1755     difftotal += map[i];
1756   for (y = 0; y < img->ysize; ++y) {
1757     i_gsamp(img, 0, img->xsize, y, line, &trans_chan, 1);
1758     for (x = 0; x < img->xsize; ++x) {
1759       line[x] = g_sat(line[x]-errp[x]/difftotal);
1760       if (line[x] < 128) {
1761         out = 0;
1762         data[y*img->xsize+x] = trans_index;
1763       }
1764       else {
1765         out = 255;
1766       }
1767       error = out - line[x];
1768       for (dx = 0; dx < mapw; ++dx) {
1769         for (dy = 0; dy < maph; ++dy) {
1770           errp[x+dx-mapo+dy*errw] += error * map[dx+mapw*dy];
1771         }
1772       }
1773     }
1774     /* shift up the error matrix */
1775     for (dy = 0; dy < maph-1; ++dy)
1776       memcpy(err+dy*errw, err+(dy+1)*errw, sizeof(*err)*errw);
1777     memset(err+(maph-1)*errw, 0, sizeof(*err)*errw);
1778   }
1779   myfree(err);
1780   myfree(line);
1781 }
1782
1783 /* builtin ordered dither maps */
1784 static unsigned char 
1785 orddith_maps[][64] =
1786 {
1787   { /* random 
1788        this is purely random - it's pretty awful
1789      */
1790      48,  72, 196, 252, 180,  92, 108,  52,
1791     228, 176,  64,   8, 236,  40,  20, 164,
1792     120, 128,  84, 116,  24,  28, 172, 220,
1793      68,   0, 188, 124, 184, 224, 192, 104,
1794     132, 100, 240, 200, 152, 160, 244,  44,
1795      96, 204, 144,  16, 140,  56, 232, 216,
1796     208,   4,  76, 212, 136, 248,  80, 168,
1797     156,  88,  32, 112, 148,  12,  36,  60,
1798   },
1799   {
1800     /* dot8
1801        perl spot.perl '($x-3.5)*($x-3.5)+($y-3.5)*($y-3.5)'
1802      */
1803     240, 232, 200, 136, 140, 192, 228, 248,
1804     220, 148, 100,  76,  80, 104, 152, 212,
1805     180, 116,  56,  32,  36,  60, 120, 176,
1806     156,  64,  28,   0,   8,  44,  88, 160,
1807     128,  92,  24,  12,   4,  40,  68, 132,
1808     184,  96,  48,  20,  16,  52, 108, 188,
1809     216, 144, 112,  72,  84, 124, 164, 224,
1810     244, 236, 196, 168, 172, 204, 208, 252,
1811   },
1812   { /* dot4
1813        perl spot.perl \
1814        'min(dist(1.5, 1.5),dist(5.5,1.5),dist(1.5,5.5),dist(5.5,5.5))'  
1815     */
1816     196,  72, 104, 220, 200,  80, 112, 224,
1817      76,   4,  24, 136,  84,   8,  32, 144,
1818     108,  28,  52, 168, 116,  36,  56, 176,
1819     216, 140, 172, 244, 228, 148, 180, 248,
1820     204,  92, 124, 236, 192,  68,  96, 208,
1821      88,  12,  44, 156,  64,   0,  16, 128,
1822     120,  40,  60, 188, 100,  20,  48, 160,
1823     232, 152, 184, 252, 212, 132, 164, 240,
1824   },
1825   { /* hline 
1826        perl spot.perl '$y-3'
1827      */
1828     160, 164, 168, 172, 176, 180, 184, 188,
1829     128, 132, 136, 140, 144, 148, 152, 156,
1830      32,  36,  40,  44,  48,  52,  56,  60,
1831       0,   4,   8,  12,  16,  20,  24,  28,
1832      64,  68,  72,  76,  80,  84,  88,  92,
1833      96, 100, 104, 108, 112, 116, 120, 124,
1834     192, 196, 200, 204, 208, 212, 216, 220,
1835     224, 228, 232, 236, 240, 244, 248, 252,
1836   },
1837   { /* vline 
1838        perl spot.perl '$x-3'
1839      */
1840     180, 100,  40,  12,  44, 104, 184, 232,
1841     204, 148,  60,  16,  64, 128, 208, 224,
1842     212, 144,  76,   8,  80, 132, 216, 244,
1843     160, 112,  68,  20,  84, 108, 172, 236,
1844     176,  96,  72,  28,  88, 152, 188, 228,
1845     200, 124,  92,   0,  32, 116, 164, 240,
1846     168, 120,  36,  24,  48, 136, 192, 248,
1847     196, 140,  52,   4,  56, 156, 220, 252,
1848   },
1849   { /* slashline 
1850        perl spot.perl '$y+$x-7'  
1851     */
1852     248, 232, 224, 192, 140,  92,  52,  28,
1853     240, 220, 196, 144, 108,  60,  12,  64,
1854     216, 180, 148, 116,  76,  20,  80, 128,
1855     204, 152, 104,  44,  16,  72, 100, 160,
1856     164,  96,  68,  24,  56, 112, 168, 176,
1857     124,  40,   8,  36,  88, 136, 184, 212,
1858      84,   4,  32, 120, 156, 188, 228, 236,
1859       0,  48, 132, 172, 200, 208, 244, 252,
1860   },
1861   { /* backline 
1862        perl spot.perl '$y-$x'
1863      */
1864       0,  32, 116, 172, 184, 216, 236, 252,
1865      56,   8,  72, 132, 136, 200, 228, 240,
1866     100,  36,  12,  40,  92, 144, 204, 220,
1867     168, 120,  60,  16,  44,  96, 156, 176,
1868     180, 164, 112,  48,  28,  52, 128, 148,
1869     208, 192, 152,  88,  84,  20,  64, 104,
1870     232, 224, 196, 140, 108,  68,  24,  76,
1871     248, 244, 212, 188, 160, 124,  80,   4,
1872   },
1873   {
1874     /* tiny
1875        good for display, bad for print
1876        hand generated
1877     */
1878       0, 128,  32, 192,   8, 136,  40, 200,
1879     224,  64, 160, 112, 232,  72, 168, 120,
1880      48, 144,  16, 208,  56, 152,  24, 216,
1881     176,  96, 240,  80, 184, 104, 248,  88,
1882      12, 140,  44, 204,   4, 132,  36, 196,
1883     236,  76, 172, 124, 228,  68, 164, 116,
1884      60, 156,  28, 220,  52, 148,  20, 212,
1885     188, 108, 252,  92, 180, 100, 244,  84,
1886   },
1887 };
1888
1889 static void
1890 transparent_ordered(i_quantize *quant, i_palidx *data, i_img *img,
1891                     i_palidx trans_index)
1892 {
1893   unsigned char *spot;
1894   i_img_dim x, y;
1895   i_sample_t *line;
1896   int trans_chan = img->channels > 2 ? 3 : 1;
1897   if (quant->tr_orddith == od_custom)
1898     spot = quant->tr_custom;
1899   else
1900     spot = orddith_maps[quant->tr_orddith];
1901
1902   line = mymalloc(img->xsize * sizeof(i_sample_t));
1903   for (y = 0; y < img->ysize; ++y) {
1904     i_gsamp(img, 0, img->xsize, y, line, &trans_chan, 1);
1905     for (x = 0; x < img->xsize; ++x) {
1906       if (line[x] < spot[(x&7)+(y&7)*8])
1907         data[x+y*img->xsize] = trans_index;
1908     }
1909   }
1910   myfree(line);
1911 }
1912