src/image/convert.cxx

   1
   2 #define JPEG_IMAGES
   3 #include "picinfo.h"
   4 #include "accel.h"
   5
   6
   7 static int im_check_for_small_amount_of_colors(PICINFO& pic, int max_colors,
   8                                                byte* pic8_buf)
   9 {
  10   unsigned long colors[256],color;
  11   int          i, num_color, j;
  12   byte         *p, *pix;
  13
  14   if (max_colors>256)
  15     max_colors = 256;
  16
  17   num_color = 0;
  18
  19   for (i=pic.w*pic.h,p=pic.pic; i>0; i--)
  20   {
  21     //construct RGB 24 bit color
  22     color  = (((u_long) *p++) << 16);
  23     color += (((u_long) *p++) << 8);
  24     color +=  *p++;
  25
  26     for(j=0;j<num_color;j++)
  27       if(color==colors[j])
  28         break;
  29
  30     if(j<num_color) //We found color -> goto next pixel
  31       continue;
  32
  33     if (num_color >= max_colors)
  34       return 0;
  35
  36     colors[num_color++]=color;
  37   }
  38
  39
  40   //Now convert pixmap with our pallete
  41   for (i=pic.w*pic.h,p=pic.pic,pix=pic8_buf; i>0; i--,pix++)
  42   {
  43     color  = (((u_long) *p++) << 16);
  44     color += (((u_long) *p++) << 8);
  45     color +=  *p++;
  46
  47     for(j=0;j<num_color;j++)
  48       if(colors[j]==color)
  49       {
  50         *pix=j;
  51         break;
  52       }
  53   }
  54
  55   //And set colormap
  56   for (i=0; i<num_color; i++) {
  57     pic.pal[i*3]   =  colors[i]>>16;
  58     pic.pal[i*3+1] = (colors[i]>>8) & 0xff;
  59     pic.pal[i*3+2] =  colors[i]     & 0xff;
  60   }
  61
  62   pic.pic=pic8_buf;
  63
  64   return 1;
  65 }
  66
  67 static int slow_quant(byte* pic24, int w, int h, byte* pic8,
  68                       byte* rm,byte *gm,byte *bm, int max_colors);
  69
  70 static int im_do_convertion_to_num_colors(PICINFO& pic, int max_colors,
  71                                                byte* pic8_buf)
  72 {
  73   byte  *redmap,*greenmap,*bluemap;
  74   int   i;
  75
  76   redmap= new byte[256];
  77   greenmap= new byte[256];
  78   bluemap= new byte[256];
  79
  80   slow_quant(pic.pic,pic.w,pic.h,pic8_buf,
  81              redmap, greenmap, bluemap, max_colors);
  82
  83   for(i=0;i<max_colors;i++)
  84   {
  85     pic.pal[i*3]  =redmap[i];
  86     pic.pal[i*3+1]=greenmap[i];
  87     pic.pal[i*3+2]=bluemap[i];
  88   }
  89
  90   delete redmap;
  91   delete greenmap;
  92   delete bluemap;
  93
  94   return 1;
  95 }
  96
  97
  98 int im_convert_true_to_pseudo(PICINFO& pic,int max_colors)
  99 {
 100   int      ret;
 101   byte     *pic8_buf,*pic24_buf;
 102   pic8_buf= new byte[pic.w*pic.h];
 103   pic24_buf=pic.pic;
 104
 105   if(im_check_for_small_amount_of_colors(pic,max_colors,pic8_buf))
 106   {
 107     delete pic24_buf;
 108     return 1;
 109   }
 110   ret=im_do_convertion_to_num_colors(pic,max_colors,pic8_buf);
 111   if(ret)
 112     delete pic24_buf;
 113   else
 114     delete pic8_buf;
 115   return ret;
 116 }
 117
 118
 119
 120 /////////////////////////////////////////////////////////////////////////
 121 /////////////////////////////////////////////////////////////////////////
 122
 123 /***************************************************************/
 124 /* The following is based on jquant2.c from version 5          */
 125 /* of the IJG JPEG software, which is                          */
 126 /* Copyright (C) 1991-1994, Thomas G. Lane.                    */
 127 /***************************************************************/
 128
 129 #define INT32         int
 130 #define MAXNUMCOLORS  256        /* maximum size of colormap */
 131
 132 #define C0_SCALE 2                /* scale R distances by this much */
 133 #define C1_SCALE 3                /* scale G distances by this much */
 134 #define C2_SCALE 1                /* and B by this much */
 135
 136 #define HIST_C0_BITS  5                /* bits of precision in R histogram */
 137 #define HIST_C1_BITS  6                /* bits of precision in G histogram */
 138 #define HIST_C2_BITS  5                /* bits of precision in B histogram */
 139
 140 /* Number of elements along histogram axes. */
 141 #define HIST_C0_ELEMS  (1<<HIST_C0_BITS)
 142 #define HIST_C1_ELEMS  (1<<HIST_C1_BITS)
 143 #define HIST_C2_ELEMS  (1<<HIST_C2_BITS)
 144
 145 /* These are the amounts to shift an input value to get a histogram index. */
 146 #define C0_SHIFT  (8-HIST_C0_BITS)
 147 #define C1_SHIFT  (8-HIST_C1_BITS)
 148 #define C2_SHIFT  (8-HIST_C2_BITS)
 149
 150
 151 typedef unsigned char JSAMPLE;
 152 typedef JSAMPLE * JSAMPROW;
 153
 154 typedef u_short histcell;        /* histogram cell; prefer an unsigned type */
 155
 156 typedef histcell * histptr;        /* for pointers to histogram cells */
 157
 158 typedef histcell hist1d[HIST_C2_ELEMS]; /* typedefs for the histogram array */
 159 typedef hist1d hist2d[HIST_C1_ELEMS];
 160 typedef hist2d hist3d[HIST_C0_ELEMS];
 161
 162 typedef short FSERROR;                /* 16 bits should be enough */
 163 typedef int LOCFSERROR;                /* use 'int' for calculation temps */
 164
 165 typedef FSERROR *FSERRPTR;        /* pointer to error array */
 166
 167 typedef struct
 168 {
 169   /* The bounds of the box (inclusive); expressed as histogram indexes */
 170   int c0min, c0max;
 171   int c1min, c1max;
 172   int c2min, c2max;
 173   /* The volume (actually 2-norm) of the box */
 174   INT32 volume;
 175   /* The number of nonzero histogram cells within this box */
 176   long colorcount;
 177 } box;
 178 typedef box * boxptr;
 179
 180 /* Local state for the IJG quantizer */
 181
 182 static hist2d * sl_histogram;        /* pointer to the 3D histogram array */
 183 static FSERRPTR sl_fserrors;        /* accumulated-errors array */
 184 static int * sl_error_limiter;        /* table for clamping the applied error */
 185 static int sl_on_odd_row;        /* flag to remember which row we are on */
 186 static JSAMPROW sl_colormap[3];        /* selected colormap */
 187 static int sl_num_colors;        /* number of selected colors */
 188
 189
 190 static void   slow_fill_histogram (byte*, int);
 191 static boxptr find_biggest_color_pop (boxptr, int);
 192 static boxptr find_biggest_volume (boxptr, int);
 193 static void   update_box (boxptr);
 194 static int    median_cut (boxptr, int, int);
 195 static void   compute_color (boxptr, int);
 196 static void   slow_select_colors (int);
 197 static int    find_nearby_colors (int, int, int, JSAMPLE []);
 198 static void   find_best_colors (int,int,int,int, JSAMPLE [], JSAMPLE []);
 199 static void   fill_inverse_cmap (int, int, int);
 200 static void   slow_map_pixels (byte*, int, int, byte*);
 201 static void   init_error_limit(void);
 202
 203
 204 /* Master control for slow quantizer. */
 205 static int slow_quant(byte* pic24, int w, int h, byte* pic8,
 206                       byte* rm,byte *gm,byte *bm, int max_colors)
 207 {
 208   size_t fs_arraysize = (w + 2) * (3 * sizeof(FSERROR));
 209
 210   /* Allocate all the temporary storage needed */
 211   if (sl_error_limiter == NULL)
 212     init_error_limit();
 213   sl_histogram = (hist2d *) malloc(sizeof(hist3d));
 214   sl_fserrors = (FSERRPTR) malloc(fs_arraysize);
 215
 216   if (! sl_error_limiter || ! sl_histogram || ! sl_fserrors) {
 217     /* we never free sl_error_limiter once acquired */
 218     if (sl_histogram) free(sl_histogram);
 219     if (sl_fserrors) free(sl_fserrors);
 220     fprintf(stderr,"Slow_quant() - failed to allocate workspace\n");
 221     return 1;
 222   }
 223
 224   sl_colormap[0] = (JSAMPROW) rm;
 225   sl_colormap[1] = (JSAMPROW) gm;
 226   sl_colormap[2] = (JSAMPROW) bm;
 227
 228   /* Compute the color histogram */
 229   slow_fill_histogram(pic24, w*h);
 230
 231   /* Select the colormap */
 232   slow_select_colors(max_colors);
 233
 234   /* Zero the histogram: now to be used as inverse color map */
 235   im_setzero((char *) sl_histogram, sizeof(hist3d));
 236
 237   /* Initialize the propagated errors to zero. */
 238   im_setzero((char *) sl_fserrors, fs_arraysize);
 239   sl_on_odd_row = FALSE;
 240
 241   /* Map the image. */
 242   slow_map_pixels(pic24, w, h, pic8);
 243
 244   /* Release working memory. */
 245   /* we never free sl_error_limiter once acquired */
 246   free(sl_histogram);
 247   free(sl_fserrors);
 248
 249   return 0;
 250 }
 251
 252
 253 static void slow_fill_histogram (byte* pic24,int numpixels)
 254 {
 255   histptr histp;
 256   hist2d * histogram = sl_histogram;
 257
 258   im_setzero((char *) histogram, sizeof(hist3d));
 259
 260   while (numpixels-- > 0)
 261   {
 262     /* get pixel value and index into the histogram */
 263     histp = & histogram[pic24[0] >> C0_SHIFT]
 264                        [pic24[1] >> C1_SHIFT]
 265                        [pic24[2] >> C2_SHIFT];
 266     /* increment, check for overflow and undo increment if so. */
 267     if (++(*histp) <= 0)
 268       (*histp)--;
 269     pic24 += 3;
 270   }
 271 }
 272
 273
 274 static boxptr find_biggest_color_pop (boxptr boxlist, int numboxes)
 275 {
 276   register boxptr boxp;
 277   register int i;
 278   register long maxc = 0;
 279   boxptr which = NULL;
 280
 281   for (i = 0, boxp = boxlist; i < numboxes; i++, boxp++) {
 282     if (boxp->colorcount > maxc && boxp->volume > 0) {
 283       which = boxp;
 284       maxc = boxp->colorcount;
 285     }
 286   }
 287   return which;
 288 }
 289
 290
 291 static boxptr find_biggest_volume (boxptr boxlist, int numboxes)
 292 {
 293   register boxptr boxp;
 294   register int i;
 295   register INT32 maxv = 0;
 296   boxptr which = NULL;
 297
 298   for (i = 0, boxp = boxlist; i < numboxes; i++, boxp++) {
 299     if (boxp->volume > maxv) {
 300       which = boxp;
 301       maxv = boxp->volume;
 302     }
 303   }
 304   return which;
 305 }
 306
 307
 308 static void update_box (boxptr boxp)
 309 {
 310   hist2d * histogram = sl_histogram;
 311   histptr histp;
 312   int c0,c1,c2;
 313   int c0min,c0max,c1min,c1max,c2min,c2max;
 314   INT32 dist0,dist1,dist2;
 315   long ccount;
 316
 317   c0min = boxp->c0min;  c0max = boxp->c0max;
 318   c1min = boxp->c1min;  c1max = boxp->c1max;
 319   c2min = boxp->c2min;  c2max = boxp->c2max;
 320
 321   if (c0max > c0min)
 322     for (c0 = c0min; c0 <= c0max; c0++)
 323       for (c1 = c1min; c1 <= c1max; c1++) {
 324         histp = & histogram[c0][c1][c2min];
 325         for (c2 = c2min; c2 <= c2max; c2++)
 326           if (*histp++ != 0) {
 327             boxp->c0min = c0min = c0;
 328             goto have_c0min;
 329           }
 330       }
 331  have_c0min:
 332   if (c0max > c0min)
 333     for (c0 = c0max; c0 >= c0min; c0--)
 334       for (c1 = c1min; c1 <= c1max; c1++) {
 335         histp = & histogram[c0][c1][c2min];
 336         for (c2 = c2min; c2 <= c2max; c2++)
 337           if (*histp++ != 0) {
 338             boxp->c0max = c0max = c0;
 339             goto have_c0max;
 340           }
 341       }
 342  have_c0max:
 343   if (c1max > c1min)
 344     for (c1 = c1min; c1 <= c1max; c1++)
 345       for (c0 = c0min; c0 <= c0max; c0++) {
 346         histp = & histogram[c0][c1][c2min];
 347         for (c2 = c2min; c2 <= c2max; c2++)
 348           if (*histp++ != 0) {
 349             boxp->c1min = c1min = c1;
 350             goto have_c1min;
 351           }
 352       }
 353  have_c1min:
 354   if (c1max > c1min)
 355     for (c1 = c1max; c1 >= c1min; c1--)
 356       for (c0 = c0min; c0 <= c0max; c0++) {
 357         histp = & histogram[c0][c1][c2min];
 358         for (c2 = c2min; c2 <= c2max; c2++)
 359           if (*histp++ != 0) {
 360             boxp->c1max = c1max = c1;
 361             goto have_c1max;
 362           }
 363       }
 364  have_c1max:
 365   if (c2max > c2min)
 366     for (c2 = c2min; c2 <= c2max; c2++)
 367       for (c0 = c0min; c0 <= c0max; c0++) {
 368         histp = & histogram[c0][c1min][c2];
 369         for (c1 = c1min; c1 <= c1max; c1++, histp += HIST_C2_ELEMS)
 370           if (*histp != 0) {
 371             boxp->c2min = c2min = c2;
 372             goto have_c2min;
 373           }
 374       }
 375  have_c2min:
 376   if (c2max > c2min)
 377     for (c2 = c2max; c2 >= c2min; c2--)
 378       for (c0 = c0min; c0 <= c0max; c0++) {
 379         histp = & histogram[c0][c1min][c2];
 380         for (c1 = c1min; c1 <= c1max; c1++, histp += HIST_C2_ELEMS)
 381           if (*histp != 0) {
 382             boxp->c2max = c2max = c2;
 383             goto have_c2max;
 384           }
 385       }
 386  have_c2max:
 387
 388   dist0 = ((c0max - c0min) << C0_SHIFT) * C0_SCALE;
 389   dist1 = ((c1max - c1min) << C1_SHIFT) * C1_SCALE;
 390   dist2 = ((c2max - c2min) << C2_SHIFT) * C2_SCALE;
 391   boxp->volume = dist0*dist0 + dist1*dist1 + dist2*dist2;
 392
 393   ccount = 0;
 394   for (c0 = c0min; c0 <= c0max; c0++)
 395     for (c1 = c1min; c1 <= c1max; c1++) {
 396       histp = & histogram[c0][c1][c2min];
 397       for (c2 = c2min; c2 <= c2max; c2++, histp++)
 398         if (*histp != 0) {
 399           ccount++;
 400         }
 401     }
 402   boxp->colorcount = ccount;
 403 }
 404
 405
 406 static int median_cut (boxptr boxlist, int numboxes, int desired_colors)
 407 {
 408   int n,lb;
 409   int c0,c1,c2,cmax;
 410   register boxptr b1,b2;
 411
 412   while (numboxes < desired_colors) {
 413     /* Select box to split.
 414      * Current algorithm: by population for first half, then by volume.
 415      */
 416     if (numboxes*2 <= desired_colors) {
 417       b1 = find_biggest_color_pop(boxlist, numboxes);
 418     } else {
 419       b1 = find_biggest_volume(boxlist, numboxes);
 420     }
 421     if (b1 == NULL)                /* no splittable boxes left! */
 422       break;
 423     b2 = &boxlist[numboxes];        /* where new box will go */
 424     /* Copy the color bounds to the new box. */
 425     b2->c0max = b1->c0max; b2->c1max = b1->c1max; b2->c2max = b1->c2max;
 426     b2->c0min = b1->c0min; b2->c1min = b1->c1min; b2->c2min = b1->c2min;
 427     /* Choose which axis to split the box on.
 428      */
 429     c0 = ((b1->c0max - b1->c0min) << C0_SHIFT) * C0_SCALE;
 430     c1 = ((b1->c1max - b1->c1min) << C1_SHIFT) * C1_SCALE;
 431     c2 = ((b1->c2max - b1->c2min) << C2_SHIFT) * C2_SCALE;
 432     cmax = c1; n = 1;
 433     if (c0 > cmax) { cmax = c0; n = 0; }
 434     if (c2 > cmax) { n = 2; }
 435     switch (n) {
 436     case 0:
 437       lb = (b1->c0max + b1->c0min) / 2;
 438       b1->c0max = lb;
 439       b2->c0min = lb+1;
 440       break;
 441     case 1:
 442       lb = (b1->c1max + b1->c1min) / 2;
 443       b1->c1max = lb;
 444       b2->c1min = lb+1;
 445       break;
 446     case 2:
 447       lb = (b1->c2max + b1->c2min) / 2;
 448       b1->c2max = lb;
 449       b2->c2min = lb+1;
 450       break;
 451     }
 452     /* Update stats for boxes */
 453     update_box(b1);
 454     update_box(b2);
 455     numboxes++;
 456   }
 457   return numboxes;
 458 }
 459
 460
 461 static void compute_color (boxptr boxp,int icolor)
 462 {
 463   /* Current algorithm: mean weighted by pixels (not colors) */
 464   /* Note it is important to get the rounding correct! */
 465   hist2d * histogram = sl_histogram;
 466   histptr histp;
 467   int c0,c1,c2;
 468   int c0min,c0max,c1min,c1max,c2min,c2max;
 469   long count;
 470   long total = 0;
 471   long c0total = 0;
 472   long c1total = 0;
 473   long c2total = 0;
 474
 475   c0min = boxp->c0min;  c0max = boxp->c0max;
 476   c1min = boxp->c1min;  c1max = boxp->c1max;
 477   c2min = boxp->c2min;  c2max = boxp->c2max;
 478
 479   for (c0 = c0min; c0 <= c0max; c0++)
 480     for (c1 = c1min; c1 <= c1max; c1++) {
 481       histp = & histogram[c0][c1][c2min];
 482       for (c2 = c2min; c2 <= c2max; c2++) {
 483         if ((count = *histp++) != 0) {
 484           total += count;
 485           c0total += ((c0 << C0_SHIFT) + ((1<<C0_SHIFT)>>1)) * count;
 486           c1total += ((c1 << C1_SHIFT) + ((1<<C1_SHIFT)>>1)) * count;
 487           c2total += ((c2 << C2_SHIFT) + ((1<<C2_SHIFT)>>1)) * count;
 488         }
 489       }
 490     }
 491
 492   sl_colormap[0][icolor] = (JSAMPLE) ((c0total + (total>>1)) / total);
 493   sl_colormap[1][icolor] = (JSAMPLE) ((c1total + (total>>1)) / total);
 494   sl_colormap[2][icolor] = (JSAMPLE) ((c2total + (total>>1)) / total);
 495 }
 496
 497
 498 static void slow_select_colors (int descolors)
 499 /* Master routine for color selection */
 500 {
 501   box boxlist[MAXNUMCOLORS];
 502   int numboxes;
 503   int i;
 504
 505   /* Initialize one box containing whole space */
 506   numboxes = 1;
 507   boxlist[0].c0min = 0;
 508   boxlist[0].c0max = 255 >> C0_SHIFT;
 509   boxlist[0].c1min = 0;
 510   boxlist[0].c1max = 255 >> C1_SHIFT;
 511   boxlist[0].c2min = 0;
 512   boxlist[0].c2max = 255 >> C2_SHIFT;
 513   /* Shrink it to actually-used volume and set its statistics */
 514   update_box(& boxlist[0]);
 515   /* Perform median-cut to produce final box list */
 516   numboxes = median_cut(boxlist, numboxes, descolors);
 517   /* Compute the representative color for each box, fill colormap */
 518   for (i = 0; i < numboxes; i++)
 519     compute_color(& boxlist[i], i);
 520   sl_num_colors = numboxes;
 521 }
 522
 523
 524 /* log2(histogram cells in update box) for each axis; this can be adjusted */
 525 #define BOX_C0_LOG  (HIST_C0_BITS-3)
 526 #define BOX_C1_LOG  (HIST_C1_BITS-3)
 527 #define BOX_C2_LOG  (HIST_C2_BITS-3)
 528
 529 #define BOX_C0_ELEMS  (1<<BOX_C0_LOG) /* # of hist cells in update box */
 530 #define BOX_C1_ELEMS  (1<<BOX_C1_LOG)
 531 #define BOX_C2_ELEMS  (1<<BOX_C2_LOG)
 532
 533 #define BOX_C0_SHIFT  (C0_SHIFT + BOX_C0_LOG)
 534 #define BOX_C1_SHIFT  (C1_SHIFT + BOX_C1_LOG)
 535 #define BOX_C2_SHIFT  (C2_SHIFT + BOX_C2_LOG)
 536
 537
 538 static int find_nearby_colors (int minc0, int minc1, int minc2,
 539                                JSAMPLE colorlist[])
 540 {
 541   int numcolors = sl_num_colors;
 542   int maxc0, maxc1, maxc2;
 543   int centerc0, centerc1, centerc2;
 544   int i, x, ncolors;
 545   INT32 minmaxdist, min_dist, max_dist, tdist;
 546   INT32 mindist[MAXNUMCOLORS];        /* min distance to colormap entry i */
 547
 548   maxc0 = minc0 + ((1 << BOX_C0_SHIFT) - (1 << C0_SHIFT));
 549   centerc0 = (minc0 + maxc0) >> 1;
 550   maxc1 = minc1 + ((1 << BOX_C1_SHIFT) - (1 << C1_SHIFT));
 551   centerc1 = (minc1 + maxc1) >> 1;
 552   maxc2 = minc2 + ((1 << BOX_C2_SHIFT) - (1 << C2_SHIFT));
 553   centerc2 = (minc2 + maxc2) >> 1;
 554
 555   minmaxdist = 0x7FFFFFFFL;
 556
 557   for (i = 0; i < numcolors; i++) {
 558     /* We compute the squared-c0-distance term, then add in the other two. */
 559     x = sl_colormap[0][i];
 560     if (x < minc0) {
 561       tdist = (x - minc0) * C0_SCALE;
 562       min_dist = tdist*tdist;
 563       tdist = (x - maxc0) * C0_SCALE;
 564       max_dist = tdist*tdist;
 565     } else if (x > maxc0) {
 566       tdist = (x - maxc0) * C0_SCALE;
 567       min_dist = tdist*tdist;
 568       tdist = (x - minc0) * C0_SCALE;
 569       max_dist = tdist*tdist;
 570     } else {
 571       /* within cell range so no contribution to min_dist */
 572       min_dist = 0;
 573       if (x <= centerc0) {
 574         tdist = (x - maxc0) * C0_SCALE;
 575         max_dist = tdist*tdist;
 576       } else {
 577         tdist = (x - minc0) * C0_SCALE;
 578         max_dist = tdist*tdist;
 579       }
 580     }
 581
 582     x = sl_colormap[1][i];
 583     if (x < minc1) {
 584       tdist = (x - minc1) * C1_SCALE;
 585       min_dist += tdist*tdist;
 586       tdist = (x - maxc1) * C1_SCALE;
 587       max_dist += tdist*tdist;
 588     } else if (x > maxc1) {
 589       tdist = (x - maxc1) * C1_SCALE;
 590       min_dist += tdist*tdist;
 591       tdist = (x - minc1) * C1_SCALE;
 592       max_dist += tdist*tdist;
 593     } else {
 594       /* within cell range so no contribution to min_dist */
 595       if (x <= centerc1) {
 596         tdist = (x - maxc1) * C1_SCALE;
 597         max_dist += tdist*tdist;
 598       } else {
 599         tdist = (x - minc1) * C1_SCALE;
 600         max_dist += tdist*tdist;
 601       }
 602     }
 603
 604     x = sl_colormap[2][i];
 605     if (x < minc2) {
 606       tdist = (x - minc2) * C2_SCALE;
 607       min_dist += tdist*tdist;
 608       tdist = (x - maxc2) * C2_SCALE;
 609       max_dist += tdist*tdist;
 610     } else if (x > maxc2) {
 611       tdist = (x - maxc2) * C2_SCALE;
 612       min_dist += tdist*tdist;
 613       tdist = (x - minc2) * C2_SCALE;
 614       max_dist += tdist*tdist;
 615     } else {
 616       /* within cell range so no contribution to min_dist */
 617       if (x <= centerc2) {
 618         tdist = (x - maxc2) * C2_SCALE;
 619         max_dist += tdist*tdist;
 620       } else {
 621         tdist = (x - minc2) * C2_SCALE;
 622         max_dist += tdist*tdist;
 623       }
 624     }
 625
 626     mindist[i] = min_dist;        /* save away the results */
 627     if (max_dist < minmaxdist)
 628       minmaxdist = max_dist;
 629   }
 630
 631   ncolors = 0;
 632   for (i = 0; i < numcolors; i++) {
 633     if (mindist[i] <= minmaxdist)
 634       colorlist[ncolors++] = (JSAMPLE) i;
 635   }
 636   return ncolors;
 637 }
 638
 639
 640 static void find_best_colors (int minc0,int minc1,int minc2,int numcolors,
 641                               JSAMPLE colorlist[], JSAMPLE bestcolor[])
 642 {
 643   int ic0, ic1, ic2;
 644   int i, icolor;
 645   register INT32 * bptr;        /* pointer into bestdist[] array */
 646   JSAMPLE * cptr;                /* pointer into bestcolor[] array */
 647   INT32 dist0, dist1;                /* initial distance values */
 648   register INT32 dist2;                /* current distance in inner loop */
 649   INT32 xx0, xx1;                /* distance increments */
 650   register INT32 xx2;
 651   INT32 inc0, inc1, inc2;        /* initial values for increments */
 652   /* This array holds the distance to the nearest-so-far color for each cell */
 653   INT32 bestdist[BOX_C0_ELEMS * BOX_C1_ELEMS * BOX_C2_ELEMS];
 654
 655   /* Initialize best-distance for each cell of the update box */
 656   bptr = bestdist;
 657   for (i = BOX_C0_ELEMS*BOX_C1_ELEMS*BOX_C2_ELEMS-1; i >= 0; i--)
 658     *bptr++ = 0x7FFFFFFFL;
 659
 660   /* Nominal steps between cell centers ("x" in Thomas article) */
 661 #define STEP_C0  ((1 << C0_SHIFT) * C0_SCALE)
 662 #define STEP_C1  ((1 << C1_SHIFT) * C1_SCALE)
 663 #define STEP_C2  ((1 << C2_SHIFT) * C2_SCALE)
 664
 665   for (i = 0; i < numcolors; i++) {
 666     icolor = colorlist[i];
 667     /* Compute (square of) distance from minc0/c1/c2 to this color */
 668     inc0 = (minc0 - (int) sl_colormap[0][icolor]) * C0_SCALE;
 669     dist0 = inc0*inc0;
 670     inc1 = (minc1 - (int) sl_colormap[1][icolor]) * C1_SCALE;
 671     dist0 += inc1*inc1;
 672     inc2 = (minc2 - (int) sl_colormap[2][icolor]) * C2_SCALE;
 673     dist0 += inc2*inc2;
 674     /* Form the initial difference increments */
 675     inc0 = inc0 * (2 * STEP_C0) + STEP_C0 * STEP_C0;
 676     inc1 = inc1 * (2 * STEP_C1) + STEP_C1 * STEP_C1;
 677     inc2 = inc2 * (2 * STEP_C2) + STEP_C2 * STEP_C2;
 678     /* Now loop over all cells in box, updating distance per Thomas method */
 679     bptr = bestdist;
 680     cptr = bestcolor;
 681     xx0 = inc0;
 682     for (ic0 = BOX_C0_ELEMS-1; ic0 >= 0; ic0--) {
 683       dist1 = dist0;
 684       xx1 = inc1;
 685       for (ic1 = BOX_C1_ELEMS-1; ic1 >= 0; ic1--) {
 686         dist2 = dist1;
 687         xx2 = inc2;
 688         for (ic2 = BOX_C2_ELEMS-1; ic2 >= 0; ic2--) {
 689           if (dist2 < *bptr) {
 690             *bptr = dist2;
 691             *cptr = (JSAMPLE) icolor;
 692           }
 693           dist2 += xx2;
 694           xx2 += 2 * STEP_C2 * STEP_C2;
 695           bptr++;
 696           cptr++;
 697         }
 698         dist1 += xx1;
 699         xx1 += 2 * STEP_C1 * STEP_C1;
 700       }
 701       dist0 += xx0;
 702       xx0 += 2 * STEP_C0 * STEP_C0;
 703     }
 704   }
 705 }
 706
 707
 708 static void fill_inverse_cmap (int c0, int c1, int c2)
 709 {
 710   hist2d * histogram = sl_histogram;
 711   int minc0, minc1, minc2;        /* lower left corner of update box */
 712   int ic0, ic1, ic2;
 713   register JSAMPLE * cptr;        /* pointer into bestcolor[] array */
 714   register histptr cachep;        /* pointer into main cache array */
 715   /* This array lists the candidate colormap indexes. */
 716   JSAMPLE colorlist[MAXNUMCOLORS];
 717   int numcolors;                /* number of candidate colors */
 718   /* This array holds the actually closest colormap index for each cell. */
 719   JSAMPLE bestcolor[BOX_C0_ELEMS * BOX_C1_ELEMS * BOX_C2_ELEMS];
 720
 721   /* Convert cell coordinates to update box ID */
 722   c0 >>= BOX_C0_LOG;
 723   c1 >>= BOX_C1_LOG;
 724   c2 >>= BOX_C2_LOG;
 725
 726   minc0 = (c0 << BOX_C0_SHIFT) + ((1 << C0_SHIFT) >> 1);
 727   minc1 = (c1 << BOX_C1_SHIFT) + ((1 << C1_SHIFT) >> 1);
 728   minc2 = (c2 << BOX_C2_SHIFT) + ((1 << C2_SHIFT) >> 1);
 729
 730   numcolors = find_nearby_colors(minc0, minc1, minc2, colorlist);
 731
 732   /* Determine the actually nearest colors. */
 733   find_best_colors(minc0, minc1, minc2, numcolors, colorlist, bestcolor);
 734
 735   /* Save the best color numbers (plus 1) in the main cache array */
 736   c0 <<= BOX_C0_LOG;                /* convert ID back to base cell indexes */
 737   c1 <<= BOX_C1_LOG;
 738   c2 <<= BOX_C2_LOG;
 739   cptr = bestcolor;
 740   for (ic0 = 0; ic0 < BOX_C0_ELEMS; ic0++) {
 741     for (ic1 = 0; ic1 < BOX_C1_ELEMS; ic1++) {
 742       cachep = & histogram[c0+ic0][c1+ic1][c2];
 743       for (ic2 = 0; ic2 < BOX_C2_ELEMS; ic2++) {
 744         *cachep++ = (histcell) (*cptr++ + 1);
 745       }
 746     }
 747   }
 748 }
 749
 750 static void slow_map_pixels (byte *pic24,int width,int height, byte* pic8)
 751 {
 752   register LOCFSERROR cur0, cur1, cur2;        /* current error or pixel value */
 753   LOCFSERROR belowerr0, belowerr1, belowerr2; /* error for pixel below cur */
 754   LOCFSERROR bpreverr0, bpreverr1, bpreverr2; /* error for below/prev col */
 755   register FSERRPTR errorptr;        /* => fserrors[] at column before current */
 756   JSAMPROW inptr;                /* => current input pixel */
 757   JSAMPROW outptr;                /* => current output pixel */
 758   histptr cachep;
 759   int dir;                        /* +1 or -1 depending on direction */
 760   int dir3;                        /* 3*dir, for advancing inptr & errorptr */
 761   int row, col;
 762   int *error_limit = sl_error_limiter;
 763   JSAMPROW colormap0 = sl_colormap[0];
 764   JSAMPROW colormap1 = sl_colormap[1];
 765   JSAMPROW colormap2 = sl_colormap[2];
 766   hist2d * histogram = sl_histogram;
 767
 768   for (row = 0; row < height; row++) {
 769
 770 /*    if ((row&0x3f) == 0) WaitCursor();
 771     ProgressMeter(0, height-1, row, "Dither");*/
 772
 773     inptr = & pic24[row * width * 3];
 774     outptr = & pic8[row * width];
 775     if (sl_on_odd_row) {
 776       /* work right to left in this row */
 777       inptr += (width-1) * 3;        /* so point to rightmost pixel */
 778       outptr += width-1;
 779       dir = -1;
 780       dir3 = -3;
 781       errorptr = sl_fserrors + (width+1)*3; /* => entry after last column */
 782       sl_on_odd_row = FALSE;        /* flip for next time */
 783     } else {
 784       /* work left to right in this row */
 785       dir = 1;
 786       dir3 = 3;
 787       errorptr = sl_fserrors;        /* => entry before first real column */
 788       sl_on_odd_row = TRUE;        /* flip for next time */
 789     }
 790     /* Preset error values: no error propagated to first pixel from left */
 791     cur0 = cur1 = cur2 = 0;
 792     /* and no error propagated to row below yet */
 793     belowerr0 = belowerr1 = belowerr2 = 0;
 794     bpreverr0 = bpreverr1 = bpreverr2 = 0;
 795
 796     for (col = width; col > 0; col--)
 797     {
 798       cur0 = (cur0 + errorptr[dir3+0] + 8) >> 4;
 799       cur1 = (cur1 + errorptr[dir3+1] + 8) >> 4;
 800       cur2 = (cur2 + errorptr[dir3+2] + 8) >> 4;
 801       cur0 = error_limit[cur0];
 802       cur1 = error_limit[cur1];
 803       cur2 = error_limit[cur2];
 804       cur0 += inptr[0];
 805       cur1 += inptr[1];
 806       cur2 += inptr[2];
 807       RANGE(cur0, 0, 255);
 808       RANGE(cur1, 0, 255);
 809       RANGE(cur2, 0, 255);
 810       /* Index into the cache with adjusted pixel value */
 811       cachep = & histogram[cur0>>C0_SHIFT][cur1>>C1_SHIFT][cur2>>C2_SHIFT];
 812       /* If we have not seen this color before, find nearest colormap */
 813       /* entry and update the cache */
 814       if (*cachep == 0)
 815         fill_inverse_cmap(cur0>>C0_SHIFT, cur1>>C1_SHIFT, cur2>>C2_SHIFT);
 816       /* Now emit the colormap index for this cell */
 817       { register int pixcode = *cachep - 1;
 818         *outptr = (JSAMPLE) pixcode;
 819         /* Compute representation error for this pixel */
 820         cur0 -= (int) colormap0[pixcode];
 821         cur1 -= (int) colormap1[pixcode];
 822         cur2 -= (int) colormap2[pixcode];
 823       }
 824       /* Compute error fractions to be propagated to adjacent pixels.
 825        * Add these into the running sums, and simultaneously shift the
 826        * next-line error sums left by 1 column.
 827        */
 828       { register LOCFSERROR bnexterr, delta;
 829
 830         bnexterr = cur0;        /* Process component 0 */
 831         delta = cur0 * 2;
 832         cur0 += delta;                /* form error * 3 */
 833         errorptr[0] = (FSERROR) (bpreverr0 + cur0);
 834         cur0 += delta;                /* form error * 5 */
 835         bpreverr0 = belowerr0 + cur0;
 836         belowerr0 = bnexterr;
 837         cur0 += delta;                /* form error * 7 */
 838         bnexterr = cur1;        /* Process component 1 */
 839         delta = cur1 * 2;
 840         cur1 += delta;                /* form error * 3 */
 841         errorptr[1] = (FSERROR) (bpreverr1 + cur1);
 842         cur1 += delta;                /* form error * 5 */
 843         bpreverr1 = belowerr1 + cur1;
 844         belowerr1 = bnexterr;
 845         cur1 += delta;                /* form error * 7 */
 846         bnexterr = cur2;        /* Process component 2 */
 847         delta = cur2 * 2;
 848         cur2 += delta;                /* form error * 3 */
 849         errorptr[2] = (FSERROR) (bpreverr2 + cur2);
 850         cur2 += delta;                /* form error * 5 */
 851         bpreverr2 = belowerr2 + cur2;
 852         belowerr2 = bnexterr;
 853         cur2 += delta;                /* form error * 7 */
 854       }
 855       /* At this point curN contains the 7/16 error value to be propagated
 856        * to the next pixel on the current line, and all the errors for the
 857        * next line have been shifted over.  We are therefore ready to move on.
 858        */
 859       inptr += dir3;                /* Advance pixel pointers to next column */
 860       outptr += dir;
 861       errorptr += dir3;                /* advance errorptr to current column */
 862     }
 863     /* Post-loop cleanup: we must unload the final error values into the
 864      * final fserrors[] entry.  Note we need not unload belowerrN because
 865      * it is for the dummy column before or after the actual array.
 866      */
 867     errorptr[0] = (FSERROR) bpreverr0; /* unload prev errs into array */
 868     errorptr[1] = (FSERROR) bpreverr1;
 869     errorptr[2] = (FSERROR) bpreverr2;
 870   }
 871 }
 872
 873
 874 static void init_error_limit ()
 875 /* Allocate and fill in the error_limiter table */
 876 /* Note this should be done only once. */
 877 {
 878   int * table;
 879   int in, out;
 880
 881   table = (int *) malloc((size_t) ((255*2+1) * sizeof(int)));
 882   if (! table) return;
 883
 884   table += 255;                /* so can index -255 .. +255 */
 885   sl_error_limiter = table;
 886
 887 #define STEPSIZE ((255+1)/16)
 888   /* Map errors 1:1 up to +- 255/16 */
 889   out = 0;
 890   for (in = 0; in < STEPSIZE; in++, out++) {
 891     table[in] = out; table[-in] = -out;
 892   }
 893   /* Map errors 1:2 up to +- 3*255/16 */
 894   for (; in < STEPSIZE*3; in++, out += (in&1) ? 0 : 1) {
 895     table[in] = out; table[-in] = -out;
 896   }
 897   /* Clamp the rest to final out value (which is (255+1)/8) */
 898   for (; in <= 255; in++) {
 899     table[in] = out; table[-in] = -out;
 900   }
 901 #undef STEPSIZE
 902 }