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