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