colorcontent.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 "allheaders.h"

/* ----------------------------------------------------------------------- *
 *      Builds an image of the color content, on a per-pixel basis,        *
 *      as a measure of the amount of divergence of each color             *
 *      component (R,G,B) from gray.                                       *
 * ----------------------------------------------------------------------- */
l_ok
pixColorContent(PIX     *pixs,
                l_int32  rwhite,
                l_int32  gwhite,
                l_int32  bwhite,
                l_int32  mingray,
                PIX    **ppixr,
                PIX    **ppixg,
                PIX    **ppixb)
{
l_int32    w, h, d, i, j, wplc, wplr, wplg, wplb;
l_int32    rval, gval, bval, rgdiff, rbdiff, gbdiff, maxval, colorval;
l_int32   *rtab, *gtab, *btab;
l_uint32   pixel;
l_uint32  *datac, *datar, *datag, *datab, *linec, *liner, *lineg, *lineb;
NUMA      *nar, *nag, *nab;
PIX       *pixc;   /* rgb */
PIX       *pixr, *pixg, *pixb;   /* 8 bpp grayscale */
PIXCMAP   *cmap;

    PROCNAME("pixColorContent");

    if (!ppixr && !ppixg && !ppixb)
        return ERROR_INT("no return val requested", procName, 1);
    if (ppixr) *ppixr = NULL;
    if (ppixg) *ppixg = NULL;
    if (ppixb) *ppixb = NULL;
    if (!pixs)
        return ERROR_INT("pixs not defined", procName, 1);
    if (mingray < 0) mingray = 0;
    pixGetDimensions(pixs, &w, &h, &d);
    if (mingray > 255)
        return ERROR_INT("mingray > 255", procName, 1);
    if (rwhite < 0 || gwhite < 0 || bwhite < 0)
        return ERROR_INT("some white vals are negative", procName, 1);
    if ((rwhite || gwhite || bwhite) && (rwhite * gwhite * bwhite == 0))
        return ERROR_INT("white vals not all zero or all nonzero", procName, 1);

    cmap = pixGetColormap(pixs);
    if (!cmap && d != 32)
        return ERROR_INT("pixs neither cmapped nor 32 bpp", procName, 1);
    if (cmap)
        pixc = pixRemoveColormap(pixs, REMOVE_CMAP_TO_FULL_COLOR);
    else
        pixc = pixClone(pixs);

    pixr = pixg = pixb = NULL;
    pixGetDimensions(pixc, &w, &h, NULL);
    if (ppixr) {
        pixr = pixCreate(w, h, 8);
        datar = pixGetData(pixr);
        wplr = pixGetWpl(pixr);
        *ppixr = pixr;
    }
    if (ppixg) {
        pixg = pixCreate(w, h, 8);
        datag = pixGetData(pixg);
        wplg = pixGetWpl(pixg);
        *ppixg = pixg;
    }
    if (ppixb) {
        pixb = pixCreate(w, h, 8);
        datab = pixGetData(pixb);
        wplb = pixGetWpl(pixb);
        *ppixb = pixb;
    }

    datac = pixGetData(pixc);
    wplc = pixGetWpl(pixc);
    if (rwhite) {  /* all white pt vals are nonzero */
        nar = numaGammaTRC(1.0, 0, rwhite);
        rtab = numaGetIArray(nar);
        nag = numaGammaTRC(1.0, 0, gwhite);
        gtab = numaGetIArray(nag);
        nab = numaGammaTRC(1.0, 0, bwhite);
        btab = numaGetIArray(nab);
    }
    for (i = 0; i < h; i++) {
        linec = datac + i * wplc;
        if (pixr)
            liner = datar + i * wplr;
        if (pixg)
            lineg = datag + i * wplg;
        if (pixb)
            lineb = datab + i * wplb;
        for (j = 0; j < w; j++) {
            pixel = linec[j];
            extractRGBValues(pixel, &rval, &gval, &bval);
            if (rwhite) {  /* color correct for white point */
                rval = rtab[rval];
                gval = gtab[gval];
                bval = btab[bval];
            }
            if (mingray > 0) {  /* dark pixels have no color value */
                maxval = L_MAX(rval, gval);
                maxval = L_MAX(maxval, bval);
                if (maxval < mingray)
                    continue;  /* colorval = 0 for each component */
            }
            rgdiff = L_ABS(rval - gval);
            rbdiff = L_ABS(rval - bval);
            gbdiff = L_ABS(gval - bval);
            if (pixr) {
                colorval = (rgdiff + rbdiff) / 2;
                SET_DATA_BYTE(liner, j, colorval);
            }
            if (pixg) {
                colorval = (rgdiff + gbdiff) / 2;
                SET_DATA_BYTE(lineg, j, colorval);
            }
            if (pixb) {
                colorval = (rbdiff + gbdiff) / 2;
                SET_DATA_BYTE(lineb, j, colorval);
            }
        }
    }

    if (rwhite) {
        numaDestroy(&nar);
        numaDestroy(&nag);
        numaDestroy(&nab);
        LEPT_FREE(rtab);
        LEPT_FREE(gtab);
        LEPT_FREE(btab);
    }
    pixDestroy(&pixc);
    return 0;
}


/* ----------------------------------------------------------------------- *
 *      Finds the 'amount' of color in an image, on a per-pixel basis,     *
 *      as a measure of the difference of the pixel color from gray.       *
 * ----------------------------------------------------------------------- */
PIX *
pixColorMagnitude(PIX     *pixs,
                  l_int32  rwhite,
                  l_int32  gwhite,
                  l_int32  bwhite,
                  l_int32  type)
{
l_int32    w, h, d, i, j, wplc, wpld;
l_int32    rval, gval, bval, rdist, gdist, bdist, colorval;
l_int32    rgdist, rbdist, gbdist, mindist, maxdist, minval, maxval;
l_int32   *rtab, *gtab, *btab;
l_uint32   pixel;
l_uint32  *datac, *datad, *linec, *lined;
NUMA      *nar, *nag, *nab;
PIX       *pixc, *pixd;
PIXCMAP   *cmap;

    PROCNAME("pixColorMagnitude");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    pixGetDimensions(pixs, &w, &h, &d);
    if (type != L_MAX_DIFF_FROM_AVERAGE_2 && type != L_MAX_MIN_DIFF_FROM_2 &&
        type != L_MAX_DIFF)
        return (PIX *)ERROR_PTR("invalid type", procName, NULL);
    if (rwhite < 0 || gwhite < 0 || bwhite < 0)
        return (PIX *)ERROR_PTR("some white vals are negative", procName, NULL);
    if ((rwhite || gwhite || bwhite) && (rwhite * gwhite * bwhite == 0))
        return (PIX *)ERROR_PTR("white vals not all zero or all nonzero",
                                procName, NULL);

    cmap = pixGetColormap(pixs);
    if (!cmap && d != 32)
        return (PIX *)ERROR_PTR("pixs not cmapped or 32 bpp", procName, NULL);
    if (cmap)
        pixc = pixRemoveColormap(pixs, REMOVE_CMAP_TO_FULL_COLOR);
    else
        pixc = pixClone(pixs);

    pixd = pixCreate(w, h, 8);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);
    datac = pixGetData(pixc);
    wplc = pixGetWpl(pixc);
    if (rwhite) {  /* all white pt vals are nonzero */
        nar = numaGammaTRC(1.0, 0, rwhite);
        rtab = numaGetIArray(nar);
        nag = numaGammaTRC(1.0, 0, gwhite);
        gtab = numaGetIArray(nag);
        nab = numaGammaTRC(1.0, 0, bwhite);
        btab = numaGetIArray(nab);
    }
    for (i = 0; i < h; i++) {
        linec = datac + i * wplc;
        lined = datad + i * wpld;
        for (j = 0; j < w; j++) {
            pixel = linec[j];
            extractRGBValues(pixel, &rval, &gval, &bval);
            if (rwhite) {  /* color correct for white point */
                rval = rtab[rval];
                gval = gtab[gval];
                bval = btab[bval];
            }
            if (type == L_MAX_DIFF_FROM_AVERAGE_2) {
                rdist = ((gval + bval ) / 2 - rval);
                rdist = L_ABS(rdist);
                gdist = ((rval + bval ) / 2 - gval);
                gdist = L_ABS(gdist);
                bdist = ((rval + gval ) / 2 - bval);
                bdist = L_ABS(bdist);
                colorval = L_MAX(rdist, gdist);
                colorval = L_MAX(colorval, bdist);
            } else if (type == L_MAX_MIN_DIFF_FROM_2) {  /* intermediate dist */
                rgdist = L_ABS(rval - gval);
                rbdist = L_ABS(rval - bval);
                gbdist = L_ABS(gval - bval);
                maxdist = L_MAX(rgdist, rbdist);
                if (gbdist >= maxdist) {
                    colorval = maxdist;
                } else {  /* gbdist is smallest or intermediate */
                    mindist = L_MIN(rgdist, rbdist);
                    colorval = L_MAX(mindist, gbdist);
                }
            } else {  /* type == L_MAX_DIFF */
                minval = L_MIN(rval, gval);
                minval = L_MIN(minval, bval);
                maxval = L_MAX(rval, gval);
                maxval = L_MAX(maxval, bval);
                colorval = maxval - minval;
            }
            SET_DATA_BYTE(lined, j, colorval);
        }
    }

    if (rwhite) {
        numaDestroy(&nar);
        numaDestroy(&nag);
        numaDestroy(&nab);
        LEPT_FREE(rtab);
        LEPT_FREE(gtab);
        LEPT_FREE(btab);
    }
    pixDestroy(&pixc);
    return pixd;
}


/* ----------------------------------------------------------------------- *
 *      Generates a mask over pixels that have sufficient color and        *
 *      are not too close to gray pixels.                                  *
 * ----------------------------------------------------------------------- */
PIX *
pixMaskOverColorPixels(PIX     *pixs,
                       l_int32  threshdiff,
                       l_int32  mindist)
{
l_int32    w, h, d, i, j, wpls, wpld, size;
l_int32    rval, gval, bval, minval, maxval;
l_uint32  *datas, *datad, *lines, *lined;
PIX       *pixc, *pixd;
PIXCMAP   *cmap;

    PROCNAME("pixMaskOverColorPixels");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    pixGetDimensions(pixs, &w, &h, &d);

    cmap = pixGetColormap(pixs);
    if (!cmap && d != 32)
        return (PIX *)ERROR_PTR("pixs not cmapped or 32 bpp", procName, NULL);
    if (cmap)
        pixc = pixRemoveColormap(pixs, REMOVE_CMAP_TO_FULL_COLOR);
    else
        pixc = pixClone(pixs);

    pixd = pixCreate(w, h, 1);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);
    datas = pixGetData(pixc);
    wpls = pixGetWpl(pixc);
    for (i = 0; i < h; i++) {
        lines = datas + i * wpls;
        lined = datad + i * wpld;
        for (j = 0; j < w; j++) {
            extractRGBValues(lines[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 - minval >= threshdiff)
                SET_DATA_BIT(lined, j);
        }
    }

    if (mindist > 1) {
        size = 2 * (mindist - 1) + 1;
        pixErodeBrick(pixd, pixd, size, size);
    }

    pixDestroy(&pixc);
    return pixd;
}


/* ----------------------------------------------------------------------- *
 *      Generates a mask over pixels that have RGB color components        *
 *      within the prescribed range (a cube in RGB color space)            *
 * ----------------------------------------------------------------------- */
PIX *
pixMaskOverColorRange(PIX     *pixs,
                      l_int32  rmin,
                      l_int32  rmax,
                      l_int32  gmin,
                      l_int32  gmax,
                      l_int32  bmin,
                      l_int32  bmax)
{
l_int32    w, h, d, i, j, wpls, wpld;
l_int32    rval, gval, bval;
l_uint32  *datas, *datad, *lines, *lined;
PIX       *pixc, *pixd;
PIXCMAP   *cmap;

    PROCNAME("pixMaskOverColorRange");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    pixGetDimensions(pixs, &w, &h, &d);

    cmap = pixGetColormap(pixs);
    if (!cmap && d != 32)
        return (PIX *)ERROR_PTR("pixs not cmapped or 32 bpp", procName, NULL);
    if (cmap)
        pixc = pixRemoveColormap(pixs, REMOVE_CMAP_TO_FULL_COLOR);
    else
        pixc = pixClone(pixs);

    pixd = pixCreate(w, h, 1);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);
    datas = pixGetData(pixc);
    wpls = pixGetWpl(pixc);
    for (i = 0; i < h; i++) {
        lines = datas + i * wpls;
        lined = datad + i * wpld;
        for (j = 0; j < w; j++) {
            extractRGBValues(lines[j], &rval, &gval, &bval);
            if (rval < rmin || rval > rmax) continue;
            if (gval < gmin || gval > gmax) continue;
            if (bval < bmin || bval > bmax) continue;
            SET_DATA_BIT(lined, j);
        }
    }

    pixDestroy(&pixc);
    return pixd;
}


/* ----------------------------------------------------------------------- *
 *   Finds the fraction of pixels with "color" that are not close to black *
 * ----------------------------------------------------------------------- */
l_ok
pixColorFraction(PIX        *pixs,
                 l_int32     darkthresh,
                 l_int32     lightthresh,
                 l_int32     diffthresh,
                 l_int32     factor,
                 l_float32  *ppixfract,
                 l_float32  *pcolorfract)
{
l_int32    i, j, w, h, wpl, rval, gval, bval, minval, maxval;
l_int32    total, npix, ncolor;
l_uint32   pixel;
l_uint32  *data, *line;

    PROCNAME("pixColorFraction");

    if (ppixfract) *ppixfract = 0.0;
    if (pcolorfract) *pcolorfract = 0.0;
    if (!ppixfract || !pcolorfract)
        return ERROR_INT("&pixfract and &colorfract not defined",
                         procName, 1);
    if (!pixs || pixGetDepth(pixs) != 32)
        return ERROR_INT("pixs not defined or not 32 bpp", procName, 1);

    pixGetDimensions(pixs, &w, &h, NULL);
    data = pixGetData(pixs);
    wpl = pixGetWpl(pixs);
    npix = ncolor = total = 0;
    for (i = 0; i < h; i += factor) {
        line = data + i * wpl;
        for (j = 0; j < w; j += factor) {
            total++;
            pixel = line[j];
            extractRGBValues(pixel, &rval, &gval, &bval);
            minval = L_MIN(rval, gval);
            minval = L_MIN(minval, bval);
            if (minval > lightthresh)  /* near white */
                continue;
            maxval = L_MAX(rval, gval);
            maxval = L_MAX(maxval, bval);
            if (maxval < darkthresh)  /* near black */
                continue;

            npix++;
            if (maxval - minval >= diffthresh)
                ncolor++;
        }
    }

    if (npix == 0) {
        L_WARNING("No pixels found for consideration\n", procName);
        return 0;
    }
    *ppixfract = (l_float32)npix / (l_float32)total;
    *pcolorfract = (l_float32)ncolor / (l_float32)npix;
    return 0;
}


/* ----------------------------------------------------------------------- *
 *     Determine if there are significant color regions in a page image    *
 * ----------------------------------------------------------------------- */
l_ok
pixFindColorRegions(PIX        *pixs,
                    PIX        *pixm,
                    l_int32     factor,
                    l_int32     lightthresh,
                    l_int32     darkthresh,
                    l_int32     mindiff,
                    l_int32     colordiff,
                    l_float32   edgefract,
                    l_float32  *pcolorfract,
                    PIX       **pcolormask1,
                    PIX       **pcolormask2,
                    PIXA       *pixadb)
{
l_int32    w, h, count, rval, gval, bval, aveval, proceed;
l_float32  ratio;
l_uint32  *carray;
BOXA      *boxa1, *boxa2;
PIX       *pix1, *pix2, *pix3, *pix4, *pix5, *pixm1, *pixm2, *pixm3;

    PROCNAME("pixFindColorRegions");

    if (pcolormask1) *pcolormask1 = NULL;
    if (pcolormask2) *pcolormask2 = NULL;
    if (!pcolorfract)
        return ERROR_INT("&colorfract not defined", procName, 1);
    *pcolorfract = 0.0;
    if (!pixs || pixGetDepth(pixs) != 32)
        return ERROR_INT("pixs not defined or not 32 bpp", procName, 1);
    if (factor < 1) factor = 1;
    if (lightthresh < 0) lightthresh = 210;  /* defaults */
    if (darkthresh < 0) darkthresh = 70;
    if (mindiff < 0) mindiff = 10;
    if (colordiff < 0) colordiff = 90;
    if (edgefract < 0.0 || edgefract > 1.0) edgefract = 0.05;

        /* Check if pixm covers most of the image.  If so, just return. */
    pixGetDimensions(pixs, &w, &h, NULL);
    if (pixm) {
        pixCountPixels(pixm, &count, NULL);
        ratio = (l_float32)count / ((l_float32)(w) * h);
        if (ratio > 0.7) {
            if (pixadb) L_INFO("pixm has big fg: %f5.2\n", procName, ratio);
            return 0;
        }
    }

        /* Get the light background color.  Use the average component value
         * and select the lightest of 10 buckets.  Require that it is
         * reddish and, using lightthresh, not too dark. */
    pixGetRankColorArray(pixs, 10, L_SELECT_AVERAGE, factor, &carray, 0, 0);
    if (!carray)
        return ERROR_INT("rank color array not made", procName, 1);
    extractRGBValues(carray[9], &rval, &gval, &bval);
    if (pixadb) L_INFO("lightest background color: (r,g,b) = (%d,%d,%d)\n",
                       procName, rval, gval, bval);
    proceed = TRUE;
    if ((rval < bval - 2) || (rval < gval - 2)) {
        if (pixadb) L_INFO("background not reddish\n", procName);
        proceed = FALSE;
    }
    aveval = (rval + gval + bval) / 3;
    if (aveval < lightthresh) {
        if (pixadb) L_INFO("background too dark\n", procName);
        proceed = FALSE;
    }
    if (pixadb) {
        pix1 = pixDisplayColorArray(carray, 10, 120, 3, 6);
        pixaAddPix(pixadb, pix1, L_INSERT);
    }
    LEPT_FREE(carray);
    if (proceed == FALSE) return 0;

        /* Make a mask pixm1 over the dark pixels in the image:
         * convert to gray using the average of the components;
         * threshold using darkthresh; do a small dilation;
         * combine with pixm. */
    pix1 = pixConvertRGBToGray(pixs, 0.33, 0.34, 0.33);
    if (pixadb) pixaAddPix(pixadb, pix1, L_COPY);
    pixm1 = pixThresholdToBinary(pix1, darkthresh);
    pixDilateBrick(pixm1, pixm1, 7, 7);
    if (pixadb) pixaAddPix(pixadb, pixm1, L_COPY);
    if (pixm) {
        pixOr(pixm1, pixm1, pixm);
        if (pixadb) pixaAddPix(pixadb, pixm1, L_COPY);
    }
    pixDestroy(&pix1);

        /* Make masks over pixels that are bluish, or greenish, or
           have a very large color saturation (max - min) value. */
    pixm2 = pixConvertRGBToBinaryArb(pixs, -1.0, 0.0, 1.0, mindiff,
                                     L_SELECT_IF_GTE);  /* b - r */
    if (pixadb) pixaAddPix(pixadb, pixm2, L_COPY);
    pix1 = pixConvertRGBToBinaryArb(pixs, -1.0, 1.0, 0.0, mindiff,
                                    L_SELECT_IF_GTE);  /* g - r */
    if (pixadb) pixaAddPix(pixadb, pix1, L_COPY);
    pixOr(pixm2, pixm2, pix1);
    pixDestroy(&pix1);
    pix1 = pixConvertRGBToGrayMinMax(pixs, L_CHOOSE_MAXDIFF);
    pix2 = pixThresholdToBinary(pix1, colordiff);
    pixInvert(pix2, pix2);
    if (pixadb) pixaAddPix(pixadb, pix2, L_COPY);
    pixOr(pixm2, pixm2, pix2);
    if (pixadb) pixaAddPix(pixadb, pixm2, L_COPY);
    pixDestroy(&pix1);
    pixDestroy(&pix2);

        /* Subtract the dark pixels represented by pixm1.
         * pixm2 now holds all the color pixels of interest  */
    pixSubtract(pixm2, pixm2, pixm1);
    pixDestroy(&pixm1);
    if (pixadb) pixaAddPix(pixadb, pixm2, L_COPY);

        /* But we're not quite finished.  Remove pixels from any component
         * that is touching the image border.  False color pixels can
         * sometimes be found there if the image is much darker near
         * the border, due to oxidation or reduced illumination.  Also
         * remove any pixels within the normalized fraction %distfract
         * of the image border. */
    pixm3 = pixRemoveBorderConnComps(pixm2, 8);
    pixDestroy(&pixm2);
    if (edgefract > 0.0) {
        pix2 = pixMakeFrameMask(w, h, edgefract, 1.0, edgefract, 1.0);
        pixAnd(pixm3, pixm3, pix2);
        pixDestroy(&pix2);
    }
    if (pixadb) pixaAddPix(pixadb, pixm3, L_COPY);

        /* Get the fraction of light color pixels */
    pixCountPixels(pixm3, &count, NULL);
    *pcolorfract = (l_float32)count / ((l_float32)(w) * h);
    if (pixadb) {
        if (count == 0)
            L_INFO("no light color pixels found\n", procName);
        else
            L_INFO("fraction of light color pixels = %5.3f\n", procName,
                   *pcolorfract);
    }

        /* Debug: extract the color pixels from pixs */
    if (pixadb && count > 0) {
            /* Use pixm3 to extract the color pixels */
        pix3 = pixCreateTemplate(pixs);
        pixSetAll(pix3);
        pixCombineMasked(pix3, pixs, pixm3);
        pixaAddPix(pixadb, pix3, L_INSERT);

            /* Use additional filtering to extract the color pixels */
        pix3 = pixCloseSafeBrick(NULL, pixm3, 15, 15);
        pixaAddPix(pixadb, pix3, L_INSERT);
        pix5 = pixCreateTemplate(pixs);
        pixSetAll(pix5);
        pixCombineMasked(pix5, pixs, pix3);
        pixaAddPix(pixadb, pix5, L_INSERT);

            /* Get the combined bounding boxes of the mask components
             * in pix3, and extract those pixels from pixs. */
        boxa1 = pixConnCompBB(pix3, 8);
        boxa2 = boxaCombineOverlaps(boxa1, NULL);
        pix4 = pixCreateTemplate(pix3);
        pixMaskBoxa(pix4, pix4, boxa2, L_SET_PIXELS);
        pixaAddPix(pixadb, pix4, L_INSERT);
        pix5 = pixCreateTemplate(pixs);
        pixSetAll(pix5);
        pixCombineMasked(pix5, pixs, pix4);
        pixaAddPix(pixadb, pix5, L_INSERT);
        boxaDestroy(&boxa1);
        boxaDestroy(&boxa2);
    }
    pixaAddPix(pixadb, pixs, L_COPY);

        /* Optional colormask returns */
    if (pcolormask2 && count > 0)
        *pcolormask2 = pixCloseSafeBrick(NULL, pixm3, 15, 15);
    if (pcolormask1 && count > 0)
        *pcolormask1 = pixm3;
    else
        pixDestroy(&pixm3);
    return 0;
}


/* ----------------------------------------------------------------------- *
 *      Finds the number of perceptually significant gray intensities      *
 *      in a grayscale image.                                              *
 * ----------------------------------------------------------------------- */
l_ok
pixNumSignificantGrayColors(PIX       *pixs,
                            l_int32    darkthresh,
                            l_int32    lightthresh,
                            l_float32  minfract,
                            l_int32    factor,
                            l_int32   *pncolors)
{
l_int32  i, w, h, count, mincount, ncolors;
NUMA    *na;

    PROCNAME("pixNumSignificantGrayColors");

    if (!pncolors)
        return ERROR_INT("&ncolors not defined", procName, 1);
    *pncolors = 0;
    if (!pixs || pixGetDepth(pixs) != 8)
        return ERROR_INT("pixs not defined or not 8 bpp", procName, 1);
    if (darkthresh < 0) darkthresh = 20;  /* defaults */
    if (lightthresh < 0) lightthresh = 236;
    if (minfract < 0.0) minfract = 0.0001;
    if (minfract > 1.0)
        return ERROR_INT("minfract > 1.0", procName, 1);
    if (minfract >= 0.001)
        L_WARNING("minfract too big; likely to underestimate ncolors\n",
                  procName);
    if (lightthresh > 255 || darkthresh >= lightthresh)
        return ERROR_INT("invalid thresholds", procName, 1);
    if (factor < 1) factor = 1;

    pixGetDimensions(pixs, &w, &h, NULL);
    mincount = (l_int32)(minfract * w * h * factor * factor);
    if ((na = pixGetGrayHistogram(pixs, factor)) == NULL)
        return ERROR_INT("na not made", procName, 1);
    ncolors = 2;  /* add in black and white */
    for (i = darkthresh; i <= lightthresh; i++) {
        numaGetIValue(na, i, &count);
        if (count >= mincount)
            ncolors++;
    }

    *pncolors = ncolors;
    numaDestroy(&na);
    return 0;
}


/* ----------------------------------------------------------------------- *
 *   Identifies images where color quantization will cause posterization   *
 *   due to the existence of many colors in low-gradient regions.          *
 * ----------------------------------------------------------------------- */
l_ok
pixColorsForQuantization(PIX      *pixs,
                         l_int32   thresh,
                         l_int32  *pncolors,
                         l_int32  *piscolor,
                         l_int32   debug)
{
l_int32    w, h, d, minside, factor;
l_float32  pixfract, colorfract;
PIX       *pixt, *pixsc, *pixg, *pixe, *pixb, *pixm;
PIXCMAP   *cmap;

    PROCNAME("pixColorsForQuantization");

    if (piscolor) *piscolor = 0;
    if (!pncolors)
        return ERROR_INT("&ncolors not defined", procName, 1);
    *pncolors = 0;
    if (!pixs)
        return ERROR_INT("pixs not defined", procName, 1);
    if ((cmap = pixGetColormap(pixs)) != NULL) {
        *pncolors = pixcmapGetCount(cmap);
        if (piscolor)
            pixcmapHasColor(cmap, piscolor);
        return 0;
    }

    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 8 && d != 32)
        return ERROR_INT("pixs not 8 or 32 bpp", procName, 1);
    if (thresh <= 0)
        thresh = 15;

        /* First test if 32 bpp has any significant color; if not,
         * convert it to gray.  Colors whose average values are within
         * 20 of black or 8 of white are ignored because they're not
         * very 'colorful'.  If less than 2.5/10000 of the pixels have
         * significant color, consider the image to be gray. */
    minside = L_MIN(w, h);
    if (d == 8) {
        pixt = pixClone(pixs);
    } else {  /* d == 32 */
        factor = L_MAX(1, minside / 400);
        pixColorFraction(pixs, 20, 248, 30, factor, &pixfract, &colorfract);
        if (pixfract * colorfract < 0.00025) {
            pixt = pixGetRGBComponent(pixs, COLOR_RED);
            d = 8;
        } else {  /* d == 32 */
            pixt = pixClone(pixs);
            if (piscolor)
                *piscolor = 1;
        }
    }

        /* If the smallest side is less than 1000, do not downscale.
         * If it is in [1000 ... 2000), downscale by 2x.  If it is >= 2000,
         * downscale by 4x.  Factors of 2 are chosen for speed.  The
         * actual resolution at which subsequent calculations take place
         * is not strongly dependent on downscaling.  */
    factor = L_MAX(1, minside / 500);
    if (factor == 1)
        pixsc = pixCopy(NULL, pixt);  /* to be sure pixs is unchanged */
    else if (factor == 2 || factor == 3)
        pixsc = pixScaleAreaMap2(pixt);
    else
        pixsc = pixScaleAreaMap(pixt, 0.25, 0.25);

        /* Basic edge mask generation procedure:
         *   ~ work on a grayscale image
         *   ~ get a 1 bpp edge mask by using an edge filter and
         *     thresholding to get fg pixels at the edges
         *   ~ for gray, dilate with a 3x3 brick Sel to get mask over
         *     all pixels within a distance of 1 pixel from the nearest
         *     edge pixel
         *   ~ for color, dilate with a 7x7 brick Sel to get mask over
         *     all pixels within a distance of 3 pixels from the nearest
         *     edge pixel  */
    if (d == 8)
        pixg = pixClone(pixsc);
    else  /* d == 32 */
        pixg = pixConvertRGBToLuminance(pixsc);
    pixe = pixSobelEdgeFilter(pixg, L_ALL_EDGES);
    pixb = pixThresholdToBinary(pixe, thresh);
    pixInvert(pixb, pixb);
    if (d == 8)
        pixm = pixMorphSequence(pixb, "d3.3", 0);
    else
        pixm = pixMorphSequence(pixb, "d7.7", 0);

        /* Mask the near-edge pixels to white, and count the colors.
         * If grayscale, don't count colors within 20 levels of
         * black or white, and only count colors with a fraction
         * of at least 1/10000 of the image pixels.
         * If color, count the number of level 4 octcubes that
         * contain at least 20 pixels.  These magic numbers are guesses
         * as to what might work, based on a small data set.  Results
         * should not be overly sensitive to their actual values. */
    if (d == 8) {
        pixSetMasked(pixg, pixm, 0xff);
        if (debug) pixWrite("junkpix8.png", pixg, IFF_PNG);
        pixNumSignificantGrayColors(pixg, 20, 236, 0.0001, 1, pncolors);
    } else {  /* d == 32 */
        pixSetMasked(pixsc, pixm, 0xffffffff);
        if (debug) pixWrite("junkpix32.png", pixsc, IFF_PNG);
        pixNumberOccupiedOctcubes(pixsc, 4, 20, -1, pncolors);
    }

    pixDestroy(&pixt);
    pixDestroy(&pixsc);
    pixDestroy(&pixg);
    pixDestroy(&pixe);
    pixDestroy(&pixb);
    pixDestroy(&pixm);
    return 0;
}


/* ----------------------------------------------------------------------- *
 *               Finds the number of unique colors in an image             *
 * ----------------------------------------------------------------------- */
l_ok
pixNumColors(PIX      *pixs,
             l_int32   factor,
             l_int32  *pncolors)
{
l_int32    w, h, d, i, j, wpl, hashsize, sum, count;
l_int32    rval, gval, bval, val;
l_int32   *inta;
l_uint32   pixel;
l_uint32  *data, *line;
PIXCMAP   *cmap;

    PROCNAME("pixNumColors");

    if (!pncolors)
        return ERROR_INT("&ncolors not defined", procName, 1);
    *pncolors = 0;
    if (!pixs)
        return ERROR_INT("pixs not defined", procName, 1);
    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 2 && d != 4 && d != 8 && d != 32)
        return ERROR_INT("d not in {2, 4, 8, 32}", procName, 1);
    if (factor < 1) factor = 1;

    data = pixGetData(pixs);
    wpl = pixGetWpl(pixs);
    sum = 0;
    if (d != 32) {  /* grayscale */
        if ((inta = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32))) == NULL)
            return ERROR_INT("calloc failure for inta", procName, 1);
        for (i = 0; i < h; i += factor) {
            line = data + i * wpl;
            for (j = 0; j < w; j += factor) {
                if (d == 8)
                    val = GET_DATA_BYTE(line, j);
                else if (d == 4)
                    val = GET_DATA_QBIT(line, j);
                else  /* d == 2 */
                    val = GET_DATA_DIBIT(line, j);
                inta[val] = 1;
            }
        }
        for (i = 0; i < 256; i++)
            if (inta[i]) sum++;
        *pncolors = sum;
        LEPT_FREE(inta);

        cmap = pixGetColormap(pixs);
        if (cmap && factor == 1) {
            count = pixcmapGetCount(cmap);
            if (sum != count)
                L_WARNING("colormap size %d differs from actual colors\n",
                          procName, count);
        }
        return 0;
    }

        /* 32 bpp rgb; quit if we get above 256 colors */
    hashsize = 5507;  /* big and prime; collisions are not likely */
    if ((inta = (l_int32 *)LEPT_CALLOC(hashsize, sizeof(l_int32))) == NULL)
        return ERROR_INT("calloc failure with hashsize", procName, 1);
    for (i = 0; i < h; i += factor) {
        line = data + i * wpl;
        for (j = 0; j < w; j += factor) {
            pixel = line[j];
            extractRGBValues(pixel, &rval, &gval, &bval);
            val = (137 * rval + 269 * gval + 353 * bval) % hashsize;
            if (inta[val] == 0) {
                inta[val] = 1;
                sum++;
                if (sum > 256) {
                    LEPT_FREE(inta);
                    return 0;
                }
            }
        }
    }

    *pncolors = sum;
    LEPT_FREE(inta);
    return 0;
}


/* ----------------------------------------------------------------------- *
 *       Find the most "populated" colors in the image (and quantize)      *
 * ----------------------------------------------------------------------- */
l_ok
pixGetMostPopulatedColors(PIX        *pixs,
                          l_int32     sigbits,
                          l_int32     factor,
                          l_int32     ncolors,
                          l_uint32  **parray,
                          PIXCMAP   **pcmap)
{
l_int32  n, i, rgbindex, rval, gval, bval;
NUMA    *nahisto, *naindex;

    PROCNAME("pixGetMostPopulatedColors");

    if (!parray && !pcmap)
        return ERROR_INT("no return val requested", procName, 1);
    if (parray) *parray = NULL;
    if (pcmap) *pcmap = NULL;
    if (!pixs || pixGetDepth(pixs) != 32)
        return ERROR_INT("pixs not defined", procName, 1);
    if (sigbits < 2 || sigbits > 6)
        return ERROR_INT("sigbits not in [2 ... 6]", procName, 1);
    if (factor < 1 || ncolors < 1)
        return ERROR_INT("factor < 1 or ncolors < 1", procName, 1);

    if ((nahisto = pixGetRGBHistogram(pixs, sigbits, factor)) == NULL)
        return ERROR_INT("nahisto not made", procName, 1);

        /* naindex contains the index into nahisto, which is the rgbindex */
    naindex = numaSortIndexAutoSelect(nahisto, L_SORT_DECREASING);
    numaDestroy(&nahisto);
    if (!naindex)
        return ERROR_INT("naindex not made", procName, 1);

    n = numaGetCount(naindex);
    ncolors = L_MIN(n, ncolors);
    if (parray) *parray = (l_uint32 *)LEPT_CALLOC(ncolors, sizeof(l_uint32));
    if (pcmap) *pcmap = pixcmapCreate(8);
    for (i = 0; i < ncolors; i++) {
        numaGetIValue(naindex, i, &rgbindex);  /* rgb index */
        getRGBFromIndex(rgbindex, sigbits, &rval, &gval, &bval);
        if (parray) composeRGBPixel(rval, gval, bval, *parray + i);
        if (pcmap) pixcmapAddColor(*pcmap, rval, gval, bval);
    }

    numaDestroy(&naindex);
    return 0;
}


PIX *
pixSimpleColorQuantize(PIX        *pixs,
                       l_int32     sigbits,
                       l_int32     factor,
                       l_int32     ncolors)
{
l_int32   w, h;
PIX      *pixd;
PIXCMAP  *cmap;

    PROCNAME("pixSimpleColorQuantize");

    if (!pixs || pixGetDepth(pixs) != 32)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (sigbits < 2 || sigbits > 4)
        return (PIX *)ERROR_PTR("sigbits not in {2,3,4}", procName, NULL);

    pixGetMostPopulatedColors(pixs, sigbits, factor, ncolors, NULL, &cmap);
    pixGetDimensions(pixs, &w, &h, NULL);
    pixd = pixCreate(w, h, 8);
    pixSetColormap(pixd, cmap);
    pixAssignToNearestColor(pixd, pixs, NULL, 4, NULL);
    return pixd;
}


/* ----------------------------------------------------------------------- *
 *            Constructs a color histogram based on rgb indices            *
 * ----------------------------------------------------------------------- */
NUMA *
pixGetRGBHistogram(PIX     *pixs,
                   l_int32  sigbits,
                   l_int32  factor)
{
l_int32     w, h, i, j, size, wpl, rval, gval, bval, npts;
l_uint32    val32, rgbindex;
l_float32  *array;
l_uint32   *data, *line, *rtab, *gtab, *btab;
NUMA       *na;

    PROCNAME("pixGetRGBHistogram");

    if (!pixs || pixGetDepth(pixs) != 32)
        return (NUMA *)ERROR_PTR("pixs not defined", procName, NULL);
    if (sigbits < 2 || sigbits > 6)
        return (NUMA *)ERROR_PTR("sigbits not in [2 ... 6]", procName, NULL);
    if (factor < 1)
        return (NUMA *)ERROR_PTR("factor < 1", procName, NULL);

        /* Get histogram size: 2^(3 * sigbits) */
    size = 1 << (3 * sigbits);  /* 64, 512, 4096, 32768, 262144 */
    na = numaMakeConstant(0, size);  /* init to all 0 */
    array = numaGetFArray(na, L_NOCOPY);

    makeRGBIndexTables(&rtab, &gtab, &btab, sigbits);

        /* Check the number of sampled pixels */
    pixGetDimensions(pixs, &w, &h, NULL);
    npts = ((w + factor - 1) / factor) * ((h + factor - 1) / factor);
    if (npts < 1000)
        L_WARNING("only sampling %d pixels\n", procName, npts);
    wpl = pixGetWpl(pixs);
    data = pixGetData(pixs);
    for (i = 0; i < h; i += factor) {
        line = data + i * wpl;
        for (j = 0; j < w; j += factor) {
            val32 = *(line + j);
            extractRGBValues(val32, &rval, &gval, &bval);
            rgbindex = rtab[rval] | gtab[gval] | btab[bval];
            array[rgbindex]++;
        }
    }

    LEPT_FREE(rtab);
    LEPT_FREE(gtab);
    LEPT_FREE(btab);
    return na;
}


l_ok
makeRGBIndexTables(l_uint32  **prtab,
                   l_uint32  **pgtab,
                   l_uint32  **pbtab,
                   l_int32     sigbits)
{
l_int32    i;
l_uint32  *rtab, *gtab, *btab;

    PROCNAME("makeRGBIndexTables");

    if (prtab) *prtab = NULL;
    if (pgtab) *pgtab = NULL;
    if (pbtab) *pbtab = NULL;
    if (!prtab || !pgtab || !pbtab)
        return ERROR_INT("not all table ptrs defined", procName, 1);
    if (sigbits < 2 || sigbits > 6)
        return ERROR_INT("sigbits not in [2 ... 6]", procName, 1);

    rtab = (l_uint32 *)LEPT_CALLOC(256, sizeof(l_uint32));
    gtab = (l_uint32 *)LEPT_CALLOC(256, sizeof(l_uint32));
    btab = (l_uint32 *)LEPT_CALLOC(256, sizeof(l_uint32));
    if (!rtab || !gtab || !btab)
        return ERROR_INT("calloc fail for tab", procName, 1);
    *prtab = rtab;
    *pgtab = gtab;
    *pbtab = btab;
    switch (sigbits) {
    case 2:
        for (i = 0; i < 256; i++) {
            rtab[i] = (i & 0xc0) >> 2;
            gtab[i] = (i & 0xc0) >> 4;
            btab[i] = (i & 0xc0) >> 6;
        }
        break;
    case 3:
        for (i = 0; i < 256; i++) {
            rtab[i] = (i & 0xe0) << 1;
            gtab[i] = (i & 0xe0) >> 2;
            btab[i] = (i & 0xe0) >> 5;
        }
        break;
    case 4:
        for (i = 0; i < 256; i++) {
            rtab[i] = (i & 0xf0) << 4;
            gtab[i] = (i & 0xf0);
            btab[i] = (i & 0xf0) >> 4;
        }
        break;
    case 5:
        for (i = 0; i < 256; i++) {
          rtab[i] = (i & 0xf8) << 7;
          gtab[i] = (i & 0xf8) << 2;
          btab[i] = (i & 0xf8) >> 3;
        }
        break;
    case 6:
        for (i = 0; i < 256; i++) {
          rtab[i] = (i & 0xfc) << 10;
          gtab[i] = (i & 0xfc) << 4;
          btab[i] = (i & 0xfc) >> 2;
        }
        break;
    default:
        L_ERROR("Illegal sigbits = %d\n", procName, sigbits);
        return ERROR_INT("sigbits not in [2 ... 6]", procName, 1);
    }

    return 0;
}


l_ok
getRGBFromIndex(l_uint32  index,
                l_int32   sigbits,
                l_int32  *prval,
                l_int32  *pgval,
                l_int32  *pbval)
{
    PROCNAME("getRGBFromIndex");

    if (prval) *prval = 0;
    if (pgval) *pgval = 0;
    if (pbval) *pbval = 0;
    if (!prval || !pgval || !pbval)
        return ERROR_INT("not all component ptrs defined", procName, 1);
    if (sigbits < 2 || sigbits > 6)
        return ERROR_INT("sigbits not in [2 ... 6]", procName, 1);

    switch (sigbits) {
    case 2:
        *prval = ((index << 2) & 0xc0) | 0x20;
        *pgval = ((index << 4) & 0xc0) | 0x20;
        *pbval = ((index << 6) & 0xc0) | 0x20;
        break;
    case 3:
        *prval = ((index >> 1) & 0xe0) | 0x10;
        *pgval = ((index << 2) & 0xe0) | 0x10;
        *pbval = ((index << 5) & 0xe0) | 0x10;
        break;
    case 4:
        *prval = ((index >> 4) & 0xf0) | 0x08;
        *pgval = (index & 0xf0) | 0x08;
        *pbval = ((index << 4) & 0xf0) | 0x08;
        break;
    case 5:
        *prval = ((index >> 7) & 0xf8) | 0x04;
        *pgval = ((index >> 2) & 0xf8) | 0x04;
        *pbval = ((index << 3) & 0xf8) | 0x04;
        break;
    case 6:
        *prval = ((index >> 10) & 0xfc) | 0x02;
        *pgval = ((index >> 4) & 0xfc) | 0x02;
        *pbval = ((index << 2) & 0xfc) | 0x02;
        break;
    default:
        L_ERROR("Illegal sigbits = %d\n", procName, sigbits);
        return ERROR_INT("sigbits not in [2 ... 6]", procName, 1);
    }

    return 0;
}


/* ----------------------------------------------------------------------- *
 *             Identify images that have highlight (red) color             *
 * ----------------------------------------------------------------------- */
l_ok
pixHasHighlightRed(PIX        *pixs,
                   l_int32     factor,
                   l_float32   fract,
                   l_float32   fthresh,
                   l_int32    *phasred,
                   l_float32  *pratio,
                   PIX       **ppixdb)
{
l_int32    w, h, count;
l_float32  ratio;
PIX       *pix1, *pix2, *pix3, *pix4;
FPIX      *fpix;

    PROCNAME("pixHasHighlightRed");

    if (pratio) *pratio = 0.0;
    if (ppixdb) *ppixdb = NULL;
    if (phasred) *phasred = 0;
    if (!pratio && !ppixdb)
        return ERROR_INT("no return val requested", procName, 1);
    if (!phasred)
        return ERROR_INT("&hasred not defined", procName, 1);
    if (!pixs || pixGetDepth(pixs) != 32)
        return ERROR_INT("pixs not defined or not 32 bpp", procName, 1);
    if (fthresh < 1.5 || fthresh > 3.5)
        L_WARNING("fthresh = %f is out of normal bounds\n", procName, fthresh);

    if (factor > 1)
        pix1 = pixScaleByIntSampling(pixs, factor);
    else
        pix1 = pixClone(pixs);

        /* Identify pixels that are either red or dark foreground */
    fpix = pixComponentFunction(pix1, 1.0, 0.0, -1.0, 0.0, 0.0, 1.0);
    pix2 = fpixThresholdToPix(fpix, fthresh);
    pixInvert(pix2, pix2);

        /* Identify pixels that are either red or light background */
    pix3 = pixGetRGBComponent(pix1, COLOR_RED);
    pix4 = pixThresholdToBinary(pix3, 130);
    pixInvert(pix4, pix4);

    pixAnd(pix4, pix4, pix2);
    pixCountPixels(pix4, &count, NULL);
    pixGetDimensions(pix4, &w, &h, NULL);
    L_INFO("count = %d, thresh = %d\n", procName, count,
           (l_int32)(fract * w * h));
    ratio = (l_float32)count / (fract * w * h);
    if (pratio) *pratio = ratio;
    if (ratio >= 1.0)
        *phasred = 1;
    if (ppixdb)
        *ppixdb = pix4;
    else
        pixDestroy(&pix4);
    pixDestroy(&pix1);
    pixDestroy(&pix2);
    pixDestroy(&pix3);
    fpixDestroy(&fpix);
    return 0;
}