Adding upstream version 1.0.2.upstream/1.0.2upstream

Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>

1 files changed, 863 insertions, 0 deletions

diff --git a/src/3rdparty/autotrace/median.c b/src/3rdparty/autotrace/median.c
new file mode 100644
index 0000000..60e8c7f
--- /dev/null
+++ b/src/3rdparty/autotrace/median.c
@@ -0,0 +1,863 @@
+/* median.c: median cut - reducing a high color bitmap to certain number of colors
+
+   Copyright (C) 2001, 2002 Martin Weber
+
+   This library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public License
+   as published by the Free Software Foundation; either version 2.1 of
+   the License, or (at your option) any later version.
+
+   This library is distributed in the hope that it will be useful, but
+   WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with this library; if not, write to the Free Software
+   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+   USA. */
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif /* Def: HAVE_CONFIG_H */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "logreport.h"
+#include "xstd.h"
+#include "quantize.h"
+
+#define MAXNUMCOLORS 256
+
+#if 0
+#define R_SCALE
+#define G_SCALE
+#define B_SCALE
+#else
+
+/* scale RGB distances by *2,*3,*1 */
+#define R_SCALE  <<1
+#define G_SCALE  *3
+#define B_SCALE
+#endif
+
+#define BITS_IN_SAMPLE	8
+
+#define R_SHIFT	(BITS_IN_SAMPLE - PRECISION_R)
+#define G_SHIFT	(BITS_IN_SAMPLE - PRECISION_G)
+#define B_SHIFT	(BITS_IN_SAMPLE - PRECISION_B)
+
+typedef struct {
+  /* The bounds of the box (inclusive); expressed as histogram indexes */
+  int Rmin, Rmax;
+  int Gmin, Gmax;
+  int Bmin, Bmax;
+  /* The volume (actually 2-norm) of the box */
+  int volume;
+  /* The number of nonzero histogram cells within this box */
+  long colorcount;
+} box, *boxptr;
+
+static void zero_histogram_rgb(Histogram histogram)
+{
+  int r, g, b;
+  for (r = 0; r < HIST_R_ELEMS; r++)
+    for (g = 0; g < HIST_G_ELEMS; g++)
+      for (b = 0; b < HIST_B_ELEMS; b++)
+        histogram[r * MR + g * MG + b] = 0;
+}
+
+static void generate_histogram_rgb(Histogram histogram, at_bitmap * image, const at_color * ignoreColor)
+{
+  unsigned char *src = image->bitmap;
+  int num_elems;
+  ColorFreq *col;
+
+  num_elems = AT_BITMAP_WIDTH(image) * AT_BITMAP_HEIGHT(image);
+  zero_histogram_rgb(histogram);
+
+  switch (AT_BITMAP_PLANES(image)) {
+  case 3:
+    while (num_elems--) {
+      /* If we have an ignorecolor, skip it. */
+      if (ignoreColor) {
+        if ((src[0] == ignoreColor->r)
+            && (src[1] == ignoreColor->g)
+            && (src[2] == ignoreColor->b)) {
+          src += 3;
+          continue;
+        }
+      }
+      col = &histogram[(src[0] >> R_SHIFT) * MR + (src[1] >> G_SHIFT) * MG + (src[2] >> B_SHIFT)];
+      (*col)++;
+      src += 3;
+    }
+    break;
+
+  case 1:
+    while (--num_elems >= 0) {
+      if (ignoreColor && src[num_elems] == ignoreColor->r)
+        continue;
+      col = &histogram[(src[num_elems] >> R_SHIFT) * MR + (src[num_elems] >> G_SHIFT) * MG + (src[num_elems] >> B_SHIFT)];
+      (*col)++;
+    }
+    break;
+  default:
+    /* To avoid compiler warning */ ;
+  }
+}
+
+static boxptr find_biggest_volume(boxptr boxlist, int numboxes)
+/* Find the splittable box with the largest (scaled) volume */
+/* Returns 0 if no splittable boxes remain */
+{
+  boxptr boxp;
+  int i;
+  int maxv = 0;
+  boxptr which = 0;
+
+  for (i = 0, boxp = boxlist; i < numboxes; i++, boxp++) {
+    if (boxp->volume > maxv) {
+      which = boxp;
+      maxv = boxp->volume;
+    }
+  }
+
+  return which;
+}
+
+static void update_box_rgb(Histogram histogram, boxptr boxp)
+/* Shrink the min/max bounds of a box to enclose only nonzero elements, */
+/* and recompute its volume and population */
+{
+  ColorFreq *histp;
+  int R, G, B;
+  int Rmin, Rmax, Gmin, Gmax, Bmin, Bmax;
+  int dist0, dist1, dist2;
+  long ccount;
+
+  Rmin = boxp->Rmin;
+  Rmax = boxp->Rmax;
+  Gmin = boxp->Gmin;
+  Gmax = boxp->Gmax;
+  Bmin = boxp->Bmin;
+  Bmax = boxp->Bmax;
+
+  if (Rmax > Rmin)
+    for (R = Rmin; R <= Rmax; R++)
+      for (G = Gmin; G <= Gmax; G++) {
+        histp = histogram + R * MR + G * MG + Bmin;
+        for (B = Bmin; B <= Bmax; B++)
+          if (*histp++ != 0) {
+            boxp->Rmin = Rmin = R;
+            goto have_Rmin;
+          }
+      }
+have_Rmin:
+  if (Rmax > Rmin)
+    for (R = Rmax; R >= Rmin; R--)
+      for (G = Gmin; G <= Gmax; G++) {
+        histp = histogram + R * MR + G * MG + Bmin;
+        for (B = Bmin; B <= Bmax; B++)
+          if (*histp++ != 0) {
+            boxp->Rmax = Rmax = R;
+            goto have_Rmax;
+          }
+      }
+have_Rmax:
+  if (Gmax > Gmin)
+    for (G = Gmin; G <= Gmax; G++)
+      for (R = Rmin; R <= Rmax; R++) {
+        histp = histogram + R * MR + G * MG + Bmin;
+        for (B = Bmin; B <= Bmax; B++)
+          if (*histp++ != 0) {
+            boxp->Gmin = Gmin = G;
+            goto have_Gmin;
+          }
+      }
+have_Gmin:
+  if (Gmax > Gmin)
+    for (G = Gmax; G >= Gmin; G--)
+      for (R = Rmin; R <= Rmax; R++) {
+        histp = histogram + R * MR + G * MG + Bmin;
+        for (B = Bmin; B <= Bmax; B++)
+          if (*histp++ != 0) {
+            boxp->Gmax = Gmax = G;
+            goto have_Gmax;
+          }
+      }
+have_Gmax:
+  if (Bmax > Bmin)
+    for (B = Bmin; B <= Bmax; B++)
+      for (R = Rmin; R <= Rmax; R++) {
+        histp = histogram + R * MR + Gmin * MG + B;
+        for (G = Gmin; G <= Gmax; G++, histp += MG)
+          if (*histp != 0) {
+            boxp->Bmin = Bmin = B;
+            goto have_Bmin;
+          }
+      }
+have_Bmin:
+  if (Bmax > Bmin)
+    for (B = Bmax; B >= Bmin; B--)
+      for (R = Rmin; R <= Rmax; R++) {
+        histp = histogram + R * MR + Gmin * MG + B;
+        for (G = Gmin; G <= Gmax; G++, histp += MG)
+          if (*histp != 0) {
+            boxp->Bmax = Bmax = B;
+            goto have_Bmax;
+          }
+      }
+have_Bmax:
+
+  /* Update box volume.
+   * We use 2-norm rather than real volume here; this biases the method
+   * against making long narrow boxes, and it has the side benefit that
+   * a box is splittable iff norm > 0.
+   * Since the differences are expressed in histogram-cell units,
+   * we have to shift back to JSAMPLE units to get consistent distances;
+   * after which, we scale according to the selected distance scale factors.
+   */
+  dist0 = Rmax - Rmin;
+  dist1 = Gmax - Gmin;
+  dist2 = Bmax - Bmin;
+  boxp->volume = dist0 * dist0 + dist1 * dist1 + dist2 * dist2;
+
+  /* Now scan remaining volume of box and compute population */
+  ccount = 0;
+  for (R = Rmin; R <= Rmax; R++)
+    for (G = Gmin; G <= Gmax; G++) {
+      histp = histogram + R * MR + G * MG + Bmin;
+      for (B = Bmin; B <= Bmax; B++, histp++)
+        if (*histp != 0) {
+          ccount++;
+        }
+    }
+
+  boxp->colorcount = ccount;
+}
+
+static int median_cut_rgb(Histogram histogram, boxptr boxlist, int numboxes, int desired_colors)
+/* Repeatedly select and split the largest box until we have enough boxes */
+{
+  int n, lb;
+  int R, G, B, cmax;
+  boxptr b1, b2;
+
+  while (numboxes < desired_colors) {
+    /* Select box to split.
+     * Current algorithm: by population for first half, then by volume.
+     */
+    b1 = find_biggest_volume(boxlist, numboxes);
+
+    if (b1 == 0)                /* no splittable boxes left! */
+      break;
+    b2 = boxlist + numboxes;    /* where new box will go */
+    /* Copy the color bounds to the new box. */
+    b2->Rmax = b1->Rmax;
+    b2->Gmax = b1->Gmax;
+    b2->Bmax = b1->Bmax;
+    b2->Rmin = b1->Rmin;
+    b2->Gmin = b1->Gmin;
+    b2->Bmin = b1->Bmin;
+    /* Choose which axis to split the box on.
+     * Current algorithm: longest scaled axis.
+     * See notes in update_box about scaling distances.
+     */
+    R = b1->Rmax - b1->Rmin;
+    G = b1->Gmax - b1->Gmin;
+    B = b1->Bmax - b1->Bmin;
+    /* We want to break any ties in favor of green, then red, blue last.
+     */
+    cmax = G;
+    n = 1;
+    if (R > cmax) {
+      cmax = R;
+      n = 0;
+    }
+    if (B > cmax) {
+      n = 2;
+    }
+    /* Choose split point along selected axis, and update box bounds.
+     * Current algorithm: split at halfway point.
+     * (Since the box has been shrunk to minimum volume,
+     * any split will produce two nonempty subboxes.)
+     * Note that lb value is max for lower box, so must be < old max.
+     */
+    switch (n) {
+    case 0:
+      lb = (b1->Rmax + b1->Rmin) / 2;
+      b1->Rmax = lb;
+      b2->Rmin = lb + 1;
+      break;
+    case 1:
+      lb = (b1->Gmax + b1->Gmin) / 2;
+      b1->Gmax = lb;
+      b2->Gmin = lb + 1;
+      break;
+    case 2:
+      lb = (b1->Bmax + b1->Bmin) / 2;
+      b1->Bmax = lb;
+      b2->Bmin = lb + 1;
+      break;
+    }
+    /* Update stats for boxes */
+    update_box_rgb(histogram, b1);
+    update_box_rgb(histogram, b2);
+    numboxes++;
+  }
+  return numboxes;
+}
+
+static void compute_color_rgb(QuantizeObj * quantobj, Histogram histogram, boxptr boxp, int icolor)
+/* Compute representative color for a box, put it in colormap[icolor] */
+{
+  /* Current algorithm: mean weighted by pixels (not colors) */
+  /* Note it is important to get the rounding correct! */
+  ColorFreq *histp;
+  int R, G, B;
+  int Rmin, Rmax;
+  int Gmin, Gmax;
+  int Bmin, Bmax;
+  unsigned long count;
+  unsigned long total = 0;
+  unsigned long Rtotal = 0;
+  unsigned long Gtotal = 0;
+  unsigned long Btotal = 0;
+
+  Rmin = boxp->Rmin;
+  Rmax = boxp->Rmax;
+  Gmin = boxp->Gmin;
+  Gmax = boxp->Gmax;
+  Bmin = boxp->Bmin;
+  Bmax = boxp->Bmax;
+
+  for (R = Rmin; R <= Rmax; R++)
+    for (G = Gmin; G <= Gmax; G++) {
+      histp = histogram + R * MR + G * MG + Bmin;
+      for (B = Bmin; B <= Bmax; B++) {
+        if ((count = *histp++) != 0) {
+          total += count;
+          Rtotal += ((R << R_SHIFT) + ((1 << R_SHIFT) >> 1)) * count;
+          Gtotal += ((G << G_SHIFT) + ((1 << G_SHIFT) >> 1)) * count;
+          Btotal += ((B << B_SHIFT) + ((1 << B_SHIFT) >> 1)) * count;
+        }
+      }
+    }
+
+  quantobj->cmap[icolor].r = (unsigned char)((Rtotal + (total >> 1)) / total);
+  quantobj->cmap[icolor].g = (unsigned char)((Gtotal + (total >> 1)) / total);
+  quantobj->cmap[icolor].b = (unsigned char)((Btotal + (total >> 1)) / total);
+  quantobj->freq[icolor] = total;
+}
+
+static void select_colors_rgb(QuantizeObj * quantobj, Histogram histogram)
+/* Master routine for color selection */
+{
+  boxptr boxlist;
+  int numboxes;
+  int desired = quantobj->desired_number_of_colors;
+  int i;
+
+  /* Allocate workspace for box list */
+  XMALLOC(boxlist, desired * sizeof(box));
+
+  /* Initialize one box containing whole space */
+  numboxes = 1;
+  boxlist[0].Rmin = 0;
+  boxlist[0].Rmax = (1 << PRECISION_R) - 1;
+  boxlist[0].Gmin = 0;
+  boxlist[0].Gmax = (1 << PRECISION_G) - 1;
+  boxlist[0].Bmin = 0;
+  boxlist[0].Bmax = (1 << PRECISION_B) - 1;
+  /* Shrink it to actually-used volume and set its statistics */
+  update_box_rgb(histogram, boxlist);
+  /* Perform median-cut to produce final box list */
+  numboxes = median_cut_rgb(histogram, boxlist, numboxes, desired);
+  quantobj->actual_number_of_colors = numboxes;
+  /* Compute the representative color for each box, fill colormap */
+  for (i = 0; i < numboxes; i++)
+    compute_color_rgb(quantobj, histogram, boxlist + i, i);
+  free(boxlist);
+}
+
+/*
+ * These routines are concerned with the time-critical task of mapping input
+ * colors to the nearest color in the selected colormap.
+ *
+ * We re-use the histogram space as an "inverse color map", essentially a
+ * cache for the results of nearest-color searches.  All colors within a
+ * histogram cell will be mapped to the same colormap entry, namely the one
+ * closest to the cell's center.  This may not be quite the closest entry to
+ * the actual input color, but it's almost as good.  A zero in the cache
+ * indicates we haven't found the nearest color for that cell yet; the array
+ * is cleared to zeroes before starting the mapping pass.  When we find the
+ * nearest color for a cell, its colormap index plus one is recorded in the
+ * cache for future use.  The pass2 scanning routines call fill_inverse_cmap
+ * when they need to use an unfilled entry in the cache.
+ *
+ * Our method of efficiently finding nearest colors is based on the "locally
+ * sorted search" idea described by Heckbert and on the incremental distance
+ * calculation described by Spencer W. Thomas in chapter III.1 of Graphics
+ * Gems II (James Arvo, ed.  Academic Press, 1991).  Thomas points out that
+ * the distances from a given colormap entry to each cell of the histogram can
+ * be computed quickly using an incremental method: the differences between
+ * distances to adjacent cells themselves differ by a constant.  This allows a
+ * fairly fast implementation of the "brute force" approach of computing the
+ * distance from every colormap entry to every histogram cell.  Unfortunately,
+ * it needs a work array to hold the best-distance-so-far for each histogram
+ * cell (because the inner loop has to be over cells, not colormap entries).
+ * The work array elements have to be ints, so the work array would need
+ * 256Kb at our recommended precision.  This is not feasible in DOS machines.
+
+[ 256*1024/4 = 65,536 ]
+
+ * To get around these problems, we apply Thomas' method to compute the
+ * nearest colors for only the cells within a small subbox of the histogram.
+ * The work array need be only as big as the subbox, so the memory usage
+ * problem is solved.  Furthermore, we need not fill subboxes that are never
+ * referenced in pass2; many images use only part of the color gamut, so a
+ * fair amount of work is saved.  An additional advantage of this
+ * approach is that we can apply Heckbert's locality criterion to quickly
+ * eliminate colormap entries that are far away from the subbox; typically
+ * three-fourths of the colormap entries are rejected by Heckbert's criterion,
+ * and we need not compute their distances to individual cells in the subbox.
+ * The speed of this approach is heavily influenced by the subbox size: too
+ * small means too much overhead, too big loses because Heckbert's criterion
+ * can't eliminate as many colormap entries.  Empirically the best subbox
+ * size seems to be about 1/512th of the histogram (1/8th in each direction).
+ *
+ * Thomas' article also describes a refined method which is asymptotically
+ * faster than the brute-force method, but it is also far more complex and
+ * cannot efficiently be applied to small subboxes.  It is therefore not
+ * useful for programs intended to be portable to DOS machines.  On machines
+ * with plenty of memory, filling the whole histogram in one shot with Thomas'
+ * refined method might be faster than the present code --- but then again,
+ * it might not be any faster, and it's certainly more complicated.
+ */
+
+/* log2(histogram cells in update box) for each axis; this can be adjusted */
+#define BOX_R_LOG  (PRECISION_R-3)
+#define BOX_G_LOG  (PRECISION_G-3)
+#define BOX_B_LOG  (PRECISION_B-3)
+
+#define BOX_R_ELEMS  (1<<BOX_R_LOG) /* # of hist cells in update box */
+#define BOX_G_ELEMS  (1<<BOX_G_LOG)
+#define BOX_B_ELEMS  (1<<BOX_B_LOG)
+
+#define BOX_R_SHIFT  (R_SHIFT + BOX_R_LOG)
+#define BOX_G_SHIFT  (G_SHIFT + BOX_G_LOG)
+#define BOX_B_SHIFT  (B_SHIFT + BOX_B_LOG)
+
+/*
+ * The next three routines implement inverse colormap filling.  They could
+ * all be folded into one big routine, but splitting them up this way saves
+ * some stack space (the mindist[] and bestdist[] arrays need not coexist)
+ * and may allow some compilers to produce better code by registerizing more
+ * inner-loop variables.
+ */
+
+static int find_nearby_colors(QuantizeObj * quantobj, int minR, int minG, int minB, int *colorlist)
+/* Locate the colormap entries close enough to an update box to be candidates
+ * for the nearest entry to some cell(s) in the update box.  The update box
+ * is specified by the center coordinates of its first cell.  The number of
+ * candidate colormap entries is returned, and their colormap indexes are
+ * placed in colorlist[].
+ * This routine uses Heckbert's "locally sorted search" criterion to select
+ * the colors that need further consideration.
+ */
+{
+  int numcolors = quantobj->actual_number_of_colors;
+  int maxR, maxG, maxB;
+  int centerR, centerG, centerB;
+  int i, x, ncolors;
+  int minmaxdist, min_dist = 0, max_dist, tdist;
+  int mindist[MAXNUMCOLORS];    /* min distance to colormap entry i */
+
+  /* Compute TRUE coordinates of update box's upper corner and center.
+   * Actually we compute the coordinates of the center of the upper-corner
+   * histogram cell, which are the upper bounds of the volume we care about.
+   * Note that since ">>" rounds down, the "center" values may be closer to
+   * min than to max; hence comparisons to them must be "<=", not "<".
+   */
+  maxR = minR + ((1 << BOX_R_SHIFT) - (1 << R_SHIFT));
+  centerR = (minR + maxR) >> 1;
+  maxG = minG + ((1 << BOX_G_SHIFT) - (1 << G_SHIFT));
+  centerG = (minG + maxG) >> 1;
+  maxB = minB + ((1 << BOX_B_SHIFT) - (1 << B_SHIFT));
+  centerB = (minB + maxB) >> 1;
+
+  /* For each color in colormap, find:
+   *  1. its minimum squared-distance to any point in the update box
+   *     (zero if color is within update box);
+   *  2. its maximum squared-distance to any point in the update box.
+   * Both of these can be found by considering only the corners of the box.
+   * We save the minimum distance for each color in mindist[];
+   * only the smallest maximum distance is of interest.
+   */
+  minmaxdist = 0x7FFFFFFFL;
+
+  for (i = 0; i < numcolors; i++) {
+    /* We compute the squared-R-distance term, then add in the other two. */
+    x = quantobj->cmap[i].r;
+    if (x < minR) {
+      tdist = (x - minR) R_SCALE;
+      min_dist = tdist * tdist;
+      tdist = (x - maxR) R_SCALE;
+      max_dist = tdist * tdist;
+    } else if (x > maxR) {
+      tdist = (x - maxR) R_SCALE;
+      min_dist = tdist * tdist;
+      tdist = (x - minR) R_SCALE;
+      max_dist = tdist * tdist;
+    } else {
+      /* within cell range so no contribution to min_dist */
+      min_dist = 0;
+      if (x <= centerR) {
+        tdist = (x - maxR) R_SCALE;
+        max_dist = tdist * tdist;
+      } else {
+        tdist = (x - minR) R_SCALE;
+        max_dist = tdist * tdist;
+      }
+    }
+
+    x = quantobj->cmap[i].g;
+    if (x < minG) {
+      tdist = (x - minG) G_SCALE;
+      min_dist += tdist * tdist;
+      tdist = (x - maxG) G_SCALE;
+      max_dist += tdist * tdist;
+    } else if (x > maxG) {
+      tdist = (x - maxG) G_SCALE;
+      min_dist += tdist * tdist;
+      tdist = (x - minG) G_SCALE;
+      max_dist += tdist * tdist;
+    } else {
+      /* within cell range so no contribution to min_dist */
+      if (x <= centerG) {
+        tdist = (x - maxG) G_SCALE;
+        max_dist += tdist * tdist;
+      } else {
+        tdist = (x - minG) G_SCALE;
+        max_dist += tdist * tdist;
+      }
+    }
+
+    x = quantobj->cmap[i].b;
+    if (x < minB) {
+      tdist = (x - minB) B_SCALE;
+      min_dist += tdist * tdist;
+      tdist = (x - maxB) B_SCALE;
+      max_dist += tdist * tdist;
+    } else if (x > maxB) {
+      tdist = (x - maxB) B_SCALE;
+      min_dist += tdist * tdist;
+      tdist = (x - minB) B_SCALE;
+      max_dist += tdist * tdist;
+    } else {
+      /* within cell range so no contribution to min_dist */
+      if (x <= centerB) {
+        tdist = (x - maxB) B_SCALE;
+        max_dist += tdist * tdist;
+      } else {
+        tdist = (x - minB) B_SCALE;
+        max_dist += tdist * tdist;
+      }
+    }
+
+    mindist[i] = min_dist;      /* save away the results */
+    if (max_dist < minmaxdist)
+      minmaxdist = max_dist;
+  }
+
+  /* Now we know that no cell in the update box is more than minmaxdist
+   * away from some colormap entry.  Therefore, only colors that are
+   * within minmaxdist of some part of the box need be considered.
+   */
+  ncolors = 0;
+  for (i = 0; i < numcolors; i++) {
+    if (mindist[i] <= minmaxdist)
+      colorlist[ncolors++] = i;
+  }
+  return ncolors;
+}
+
+static void find_best_colors(QuantizeObj * quantobj, int minR, int minG, int minB, int numcolors, int *colorlist, int *bestcolor)
+/* Find the closest colormap entry for each cell in the update box,
+  given the list of candidate colors prepared by find_nearby_colors.
+  Return the indexes of the closest entries in the bestcolor[] array.
+  This routine uses Thomas' incremental distance calculation method to
+  find the distance from a colormap entry to successive cells in the box.
+ */
+{
+  int iR, iG, iB;
+  int i, icolor;
+  int *bptr;                    /* pointer into bestdist[] array */
+  int *cptr;                    /* pointer into bestcolor[] array */
+  int dist0, dist1;             /* initial distance values */
+  int dist2;                    /* current distance in inner loop */
+  int xx0, xx1;                 /* distance increments */
+  int xx2;
+  int inR, inG, inB;            /* initial values for increments */
+
+  /* This array holds the distance to the nearest-so-far color for each cell */
+  int bestdist[BOX_R_ELEMS * BOX_G_ELEMS * BOX_B_ELEMS];
+
+  /* Initialize best-distance for each cell of the update box */
+  bptr = bestdist;
+  for (i = BOX_R_ELEMS * BOX_G_ELEMS * BOX_B_ELEMS - 1; i >= 0; i--)
+    *bptr++ = 0x7FFFFFFFL;
+
+  /* For each color selected by find_nearby_colors,
+   * compute its distance to the center of each cell in the box.
+   * If that's less than best-so-far, update best distance and color number.
+   */
+
+  /* Nominal steps between cell centers ("x" in Thomas article) */
+#define STEP_R  ((1 << R_SHIFT) R_SCALE)
+#define STEP_G  ((1 << G_SHIFT) G_SCALE)
+#define STEP_B  ((1 << B_SHIFT) B_SCALE)
+
+  for (i = 0; i < numcolors; i++) {
+    icolor = colorlist[i];
+    /* Compute (square of) distance from minR/G/B to this color */
+    inR = (minR - quantobj->cmap[icolor].r) R_SCALE;
+    dist0 = inR * inR;
+    inG = (minG - quantobj->cmap[icolor].g) G_SCALE;
+    dist0 += inG * inG;
+    inB = (minB - quantobj->cmap[icolor].b) B_SCALE;
+    dist0 += inB * inB;
+    /* Form the initial difference increments */
+    inR = inR * (2 * STEP_R) + STEP_R * STEP_R;
+    inG = inG * (2 * STEP_G) + STEP_G * STEP_G;
+    inB = inB * (2 * STEP_B) + STEP_B * STEP_B;
+    /* Now loop over all cells in box, updating distance per Thomas method */
+    bptr = bestdist;
+    cptr = bestcolor;
+    xx0 = inR;
+    for (iR = BOX_R_ELEMS - 1; iR >= 0; iR--) {
+      dist1 = dist0;
+      xx1 = inG;
+      for (iG = BOX_G_ELEMS - 1; iG >= 0; iG--) {
+        dist2 = dist1;
+        xx2 = inB;
+        for (iB = BOX_B_ELEMS - 1; iB >= 0; iB--) {
+          if (dist2 < *bptr) {
+            *bptr = dist2;
+            *cptr = icolor;
+          }
+          dist2 += xx2;
+          xx2 += 2 * STEP_B * STEP_B;
+          bptr++;
+          cptr++;
+        }
+        dist1 += xx1;
+        xx1 += 2 * STEP_G * STEP_G;
+      }
+      dist0 += xx0;
+      xx0 += 2 * STEP_R * STEP_R;
+    }
+  }
+}
+
+static void fill_inverse_cmap_rgb(QuantizeObj * quantobj, Histogram histogram, int R, int G, int B)
+/* Fill the inverse-colormap entries in the update box that contains
+ histogram cell R/G/B.  (Only that one cell MUST be filled, but
+ we can fill as many others as we wish.) */
+{
+  int minR, minG, minB;         /* lower left corner of update box */
+  int iR, iG, iB;
+  int *cptr;                    /* pointer into bestcolor[] array */
+  ColorFreq *cachep;            /* pointer into main cache array */
+  /* This array lists the candidate colormap indexes. */
+  int colorlist[MAXNUMCOLORS];
+  int numcolors;                /* number of candidate colors */
+  /* This array holds the actually closest colormap index for each cell. */
+  int bestcolor[BOX_R_ELEMS * BOX_G_ELEMS * BOX_B_ELEMS];
+
+  /* Convert cell coordinates to update box ID */
+  R >>= BOX_R_LOG;
+  G >>= BOX_G_LOG;
+  B >>= BOX_B_LOG;
+
+  /* Compute TRUE coordinates of update box's origin corner.
+   * Actually we compute the coordinates of the center of the corner
+   * histogram cell, which are the lower bounds of the volume we care about.
+   */
+  minR = (R << BOX_R_SHIFT) + ((1 << R_SHIFT) >> 1);
+  minG = (G << BOX_G_SHIFT) + ((1 << G_SHIFT) >> 1);
+  minB = (B << BOX_B_SHIFT) + ((1 << B_SHIFT) >> 1);
+
+  /* Determine which colormap entries are close enough to be candidates
+   * for the nearest entry to some cell in the update box.
+   */
+  numcolors = find_nearby_colors(quantobj, minR, minG, minB, colorlist);
+
+  /* Determine the actually nearest colors. */
+  find_best_colors(quantobj, minR, minG, minB, numcolors, colorlist, bestcolor);
+
+  /* Save the best color numbers (plus 1) in the main cache array */
+  R <<= BOX_R_LOG;              /* convert ID back to base cell indexes */
+  G <<= BOX_G_LOG;
+  B <<= BOX_B_LOG;
+  cptr = bestcolor;
+  for (iR = 0; iR < BOX_R_ELEMS; iR++) {
+    for (iG = 0; iG < BOX_G_ELEMS; iG++) {
+      cachep = &histogram[(R + iR) * MR + (G + iG) * MG + B];
+      for (iB = 0; iB < BOX_B_ELEMS; iB++) {
+        *cachep++ = (*cptr++) + 1;
+      }
+    }
+  }
+}
+
+/*  This is pass 1  */
+static void median_cut_pass1_rgb(QuantizeObj * quantobj, at_bitmap * image, const at_color * ignoreColor)
+{
+  generate_histogram_rgb(quantobj->histogram, image, ignoreColor);
+  select_colors_rgb(quantobj, quantobj->histogram);
+}
+
+/* Map some rows of pixels to the output colormapped representation. */
+static void median_cut_pass2_rgb(QuantizeObj * quantobj, at_bitmap * image, const at_color * bgColor)
+ /* This version performs no dithering */
+{
+  Histogram histogram = quantobj->histogram;
+  ColorFreq *cachep;
+  int R, G, B;
+  int origR, origG, origB;
+  int row, col;
+  int spp = AT_BITMAP_PLANES(image);
+  int width = AT_BITMAP_WIDTH(image);
+  int height = AT_BITMAP_HEIGHT(image);
+  unsigned char *src, *dest;
+  at_color bg_color = { 0xff, 0xff, 0xff };
+
+  zero_histogram_rgb(histogram);
+
+  if (bgColor) {
+    /* Find the nearest colormap entry for the background color. */
+    R = bgColor->r >> R_SHIFT;
+    G = bgColor->g >> G_SHIFT;
+    B = bgColor->b >> B_SHIFT;
+    cachep = &histogram[R * MR + G * MG + B];
+    if (*cachep == 0)
+      fill_inverse_cmap_rgb(quantobj, histogram, R, G, B);
+    bg_color = quantobj->cmap[*cachep - 1];
+  }
+
+  src = dest = image->bitmap;
+  if (spp == 3) {
+    for (row = 0; row < height; row++) {
+      for (col = 0; col < width; col++) {
+        /* get pixel value and index into the cache */
+        origR = (*src++);
+        origG = (*src++);
+        origB = (*src++);
+
+        /*
+           if (origR > 253 && origG > 253 && origB > 253)
+           {
+           (*dest++) = 255; (*dest++) = 255; (*dest++) = 255;
+           continue;
+           }
+         */
+
+        /* get pixel value and index into the cache */
+        R = origR >> R_SHIFT;
+        G = origG >> G_SHIFT;
+        B = origB >> B_SHIFT;
+        cachep = &histogram[R * MR + G * MG + B];
+        /* If we have not seen this color before, find nearest
+           colormap entry and update the cache */
+        if (*cachep == 0) {
+          fill_inverse_cmap_rgb(quantobj, histogram, R, G, B);
+        }
+        /* Now emit the colormap index for this cell */
+        dest[0] = quantobj->cmap[*cachep - 1].r;
+        dest[1] = quantobj->cmap[*cachep - 1].g;
+        dest[2] = quantobj->cmap[*cachep - 1].b;
+
+        /* If the colormap entry for this pixel is the same as the
+           background's colormap entry, set the pixel to the
+           background color. */
+        if (bgColor && (dest[0] == bg_color.r && dest[1] == bg_color.g && dest[2] == bg_color.b)) {
+          dest[0] = bgColor->r;
+          dest[1] = bgColor->g;
+          dest[2] = bgColor->b;
+        }
+        dest += 3;
+      }
+    }
+  } else if (spp == 1) {
+    long idx = width * height;
+    while (--idx >= 0) {
+      origR = src[idx];
+      R = origR >> R_SHIFT;
+      G = origR >> G_SHIFT;
+      B = origR >> B_SHIFT;
+      cachep = &histogram[R * MR + G * MG + B];
+      if (*cachep == 0)
+        fill_inverse_cmap_rgb(quantobj, histogram, R, G, B);
+
+      dest[idx] = quantobj->cmap[*cachep - 1].r;
+
+      /* If the colormap entry for this pixel is the same as the
+         background's colormap entry, set the pixel to the
+         background color. */
+      if (bgColor && dest[idx] == bg_color.r)
+        dest[idx] = bgColor->r;
+    }
+  }
+}
+
+static QuantizeObj *initialize_median_cut(int num_colors)
+{
+  QuantizeObj *quantobj;
+
+  /* Initialize the data structures */
+  XMALLOC(quantobj, sizeof(QuantizeObj));
+
+  XMALLOC(quantobj->histogram, sizeof(ColorFreq) * HIST_R_ELEMS * HIST_G_ELEMS * HIST_B_ELEMS);
+  quantobj->desired_number_of_colors = num_colors;
+
+  return quantobj;
+}
+
+void quantize(at_bitmap * image, long ncolors, const at_color * bgColor, QuantizeObj ** iQuant, at_exception_type * exp)
+{
+  QuantizeObj *quantobj;
+  unsigned int spp = AT_BITMAP_PLANES(image);
+
+  if (spp != 3 && spp != 1) {
+    LOG("quantize: %u-plane images are not supported", spp);
+    at_exception_fatal(exp, "quantize: wrong plane images are passed");
+    return;
+  }
+
+  /* If a pointer was sent in, let's use it. */
+  if (iQuant) {
+    if (*iQuant == NULL) {
+      quantobj = initialize_median_cut(ncolors);
+      median_cut_pass1_rgb(quantobj, image, bgColor);
+      *iQuant = quantobj;
+    } else
+      quantobj = *iQuant;
+  } else {
+    quantobj = initialize_median_cut(ncolors);
+    median_cut_pass1_rgb(quantobj, image, NULL);
+  }
+
+  median_cut_pass2_rgb(quantobj, image, bgColor);
+
+  if (iQuant == NULL)
+    quantize_object_free(quantobj);
+}
+
+void quantize_object_free(QuantizeObj * quantobj)
+{
+  free(quantobj->histogram);
+  free(quantobj);
+}