gimpimage-convert-indexed.c 139 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
/* GIMP - The GNU Image Manipulation Program
 * Copyright (C) 1995 Spencer Kimball and Peter Mattis
 *
 * gimpimage-convert-indexed.c
 * Copyright (C) 1997-2004 Adam D. Moss <adam@gimp.org>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program 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 General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
18
 * along with this program.  If not, see <https://www.gnu.org/licenses/>.
19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138
 */

/*
 * 2005-09-04 - Switch 'positional' dither matrix to a 32x32 Bayer,
 *  which generates results that compress somewhat better (and may look
 *  worse or better depending on what you enjoy...).  [adam@gimp.org]
 *
 * 2004-12-12 - Use a slower but much nicer technique for finding the
 *  two best colors to dither between when using fixed/positional
 *  dither methods.  Makes positional dither much less lame.  [adam@gimp.org]
 *
 * 2002-02-10 - Quantizer version 3.0 (the rest of the commit started
 *  a year ago -- whoops).  Divide colors within CIE L*a*b* space using
 *  CPercep module (cpercep.[ch]), color-match and dither likewise,
 *  change the underlying box selection criteria and division point
 *  logic, bump luminance precision upwards, etc.etc.  Generally
 *  chooses a much richer color set, especially for low numbers of
 *  colors.  n.b.: Less luminance-sloppy in straight remapping which is
 *  good for color but a bit worse for high-frequency detail (that's
 *  partly what fs-dithering is for -- use it).  [adam@gimp.org]
 *
 * 2001-03-25 - Define accessor function/macro for histogram reads and
 *  writes.  This slows us down a little because we avoid some of the
 *  dirty tricks we used when we knew that the histogram was a straight
 *  3d array, so I've recovered some of the speed loss by implementing
 *  a 5d accessor function with good locality of reference.  This change
 *  is the first step towards quantizing in a more interesting colorspace
 *  than frumpy old RGB.  [Adam]
 *
 * 2000/01/30 - Use palette_selector instead of option_menu for custom
 *  palette. Use libgimp callback functions.  [Sven]
 *
 * 99/09/01 - Created a low-bleed FS-dither option.  [Adam]
 *
 * 99/08/29 - Deterministic color dithering to arbitrary palettes.
 *  Ideal for animations that are going to be delta-optimized or simply
 *  don't want to look 'busy' in static areas.  Also a bunch of bugfixes
 *  and tweaks.  [Adam]
 *
 * 99/08/28 - Deterministic alpha dithering over layers, reduced bleeding
 *  of transparent values into opaque values, added optional stage to
 *  remove duplicate or unused color entries from final colormap. [Adam]
 *
 * 99/02/24 - Many revisions to the box-cut quantizer used in RGB->INDEXED
 *  conversion.  Box to be cut is chosen on the basis of possessing an axis
 *  with the largest sum of weighted perceptible error, rather than based on
 *  volume or population.  The box is split along this axis rather than its
 *  longest axis, at the point of error mean rather than simply at its centre.
 *  Error-limiting in the F-S dither has been disabled - it may become optional
 *  again later.  If you're convinced that you have an image where the old
 *  dither looks better, let me know.  [Adam]
 *
 * 99/01/10 - Hourglass... [Adam]
 *
 * 98/07/25 - Convert-to-indexed now remembers the last invocation's
 *  settings.  Also, GRAY->INDEXED is more flexible.  [Adam]
 *
 * 98/07/05 - Sucked the warning about quantizing to too many colors into
 *  a text widget embedded in the dialog, improved intelligence of dialog
 *  to default 'custom palette' selection to 'Web' if available, and
 *  in this case not bother to present the native WWW-palette radio
 *  button.  [Adam]
 *
 * 98/04/13 - avoid a division by zero when converting an empty gray-scale
 *  image (who would like to do such a thing anyway??)  [Sven ]
 *
 * 98/03/23 - fixed a longstanding fencepost - hopefully the *right*
 *  way, *again*.  [Adam]
 *
 * 97/11/14 - added a proper pdb interface and support for dithering
 *  to custom palettes (based on a patch by Eric Hernes) [Yosh]
 *
 * 97/11/04 - fixed the accidental use of the color-counting case
 *  when palette_type is WEB or MONO. [Adam]
 *
 * 97/10/25 - color-counting implemented (could use some hashing, but
 *  performance actually seems okay) - now RGB->INDEXED conversion isn't
 *  destructive if it doesn't have to be. [Adam]
 *
 * 97/10/14 - fixed divide-by-zero when converting a completely transparent
 *  RGB image to indexed. [Adam]
 *
 * 97/07/01 - started todo/revision log.  Put code back in to
 *  eliminate full-alpha pixels from RGB histogram.
 *  [Adam D. Moss - adam@gimp.org]
 */

  /* TODO for Convert:
   *
   * . Tweak, tweak, tweak.  Old RGB code was tuned muchly.
   *
   * . Re-enable Heckbert locality for matching, benchmark it
   *
   * . Try faster fixed-point sRGB<->L*a*b* pixel conversion (see cpercep.c)
   *
   * . Use palette of another open INDEXED image?
   *
   * . Do error-splitting trick for GREY->INDEXED (hardly worth it)
   */

  /* CODE READABILITY BUGS:
   *
   * . Most uses of variants of the R,G,B variable naming convention
   *   are referring to L*a*b* co-ordinates, not RGB co-ordinates!
   *
   * . Each said variable is usually one of the following, but it is
   *   rarely clear which one:
   *     - (assumed sRGB) raw non-linear 8-bit RGB co-ordinates
   *     - 'full'-precision (unshifted) 8-bit L*a*b* co-ordinates
   *     - box-space (reduced-precision shifted L*a*b*) co-ordinates
   */

#include "config.h"

#include <stdlib.h>
#include <stdio.h>
#include <string.h>

#include <cairo.h>
#include <gdk-pixbuf/gdk-pixbuf.h>
139
#define GEGL_ITERATOR2_API
140 141 142 143 144 145 146
#include <gegl.h>

#include "libgimpcolor/gimpcolor.h"
#include "libgimpmath/gimpmath.h"

#include "core-types.h"

147
#include "gegl/gimp-babl.h"
148 149 150 151 152 153 154 155 156 157 158 159
#include "gegl/gimp-gegl-utils.h"

#include "gimp.h"
#include "gimpcontainer.h"
#include "gimpdrawable.h"
#include "gimperror.h"
#include "gimpimage.h"
#include "gimpimage-color-profile.h"
#include "gimpimage-colormap.h"
#include "gimpimage-undo.h"
#include "gimpimage-undo-push.h"
#include "gimplayer.h"
160
#include "gimpobjectqueue.h"
161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487
#include "gimppalette.h"
#include "gimpprogress.h"

#include "text/gimptextlayer.h"

#include "gimpimage-convert-fsdither.h"
#include "gimpimage-convert-data.h"
#include "gimpimage-convert-indexed.h"

#include "gimp-intl.h"


/* basic memory/quality tradeoff */
#define PRECISION_R 8
#define PRECISION_G 6
#define PRECISION_B 6

#define HIST_R_ELEMS (1<<PRECISION_R)
#define HIST_G_ELEMS (1<<PRECISION_G)
#define HIST_B_ELEMS (1<<PRECISION_B)

#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)

/* we've stretched our non-cubic L*a*b* volume to touch the
 * faces of the logical cube we've allocated for it, so re-scale
 * again in inverse proportion to get back to linear proportions.
 */
#define R_SCALE 13              /*  scale R (L*) distances by this much  */
#define G_SCALE 24              /*  scale G (a*) distances by this much  */
#define B_SCALE 26              /*  and B (b*) by this much              */


typedef struct _Color Color;
typedef struct _QuantizeObj QuantizeObj;

typedef void (* Pass1Func)     (QuantizeObj *quantize_obj);
typedef void (* Pass2InitFunc) (QuantizeObj *quantize_obj);
typedef void (* Pass2Func)     (QuantizeObj *quantize_obj,
                                GimpLayer   *layer,
                                GeglBuffer  *new_buffer);
typedef void (* CleanupFunc)   (QuantizeObj *quantize_obj);

typedef gulong ColorFreq;
typedef ColorFreq * CFHistogram;

typedef enum { AXIS_UNDEF, AXIS_RED, AXIS_BLUE, AXIS_GREEN } AxisType;

typedef double etype;


/*
  We provide two different histogram access interfaces.  HIST_LIN()
  accesses the histogram in histogram-native space, taking absolute
  histogram co-ordinates.  HIST_RGB() accesses the histogram in RGB
  space.  This latter takes unsigned 8-bit co-ordinates, internally
  converts those co-ordinates to histogram-native space and returns
  the access pointer to the corresponding histogram cell.

  Using these two interfaces we can import RGB data into a more
  interesting space and efficiently work in the latter space until
  it is time to output the quantized values in RGB again.  For
  this final conversion we implement the function lin_to_rgb().

  We effectively pull our three-dimensional space into five dimensions
  such that the most-entropic bits lay in the lowest bits of the resulting
  array index.  This gives significantly better locality of reference
  and hence a small speedup despite the extra work involved in calculating
  the index.

  Why not six dimensions?  The expansion of dimensionality is good for random
  access such as histogram population and the query pattern typical of
  dithering but we have some code which iterates in a scanning manner, for
  which the expansion is suboptimal.  The compromise is to leave the B
  dimension unmolested in the lower-order bits of the index, since this is
  the dimension most commonly iterated through in the inner loop of the
  scans.

  --adam

  RhGhRlGlB
*/
#define VOL_GBITS  (PRECISION_G)
#define VOL_BBITS  (PRECISION_B)
#define VOL_RBITS  (PRECISION_R)
#define VOL_GBITSh (VOL_GBITS - 3)
#define VOL_GBITSl (VOL_GBITS - VOL_GBITSh)
#define VOL_BBITSh (VOL_BBITS - 4)
#define VOL_BBITSl (VOL_BBITS - VOL_BBITSh)
#define VOL_RBITSh (VOL_RBITS - 3)
#define VOL_RBITSl (VOL_RBITS - VOL_RBITSh)
#define VOL_GMASKh (((1<<VOL_GBITSh)-1) << VOL_GBITSl)
#define VOL_GMASKl ((1<<VOL_GBITSl)-1)
#define VOL_BMASKh (((1<<VOL_BBITSh)-1) << VOL_BBITSl)
#define VOL_BMASKl ((1<<VOL_BBITSl)-1)
#define VOL_RMASKh (((1<<VOL_RBITSh)-1) << VOL_RBITSl)
#define VOL_RMASKl ((1<<VOL_RBITSl)-1)
/* The 5d compromise thing. */
#define REF_FUNC(r,g,b) \
( \
 (((r) & VOL_RMASKh) << (VOL_BBITS + VOL_GBITS)) | \
 (((r) & VOL_RMASKl) << (VOL_GBITSl + VOL_BBITS)) | \
 (((g) & VOL_GMASKh) << (VOL_RBITSl + VOL_BBITS)) | \
 (((g) & VOL_GMASKl) << (VOL_BBITS)) | \
 (b) \
)
/* The full-on 6d thing. */
/*
#define REF_FUNC(r,g,b) \
( \
 (((r) & VOL_RMASKh) << (VOL_BBITS + VOL_GBITS)) | \
 (((r) & VOL_RMASKl) << (VOL_GBITSl + VOL_BBITSl)) | \
 (((g) & VOL_GMASKh) << (VOL_RBITSl + VOL_BBITS)) | \
 (((g) & VOL_GMASKl) << (VOL_BBITSl)) | \
 (((b) & VOL_BMASKh) << (VOL_RBITSl + VOL_GBITSl)) | \
 ((b) & VOL_BMASKl) \
)
*/
/* The boring old 3d thing. */
/*
#define REF_FUNC(r,g,b) (((r)<<(PRECISION_G+PRECISION_B)) | ((g)<<(PRECISION_B)) | (b))
*/

/* You even get to choose whether you want the accessor function
   implemented as a macro or an inline function.  Don't say I never
   give you anything. */
/*
#define HIST_LIN(hist_ptr,r,g,b) (&(hist_ptr)[REF_FUNC((r),(g),(b))])
*/
static inline ColorFreq *
HIST_LIN (ColorFreq  *hist_ptr,
          const gint  r,
          const gint  g,
          const gint  b)
{
  return (&(hist_ptr) [REF_FUNC (r, g, b)]);
}


#define LOWA   (-86.181F)
#define LOWB  (-107.858F)
#define HIGHA   (98.237F)
#define HIGHB   (94.480F)

#if 1
#define LRAT (2.55F)
#define ARAT (255.0F / (HIGHA - LOWA))
#define BRAT (255.0F / (HIGHB - LOWB))
#else
#define LRAT (1.0F)
#define ARAT (1.0F)
#define BRAT (1.0F)
#endif

static const Babl *rgb_to_lab_fish = NULL;
static const Babl *lab_to_rgb_fish = NULL;

static inline void
rgb_to_unshifted_lin (const guchar  r,
                      const guchar  g,
                      const guchar  b,
                      gint         *hr,
                      gint         *hg,
                      gint         *hb)
{
  gint   or, og, ob;
  gfloat rgb[3] = { r / 255.0, g / 255.0, b / 255.0 };
  gfloat lab[3];

  babl_process (rgb_to_lab_fish, rgb, lab, 1);

  /* fprintf(stderr, " %d-%d-%d -> %0.3f,%0.3f,%0.3f ", r, g, b, sL, sa, sb);*/

  or = RINT(lab[0] * LRAT);
  og = RINT((lab[1] - LOWA) * ARAT);
  ob = RINT((lab[2] - LOWB) * BRAT);

  *hr = CLAMP(or, 0, 255);
  *hg = CLAMP(og, 0, 255);
  *hb = CLAMP(ob, 0, 255);

  /*  fprintf(stderr, " %d:%d:%d ", *hr, *hg, *hb); */
}


static inline void
rgb_to_lin (const guchar  r,
            const guchar  g,
            const guchar  b,
            gint         *hr,
            gint         *hg,
            gint         *hb)
{
  gint or, og, ob;

  /*
  double sL, sa, sb;
  {
    double low_l = 999.0, low_a = 999.9, low_b = 999.0;
    double high_l = -999.0, high_a = -999.0, high_b = -999.0;

    int r,g,b;

    for (r=0; r<256; r++)
      for (g=0; g<256; g++)
        for (b=0; b<256; b++)
          {
            cpercep_rgb_to_space(r,g,b, &sL, &sa, &sb);

            if (sL > high_l)
              high_l = sL;
            if (sL < low_l)
              low_l = sL;
            if (sa > high_a)
              high_a = sa;
            if (sa < low_a)
              low_a = sa;
            if (sb > high_b)
              high_b = sb;
            if (sb < low_b)
              low_b = sb;
          }

    fprintf(stderr, " [L: %0.3f -> %0.3f / a: %0.3f -> %0.3f / b: %0.3f -> %0.3f]\t", low_l, high_l, low_a, high_a, low_b, high_b);

    exit(-1);
  }
  */

  rgb_to_unshifted_lin (r, g, b, &or, &og, &ob);

#if 0
#define RSDF(r) ((r) >= ((HIST_R_ELEMS-1) << R_SHIFT) ? HIST_R_ELEMS-1 : \
                 ((r) + ((1<<R_SHIFT)>>1) ) >> R_SHIFT)
#define GSDF(g) ((g) >= ((HIST_G_ELEMS-1) << G_SHIFT) ? HIST_G_ELEMS-1 : \
                 ((g) + ((1<<G_SHIFT)>>1) ) >> G_SHIFT)
#define BSDF(b) ((b) >= ((HIST_B_ELEMS-1) << B_SHIFT) ? HIST_B_ELEMS-1 : \
                 ((b) + ((1<<B_SHIFT)>>1) ) >> B_SHIFT)
#else
#define RSDF(r) ((r) >> R_SHIFT)
#define GSDF(g) ((g) >> G_SHIFT)
#define BSDF(b) ((b) >> B_SHIFT)
#endif

  or = RSDF (or);
  og = GSDF (og);
  ob = BSDF (ob);

  *hr = or;
  *hg = og;
  *hb = ob;
}


static inline ColorFreq *
HIST_RGB (ColorFreq  *hist_ptr,
          const gint  r,
          const gint  g,
          const gint  b)
{
  gint hr, hg, hb;

  rgb_to_lin (r, g, b, &hr, &hg, &hb);

  return HIST_LIN (hist_ptr, hr, hg, hb);
}


static inline void
lin_to_rgb (const gdouble  hr,
            const gdouble  hg,
            const gdouble  hb,
            guchar        *r,
            guchar        *g,
            guchar        *b)
{
  gfloat  rgb[3];
  gfloat  lab[3];
  gdouble ir, ig, ib;

  ir = ((gdouble) (hr)) * 255.0F / (gdouble) (HIST_R_ELEMS - 1);
  ig = ((gdouble)( hg)) * 255.0F / (gdouble) (HIST_G_ELEMS - 1);
  ib = ((gdouble)( hb)) * 255.0F / (gdouble) (HIST_B_ELEMS - 1);

  ir = ir / LRAT;
  ig = (ig / ARAT) + LOWA;
  ib = (ib / BRAT) + LOWB;

  lab[0] = ir;
  lab[1] = ig;
  lab[2] = ib;

  babl_process (lab_to_rgb_fish, lab, rgb, 1);

  *r = RINT (CLAMP (rgb[0] * 255, 0.0F, 255.0F));
  *g = RINT (CLAMP (rgb[1] * 255, 0.0F, 255.0F));
  *b = RINT (CLAMP (rgb[2] * 255, 0.0F, 255.0F));
}



struct _Color
{
  gint red;
  gint green;
  gint blue;
};

struct _QuantizeObj
{
  Pass1Func     first_pass;       /* first pass over image data creates colormap  */
  Pass2InitFunc second_pass_init; /* Initialize data which persists over invocations */
  Pass2Func     second_pass;      /* second pass maps from image data to colormap */
  CleanupFunc   delete_func;      /* function to clean up data associated with private */

  GimpPalette  *custom_palette;           /* The custom palette, if any        */

  gint          desired_number_of_colors; /* Number of colors we will allow    */
  gint          actual_number_of_colors;  /* Number of colors actually needed  */
  Color         cmap[256];                /* colormap created by quantization  */
  Color         clin[256];                /* .. converted back to linear space */
  gulong        index_used_count[256];    /* how many times an index was used  */
  CFHistogram   histogram;                /* holds the histogram               */

488
  gboolean      want_dither_alpha;
489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523
  gint          error_freedom;            /* 0=much bleed, 1=controlled bleed */

  GimpProgress *progress;
};

typedef struct
{
  /*  The bounds of the box (inclusive); expressed as histogram indexes  */
  gint    Rmin, Rmax;
  gint    Rhalferror;
  gint    Gmin, Gmax;
  gint    Ghalferror;
  gint    Bmin, Bmax;
  gint    Bhalferror;

  /*  The volume (actually 2-norm) of the box  */
  gint    volume;

  /*  The number of nonzero histogram cells within this box  */
  glong   colorcount;

  /* The sum of the weighted error within this box */
  guint64 error;
  /* The sum of the unweighted error within this box */
  guint64 rerror;
  guint64 gerror;
  guint64 berror;

} box, *boxptr;


static void          zero_histogram_gray     (CFHistogram   histogram);
static void          zero_histogram_rgb      (CFHistogram   histogram);
static void          generate_histogram_gray (CFHistogram   hostogram,
                                              GimpLayer    *layer,
524
                                              gboolean      dither_alpha);
525 526 527
static void          generate_histogram_rgb  (CFHistogram   histogram,
                                              GimpLayer    *layer,
                                              gint          col_limit,
528
                                              gboolean      dither_alpha,
529
                                              GimpProgress *progress);
530 531

static QuantizeObj * initialize_median_cut   (GimpImageBaseType      old_type,
532
                                              gint                   max_colors,
533 534 535
                                              GimpConvertDitherType  dither_type,
                                              GimpConvertPaletteType palette_type,
                                              GimpPalette           *custom_palette,
536
                                              gboolean               dither_alpha,
537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634
                                              GimpProgress          *progress);

static void          compute_color_lin8      (QuantizeObj           *quantobj,
                                              CFHistogram            histogram,
                                              boxptr                 boxp,
                                              const int              icolor);


static guchar    found_cols[MAXNUMCOLORS][3];
static gint      num_found_cols;
static gboolean  needs_quantize;


/**********************************************************/
typedef struct
{
  glong  used_count;
  guchar initial_index;
} PalEntry;

static int
mapping_compare (const void *a,
                 const void *b)
{
  PalEntry *m1 = (PalEntry *) a;
  PalEntry *m2 = (PalEntry *) b;

  return (m2->used_count - m1->used_count);
}

/* FWIW, the make_remap_table() and mapping_compare() function source
 * and PalEntry may be re-used under the XFree86-style license.
 * <adam@gimp.org>
 */
static void
make_remap_table (const guchar  old_palette[],
                  guchar        new_palette[],
                  const gulong  index_used_count[],
                  guchar        remap_table[],
                  gint         *num_entries)
{
  gint      i, j, k;
  guchar    temppal[256 * 3];
  gulong    tempuse[256];
  gulong    transmap[256];
  PalEntry *palentries;
  gint      used = 0;

  memset (temppal, 0, 256 * 3);
  memset (tempuse, 0, 256 * sizeof (gulong));
  memset (transmap, 255, 256 * sizeof (gulong));

  /* First pass - only collect entries which are marked as being used
   * at all in index_used_count.
   */
  for (i = 0; i < *num_entries; i++)
    {
      if (index_used_count[i])
        {
          temppal[used*3 + 0] = old_palette[i*3 + 0];
          temppal[used*3 + 1] = old_palette[i*3 + 1];
          temppal[used*3 + 2] = old_palette[i*3 + 2];

          tempuse[used] = index_used_count[i];
          transmap[i] = used;

          used++;
        }
    }

  /* Second pass - remove duplicates. (O(n^3), could do better!) */
  for (i = 0; i < used; i++)
    {
      for (j = 0; j < i; j++)
        {
          if ((temppal[i*3 + 1] == temppal[j*3 + 1]) &&
              (temppal[i*3 + 0] == temppal[j*3 + 0]) &&
              (temppal[i*3 + 2] == temppal[j*3 + 2]) &&
              tempuse[j] &&
              tempuse[i])
            {
              /* Move the 'used' tally from one to the other. */
              tempuse[i] += tempuse[j];
              /* zero one of them, deactivating its entry. */
              tempuse[j] = 0;

              /* change all mappings from this dead index to the live
               * one.
               */
              for (k = 0; k < *num_entries; k++)
                {
                  if (index_used_count[k] && (transmap[k] == j))
                    transmap[k] = i;
                }
            }
        }
    }

635
  /* Third pass - rank all used indices to the beginning of the
636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702
   * palette.
   */
  palentries = g_new (PalEntry, used);

  for (i = 0; i < used; i++)
    {
      palentries[i].initial_index = i;
      palentries[i].used_count    = tempuse[i];
    }

  qsort (palentries, used, sizeof (PalEntry), &mapping_compare);

  for (i = 0; i < *num_entries; i++)
    {
      if (index_used_count[i])
        {
          for (j = 0; j < used; j++)
            {
              if ((transmap[i] == palentries[j].initial_index)
                  && (palentries[j].used_count))
                {
                  remap_table[i] = j;
                  break;
                }
            }
        }
    }
  for (i = 0; i < *num_entries; i++)
    {
      if (index_used_count[i])
        {
          new_palette[remap_table[i] * 3 + 0] = old_palette[i * 3 + 0];
          new_palette[remap_table[i] * 3 + 1] = old_palette[i * 3 + 1];
          new_palette[remap_table[i] * 3 + 2] = old_palette[i * 3 + 2];
        }
    }

  *num_entries = 0;

  for (j = 0; j < used; j++)
    {
      if (palentries[j].used_count)
        {
          (*num_entries)++;
        }
    }

  g_free (palentries);
}

static void
remap_indexed_layer (GimpLayer    *layer,
                     const guchar *remap_table,
                     gint          num_entries)
{
  GeglBufferIterator *iter;
  const Babl         *format;
  gint                bpp;
  gboolean            has_alpha;

  format  = gimp_drawable_get_format (GIMP_DRAWABLE (layer));

  bpp       = babl_format_get_bytes_per_pixel (format);
  has_alpha = babl_format_has_alpha (format);

  iter = gegl_buffer_iterator_new (gimp_drawable_get_buffer (GIMP_DRAWABLE (layer)),
                                   NULL, 0, NULL,
703
                                   GEGL_ACCESS_READWRITE, GEGL_ABYSS_NONE, 1);
704 705 706

  while (gegl_buffer_iterator_next (iter))
    {
707
      guchar *data   = iter->items[0].data;
708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752
      gint    length = iter->length;

      if (has_alpha)
        {
          while (length--)
            {
              if (data[ALPHA_I])
                data[INDEXED] = remap_table[data[INDEXED]];

              data += bpp;
            }
        }
      else
        {
          while (length--)
            {
              data[INDEXED] = remap_table[data[INDEXED]];

              data += bpp;
            }
        }
    }
}

static gint
color_quicksort (const void *c1,
                 const void *c2)
{
  Color *color1 = (Color *) c1;
  Color *color2 = (Color *) c2;

  gdouble v1 = GIMP_RGB_LUMINANCE (color1->red, color1->green, color1->blue);
  gdouble v2 = GIMP_RGB_LUMINANCE (color2->red, color2->green, color2->blue);

  if (v1 < v2)
    return -1;
  else if (v1 > v2)
    return 1;
  else
    return 0;
}

gboolean
gimp_image_convert_indexed (GimpImage               *image,
                            GimpConvertPaletteType   palette_type,
753 754 755 756 757
                            gint                     max_colors,
                            gboolean                 remove_duplicates,
                            GimpConvertDitherType    dither_type,
                            gboolean                 dither_alpha,
                            gboolean                 dither_text_layers,
758 759 760 761
                            GimpPalette             *custom_palette,
                            GimpProgress            *progress,
                            GError                 **error)
{
762 763 764
  QuantizeObj       *quantobj     = NULL;
  GimpObjectQueue   *queue        = NULL;
  GimpProgress      *sub_progress = NULL;
765 766 767
  GimpImageBaseType  old_type;
  GList             *all_layers;
  GList             *list;
768
  GimpColorProfile  *src_profile  = NULL;
769 770 771 772
  GimpColorProfile  *dest_profile = NULL;

  g_return_val_if_fail (GIMP_IS_IMAGE (image), FALSE);
  g_return_val_if_fail (gimp_image_get_base_type (image) != GIMP_INDEXED, FALSE);
773 774 775
  g_return_val_if_fail (gimp_babl_is_valid (GIMP_INDEXED,
                                            gimp_image_get_precision (image)),
                        FALSE);
776 777 778 779 780 781 782 783
  g_return_val_if_fail (custom_palette == NULL ||
                        GIMP_IS_PALETTE (custom_palette), FALSE);
  g_return_val_if_fail (custom_palette == NULL ||
                        gimp_palette_get_n_colors (custom_palette) <= 256,
                        FALSE);
  g_return_val_if_fail (progress == NULL || GIMP_IS_PROGRESS (progress), FALSE);
  g_return_val_if_fail (error == NULL || *error == NULL, FALSE);

784
  if (palette_type == GIMP_CONVERT_PALETTE_CUSTOM)
785 786
    {
      if (! custom_palette)
787
        palette_type = GIMP_CONVERT_PALETTE_MONO;
788 789 790 791

      if (gimp_palette_get_n_colors (custom_palette) == 0)
        {
          g_set_error_literal (error, GIMP_ERROR, GIMP_FAILED,
792
                               _("Cannot convert image: palette is empty."));
793 794 795 796 797 798 799 800 801 802 803 804 805
          return FALSE;
        }
    }

  gimp_set_busy (image->gimp);

  all_layers = gimp_image_get_layer_list (image);

  g_object_freeze_notify (G_OBJECT (image));

  gimp_image_undo_group_start (image, GIMP_UNDO_GROUP_IMAGE_CONVERT,
                               C_("undo-type", "Convert Image to Indexed"));

806 807
  src_profile = gimp_color_managed_get_color_profile (GIMP_COLOR_MANAGED (image));

808 809 810 811 812 813 814 815
  /*  Push the image type to the stack  */
  gimp_image_undo_push_image_type (image, NULL);

  /*  Set the new base type  */
  old_type = gimp_image_get_base_type (image);

  g_object_set (image, "base-type", GIMP_INDEXED, NULL);

816 817
  /* when converting from GRAY, always convert to the new type's
   * builtin profile
818 819
   */
  if (old_type == GIMP_GRAY)
820
    dest_profile = gimp_image_get_builtin_color_profile (image);
821 822 823 824 825 826 827 828 829 830 831

  /*  Build histogram if necessary.  */
  rgb_to_lab_fish = babl_fish (babl_format ("R'G'B' float"),
                               babl_format ("CIE Lab float"));
  lab_to_rgb_fish = babl_fish (babl_format ("CIE Lab float"),
                               babl_format ("R'G'B' float"));

  /* don't dither if the input is grayscale and we are simply mapping
   * every color
   */
  if (old_type     == GIMP_GRAY &&
832
      max_colors   == 256       &&
833
      palette_type == GIMP_CONVERT_PALETTE_GENERATE)
834
    {
835
      dither_type = GIMP_CONVERT_DITHER_NONE;
836 837
    }

838 839 840 841 842 843 844 845
  if (progress)
    {
      queue        = gimp_object_queue_new (progress);
      sub_progress = GIMP_PROGRESS (queue);

      gimp_object_queue_push_list (queue, all_layers);
    }

846
  quantobj = initialize_median_cut (old_type, max_colors, dither_type,
847
                                    palette_type, custom_palette,
848
                                    dither_alpha,
849
                                    sub_progress);
850

851
  if (palette_type == GIMP_CONVERT_PALETTE_GENERATE)
852 853 854 855 856 857 858 859 860 861 862 863 864 865
    {
      if (old_type == GIMP_GRAY)
        zero_histogram_gray (quantobj->histogram);
      else
        zero_histogram_rgb (quantobj->histogram);

      /* To begin, assume that there are fewer colors in the image
       *  than the user actually asked for.  In that case, we don't
       *  need to quantize or color-dither.
       */
      needs_quantize = FALSE;
      num_found_cols = 0;

      /*  Build the histogram  */
866
      for (list = all_layers;
867
           list;
868
           list = g_list_next (list))
869 870 871
        {
          GimpLayer *layer = list->data;

872 873 874
          if (queue)
            gimp_object_queue_pop (queue);

875 876 877
          if (old_type == GIMP_GRAY)
            {
              generate_histogram_gray (quantobj->histogram,
878
                                       layer, dither_alpha);
879 880 881 882 883 884 885 886
            }
          else
            {
              /* Note: generate_histogram_rgb may set needs_quantize
               * if the image contains more colors than the limit
               * specified by the user.
               */
              generate_histogram_rgb (quantobj->histogram,
887
                                      layer, max_colors, dither_alpha,
888
                                      sub_progress);
889 890 891 892 893 894 895 896 897 898
            }
        }
    }

  if (progress)
    gimp_progress_set_text_literal (progress,
                                    _("Converting to indexed colors (stage 2)"));

  if (old_type == GIMP_RGB &&
      ! needs_quantize     &&
899
      palette_type == GIMP_CONVERT_PALETTE_GENERATE)
900 901 902 903 904 905 906 907 908 909 910 911 912 913 914
    {
      gint i;

      /*  If this is an RGB image, and the user wanted a custom-built
       *  generated palette, and this image has no more colors than
       *  the user asked for, we don't need the first pass
       *  (quantization).
       *
       *  There's also no point in dithering, since there's no error
       *  to spread.  So we destroy the old quantobj and make a new
       *  one with the remapping function set to a special LUT-based
       *  no-dither remapper.
       */

      quantobj->delete_func (quantobj);
915
      quantobj = initialize_median_cut (old_type, max_colors,
916
                                        GIMP_CONVERT_DITHER_NODESTRUCT,
917 918
                                        palette_type,
                                        custom_palette,
919
                                        dither_alpha,
920
                                        sub_progress);
921 922 923 924 925 926 927 928 929 930 931 932 933 934 935
      /* We can skip the first pass (palette creation) */

      quantobj->actual_number_of_colors = num_found_cols;
      for (i = 0; i < num_found_cols; i++)
        {
          quantobj->cmap[i].red   = found_cols[i][0];
          quantobj->cmap[i].green = found_cols[i][1];
          quantobj->cmap[i].blue  = found_cols[i][2];
        }
    }
  else
    {
      quantobj->first_pass (quantobj);
    }

936
  if (palette_type == GIMP_CONVERT_PALETTE_GENERATE)
937 938 939 940 941
    qsort (quantobj->cmap,
           quantobj->actual_number_of_colors, sizeof (Color),
           color_quicksort);

  if (progress)
942 943 944 945 946 947 948
    {
      gimp_progress_set_text_literal (progress,
                                      _("Converting to indexed colors (stage 3)"));

      gimp_object_queue_clear (queue);
      gimp_object_queue_push_list (queue, all_layers);
    }
949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973

  /* Initialise data which must persist across indexed layer iterations */
  if (quantobj->second_pass_init)
    quantobj->second_pass_init (quantobj);

  /*  Set the generated palette on the image, we need it to
   *  convert the layers. We optionally remove duplicate entries
   *  after the layer conversion.
   */
  {
    guchar colormap[GIMP_IMAGE_COLORMAP_SIZE];
    gint   i, j;

    for (i = 0, j = 0; i < quantobj->actual_number_of_colors; i++)
      {
        colormap[j++] = quantobj->cmap[i].red;
        colormap[j++] = quantobj->cmap[i].green;
        colormap[j++] = quantobj->cmap[i].blue;
      }

    gimp_image_set_colormap (image, colormap,
                             quantobj->actual_number_of_colors, TRUE);
  }

  /*  Convert all layers  */
974
  for (list = all_layers;
975
       list;
976
       list = g_list_next (list))
977 978 979 980
    {
      GimpLayer *layer = list->data;
      gboolean   quantize;

981 982 983
      if (queue)
        gimp_object_queue_pop (queue);

984
      if (gimp_item_is_text_layer (GIMP_ITEM (layer)))
985
        quantize = dither_text_layers;
986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010
      else
        quantize = TRUE;

      if (quantize)
        {
          GeglBuffer *new_buffer;
          gboolean    has_alpha;

          has_alpha = gimp_drawable_has_alpha (GIMP_DRAWABLE (layer));

          new_buffer =
            gegl_buffer_new (GEGL_RECTANGLE (0, 0,
                                             gimp_item_get_width  (GIMP_ITEM (layer)),
                                             gimp_item_get_height (GIMP_ITEM (layer))),
                             gimp_image_get_layer_format (image,
                                                          has_alpha));

          quantobj->second_pass (quantobj, layer, new_buffer);

          gimp_drawable_set_buffer (GIMP_DRAWABLE (layer), TRUE, NULL,
                                    new_buffer);
          g_object_unref (new_buffer);
        }
      else
        {
1011 1012
          gimp_drawable_convert_type (GIMP_DRAWABLE (layer), image,
                                      GIMP_INDEXED,
1013
                                      gimp_drawable_get_precision (GIMP_DRAWABLE (layer)),
1014
                                      gimp_drawable_has_alpha (GIMP_DRAWABLE (layer)),
1015
                                      src_profile,
1016
                                      dest_profile,
1017
                                      GEGL_DITHER_NONE, GEGL_DITHER_NONE,
1018
                                      TRUE, sub_progress);
1019 1020 1021 1022
        }
    }

  /*  Set the final palette on the image  */
1023 1024 1025
  if (remove_duplicates &&
      (palette_type != GIMP_CONVERT_PALETTE_GENERATE) &&
      (palette_type != GIMP_CONVERT_PALETTE_MONO))
1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066
    {
      guchar colormap[GIMP_IMAGE_COLORMAP_SIZE];
      gint   i, j;
      guchar old_palette[256 * 3];
      guchar new_palette[256 * 3];
      guchar remap_table[256];
      gint   num_entries;

      for (i = 0, j = 0; i < quantobj->actual_number_of_colors; i++)
        {
          old_palette[j++] = quantobj->cmap[i].red;
          old_palette[j++] = quantobj->cmap[i].green;
          old_palette[j++] = quantobj->cmap[i].blue;
        }

      num_entries = quantobj->actual_number_of_colors;

      /* Generate a remapping table */
      make_remap_table (old_palette, new_palette,
                        quantobj->index_used_count,
                        remap_table, &num_entries);

      /*  Convert all layers  */
      for (list = all_layers; list; list = g_list_next (list))
        {
          remap_indexed_layer (list->data, remap_table, num_entries);
        }

      for (i = 0, j = 0; i < num_entries; i++)
        {
          colormap[j] = new_palette[j]; j++;
          colormap[j] = new_palette[j]; j++;
          colormap[j] = new_palette[j]; j++;
        }

      gimp_image_set_colormap (image, colormap, num_entries, TRUE);
    }

  /*  When converting from GRAY, set the new profile.
   */
  if (old_type == GIMP_GRAY)
1067
    gimp_image_set_color_profile (image, dest_profile, NULL);
1068 1069 1070 1071 1072 1073 1074 1075 1076 1077

  /*  Delete the quantizer object, if there is one */
  if (quantobj)
    quantobj->delete_func (quantobj);

  gimp_image_undo_group_end (image);

  gimp_image_mode_changed (image);
  g_object_thaw_notify (G_OBJECT (image));

1078 1079
  g_clear_object (&queue);

1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111
  g_list_free (all_layers);

  gimp_unset_busy (image->gimp);

  return TRUE;
}

/*
 *  Indexed color conversion machinery
 */

static void
zero_histogram_gray (CFHistogram histogram)
{
  gint i;

  for (i = 0; i < 256; i++)
    histogram[i] = 0;
}


static void
zero_histogram_rgb (CFHistogram histogram)
{
  memset (histogram, 0,
          HIST_R_ELEMS * HIST_G_ELEMS * HIST_B_ELEMS * sizeof (ColorFreq));
}


static void
generate_histogram_gray (CFHistogram  histogram,
                         GimpLayer   *layer,
1112
                         gboolean     dither_alpha)
1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128
{
  GeglBufferIterator *iter;
  const Babl         *format;
  gint                bpp;
  gboolean            has_alpha;

  format = gimp_drawable_get_format (GIMP_DRAWABLE (layer));

  g_return_if_fail (format == babl_format ("Y' u8") ||
                    format == babl_format ("Y'A u8"));

  bpp       = babl_format_get_bytes_per_pixel (format);
  has_alpha = babl_format_has_alpha (format);

  iter = gegl_buffer_iterator_new (gimp_drawable_get_buffer (GIMP_DRAWABLE (layer)),
                                   NULL, 0, format,
1129
                                   GEGL_ACCESS_READ, GEGL_ABYSS_NONE, 1);
1130 1131 1132

  while (gegl_buffer_iterator_next (iter))
    {
1133
      const guchar *data   = iter->items[0].data;
1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162
      gint          length = iter->length;

      if (has_alpha)
        {
          while (length--)
            {
              if (data[ALPHA_G] > 127)
                histogram[*data]++;

              data += bpp;
            }
        }
      else
        {
          while (length--)
            {
              histogram[*data]++;

              data += bpp;
            }
        }
    }
}


static void
generate_histogram_rgb (CFHistogram   histogram,
                        GimpLayer    *layer,
                        gint          col_limit,
1163
                        gboolean      dither_alpha,
1164
                        GimpProgress *progress)
1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195
{
  GeglBufferIterator *iter;
  const Babl         *format;
  GeglRectangle      *roi;
  ColorFreq          *colfreq;
  gint                nfc_iter;
  gint                row, col, coledge;
  gint                offsetx, offsety;
  glong               layer_size;
  glong               total_size = 0;
  gint                count      = 0;
  gint                bpp;
  gboolean            has_alpha;

  format = gimp_drawable_get_format (GIMP_DRAWABLE (layer));

  g_return_if_fail (format == babl_format ("R'G'B' u8") ||
                    format == babl_format ("R'G'B'A u8"));

  bpp       = babl_format_get_bytes_per_pixel (format);
  has_alpha = babl_format_has_alpha (format);

  gimp_item_get_offset (GIMP_ITEM (layer), &offsetx, &offsety);

  layer_size = (gimp_item_get_width  (GIMP_ITEM (layer)) *
                gimp_item_get_height (GIMP_ITEM (layer)));

  /*  g_printerr ("col_limit = %d, nfc = %d\n", col_limit, num_found_cols); */

  iter = gegl_buffer_iterator_new (gimp_drawable_get_buffer (GIMP_DRAWABLE (layer)),
                                   NULL, 0, format,
1196 1197
                                   GEGL_ACCESS_READ, GEGL_ABYSS_NONE, 1);
  roi = &iter->items[0].roi;
1198 1199 1200 1201 1202 1203

  if (progress)
    gimp_progress_set_value (progress, 0.0);

  while (gegl_buffer_iterator_next (iter))
    {
1204
      const guchar *data   = iter->items[0].data;
1205 1206 1207 1208 1209 1210 1211 1212
      gint          length = iter->length;

      total_size += length;

      /* g_printerr (" [%d,%d - %d,%d]", srcPR.x, src_roi->y, offsetx, offsety); */

      if (needs_quantize)
        {
1213
          if (dither_alpha)
1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269 1270 1271 1272 1273 1274 1275 1276
            {
              /* if alpha-dithering,
                 we need to be deterministic w.r.t. offsets */

              col = roi->x + offsetx;
              coledge = col + roi->width;
              row = roi->y + offsety;

              while (length--)
                {
                  gboolean transparent = FALSE;

                  if (has_alpha &&
                      data[ALPHA] <
                      DM[col & DM_WIDTHMASK][row & DM_HEIGHTMASK])
                    transparent = TRUE;

                  if (! transparent)
                    {
                      colfreq = HIST_RGB (histogram,
                                          data[RED],
                                          data[GREEN],
                                          data[BLUE]);
                      (*colfreq)++;
                    }

                  col++;
                  if (col == coledge)
                    {
                      col = roi->x + offsetx;
                      row++;
                    }

                  data += bpp;
                }
            }
          else
            {
              while (length--)
                {
                  if ((has_alpha && ((data[ALPHA] > 127)))
                      || (!has_alpha))
                    {
                      colfreq = HIST_RGB (histogram,
                                          data[RED],
                                          data[GREEN],
                                          data[BLUE]);
                      (*colfreq)++;
                    }

                  data += bpp;
                }
            }
        }
      else
        {
          /* if alpha-dithering, we need to be deterministic w.r.t. offsets */
          col = roi->x + offsetx;
          coledge = col + roi->width;
          row = roi->y + offsety;

          while (length--)
            {
1277
              gboolean transparent = FALSE;
1278

1279 1280 1281 1282 1283
              if (has_alpha)
                {
                  if (dither_alpha)
                    {
                      if (data[ALPHA] <
1284
                          DM[col & DM_WIDTHMASK][row & DM_HEIGHTMASK])
1285 1286 1287 1288 1289 1290 1291 1292 1293 1294
                        transparent = TRUE;
                    }
                  else
                    {
                      if (data[ALPHA] <= 127)
                        transparent = TRUE;
                    }
                }

              if (! transparent)
1295 1296 1297 1298 1299 1300 1301 1302 1303 1304 1305 1306 1307 1308 1309 1310 1311 1312 1313 1314 1315 1316 1317 1318 1319 1320 1321 1322 1323 1324 1325 1326 1327 1328 1329 1330 1331 1332 1333 1334 1335 1336 1337 1338 1339 1340 1341 1342 1343 1344 1345 1346 1347 1348 1349 1350 1351 1352 1353 1354 1355
                {
                  colfreq = HIST_RGB (histogram,
                                      data[RED],
                                      data[GREEN],
                                      data[BLUE]);
                  (*colfreq)++;

                  if (!needs_quantize)
                    {
                      for (nfc_iter = 0;
                           nfc_iter < num_found_cols;
                           nfc_iter++)
                        {
                          if ((data[RED]   == found_cols[nfc_iter][0]) &&
                              (data[GREEN] == found_cols[nfc_iter][1]) &&
                              (data[BLUE]  == found_cols[nfc_iter][2]))
                            goto already_found;
                        }

                      /* Color was not in the table of
                       * existing colors
                       */

                      num_found_cols++;

                      if (num_found_cols > col_limit)
                        {
                          /* There are more colors in the image than
                           *  were allowed.  We switch to plain
                           *  histogram calculation with a view to
                           *  quantizing at a later stage.
                           */
                          needs_quantize = TRUE;
                          /* g_print ("\nmax colors exceeded - needs quantize.\n");*/
                          goto already_found;
                        }
                      else
                        {
                          /* Remember the new color we just found.
                           */
                          found_cols[num_found_cols-1][0] = data[RED];
                          found_cols[num_found_cols-1][1] = data[GREEN];
                          found_cols[num_found_cols-1][2] = data[BLUE];
                        }
                    }
                }
            already_found:

              col++;
              if (col == coledge)
                {
                  col = roi->x + offsetx;
                  row++;
                }

              data += bpp;
            }
        }

      if (progress && (count % 16 == 0))
        gimp_progress_set_value (progress,
1356
                                 (gdouble) total_size / (gdouble) layer_size);
1357 1358 1359 1360 1361 1362 1363 1364 1365 1366 1367 1368 1369 1370 1371 1372 1373 1374 1375 1376 1377 1378 1379 1380 1381 1382 1383 1384 1385 1386 1387 1388 1389 1390 1391 1392 1393 1394 1395 1396 1397 1398 1399 1400 1401 1402 1403 1404 1405 1406 1407 1408 1409 1410 1411 1412 1413 1414 1415 1416 1417 1418 1419 1420 1421 1422 1423 1424 1425 1426 1427 1428 1429 1430 1431 1432 1433 1434 1435 1436 1437 1438 1439 1440 1441 1442 1443 1444 1445 1446 1447 1448 1449 1450 1451 1452 1453 1454 1455 1456 1457 1458 1459 1460 1461 1462 1463 1464 1465 1466 1467 1468 1469 1470 1471 1472 1473 1474 1475 1476 1477 1478 1479 1480 1481 1482 1483 1484 1485 1486 1487 1488 1489 1490 1491 1492 1493 1494 1495 1496 1497 1498 1499 1500 1501 1502 1503 1504 1505 1506 1507 1508 1509 1510 1511 1512 1513 1514 1515 1516 1517 1518 1519 1520 1521 1522 1523 1524 1525 1526 1527 1528 1529 1530 1531 1532 1533 1534 1535 1536 1537 1538 1539 1540 1541 1542 1543 1544 1545 1546 1547 1548 1549 1550 1551 1552 1553 1554 1555 1556 1557 1558 1559 1560 1561 1562 1563 1564 1565 1566 1567 1568 1569 1570 1571 1572 1573 1574 1575 1576 1577 1578 1579 1580 1581 1582 1583 1584 1585 1586 1587 1588 1589 1590 1591 1592 1593 1594 1595 1596 1597 1598 1599 1600 1601 1602 1603 1604 1605 1606 1607 1608 1609 1610 1611 1612 1613 1614 1615 1616 1617 1618 1619 1620 1621 1622 1623 1624 1625 1626 1627 1628 1629 1630 1631 1632 1633 1634 1635 1636 1637 1638 1639 1640 1641 1642 1643 1644 1645 1646 1647 1648 1649 1650 1651 1652 1653 1654 1655 1656 1657 1658 1659 1660 1661 1662 1663 1664 1665 1666 1667 1668 1669 1670 1671 1672 1673 1674 1675 1676 1677 1678 1679 1680 1681 1682 1683 1684 1685 1686 1687 1688 1689 1690 1691 1692 1693 1694 1695 1696 1697 1698 1699 1700 1701 1702 1703 1704 1705 1706 1707 1708 1709 1710 1711 1712 1713 1714 1715 1716 1717 1718 1719 1720 1721 1722 1723 1724 1725 1726 1727 1728 1729 1730 1731 1732 1733 1734 1735 1736 1737 1738 1739 1740 1741 1742 1743 1744 1745 1746 1747 1748 1749 1750 1751 1752 1753 1754 1755 1756 1757 1758 1759 1760 1761 1762 1763 1764 1765 1766 1767 1768 1769 1770 1771 1772 1773 1774 1775 1776 1777 1778 1779 1780 1781 1782 1783 1784 1785 1786 1787 1788 1789 1790 1791 1792 1793 1794 1795 1796 1797 1798 1799 1800 1801 1802 1803 1804 1805 1806 1807 1808 1809 1810 1811 1812 1813 1814 1815 1816 1817 1818 1819 1820 1821 1822 1823 1824 1825 1826 1827 1828 1829 1830 1831 1832 1833 1834 1835 1836 1837 1838 1839 1840 1841 1842 1843 1844 1845 1846 1847 1848 1849 1850 1851 1852 1853 1854 1855 1856 1857 1858 1859 1860 1861 1862 1863 1864 1865 1866 1867 1868 1869 1870 1871 1872 1873 1874 1875 1876 1877 1878 1879 1880 1881 1882 1883 1884 1885 1886 1887 1888 1889 1890 1891 1892 1893 1894 1895 1896 1897 1898 1899 1900 1901 1902 1903 1904 1905 1906 1907 1908 1909 1910 1911 1912 1913 1914 1915 1916 1917 1918 1919 1920 1921 1922 1923 1924
    }

/*  g_print ("O: col_limit = %d, nfc = %d\n", col_limit, num_found_cols);*/
}



static boxptr
find_split_candidate (const boxptr  boxlist,
                      const gint    numboxes,
                      AxisType     *which_axis,
                      const gint    desired_colors)
{
  boxptr  boxp;
  gint    i;
  etype   maxc = 0;
  boxptr  which = NULL;
  gdouble Lbias;

  *which_axis = AXIS_UNDEF;

  /* we only perform the initial L-split bias /at all/ if the final
     number of desired colors is quite low, otherwise it all comes
     out in the wash anyway and this initial bias generally only hurts
     us in the long run. */
  if (desired_colors <= 16)
    {
#define BIAS_FACTOR 2.66F
#define BIAS_NUMBER 2 /* 0 */

      /* we bias towards splitting across L* for first few colors */
      Lbias = (numboxes > BIAS_NUMBER) ? 1.0F : ((gdouble) (BIAS_NUMBER + 1) -
                                                 ((gdouble) numboxes)) /
        ((gdouble) BIAS_NUMBER / BIAS_FACTOR);
      /*Lbias = 1.0;
        fprintf(stderr, " [[%d]] ", numboxes);
        fprintf(stderr, "Using ramped L-split bias.\n");
        fprintf(stderr, "R\n");
      */
    }
  else
    Lbias = 1.0F;

  for (i = 0, boxp = boxlist; i < numboxes; i++, boxp++)
    {
      if (boxp->volume > 0)
        {
#ifndef _MSC_VER
          etype rpe = (double)((boxp->rerror) * R_SCALE * R_SCALE);
          etype gpe = (double)((boxp->gerror) * G_SCALE * G_SCALE);
          etype bpe = (double)((boxp->berror) * B_SCALE * B_SCALE);
#else
          /*
           * Sorry about the mess, otherwise would get :
           * error C2520: conversion from unsigned __int64 to double
           *              not implemented, use signed __int64
           */
          etype rpe = (double)(((__int64)boxp->rerror) * R_SCALE * R_SCALE);
          etype gpe = (double)(((__int64)boxp->gerror) * G_SCALE * G_SCALE);
          etype bpe = (double)(((__int64)boxp->berror) * B_SCALE * B_SCALE);
#endif

          if (Lbias * rpe > maxc &&
              boxp->Rmin < boxp->Rmax)
            {
              which = boxp;
              maxc  = Lbias * rpe;
              *which_axis = AXIS_RED;
            }

          if (gpe > maxc &&
              boxp->Gmin < boxp->Gmax)
            {
              which = boxp;
              maxc  = gpe;
              *which_axis = AXIS_GREEN;
            }

          if (bpe > maxc &&
              boxp->Bmin < boxp->Bmax)
            {
              which = boxp;
              maxc  = bpe;
              *which_axis = AXIS_BLUE;
            }
        }
    }

  /*  fprintf(stderr, " %f,%p ", maxc, which); */
  /*  fprintf(stderr, " %llu ", maxc); */

  return which;
}


/* Find the splittable box with the largest (scaled) volume Returns
 * NULL if no splittable boxes remain
 */
static boxptr
find_biggest_volume (const boxptr boxlist,
                     const gint   numboxes)
{
  boxptr boxp;
  gint   i;
  gint   maxv = 0;
  boxptr which = NULL;

  for (i = 0, boxp = boxlist; i < numboxes; i++, boxp++)
    {
      if (boxp->volume > maxv)
        {
          which = boxp;
          maxv = boxp->volume;
        }
    }

  return which;
}


/* Shrink the min/max bounds of a box to enclose only nonzero
 * elements, and recompute its volume and population
 */
static void
update_box_gray (const CFHistogram histogram,
                 boxptr            boxp)
{
  gint      i, min, max, dist;
  ColorFreq ccount;

  min = boxp->Rmin;
  max = boxp->Rmax;

  if (max > min)
    for (i = min; i <= max; i++)
      {
        if (histogram[i] != 0)
          {
            boxp->Rmin = min = i;
            break;
          }
      }

  if (max > min)
    for (i = max; i >= min; i--)
      {
        if (histogram[i] != 0)
          {
            boxp->Rmax = max = i;
            break;
          }
      }

  /* 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.
   */
  dist = max - min;
  boxp->volume = dist * dist;

  /* Now scan remaining volume of box and compute population */
  ccount = 0;
  for (i = min; i <= max; i++)
    if (histogram[i] != 0)
      ccount++;

  boxp->colorcount = ccount;
}


/* Shrink the min/max bounds of a box to enclose only nonzero
 * elements, and recompute its volume, population and error
 */
static void
update_box_rgb (const CFHistogram histogram,
                boxptr            boxp,
                const gint        cells_remaining)
{
  gint      R, G, B;
  gint      Rmin, Rmax, Gmin, Gmax, Bmin, Bmax;
  gint      dist0, dist1, dist2;
  ColorFreq ccount;
  /*
  guint64 tempRerror;
  guint64 tempGerror;
  guint64 tempBerror;
  */
  QuantizeObj dummyqo;
  box         dummybox;

  /* fprintf(stderr, "U"); */

  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++)
        {
          for (B = Bmin; B <= Bmax; B++)
            {
              if (*HIST_LIN (histogram, R, G, B) != 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++)
        {
          for (B = Bmin; B <= Bmax; B++)
            {
              if (*HIST_LIN (histogram, R, G, B) != 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++)
        {
          for (B = Bmin; B <= Bmax; B++)
            {
              if (*HIST_LIN (histogram, R, G, B) != 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++)
        {
          for (B = Bmin; B <= Bmax; B++)
            {
              if (*HIST_LIN (histogram, R, G, B) != 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++)
        {
          for (G = Gmin; G <= Gmax; G++)
            {
              if (*HIST_LIN (histogram, R, G, B) != 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++)
        {
          for (G = Gmin; G <= Gmax; G++)
            {
              if (*HIST_LIN (histogram, R, G, B) != 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. (ADM: note: this isn't true.)
   * 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 = ((1 + Rmax - Rmin) << R_SHIFT) * R_SCALE;
  dist1 = ((1 + Gmax - Gmin) << G_SHIFT) * G_SCALE;
  dist2 = ((1 + Bmax - Bmin) << B_SHIFT) * B_SCALE;
  boxp->volume = dist0*dist0 + dist1*dist1 + dist2*dist2;
  /*  boxp->volume = dist0 * dist1 * dist2; */

  compute_color_lin8(&dummyqo, histogram, boxp, 0);

  /*printf("(%d %d %d)\n", dummyqo.cmap[0].red,dummyqo.cmap[0].green,dummyqo.cmap[0].blue);
    fflush(stdout);*/

  /* Now scan remaining volume of box and compute population */
  ccount = 0;
  boxp->error = 0;
  boxp->rerror = 0;
  boxp->gerror = 0;
  boxp->berror = 0;
  for (R = Rmin; R <= Rmax; R++)
    {
      for (G = Gmin; G <= Gmax; G++)
        {
          for (B = Bmin; B <= Bmax; B++)
            {
              ColorFreq freq_here;

              freq_here = *HIST_LIN (histogram, R, G, B);

              if (freq_here != 0)
                {
                  int ge, be, re;

                  dummybox.Rmin = dummybox.Rmax = R;
                  dummybox.Gmin = dummybox.Gmax = G;
                  dummybox.Bmin = dummybox.Bmax = B;
                  compute_color_lin8(&dummyqo, histogram, &dummybox, 1);

                  re = dummyqo.cmap[0].red   - dummyqo.cmap[1].red;
                  ge = dummyqo.cmap[0].green - dummyqo.cmap[1].green;
                  be = dummyqo.cmap[0].blue  - dummyqo.cmap[1].blue;

                  boxp->rerror += freq_here * (re) * (re);
                  boxp->gerror += freq_here * (ge) * (ge);
                  boxp->berror += freq_here * (be) * (be);

                  ccount += freq_here;
                }
            }
        }
    }

#if 0
  fg d;flg fd;kg fld;gflkfld
  /* Scan again, taking note of halfway error point for red axis */
  tempRerror = 0;
  boxp->Rhalferror = Rmin;
#warning r<=?
  for (R = Rmin; R <= Rmax; R++)
    {
      for (G = Gmin; G <= Gmax; G++)
        {
          for (B = Bmin; B <= Bmax; B++)
            {
              ColorFreq freq_here;
              freq_here = *HIST_LIN(histogram, R, G, B);
              if (freq_here != 0)
                {
                  int re;
                  int idist;
                  double dist;

                  dummybox.Rmin = dummybox.Rmax = R;
                  dummybox.Gmin = dummybox.Gmax = G;
                  dummybox.Bmin = dummybox.Bmax = B;
                  compute_color_lin8(&dummyqo, histogram, &dummybox, 1);

                  re = dummyqo.cmap[0].red   - dummyqo.cmap[1].red;

                  tempRerror += freq_here * (re) * (re);

                  if (tempRerror*2 >= boxp->rerror)
                    goto green_axisscan;
                  else
                    boxp->Rhalferror = R;
                }
            }
        }
    }
  fprintf(stderr, " D:");
 green_axisscan:

  fprintf(stderr, "<%d: %llu/%llu> ", R, tempRerror, boxp->rerror);
  /* Scan again, taking note of halfway error point for green axis */
  tempGerror = 0;
  boxp->Ghalferror = Gmin;
#warning G<=?
  for (G = Gmin; G <= Gmax; G++)
    {
      for (R = Rmin; R <= Rmax; R++)
        {
          for (B = Bmin; B <= Bmax; B++)
            {
              ColorFreq freq_here;
              freq_here = *HIST_LIN(histogram, R, G, B);
              if (freq_here != 0)
                {
                  int ge;
                  dummybox.Rmin = dummybox.Rmax = R;
                  dummybox.Gmin = dummybox.Gmax = G;
                  dummybox.Bmin = dummybox.Bmax = B;
                  compute_color_lin8(&dummyqo, histogram, &dummybox, 1);

                  ge = dummyqo.cmap[0].green - dummyqo.cmap[1].green;

                  tempGerror += freq_here * (ge) * (ge);

                  if (tempGerror*2 >= boxp->gerror)
                    goto blue_axisscan;
                  else
                    boxp->Ghalferror = G;
                }
            }
        }
    }

 blue_axisscan:
  /* Scan again, taking note of halfway error point for blue axis */
  tempBerror = 0;
  boxp->Bhalferror = Bmin;
#warning B<=?
  for (B = Bmin; B <= Bmax; B++)
    {
      for (R = Rmin; R <= Rmax; R++)
        {
          for (G = Gmin; G <= Gmax; G++)
            {
              ColorFreq freq_here;
              freq_here = *HIST_LIN(histogram, R, G, B);
              if (freq_here != 0)
                {
                  int be;
                  dummybox.Rmin = dummybox.Rmax = R;
                  dummybox.Gmin = dummybox.Gmax = G;
                  dummybox.Bmin = dummybox.Bmax = B;
                  compute_color_lin8(&dummyqo, histogram, &dummybox, 1);

                  be = dummyqo.cmap[0].blue  - dummyqo.cmap[1].blue;

                  tempBerror += freq_here * (be) * (be);

                  if (tempBerror*2 >= boxp->berror)
                    goto finished_axesscan;
                  else
                    boxp->Bhalferror = B;
                }
            }
        }
    }
 finished_axesscan:
#else

  boxp->Rhalferror = Rmin + (Rmax - Rmin + 1) / 2;
  boxp->Ghalferror = Gmin + (Gmax - Gmin + 1) / 2;
  boxp->Bhalferror = Bmin + (Bmax - Bmin + 1) / 2;

  if (dist0 && dist1 && dist2)
    {
      AxisType longest_ax      = AXIS_UNDEF;
      gint     longest_length  = 0;
      gint     longest_length2 = 0;
      gint     ratio;

      /*
        fprintf(stderr, "[%d,%d,%d=%d,%d,%d] ",
        (Rmax - Rmin), (Gmax - Gmin), (Bmax - Bmin),
        dist0, dist1, dist2);
      */

      if (dist0 >= longest_length)
        {
          longest_length2 = longest_length;
          longest_length = dist0;
          longest_ax = AXIS_RED;
        }
      else if (dist0 >= longest_length2)
        {
          longest_length2 = dist0;
        }

      if (dist1 >= longest_length)
        {
          longest_length2 = longest_length;
          longest_length = dist1;
          longest_ax = AXIS_GREEN;
        }
      else if (dist1 >= longest_length2)
        {
          longest_length2 = dist1;
        }

      if (dist2 >= longest_length)
        {
          longest_length2 = longest_length;
          longest_length = dist2;
          longest_ax = AXIS_BLUE;
        }
      else if (dist2 >= longest_length2)
        {
          longest_length2 = dist2;
        }

      if (longest_length2 == 0)
        longest_length2 = 1;

      ratio = (longest_length + longest_length2/2) / longest_length2;
      /* fprintf(stderr, " ratio:(%d/%d)=%d ", longest_length, longest_length2, ratio);
         fprintf(stderr, "C%d ", cells_remaining); */

      if (ratio > cells_remaining + 1)
        ratio = cells_remaining + 1;

      if (ratio > 2)
        {
          switch (longest_ax)
            {
            case AXIS_RED:
              if (Rmin + (Rmax - Rmin + ratio / 2) / ratio < Rmax)
                {
                  /* fprintf(stderr, "FR%d \007\n",ratio);*/
                  boxp->Rhalferror = Rmin +  (Rmax - Rmin + ratio / 2) / ratio;
                }
              break;
            case AXIS_GREEN:
              if (Gmin + (Gmax - Gmin + ratio / 2) / ratio < Gmax)
                {
                  /* fprintf(stderr, "FG%d \007\n",ratio);*/
                  boxp->Ghalferror = Gmin + (Gmax - Gmin + ratio / 2) / ratio;
                }
              break;
            case AXIS_BLUE:
              if (Bmin + (Bmax - Bmin + ratio / 2) / ratio < Bmax)
                {
                  /* fprintf(stderr, "FB%d \007\n",ratio);*/
                  boxp->Bhalferror = Bmin + (Bmax - Bmin + ratio / 2) / ratio;
                }
              break;
            default:
              g_warning ("GRR, UNDEF LONGEST AXIS\007\n");
            }
        }
    }

  if (boxp->Rhalferror == Rmax)
    boxp->Rhalferror = Rmin;
  if (boxp->Ghalferror == Gmax)
    boxp->Ghalferror = Gmin;
  if (boxp->Bhalferror == Bmax)
    boxp->Bhalferror = Bmin;

  /*
  boxp->Rhalferror = RSDF(dummyqo.cmap[0].red);
  boxp->Ghalferror = GSDF(dummyqo.cmap[0].green);
  boxp->Bhalferror = BSDF(dummyqo.cmap[0].blue);
  */

  /*
  boxp->Rhalferror = (RSDF(dummyqo.cmap[0].red) + (Rmin+Rmax)/2)/2;
  boxp->Ghalferror = (GSDF(dummyqo.cmap[0].green) + (Gmin+Gmax)/2)/2;
  boxp->Bhalferror = (BSDF(dummyqo.cmap[0].blue) + (Bmin+Bmax)/2)/2;
  */


#endif
  /*
  fprintf(stderr, " %d,%d", dummyqo.cmap[0].blue, boxp->Bmax);

1925 1926 1927 1928 1929 1930
  gimp_assert (boxp->Rhalferror >= boxp->Rmin);
  gimp_assert (boxp->Rhalferror < boxp->Rmax);
  gimp_assert (boxp->Ghalferror >= boxp->Gmin);
  gimp_assert (boxp->Ghalferror < boxp->Gmax);
  gimp_assert (boxp->Bhalferror >= boxp->Bmin);
  gimp_assert (boxp->Bhalferror < boxp->Bmax);*/
1931 1932 1933 1934 1935 1936 1937 1938 1939 1940 1941 1942 1943 1944 1945 1946 1947 1948 1949 1950 1951 1952 1953 1954 1955 1956 1957 1958 1959 1960 1961 1962 1963 1964 1965 1966 1967 1968 1969 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022 2023 2024 2025 2026 2027 2028 2029

  /*boxp->error = (sqrt((double)(boxp->error/ccount)));*/
  /*  boxp->rerror = (sqrt((double)((boxp->rerror)/ccount)));
  boxp->gerror = (sqrt((double)((boxp->gerror)/ccount)));
  boxp->berror = (sqrt((double)((boxp->berror)/ccount)));*/
  /*printf(":%lld / %ld: ", boxp->error, ccount);
  printf("(%d-%d-%d)(%d-%d-%d)(%d-%d-%d)\n",
         Rmin, boxp->Rhalferror, Rmax,
         Gmin, boxp->Ghalferror, Gmax,
         Bmin, boxp->Bhalferror, Bmax
         );
         fflush(stdout);*/

  boxp->colorcount = ccount;
}


/* Repeatedly select and split the largest box until we have enough
 * boxes
 */
static gint
median_cut_gray (CFHistogram histogram,
                 boxptr      boxlist,
                 gint        numboxes,
                 gint        desired_colors)
{
  gint   lb;
  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 == NULL)           /* 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->Rmin = b1->Rmin;

      /* 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.
       */
      lb = (b1->Rmax + b1->Rmin) / 2;
      b1->Rmax = lb;
      b2->Rmin = lb + 1;

      /* Update stats for boxes */
      update_box_gray (histogram, b1);
      update_box_gray (histogram, b2);
      numboxes++;
    }

  return numboxes;
}

/* Repeatedly select and split the largest box until we have enough
 * boxes
 */
static gint
median_cut_rgb (CFHistogram   histogram,
                boxptr        boxlist,
                gint          numboxes,
                gint          desired_colors,
                GimpProgress *progress)
{
  gint     lb;
  boxptr   b1, b2;
  AxisType which_axis;

  while (numboxes < desired_colors)
    {
      b1 = find_split_candidate (boxlist, numboxes, &which_axis, desired_colors);

      if (b1 == NULL)           /* 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 split point along selected axis, and update box bounds.
       * Note that lb value is max for lower box, so must be < old max.
       */
      switch (which_axis)
        {
        case AXIS_RED:
          lb = b1->Rhalferror;/* *0 + (b1->Rmax + b1->Rmin) / 2; */
          b1->Rmax = lb;
          b2->Rmin = lb+1;
2030 2031
          g_return_val_if_fail (b1->Rmax >= b1->Rmin, numboxes);
          g_return_val_if_fail (b2->Rmax >= b2->Rmin, numboxes);
2032 2033 2034 2035 2036
          break;
        case AXIS_GREEN:
          lb = b1->Ghalferror;/* *0 + (b1->Gmax + b1->Gmin) / 2; */
          b1->Gmax = lb;
          b2->Gmin = lb+1;
2037 2038
          g_return_val_if_fail (b1->Gmax >= b1->Gmin, numboxes);
          g_return_val_if_fail (b2->Gmax >= b2->Gmin, numboxes);
2039 2040 2041 2042 2043
          break;
        case AXIS_BLUE:
          lb = b1->Bhalferror;/* *0 + (b1->Bmax + b1->Bmin) / 2; */
          b1->Bmax = lb;
          b2->Bmin = lb+1;
2044 2045
          g_return_val_if_fail (b1->Bmax >= b1->Bmin, numboxes);
          g_return_val_if_fail (b2->Bmax >= b2->Bmin, numboxes);
2046 2047 2048 2049 2050 2051 2052 2053 2054 2055 2056 2057 2058 2059 2060 2061 2062 2063 2064 2065 2066 2067 2068 2069 2070 2071 2072 2073 2074 2075 2076 2077 2078 2079 2080 2081 2082 2083 2084 2085 2086 2087 2088 2089 2090 2091 2092 2093 2094 2095 2096 2097 2098 2099 2100 2101 2102 2103 2104 2105 2106 2107 2108 2109 2110 2111 2112 2113 2114 2115 2116 2117 2118 2119 2120 2121 2122 2123 2124 2125 2126 2127 2128 2129 2130 2131 2132 2133 2134 2135 2136 2137 2138 2139 2140 2141 2142 2143 2144 2145 2146 2147 2148 2149 2150 2151 2152 2153 2154 2155 2156 2157 2158 2159 2160 2161 2162 2163 2164 2165 2166 2167 2168 2169 2170 2171 2172 2173 2174 2175 2176 2177 2178 2179 2180 2181 2182 2183 2184 2185 2186 2187 2188 2189 2190 2191 2192 2193 2194 2195 2196 2197 2198 2199 2200 2201 2202 2203 2204 2205 2206 2207 2208 2209 2210 2211 2212 2213 2214 2215 2216 2217 2218 2219 2220 2221 2222 2223 2224 2225 2226 2227 2228 2229 2230 2231 2232 2233 2234 2235 2236 2237 2238 2239 2240 2241 2242 2243 2244 2245 2246 2247 2248 2249 2250 2251 2252 2253 2254 2255 2256 2257 2258 2259 2260 2261 2262 2263 2264 2265 2266 2267 2268 2269 2270 2271 2272 2273 2274 2275 2276 2277 2278 2279 2280 2281 2282 2283 2284 2285 2286 2287 2288 2289 2290 2291 2292 2293 2294 2295 2296 2297 2298 2299 2300 2301 2302 2303 2304 2305 2306 2307 2308 2309 2310 2311 2312 2313 2314 2315 2316 2317 2318 2319 2320 2321 2322 2323 2324 2325 2326 2327 2328 2329 2330 2331 2332 2333 2334 2335 2336 2337 2338 2339 2340 2341 2342 2343 2344 2345 2346 2347 2348 2349 2350 2351 2352 2353 2354 2355 2356 2357 2358 2359 2360 2361 2362 2363 2364 2365 2366 2367 2368 2369 2370 2371 2372 2373 2374 2375 2376 2377 2378 2379 2380 2381 2382 2383 2384 2385 2386 2387 2388 2389 2390 2391 2392 2393 2394 2395 2396 2397 2398 2399 2400 2401 2402 2403 2404 2405 2406 2407 2408 2409 2410 2411 2412 2413 2414 2415 2416 2417 2418 2419 2420 2421 2422 2423 2424 2425 2426 2427 2428 2429 2430 2431 2432 2433 2434 2435 2436 2437 2438 2439 2440 2441 2442 2443 2444 2445 2446 2447 2448 2449 2450 2451 2452 2453 2454 2455 2456 2457 2458 2459 2460 2461 2462 2463 2464 2465 2466 2467 2468 2469 2470 2471 2472 2473 2474 2475 2476 2477 2478 2479 2480 2481 2482 2483 2484 2485 2486 2487 2488 2489 2490 2491 2492 2493 2494 2495 2496 2497 2498 2499 2500 2501 2502 2503 2504 2505 2506 2507 2508 2509 2510 2511 2512 2513 2514 2515 2516 2517 2518 2519 2520 2521 2522 2523 2524 2525 2526 2527 2528 2529 2530 2531 2532 2533 2534 2535 2536 2537 2538 2539 2540 2541 2542 2543 2544 2545 2546 2547 2548 2549 2550 2551 2552 2553 2554 2555 2556 2557 2558 2559 2560 2561 2562 2563 2564 2565 2566 2567 2568 2569 2570 2571 2572 2573 2574 2575 2576 2577 2578 2579 2580 2581 2582 2583 2584 2585 2586 2587 2588 2589 2590 2591 2592 2593 2594 2595 2596 2597 2598 2599 2600 2601 2602 2603 2604 2605 2606 2607 2608 2609 2610 2611 2612 2613 2614 2615 2616 2617 2618 2619 2620 2621 2622 2623 2624 2625 2626 2627 2628 2629 2630 2631 2632 2633 2634 2635 2636 2637 2638 2639 2640 2641 2642 2643 2644 2645 2646 2647 2648 2649 2650 2651 2652 2653 2654 2655 2656 2657 2658 2659 2660 2661 2662 2663 2664 2665 2666 2667 2668 2669 2670 2671 2672 2673 2674 2675 2676 2677 2678 2679 2680 2681 2682 2683 2684 2685 2686 2687 2688 2689 2690 2691 2692 2693 2694 2695 2696 2697 2698 2699 2700 2701 2702 2703 2704 2705 2706 2707 2708 2709 2710 2711 2712 2713 2714 2715 2716 2717 2718 2719 2720 2721 2722 2723 2724 2725 2726 2727 2728 2729 2730 2731 2732 2733 2734 2735 2736 2737 2738 2739 2740 2741 2742 2743 2744 2745 2746 2747 2748 2749 2750 2751 2752 2753 2754 2755 2756 2757 2758 2759 2760 2761 2762 2763 2764 2765 2766 2767 2768 2769 2770 2771 2772 2773 2774 2775 2776 2777 2778 2779 2780 2781 2782 2783 2784 2785 2786 2787 2788 2789 2790 2791 2792 2793 2794 2795 2796 2797 2798 2799 2800 2801 2802 2803 2804 2805 2806 2807 2808 2809 2810 2811 2812 2813 2814 2815 2816 2817 2818 2819 2820 2821 2822 2823 2824 2825 2826 2827 2828 2829 2830 2831 2832 2833 2834
          break;
        default:
          g_error ("Uh-oh.");
        }
      /* Update stats for boxes */
      numboxes++;

      if (progress && (numboxes % 16 == 0))
        gimp_progress_set_value (progress, (gdouble) numboxes / desired_colors);

      update_box_rgb (histogram, b1, desired_colors - numboxes);
      update_box_rgb (histogram, b2, desired_colors - numboxes);
    }

  return numboxes;
}


/* Compute representative color for a box, put it in colormap[icolor]
 */
static void
compute_color_gray (QuantizeObj *quantobj,
                    CFHistogram  histogram,
                    boxptr       boxp,
                    int          icolor)
{
  gint    i, min, max;
  guint64 count;
  guint64 total;
  guint64 gtotal;

  min = boxp->Rmin;
  max = boxp->Rmax;

  total = 0;
  gtotal = 0;

  for (i = min; i <= max; i++)
    {
      count = histogram[i];
      if (count != 0)
        {
          total += count;
          gtotal += i * count;
        }
    }

  if (total != 0)
    {
      quantobj->cmap[icolor].red   =
      quantobj->cmap[icolor].green =
      quantobj->cmap[icolor].blue  = (gtotal + (total >> 1)) / total;
    }
  else
    {
      /* The only situation where total==0 is if the image was null or
       *  all-transparent.  In that case we just put a dummy value in
       *  the colormap.
       */
      quantobj->cmap[icolor].red   =
      quantobj->cmap[icolor].green =
      quantobj->cmap[icolor].blue  = 0;
    }
}


/* Compute representative color for a box, put it in colormap[icolor]
 */
static void
compute_color_rgb (QuantizeObj *quantobj,
                   CFHistogram  histogram,
                   boxptr       boxp,
                   int          icolor)
{
  /* Current algorithm: mean weighted by pixels (not colors) */
  /* Note it is important to get the rounding correct! */
  gint      R, G, B;
  gint      Rmin, Rmax;
  gint      Gmin, Gmax;
  gint      Bmin, Bmax;
  ColorFreq total  = 0;
  ColorFreq Rtotal = 0;
  ColorFreq Gtotal = 0;
  ColorFreq 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++)
      {
        for (B = Bmin; B <= Bmax; B++)
          {
            ColorFreq this_freq = *HIST_LIN (histogram, R, G, B);

            if (this_freq != 0)
              {
                total += this_freq;
                Rtotal += R * this_freq;
                Gtotal += G * this_freq;
                Btotal += B * this_freq;
              }
          }
      }

  if (total > 0)
    {
      guchar red, green, blue;

      lin_to_rgb (/*(Rtotal + (total>>1)) / total,
                    (Gtotal + (total>>1)) / total,
                    (Btotal + (total>>1)) / total,*/
                  (double)Rtotal / (double)total,
                  (double)Gtotal / (double)total,
                  (double)Btotal / (double)total,
                  &red, &green, &blue);

      quantobj->cmap[icolor].red   = red;
      quantobj->cmap[icolor].green = green;
      quantobj->cmap[icolor].blue  = blue;
    }
  else
    {
      /* The only situation where total==0 is if the image was null or
       *  all-transparent.  In that case we just put a dummy value in
       *  the colormap.
       */
      quantobj->cmap[icolor].red   = 0;
      quantobj->cmap[icolor].green = 0;
      quantobj->cmap[icolor].blue  = 0;
    }
}


/* Compute representative color for a box, put it in colormap[icolor]
 */
static void
compute_color_lin8 (QuantizeObj *quantobj,
                    CFHistogram  histogram,
                    boxptr       boxp,
                    const gint   icolor)
{
  /* Current algorithm: mean weighted by pixels (not colors) */
  /* Note it is important to get the rounding correct! */
  gint      R, G, B;
  gint      Rmin, Rmax;
  gint      Gmin, Gmax;
  gint      Bmin, Bmax;
  ColorFreq total  = 0;
  ColorFreq Rtotal = 0;
  ColorFreq Gtotal = 0;
  ColorFreq 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++)
      {
        for (B = Bmin; B <= Bmax; B++)
          {
            ColorFreq this_freq = *HIST_LIN (histogram, R, G, B);

            if (this_freq != 0)
              {
                Rtotal += R * this_freq;
                Gtotal += G * this_freq;
                Btotal += B * this_freq;
                total += this_freq;
              }
          }
      }

  if (total != 0)
    {
      quantobj->cmap[icolor].red   = ((Rtotal << R_SHIFT) + (total>>1)) / total;
      quantobj->cmap[icolor].green = ((Gtotal << G_SHIFT) + (total>>1)) / total;
      quantobj->cmap[icolor].blue  = ((Btotal << B_SHIFT) + (total>>1)) / total;
    }
  else
    {
      /* The only situation where total==0 is if the image was null or
       *  all-transparent.  In that case we just put a dummy value in
       *  the colormap.
       */
      g_warning ("eep.");
      quantobj->cmap[icolor].red   = 0;
      quantobj->cmap[icolor].green = 128;
      quantobj->cmap[icolor].blue  = 128;
    }
}


/* Master routine for color selection
 */
static void
select_colors_gray (QuantizeObj *quantobj,
                    CFHistogram  histogram)
{
  boxptr boxlist;
  gint   numboxes;
  gint   desired = quantobj->desired_number_of_colors;
  gint   i;

  /* Allocate workspace for box list */
  boxlist = g_new (box, desired);

  /* Initialize one box containing whole space */
  numboxes = 1;
  boxlist[0].Rmin = 0;
  boxlist[0].Rmax = 255;
  /* Shrink it to actually-used volume and set its statistics */
  update_box_gray (histogram, boxlist);
  /* Perform median-cut to produce final box list */
  numboxes = median_cut_gray (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_gray (quantobj, histogram, boxlist + i, i);
}


/* Master routine for color selection
 */
static void
select_colors_rgb (QuantizeObj *quantobj,
                   CFHistogram  histogram)
{
  boxptr boxlist;
  gint   numboxes;
  gint   desired = quantobj->desired_number_of_colors;
  gint  i;

  /* Allocate workspace for box list */
  boxlist = g_new (box, desired);

  /* Initialize one box containing whole space */
  numboxes = 1;
  boxlist[0].Rmin = 0;
  boxlist[0].Rmax = HIST_R_ELEMS - 1;
  boxlist[0].Gmin = 0;
  boxlist[0].Gmax = HIST_G_ELEMS - 1;
  boxlist[0].Bmin = 0;
  boxlist[0].Bmax = HIST_B_ELEMS - 1;
  /* Shrink it to actually-used volume and set its statistics */
  update_box_rgb (histogram, &boxlist[0], quantobj->desired_number_of_colors);
  /* Perform median-cut to produce final box list */
  numboxes = median_cut_rgb (histogram, boxlist, numboxes, desired,
                             quantobj->progress);

  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);
    }

  g_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.
 *
 * 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)*/

/*adam*/
#define BOX_R_LOG 0
#define BOX_G_LOG 0
#define BOX_B_LOG 0

#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.
 */

/* 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.
 */
static gint
find_nearby_colors (QuantizeObj *quantobj,
                    int          minR,
                    int          minG,
                    int          minB,
                    int          colorlist[])
{
  int numcolors = quantobj->actual_number_of_colors;
  int maxR, maxG, maxB;
  int centerR, centerG, centerB;
  int i, x, ncolors;
  int minmaxdist, min_dist, 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) >> 1;
  maxG = minG + ((1 << BOX_G_SHIFT) - (1 << G_SHIFT));
  centerG = (minG + maxG + 1) >> 1;
  maxB = minB + ((1 << BOX_B_SHIFT) - (1 << B_SHIFT));
  centerB = (minB + maxB + 1) >> 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->clin[i].red;
      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->clin[i].green;
      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->clin[i].blue;
      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;
}


/* 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.
 */
static void
find_best_colors (QuantizeObj *quantobj,
                  gint         minR,
                  gint         minG,
                  gint         minB,
                  gint         numcolors,
                  gint         colorlist[],
                  gint         bestcolor[])
{
  gint  iR, iG, iB;
  gint  i, icolor;
  gint *bptr;           /* pointer into bestdist[] array */
  gint *cptr;           /* pointer into bestcolor[] array */
  gint  dist0, dist1;     /* initial distance values */
  gint  dist2;            /* current distance in inner loop */
  gint  xx0, xx1;         /* distance increments */
  gint  xx2;
  gint  inR, inG, inB;    /* initial values for increments */

  /* This array holds the distance to the nearest-so-far color for each cell */
  gint  bestdist[BOX_R_ELEMS * BOX_G_ELEMS * BOX_B_ELEMS] = { 0, };

  /* 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->clin[icolor].red) * R_SCALE;
      dist0 = inR*inR;
      /* special-case for L*==0: chroma diffs irrelevant */
      /*    if (minR > 0 || quantobj->clin[icolor].red > 0) */
      {
        inG = (minG - quantobj->clin[icolor].green) * G_SCALE;
        dist0 += inG*inG;
        inB = (minB - quantobj->clin[icolor].blue) * B_SCALE;
        dist0 += inB*inB;
      }
      /*    else
            {
            inG = 0;
            inB = 0;
            } */
      /* 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;
        }
    }
}


/* 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.)
 */
static void
fill_inverse_cmap_gray (QuantizeObj *quantobj,
                        CFHistogram  histogram,
                        gint         pixel)
{
  Color *cmap = quantobj->cmap;
  glong  mindist;
  gint   mindisti;
  gint   i;

  mindist  = 65536;
  mindisti = -1;

  for (i = 0; i < quantobj->actual_number_of_colors; i++)
    {
      glong dist = ABS (pixel - cmap[i].red);

      if (dist < mindist)
        {
          mindist  = dist;
          mindisti = i;
        }
    }

  if (i >= 0)
    histogram[pixel] = mindisti + 1;
}


/* 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.)
 */
static void
fill_inverse_cmap_rgb (QuantizeObj *quantobj,
                       CFHistogram  histogram,
                       gint         R,
                       gint         G,
                       gint         B)
{
  gint  minR, minG, minB; /* lower left corner of update box */
  gint  iR, iG, iB;
  gint *cptr;           /* pointer into bestcolor[] array */
  /* This array lists the candidate colormap indexes. */
  gint  colorlist[MAXNUMCOLORS];
  gint  numcolors;                /* number of candidate colors */
  /* This array holds the actually closest colormap index for each cell. */
  gint  bestcolor[BOX_R_ELEMS * BOX_G_ELEMS * BOX_B_ELEMS] = { 0, };

  /* 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++)
        {
          for (iB = 0; iB < BOX_B_ELEMS; iB++)
            {
              *HIST_LIN (histogram, R + iR, G + iG, B + iB) = (*cptr++) + 1;
            }
        }
    }
}


/*  This is pass 1  */

static void
median_cut_pass1_gray (QuantizeObj *quantobj)
{
  select_colors_gray (quantobj, quantobj->histogram);
}


static void
median_cut_pass1_rgb (QuantizeObj *quantobj)
{
  select_colors_rgb (quantobj, quantobj->histogram);
}


static void
monopal_pass1 (QuantizeObj *quantobj)
{