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