colorquant2.c

Go to the documentation of this file.

/*====================================================================*
 -  Copyright (C) 2001 Leptonica.  All rights reserved.
 -
 -  Redistribution and use in source and binary forms, with or without
 -  modification, are permitted provided that the following conditions
 -  are met:
 -  1. Redistributions of source code must retain the above copyright
 -     notice, this list of conditions and the following disclaimer.
 -  2. Redistributions in binary form must reproduce the above
 -     copyright notice, this list of conditions and the following
 -     disclaimer in the documentation and/or other materials
 -     provided with the distribution.
 -
 -  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 -  ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 -  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 -  A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL ANY
 -  CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 -  EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 -  PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
 -  PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
 -  OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
 -  NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 -  SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *====================================================================*/

#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                             *
 *------------------------------------------------------------------------*/
PIX *
pixMedianCutQuant(PIX     *pixs,
                  l_int32  ditherflag)
{
    return pixMedianCutQuantGeneral(pixs, ditherflag,
                                    0, 256, DEFAULT_SIG_BITS, 1, 1);
}


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("\n  Pixel fraction neither white nor black = %6.3f"
                   "\n  Color fraction of those pixels = %6.3f"
                   "\n  Quantizing in gray\n",
                   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);
        LEPT_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!\n", procName);
            break;
        }
        if (vbox1->vol > 1)
            vbox1->sortparam = vbox1->npix;
        LEPT_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!\n", 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!\n", procName);
            break;
        }
        if (vbox1->vol > 1)
            vbox1->sortparam = vbox1->npix * vbox1->vol;
        LEPT_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!\n", 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);
    LEPT_FREE(histo);
    return pixd;
}


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("\n  Pixel fraction neither white nor black = %6.3f"
                      "\n  Color fraction of those pixels = %6.3f"
                      "\n  Quantizing in gray\n",
                      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 *)LEPT_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("Too many colors: extra = %d\n", procName, -unused);
    if (unused > 0) {  /* fill in with black; these won't be used */
        L_INFO("%d unused colors\n", 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);
    LEPT_FREE(lut);
    return pixd;
}


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("ncolor too small; setting to %d\n", procName, maxncolors);
        ncolor = maxncolors;
    }
    if (ngray < maxncolors) {
        L_WARNING("ngray too small; setting to %d\n", 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                    *
 *------------------------------------------------------------------------*/
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 *)LEPT_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                           *
 *------------------------------------------------------------------------*/
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;
}


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, success;
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 */
        success = TRUE;
        bufu8r = bufu8g = bufu8b = NULL;
        buf1r = buf1g = buf1b = buf2r = buf2g = buf2b = NULL;
        bufu8r = (l_uint8 *)LEPT_CALLOC(w, sizeof(l_uint8));
        bufu8g = (l_uint8 *)LEPT_CALLOC(w, sizeof(l_uint8));
        bufu8b = (l_uint8 *)LEPT_CALLOC(w, sizeof(l_uint8));
        buf1r = (l_int32 *)LEPT_CALLOC(w, sizeof(l_int32));
        buf1g = (l_int32 *)LEPT_CALLOC(w, sizeof(l_int32));
        buf1b = (l_int32 *)LEPT_CALLOC(w, sizeof(l_int32));
        buf2r = (l_int32 *)LEPT_CALLOC(w, sizeof(l_int32));
        buf2g = (l_int32 *)LEPT_CALLOC(w, sizeof(l_int32));
        buf2b = (l_int32 *)LEPT_CALLOC(w, sizeof(l_int32));
        if (!bufu8r || !bufu8g || !bufu8b || !buf1r || !buf1g ||
            !buf1b || !buf2r || !buf2g || !buf2b) {
            L_ERROR("buffer not made\n", procName);
            success = FALSE;
            goto buffer_cleanup;
        }

            /* 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 {
                        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 {
                        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 {
                        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]);
        }

buffer_cleanup:
        LEPT_FREE(bufu8r);
        LEPT_FREE(bufu8g);
        LEPT_FREE(bufu8b);
        LEPT_FREE(buf1r);
        LEPT_FREE(buf1g);
        LEPT_FREE(buf1b);
        LEPT_FREE(buf2r);
        LEPT_FREE(buf2g);
        LEPT_FREE(buf2b);
        if (!success) pixDestroy(&pixd);
    }

    return pixd;
}


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


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


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 (pvbox1) *pvbox1 = NULL;
    if (pvbox2) *pvbox2 = NULL;
    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);

    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. */
    vbox1 = vbox2 = NULL;
    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;
            }
        }
    }
    *pvbox1 = vbox1;
    *pvbox2 = vbox2;
    if (!vbox1)
        return ERROR_INT("vbox1 not made; shouldn't happen", procName, 1);
    if (!vbox2)
        return ERROR_INT("vbox2 not made; shouldn't happen", procName, 1);
    vbox1->npix = vboxGetCount(vbox1, histo, sigbits);
    vbox2->npix = vboxGetCount(vbox2, histo, sigbits);
    vbox1->vol = vboxGetVolume(vbox1);
    vbox2->vol = vboxGetVolume(vbox2);

    return 0;
}


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);
        LEPT_FREE(vbox);
        index++;
    }

    return cmap;
}


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


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


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


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