Home | History | Annotate | Download | only in libjpeg_turbo
      1 /*
      2  * jquant2.c
      3  *
      4  * This file was part of the Independent JPEG Group's software:
      5  * Copyright (C) 1991-1996, Thomas G. Lane.
      6  * libjpeg-turbo Modifications:
      7  * Copyright (C) 2009, D. R. Commander.
      8  * For conditions of distribution and use, see the accompanying README file.
      9  *
     10  * This file contains 2-pass color quantization (color mapping) routines.
     11  * These routines provide selection of a custom color map for an image,
     12  * followed by mapping of the image to that color map, with optional
     13  * Floyd-Steinberg dithering.
     14  * It is also possible to use just the second pass to map to an arbitrary
     15  * externally-given color map.
     16  *
     17  * Note: ordered dithering is not supported, since there isn't any fast
     18  * way to compute intercolor distances; it's unclear that ordered dither's
     19  * fundamental assumptions even hold with an irregularly spaced color map.
     20  */
     21 
     22 #define JPEG_INTERNALS
     23 #include "jinclude.h"
     24 #include "jpeglib.h"
     25 
     26 #ifdef QUANT_2PASS_SUPPORTED
     27 
     28 
     29 /*
     30  * This module implements the well-known Heckbert paradigm for color
     31  * quantization.  Most of the ideas used here can be traced back to
     32  * Heckbert's seminal paper
     33  *   Heckbert, Paul.  "Color Image Quantization for Frame Buffer Display",
     34  *   Proc. SIGGRAPH '82, Computer Graphics v.16 #3 (July 1982), pp 297-304.
     35  *
     36  * In the first pass over the image, we accumulate a histogram showing the
     37  * usage count of each possible color.  To keep the histogram to a reasonable
     38  * size, we reduce the precision of the input; typical practice is to retain
     39  * 5 or 6 bits per color, so that 8 or 4 different input values are counted
     40  * in the same histogram cell.
     41  *
     42  * Next, the color-selection step begins with a box representing the whole
     43  * color space, and repeatedly splits the "largest" remaining box until we
     44  * have as many boxes as desired colors.  Then the mean color in each
     45  * remaining box becomes one of the possible output colors.
     46  *
     47  * The second pass over the image maps each input pixel to the closest output
     48  * color (optionally after applying a Floyd-Steinberg dithering correction).
     49  * This mapping is logically trivial, but making it go fast enough requires
     50  * considerable care.
     51  *
     52  * Heckbert-style quantizers vary a good deal in their policies for choosing
     53  * the "largest" box and deciding where to cut it.  The particular policies
     54  * used here have proved out well in experimental comparisons, but better ones
     55  * may yet be found.
     56  *
     57  * In earlier versions of the IJG code, this module quantized in YCbCr color
     58  * space, processing the raw upsampled data without a color conversion step.
     59  * This allowed the color conversion math to be done only once per colormap
     60  * entry, not once per pixel.  However, that optimization precluded other
     61  * useful optimizations (such as merging color conversion with upsampling)
     62  * and it also interfered with desired capabilities such as quantizing to an
     63  * externally-supplied colormap.  We have therefore abandoned that approach.
     64  * The present code works in the post-conversion color space, typically RGB.
     65  *
     66  * To improve the visual quality of the results, we actually work in scaled
     67  * RGB space, giving G distances more weight than R, and R in turn more than
     68  * B.  To do everything in integer math, we must use integer scale factors.
     69  * The 2/3/1 scale factors used here correspond loosely to the relative
     70  * weights of the colors in the NTSC grayscale equation.
     71  * If you want to use this code to quantize a non-RGB color space, you'll
     72  * probably need to change these scale factors.
     73  */
     74 
     75 #define R_SCALE 2		/* scale R distances by this much */
     76 #define G_SCALE 3		/* scale G distances by this much */
     77 #define B_SCALE 1		/* and B by this much */
     78 
     79 static const int c_scales[3]={R_SCALE, G_SCALE, B_SCALE};
     80 #define C0_SCALE c_scales[rgb_red[cinfo->out_color_space]]
     81 #define C1_SCALE c_scales[rgb_green[cinfo->out_color_space]]
     82 #define C2_SCALE c_scales[rgb_blue[cinfo->out_color_space]]
     83 
     84 /*
     85  * First we have the histogram data structure and routines for creating it.
     86  *
     87  * The number of bits of precision can be adjusted by changing these symbols.
     88  * We recommend keeping 6 bits for G and 5 each for R and B.
     89  * If you have plenty of memory and cycles, 6 bits all around gives marginally
     90  * better results; if you are short of memory, 5 bits all around will save
     91  * some space but degrade the results.
     92  * To maintain a fully accurate histogram, we'd need to allocate a "long"
     93  * (preferably unsigned long) for each cell.  In practice this is overkill;
     94  * we can get by with 16 bits per cell.  Few of the cell counts will overflow,
     95  * and clamping those that do overflow to the maximum value will give close-
     96  * enough results.  This reduces the recommended histogram size from 256Kb
     97  * to 128Kb, which is a useful savings on PC-class machines.
     98  * (In the second pass the histogram space is re-used for pixel mapping data;
     99  * in that capacity, each cell must be able to store zero to the number of
    100  * desired colors.  16 bits/cell is plenty for that too.)
    101  * Since the JPEG code is intended to run in small memory model on 80x86
    102  * machines, we can't just allocate the histogram in one chunk.  Instead
    103  * of a true 3-D array, we use a row of pointers to 2-D arrays.  Each
    104  * pointer corresponds to a C0 value (typically 2^5 = 32 pointers) and
    105  * each 2-D array has 2^6*2^5 = 2048 or 2^6*2^6 = 4096 entries.  Note that
    106  * on 80x86 machines, the pointer row is in near memory but the actual
    107  * arrays are in far memory (same arrangement as we use for image arrays).
    108  */
    109 
    110 #define MAXNUMCOLORS  (MAXJSAMPLE+1) /* maximum size of colormap */
    111 
    112 /* These will do the right thing for either R,G,B or B,G,R color order,
    113  * but you may not like the results for other color orders.
    114  */
    115 #define HIST_C0_BITS  5		/* bits of precision in R/B histogram */
    116 #define HIST_C1_BITS  6		/* bits of precision in G histogram */
    117 #define HIST_C2_BITS  5		/* bits of precision in B/R histogram */
    118 
    119 /* Number of elements along histogram axes. */
    120 #define HIST_C0_ELEMS  (1<<HIST_C0_BITS)
    121 #define HIST_C1_ELEMS  (1<<HIST_C1_BITS)
    122 #define HIST_C2_ELEMS  (1<<HIST_C2_BITS)
    123 
    124 /* These are the amounts to shift an input value to get a histogram index. */
    125 #define C0_SHIFT  (BITS_IN_JSAMPLE-HIST_C0_BITS)
    126 #define C1_SHIFT  (BITS_IN_JSAMPLE-HIST_C1_BITS)
    127 #define C2_SHIFT  (BITS_IN_JSAMPLE-HIST_C2_BITS)
    128 
    129 
    130 typedef UINT16 histcell;	/* histogram cell; prefer an unsigned type */
    131 
    132 typedef histcell FAR * histptr;	/* for pointers to histogram cells */
    133 
    134 typedef histcell hist1d[HIST_C2_ELEMS]; /* typedefs for the array */
    135 typedef hist1d FAR * hist2d;	/* type for the 2nd-level pointers */
    136 typedef hist2d * hist3d;	/* type for top-level pointer */
    137 
    138 
    139 /* Declarations for Floyd-Steinberg dithering.
    140  *
    141  * Errors are accumulated into the array fserrors[], at a resolution of
    142  * 1/16th of a pixel count.  The error at a given pixel is propagated
    143  * to its not-yet-processed neighbors using the standard F-S fractions,
    144  *		...	(here)	7/16
    145  *		3/16	5/16	1/16
    146  * We work left-to-right on even rows, right-to-left on odd rows.
    147  *
    148  * We can get away with a single array (holding one row's worth of errors)
    149  * by using it to store the current row's errors at pixel columns not yet
    150  * processed, but the next row's errors at columns already processed.  We
    151  * need only a few extra variables to hold the errors immediately around the
    152  * current column.  (If we are lucky, those variables are in registers, but
    153  * even if not, they're probably cheaper to access than array elements are.)
    154  *
    155  * The fserrors[] array has (#columns + 2) entries; the extra entry at
    156  * each end saves us from special-casing the first and last pixels.
    157  * Each entry is three values long, one value for each color component.
    158  *
    159  * Note: on a wide image, we might not have enough room in a PC's near data
    160  * segment to hold the error array; so it is allocated with alloc_large.
    161  */
    162 
    163 #if BITS_IN_JSAMPLE == 8
    164 typedef INT16 FSERROR;		/* 16 bits should be enough */
    165 typedef int LOCFSERROR;		/* use 'int' for calculation temps */
    166 #else
    167 typedef INT32 FSERROR;		/* may need more than 16 bits */
    168 typedef INT32 LOCFSERROR;	/* be sure calculation temps are big enough */
    169 #endif
    170 
    171 typedef FSERROR FAR *FSERRPTR;	/* pointer to error array (in FAR storage!) */
    172 
    173 
    174 /* Private subobject */
    175 
    176 typedef struct {
    177   struct jpeg_color_quantizer pub; /* public fields */
    178 
    179   /* Space for the eventually created colormap is stashed here */
    180   JSAMPARRAY sv_colormap;	/* colormap allocated at init time */
    181   int desired;			/* desired # of colors = size of colormap */
    182 
    183   /* Variables for accumulating image statistics */
    184   hist3d histogram;		/* pointer to the histogram */
    185 
    186   boolean needs_zeroed;		/* TRUE if next pass must zero histogram */
    187 
    188   /* Variables for Floyd-Steinberg dithering */
    189   FSERRPTR fserrors;		/* accumulated errors */
    190   boolean on_odd_row;		/* flag to remember which row we are on */
    191   int * error_limiter;		/* table for clamping the applied error */
    192 } my_cquantizer;
    193 
    194 typedef my_cquantizer * my_cquantize_ptr;
    195 
    196 
    197 /*
    198  * Prescan some rows of pixels.
    199  * In this module the prescan simply updates the histogram, which has been
    200  * initialized to zeroes by start_pass.
    201  * An output_buf parameter is required by the method signature, but no data
    202  * is actually output (in fact the buffer controller is probably passing a
    203  * NULL pointer).
    204  */
    205 
    206 METHODDEF(void)
    207 prescan_quantize (j_decompress_ptr cinfo, JSAMPARRAY input_buf,
    208 		  JSAMPARRAY output_buf, int num_rows)
    209 {
    210   my_cquantize_ptr cquantize = (my_cquantize_ptr) cinfo->cquantize;
    211   register JSAMPROW ptr;
    212   register histptr histp;
    213   register hist3d histogram = cquantize->histogram;
    214   int row;
    215   JDIMENSION col;
    216   JDIMENSION width = cinfo->output_width;
    217 
    218   for (row = 0; row < num_rows; row++) {
    219     ptr = input_buf[row];
    220     for (col = width; col > 0; col--) {
    221       /* get pixel value and index into the histogram */
    222       histp = & histogram[GETJSAMPLE(ptr[0]) >> C0_SHIFT]
    223 			 [GETJSAMPLE(ptr[1]) >> C1_SHIFT]
    224 			 [GETJSAMPLE(ptr[2]) >> C2_SHIFT];
    225       /* increment, check for overflow and undo increment if so. */
    226       if (++(*histp) <= 0)
    227 	(*histp)--;
    228       ptr += 3;
    229     }
    230   }
    231 }
    232 
    233 
    234 /*
    235  * Next we have the really interesting routines: selection of a colormap
    236  * given the completed histogram.
    237  * These routines work with a list of "boxes", each representing a rectangular
    238  * subset of the input color space (to histogram precision).
    239  */
    240 
    241 typedef struct {
    242   /* The bounds of the box (inclusive); expressed as histogram indexes */
    243   int c0min, c0max;
    244   int c1min, c1max;
    245   int c2min, c2max;
    246   /* The volume (actually 2-norm) of the box */
    247   INT32 volume;
    248   /* The number of nonzero histogram cells within this box */
    249   long colorcount;
    250 } box;
    251 
    252 typedef box * boxptr;
    253 
    254 
    255 LOCAL(boxptr)
    256 find_biggest_color_pop (boxptr boxlist, int numboxes)
    257 /* Find the splittable box with the largest color population */
    258 /* Returns NULL if no splittable boxes remain */
    259 {
    260   register boxptr boxp;
    261   register int i;
    262   register long maxc = 0;
    263   boxptr which = NULL;
    264 
    265   for (i = 0, boxp = boxlist; i < numboxes; i++, boxp++) {
    266     if (boxp->colorcount > maxc && boxp->volume > 0) {
    267       which = boxp;
    268       maxc = boxp->colorcount;
    269     }
    270   }
    271   return which;
    272 }
    273 
    274 
    275 LOCAL(boxptr)
    276 find_biggest_volume (boxptr boxlist, int numboxes)
    277 /* Find the splittable box with the largest (scaled) volume */
    278 /* Returns NULL if no splittable boxes remain */
    279 {
    280   register boxptr boxp;
    281   register int i;
    282   register INT32 maxv = 0;
    283   boxptr which = NULL;
    284 
    285   for (i = 0, boxp = boxlist; i < numboxes; i++, boxp++) {
    286     if (boxp->volume > maxv) {
    287       which = boxp;
    288       maxv = boxp->volume;
    289     }
    290   }
    291   return which;
    292 }
    293 
    294 
    295 LOCAL(void)
    296 update_box (j_decompress_ptr cinfo, boxptr boxp)
    297 /* Shrink the min/max bounds of a box to enclose only nonzero elements, */
    298 /* and recompute its volume and population */
    299 {
    300   my_cquantize_ptr cquantize = (my_cquantize_ptr) cinfo->cquantize;
    301   hist3d histogram = cquantize->histogram;
    302   histptr histp;
    303   int c0,c1,c2;
    304   int c0min,c0max,c1min,c1max,c2min,c2max;
    305   INT32 dist0,dist1,dist2;
    306   long ccount;
    307 
    308   c0min = boxp->c0min;  c0max = boxp->c0max;
    309   c1min = boxp->c1min;  c1max = boxp->c1max;
    310   c2min = boxp->c2min;  c2max = boxp->c2max;
    311 
    312   if (c0max > c0min)
    313     for (c0 = c0min; c0 <= c0max; c0++)
    314       for (c1 = c1min; c1 <= c1max; c1++) {
    315 	histp = & histogram[c0][c1][c2min];
    316 	for (c2 = c2min; c2 <= c2max; c2++)
    317 	  if (*histp++ != 0) {
    318 	    boxp->c0min = c0min = c0;
    319 	    goto have_c0min;
    320 	  }
    321       }
    322  have_c0min:
    323   if (c0max > c0min)
    324     for (c0 = c0max; c0 >= c0min; c0--)
    325       for (c1 = c1min; c1 <= c1max; c1++) {
    326 	histp = & histogram[c0][c1][c2min];
    327 	for (c2 = c2min; c2 <= c2max; c2++)
    328 	  if (*histp++ != 0) {
    329 	    boxp->c0max = c0max = c0;
    330 	    goto have_c0max;
    331 	  }
    332       }
    333  have_c0max:
    334   if (c1max > c1min)
    335     for (c1 = c1min; c1 <= c1max; c1++)
    336       for (c0 = c0min; c0 <= c0max; c0++) {
    337 	histp = & histogram[c0][c1][c2min];
    338 	for (c2 = c2min; c2 <= c2max; c2++)
    339 	  if (*histp++ != 0) {
    340 	    boxp->c1min = c1min = c1;
    341 	    goto have_c1min;
    342 	  }
    343       }
    344  have_c1min:
    345   if (c1max > c1min)
    346     for (c1 = c1max; c1 >= c1min; c1--)
    347       for (c0 = c0min; c0 <= c0max; c0++) {
    348 	histp = & histogram[c0][c1][c2min];
    349 	for (c2 = c2min; c2 <= c2max; c2++)
    350 	  if (*histp++ != 0) {
    351 	    boxp->c1max = c1max = c1;
    352 	    goto have_c1max;
    353 	  }
    354       }
    355  have_c1max:
    356   if (c2max > c2min)
    357     for (c2 = c2min; c2 <= c2max; c2++)
    358       for (c0 = c0min; c0 <= c0max; c0++) {
    359 	histp = & histogram[c0][c1min][c2];
    360 	for (c1 = c1min; c1 <= c1max; c1++, histp += HIST_C2_ELEMS)
    361 	  if (*histp != 0) {
    362 	    boxp->c2min = c2min = c2;
    363 	    goto have_c2min;
    364 	  }
    365       }
    366  have_c2min:
    367   if (c2max > c2min)
    368     for (c2 = c2max; c2 >= c2min; c2--)
    369       for (c0 = c0min; c0 <= c0max; c0++) {
    370 	histp = & histogram[c0][c1min][c2];
    371 	for (c1 = c1min; c1 <= c1max; c1++, histp += HIST_C2_ELEMS)
    372 	  if (*histp != 0) {
    373 	    boxp->c2max = c2max = c2;
    374 	    goto have_c2max;
    375 	  }
    376       }
    377  have_c2max:
    378 
    379   /* Update box volume.
    380    * We use 2-norm rather than real volume here; this biases the method
    381    * against making long narrow boxes, and it has the side benefit that
    382    * a box is splittable iff norm > 0.
    383    * Since the differences are expressed in histogram-cell units,
    384    * we have to shift back to JSAMPLE units to get consistent distances;
    385    * after which, we scale according to the selected distance scale factors.
    386    */
    387   dist0 = ((c0max - c0min) << C0_SHIFT) * C0_SCALE;
    388   dist1 = ((c1max - c1min) << C1_SHIFT) * C1_SCALE;
    389   dist2 = ((c2max - c2min) << C2_SHIFT) * C2_SCALE;
    390   boxp->volume = dist0*dist0 + dist1*dist1 + dist2*dist2;
    391 
    392   /* Now scan remaining volume of box and compute population */
    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 LOCAL(int)
    407 median_cut (j_decompress_ptr cinfo, boxptr boxlist, int numboxes,
    408 	    int desired_colors)
    409 /* Repeatedly select and split the largest box until we have enough boxes */
    410 {
    411   int n,lb;
    412   int c0,c1,c2,cmax;
    413   register boxptr b1,b2;
    414 
    415   while (numboxes < desired_colors) {
    416     /* Select box to split.
    417      * Current algorithm: by population for first half, then by volume.
    418      */
    419     if (numboxes*2 <= desired_colors) {
    420       b1 = find_biggest_color_pop(boxlist, numboxes);
    421     } else {
    422       b1 = find_biggest_volume(boxlist, numboxes);
    423     }
    424     if (b1 == NULL)		/* no splittable boxes left! */
    425       break;
    426     b2 = &boxlist[numboxes];	/* where new box will go */
    427     /* Copy the color bounds to the new box. */
    428     b2->c0max = b1->c0max; b2->c1max = b1->c1max; b2->c2max = b1->c2max;
    429     b2->c0min = b1->c0min; b2->c1min = b1->c1min; b2->c2min = b1->c2min;
    430     /* Choose which axis to split the box on.
    431      * Current algorithm: longest scaled axis.
    432      * See notes in update_box about scaling distances.
    433      */
    434     c0 = ((b1->c0max - b1->c0min) << C0_SHIFT) * C0_SCALE;
    435     c1 = ((b1->c1max - b1->c1min) << C1_SHIFT) * C1_SCALE;
    436     c2 = ((b1->c2max - b1->c2min) << C2_SHIFT) * C2_SCALE;
    437     /* We want to break any ties in favor of green, then red, blue last.
    438      * This code does the right thing for R,G,B or B,G,R color orders only.
    439      */
    440     if (rgb_red[cinfo->out_color_space] == 0) {
    441       cmax = c1; n = 1;
    442       if (c0 > cmax) { cmax = c0; n = 0; }
    443       if (c2 > cmax) { n = 2; }
    444     }
    445     else {
    446       cmax = c1; n = 1;
    447       if (c2 > cmax) { cmax = c2; n = 2; }
    448       if (c0 > cmax) { n = 0; }
    449     }
    450     /* Choose split point along selected axis, and update box bounds.
    451      * Current algorithm: split at halfway point.
    452      * (Since the box has been shrunk to minimum volume,
    453      * any split will produce two nonempty subboxes.)
    454      * Note that lb value is max for lower box, so must be < old max.
    455      */
    456     switch (n) {
    457     case 0:
    458       lb = (b1->c0max + b1->c0min) / 2;
    459       b1->c0max = lb;
    460       b2->c0min = lb+1;
    461       break;
    462     case 1:
    463       lb = (b1->c1max + b1->c1min) / 2;
    464       b1->c1max = lb;
    465       b2->c1min = lb+1;
    466       break;
    467     case 2:
    468       lb = (b1->c2max + b1->c2min) / 2;
    469       b1->c2max = lb;
    470       b2->c2min = lb+1;
    471       break;
    472     }
    473     /* Update stats for boxes */
    474     update_box(cinfo, b1);
    475     update_box(cinfo, b2);
    476     numboxes++;
    477   }
    478   return numboxes;
    479 }
    480 
    481 
    482 LOCAL(void)
    483 compute_color (j_decompress_ptr cinfo, boxptr boxp, int icolor)
    484 /* Compute representative color for a box, put it in colormap[icolor] */
    485 {
    486   /* Current algorithm: mean weighted by pixels (not colors) */
    487   /* Note it is important to get the rounding correct! */
    488   my_cquantize_ptr cquantize = (my_cquantize_ptr) cinfo->cquantize;
    489   hist3d histogram = cquantize->histogram;
    490   histptr histp;
    491   int c0,c1,c2;
    492   int c0min,c0max,c1min,c1max,c2min,c2max;
    493   long count;
    494   long total = 0;
    495   long c0total = 0;
    496   long c1total = 0;
    497   long c2total = 0;
    498 
    499   c0min = boxp->c0min;  c0max = boxp->c0max;
    500   c1min = boxp->c1min;  c1max = boxp->c1max;
    501   c2min = boxp->c2min;  c2max = boxp->c2max;
    502 
    503   for (c0 = c0min; c0 <= c0max; c0++)
    504     for (c1 = c1min; c1 <= c1max; c1++) {
    505       histp = & histogram[c0][c1][c2min];
    506       for (c2 = c2min; c2 <= c2max; c2++) {
    507 	if ((count = *histp++) != 0) {
    508 	  total += count;
    509 	  c0total += ((c0 << C0_SHIFT) + ((1<<C0_SHIFT)>>1)) * count;
    510 	  c1total += ((c1 << C1_SHIFT) + ((1<<C1_SHIFT)>>1)) * count;
    511 	  c2total += ((c2 << C2_SHIFT) + ((1<<C2_SHIFT)>>1)) * count;
    512 	}
    513       }
    514     }
    515 
    516   cinfo->colormap[0][icolor] = (JSAMPLE) ((c0total + (total>>1)) / total);
    517   cinfo->colormap[1][icolor] = (JSAMPLE) ((c1total + (total>>1)) / total);
    518   cinfo->colormap[2][icolor] = (JSAMPLE) ((c2total + (total>>1)) / total);
    519 }
    520 
    521 
    522 LOCAL(void)
    523 select_colors (j_decompress_ptr cinfo, int desired_colors)
    524 /* Master routine for color selection */
    525 {
    526   boxptr boxlist;
    527   int numboxes;
    528   int i;
    529 
    530   /* Allocate workspace for box list */
    531   boxlist = (boxptr) (*cinfo->mem->alloc_small)
    532     ((j_common_ptr) cinfo, JPOOL_IMAGE, desired_colors * SIZEOF(box));
    533   /* Initialize one box containing whole space */
    534   numboxes = 1;
    535   boxlist[0].c0min = 0;
    536   boxlist[0].c0max = MAXJSAMPLE >> C0_SHIFT;
    537   boxlist[0].c1min = 0;
    538   boxlist[0].c1max = MAXJSAMPLE >> C1_SHIFT;
    539   boxlist[0].c2min = 0;
    540   boxlist[0].c2max = MAXJSAMPLE >> C2_SHIFT;
    541   /* Shrink it to actually-used volume and set its statistics */
    542   update_box(cinfo, & boxlist[0]);
    543   /* Perform median-cut to produce final box list */
    544   numboxes = median_cut(cinfo, boxlist, numboxes, desired_colors);
    545   /* Compute the representative color for each box, fill colormap */
    546   for (i = 0; i < numboxes; i++)
    547     compute_color(cinfo, & boxlist[i], i);
    548   cinfo->actual_number_of_colors = numboxes;
    549   TRACEMS1(cinfo, 1, JTRC_QUANT_SELECTED, numboxes);
    550 }
    551 
    552 
    553 /*
    554  * These routines are concerned with the time-critical task of mapping input
    555  * colors to the nearest color in the selected colormap.
    556  *
    557  * We re-use the histogram space as an "inverse color map", essentially a
    558  * cache for the results of nearest-color searches.  All colors within a
    559  * histogram cell will be mapped to the same colormap entry, namely the one
    560  * closest to the cell's center.  This may not be quite the closest entry to
    561  * the actual input color, but it's almost as good.  A zero in the cache
    562  * indicates we haven't found the nearest color for that cell yet; the array
    563  * is cleared to zeroes before starting the mapping pass.  When we find the
    564  * nearest color for a cell, its colormap index plus one is recorded in the
    565  * cache for future use.  The pass2 scanning routines call fill_inverse_cmap
    566  * when they need to use an unfilled entry in the cache.
    567  *
    568  * Our method of efficiently finding nearest colors is based on the "locally
    569  * sorted search" idea described by Heckbert and on the incremental distance
    570  * calculation described by Spencer W. Thomas in chapter III.1 of Graphics
    571  * Gems II (James Arvo, ed.  Academic Press, 1991).  Thomas points out that
    572  * the distances from a given colormap entry to each cell of the histogram can
    573  * be computed quickly using an incremental method: the differences between
    574  * distances to adjacent cells themselves differ by a constant.  This allows a
    575  * fairly fast implementation of the "brute force" approach of computing the
    576  * distance from every colormap entry to every histogram cell.  Unfortunately,
    577  * it needs a work array to hold the best-distance-so-far for each histogram
    578  * cell (because the inner loop has to be over cells, not colormap entries).
    579  * The work array elements have to be INT32s, so the work array would need
    580  * 256Kb at our recommended precision.  This is not feasible in DOS machines.
    581  *
    582  * To get around these problems, we apply Thomas' method to compute the
    583  * nearest colors for only the cells within a small subbox of the histogram.
    584  * The work array need be only as big as the subbox, so the memory usage
    585  * problem is solved.  Furthermore, we need not fill subboxes that are never
    586  * referenced in pass2; many images use only part of the color gamut, so a
    587  * fair amount of work is saved.  An additional advantage of this
    588  * approach is that we can apply Heckbert's locality criterion to quickly
    589  * eliminate colormap entries that are far away from the subbox; typically
    590  * three-fourths of the colormap entries are rejected by Heckbert's criterion,
    591  * and we need not compute their distances to individual cells in the subbox.
    592  * The speed of this approach is heavily influenced by the subbox size: too
    593  * small means too much overhead, too big loses because Heckbert's criterion
    594  * can't eliminate as many colormap entries.  Empirically the best subbox
    595  * size seems to be about 1/512th of the histogram (1/8th in each direction).
    596  *
    597  * Thomas' article also describes a refined method which is asymptotically
    598  * faster than the brute-force method, but it is also far more complex and
    599  * cannot efficiently be applied to small subboxes.  It is therefore not
    600  * useful for programs intended to be portable to DOS machines.  On machines
    601  * with plenty of memory, filling the whole histogram in one shot with Thomas'
    602  * refined method might be faster than the present code --- but then again,
    603  * it might not be any faster, and it's certainly more complicated.
    604  */
    605 
    606 
    607 /* log2(histogram cells in update box) for each axis; this can be adjusted */
    608 #define BOX_C0_LOG  (HIST_C0_BITS-3)
    609 #define BOX_C1_LOG  (HIST_C1_BITS-3)
    610 #define BOX_C2_LOG  (HIST_C2_BITS-3)
    611 
    612 #define BOX_C0_ELEMS  (1<<BOX_C0_LOG) /* # of hist cells in update box */
    613 #define BOX_C1_ELEMS  (1<<BOX_C1_LOG)
    614 #define BOX_C2_ELEMS  (1<<BOX_C2_LOG)
    615 
    616 #define BOX_C0_SHIFT  (C0_SHIFT + BOX_C0_LOG)
    617 #define BOX_C1_SHIFT  (C1_SHIFT + BOX_C1_LOG)
    618 #define BOX_C2_SHIFT  (C2_SHIFT + BOX_C2_LOG)
    619 
    620 
    621 /*
    622  * The next three routines implement inverse colormap filling.  They could
    623  * all be folded into one big routine, but splitting them up this way saves
    624  * some stack space (the mindist[] and bestdist[] arrays need not coexist)
    625  * and may allow some compilers to produce better code by registerizing more
    626  * inner-loop variables.
    627  */
    628 
    629 LOCAL(int)
    630 find_nearby_colors (j_decompress_ptr cinfo, int minc0, int minc1, int minc2,
    631 		    JSAMPLE colorlist[])
    632 /* Locate the colormap entries close enough to an update box to be candidates
    633  * for the nearest entry to some cell(s) in the update box.  The update box
    634  * is specified by the center coordinates of its first cell.  The number of
    635  * candidate colormap entries is returned, and their colormap indexes are
    636  * placed in colorlist[].
    637  * This routine uses Heckbert's "locally sorted search" criterion to select
    638  * the colors that need further consideration.
    639  */
    640 {
    641   int numcolors = cinfo->actual_number_of_colors;
    642   int maxc0, maxc1, maxc2;
    643   int centerc0, centerc1, centerc2;
    644   int i, x, ncolors;
    645   INT32 minmaxdist, min_dist, max_dist, tdist;
    646   INT32 mindist[MAXNUMCOLORS];	/* min distance to colormap entry i */
    647 
    648   /* Compute true coordinates of update box's upper corner and center.
    649    * Actually we compute the coordinates of the center of the upper-corner
    650    * histogram cell, which are the upper bounds of the volume we care about.
    651    * Note that since ">>" rounds down, the "center" values may be closer to
    652    * min than to max; hence comparisons to them must be "<=", not "<".
    653    */
    654   maxc0 = minc0 + ((1 << BOX_C0_SHIFT) - (1 << C0_SHIFT));
    655   centerc0 = (minc0 + maxc0) >> 1;
    656   maxc1 = minc1 + ((1 << BOX_C1_SHIFT) - (1 << C1_SHIFT));
    657   centerc1 = (minc1 + maxc1) >> 1;
    658   maxc2 = minc2 + ((1 << BOX_C2_SHIFT) - (1 << C2_SHIFT));
    659   centerc2 = (minc2 + maxc2) >> 1;
    660 
    661   /* For each color in colormap, find:
    662    *  1. its minimum squared-distance to any point in the update box
    663    *     (zero if color is within update box);
    664    *  2. its maximum squared-distance to any point in the update box.
    665    * Both of these can be found by considering only the corners of the box.
    666    * We save the minimum distance for each color in mindist[];
    667    * only the smallest maximum distance is of interest.
    668    */
    669   minmaxdist = 0x7FFFFFFFL;
    670 
    671   for (i = 0; i < numcolors; i++) {
    672     /* We compute the squared-c0-distance term, then add in the other two. */
    673     x = GETJSAMPLE(cinfo->colormap[0][i]);
    674     if (x < minc0) {
    675       tdist = (x - minc0) * C0_SCALE;
    676       min_dist = tdist*tdist;
    677       tdist = (x - maxc0) * C0_SCALE;
    678       max_dist = tdist*tdist;
    679     } else if (x > maxc0) {
    680       tdist = (x - maxc0) * C0_SCALE;
    681       min_dist = tdist*tdist;
    682       tdist = (x - minc0) * C0_SCALE;
    683       max_dist = tdist*tdist;
    684     } else {
    685       /* within cell range so no contribution to min_dist */
    686       min_dist = 0;
    687       if (x <= centerc0) {
    688 	tdist = (x - maxc0) * C0_SCALE;
    689 	max_dist = tdist*tdist;
    690       } else {
    691 	tdist = (x - minc0) * C0_SCALE;
    692 	max_dist = tdist*tdist;
    693       }
    694     }
    695 
    696     x = GETJSAMPLE(cinfo->colormap[1][i]);
    697     if (x < minc1) {
    698       tdist = (x - minc1) * C1_SCALE;
    699       min_dist += tdist*tdist;
    700       tdist = (x - maxc1) * C1_SCALE;
    701       max_dist += tdist*tdist;
    702     } else if (x > maxc1) {
    703       tdist = (x - maxc1) * C1_SCALE;
    704       min_dist += tdist*tdist;
    705       tdist = (x - minc1) * C1_SCALE;
    706       max_dist += tdist*tdist;
    707     } else {
    708       /* within cell range so no contribution to min_dist */
    709       if (x <= centerc1) {
    710 	tdist = (x - maxc1) * C1_SCALE;
    711 	max_dist += tdist*tdist;
    712       } else {
    713 	tdist = (x - minc1) * C1_SCALE;
    714 	max_dist += tdist*tdist;
    715       }
    716     }
    717 
    718     x = GETJSAMPLE(cinfo->colormap[2][i]);
    719     if (x < minc2) {
    720       tdist = (x - minc2) * C2_SCALE;
    721       min_dist += tdist*tdist;
    722       tdist = (x - maxc2) * C2_SCALE;
    723       max_dist += tdist*tdist;
    724     } else if (x > maxc2) {
    725       tdist = (x - maxc2) * C2_SCALE;
    726       min_dist += tdist*tdist;
    727       tdist = (x - minc2) * C2_SCALE;
    728       max_dist += tdist*tdist;
    729     } else {
    730       /* within cell range so no contribution to min_dist */
    731       if (x <= centerc2) {
    732 	tdist = (x - maxc2) * C2_SCALE;
    733 	max_dist += tdist*tdist;
    734       } else {
    735 	tdist = (x - minc2) * C2_SCALE;
    736 	max_dist += tdist*tdist;
    737       }
    738     }
    739 
    740     mindist[i] = min_dist;	/* save away the results */
    741     if (max_dist < minmaxdist)
    742       minmaxdist = max_dist;
    743   }
    744 
    745   /* Now we know that no cell in the update box is more than minmaxdist
    746    * away from some colormap entry.  Therefore, only colors that are
    747    * within minmaxdist of some part of the box need be considered.
    748    */
    749   ncolors = 0;
    750   for (i = 0; i < numcolors; i++) {
    751     if (mindist[i] <= minmaxdist)
    752       colorlist[ncolors++] = (JSAMPLE) i;
    753   }
    754   return ncolors;
    755 }
    756 
    757 
    758 LOCAL(void)
    759 find_best_colors (j_decompress_ptr cinfo, int minc0, int minc1, int minc2,
    760 		  int numcolors, JSAMPLE colorlist[], JSAMPLE bestcolor[])
    761 /* Find the closest colormap entry for each cell in the update box,
    762  * given the list of candidate colors prepared by find_nearby_colors.
    763  * Return the indexes of the closest entries in the bestcolor[] array.
    764  * This routine uses Thomas' incremental distance calculation method to
    765  * find the distance from a colormap entry to successive cells in the box.
    766  */
    767 {
    768   int ic0, ic1, ic2;
    769   int i, icolor;
    770   register INT32 * bptr;	/* pointer into bestdist[] array */
    771   JSAMPLE * cptr;		/* pointer into bestcolor[] array */
    772   INT32 dist0, dist1;		/* initial distance values */
    773   register INT32 dist2;		/* current distance in inner loop */
    774   INT32 xx0, xx1;		/* distance increments */
    775   register INT32 xx2;
    776   INT32 inc0, inc1, inc2;	/* initial values for increments */
    777   /* This array holds the distance to the nearest-so-far color for each cell */
    778   INT32 bestdist[BOX_C0_ELEMS * BOX_C1_ELEMS * BOX_C2_ELEMS];
    779 
    780   /* Initialize best-distance for each cell of the update box */
    781   bptr = bestdist;
    782   for (i = BOX_C0_ELEMS*BOX_C1_ELEMS*BOX_C2_ELEMS-1; i >= 0; i--)
    783     *bptr++ = 0x7FFFFFFFL;
    784 
    785   /* For each color selected by find_nearby_colors,
    786    * compute its distance to the center of each cell in the box.
    787    * If that's less than best-so-far, update best distance and color number.
    788    */
    789 
    790   /* Nominal steps between cell centers ("x" in Thomas article) */
    791 #define STEP_C0  ((1 << C0_SHIFT) * C0_SCALE)
    792 #define STEP_C1  ((1 << C1_SHIFT) * C1_SCALE)
    793 #define STEP_C2  ((1 << C2_SHIFT) * C2_SCALE)
    794 
    795   for (i = 0; i < numcolors; i++) {
    796     icolor = GETJSAMPLE(colorlist[i]);
    797     /* Compute (square of) distance from minc0/c1/c2 to this color */
    798     inc0 = (minc0 - GETJSAMPLE(cinfo->colormap[0][icolor])) * C0_SCALE;
    799     dist0 = inc0*inc0;
    800     inc1 = (minc1 - GETJSAMPLE(cinfo->colormap[1][icolor])) * C1_SCALE;
    801     dist0 += inc1*inc1;
    802     inc2 = (minc2 - GETJSAMPLE(cinfo->colormap[2][icolor])) * C2_SCALE;
    803     dist0 += inc2*inc2;
    804     /* Form the initial difference increments */
    805     inc0 = inc0 * (2 * STEP_C0) + STEP_C0 * STEP_C0;
    806     inc1 = inc1 * (2 * STEP_C1) + STEP_C1 * STEP_C1;
    807     inc2 = inc2 * (2 * STEP_C2) + STEP_C2 * STEP_C2;
    808     /* Now loop over all cells in box, updating distance per Thomas method */
    809     bptr = bestdist;
    810     cptr = bestcolor;
    811     xx0 = inc0;
    812     for (ic0 = BOX_C0_ELEMS-1; ic0 >= 0; ic0--) {
    813       dist1 = dist0;
    814       xx1 = inc1;
    815       for (ic1 = BOX_C1_ELEMS-1; ic1 >= 0; ic1--) {
    816 	dist2 = dist1;
    817 	xx2 = inc2;
    818 	for (ic2 = BOX_C2_ELEMS-1; ic2 >= 0; ic2--) {
    819 	  if (dist2 < *bptr) {
    820 	    *bptr = dist2;
    821 	    *cptr = (JSAMPLE) icolor;
    822 	  }
    823 	  dist2 += xx2;
    824 	  xx2 += 2 * STEP_C2 * STEP_C2;
    825 	  bptr++;
    826 	  cptr++;
    827 	}
    828 	dist1 += xx1;
    829 	xx1 += 2 * STEP_C1 * STEP_C1;
    830       }
    831       dist0 += xx0;
    832       xx0 += 2 * STEP_C0 * STEP_C0;
    833     }
    834   }
    835 }
    836 
    837 
    838 LOCAL(void)
    839 fill_inverse_cmap (j_decompress_ptr cinfo, int c0, int c1, int c2)
    840 /* Fill the inverse-colormap entries in the update box that contains */
    841 /* histogram cell c0/c1/c2.  (Only that one cell MUST be filled, but */
    842 /* we can fill as many others as we wish.) */
    843 {
    844   my_cquantize_ptr cquantize = (my_cquantize_ptr) cinfo->cquantize;
    845   hist3d histogram = cquantize->histogram;
    846   int minc0, minc1, minc2;	/* lower left corner of update box */
    847   int ic0, ic1, ic2;
    848   register JSAMPLE * cptr;	/* pointer into bestcolor[] array */
    849   register histptr cachep;	/* pointer into main cache array */
    850   /* This array lists the candidate colormap indexes. */
    851   JSAMPLE colorlist[MAXNUMCOLORS];
    852   int numcolors;		/* number of candidate colors */
    853   /* This array holds the actually closest colormap index for each cell. */
    854   JSAMPLE bestcolor[BOX_C0_ELEMS * BOX_C1_ELEMS * BOX_C2_ELEMS];
    855 
    856   /* Convert cell coordinates to update box ID */
    857   c0 >>= BOX_C0_LOG;
    858   c1 >>= BOX_C1_LOG;
    859   c2 >>= BOX_C2_LOG;
    860 
    861   /* Compute true coordinates of update box's origin corner.
    862    * Actually we compute the coordinates of the center of the corner
    863    * histogram cell, which are the lower bounds of the volume we care about.
    864    */
    865   minc0 = (c0 << BOX_C0_SHIFT) + ((1 << C0_SHIFT) >> 1);
    866   minc1 = (c1 << BOX_C1_SHIFT) + ((1 << C1_SHIFT) >> 1);
    867   minc2 = (c2 << BOX_C2_SHIFT) + ((1 << C2_SHIFT) >> 1);
    868 
    869   /* Determine which colormap entries are close enough to be candidates
    870    * for the nearest entry to some cell in the update box.
    871    */
    872   numcolors = find_nearby_colors(cinfo, minc0, minc1, minc2, colorlist);
    873 
    874   /* Determine the actually nearest colors. */
    875   find_best_colors(cinfo, minc0, minc1, minc2, numcolors, colorlist,
    876 		   bestcolor);
    877 
    878   /* Save the best color numbers (plus 1) in the main cache array */
    879   c0 <<= BOX_C0_LOG;		/* convert ID back to base cell indexes */
    880   c1 <<= BOX_C1_LOG;
    881   c2 <<= BOX_C2_LOG;
    882   cptr = bestcolor;
    883   for (ic0 = 0; ic0 < BOX_C0_ELEMS; ic0++) {
    884     for (ic1 = 0; ic1 < BOX_C1_ELEMS; ic1++) {
    885       cachep = & histogram[c0+ic0][c1+ic1][c2];
    886       for (ic2 = 0; ic2 < BOX_C2_ELEMS; ic2++) {
    887 	*cachep++ = (histcell) (GETJSAMPLE(*cptr++) + 1);
    888       }
    889     }
    890   }
    891 }
    892 
    893 
    894 /*
    895  * Map some rows of pixels to the output colormapped representation.
    896  */
    897 
    898 METHODDEF(void)
    899 pass2_no_dither (j_decompress_ptr cinfo,
    900 		 JSAMPARRAY input_buf, JSAMPARRAY output_buf, int num_rows)
    901 /* This version performs no dithering */
    902 {
    903   my_cquantize_ptr cquantize = (my_cquantize_ptr) cinfo->cquantize;
    904   hist3d histogram = cquantize->histogram;
    905   register JSAMPROW inptr, outptr;
    906   register histptr cachep;
    907   register int c0, c1, c2;
    908   int row;
    909   JDIMENSION col;
    910   JDIMENSION width = cinfo->output_width;
    911 
    912   for (row = 0; row < num_rows; row++) {
    913     inptr = input_buf[row];
    914     outptr = output_buf[row];
    915     for (col = width; col > 0; col--) {
    916       /* get pixel value and index into the cache */
    917       c0 = GETJSAMPLE(*inptr++) >> C0_SHIFT;
    918       c1 = GETJSAMPLE(*inptr++) >> C1_SHIFT;
    919       c2 = GETJSAMPLE(*inptr++) >> C2_SHIFT;
    920       cachep = & histogram[c0][c1][c2];
    921       /* If we have not seen this color before, find nearest colormap entry */
    922       /* and update the cache */
    923       if (*cachep == 0)
    924 	fill_inverse_cmap(cinfo, c0,c1,c2);
    925       /* Now emit the colormap index for this cell */
    926       *outptr++ = (JSAMPLE) (*cachep - 1);
    927     }
    928   }
    929 }
    930 
    931 
    932 METHODDEF(void)
    933 pass2_fs_dither (j_decompress_ptr cinfo,
    934 		 JSAMPARRAY input_buf, JSAMPARRAY output_buf, int num_rows)
    935 /* This version performs Floyd-Steinberg dithering */
    936 {
    937   my_cquantize_ptr cquantize = (my_cquantize_ptr) cinfo->cquantize;
    938   hist3d histogram = cquantize->histogram;
    939   register LOCFSERROR cur0, cur1, cur2;	/* current error or pixel value */
    940   LOCFSERROR belowerr0, belowerr1, belowerr2; /* error for pixel below cur */
    941   LOCFSERROR bpreverr0, bpreverr1, bpreverr2; /* error for below/prev col */
    942   register FSERRPTR errorptr;	/* => fserrors[] at column before current */
    943   JSAMPROW inptr;		/* => current input pixel */
    944   JSAMPROW outptr;		/* => current output pixel */
    945   histptr cachep;
    946   int dir;			/* +1 or -1 depending on direction */
    947   int dir3;			/* 3*dir, for advancing inptr & errorptr */
    948   int row;
    949   JDIMENSION col;
    950   JDIMENSION width = cinfo->output_width;
    951   JSAMPLE *range_limit = cinfo->sample_range_limit;
    952   int *error_limit = cquantize->error_limiter;
    953   JSAMPROW colormap0 = cinfo->colormap[0];
    954   JSAMPROW colormap1 = cinfo->colormap[1];
    955   JSAMPROW colormap2 = cinfo->colormap[2];
    956   SHIFT_TEMPS
    957 
    958   for (row = 0; row < num_rows; row++) {
    959     inptr = input_buf[row];
    960     outptr = output_buf[row];
    961     if (cquantize->on_odd_row) {
    962       /* work right to left in this row */
    963       inptr += (width-1) * 3;	/* so point to rightmost pixel */
    964       outptr += width-1;
    965       dir = -1;
    966       dir3 = -3;
    967       errorptr = cquantize->fserrors + (width+1)*3; /* => entry after last column */
    968       cquantize->on_odd_row = FALSE; /* flip for next time */
    969     } else {
    970       /* work left to right in this row */
    971       dir = 1;
    972       dir3 = 3;
    973       errorptr = cquantize->fserrors; /* => entry before first real column */
    974       cquantize->on_odd_row = TRUE; /* flip for next time */
    975     }
    976     /* Preset error values: no error propagated to first pixel from left */
    977     cur0 = cur1 = cur2 = 0;
    978     /* and no error propagated to row below yet */
    979     belowerr0 = belowerr1 = belowerr2 = 0;
    980     bpreverr0 = bpreverr1 = bpreverr2 = 0;
    981 
    982     for (col = width; col > 0; col--) {
    983       /* curN holds the error propagated from the previous pixel on the
    984        * current line.  Add the error propagated from the previous line
    985        * to form the complete error correction term for this pixel, and
    986        * round the error term (which is expressed * 16) to an integer.
    987        * RIGHT_SHIFT rounds towards minus infinity, so adding 8 is correct
    988        * for either sign of the error value.
    989        * Note: errorptr points to *previous* column's array entry.
    990        */
    991       cur0 = RIGHT_SHIFT(cur0 + errorptr[dir3+0] + 8, 4);
    992       cur1 = RIGHT_SHIFT(cur1 + errorptr[dir3+1] + 8, 4);
    993       cur2 = RIGHT_SHIFT(cur2 + errorptr[dir3+2] + 8, 4);
    994       /* Limit the error using transfer function set by init_error_limit.
    995        * See comments with init_error_limit for rationale.
    996        */
    997       cur0 = error_limit[cur0];
    998       cur1 = error_limit[cur1];
    999       cur2 = error_limit[cur2];
   1000       /* Form pixel value + error, and range-limit to 0..MAXJSAMPLE.
   1001        * The maximum error is +- MAXJSAMPLE (or less with error limiting);
   1002        * this sets the required size of the range_limit array.
   1003        */
   1004       cur0 += GETJSAMPLE(inptr[0]);
   1005       cur1 += GETJSAMPLE(inptr[1]);
   1006       cur2 += GETJSAMPLE(inptr[2]);
   1007       cur0 = GETJSAMPLE(range_limit[cur0]);
   1008       cur1 = GETJSAMPLE(range_limit[cur1]);
   1009       cur2 = GETJSAMPLE(range_limit[cur2]);
   1010       /* Index into the cache with adjusted pixel value */
   1011       cachep = & histogram[cur0>>C0_SHIFT][cur1>>C1_SHIFT][cur2>>C2_SHIFT];
   1012       /* If we have not seen this color before, find nearest colormap */
   1013       /* entry and update the cache */
   1014       if (*cachep == 0)
   1015 	fill_inverse_cmap(cinfo, cur0>>C0_SHIFT,cur1>>C1_SHIFT,cur2>>C2_SHIFT);
   1016       /* Now emit the colormap index for this cell */
   1017       { register int pixcode = *cachep - 1;
   1018 	*outptr = (JSAMPLE) pixcode;
   1019 	/* Compute representation error for this pixel */
   1020 	cur0 -= GETJSAMPLE(colormap0[pixcode]);
   1021 	cur1 -= GETJSAMPLE(colormap1[pixcode]);
   1022 	cur2 -= GETJSAMPLE(colormap2[pixcode]);
   1023       }
   1024       /* Compute error fractions to be propagated to adjacent pixels.
   1025        * Add these into the running sums, and simultaneously shift the
   1026        * next-line error sums left by 1 column.
   1027        */
   1028       { register LOCFSERROR bnexterr, delta;
   1029 
   1030 	bnexterr = cur0;	/* Process component 0 */
   1031 	delta = cur0 * 2;
   1032 	cur0 += delta;		/* form error * 3 */
   1033 	errorptr[0] = (FSERROR) (bpreverr0 + cur0);
   1034 	cur0 += delta;		/* form error * 5 */
   1035 	bpreverr0 = belowerr0 + cur0;
   1036 	belowerr0 = bnexterr;
   1037 	cur0 += delta;		/* form error * 7 */
   1038 	bnexterr = cur1;	/* Process component 1 */
   1039 	delta = cur1 * 2;
   1040 	cur1 += delta;		/* form error * 3 */
   1041 	errorptr[1] = (FSERROR) (bpreverr1 + cur1);
   1042 	cur1 += delta;		/* form error * 5 */
   1043 	bpreverr1 = belowerr1 + cur1;
   1044 	belowerr1 = bnexterr;
   1045 	cur1 += delta;		/* form error * 7 */
   1046 	bnexterr = cur2;	/* Process component 2 */
   1047 	delta = cur2 * 2;
   1048 	cur2 += delta;		/* form error * 3 */
   1049 	errorptr[2] = (FSERROR) (bpreverr2 + cur2);
   1050 	cur2 += delta;		/* form error * 5 */
   1051 	bpreverr2 = belowerr2 + cur2;
   1052 	belowerr2 = bnexterr;
   1053 	cur2 += delta;		/* form error * 7 */
   1054       }
   1055       /* At this point curN contains the 7/16 error value to be propagated
   1056        * to the next pixel on the current line, and all the errors for the
   1057        * next line have been shifted over.  We are therefore ready to move on.
   1058        */
   1059       inptr += dir3;		/* Advance pixel pointers to next column */
   1060       outptr += dir;
   1061       errorptr += dir3;		/* advance errorptr to current column */
   1062     }
   1063     /* Post-loop cleanup: we must unload the final error values into the
   1064      * final fserrors[] entry.  Note we need not unload belowerrN because
   1065      * it is for the dummy column before or after the actual array.
   1066      */
   1067     errorptr[0] = (FSERROR) bpreverr0; /* unload prev errs into array */
   1068     errorptr[1] = (FSERROR) bpreverr1;
   1069     errorptr[2] = (FSERROR) bpreverr2;
   1070   }
   1071 }
   1072 
   1073 
   1074 /*
   1075  * Initialize the error-limiting transfer function (lookup table).
   1076  * The raw F-S error computation can potentially compute error values of up to
   1077  * +- MAXJSAMPLE.  But we want the maximum correction applied to a pixel to be
   1078  * much less, otherwise obviously wrong pixels will be created.  (Typical
   1079  * effects include weird fringes at color-area boundaries, isolated bright
   1080  * pixels in a dark area, etc.)  The standard advice for avoiding this problem
   1081  * is to ensure that the "corners" of the color cube are allocated as output
   1082  * colors; then repeated errors in the same direction cannot cause cascading
   1083  * error buildup.  However, that only prevents the error from getting
   1084  * completely out of hand; Aaron Giles reports that error limiting improves
   1085  * the results even with corner colors allocated.
   1086  * A simple clamping of the error values to about +- MAXJSAMPLE/8 works pretty
   1087  * well, but the smoother transfer function used below is even better.  Thanks
   1088  * to Aaron Giles for this idea.
   1089  */
   1090 
   1091 LOCAL(void)
   1092 init_error_limit (j_decompress_ptr cinfo)
   1093 /* Allocate and fill in the error_limiter table */
   1094 {
   1095   my_cquantize_ptr cquantize = (my_cquantize_ptr) cinfo->cquantize;
   1096   int * table;
   1097   int in, out;
   1098 
   1099   table = (int *) (*cinfo->mem->alloc_small)
   1100     ((j_common_ptr) cinfo, JPOOL_IMAGE, (MAXJSAMPLE*2+1) * SIZEOF(int));
   1101   table += MAXJSAMPLE;		/* so can index -MAXJSAMPLE .. +MAXJSAMPLE */
   1102   cquantize->error_limiter = table;
   1103 
   1104 #define STEPSIZE ((MAXJSAMPLE+1)/16)
   1105   /* Map errors 1:1 up to +- MAXJSAMPLE/16 */
   1106   out = 0;
   1107   for (in = 0; in < STEPSIZE; in++, out++) {
   1108     table[in] = out; table[-in] = -out;
   1109   }
   1110   /* Map errors 1:2 up to +- 3*MAXJSAMPLE/16 */
   1111   for (; in < STEPSIZE*3; in++, out += (in&1) ? 0 : 1) {
   1112     table[in] = out; table[-in] = -out;
   1113   }
   1114   /* Clamp the rest to final out value (which is (MAXJSAMPLE+1)/8) */
   1115   for (; in <= MAXJSAMPLE; in++) {
   1116     table[in] = out; table[-in] = -out;
   1117   }
   1118 #undef STEPSIZE
   1119 }
   1120 
   1121 
   1122 /*
   1123  * Finish up at the end of each pass.
   1124  */
   1125 
   1126 METHODDEF(void)
   1127 finish_pass1 (j_decompress_ptr cinfo)
   1128 {
   1129   my_cquantize_ptr cquantize = (my_cquantize_ptr) cinfo->cquantize;
   1130 
   1131   /* Select the representative colors and fill in cinfo->colormap */
   1132   cinfo->colormap = cquantize->sv_colormap;
   1133   select_colors(cinfo, cquantize->desired);
   1134   /* Force next pass to zero the color index table */
   1135   cquantize->needs_zeroed = TRUE;
   1136 }
   1137 
   1138 
   1139 METHODDEF(void)
   1140 finish_pass2 (j_decompress_ptr cinfo)
   1141 {
   1142   /* no work */
   1143 }
   1144 
   1145 
   1146 /*
   1147  * Initialize for each processing pass.
   1148  */
   1149 
   1150 METHODDEF(void)
   1151 start_pass_2_quant (j_decompress_ptr cinfo, boolean is_pre_scan)
   1152 {
   1153   my_cquantize_ptr cquantize = (my_cquantize_ptr) cinfo->cquantize;
   1154   hist3d histogram = cquantize->histogram;
   1155   int i;
   1156 
   1157   /* Only F-S dithering or no dithering is supported. */
   1158   /* If user asks for ordered dither, give him F-S. */
   1159   if (cinfo->dither_mode != JDITHER_NONE)
   1160     cinfo->dither_mode = JDITHER_FS;
   1161 
   1162   if (is_pre_scan) {
   1163     /* Set up method pointers */
   1164     cquantize->pub.color_quantize = prescan_quantize;
   1165     cquantize->pub.finish_pass = finish_pass1;
   1166     cquantize->needs_zeroed = TRUE; /* Always zero histogram */
   1167   } else {
   1168     /* Set up method pointers */
   1169     if (cinfo->dither_mode == JDITHER_FS)
   1170       cquantize->pub.color_quantize = pass2_fs_dither;
   1171     else
   1172       cquantize->pub.color_quantize = pass2_no_dither;
   1173     cquantize->pub.finish_pass = finish_pass2;
   1174 
   1175     /* Make sure color count is acceptable */
   1176     i = cinfo->actual_number_of_colors;
   1177     if (i < 1)
   1178       ERREXIT1(cinfo, JERR_QUANT_FEW_COLORS, 1);
   1179     if (i > MAXNUMCOLORS)
   1180       ERREXIT1(cinfo, JERR_QUANT_MANY_COLORS, MAXNUMCOLORS);
   1181 
   1182     if (cinfo->dither_mode == JDITHER_FS) {
   1183       size_t arraysize = (size_t) ((cinfo->output_width + 2) *
   1184 				   (3 * SIZEOF(FSERROR)));
   1185       /* Allocate Floyd-Steinberg workspace if we didn't already. */
   1186       if (cquantize->fserrors == NULL)
   1187 	cquantize->fserrors = (FSERRPTR) (*cinfo->mem->alloc_large)
   1188 	  ((j_common_ptr) cinfo, JPOOL_IMAGE, arraysize);
   1189       /* Initialize the propagated errors to zero. */
   1190       jzero_far((void FAR *) cquantize->fserrors, arraysize);
   1191       /* Make the error-limit table if we didn't already. */
   1192       if (cquantize->error_limiter == NULL)
   1193 	init_error_limit(cinfo);
   1194       cquantize->on_odd_row = FALSE;
   1195     }
   1196 
   1197   }
   1198   /* Zero the histogram or inverse color map, if necessary */
   1199   if (cquantize->needs_zeroed) {
   1200     for (i = 0; i < HIST_C0_ELEMS; i++) {
   1201       jzero_far((void FAR *) histogram[i],
   1202 		HIST_C1_ELEMS*HIST_C2_ELEMS * SIZEOF(histcell));
   1203     }
   1204     cquantize->needs_zeroed = FALSE;
   1205   }
   1206 }
   1207 
   1208 
   1209 /*
   1210  * Switch to a new external colormap between output passes.
   1211  */
   1212 
   1213 METHODDEF(void)
   1214 new_color_map_2_quant (j_decompress_ptr cinfo)
   1215 {
   1216   my_cquantize_ptr cquantize = (my_cquantize_ptr) cinfo->cquantize;
   1217 
   1218   /* Reset the inverse color map */
   1219   cquantize->needs_zeroed = TRUE;
   1220 }
   1221 
   1222 
   1223 /*
   1224  * Module initialization routine for 2-pass color quantization.
   1225  */
   1226 
   1227 GLOBAL(void)
   1228 jinit_2pass_quantizer (j_decompress_ptr cinfo)
   1229 {
   1230   my_cquantize_ptr cquantize;
   1231   int i;
   1232 
   1233   cquantize = (my_cquantize_ptr)
   1234     (*cinfo->mem->alloc_small) ((j_common_ptr) cinfo, JPOOL_IMAGE,
   1235 				SIZEOF(my_cquantizer));
   1236   cinfo->cquantize = (struct jpeg_color_quantizer *) cquantize;
   1237   cquantize->pub.start_pass = start_pass_2_quant;
   1238   cquantize->pub.new_color_map = new_color_map_2_quant;
   1239   cquantize->fserrors = NULL;	/* flag optional arrays not allocated */
   1240   cquantize->error_limiter = NULL;
   1241 
   1242   /* Make sure jdmaster didn't give me a case I can't handle */
   1243   if (cinfo->out_color_components != 3)
   1244     ERREXIT(cinfo, JERR_NOTIMPL);
   1245 
   1246   /* Allocate the histogram/inverse colormap storage */
   1247   cquantize->histogram = (hist3d) (*cinfo->mem->alloc_small)
   1248     ((j_common_ptr) cinfo, JPOOL_IMAGE, HIST_C0_ELEMS * SIZEOF(hist2d));
   1249   for (i = 0; i < HIST_C0_ELEMS; i++) {
   1250     cquantize->histogram[i] = (hist2d) (*cinfo->mem->alloc_large)
   1251       ((j_common_ptr) cinfo, JPOOL_IMAGE,
   1252        HIST_C1_ELEMS*HIST_C2_ELEMS * SIZEOF(histcell));
   1253   }
   1254   cquantize->needs_zeroed = TRUE; /* histogram is garbage now */
   1255 
   1256   /* Allocate storage for the completed colormap, if required.
   1257    * We do this now since it is FAR storage and may affect
   1258    * the memory manager's space calculations.
   1259    */
   1260   if (cinfo->enable_2pass_quant) {
   1261     /* Make sure color count is acceptable */
   1262     int desired = cinfo->desired_number_of_colors;
   1263     /* Lower bound on # of colors ... somewhat arbitrary as long as > 0 */
   1264     if (desired < 8)
   1265       ERREXIT1(cinfo, JERR_QUANT_FEW_COLORS, 8);
   1266     /* Make sure colormap indexes can be represented by JSAMPLEs */
   1267     if (desired > MAXNUMCOLORS)
   1268       ERREXIT1(cinfo, JERR_QUANT_MANY_COLORS, MAXNUMCOLORS);
   1269     cquantize->sv_colormap = (*cinfo->mem->alloc_sarray)
   1270       ((j_common_ptr) cinfo,JPOOL_IMAGE, (JDIMENSION) desired, (JDIMENSION) 3);
   1271     cquantize->desired = desired;
   1272   } else
   1273     cquantize->sv_colormap = NULL;
   1274 
   1275   /* Only F-S dithering or no dithering is supported. */
   1276   /* If user asks for ordered dither, give him F-S. */
   1277   if (cinfo->dither_mode != JDITHER_NONE)
   1278     cinfo->dither_mode = JDITHER_FS;
   1279 
   1280   /* Allocate Floyd-Steinberg workspace if necessary.
   1281    * This isn't really needed until pass 2, but again it is FAR storage.
   1282    * Although we will cope with a later change in dither_mode,
   1283    * we do not promise to honor max_memory_to_use if dither_mode changes.
   1284    */
   1285   if (cinfo->dither_mode == JDITHER_FS) {
   1286     cquantize->fserrors = (FSERRPTR) (*cinfo->mem->alloc_large)
   1287       ((j_common_ptr) cinfo, JPOOL_IMAGE,
   1288        (size_t) ((cinfo->output_width + 2) * (3 * SIZEOF(FSERROR))));
   1289     /* Might as well create the error-limiting table too. */
   1290     init_error_limit(cinfo);
   1291   }
   1292 }
   1293 
   1294 #endif /* QUANT_2PASS_SUPPORTED */
   1295