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