colorquant2.c

Go to the documentation of this file.

/*====================================================================*
 -  Copyright (C) 2001 Leptonica.  All rights reserved.
 -  This software is distributed in the hope that it will be
 -  useful, but with NO WARRANTY OF ANY KIND.
 -  No author or distributor accepts responsibility to anyone for the
 -  consequences of using this software, or for whether it serves any
 -  particular purpose or works at all, unless he or she says so in
 -  writing.  Everyone is granted permission to copy, modify and
 -  redistribute this source code, for commercial or non-commercial
 -  purposes, with the following restrictions: (1) the origin of this
 -  source code must not be misrepresented; (2) modified versions must
 -  be plainly marked as such; and (3) this notice may not be removed
 -  or altered from any source or modified source distribution.
 *====================================================================*/

/*
 *  colorquant2.c
 *                     
 *  Modified median cut color quantization
 *
 *      High level
 *          PIX              *pixMedianCutQuant()
 *          PIX              *pixMedianCutQuantGeneral()
 *          PIX              *pixMedianCutQuantMixed()
 *          PIX              *pixFewColorsMedianCutQuantMixed()
 *
 *      Median cut indexed histogram
 *          l_int32          *pixMedianCutHisto()
 *
 *      Static helpers
 *          static PIXCMAP   *pixcmapGenerateFromHisto()
 *          static PIX       *pixQuantizeWithColormap()
 *          static void       getColorIndexMedianCut()
 *          static L_BOX3D   *pixGetColorRegion()
 *          static l_int32    medianCutApply()
 *          static PIXCMAP   *pixcmapGenerateFromMedianCuts()
 *          static l_int32    vboxGetAverageColor()
 *          static l_int32    vboxGetCount()
 *          static l_int32    vboxGetVolume()
 *          static L_BOX3D   *box3dCreate();
 *          static L_BOX3D   *box3dCopy();
 *
 *   Paul Heckbert published the median cut algorithm, "Color Image
 *   Quantization for Frame Buffer Display," in Proc. SIGGRAPH '82,
 *   Boston, July 1982, pp. 297-307.  A copy of the paper without
 *   figures can be found on the web.
 *
 *   Median cut starts with either the full color space or the occupied
 *   region of color space.  If you're not dithering, the occupied region
 *   can be used, but with dithering, pixels can end up in any place
 *   in the color space, so you must represent the entire color space in
 *   the final colormap.
 *
 *   Color components are quantized to typically 5 or 6 significant
 *   bits (for each of r, g and b).   Call a 3D region of color
 *   space a 'vbox'.  Any color in this quantized space is represented
 *   by an element of a linear histogram array, indexed by rgb value.
 *   The initial region is then divided into two regions that have roughly
 *   equal pixel occupancy (hence the name "median cut").  Subdivision
 *   continues until the requisite number of vboxes has been generated.
 *
 *   But the devil is in the details of the subdivision process.
 *   Here are some choices that you must make:
 *     (1) Along which axis to subdivide?
 *     (2) Which box to put the bin with the median pixel?
 *     (3) How to order the boxes for subdivision?
 *     (4) How to adequately handle boxes with very small numbers of pixels?
 *     (5) How to prevent a little-represented but highly visible color
 *         from being masked out by other colors in its vbox.
 *
 *   Taking these in order:
 *     (1) Heckbert suggests using either the largest vbox side, or the vbox
 *         side with the largest variance in pixel occupancy.  We choose
 *         to divide based on the largest vbox side.
 *     (2) Suppose you've chosen a side.  Then you have a histogram
 *         of pixel occupancy in 2D slices of the vbox.  One of those
 *         slices includes the median pixel.  Suppose there are L bins
 *         to the left (smaller index) and R bins to the right.  Then
 *         this slice (or bin) should be assigned to the box containing
 *         the smaller of L and R.  This both shortens the larger
 *         of the subdivided dimensions and helps a low-count color
 *         far from the subdivision boundary to better express itself.
 *     (2a) One can also ask if the boundary should be moved even
 *         farther into the longer side.  This is feasable if we have
 *         a method for doing extra subdivisions on the high count
 *         vboxes.  And we do (see (3)).
 *     (3) To make sure that the boxes are subdivided toward equal
 *         occupancy, use an occupancy-sorted priority queue, rather
 *         than a simple queue.
 *     (4) With a priority queue, boxes with small number of pixels
 *         won't be repeatedly subdivided.  This is good.
 *     (5) Use of a priority queue allows tricks such as in (2a) to let
 *         small occupancy clusters be better expressed.  In addition,
 *         rather than splitting near the median, small occupancy colors
 *         are best reproduced by cutting half-way into the longer side.
 *
 *   However, serious problems can arise with dithering if a priority
 *   queue is used based on population alone.  If the picture has
 *   large regions of nearly constant color, some vboxes can be very
 *   large and have a sizeable population (but not big enough to get to
 *   the head of the queue).  If one of these large, occupied vboxes
 *   is near in color to a nearly constant color region of the
 *   image, dithering can inject pixels from the large vbox into
 *   the nearly uniform region.  These pixels can be very far away
 *   in color, and the oscillations are highly visible.  To prevent
 *   this, we can take either or both of these actions:
 *
 *     (1) Subdivide a fraction (< 1.0) based on population, and
 *         do the rest of the subdivision based on the product of
 *         the vbox volume and its population.  By using the product,
 *         we avoid further subdivision of nearly empty vboxes, and
 *         directly target large vboxes with significant population.
 *
 *     (2) Threshold the excess color transferred in dithering to
 *         neighboring pixels.
 *
 *   Doing either of these will stop the most annoying oscillations
 *   in dithering.  Furthermore, by doing (1), we also improve the
 *   rendering of regions of nearly constant color, both with and
 *   without dithering.  It turns out that the image quality is
 *   not sensitive to the value of the parameter in (1); values
 *   between 0.3 and 0.9 give very good results.
 *
 *   Here's the lesson: subdivide the color space into vboxes such
 *   that (1) the most populated vboxes that can be further
 *   subdivided (i.e., that occupy more than one quantum volume
 *   in color space) all have approximately the same population,
 *   and (2) all large vboxes have no significant population.
 *   If these conditions are met, the quantization will be excellent.
 *
 *   Once the subdivision has been made, the colormap is generated,
 *   with one color for each vbox and using the average color in the vbox.
 *   At the same time, the histogram array is converted to an inverse
 *   colormap table, storing the colormap index in every cell in the
 *   vbox.  Finally, using both the colormap and the inverse colormap,
 *   a colormapped pix is quickly generated from the original rgb pix.
 *
 *   In the present implementation, subdivided regions of colorspace
 *   that are not occupied are retained, but not further subdivided.
 *   This is required for our inverse colormap lookup table for
 *   dithering, because dithered pixels may fall into these unoccupied
 *   regions.  For such empty regions, we use the center as the rgb
 *   colormap value.
 *
 *   This variation on median cut can be referred to as "Modified Median
 *   Cut" quantization, or MMCQ.  Overall, the undithered MMCQ gives
 *   comparable results to the two-pass Octcube Quantizer (OQ).
 *   Comparing the two methods on the test24.jpg painting, we see:
 *
 *     (1) For rendering spot color (the various reds and pinks in
 *         the image), MMCQ is not as good as OQ.
 *
 *     (2) For rendering majority color regions, MMCQ does a better
 *         job of avoiding posterization.  That is, it does better
 *         dividing the color space up in the most heavily populated regions.
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "allheaders.h"

    /* Median cut 3-d volume element.  Sort on first element, which
     * can be the number of pixels, the volume or a combination
     * of these.   */
struct L_Box3d
{
    l_float32        sortparam;  /* parameter on which to sort the vbox */
    l_int32          npix;       /* number of pixels in the vbox        */
    l_int32          vol;        /* quantized volume of vbox            */
    l_int32          r1;         /* min r index in the vbox             */
    l_int32          r2;         /* max r index in the vbox             */
    l_int32          g1;         /* min g index in the vbox             */
    l_int32          g2;         /* max g index in the vbox             */
    l_int32          b1;         /* min b index in the vbox             */
    l_int32          b2;         /* max b index in the vbox             */
};
typedef struct L_Box3d  L_BOX3D;

    /* Static median cut helper functions */
static PIXCMAP *pixcmapGenerateFromHisto(PIX *pixs, l_int32 depth,
                                         l_int32 *histo, l_int32 histosize,
                                         l_int32 sigbits);
static PIX *pixQuantizeWithColormap(PIX *pixs, l_int32 ditherflag,
                                    l_int32 outdepth,
                                    PIXCMAP *cmap, l_int32 *indexmap,
                                    l_int32 mapsize, l_int32 sigbits);
static void getColorIndexMedianCut(l_uint32 pixel, l_int32 rshift,
                                   l_uint32 mask, l_int32 sigbits,
                                   l_int32 *pindex);
static L_BOX3D *pixGetColorRegion(PIX *pixs, l_int32 sigbits,
                                  l_int32 subsample);
static l_int32 medianCutApply(l_int32 *histo, l_int32 sigbits,
                              L_BOX3D *vbox, L_BOX3D **pvbox1,
                              L_BOX3D **pvbox2);
static PIXCMAP *pixcmapGenerateFromMedianCuts(L_HEAP *lh, l_int32 *histo,
                                              l_int32 sigbits);
static l_int32 vboxGetAverageColor(L_BOX3D *vbox, l_int32 *histo,
                                   l_int32 sigbits, l_int32 index,
                                   l_int32  *prval, l_int32 *pgval,
                                   l_int32  *pbval);
static l_int32 vboxGetCount(L_BOX3D *vbox, l_int32 *histo, l_int32 sigbits);
static l_int32 vboxGetVolume(L_BOX3D *vbox);
static L_BOX3D *box3dCreate(l_int32 r1, l_int32 r2, l_int32 g1,
                            l_int32 g2, l_int32 b1, l_int32 b2);
static L_BOX3D *box3dCopy(L_BOX3D *vbox);


    /* 5 significant bits for each component is generally satisfactory */
static const l_int32  DEFAULT_SIG_BITS = 5;
static const l_int32  MAX_ITERS_ALLOWED = 5000;  /* prevents infinite looping */

    /* Specify fraction of vboxes made that are sorted on population alone.
     * The remaining vboxes are sorted on (population * vbox-volume).  */
static const l_float32  FRACT_BY_POPULATION = 0.85;

    /* To get the max value of 'dif' in the dithering color transfer,
     * divide DIF_CAP by 8. */
static const l_int32  DIF_CAP = 100;


#ifndef   NO_CONSOLE_IO
#define   DEBUG_MC_COLORS       0
#define   DEBUG_SPLIT_AXES      0
#endif   /* ~NO_CONSOLE_IO */


/*------------------------------------------------------------------------*
 *                                 High level                             *
 *------------------------------------------------------------------------*/
/*!
 *  pixMedianCutQuant()
 *
 *      Input:  pixs  (32 bpp; rgb color)
 *              ditherflag (1 for dither; 0 for no dither)
 *      Return: pixd (8 bit with colormap), or null on error
 *
 *  Notes:
 *      (1) Simple interface.  See pixMedianCutQuantGeneral() for
 *          use of defaulted parameters.
 */ 
PIX *
pixMedianCutQuant(PIX     *pixs,
                  l_int32  ditherflag)
{
    return pixMedianCutQuantGeneral(pixs, ditherflag,
                                    0, 256, DEFAULT_SIG_BITS, 1, 1);
}


/*!
 *  pixMedianCutQuantGeneral()
 *
 *      Input:  pixs  (32 bpp; rgb color)
 *              ditherflag (1 for dither; 0 for no dither)
 *              outdepth (output depth; valid: 0, 1, 2, 4, 8)
 *              maxcolors (between 2 and 256)
 *              sigbits (valid: 5 or 6; use 0 for default)
 *              maxsub (max subsampling, integer; use 0 for default;
 *                      1 for no subsampling)
 *              checkbw (1 to check if color content is very small,
 *                       0 to assume there is sufficient color)
 *      Return: pixd (8 bit with colormap), or null on error
 *
 *  Notes:
 *      (1) @maxcolors must be in the range [2 ... 256].
 *      (2) Use @outdepth = 0 to have the output depth computed as the
 *          minimum required to hold the actual colors found, given
 *          the @maxcolors constraint.
 *      (3) Use @outdepth = 1, 2, 4 or 8 to specify the output depth.
 *          In that case, @maxcolors must not exceed 2^(outdepth).
 *      (4) If there are fewer quantized colors in the image than @maxcolors,
 *          the colormap is simply generated from those colors.
 *      (5) @maxsub is the maximum allowed subsampling to be used in the
 *          computation of the color histogram and region of occupied
 *          color space.  The subsampling is chosen internally for
 *          efficiency, based on the image size, but this parameter
 *          limits it.  Use @maxsub = 0 for the internal default, which is the
 *          maximum allowed subsampling.  Use @maxsub = 1 to prevent
 *          subsampling.  In general use @maxsub >= 1 to specify the
 *          maximum subsampling to be allowed, where the actual subsampling
 *          will be the minimum of this value and the internally
 *          determined default value.
 *      (6) If the image appears gray because either most of the pixels
 *          are gray or most of the pixels are essentially black or white,
 *          the image is trivially quantized with a grayscale colormap.  The
 *          reason is that median cut divides the color space into rectangular
 *          regions, and it does a very poor job if all the pixels are
 *          near the diagonal of the color space cube.
 */ 
PIX *
pixMedianCutQuantGeneral(PIX     *pixs,
                         l_int32  ditherflag,
                         l_int32  outdepth,
                         l_int32  maxcolors,
                         l_int32  sigbits,
                         l_int32  maxsub,
                         l_int32  checkbw)
{
l_int32    i, subsample, histosize, smalln, ncolors, niters, popcolors;
l_int32    w, h, minside, factor, index, rval, gval, bval;
l_int32   *histo;
l_float32  pixfract, colorfract;
L_BOX3D   *vbox, *vbox1, *vbox2;
L_HEAP    *lh, *lhs;
PIX       *pixd;
PIXCMAP   *cmap;

    PROCNAME("pixMedianCutQuantGeneral");

    if (!pixs || pixGetDepth(pixs) != 32)
        return (PIX *)ERROR_PTR("pixs undefined or not 32 bpp", procName, NULL);
    if (maxcolors < 2 || maxcolors > 256)
        return (PIX *)ERROR_PTR("maxcolors not in [2...256]", procName, NULL);
    if (outdepth != 0 && outdepth != 1 && outdepth != 2 && outdepth != 4 &&
        outdepth != 8)
        return (PIX *)ERROR_PTR("outdepth not in {0,1,2,4,8}", procName, NULL);
    if (outdepth > 0 && (maxcolors > (1 << outdepth)))
        return (PIX *)ERROR_PTR("maxcolors > 2^(outdepth)", procName, NULL);
    if (sigbits == 0)
        sigbits = DEFAULT_SIG_BITS;
    else if (sigbits < 5 || sigbits > 6)
        return (PIX *)ERROR_PTR("sigbits not 5 or 6", procName, NULL);
    if (maxsub <= 0)
        maxsub = 10;  /* default will prevail for 10^7 pixels or less */

        /* Determine if the image has sufficient color content.
         * If pixfract << 1, most pixels are close to black or white.
         * If colorfract << 1, the pixels that are not near
         *   black or white have very little color.
         * If with little color, quantize with a grayscale colormap. */
    pixGetDimensions(pixs, &w, &h, NULL);
    if (checkbw) {
        minside = L_MIN(w, h);
        factor = L_MAX(1, minside / 400);
        pixColorFraction(pixs, 20, 244, 20, factor, &pixfract, &colorfract);
        if (pixfract * colorfract < 0.00025) {
            L_INFO_FLOAT2("\n  Pixel fraction neither white nor black = %6.3f"
                          "\n  Color fraction of those pixels = %6.3f"
                          "\n  Quantizing in gray",
                          procName, pixfract, colorfract);
            return pixConvertTo8(pixs, 1);
        }
    }

        /* Compute the color space histogram.  Default sampling
         * is about 10^5 pixels.  */
    if (maxsub == 1)
        subsample = 1;
    else {
        subsample = (l_int32)(sqrt((l_float64)(w * h) / 100000.));
        subsample = L_MAX(1, L_MIN(maxsub, subsample));
    }
    histo = pixMedianCutHisto(pixs, sigbits, subsample);
    histosize = 1 << (3 * sigbits);

        /* See if the number of quantized colors is less than maxcolors */
    ncolors = 0;
    smalln = TRUE;
    for (i = 0; i < histosize; i++) {
        if (histo[i])
            ncolors++;
        if (ncolors > maxcolors) {
            smalln = FALSE;
            break;
        }
    }
    if (smalln) {  /* finish up now */
        if (outdepth == 0) {
            if (ncolors <= 2)
                outdepth = 1;
            else if (ncolors <= 4)
                outdepth = 2;
            else if (ncolors <= 16)
                outdepth = 4;
            else 
                outdepth = 8;
        }
        cmap = pixcmapGenerateFromHisto(pixs, outdepth,
                                        histo, histosize, sigbits);
        pixd = pixQuantizeWithColormap(pixs, ditherflag, outdepth, cmap,
                                       histo, histosize, sigbits);
        FREE(histo);
        return pixd;
    }

        /* Initial vbox: minimum region in colorspace occupied by pixels */
    if (ditherflag || subsample > 1)  /* use full color space */
        vbox = box3dCreate(0, (1 << sigbits) - 1,
                           0, (1 << sigbits) - 1,
                           0, (1 << sigbits) - 1);
    else
        vbox = pixGetColorRegion(pixs, sigbits, subsample);
    vbox->npix = vboxGetCount(vbox, histo, sigbits);
    vbox->vol = vboxGetVolume(vbox);

        /* For a fraction 'popcolors' of the desired 'maxcolors',
         * generate median cuts based on population, putting
         * everything on a priority queue sorted by population. */
    lh = lheapCreate(0, L_SORT_DECREASING);
    lheapAdd(lh, vbox);
    ncolors = 1;
    niters = 0;
    popcolors = (l_int32)(FRACT_BY_POPULATION * maxcolors);
    while (1) {
        vbox = (L_BOX3D *)lheapRemove(lh);
        if (vboxGetCount(vbox, histo, sigbits) == 0)  { /* just put it back */
            lheapAdd(lh, vbox);
            continue;
        }
        medianCutApply(histo, sigbits, vbox, &vbox1, &vbox2);
        if (!vbox1) {
            L_WARNING("vbox1 not defined; shouldn't happen!", procName);
            break;
        }
        if (vbox1->vol > 1)
            vbox1->sortparam = vbox1->npix;
        FREE(vbox);
        lheapAdd(lh, vbox1);
        if (vbox2) {  /* vbox2 can be NULL */
            if (vbox2->vol > 1)
                vbox2->sortparam = vbox2->npix;
            lheapAdd(lh, vbox2);
            ncolors++;
        }
        if (ncolors >= popcolors)
            break;
        if (niters++ > MAX_ITERS_ALLOWED) {
            L_WARNING("infinite loop; perhaps too few pixels!", procName);
            break;
        }
    }

        /* Re-sort by the product of pixel occupancy times the size
         * in color space. */
    lhs = lheapCreate(0, L_SORT_DECREASING);
    while ((vbox = (L_BOX3D *)lheapRemove(lh))) {
        vbox->sortparam = vbox->npix * vbox->vol;
        lheapAdd(lhs, vbox);
    }
    lheapDestroy(&lh, TRUE);

        /* For the remaining (maxcolors - popcolors), generate the
         * median cuts using the (npix * vol) sorting. */
    while (1) {
        vbox = (L_BOX3D *)lheapRemove(lhs);
        if (vboxGetCount(vbox, histo, sigbits) == 0)  { /* just put it back */
            lheapAdd(lhs, vbox);
            continue;
        }
        medianCutApply(histo, sigbits, vbox, &vbox1, &vbox2);
        if (!vbox1) {
            L_WARNING("vbox1 not defined; shouldn't happen!", procName);
            break;
        }
        if (vbox1->vol > 1)
            vbox1->sortparam = vbox1->npix * vbox1->vol;
        FREE(vbox);
        lheapAdd(lhs, vbox1);
        if (vbox2) {  /* vbox2 can be NULL */
            if (vbox2->vol > 1)
                vbox2->sortparam = vbox2->npix * vbox2->vol;
            lheapAdd(lhs, vbox2);
            ncolors++;
        }
        if (ncolors >= maxcolors)
            break;
        if (niters++ > MAX_ITERS_ALLOWED) {
            L_WARNING("infinite loop; perhaps too few pixels!", procName);
            break;
        }
    }

        /* Re-sort by pixel occupancy.  This is not necessary,
         * but it makes a more useful listing.  */
    lh = lheapCreate(0, L_SORT_DECREASING);
    while ((vbox = (L_BOX3D *)lheapRemove(lhs))) {
        vbox->sortparam = vbox->npix;
/*        vbox->sortparam = vbox->npix * vbox->vol; */
        lheapAdd(lh, vbox);
    }
    lheapDestroy(&lhs, TRUE);

        /* Generate colormap from median cuts and quantize pixd */
    cmap = pixcmapGenerateFromMedianCuts(lh, histo, sigbits);
    if (outdepth == 0) {
        ncolors = pixcmapGetCount(cmap);
        if (ncolors <= 2)
            outdepth = 1;
        else if (ncolors <= 4)
            outdepth = 2;
        else if (ncolors <= 16)
            outdepth = 4;
        else 
            outdepth = 8;
    }
    pixd = pixQuantizeWithColormap(pixs, ditherflag, outdepth, cmap,
                                   histo, histosize, sigbits);

        /* Force darkest color to black if each component <= 4 */
    pixcmapGetRankIntensity(cmap, 0.0, &index);
    pixcmapGetColor(cmap, index, &rval, &gval, &bval);
    if (rval < 5 && gval < 5 && bval < 5)
        pixcmapResetColor(cmap, index, 0, 0, 0);

        /* Force lightest color to white if each component >= 252 */
    pixcmapGetRankIntensity(cmap, 1.0, &index);
    pixcmapGetColor(cmap, index, &rval, &gval, &bval);
    if (rval > 251 && gval > 251 && bval > 251)
        pixcmapResetColor(cmap, index, 255, 255, 255);

    lheapDestroy(&lh, TRUE);
    FREE(histo);
    return pixd;
}


/*!
 *  pixMedianCutQuantMixed()
 *
 *      Input:  pixs  (32 bpp; rgb color)
 *              ncolor (maximum number of colors assigned to pixels with
 *                      significant color)
 *              ngray (number of gray colors to be used; must be >= 2)
 *              darkthresh (threshold near black; if the lightest component
 *                          is below this, the pixel is not considered to
 *                          be gray or color; uses 0 for default)
 *              lightthresh (threshold near white; if the darkest component
 *                           is above this, the pixel is not considered to
 *                           be gray or color; use 0 for default)
 *              diffthresh (thresh for the max difference between component
 *                          values; for differences below this, the pixel
 *                          is considered to be gray; use 0 for default)
 *      Return: pixd (8 bpp cmapped), or null on error
 *
 *  Notes:
 *      (1) ncolor + ngray must not exceed 255.
 *      (2) The method makes use of pixMedianCutQuantGeneral() with
 *          minimal addition.
 *          (a) Preprocess the image, setting all pixels with little color
 *              to black, and populating an auxiliary 8 bpp image with the
 *              expected colormap values corresponding to the set of
 *              quantized gray values.
 *          (b) Color quantize the altered input image to n + 1 colors.
 *          (c) Augment the colormap with the gray indices, and
 *              substitute the gray quantized values from the auxiliary
 *              image for those in the color quantized output that had
 *              been quantized as black.
 *      (3) Median cut color quantization is relatively poor for grayscale
 *          images with many colors, when compared to octcube quantization.
 *          Thus, for images with both gray and color, it is important
 *          to quantize the gray pixels by another method.  Here, we
 *          are conservative in detecting color, preferring to use
 *          a few extra bits to encode colorful pixels that push them
 *          to gray.  This is particularly reasonable with this function,
 *          because it handles the gray and color pixels separately,
 *          using median cut color quantization for the color pixels
 *          and equal-bin grayscale quantization for the non-color pixels.
 */
PIX *
pixMedianCutQuantMixed(PIX     *pixs,
                       l_int32  ncolor,
                       l_int32  ngray,
                       l_int32  darkthresh,
                       l_int32  lightthresh,
                       l_int32  diffthresh)
{
l_int32    i, j, w, h, wplc, wplg, wpld, nc, unused, iscolor, factor, minside;
l_int32    rval, gval, bval, minval, maxval, val, grayval;
l_float32  pixfract, colorfract;
l_int32   *lut;
l_uint32  *datac, *datag, *datad, *linec, *lineg, *lined;
PIX       *pixc, *pixg, *pixd;
PIXCMAP   *cmap;

    PROCNAME("pixMedianCutQuantMixed");

    if (!pixs || pixGetDepth(pixs) != 32)
        return (PIX *)ERROR_PTR("pixs undefined or not 32 bpp", procName, NULL);
    if (ngray < 2)
        return (PIX *)ERROR_PTR("ngray < 2", procName, NULL);
    if (ncolor + ngray > 255)
        return (PIX *)ERROR_PTR("ncolor + ngray > 255", procName, NULL);
    if (darkthresh <= 0) darkthresh = 20;
    if (lightthresh <= 0) lightthresh = 244;
    if (diffthresh <= 0) diffthresh = 20;

        /* First check if this should be quantized in gray.
         * Use a more sensitive parameter for detecting color than with
         * pixMedianCutQuantGeneral(), because this function can handle
         * gray pixels well.  */
    pixGetDimensions(pixs, &w, &h, NULL);
    minside = L_MIN(w, h);
    factor = L_MAX(1, minside / 400);
    pixColorFraction(pixs, darkthresh, lightthresh, diffthresh, factor,
                     &pixfract, &colorfract);
    if (pixfract * colorfract < 0.0001) {
        L_INFO_FLOAT2("\n  Pixel fraction neither white nor black = %6.3f"
                      "\n  Color fraction of those pixels = %6.3f"
                      "\n  Quantizing in gray",
                      procName, pixfract, colorfract);
        pixg = pixConvertTo8(pixs, 0);
        pixd = pixThresholdOn8bpp(pixg, ngray, 1);
        pixDestroy(&pixg);
        return pixd;
    }

        /* OK, there is color in the image.
         * Preprocess to handle the gray pixels.  Set the color pixels in pixc
         * to black, and store their (eventual) colormap indices in pixg.*/
    pixc = pixCopy(NULL, pixs);
    pixg = pixCreate(w, h, 8);  /* color pixels will remain 0 here */
    datac = pixGetData(pixc);
    datag = pixGetData(pixg);
    wplc = pixGetWpl(pixc);
    wplg = pixGetWpl(pixg);
    lut = (l_int32 *)CALLOC(256, sizeof(l_int32));
    for (i = 0; i < 256; i++)
        lut[i] = ncolor + 1 + (i * (ngray - 1) + 128) / 255;
    for (i = 0; i < h; i++) {
        linec = datac + i * wplc;
        lineg = datag + i * wplg;
        for (j = 0; j < w; j++) {
            iscolor = FALSE;
            extractRGBValues(linec[j], &rval, &gval, &bval);
            minval = L_MIN(rval, gval);
            minval = L_MIN(minval, bval);
            maxval = L_MAX(rval, gval);
            maxval = L_MAX(maxval, bval);
            if (maxval >= darkthresh &&
                minval <= lightthresh &&
                maxval - minval >= diffthresh) {
                iscolor = TRUE;
            }
            if (!iscolor) {
                linec[j] = 0x0;  /* set to black */
                grayval = (maxval + minval) / 2;
                SET_DATA_BYTE(lineg, j, lut[grayval]);
            }
        }
    }

        /* Median cut on color pixels plus black */
    pixd = pixMedianCutQuantGeneral(pixc, FALSE, 8, ncolor + 1,
                                    DEFAULT_SIG_BITS, 1, 0);

        /* Augment the colormap with gray values.  The new cmap
         * indices should agree with the values previously stored in pixg. */
    cmap = pixGetColormap(pixd);
    nc = pixcmapGetCount(cmap);
    unused = ncolor  + 1 - nc;
    if (unused < 0)
        L_ERROR_INT("Too many colors: extra = %d", procName, -unused);
    if (unused > 0) {  /* fill in with black; these won't be used */
        L_INFO_INT("%d unused colors", procName, unused);
        for (i = 0; i < unused; i++)
            pixcmapAddColor(cmap, 0, 0, 0);
    }
    for (i = 0; i < ngray; i++) {
        grayval = (255 * i) / (ngray - 1);
        pixcmapAddColor(cmap, grayval, grayval, grayval);
    }

        /* Substitute cmap indices for the gray pixels into pixd */
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);
    for (i = 0; i < h; i++) {
        lined = datad + i * wpld;
        lineg = datag + i * wplg;
        for (j = 0; j < w; j++) {
            val = GET_DATA_BYTE(lineg, j);  /* if 0, it's a color pixel */
            if (val)
                SET_DATA_BYTE(lined, j, val);
        }
    }

    pixDestroy(&pixc);
    pixDestroy(&pixg);
    FREE(lut);
    return pixd;
}


/*!
 *  pixFewColorsMedianCutQuantMixed()
 *
 *      Input:  pixs (32 bpp rgb)
 *              ncolor (number of colors to be assigned to pixels with
 *                       significant color)
 *              ngray (number of gray colors to be used; must be >= 2)
 *              maxncolors (maximum number of colors to be returned
 *                         from pixColorsForQuantization(); use 0 for default)
 *              darkthresh (threshold near black; if the lightest component
 *                          is below this, the pixel is not considered to
 *                          be gray or color; use 0 for default)
 *              lightthresh (threshold near white; if the darkest component
 *                           is above this, the pixel is not considered to
 *                           be gray or color; use 0 for default)
 *              diffthresh (thresh for the max difference between component
 *                          values; for differences below this, the pixel
 *                          is considered to be gray; use 0 for default)
 *                          considered gray; use 0 for default)
 *      Return: pixd (8 bpp, median cut quantized for pixels that are
 *                    not gray; gray pixels are quantized separately
 *                    over the full gray range); null if too many colors
 *                    or on error
 *
 *  Notes:
 *      (1) This is the "few colors" version of pixMedianCutQuantMixed().
 *          It fails (returns NULL) if it finds more than maxncolors, but
 *          otherwise it gives the same result.
 *      (2) Recommended input parameters are:
 *              @maxncolors:  20
 *              @darkthresh:  20
 *              @lightthresh: 244
 *              @diffthresh:  15  (any higher can miss colors differing
 *                                 slightly from gray)
 *      (3) Both ncolor and ngray should be at least equal to maxncolors.
 *          If they're not, they are automatically increased, and a
 *          warning is given.
 *      (4) If very little color content is found, the input is
 *          converted to gray and quantized in equal intervals.
 *      (5) This can be useful for quantizing orthographically generated
 *          images such as color maps, where there may be more than 256 colors
 *          because of aliasing or jpeg artifacts on text or lines, but
 *          there are a relatively small number of solid colors.
 *      (6) Example of usage:
 *             // Try to quantize, using default values for mixed med cut
 *             Pix *pixq = pixFewColorsMedianCutQuantMixed(pixs, 100, 20,
 *                             0, 0, 0, 0);
 *             if (!pixq)  // too many colors; don't quantize
 *                 pixq = pixClone(pixs);
 */
PIX *
pixFewColorsMedianCutQuantMixed(PIX       *pixs,
                                l_int32    ncolor,
                                l_int32    ngray,
                                l_int32    maxncolors,
                                l_int32    darkthresh,
                                l_int32    lightthresh,
                                l_int32    diffthresh)
{
l_int32  ncolors, iscolor;
PIX     *pixg, *pixd;

    PROCNAME("pixFewColorsMedianCutQuantMixed");

    if (!pixs || pixGetDepth(pixs) != 32)
        return (PIX *)ERROR_PTR("pixs undefined or not 32 bpp", procName, NULL);
    if (maxncolors <= 0) maxncolors = 20;
    if (darkthresh <= 0) darkthresh = 20;
    if (lightthresh <= 0) lightthresh = 244;
    if (diffthresh <= 0) diffthresh = 15;
    if (ncolor < maxncolors) {
        L_WARNING_INT("ncolor too small; setting to %d", procName, maxncolors);
        ncolor = maxncolors;
    }
    if (ngray < maxncolors) {
        L_WARNING_INT("ngray too small; setting to %d", procName, maxncolors);
        ngray = maxncolors;
    }

        /* Estimate the color content and the number of colors required */
    pixColorsForQuantization(pixs, 15, &ncolors, &iscolor, 0);

        /* Note that maxncolors applies to all colors required to quantize,
         * both gray and colorful */
    if (ncolors > maxncolors)
        return (PIX *)ERROR_PTR("too many colors", procName, NULL);

        /* If no color, return quantized gray pix */
    if (!iscolor) {
        pixg = pixConvertTo8(pixs, 0);
        pixd = pixThresholdOn8bpp(pixg, ngray, 1);
        pixDestroy(&pixg);
        return pixd;
    }

        /* Use the mixed gray/color quantizer */
    return pixMedianCutQuantMixed(pixs, ncolor, ngray, darkthresh,
                                  lightthresh, diffthresh);
}



/*------------------------------------------------------------------------*
 *                        Median cut indexed histogram                    *
 *------------------------------------------------------------------------*/
/*!
 *  pixMedianCutHisto()
 *
 *      Input:  pixs  (32 bpp; rgb color)
 *              sigbits (valid: 5 or 6)
 *              subsample (integer > 0)
 *      Return: histo (1-d array, giving the number of pixels in
 *                     each quantized region of color space), or null on error
 *
 *  Notes:
 *      (1) Array is indexed by (3 * sigbits) bits.  The array size
 *          is 2^(3 * sigbits).
 *      (2) Indexing into the array from rgb uses red sigbits as
 *          most significant and blue as least.
 */
l_int32 *
pixMedianCutHisto(PIX     *pixs,
                  l_int32  sigbits,
                  l_int32  subsample)
{
l_int32    i, j, w, h, wpl, rshift, index, histosize;
l_int32   *histo;
l_uint32   mask, pixel;
l_uint32  *data, *line;

    PROCNAME("pixMedianCutHisto");

    if (!pixs)
        return (l_int32 *)ERROR_PTR("pixs not defined", procName, NULL);
    if (pixGetDepth(pixs) != 32)
        return (l_int32 *)ERROR_PTR("pixs not 32 bpp", procName, NULL);
    if (sigbits < 5 || sigbits > 6)
        return (l_int32 *)ERROR_PTR("sigbits not 5 or 6", procName, NULL);
    if (subsample <= 0)
        return (l_int32 *)ERROR_PTR("subsample not > 0", procName, NULL);

    histosize = 1 << (3 * sigbits);
    if ((histo = (l_int32 *)CALLOC(histosize, sizeof(l_int32))) == NULL)
        return (l_int32 *)ERROR_PTR("histo not made", procName, NULL);

    rshift = 8 - sigbits;
    mask = 0xff >> rshift;
    pixGetDimensions(pixs, &w, &h, NULL);
    data = pixGetData(pixs);
    wpl = pixGetWpl(pixs);
    for (i = 0; i < h; i += subsample) {
        line = data + i * wpl; 
        for (j = 0; j < w; j += subsample) {
            pixel = line[j];
            getColorIndexMedianCut(pixel, rshift, mask, sigbits, &index);
            histo[index]++;
        }
    }

    return histo;
}


/*------------------------------------------------------------------------*
 *                               Static helpers                           *
 *------------------------------------------------------------------------*/
/*!
 *  pixcmapGenerateFromHisto()
 *
 *      Input:  pixs  (32 bpp; rgb color)
 *              depth (of colormap)
 *              histo
 *              histosize
 *              sigbits
 *      Return: colormap, or null on error
 *
 *  Notes:
 *      (1) This is used when the number of colors in the histo
 *          is not greater than maxcolors.
 *      (2) As a side-effect, the histo becomes an inverse colormap,
 *          labeling the cmap indices for each existing color.
 */
static PIXCMAP *
pixcmapGenerateFromHisto(PIX      *pixs,
                         l_int32   depth,
                         l_int32  *histo,
                         l_int32   histosize,
                         l_int32   sigbits)
{
l_int32   i, index, shift, rval, gval, bval;
l_uint32  mask;
PIXCMAP  *cmap;

    PROCNAME("pixcmapGenerateFromHisto");

    if (!pixs)
        return (PIXCMAP *)ERROR_PTR("pixs not defined", procName, NULL);
    if (pixGetDepth(pixs) != 32)
        return (PIXCMAP *)ERROR_PTR("pixs not 32 bpp", procName, NULL);
    if (!histo)
        return (PIXCMAP *)ERROR_PTR("histo not defined", procName, NULL);

        /* Capture the rgb values of each occupied cube in the histo,
         * and re-label the histo value with the colormap index. */
    cmap = pixcmapCreate(depth);
    shift = 8 - sigbits;
    mask = 0xff >> shift;
    for (i = 0, index = 0; i < histosize; i++) {
        if (histo[i]) {
            rval = (i >> (2 * sigbits)) << shift;
            gval = ((i >> sigbits) & mask) << shift;
            bval = (i & mask) << shift;
            pixcmapAddColor(cmap, rval, gval, bval);
            histo[i] = index++;
        }
    }

    return cmap;
}


/*!
 *  pixQuantizeWithColormap()
 *
 *      Input:  pixs  (32 bpp; rgb color)
 *              ditherflag (1 for dither; 0 for no dither)
 *              outdepth
 *              cmap
 *              indexmap
 *              histosize
 *              sigbits
 *      Return: pixd (quantized to colormap), or null on error
 *
 *  Notes:
 *      (1) The indexmap is a LUT that takes the rgb indices of the
 *          pixel and returns the index into the colormap.
 *      (2) If ditherflag is 1, @outdepth is ignored and the output
 *          depth is set to 8.
 */
static PIX *
pixQuantizeWithColormap(PIX      *pixs,
                        l_int32   ditherflag,
                        l_int32   outdepth,
                        PIXCMAP  *cmap,
                        l_int32  *indexmap,
                        l_int32   mapsize,
                        l_int32   sigbits)
{
l_uint8   *bufu8r, *bufu8g, *bufu8b;
l_int32    i, j, w, h, wpls, wpld, rshift, index, cmapindex;
l_int32    rval, gval, bval, rc, gc, bc;
l_int32    dif, val1, val2, val3;
l_int32   *buf1r, *buf1g, *buf1b, *buf2r, *buf2g, *buf2b;
l_uint32  *datas, *datad, *lines, *lined;
l_uint32   mask, pixel;
PIX       *pixd;

    PROCNAME("pixQuantizeWithColormap");

    if (!pixs || pixGetDepth(pixs) != 32)
        return (PIX *)ERROR_PTR("pixs not 32 bpp", procName, NULL);
    if (!cmap)
        return (PIX *)ERROR_PTR("cmap not defined", procName, NULL);
    if (!indexmap)
        return (PIX *)ERROR_PTR("indexmap not defined", procName, NULL);
    if (ditherflag)
        outdepth = 8;

    pixGetDimensions(pixs, &w, &h, NULL);
    pixd = pixCreate(w, h, outdepth);
    pixSetColormap(pixd, cmap);
    pixCopyResolution(pixd, pixs);
    pixCopyInputFormat(pixd, pixs);
    datas = pixGetData(pixs);
    datad = pixGetData(pixd);
    wpls = pixGetWpl(pixs);
    wpld = pixGetWpl(pixd);

    rshift = 8 - sigbits;
    mask = 0xff >> rshift;
    if (ditherflag == 0) {
        for (i = 0; i < h; i++) {
            lines = datas + i * wpls;
            lined = datad + i * wpld;
            if (outdepth == 1) {
                for (j = 0; j < w; j++) {
                    pixel = lines[j];
                    getColorIndexMedianCut(pixel, rshift, mask,
                                           sigbits, &index);
                    if (indexmap[index])
                        SET_DATA_BIT(lined, j);
                }
            }
            else if (outdepth == 2) {
                for (j = 0; j < w; j++) {
                    pixel = lines[j];
                    getColorIndexMedianCut(pixel, rshift, mask,
                                           sigbits, &index);
                    SET_DATA_DIBIT(lined, j, indexmap[index]);
                }
            }
            else if (outdepth == 4) {
                for (j = 0; j < w; j++) {
                    pixel = lines[j];
                    getColorIndexMedianCut(pixel, rshift, mask,
                                           sigbits, &index);
                    SET_DATA_QBIT(lined, j, indexmap[index]);
                }
            }
            else {  /* outdepth == 8 */
                for (j = 0; j < w; j++) {
                    pixel = lines[j];
                    getColorIndexMedianCut(pixel, rshift, mask,
                                           sigbits, &index);
                    SET_DATA_BYTE(lined, j, indexmap[index]);
                }
            }
        }
    }
    else {  /* ditherflag == 1 */
        bufu8r = (l_uint8 *)CALLOC(w, sizeof(l_uint8));
        bufu8g = (l_uint8 *)CALLOC(w, sizeof(l_uint8));
        bufu8b = (l_uint8 *)CALLOC(w, sizeof(l_uint8));
        buf1r = (l_int32 *)CALLOC(w, sizeof(l_int32));
        buf1g = (l_int32 *)CALLOC(w, sizeof(l_int32));
        buf1b = (l_int32 *)CALLOC(w, sizeof(l_int32));
        buf2r = (l_int32 *)CALLOC(w, sizeof(l_int32));
        buf2g = (l_int32 *)CALLOC(w, sizeof(l_int32));
        buf2b = (l_int32 *)CALLOC(w, sizeof(l_int32));
        if (!bufu8r || !bufu8g || !bufu8b)
            return (PIX *)ERROR_PTR("uint8 line buf not made", procName, NULL);
        if (!buf1r || !buf1g || !buf1b || !buf2r || !buf2g || !buf2b)
            return (PIX *)ERROR_PTR("mono line buf not made", procName, NULL);

            /* Start by priming buf2; line 1 is above line 2 */
        pixGetRGBLine(pixs, 0, bufu8r, bufu8g, bufu8b);
        for (j = 0; j < w; j++) {
            buf2r[j] = 64 * bufu8r[j];
            buf2g[j] = 64 * bufu8g[j];
            buf2b[j] = 64 * bufu8b[j];
        }

        for (i = 0; i < h - 1; i++) {
                /* Swap data 2 --> 1, and read in new line 2 */
            memcpy(buf1r, buf2r, 4 * w);
            memcpy(buf1g, buf2g, 4 * w);
            memcpy(buf1b, buf2b, 4 * w);
            pixGetRGBLine(pixs, i + 1, bufu8r, bufu8g, bufu8b);
            for (j = 0; j < w; j++) {
                buf2r[j] = 64 * bufu8r[j];
                buf2g[j] = 64 * bufu8g[j];
                buf2b[j] = 64 * bufu8b[j];
            }

                /* Dither */
            lined = datad + i * wpld;
            for (j = 0; j < w - 1; j++) {
                rval = buf1r[j] / 64;
                gval = buf1g[j] / 64;
                bval = buf1b[j] / 64;
                index = ((rval >> rshift) << (2 * sigbits)) +
                        ((gval >> rshift) << sigbits) + (bval >> rshift);
                cmapindex = indexmap[index];
                SET_DATA_BYTE(lined, j, cmapindex);
                pixcmapGetColor(cmap, cmapindex, &rc, &gc, &bc);

                dif = buf1r[j] / 8 - 8 * rc;
                if (dif > DIF_CAP) dif = DIF_CAP;
                if (dif < -DIF_CAP) dif = -DIF_CAP;
                if (dif != 0) {
                    val1 = buf1r[j + 1] + 3 * dif;
                    val2 = buf2r[j] + 3 * dif;
                    val3 = buf2r[j + 1] + 2 * dif;
                    if (dif > 0) {
                        buf1r[j + 1] = L_MIN(16383, val1);
                        buf2r[j] = L_MIN(16383, val2);
                        buf2r[j + 1] = L_MIN(16383, val3);
                    }
                    else if (dif < 0) {
                        buf1r[j + 1] = L_MAX(0, val1);
                        buf2r[j] = L_MAX(0, val2);
                        buf2r[j + 1] = L_MAX(0, val3);
                    }
                }

                dif = buf1g[j] / 8 - 8 * gc;
                if (dif > DIF_CAP) dif = DIF_CAP;
                if (dif < -DIF_CAP) dif = -DIF_CAP;
                if (dif != 0) {
                    val1 = buf1g[j + 1] + 3 * dif;
                    val2 = buf2g[j] + 3 * dif;
                    val3 = buf2g[j + 1] + 2 * dif;
                    if (dif > 0) {
                        buf1g[j + 1] = L_MIN(16383, val1);
                        buf2g[j] = L_MIN(16383, val2);
                        buf2g[j + 1] = L_MIN(16383, val3);
                    }
                    else if (dif < 0) {
                        buf1g[j + 1] = L_MAX(0, val1);
                        buf2g[j] = L_MAX(0, val2);
                        buf2g[j + 1] = L_MAX(0, val3);
                    }
                }

                dif = buf1b[j] / 8 - 8 * bc;
                if (dif > DIF_CAP) dif = DIF_CAP;
                if (dif < -DIF_CAP) dif = -DIF_CAP;
                if (dif != 0) {
                    val1 = buf1b[j + 1] + 3 * dif;
                    val2 = buf2b[j] + 3 * dif;
                    val3 = buf2b[j + 1] + 2 * dif;
                    if (dif > 0) {
                        buf1b[j + 1] = L_MIN(16383, val1);
                        buf2b[j] = L_MIN(16383, val2);
                        buf2b[j + 1] = L_MIN(16383, val3);
                    }
                    else if (dif < 0) {
                        buf1b[j + 1] = L_MAX(0, val1);
                        buf2b[j] = L_MAX(0, val2);
                        buf2b[j + 1] = L_MAX(0, val3);
                    }
                }
            }

                /* Get last pixel in row; no downward propagation */
            rval = buf1r[w - 1] / 64;
            gval = buf1g[w - 1] / 64;
            bval = buf1b[w - 1] / 64;
            index = ((rval >> rshift) << (2 * sigbits)) +
                    ((gval >> rshift) << sigbits) + (bval >> rshift);
            SET_DATA_BYTE(lined, w - 1, indexmap[index]);
        }

            /* Get last row of pixels; no leftward propagation */
        lined = datad + (h - 1) * wpld;
        for (j = 0; j < w; j++) {
            rval = buf2r[j] / 64;
            gval = buf2g[j] / 64;
            bval = buf2b[j] / 64;
            index = ((rval >> rshift) << (2 * sigbits)) +
                    ((gval >> rshift) << sigbits) + (bval >> rshift);
            SET_DATA_BYTE(lined, j, indexmap[index]);
        }

        FREE(bufu8r);
        FREE(bufu8g);
        FREE(bufu8b);
        FREE(buf1r);
        FREE(buf1g);
        FREE(buf1b);
        FREE(buf2r);
        FREE(buf2g);
        FREE(buf2b);
    }

    return pixd;
}


/*!
 *  getColorIndexMedianCut()
 *
 *      Input:  pixel (32 bit rgb)
 *              rshift (of component: 8 - sigbits)
 *              mask (over sigbits)
 *              sigbits
 *              &index (<return> rgb index value)
 *      Return: void
 *
 *  Notes:
 *      (1) This is used on each pixel in the source image.  No checking
 *          is done on input values.
 */
static void
getColorIndexMedianCut(l_uint32  pixel,
                       l_int32   rshift,
                       l_uint32  mask,
                       l_int32   sigbits,
                       l_int32  *pindex)
{
l_int32  rval, gval, bval;

    rval = pixel >> (24 + rshift);
    gval = (pixel >> (16 + rshift)) & mask;
    bval = (pixel >> (8 + rshift)) & mask;
    *pindex = (rval << (2 * sigbits)) + (gval << sigbits) + bval;
    return;
}


/*!
 *  pixGetColorRegion()
 *
 *      Input:  pixs  (32 bpp; rgb color)
 *              sigbits (valid: 5, 6)
 *              subsample (integer > 0)
 *      Return: vbox (minimum 3D box in color space enclosing all pixels),
 *              or null on error
 *
 *  Notes:
 *      (1) Computes the minimum 3D box in color space enclosing all
 *          pixels in the image.
 */
static L_BOX3D *
pixGetColorRegion(PIX     *pixs,
                  l_int32  sigbits,
                  l_int32  subsample)
{
l_int32    rmin, rmax, gmin, gmax, bmin, bmax, rval, gval, bval;
l_int32    w, h, wpl, i, j, rshift;
l_uint32   mask, pixel;
l_uint32  *data, *line;

    PROCNAME("pixGetColorRegion");

    if (!pixs)
        return (L_BOX3D *)ERROR_PTR("pixs not defined", procName, NULL);

    rmin = gmin = bmin = 1000000;
    rmax = gmax = bmax = 0;
    rshift = 8 - sigbits;
    mask = 0xff >> rshift;
    pixGetDimensions(pixs, &w, &h, NULL);
    data = pixGetData(pixs);
    wpl = pixGetWpl(pixs);
    for (i = 0; i < h; i += subsample) {
        line = data + i * wpl;
        for (j = 0; j < w; j += subsample) {
            pixel = line[j];
            rval = pixel >> (24 + rshift);
            gval = (pixel >> (16 + rshift)) & mask;
            bval = (pixel >> (8 + rshift)) & mask;
            if (rval < rmin)
                rmin = rval;
            else if (rval > rmax)
                rmax = rval;
            if (gval < gmin)
                gmin = gval;
            else if (gval > gmax)
                gmax = gval;
            if (bval < bmin)
                bmin = bval;
            else if (bval > bmax)
                bmax = bval;
        }
    }

    return box3dCreate(rmin, rmax, gmin, gmax, bmin, bmax);
}


/*!
 *  medianCutApply()
 *
 *      Input:  histo  (array; in rgb colorspace)
 *              sigbits
 *              vbox (input 3D box)
 *              &vbox1, vbox2 (<return> vbox split in two parts)
 *      Return: 0 if OK, 1 on error
 */
static l_int32
medianCutApply(l_int32   *histo,
               l_int32    sigbits,
               L_BOX3D   *vbox,
               L_BOX3D  **pvbox1,
               L_BOX3D  **pvbox2)
{
l_int32   i, j, k, sum, rw, gw, bw, maxw, index;
l_int32   total, left, right;
l_int32   partialsum[128];
L_BOX3D  *vbox1, *vbox2;

    PROCNAME("medianCutApply");

    if (!histo) 
        return ERROR_INT("histo not defined", procName, 1);
    if (!vbox) 
        return ERROR_INT("vbox not defined", procName, 1);
    if (!pvbox1 || !pvbox2) 
        return ERROR_INT("&vbox1 and &vbox2 not both defined", procName, 1);

    *pvbox1 = *pvbox2 = NULL;
    if (vboxGetCount(vbox, histo, sigbits) == 0)
        return ERROR_INT("no pixels in vbox", procName, 1);

        /* If the vbox occupies just one element in color space, it can't
         * be split.  Leave the 'sortparam' field at 0, so that it goes to
         * the tail of the priority queue and stays there, thereby avoiding
         * an infinite loop (take off, put back on the head) if it
         * happens to be the most populous box! */
    rw = vbox->r2 - vbox->r1 + 1;
    gw = vbox->g2 - vbox->g1 + 1;
    bw = vbox->b2 - vbox->b1 + 1;
    if (rw == 1 && gw == 1 && bw == 1) {
        *pvbox1 = box3dCopy(vbox);
        return 0;
    }

        /* Select the longest axis for splitting */
    maxw = L_MAX(rw, gw);
    maxw = L_MAX(maxw, bw);
#if  DEBUG_SPLIT_AXES
    if (rw == maxw)
        fprintf(stderr, "red split\n");
    else if (gw == maxw)
        fprintf(stderr, "green split\n");
    else
        fprintf(stderr, "blue split\n");
#endif  /* DEBUG_SPLIT_AXES */

        /* Find the partial sum arrays along the selected axis. */
    total = 0;
    if (maxw == rw) {
        for (i = vbox->r1; i <= vbox->r2; i++) {
            sum = 0;
            for (j = vbox->g1; j <= vbox->g2; j++) {
                for (k = vbox->b1; k <= vbox->b2; k++) {
                    index = (i << (2 * sigbits)) + (j << sigbits) + k;
                    sum += histo[index];
                }
            }
            total += sum;
            partialsum[i] = total;
        }
    }
    else if (maxw == gw) {
        for (i = vbox->g1; i <= vbox->g2; i++) {
            sum = 0;
            for (j = vbox->r1; j <= vbox->r2; j++) {
                for (k = vbox->b1; k <= vbox->b2; k++) {
                    index = (i << sigbits) + (j << (2 * sigbits)) + k;
                    sum += histo[index];
                }
            }
            total += sum;
            partialsum[i] = total;
        }
    }
    else {  /* maxw == bw */
        for (i = vbox->b1; i <= vbox->b2; i++) {
            sum = 0;
            for (j = vbox->r1; j <= vbox->r2; j++) {
                for (k = vbox->g1; k <= vbox->g2; k++) {
                    index = i + (j << (2 * sigbits)) + (k << sigbits);
                    sum += histo[index];
                }
            }
            total += sum;
            partialsum[i] = total;
        }
    }

        /* Determine the cut planes, making sure that two vboxes
         * are always produced.  Generate the two vboxes and compute
         * the sum in each of them.  Choose the cut plane within
         * the greater of the (left, right) sides of the bin in which
         * the median pixel resides.  Here's the surprise: go halfway
         * into that side.  By doing that, you technically move away
         * from "median cut," but in the process a significant number
         * of low-count vboxes are produced, allowing much better
         * reproduction of low-count spot colors. */
    if (maxw == rw) {
        for (i = vbox->r1; i <= vbox->r2; i++) {
            if (partialsum[i] > total / 2) {
                vbox1 = box3dCopy(vbox);
                vbox2 = box3dCopy(vbox);
                left = i - vbox->r1;
                right = vbox->r2 - i;
                if (left <= right)
                    vbox1->r2 = L_MIN(vbox->r2 - 1, i + right / 2);
                else  /* left > right */
                    vbox1->r2 = L_MAX(vbox->r1, i - 1 - left / 2);
                vbox2->r1 = vbox1->r2 + 1;
                break;
            }
        }
    }
    else if (maxw == gw) {
        for (i = vbox->g1; i <= vbox->g2; i++) {
            if (partialsum[i] > total / 2) {
                vbox1 = box3dCopy(vbox);
                vbox2 = box3dCopy(vbox);
                left = i - vbox->g1;
                right = vbox->g2 - i;
                if (left <= right)
                    vbox1->g2 = L_MIN(vbox->g2 - 1, i + right / 2);
                else  /* left > right */
                    vbox1->g2 = L_MAX(vbox->g1, i - 1 - left / 2);
                vbox2->g1 = vbox1->g2 + 1;
                break;
            }
        }
    }
    else {  /* maxw == bw */
        for (i = vbox->b1; i <= vbox->b2; i++) {
            if (partialsum[i] > total / 2) {
                vbox1 = box3dCopy(vbox);
                vbox2 = box3dCopy(vbox);
                left = i - vbox->b1;
                right = vbox->b2 - i;
                if (left <= right)
                    vbox1->b2 = L_MIN(vbox->b2 - 1, i + right / 2);
                else  /* left > right */
                    vbox1->b2 = L_MAX(vbox->b1, i - 1 - left / 2);
                vbox2->b1 = vbox1->b2 + 1;
                break;
            }
        }
    }
    vbox1->npix = vboxGetCount(vbox1, histo, sigbits);
    vbox2->npix = vboxGetCount(vbox2, histo, sigbits);
    vbox1->vol = vboxGetVolume(vbox1);
    vbox2->vol = vboxGetVolume(vbox2);
    *pvbox1 = vbox1;
    *pvbox2 = vbox2;

    return 0;
}


/*!
 *  pixcmapGenerateFromMedianCuts()
 *
 *      Input:  lh (priority queue of pointers to vboxes)
 *              histo
 *              sigbits (valid: 5 or 6)
 *      Return: cmap, or null on error
 *
 *  Notes:
 *      (1) Each vbox in the heap represents a color in the colormap.
 *      (2) As a side-effect, the histo becomes an inverse colormap,
 *          where the part of the array correpsonding to each vbox
 *          is labeled with the cmap index for that vbox.  Then
 *          for each rgb pixel, the colormap index is found directly
 *          by mapping the rgb value to the histo array index.
 */
static PIXCMAP *
pixcmapGenerateFromMedianCuts(L_HEAP   *lh,
                              l_int32  *histo,
                              l_int32   sigbits)
{
l_int32   index, rval, gval, bval;
L_BOX3D  *vbox;
PIXCMAP  *cmap;

    PROCNAME("pixcmapGenerateFromMedianCuts");

    if (!lh) 
        return (PIXCMAP *)ERROR_PTR("lh not defined", procName, NULL);
    if (!histo) 
        return (PIXCMAP *)ERROR_PTR("histo not defined", procName, NULL);

    rval = gval = bval = 0;  /* make compiler happy */
    cmap = pixcmapCreate(8);
    index = 0;
    while (lheapGetCount(lh) > 0) {
        vbox = (L_BOX3D *)lheapRemove(lh);
        vboxGetAverageColor(vbox, histo, sigbits, index, &rval, &gval, &bval);
        pixcmapAddColor(cmap, rval, gval, bval);
        FREE(vbox);
        index++;
    }

    return cmap;
}


/*!
 *  vboxGetAverageColor()
 *
 *      Input:  vbox (3d region of color space for one quantized color)
 *              histo
 *              sigbits (valid: 5 or 6)
 *              index (if >= 0, assign to all colors in histo in this vbox)
 *              &rval, &gval, &bval (<returned> average color)
 *      Return: cmap, or null on error
 *
 *  Notes:
 *      (1) The vbox represents one color in the colormap.
 *      (2) If index >= 0, as a side-effect, all array elements in
 *          the histo corresponding to the vbox are labeled with this
 *          cmap index for that vbox.  Otherwise, the histo array
 *          is not changed.
 */
static l_int32
vboxGetAverageColor(L_BOX3D  *vbox,
                    l_int32  *histo,
                    l_int32   sigbits,
                    l_int32   index,
                    l_int32  *prval,
                    l_int32  *pgval,
                    l_int32  *pbval)
{
l_int32  i, j, k, ntot, mult, histoindex, rsum, gsum, bsum;

    PROCNAME("vboxGetAverageColor");

    if (!vbox) 
        return ERROR_INT("vbox not defined", procName, 1);
    if (!histo) 
        return ERROR_INT("histo not defined", procName, 1);
    if (!prval || !pgval || !pbval)
        return ERROR_INT("&p*val not all defined", procName, 1);

    *prval = *pgval = *pbval = 0;
    ntot = 0;
    mult = 1 << (8 - sigbits);
    rsum = gsum = bsum = 0;
    for (i = vbox->r1; i <= vbox->r2; i++) {
        for (j = vbox->g1; j <= vbox->g2; j++) {
            for (k = vbox->b1; k <= vbox->b2; k++) {
                 histoindex = (i << (2 * sigbits)) + (j << sigbits) + k;
                 ntot += histo[histoindex];
                 rsum += (l_int32)(histo[histoindex] * (i + 0.5) * mult);
                 gsum += (l_int32)(histo[histoindex] * (j + 0.5) * mult);
                 bsum += (l_int32)(histo[histoindex] * (k + 0.5) * mult);
                 if (index >= 0)
                     histo[histoindex] = index;
            }
        }
    }

    if (ntot == 0) {
        *prval = mult * (vbox->r1 + vbox->r2 + 1) / 2;
        *pgval = mult * (vbox->g1 + vbox->g2 + 1) / 2;
        *pbval = mult * (vbox->b1 + vbox->b2 + 1) / 2;
    }
    else {
        *prval = rsum / ntot;
        *pgval = gsum / ntot;
        *pbval = bsum / ntot;
    }

#if  DEBUG_MC_COLORS
    fprintf(stderr, "ntot[%d] = %d: [%d, %d, %d], (%d, %d, %d)\n",
            index, ntot, vbox->r2 - vbox->r1 + 1,
            vbox->g2 - vbox->g1 + 1, vbox->b2 - vbox->b1 + 1,
            *prval, *pgval, *pbval);
#endif  /* DEBUG_MC_COLORS */

    return 0;
}


/*!
 *  vboxGetCount()
 *
 *      Input:  vbox (3d region of color space for one quantized color)
 *              histo
 *              sigbits (valid: 5 or 6)
 *      Return: number of image pixels in this region, or 0 on error
 */
static l_int32
vboxGetCount(L_BOX3D  *vbox,
             l_int32  *histo,
             l_int32   sigbits)
{
l_int32  i, j, k, npix, index;

    PROCNAME("vboxGetCount");

    if (!vbox) 
        return ERROR_INT("vbox not defined", procName, 0);
    if (!histo) 
        return ERROR_INT("histo not defined", procName, 0);

    npix = 0;
    for (i = vbox->r1; i <= vbox->r2; i++) {
        for (j = vbox->g1; j <= vbox->g2; j++) {
            for (k = vbox->b1; k <= vbox->b2; k++) {
                 index = (i << (2 * sigbits)) + (j << sigbits) + k;
                 npix += histo[index];
            }
        }
    }

    return npix;
}


/*!
 *  vboxGetVolume()
 *
 *      Input:  vbox (3d region of color space for one quantized color)
 *      Return: quantized volume of vbox, or 0 on error
 */
static l_int32
vboxGetVolume(L_BOX3D  *vbox)
{
    PROCNAME("vboxGetVolume");

    if (!vbox) 
        return ERROR_INT("vbox not defined", procName, 0);

    return ((vbox->r2 - vbox->r1 + 1) * (vbox->g2 - vbox->g1 + 1) *
            (vbox->b2 - vbox->b1 + 1));
}


static L_BOX3D *
box3dCreate(l_int32  r1,
            l_int32  r2,
            l_int32  g1,
            l_int32  g2,
            l_int32  b1,
            l_int32  b2)
{
L_BOX3D  *vbox;

    vbox = (L_BOX3D *)CALLOC(1, sizeof(L_BOX3D));
    vbox->r1 = r1;
    vbox->r2 = r2;
    vbox->g1 = g1;
    vbox->g2 = g2;
    vbox->b1 = b1;
    vbox->b2 = b2;
    return vbox;
}


/*
 *  Note: don't copy the sortparam.
 */
static L_BOX3D *
box3dCopy(L_BOX3D  *vbox)
{
L_BOX3D  *vboxc;

    PROCNAME("box3dCopy");

    if (!vbox) 
        return (L_BOX3D *)ERROR_PTR("vbox not defined", procName, NULL);

    vboxc = box3dCreate(vbox->r1, vbox->r2, vbox->g1, vbox->g2,
                        vbox->b1, vbox->b2);
    vboxc->npix = vbox->npix;
    vboxc->vol = vbox->vol;
    return vboxc;
}