Leolo's guassian2 patch
[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       goto fail;
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     goto fail;
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  fail:
1464   CF_CLEANUP;
1465
1466   return 0;
1467 }
1468 /* Prescan finds the boxes in the image that have the highest number of colors 
1469    and that result is used as the initial value for the vectores */
1470
1471
1472 static void prescan(i_img **imgs,int count, int cnum, cvec *clr, i_sample_t *line) {
1473   int i,k,j;
1474   i_img_dim x,y;
1475   i_sample_t *val;
1476   const int *chans;
1477
1478   pbox prebox[512];
1479   for(i=0;i<512;i++) {
1480     prebox[i].boxnum=i;
1481     prebox[i].pixcnt=0;
1482     prebox[i].cand=1;
1483   }
1484
1485   /* process each image */
1486   for (i = 0; i < count; ++i) {
1487     i_img *im = imgs[i];
1488     chans = im->channels >= 3 ? NULL : gray_samples;
1489     for(y=0;y<im->ysize;y++) {
1490       i_gsamp(im, 0, im->xsize, y, line, chans, 3);
1491       val = line;
1492       for(x=0;x<im->xsize;x++) {
1493         prebox[pixbox_ch(val)].pixcnt++;
1494       }
1495     }
1496   }
1497
1498   for(i=0;i<512;i++) prebox[i].pdc=prebox[i].pixcnt;
1499   qsort(prebox,512,sizeof(pbox),(cmpfunc)pboxcmp);
1500
1501   for(i=0;i<cnum;i++) {
1502     /*      printf("Color %d\n",i); 
1503             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); } 
1504             printf("\n\n"); */
1505     reorder(prebox);
1506   }
1507   
1508   /*    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); } */
1509   
1510   k=0;
1511   j=1;
1512   i=0;
1513   while(i<cnum) {
1514     /*    printf("prebox[%d].cand=%d\n",k,prebox[k].cand); */
1515     if (clr[i].fixed) { i++; continue; } /* reserved go to next */
1516     if (j>=prebox[k].cand) { k++; j=1; } else {
1517       if (prebox[k].cand == 2) boxcenter(prebox[k].boxnum,&(clr[i]));
1518       else boxrand(prebox[k].boxnum,&(clr[i]));
1519       /*      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); */
1520       j++;
1521       i++;
1522     }
1523   }
1524 }
1525   
1526
1527 static void reorder(pbox prescan[512]) {
1528   int nidx;
1529   pbox c;
1530
1531   nidx=0;
1532   c=prescan[0];
1533   
1534   c.cand++;
1535   c.pdc=c.pixcnt/(c.cand*c.cand); 
1536   /*  c.pdc=c.pixcnt/c.cand; */
1537   while(nidx < 511 && c.pdc < prescan[nidx+1].pdc) {
1538     prescan[nidx]=prescan[nidx+1];
1539     nidx++;
1540   }
1541   prescan[nidx]=c;
1542 }
1543
1544 static int
1545 pboxcmp(const pbox *a,const pbox *b) {
1546   if (a->pixcnt > b->pixcnt) return -1;
1547   if (a->pixcnt < b->pixcnt) return 1;
1548   return 0;
1549 }
1550
1551 static void
1552 boxcenter(int box,cvec *cv) {
1553   cv->r=15+((box&448)>>1);
1554   cv->g=15+((box&56)<<2);
1555   cv->b=15+((box&7)<<5);
1556 }
1557
1558 static void
1559 bbox(int box,int *r0,int *r1,int *g0,int *g1,int *b0,int *b1) {
1560   *r0=(box&448)>>1;
1561   *r1=(*r0)|31;
1562   *g0=(box&56)<<2;
1563   *g1=(*g0)|31;
1564   *b0=(box&7)<<5;
1565   *b1=(*b0)|31;
1566 }
1567
1568 static void
1569 boxrand(int box,cvec *cv) {
1570   cv->r=6+(rand()%25)+((box&448)>>1);
1571   cv->g=6+(rand()%25)+((box&56)<<2);
1572   cv->b=6+(rand()%25)+((box&7)<<5);
1573 }
1574
1575 static float
1576 frandn(void) {
1577
1578   float u1,u2,w;
1579   
1580   w=1;
1581   
1582   while (w >= 1 || w == 0) {
1583     u1 = 2 * frand() - 1;
1584     u2 = 2 * frand() - 1;
1585     w = u1*u1 + u2*u2;
1586   }
1587   
1588   w = sqrt((-2*log(w))/w);
1589   return u1*w;
1590 }
1591
1592 /* Create hash index */
1593 static
1594 void
1595 cr_hashindex(cvec clr[256],int cnum,hashbox hb[512]) {
1596   
1597   int bx,mind,cd,cumcnt,i;
1598 /*  printf("indexing... \n");*/
1599   
1600   cumcnt=0;
1601   for(bx=0; bx<512; bx++) {
1602     mind=196608;
1603     for(i=0; i<cnum; i++) { 
1604       cd = maxdist(bx,&clr[i]);
1605       if (cd < mind) { mind=cd; } 
1606     }
1607     
1608     hb[bx].cnt=0;
1609     for(i=0;i<cnum;i++) if (mindist(bx,&clr[i])<mind) hb[bx].vec[hb[bx].cnt++]=i;
1610     /*printf("box %d -> approx -> %d\n",bx,hb[bx].cnt); */
1611     /*  statbox(bx,cnum,clr); */
1612     cumcnt+=hb[bx].cnt;
1613   }
1614   
1615 /*  printf("Average search space: %d\n",cumcnt/512); */
1616 }
1617
1618 static int
1619 maxdist(int boxnum,cvec *cv) {
1620   int r0,r1,g0,g1,b0,b1;
1621   int r,g,b,mr,mg,mb;
1622
1623   r=cv->r;
1624   g=cv->g;
1625   b=cv->b;
1626   
1627   bbox(boxnum,&r0,&r1,&g0,&g1,&b0,&b1);
1628
1629   mr=i_max(abs(b-b0),abs(b-b1));
1630   mg=i_max(abs(g-g0),abs(g-g1));
1631   mb=i_max(abs(r-r0),abs(r-r1));
1632   
1633   return PWR2(mr)+PWR2(mg)+PWR2(mb);
1634 }
1635
1636 static int
1637 mindist(int boxnum,cvec *cv) {
1638   int r0,r1,g0,g1,b0,b1;
1639   int r,g,b,mr,mg,mb;
1640
1641   r=cv->r;
1642   g=cv->g;
1643   b=cv->b;
1644   
1645   bbox(boxnum,&r0,&r1,&g0,&g1,&b0,&b1);
1646
1647   /*  printf("box %d, (%d,%d,%d)-(%d,%d,%d) vec (%d,%d,%d) ",boxnum,r0,g0,b0,r1,g1,b1,r,g,b); */
1648
1649   if (r0<=r && r<=r1 && g0<=g && g<=g1 && b0<=b && b<=b1) return 0;
1650
1651   mr=i_min(abs(b-b0),abs(b-b1));
1652   mg=i_min(abs(g-g0),abs(g-g1));
1653   mb=i_min(abs(r-r0),abs(r-r1));
1654   
1655   mr=PWR2(mr);
1656   mg=PWR2(mg);
1657   mb=PWR2(mb);
1658
1659   if (r0<=r && r<=r1 && g0<=g && g<=g1) return mb;
1660   if (r0<=r && r<=r1 && b0<=b && b<=b1) return mg;
1661   if (b0<=b && b<=b1 && g0<=g && g<=g1) return mr;
1662
1663   if (r0<=r && r<=r1) return mg+mb;
1664   if (g0<=g && g<=g1) return mr+mb;
1665   if (b0<=b && b<=b1) return mg+mr;
1666
1667   return mr+mg+mb;
1668 }
1669
1670 static void transparent_threshold(i_quantize *, i_palidx *, i_img *, i_palidx);
1671 static void transparent_errdiff(i_quantize *, i_palidx *, i_img *, i_palidx);
1672 static void transparent_ordered(i_quantize *, i_palidx *, i_img *, i_palidx);
1673
1674 /*
1675 =item i_quant_transparent(C<quant>, C<data>, C<img>, C<trans_index>)
1676
1677 =category Image quantization
1678
1679 Dither the alpha channel on C<img> into the palette indexes in
1680 C<data>.  Pixels to be transparent are replaced with C<trans_pixel>.
1681
1682 The method used depends on the tr_* members of C<quant>.
1683
1684 =cut
1685 */
1686
1687 void 
1688 i_quant_transparent(i_quantize *quant, i_palidx *data, i_img *img,
1689                        i_palidx trans_index)
1690 {
1691   switch (quant->transp) {
1692   case tr_none:
1693     break;
1694     
1695   default:
1696     quant->tr_threshold = 128;
1697     /* fall through */
1698   case tr_threshold:
1699     transparent_threshold(quant, data, img, trans_index);
1700     break;
1701     
1702   case tr_errdiff:
1703     transparent_errdiff(quant, data, img, trans_index);
1704     break;
1705
1706   case tr_ordered:
1707     transparent_ordered(quant, data, img, trans_index);
1708     break;
1709   }
1710 }
1711
1712 static void
1713 transparent_threshold(i_quantize *quant, i_palidx *data, i_img *img,
1714                       i_palidx trans_index)
1715 {
1716   i_img_dim x, y;
1717   i_sample_t *line = mymalloc(img->xsize * sizeof(i_sample_t));
1718   int trans_chan = img->channels > 2 ? 3 : 1;
1719   
1720   for (y = 0; y < img->ysize; ++y) {
1721     i_gsamp(img, 0, img->xsize, y, line, &trans_chan, 1);
1722     for (x = 0; x < img->xsize; ++x) {
1723       if (line[x] < quant->tr_threshold)
1724         data[y*img->xsize+x] = trans_index;
1725     }
1726   }
1727   myfree(line);
1728 }
1729
1730 static void
1731 transparent_errdiff(i_quantize *quant, i_palidx *data, i_img *img,
1732                     i_palidx trans_index)
1733 {
1734   int *map;
1735   int index;
1736   int mapw, maph, mapo;
1737   int errw, *err, *errp;
1738   int difftotal, out, error;
1739   i_img_dim x, y, dx, dy;
1740   int i;
1741   i_sample_t *line;
1742   int trans_chan = img->channels > 2 ? 3 : 1;
1743
1744   /* no custom map for transparency (yet) */
1745   index = quant->tr_errdiff & ed_mask;
1746   if (index >= ed_custom) index = ed_floyd;
1747   map = maps[index].map;
1748   mapw = maps[index].width;
1749   maph = maps[index].height;
1750   mapo = maps[index].orig;
1751
1752   errw = img->xsize+mapw-1;
1753   err = mymalloc(sizeof(*err) * maph * errw);
1754   errp = err+mapo;
1755   memset(err, 0, sizeof(*err) * maph * errw);
1756
1757   line = mymalloc(img->xsize * sizeof(i_sample_t));
1758   difftotal = 0;
1759   for (i = 0; i < maph * mapw; ++i)
1760     difftotal += map[i];
1761   for (y = 0; y < img->ysize; ++y) {
1762     i_gsamp(img, 0, img->xsize, y, line, &trans_chan, 1);
1763     for (x = 0; x < img->xsize; ++x) {
1764       line[x] = g_sat(line[x]-errp[x]/difftotal);
1765       if (line[x] < 128) {
1766         out = 0;
1767         data[y*img->xsize+x] = trans_index;
1768       }
1769       else {
1770         out = 255;
1771       }
1772       error = out - line[x];
1773       for (dx = 0; dx < mapw; ++dx) {
1774         for (dy = 0; dy < maph; ++dy) {
1775           errp[x+dx-mapo+dy*errw] += error * map[dx+mapw*dy];
1776         }
1777       }
1778     }
1779     /* shift up the error matrix */
1780     for (dy = 0; dy < maph-1; ++dy)
1781       memcpy(err+dy*errw, err+(dy+1)*errw, sizeof(*err)*errw);
1782     memset(err+(maph-1)*errw, 0, sizeof(*err)*errw);
1783   }
1784   myfree(err);
1785   myfree(line);
1786 }
1787
1788 /* builtin ordered dither maps */
1789 static unsigned char 
1790 orddith_maps[][64] =
1791 {
1792   { /* random 
1793        this is purely random - it's pretty awful
1794      */
1795      48,  72, 196, 252, 180,  92, 108,  52,
1796     228, 176,  64,   8, 236,  40,  20, 164,
1797     120, 128,  84, 116,  24,  28, 172, 220,
1798      68,   0, 188, 124, 184, 224, 192, 104,
1799     132, 100, 240, 200, 152, 160, 244,  44,
1800      96, 204, 144,  16, 140,  56, 232, 216,
1801     208,   4,  76, 212, 136, 248,  80, 168,
1802     156,  88,  32, 112, 148,  12,  36,  60,
1803   },
1804   {
1805     /* dot8
1806        perl spot.perl '($x-3.5)*($x-3.5)+($y-3.5)*($y-3.5)'
1807      */
1808     240, 232, 200, 136, 140, 192, 228, 248,
1809     220, 148, 100,  76,  80, 104, 152, 212,
1810     180, 116,  56,  32,  36,  60, 120, 176,
1811     156,  64,  28,   0,   8,  44,  88, 160,
1812     128,  92,  24,  12,   4,  40,  68, 132,
1813     184,  96,  48,  20,  16,  52, 108, 188,
1814     216, 144, 112,  72,  84, 124, 164, 224,
1815     244, 236, 196, 168, 172, 204, 208, 252,
1816   },
1817   { /* dot4
1818        perl spot.perl \
1819        'min(dist(1.5, 1.5),dist(5.5,1.5),dist(1.5,5.5),dist(5.5,5.5))'  
1820     */
1821     196,  72, 104, 220, 200,  80, 112, 224,
1822      76,   4,  24, 136,  84,   8,  32, 144,
1823     108,  28,  52, 168, 116,  36,  56, 176,
1824     216, 140, 172, 244, 228, 148, 180, 248,
1825     204,  92, 124, 236, 192,  68,  96, 208,
1826      88,  12,  44, 156,  64,   0,  16, 128,
1827     120,  40,  60, 188, 100,  20,  48, 160,
1828     232, 152, 184, 252, 212, 132, 164, 240,
1829   },
1830   { /* hline 
1831        perl spot.perl '$y-3'
1832      */
1833     160, 164, 168, 172, 176, 180, 184, 188,
1834     128, 132, 136, 140, 144, 148, 152, 156,
1835      32,  36,  40,  44,  48,  52,  56,  60,
1836       0,   4,   8,  12,  16,  20,  24,  28,
1837      64,  68,  72,  76,  80,  84,  88,  92,
1838      96, 100, 104, 108, 112, 116, 120, 124,
1839     192, 196, 200, 204, 208, 212, 216, 220,
1840     224, 228, 232, 236, 240, 244, 248, 252,
1841   },
1842   { /* vline 
1843        perl spot.perl '$x-3'
1844      */
1845     180, 100,  40,  12,  44, 104, 184, 232,
1846     204, 148,  60,  16,  64, 128, 208, 224,
1847     212, 144,  76,   8,  80, 132, 216, 244,
1848     160, 112,  68,  20,  84, 108, 172, 236,
1849     176,  96,  72,  28,  88, 152, 188, 228,
1850     200, 124,  92,   0,  32, 116, 164, 240,
1851     168, 120,  36,  24,  48, 136, 192, 248,
1852     196, 140,  52,   4,  56, 156, 220, 252,
1853   },
1854   { /* slashline 
1855        perl spot.perl '$y+$x-7'  
1856     */
1857     248, 232, 224, 192, 140,  92,  52,  28,
1858     240, 220, 196, 144, 108,  60,  12,  64,
1859     216, 180, 148, 116,  76,  20,  80, 128,
1860     204, 152, 104,  44,  16,  72, 100, 160,
1861     164,  96,  68,  24,  56, 112, 168, 176,
1862     124,  40,   8,  36,  88, 136, 184, 212,
1863      84,   4,  32, 120, 156, 188, 228, 236,
1864       0,  48, 132, 172, 200, 208, 244, 252,
1865   },
1866   { /* backline 
1867        perl spot.perl '$y-$x'
1868      */
1869       0,  32, 116, 172, 184, 216, 236, 252,
1870      56,   8,  72, 132, 136, 200, 228, 240,
1871     100,  36,  12,  40,  92, 144, 204, 220,
1872     168, 120,  60,  16,  44,  96, 156, 176,
1873     180, 164, 112,  48,  28,  52, 128, 148,
1874     208, 192, 152,  88,  84,  20,  64, 104,
1875     232, 224, 196, 140, 108,  68,  24,  76,
1876     248, 244, 212, 188, 160, 124,  80,   4,
1877   },
1878   {
1879     /* tiny
1880        good for display, bad for print
1881        hand generated
1882     */
1883       0, 128,  32, 192,   8, 136,  40, 200,
1884     224,  64, 160, 112, 232,  72, 168, 120,
1885      48, 144,  16, 208,  56, 152,  24, 216,
1886     176,  96, 240,  80, 184, 104, 248,  88,
1887      12, 140,  44, 204,   4, 132,  36, 196,
1888     236,  76, 172, 124, 228,  68, 164, 116,
1889      60, 156,  28, 220,  52, 148,  20, 212,
1890     188, 108, 252,  92, 180, 100, 244,  84,
1891   },
1892 };
1893
1894 static void
1895 transparent_ordered(i_quantize *quant, i_palidx *data, i_img *img,
1896                     i_palidx trans_index)
1897 {
1898   unsigned char *spot;
1899   i_img_dim x, y;
1900   i_sample_t *line;
1901   int trans_chan = img->channels > 2 ? 3 : 1;
1902   if (quant->tr_orddith == od_custom)
1903     spot = quant->tr_custom;
1904   else
1905     spot = orddith_maps[quant->tr_orddith];
1906
1907   line = mymalloc(img->xsize * sizeof(i_sample_t));
1908   for (y = 0; y < img->ysize; ++y) {
1909     i_gsamp(img, 0, img->xsize, y, line, &trans_chan, 1);
1910     for (x = 0; x < img->xsize; ++x) {
1911       if (line[x] < spot[(x&7)+(y&7)*8])
1912         data[x+y*img->xsize] = trans_index;
1913     }
1914   }
1915   myfree(line);
1916 }
1917