quant.c

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