grayquant.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"

static void ditherToBinaryLow(l_uint32 *datad, l_int32 w, l_int32 h,
                              l_int32 wpld, l_uint32 *datas, l_int32 wpls,
                              l_uint32 *bufs1, l_uint32 *bufs2,
                              l_int32 lowerclip, l_int32 upperclip);
static void thresholdToBinaryLow(l_uint32 *datad, l_int32 w, l_int32 h,
                                 l_int32 wpld, l_uint32 *datas, l_int32 d,
                                 l_int32 wpls, l_int32 thresh);
static void ditherTo2bppLow(l_uint32 *datad, l_int32 w, l_int32 h, l_int32 wpld,
                            l_uint32 *datas, l_int32 wpls, l_uint32 *bufs1,
                            l_uint32 *bufs2, l_int32 *tabval, l_int32 *tab38,
                            l_int32   *tab14);
static void ditherTo2bppLineLow(l_uint32 *lined, l_int32 w, l_uint32 *bufs1,
                                l_uint32 *bufs2, l_int32 *tabval,
                                l_int32 *tab38, l_int32 *tab14,
                                l_int32 lastlineflag);
static l_int32 make8To2DitherTables(l_int32 **ptabval, l_int32 **ptab38,
                                    l_int32 **ptab14, l_int32 cliptoblack,
                                    l_int32 cliptowhite);
static void thresholdTo2bppLow(l_uint32 *datad, l_int32 h, l_int32 wpld,
                               l_uint32 *datas, l_int32 wpls, l_int32 *tab);
static void thresholdTo4bppLow(l_uint32 *datad, l_int32 h, l_int32 wpld,
                               l_uint32 *datas, l_int32 wpls, l_int32 *tab);
static l_int32 *makeGrayQuantTargetTable(l_int32 nlevels, l_int32 depth);
static l_int32 makeGrayQuantColormapArb(PIX *pixs, l_int32 *tab,
                                        l_int32 outdepth, PIXCMAP **pcmap);
static l_int32 numaFillCmapFromHisto(NUMA *na, PIXCMAP *cmap,
                                     l_float32 minfract, l_int32 maxsize,
                                     l_int32 **plut);

#ifndef  NO_CONSOLE_IO
#define DEBUG_UNROLLING 0
#endif   /* ~NO_CONSOLE_IO */

/*------------------------------------------------------------------*
 *             Binarization by Floyd-Steinberg dithering            *
 *------------------------------------------------------------------*/
PIX *
pixDitherToBinary(PIX  *pixs)
{
    PROCNAME("pixDitherToBinary");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (pixGetDepth(pixs) != 8)
        return (PIX *)ERROR_PTR("must be 8 bpp for dithering", procName, NULL);

    return pixDitherToBinarySpec(pixs, DEFAULT_CLIP_LOWER_1,
                                 DEFAULT_CLIP_UPPER_1);
}


PIX *
pixDitherToBinarySpec(PIX     *pixs,
                      l_int32  lowerclip,
                      l_int32  upperclip)
{
l_int32    w, h, d, wplt, wpld;
l_uint32  *datat, *datad;
l_uint32  *bufs1, *bufs2;
PIX       *pixt, *pixd;

    PROCNAME("pixDitherToBinarySpec");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 8)
        return (PIX *)ERROR_PTR("must be 8 bpp for dithering", procName, NULL);
    if (lowerclip < 0 || lowerclip > 255)
        return (PIX *)ERROR_PTR("invalid value for lowerclip", procName, NULL);
    if (upperclip < 0 || upperclip > 255)
        return (PIX *)ERROR_PTR("invalid value for upperclip", procName, NULL);

    if ((pixd = pixCreate(w, h, 1)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    pixCopyResolution(pixd, pixs);
    pixCopyInputFormat(pixd, pixs);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);

        /* Remove colormap if it exists */
    if ((pixt = pixRemoveColormap(pixs, REMOVE_CMAP_TO_GRAYSCALE)) == NULL) {
        pixDestroy(&pixd);
        return (PIX *)ERROR_PTR("pixt not made", procName, NULL);
    }
    datat = pixGetData(pixt);
    wplt = pixGetWpl(pixt);

        /* Two line buffers, 1 for current line and 2 for next line */
    bufs1 = (l_uint32 *)LEPT_CALLOC(wplt, sizeof(l_uint32));
    bufs2 = (l_uint32 *)LEPT_CALLOC(wplt, sizeof(l_uint32));
    if (!bufs1 || !bufs2) {
        LEPT_FREE(bufs1);
        LEPT_FREE(bufs2);
        pixDestroy(&pixd);
        pixDestroy(&pixt);
        return (PIX *)ERROR_PTR("bufs1, bufs2 not both made", procName, NULL);
    }

    ditherToBinaryLow(datad, w, h, wpld, datat, wplt, bufs1, bufs2,
                      lowerclip, upperclip);

    LEPT_FREE(bufs1);
    LEPT_FREE(bufs2);
    pixDestroy(&pixt);
    return pixd;
}


static void
ditherToBinaryLow(l_uint32  *datad,
                  l_int32    w,
                  l_int32    h,
                  l_int32    wpld,
                  l_uint32  *datas,
                  l_int32    wpls,
                  l_uint32  *bufs1,
                  l_uint32  *bufs2,
                  l_int32    lowerclip,
                  l_int32    upperclip)
{
l_int32    i;
l_uint32  *lined;

        /* do all lines except last line */
    memcpy(bufs2, datas, 4 * wpls);  /* prime the buffer */
    for (i = 0; i < h - 1; i++) {
        memcpy(bufs1, bufs2, 4 * wpls);
        memcpy(bufs2, datas + (i + 1) * wpls, 4 * wpls);
        lined = datad + i * wpld;
        ditherToBinaryLineLow(lined, w, bufs1, bufs2, lowerclip, upperclip, 0);
    }

        /* do last line */
    memcpy(bufs1, bufs2, 4 * wpls);
    lined = datad + (h - 1) * wpld;
    ditherToBinaryLineLow(lined, w, bufs1, bufs2, lowerclip, upperclip, 1);
}


void
ditherToBinaryLineLow(l_uint32  *lined,
                      l_int32    w,
                      l_uint32  *bufs1,
                      l_uint32  *bufs2,
                      l_int32    lowerclip,
                      l_int32    upperclip,
                      l_int32    lastlineflag)
{
l_int32   j;
l_int32   oval, eval;
l_uint8   fval1, fval2, rval, bval, dval;

    if (lastlineflag == 0) {
        for (j = 0; j < w - 1; j++) {
            oval = GET_DATA_BYTE(bufs1, j);
            if (oval > 127) {   /* binarize to OFF */
                if ((eval = 255 - oval) > upperclip) {
                        /* subtract from neighbors */
                    fval1 = (3 * eval) / 8;
                    fval2 = eval / 4;
                    rval = GET_DATA_BYTE(bufs1, j + 1);
                    rval = L_MAX(0, rval - fval1);
                    SET_DATA_BYTE(bufs1, j + 1, rval);
                    bval = GET_DATA_BYTE(bufs2, j);
                    bval = L_MAX(0, bval - fval1);
                    SET_DATA_BYTE(bufs2, j, bval);
                    dval = GET_DATA_BYTE(bufs2, j + 1);
                    dval = L_MAX(0, dval - fval2);
                    SET_DATA_BYTE(bufs2, j + 1, dval);
                }
            } else {   /* oval <= 127; binarize to ON  */
                SET_DATA_BIT(lined, j);   /* ON pixel */
                if (oval > lowerclip) {
                        /* add to neighbors */
                    fval1 = (3 * oval) / 8;
                    fval2 = oval / 4;
                    rval = GET_DATA_BYTE(bufs1, j + 1);
                    rval = L_MIN(255, rval + fval1);
                    SET_DATA_BYTE(bufs1, j + 1, rval);
                    bval = GET_DATA_BYTE(bufs2, j);
                    bval = L_MIN(255, bval + fval1);
                    SET_DATA_BYTE(bufs2, j, bval);
                    dval = GET_DATA_BYTE(bufs2, j + 1);
                    dval = L_MIN(255, dval + fval2);
                    SET_DATA_BYTE(bufs2, j + 1, dval);
                }
            }
        }

            /* do last column: j = w - 1 */
        oval = GET_DATA_BYTE(bufs1, j);
        if (oval > 127) {  /* binarize to OFF */
            if ((eval = 255 - oval) > upperclip) {
                    /* subtract from neighbors */
                fval1 = (3 * eval) / 8;
                bval = GET_DATA_BYTE(bufs2, j);
                bval = L_MAX(0, bval - fval1);
                SET_DATA_BYTE(bufs2, j, bval);
            }
        } else {  /*oval <= 127; binarize to ON */
            SET_DATA_BIT(lined, j);   /* ON pixel */
            if (oval > lowerclip) {
                    /* add to neighbors */
                fval1 = (3 * oval) / 8;
                bval = GET_DATA_BYTE(bufs2, j);
                bval = L_MIN(255, bval + fval1);
                SET_DATA_BYTE(bufs2, j, bval);
            }
        }
    } else {   /* lastlineflag == 1 */
        for (j = 0; j < w - 1; j++) {
            oval = GET_DATA_BYTE(bufs1, j);
            if (oval > 127) {   /* binarize to OFF */
                if ((eval = 255 - oval) > upperclip) {
                        /* subtract from neighbors */
                    fval1 = (3 * eval) / 8;
                    rval = GET_DATA_BYTE(bufs1, j + 1);
                    rval = L_MAX(0, rval - fval1);
                    SET_DATA_BYTE(bufs1, j + 1, rval);
                }
            } else {   /* oval <= 127; binarize to ON  */
                SET_DATA_BIT(lined, j);   /* ON pixel */
                if (oval > lowerclip) {
                        /* add to neighbors */
                    fval1 = (3 * oval) / 8;
                    rval = GET_DATA_BYTE(bufs1, j + 1);
                    rval = L_MIN(255, rval + fval1);
                    SET_DATA_BYTE(bufs1, j + 1, rval);
                }
            }
        }

            /* do last pixel: (i, j) = (h - 1, w - 1) */
        oval = GET_DATA_BYTE(bufs1, j);
        if (oval < 128)
            SET_DATA_BIT(lined, j);   /* ON pixel */
    }
}


/*------------------------------------------------------------------*
 *       Simple (pixelwise) binarization with fixed threshold       *
 *------------------------------------------------------------------*/
PIX *
pixThresholdToBinary(PIX     *pixs,
                     l_int32  thresh)
{
l_int32    d, w, h, wplt, wpld;
l_uint32  *datat, *datad;
PIX       *pixt, *pixd;

    PROCNAME("pixThresholdToBinary");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 4 && d != 8)
        return (PIX *)ERROR_PTR("pixs must be 4 or 8 bpp", procName, NULL);
    if (thresh < 0)
        return (PIX *)ERROR_PTR("thresh must be non-negative", procName, NULL);
    if (d == 4 && thresh > 16)
        return (PIX *)ERROR_PTR("4 bpp thresh not in {0-16}", procName, NULL);
    if (d == 8 && thresh > 256)
        return (PIX *)ERROR_PTR("8 bpp thresh not in {0-256}", procName, NULL);

    if ((pixd = pixCreate(w, h, 1)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    pixCopyResolution(pixd, pixs);
    pixCopyInputFormat(pixd, pixs);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);

        /* Remove colormap if it exists.  If there is a colormap,
         * pixt will be 8 bpp regardless of the depth of pixs. */
    pixt = pixRemoveColormap(pixs, REMOVE_CMAP_TO_GRAYSCALE);
    datat = pixGetData(pixt);
    wplt = pixGetWpl(pixt);
    if (pixGetColormap(pixs) && d == 4) {  /* promoted to 8 bpp */
        d = 8;
        thresh *= 16;
    }

    thresholdToBinaryLow(datad, w, h, wpld, datat, d, wplt, thresh);
    pixDestroy(&pixt);
    return pixd;
}


static void
thresholdToBinaryLow(l_uint32  *datad,
                     l_int32    w,
                     l_int32    h,
                     l_int32    wpld,
                     l_uint32  *datas,
                     l_int32    d,
                     l_int32    wpls,
                     l_int32    thresh)
{
l_int32    i;
l_uint32  *lines, *lined;

    for (i = 0; i < h; i++) {
        lines = datas + i * wpls;
        lined = datad + i * wpld;
        thresholdToBinaryLineLow(lined, w, lines, d, thresh);
    }
}


/*
 *  thresholdToBinaryLineLow()
 *
 */
void
thresholdToBinaryLineLow(l_uint32  *lined,
                         l_int32    w,
                         l_uint32  *lines,
                         l_int32    d,
                         l_int32    thresh)
{
l_int32  j, k, gval, scount, dcount;
l_uint32 sword, dword;

    PROCNAME("thresholdToBinaryLineLow");

    switch (d)
    {
    case 4:
            /* Unrolled as 4 source words, 1 dest word */
        for (j = 0, scount = 0, dcount = 0; j + 31 < w; j += 32) {
            dword = 0;
            for (k = 0; k < 4; k++) {
                sword = lines[scount++];
                dword <<= 8;
                gval = (sword >> 28) & 0xf;
                    /* Trick used here and below: if gval < thresh then
                     * gval - thresh < 0, so its high-order bit is 1, and
                     * ((gval - thresh) >> 31) & 1 == 1; likewise, if
                     * gval >= thresh, then ((gval - thresh) >> 31) & 1 == 0
                     * Doing it this way avoids a random (and thus easily
                     * mispredicted) branch on each pixel. */
                dword |= ((gval - thresh) >> 24) & 128;
                gval = (sword >> 24) & 0xf;
                dword |= ((gval - thresh) >> 25) & 64;
                gval = (sword >> 20) & 0xf;
                dword |= ((gval - thresh) >> 26) & 32;
                gval = (sword >> 16) & 0xf;
                dword |= ((gval - thresh) >> 27) & 16;
                gval = (sword >> 12) & 0xf;
                dword |= ((gval - thresh) >> 28) & 8;
                gval = (sword >> 8) & 0xf;
                dword |= ((gval - thresh) >> 29) & 4;
                gval = (sword >> 4) & 0xf;
                dword |= ((gval - thresh) >> 30) & 2;
                gval = sword & 0xf;
                dword |= ((gval - thresh) >> 31) & 1;
            }
            lined[dcount++] = dword;
        }

        if (j < w) {
          dword = 0;
          for (; j < w; j++) {
              if ((j & 7) == 0) {
                  sword = lines[scount++];
              }
              gval = (sword >> 28) & 0xf;
              sword <<= 4;
              dword |= (((gval - thresh) >> 31) & 1) << (31 - (j & 31));
          }
          lined[dcount] = dword;
        }
#if DEBUG_UNROLLING
#define CHECK_BIT(a, b, c) if (GET_DATA_BIT(a, b) != c) { \
    fprintf(stderr, "Error: mismatch at %d/%d(%d), %d vs %d\n", \
            j, w, d, GET_DATA_BIT(a, b), c); }
        for (j = 0; j < w; j++) {
            gval = GET_DATA_QBIT(lines, j);
            CHECK_BIT(lined, j, gval < thresh ? 1 : 0);
        }
#endif
        break;
    case 8:
            /* Unrolled as 8 source words, 1 dest word */
        for (j = 0, scount = 0, dcount = 0; j + 31 < w; j += 32) {
            dword = 0;
            for (k = 0; k < 8; k++) {
                sword = lines[scount++];
                dword <<= 4;
                gval = (sword >> 24) & 0xff;
                dword |= ((gval - thresh) >> 28) & 8;
                gval = (sword >> 16) & 0xff;
                dword |= ((gval - thresh) >> 29) & 4;
                gval = (sword >> 8) & 0xff;
                dword |= ((gval - thresh) >> 30) & 2;
                gval = sword & 0xff;
                dword |= ((gval - thresh) >> 31) & 1;
            }
            lined[dcount++] = dword;
        }

        if (j < w) {
            dword = 0;
            for (; j < w; j++) {
                if ((j & 3) == 0) {
                    sword = lines[scount++];
                }
                gval = (sword >> 24) & 0xff;
                sword <<= 8;
                dword |= (l_uint64)(((gval - thresh) >> 31) & 1)
                             << (31 - (j & 31));
            }
            lined[dcount] = dword;
        }
#if DEBUG_UNROLLING
        for (j = 0; j < w; j++) {
            gval = GET_DATA_BYTE(lines, j);
            CHECK_BIT(lined, j, gval < thresh ? 1 : 0);
        }
#undef CHECK_BIT
#endif
        break;
    default:
        L_ERROR("src depth not 4 or 8 bpp\n", procName);
        break;
    }
}


/*------------------------------------------------------------------*
 *                Binarization with variable threshold              *
 *------------------------------------------------------------------*/
PIX *
pixVarThresholdToBinary(PIX  *pixs,
                        PIX  *pixg)
{
l_int32    i, j, vals, valg, w, h, d, wpls, wplg, wpld;
l_uint32  *datas, *datag, *datad, *lines, *lineg, *lined;
PIX       *pixd;

    PROCNAME("pixVarThresholdToBinary");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (!pixg)
        return (PIX *)ERROR_PTR("pixg not defined", procName, NULL);
    if (!pixSizesEqual(pixs, pixg))
        return (PIX *)ERROR_PTR("pix sizes not equal", procName, NULL);
    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 8)
        return (PIX *)ERROR_PTR("pixs must be 8 bpp", procName, NULL);

    pixd = pixCreate(w, h, 1);
    pixCopyResolution(pixd, pixs);
    pixCopyInputFormat(pixd, pixs);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);
    datas = pixGetData(pixs);
    wpls = pixGetWpl(pixs);
    datag = pixGetData(pixg);
    wplg = pixGetWpl(pixg);
    for (i = 0; i < h; i++) {
        lines = datas + i * wpls;
        lineg = datag + i * wplg;
        lined = datad + i * wpld;
        for (j = 0; j < w; j++) {
            vals = GET_DATA_BYTE(lines, j);
            valg = GET_DATA_BYTE(lineg, j);
            if (vals < valg)
                SET_DATA_BIT(lined, j);
        }
    }

    return pixd;
}


/*------------------------------------------------------------------*
 *                  Binarization by adaptive mapping                *
 *------------------------------------------------------------------*/
PIX *
pixAdaptThresholdToBinary(PIX       *pixs,
                          PIX       *pixm,
                          l_float32  gamma)
{
    PROCNAME("pixAdaptThresholdToBinary");

    if (!pixs || pixGetDepth(pixs) != 8)
        return (PIX *)ERROR_PTR("pixs undefined or not 8 bpp", procName, NULL);

    return pixAdaptThresholdToBinaryGen(pixs, pixm, gamma, 50, 170, 200);
}


PIX *
pixAdaptThresholdToBinaryGen(PIX       *pixs,
                             PIX       *pixm,
                             l_float32  gamma,
                             l_int32    blackval,
                             l_int32    whiteval,
                             l_int32    thresh)
{
PIX  *pix1, *pixd;

    PROCNAME("pixAdaptThresholdToBinaryGen");

    if (!pixs || pixGetDepth(pixs) != 8)
        return (PIX *)ERROR_PTR("pixs undefined or not 8 bpp", procName, NULL);

    pix1 = pixBackgroundNormSimple(pixs, pixm, NULL);
    pixGammaTRC(pix1, pix1, gamma, blackval, whiteval);
    pixd = pixThresholdToBinary(pix1, thresh);
    pixDestroy(&pix1);
    return pixd;
}


/*--------------------------------------------------------------------*
 *       Generate a binary mask from pixels of particular value(s)    *
 *--------------------------------------------------------------------*/
PIX *
pixGenerateMaskByValue(PIX     *pixs,
                       l_int32  val,
                       l_int32  usecmap)
{
l_int32    i, j, w, h, d, wplg, wpld;
l_uint32  *datag, *datad, *lineg, *lined;
PIX       *pixg, *pixd;

    PROCNAME("pixGenerateMaskByValue");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    d = pixGetDepth(pixs);
    if (d != 2 && d != 4 && d != 8)
        return (PIX *)ERROR_PTR("not 2, 4 or 8 bpp", procName, NULL);

    if (!usecmap && pixGetColormap(pixs))
        pixg = pixRemoveColormap(pixs, REMOVE_CMAP_TO_GRAYSCALE);
    else
        pixg = pixClone(pixs);
    pixGetDimensions(pixg, &w, &h, &d);
    if (d == 8 && (val < 0 || val > 255)) {
        pixDestroy(&pixg);
        return (PIX *)ERROR_PTR("val out of 8 bpp range", procName, NULL);
    }
    if (d == 4 && (val < 0 || val > 15)) {
        pixDestroy(&pixg);
        return (PIX *)ERROR_PTR("val out of 4 bpp range", procName, NULL);
    }
    if (d == 2 && (val < 0 || val > 3)) {
        pixDestroy(&pixg);
        return (PIX *)ERROR_PTR("val out of 2 bpp range", procName, NULL);
    }

    pixd = pixCreate(w, h, 1);
    pixCopyResolution(pixd, pixg);
    pixCopyInputFormat(pixd, pixs);
    datag = pixGetData(pixg);
    wplg = pixGetWpl(pixg);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);
    for (i = 0; i < h; i++) {
        lineg = datag + i * wplg;
        lined = datad + i * wpld;
        for (j = 0; j < w; j++) {
            if (d == 8) {
                if (GET_DATA_BYTE(lineg, j) == val)
                    SET_DATA_BIT(lined, j);
            } else if (d == 4) {
                if (GET_DATA_QBIT(lineg, j) == val)
                    SET_DATA_BIT(lined, j);
            } else {  /* d == 2 */
                if (GET_DATA_DIBIT(lineg, j) == val)
                    SET_DATA_BIT(lined, j);
            }
        }
    }

    pixDestroy(&pixg);
    return pixd;
}


PIX *
pixGenerateMaskByBand(PIX     *pixs,
                      l_int32  lower,
                      l_int32  upper,
                      l_int32  inband,
                      l_int32  usecmap)
{
l_int32    i, j, w, h, d, wplg, wpld, val;
l_uint32  *datag, *datad, *lineg, *lined;
PIX       *pixg, *pixd;

    PROCNAME("pixGenerateMaskByBand");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    d = pixGetDepth(pixs);
    if (d != 2 && d != 4 && d != 8)
        return (PIX *)ERROR_PTR("not 2, 4 or 8 bpp", procName, NULL);
    if (lower < 0 || lower > upper)
        return (PIX *)ERROR_PTR("lower < 0 or lower > upper!", procName, NULL);

    if (!usecmap && pixGetColormap(pixs))
        pixg = pixRemoveColormap(pixs, REMOVE_CMAP_TO_GRAYSCALE);
    else
        pixg = pixClone(pixs);
    pixGetDimensions(pixg, &w, &h, &d);
    if (d == 8 && upper > 255) {
        pixDestroy(&pixg);
        return (PIX *)ERROR_PTR("d == 8 and upper > 255", procName, NULL);
    }
    if (d == 4 && upper > 15) {
        pixDestroy(&pixg);
        return (PIX *)ERROR_PTR("d == 4 and upper > 15", procName, NULL);
    }
    if (d == 2 && upper > 3) {
        pixDestroy(&pixg);
        return (PIX *)ERROR_PTR("d == 2 and upper > 3", procName, NULL);
    }

    pixd = pixCreate(w, h, 1);
    pixCopyResolution(pixd, pixg);
    pixCopyInputFormat(pixd, pixs);
    datag = pixGetData(pixg);
    wplg = pixGetWpl(pixg);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);
    for (i = 0; i < h; i++) {
        lineg = datag + i * wplg;
        lined = datad + i * wpld;
        for (j = 0; j < w; j++) {
            if (d == 8)
                val = GET_DATA_BYTE(lineg, j);
            else if (d == 4)
                val = GET_DATA_QBIT(lineg, j);
            else  /* d == 2 */
                val = GET_DATA_DIBIT(lineg, j);
            if (inband) {
                if (val >= lower && val <= upper)
                    SET_DATA_BIT(lined, j);
            } else {  /* out of band */
                if (val < lower || val > upper)
                    SET_DATA_BIT(lined, j);
            }
        }
    }

    pixDestroy(&pixg);
    return pixd;
}


/*------------------------------------------------------------------*
 *                Thresholding to 2 bpp by dithering                *
 *------------------------------------------------------------------*/
PIX *
pixDitherTo2bpp(PIX     *pixs,
                l_int32  cmapflag)
{
    PROCNAME("pixDitherTo2bpp");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (pixGetDepth(pixs) != 8)
        return (PIX *)ERROR_PTR("must be 8 bpp for dithering", procName, NULL);

    return pixDitherTo2bppSpec(pixs, DEFAULT_CLIP_LOWER_2,
                               DEFAULT_CLIP_UPPER_2, cmapflag);
}


PIX *
pixDitherTo2bppSpec(PIX     *pixs,
                    l_int32  lowerclip,
                    l_int32  upperclip,
                    l_int32  cmapflag)
{
l_int32    w, h, d, wplt, wpld;
l_int32   *tabval, *tab38, *tab14;
l_uint32  *datat, *datad;
l_uint32  *bufs1, *bufs2;
PIX       *pixt, *pixd;
PIXCMAP   *cmap;

    PROCNAME("pixDitherTo2bppSpec");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 8)
        return (PIX *)ERROR_PTR("must be 8 bpp for dithering", procName, NULL);
    if (lowerclip < 0 || lowerclip > 255)
        return (PIX *)ERROR_PTR("invalid value for lowerclip", procName, NULL);
    if (upperclip < 0 || upperclip > 255)
        return (PIX *)ERROR_PTR("invalid value for upperclip", procName, NULL);

    if ((pixd = pixCreate(w, h, 2)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    pixCopyResolution(pixd, pixs);
    pixCopyInputFormat(pixd, pixs);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);

        /* If there is a colormap, remove it */
    pixt = pixRemoveColormap(pixs, REMOVE_CMAP_TO_GRAYSCALE);
    datat = pixGetData(pixt);
    wplt = pixGetWpl(pixt);

        /* Two line buffers, 1 for current line and 2 for next line */
    bufs1 = (l_uint32 *)LEPT_CALLOC(wplt, sizeof(l_uint32));
    bufs2 = (l_uint32 *)LEPT_CALLOC(wplt, sizeof(l_uint32));
    if (!bufs1 || !bufs2) {
        LEPT_FREE(bufs1);
        LEPT_FREE(bufs2);
        pixDestroy(&pixd);
        pixDestroy(&pixt);
        return (PIX *)ERROR_PTR("bufs1, bufs2 not both made", procName, NULL);
    }

        /* 3 lookup tables: 2-bit value, (3/8)excess, and (1/4)excess */
    make8To2DitherTables(&tabval, &tab38, &tab14, lowerclip, upperclip);

    ditherTo2bppLow(datad, w, h, wpld, datat, wplt, bufs1, bufs2,
                    tabval, tab38, tab14);

    if (cmapflag) {
        cmap = pixcmapCreateLinear(2, 4);
        pixSetColormap(pixd, cmap);
    }

    LEPT_FREE(bufs1);
    LEPT_FREE(bufs2);
    LEPT_FREE(tabval);
    LEPT_FREE(tab38);
    LEPT_FREE(tab14);
    pixDestroy(&pixt);
    return pixd;
}


static void
ditherTo2bppLow(l_uint32  *datad,
                l_int32    w,
                l_int32    h,
                l_int32    wpld,
                l_uint32  *datas,
                l_int32    wpls,
                l_uint32  *bufs1,
                l_uint32  *bufs2,
                l_int32   *tabval,
                l_int32   *tab38,
                l_int32   *tab14)
{
l_int32      i;
l_uint32    *lined;

        /* do all lines except last line */
    memcpy(bufs2, datas, 4 * wpls);  /* prime the buffer */
    for (i = 0; i < h - 1; i++) {
        memcpy(bufs1, bufs2, 4 * wpls);
        memcpy(bufs2, datas + (i + 1) * wpls, 4 * wpls);
        lined = datad + i * wpld;
        ditherTo2bppLineLow(lined, w, bufs1, bufs2, tabval, tab38, tab14, 0);
    }

        /* do last line */
    memcpy(bufs1, bufs2, 4 * wpls);
    lined = datad + (h - 1) * wpld;
    ditherTo2bppLineLow(lined, w, bufs1, bufs2, tabval, tab38, tab14, 1);
}


static void
ditherTo2bppLineLow(l_uint32  *lined,
                    l_int32    w,
                    l_uint32  *bufs1,
                    l_uint32  *bufs2,
                    l_int32   *tabval,
                    l_int32   *tab38,
                    l_int32   *tab14,
                    l_int32    lastlineflag)
{
l_int32  j;
l_int32  oval, tab38val, tab14val;
l_uint8  rval, bval, dval;

    if (lastlineflag == 0) {
        for (j = 0; j < w - 1; j++) {
            oval = GET_DATA_BYTE(bufs1, j);
            SET_DATA_DIBIT(lined, j, tabval[oval]);
            rval = GET_DATA_BYTE(bufs1, j + 1);
            bval = GET_DATA_BYTE(bufs2, j);
            dval = GET_DATA_BYTE(bufs2, j + 1);
            tab38val = tab38[oval];
            tab14val = tab14[oval];
            if (tab38val < 0) {
                rval = L_MAX(0, rval + tab38val);
                bval = L_MAX(0, bval + tab38val);
                dval = L_MAX(0, dval + tab14val);
            } else {
                rval = L_MIN(255, rval + tab38val);
                bval = L_MIN(255, bval + tab38val);
                dval = L_MIN(255, dval + tab14val);
            }
            SET_DATA_BYTE(bufs1, j + 1, rval);
            SET_DATA_BYTE(bufs2, j, bval);
            SET_DATA_BYTE(bufs2, j + 1, dval);
        }

            /* do last column: j = w - 1 */
        oval = GET_DATA_BYTE(bufs1, j);
        SET_DATA_DIBIT(lined, j, tabval[oval]);
        bval = GET_DATA_BYTE(bufs2, j);
        tab38val = tab38[oval];
        if (tab38val < 0)
            bval = L_MAX(0, bval + tab38val);
        else
            bval = L_MIN(255, bval + tab38val);
        SET_DATA_BYTE(bufs2, j, bval);
    } else {   /* lastlineflag == 1 */
        for (j = 0; j < w - 1; j++) {
            oval = GET_DATA_BYTE(bufs1, j);
            SET_DATA_DIBIT(lined, j, tabval[oval]);
            rval = GET_DATA_BYTE(bufs1, j + 1);
            tab38val = tab38[oval];
            if (tab38val < 0)
                rval = L_MAX(0, rval + tab38val);
            else
                rval = L_MIN(255, rval + tab38val);
            SET_DATA_BYTE(bufs1, j + 1, rval);
        }

            /* do last pixel: (i, j) = (h - 1, w - 1) */
        oval = GET_DATA_BYTE(bufs1, j);
        SET_DATA_DIBIT(lined, j, tabval[oval]);
    }
}


static l_int32
make8To2DitherTables(l_int32 **ptabval,
                     l_int32 **ptab38,
                     l_int32 **ptab14,
                     l_int32   cliptoblack,
                     l_int32   cliptowhite)
{
l_int32   i;
l_int32  *tabval, *tab38, *tab14;

        /* 3 lookup tables: 2-bit value, (3/8)excess, and (1/4)excess */
    tabval = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32));
    tab38 = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32));
    tab14 = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32));
    *ptabval = tabval;
    *ptab38 = tab38;
    *ptab14 = tab14;

    for (i = 0; i < 256; i++) {
        if (i <= cliptoblack) {
            tabval[i] = 0;
            tab38[i] = 0;
            tab14[i] = 0;
        } else if (i < 43) {
            tabval[i] = 0;
            tab38[i] = (3 * i + 4) / 8;
            tab14[i] = (i + 2) / 4;
        } else if (i < 85) {
            tabval[i] = 1;
            tab38[i] = (3 * (i - 85) - 4) / 8;
            tab14[i] = ((i - 85) - 2) / 4;
        } else if (i < 128) {
            tabval[i] = 1;
            tab38[i] = (3 * (i - 85) + 4) / 8;
            tab14[i] = ((i - 85) + 2) / 4;
        } else if (i < 170) {
            tabval[i] = 2;
            tab38[i] = (3 * (i - 170) - 4) / 8;
            tab14[i] = ((i - 170) - 2) / 4;
        } else if (i < 213) {
            tabval[i] = 2;
            tab38[i] = (3 * (i - 170) + 4) / 8;
            tab14[i] = ((i - 170) + 2) / 4;
        } else if (i < 255 - cliptowhite) {
            tabval[i] = 3;
            tab38[i] = (3 * (i - 255) - 4) / 8;
            tab14[i] = ((i - 255) - 2) / 4;
        } else {  /* i >= 255 - cliptowhite */
            tabval[i] = 3;
            tab38[i] = 0;
            tab14[i] = 0;
        }
    }

    return 0;
}


/*--------------------------------------------------------------------*
 *  Simple (pixelwise) thresholding to 2 bpp with optional colormap   *
 *--------------------------------------------------------------------*/
PIX *
pixThresholdTo2bpp(PIX     *pixs,
                   l_int32  nlevels,
                   l_int32  cmapflag)
{
l_int32   *qtab;
l_int32    w, h, d, wplt, wpld;
l_uint32  *datat, *datad;
PIX       *pixt, *pixd;
PIXCMAP   *cmap;

    PROCNAME("pixThresholdTo2bpp");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 8)
        return (PIX *)ERROR_PTR("pixs not 8 bpp", procName, NULL);
    if (nlevels < 2 || nlevels > 4)
        return (PIX *)ERROR_PTR("nlevels not in {2, 3, 4}", procName, NULL);

    if ((pixd = pixCreate(w, h, 2)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    pixCopyResolution(pixd, pixs);
    pixCopyInputFormat(pixd, pixs);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);

    if (cmapflag) {   /* hold out (4 - nlevels) cmap entries */
        cmap = pixcmapCreateLinear(2, nlevels);
        pixSetColormap(pixd, cmap);
    }

        /* If there is a colormap in the src, remove it */
    pixt = pixRemoveColormap(pixs, REMOVE_CMAP_TO_GRAYSCALE);
    datat = pixGetData(pixt);
    wplt = pixGetWpl(pixt);

        /* Make the appropriate table */
    if (cmapflag)
        qtab = makeGrayQuantIndexTable(nlevels);
    else
        qtab = makeGrayQuantTargetTable(4, 2);

    thresholdTo2bppLow(datad, h, wpld, datat, wplt, qtab);

    LEPT_FREE(qtab);
    pixDestroy(&pixt);
    return pixd;
}


static void
thresholdTo2bppLow(l_uint32  *datad,
                   l_int32    h,
                   l_int32    wpld,
                   l_uint32  *datas,
                   l_int32    wpls,
                   l_int32   *tab)
{
l_uint8    sval1, sval2, sval3, sval4, dval;
l_int32    i, j, k;
l_uint32  *lines, *lined;

    for (i = 0; i < h; i++) {
        lines = datas + i * wpls;
        lined = datad + i * wpld;
        for (j = 0; j < wpls; j++) {
            k = 4 * j;
            sval1 = GET_DATA_BYTE(lines, k);
            sval2 = GET_DATA_BYTE(lines, k + 1);
            sval3 = GET_DATA_BYTE(lines, k + 2);
            sval4 = GET_DATA_BYTE(lines, k + 3);
            dval = (tab[sval1] << 6) | (tab[sval2] << 4) |
                   (tab[sval3] << 2) | tab[sval4];
            SET_DATA_BYTE(lined, j, dval);
        }
    }
}


/*----------------------------------------------------------------------*
 *               Simple (pixelwise) thresholding to 4 bpp               *
 *----------------------------------------------------------------------*/
PIX *
pixThresholdTo4bpp(PIX     *pixs,
                   l_int32  nlevels,
                   l_int32  cmapflag)
{
l_int32   *qtab;
l_int32    w, h, d, wplt, wpld;
l_uint32  *datat, *datad;
PIX       *pixt, *pixd;
PIXCMAP   *cmap;

    PROCNAME("pixThresholdTo4bpp");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 8)
        return (PIX *)ERROR_PTR("pixs not 8 bpp", procName, NULL);
    if (nlevels < 2 || nlevels > 16)
        return (PIX *)ERROR_PTR("nlevels not in [2,...,16]", procName, NULL);

    if ((pixd = pixCreate(w, h, 4)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    pixCopyResolution(pixd, pixs);
    pixCopyInputFormat(pixd, pixs);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);

    if (cmapflag) {   /* hold out (16 - nlevels) cmap entries */
        cmap = pixcmapCreateLinear(4, nlevels);
        pixSetColormap(pixd, cmap);
    }

        /* If there is a colormap in the src, remove it */
    pixt = pixRemoveColormap(pixs, REMOVE_CMAP_TO_GRAYSCALE);
    datat = pixGetData(pixt);
    wplt = pixGetWpl(pixt);

        /* Make the appropriate table */
    if (cmapflag)
        qtab = makeGrayQuantIndexTable(nlevels);
    else
        qtab = makeGrayQuantTargetTable(16, 4);

    thresholdTo4bppLow(datad, h, wpld, datat, wplt, qtab);

    LEPT_FREE(qtab);
    pixDestroy(&pixt);
    return pixd;
}


static void
thresholdTo4bppLow(l_uint32  *datad,
                   l_int32    h,
                   l_int32    wpld,
                   l_uint32  *datas,
                   l_int32    wpls,
                   l_int32   *tab)
{
l_uint8    sval1, sval2, sval3, sval4;
l_uint16   dval;
l_int32    i, j, k;
l_uint32  *lines, *lined;

    for (i = 0; i < h; i++) {
        lines = datas + i * wpls;
        lined = datad + i * wpld;
        for (j = 0; j < wpls; j++) {
            k = 4 * j;
            sval1 = GET_DATA_BYTE(lines, k);
            sval2 = GET_DATA_BYTE(lines, k + 1);
            sval3 = GET_DATA_BYTE(lines, k + 2);
            sval4 = GET_DATA_BYTE(lines, k + 3);
            dval = (tab[sval1] << 12) | (tab[sval2] << 8) |
                   (tab[sval3] << 4) | tab[sval4];
            SET_DATA_TWO_BYTES(lined, j, dval);
        }
    }
}


/*----------------------------------------------------------------------*
 *    Simple (pixelwise) thresholding on 8 bpp with optional colormap   *
 *----------------------------------------------------------------------*/
PIX *
pixThresholdOn8bpp(PIX     *pixs,
                   l_int32  nlevels,
                   l_int32  cmapflag)
{
l_int32   *qtab;  /* quantization table */
l_int32    i, j, w, h, wpld, val, newval;
l_uint32  *datad, *lined;
PIX       *pixd;
PIXCMAP   *cmap;

    PROCNAME("pixThresholdOn8bpp");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (pixGetDepth(pixs) != 8)
        return (PIX *)ERROR_PTR("pixs not 8 bpp", procName, NULL);
    if (nlevels < 2 || nlevels > 256)
        return (PIX *)ERROR_PTR("nlevels not in [2,...,256]", procName, NULL);

        /* Get a new pixd; if there is a colormap in the src, remove it */
    if (pixGetColormap(pixs))
        pixd = pixRemoveColormap(pixs, REMOVE_CMAP_TO_GRAYSCALE);
    else
        pixd = pixCopy(NULL, pixs);

    if (cmapflag) {   /* hold out (256 - nlevels) cmap entries */
        cmap = pixcmapCreateLinear(8, nlevels);
        pixSetColormap(pixd, cmap);
    }

    if (cmapflag)
        qtab = makeGrayQuantIndexTable(nlevels);
    else
        qtab = makeGrayQuantTargetTable(nlevels, 8);

    pixGetDimensions(pixd, &w, &h, NULL);
    pixCopyResolution(pixd, pixs);
    pixCopyInputFormat(pixd, pixs);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);
    for (i = 0; i < h; i++) {
        lined = datad + i * wpld;
        for (j = 0; j < w; j++) {
            val = GET_DATA_BYTE(lined, j);
            newval = qtab[val];
            SET_DATA_BYTE(lined, j, newval);
        }
    }

    LEPT_FREE(qtab);
    return pixd;
}


/*----------------------------------------------------------------------*
 *    Arbitrary (pixelwise) thresholding from 8 bpp to 2, 4 or 8 bpp    *
 *----------------------------------------------------------------------*/
PIX *
pixThresholdGrayArb(PIX         *pixs,
                    const char  *edgevals,
                    l_int32      outdepth,
                    l_int32      use_average,
                    l_int32      setblack,
                    l_int32      setwhite)
{
l_int32   *qtab;
l_int32    w, h, d, i, j, n, wplt, wpld, val, newval;
l_uint32  *datat, *datad, *linet, *lined;
NUMA      *na;
PIX       *pixt, *pixd;
PIXCMAP   *cmap;

    PROCNAME("pixThresholdGrayArb");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 8)
        return (PIX *)ERROR_PTR("pixs not 8 bpp", procName, NULL);
    if (!edgevals)
        return (PIX *)ERROR_PTR("edgevals not defined", procName, NULL);
    if (outdepth != 0 && outdepth != 2 && outdepth != 4 && outdepth != 8)
        return (PIX *)ERROR_PTR("invalid outdepth", procName, NULL);

        /* Parse and sort (if required) the bin edge values */
    na = parseStringForNumbers(edgevals, " \t\n,");
    n = numaGetCount(na);
    if (n > 255) {
        numaDestroy(&na);
        return (PIX *)ERROR_PTR("more than 256 levels", procName, NULL);
    }
    if (outdepth == 0) {
        if (n <= 3)
            outdepth = 2;
        else if (n <= 15)
            outdepth = 4;
        else
            outdepth = 8;
    } else if (n + 1 > (1 << outdepth)) {
        L_WARNING("outdepth too small; setting to 8 bpp\n", procName);
        outdepth = 8;
    }
    numaSort(na, na, L_SORT_INCREASING);

        /* Make the quantization LUT and the colormap */
    makeGrayQuantTableArb(na, outdepth, &qtab, &cmap);
    if (use_average) {  /* use the average value in each bin */
        pixcmapDestroy(&cmap);
        makeGrayQuantColormapArb(pixs, qtab, outdepth, &cmap);
    }
    pixcmapSetBlackAndWhite(cmap, setblack, setwhite);
    numaDestroy(&na);

    if ((pixd = pixCreate(w, h, outdepth)) == NULL) {
        LEPT_FREE(qtab);
        pixcmapDestroy(&cmap);
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    }
    pixCopyResolution(pixd, pixs);
    pixCopyInputFormat(pixd, pixs);
    pixSetColormap(pixd, cmap);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);

        /* If there is a colormap in the src, remove it */
    pixt = pixRemoveColormap(pixs, REMOVE_CMAP_TO_GRAYSCALE);
    datat = pixGetData(pixt);
    wplt = pixGetWpl(pixt);

    if (outdepth == 2) {
        thresholdTo2bppLow(datad, h, wpld, datat, wplt, qtab);
    } else if (outdepth == 4) {
        thresholdTo4bppLow(datad, h, wpld, datat, wplt, qtab);
    } else {
        for (i = 0; i < h; i++) {
            lined = datad + i * wpld;
            linet = datat + i * wplt;
            for (j = 0; j < w; j++) {
                val = GET_DATA_BYTE(linet, j);
                newval = qtab[val];
                SET_DATA_BYTE(lined, j, newval);
            }
        }
    }

    LEPT_FREE(qtab);
    pixDestroy(&pixt);
    return pixd;
}


/*----------------------------------------------------------------------*
 *     Quantization tables for linear thresholds of grayscale images    *
 *----------------------------------------------------------------------*/
l_int32 *
makeGrayQuantIndexTable(l_int32  nlevels)
{
l_int32   *tab;
l_int32    i, j, thresh;

    PROCNAME("makeGrayQuantIndexTable");

    if ((tab = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32))) == NULL)
        return (l_int32 *)ERROR_PTR("calloc fail for tab", procName, NULL);
    for (i = 0; i < 256; i++) {
        for (j = 0; j < nlevels; j++) {
            thresh = 255 * (2 * j + 1) / (2 * nlevels - 2);
            if (i <= thresh) {
                tab[i] = j;
/*                fprintf(stderr, "tab[%d] = %d\n", i, j); */
                break;
            }
        }
    }
    return tab;
}


static l_int32 *
makeGrayQuantTargetTable(l_int32  nlevels,
                         l_int32  depth)
{
l_int32   *tab;
l_int32    i, j, thresh, maxval, quantval;

    PROCNAME("makeGrayQuantTargetTable");

    if ((tab = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32))) == NULL)
        return (l_int32 *)ERROR_PTR("calloc fail for tab", procName, NULL);

    maxval = (1 << depth) - 1;
    if (depth < 8)
        nlevels = 1 << depth;
    for (i = 0; i < 256; i++) {
        for (j = 0; j < nlevels; j++) {
            thresh = 255 * (2 * j + 1) / (2 * nlevels - 2);
            if (i <= thresh) {
                quantval = maxval * j / (nlevels - 1);
                tab[i] = quantval;
/*                fprintf(stderr, "tab[%d] = %d\n", i, tab[i]); */
                break;
            }
        }
    }
    return tab;
}


/*----------------------------------------------------------------------*
 *   Quantization table for arbitrary thresholding of grayscale images  *
 *----------------------------------------------------------------------*/
l_ok
makeGrayQuantTableArb(NUMA      *na,
                      l_int32    outdepth,
                      l_int32  **ptab,
                      PIXCMAP  **pcmap)
{
l_int32   i, j, n, jstart, ave, val;
l_int32  *tab;
PIXCMAP  *cmap;

    PROCNAME("makeGrayQuantTableArb");

    if (!ptab)
        return ERROR_INT("&tab not defined", procName, 1);
    *ptab = NULL;
    if (!pcmap)
        return ERROR_INT("&cmap not defined", procName, 1);
    *pcmap = NULL;
    if (!na)
        return ERROR_INT("na not defined", procName, 1);
    n = numaGetCount(na);
    if (n + 1 > (1 << outdepth))
        return ERROR_INT("more bins than cmap levels", procName, 1);

    if ((cmap = pixcmapCreate(outdepth)) == NULL)
        return ERROR_INT("cmap not made", procName, 1);
    tab = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32));
    *ptab = tab;
    *pcmap = cmap;

        /* First n bins */
    jstart = 0;
    for (i = 0; i < n; i++) {
        numaGetIValue(na, i, &val);
        ave = (jstart + val) / 2;
        pixcmapAddColor(cmap, ave, ave, ave);
        for (j = jstart; j < val; j++)
            tab[j] = i;
        jstart = val;
    }

        /* Last bin */
    ave = (jstart + 255) / 2;
    pixcmapAddColor(cmap, ave, ave, ave);
    for (j = jstart; j < 256; j++)
        tab[j] = n;

    return 0;
}


static l_int32
makeGrayQuantColormapArb(PIX       *pixs,
                         l_int32   *tab,
                         l_int32    outdepth,
                         PIXCMAP  **pcmap)
{
l_int32    i, j, index, w, h, d, nbins, wpl, factor, val;
l_int32   *bincount, *binave, *binstart;
l_uint32  *line, *data;

    PROCNAME("makeGrayQuantColormapArb");

    if (!pcmap)
        return ERROR_INT("&cmap not defined", procName, 1);
    *pcmap = NULL;
    if (!pixs)
        return ERROR_INT("pixs not defined", procName, 1);
    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 8)
        return ERROR_INT("pixs not 8 bpp", procName, 1);
    if (!tab)
        return ERROR_INT("tab not defined", procName, 1);
    nbins = tab[255] + 1;
    if (nbins > (1 << outdepth))
        return ERROR_INT("more bins than cmap levels", procName, 1);

        /* Find the count and weighted count for each bin */
    if ((bincount = (l_int32 *)LEPT_CALLOC(nbins, sizeof(l_int32))) == NULL)
        return ERROR_INT("calloc fail for bincount", procName, 1);
    if ((binave = (l_int32 *)LEPT_CALLOC(nbins, sizeof(l_int32))) == NULL) {
        LEPT_FREE(bincount);
        return ERROR_INT("calloc fail for binave", procName, 1);
    }
    factor = (l_int32)(sqrt((l_float64)(w * h) / 30000.) + 0.5);
    factor = L_MAX(1, factor);
    data = pixGetData(pixs);
    wpl = pixGetWpl(pixs);
    for (i = 0; i < h; i += factor) {
        line = data + i * wpl;
        for (j = 0; j < w; j += factor) {
            val = GET_DATA_BYTE(line, j);
            bincount[tab[val]]++;
            binave[tab[val]] += val;
        }
    }

        /* Find the smallest gray values in each bin */
    binstart = (l_int32 *)LEPT_CALLOC(nbins, sizeof(l_int32));
    for (i = 1, index = 1; i < 256; i++) {
        if (tab[i] < index) continue;
        if (tab[i] == index)
            binstart[index++] = i;
    }

        /* Get the averages.  If there are no samples in a bin, use
         * the center value of the bin. */
    *pcmap = pixcmapCreate(outdepth);
    for (i = 0; i < nbins; i++) {
        if (bincount[i]) {
            val = binave[i] / bincount[i];
        } else {  /* no samples in the bin */
            if (i < nbins - 1)
                val = (binstart[i] + binstart[i + 1]) / 2;
            else  /* last bin */
                val = (binstart[i] + 255) / 2;
        }
        pixcmapAddColor(*pcmap, val, val, val);
    }

    LEPT_FREE(bincount);
    LEPT_FREE(binave);
    LEPT_FREE(binstart);
    return 0;
}


/*--------------------------------------------------------------------*
 *                 Thresholding from 32 bpp rgb to 1 bpp              *
 *--------------------------------------------------------------------*/
PIX *
pixGenerateMaskByBand32(PIX       *pixs,
                        l_uint32   refval,
                        l_int32    delm,
                        l_int32    delp,
                        l_float32  fractm,
                        l_float32  fractp)
{
l_int32    i, j, w, h, d, wpls, wpld;
l_int32    rref, gref, bref, rval, gval, bval;
l_int32    rmin, gmin, bmin, rmax, gmax, bmax;
l_uint32   pixel;
l_uint32  *datas, *datad, *lines, *lined;
PIX       *pixd;

    PROCNAME("pixGenerateMaskByBand32");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 32)
        return (PIX *)ERROR_PTR("not 32 bpp", procName, NULL);
    if (delm < 0 || delp < 0)
        return (PIX *)ERROR_PTR("delm and delp must be >= 0", procName, NULL);
    if (fractm < 0.0 || fractm > 1.0 || fractp < 0.0 || fractp > 1.0)
        return (PIX *)ERROR_PTR("fractm and/or fractp invalid", procName, NULL);

    extractRGBValues(refval, &rref, &gref, &bref);
    if (fractm == 0.0 && fractp == 0.0) {
        rmin = rref - delm;
        gmin = gref - delm;
        bmin = bref - delm;
        rmax = rref + delm;
        gmax = gref + delm;
        bmax = bref + delm;
    } else if (delm == 0 && delp == 0) {
        rmin = (l_int32)((1.0 - fractm) * rref);
        gmin = (l_int32)((1.0 - fractm) * gref);
        bmin = (l_int32)((1.0 - fractm) * bref);
        rmax = rref + (l_int32)(fractp * (255 - rref));
        gmax = gref + (l_int32)(fractp * (255 - gref));
        bmax = bref + (l_int32)(fractp * (255 - bref));
    } else {
        L_ERROR("bad input: either (delm, delp) or (fractm, fractp) "
                "must be 0\n", procName);
        return NULL;
    }

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

    return pixd;
}


PIX *
pixGenerateMaskByDiscr32(PIX      *pixs,
                         l_uint32  refval1,
                         l_uint32  refval2,
                         l_int32   distflag)
{
l_int32    i, j, w, h, d, wpls, wpld;
l_int32    rref1, gref1, bref1, rref2, gref2, bref2, rval, gval, bval;
l_uint32   pixel, dist1, dist2;
l_uint32  *datas, *datad, *lines, *lined;
PIX       *pixd;

    PROCNAME("pixGenerateMaskByDiscr32");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 32)
        return (PIX *)ERROR_PTR("not 32 bpp", procName, NULL);
    if (distflag != L_MANHATTAN_DISTANCE && distflag != L_EUCLIDEAN_DISTANCE)
        return (PIX *)ERROR_PTR("invalid distflag", procName, NULL);

    extractRGBValues(refval1, &rref1, &gref1, &bref1);
    extractRGBValues(refval2, &rref2, &gref2, &bref2);
    pixd = pixCreate(w, h, 1);
    pixCopyResolution(pixd, pixs);
    pixCopyInputFormat(pixd, pixs);
    datas = pixGetData(pixs);
    wpls = pixGetWpl(pixs);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);
    for (i = 0; i < h; i++) {
        lines = datas + i * wpls;
        lined = datad + i * wpld;
        for (j = 0; j < w; j++) {
            pixel = lines[j];
            extractRGBValues(pixel, &rval, &gval, &bval);
            if (distflag == L_MANHATTAN_DISTANCE) {
                dist1 = L_ABS(rref1 - rval);
                dist2 = L_ABS(rref2 - rval);
                dist1 += L_ABS(gref1 - gval);
                dist2 += L_ABS(gref2 - gval);
                dist1 += L_ABS(bref1 - bval);
                dist2 += L_ABS(bref2 - bval);
            } else {
                dist1 = (rref1 - rval) * (rref1 - rval);
                dist2 = (rref2 - rval) * (rref2 - rval);
                dist1 += (gref1 - gval) * (gref1 - gval);
                dist2 += (gref2 - gval) * (gref2 - gval);
                dist1 += (bref1 - bval) * (bref1 - bval);
                dist2 += (bref2 - bval) * (bref2 - bval);
            }
            if (dist1 < dist2)
                SET_DATA_BIT(lined, j);
        }
    }

    return pixd;
}


/*----------------------------------------------------------------------*
 *                Histogram-based grayscale quantization                *
 *----------------------------------------------------------------------*/
PIX *
pixGrayQuantFromHisto(PIX       *pixd,
                      PIX       *pixs,
                      PIX       *pixm,
                      l_float32  minfract,
                      l_int32    maxsize)
{
l_int32    w, h, wd, hd, wm, hm, wpls, wplm, wpld;
l_int32    nc, nestim, i, j, vals, vald;
l_int32   *lut;
l_uint32  *datas, *datam, *datad, *lines, *linem, *lined;
NUMA      *na;
PIX       *pixmr;  /* resized mask */
PIXCMAP   *cmap;

    PROCNAME("pixGrayQuantFromHisto");

    if (!pixs || pixGetDepth(pixs) != 8)
        return (PIX *)ERROR_PTR("pixs undefined or not 8 bpp", procName, NULL);
    if (minfract < 0.01) {
        L_WARNING("minfract < 0.01; setting to 0.05\n", procName);
        minfract = 0.05;
    }
    if (maxsize < 2) {
        L_WARNING("maxsize < 2; setting to 10\n", procName);
        maxsize = 10;
    }
    if ((pixd && !pixm) || (!pixd && pixm))
        return (PIX *)ERROR_PTR("(pixd,pixm) not defined together",
                                procName, NULL);
    pixGetDimensions(pixs, &w, &h, NULL);
    if (pixd) {
        if (pixGetDepth(pixm) != 1)
            return (PIX *)ERROR_PTR("pixm not 1 bpp", procName, NULL);
        if ((cmap = pixGetColormap(pixd)) == NULL)
            return (PIX *)ERROR_PTR("pixd not cmapped", procName, NULL);
        pixGetDimensions(pixd, &wd, &hd, NULL);
        if (w != wd || h != hd)
            return (PIX *)ERROR_PTR("pixs, pixd sizes differ", procName, NULL);
        nc = pixcmapGetCount(cmap);
        nestim = nc + (l_int32)(1.5 * 255 / maxsize);
        fprintf(stderr, "nestim = %d\n", nestim);
        if (nestim > 255) {
            L_ERROR("Estimate %d colors!\n", procName, nestim);
            return (PIX *)ERROR_PTR("probably too many colors", procName, NULL);
        }
        pixGetDimensions(pixm, &wm, &hm, NULL);
        if (w != wm || h != hm) {  /* resize the mask */
            L_WARNING("mask and dest sizes not equal\n", procName);
            pixmr = pixCreateNoInit(w, h, 1);
            pixRasterop(pixmr, 0, 0, wm, hm, PIX_SRC, pixm, 0, 0);
            pixRasterop(pixmr, wm, 0, w - wm, h, PIX_SET, NULL, 0, 0);
            pixRasterop(pixmr, 0, hm, wm, h - hm, PIX_SET, NULL, 0, 0);
        } else {
            pixmr = pixClone(pixm);
        }
    } else {
        pixd = pixCreateTemplate(pixs);
        cmap = pixcmapCreate(8);
        pixSetColormap(pixd, cmap);
    }
    pixCopyResolution(pixd, pixs);
    pixCopyInputFormat(pixd, pixs);

        /* Use original mask, if it exists, to select gray pixels */
    na = pixGetGrayHistogramMasked(pixs, pixm, 0, 0, 1);

        /* Fill out the cmap with gray colors, and generate the lut
         * for pixel assignment.  Issue a warning on failure.  */
    if (numaFillCmapFromHisto(na, cmap, minfract, maxsize, &lut))
        L_ERROR("ran out of colors in cmap!\n", procName);
    numaDestroy(&na);

        /* Assign the gray pixels to their cmap indices */
    datas = pixGetData(pixs);
    datad = pixGetData(pixd);
    wpls = pixGetWpl(pixs);
    wpld = pixGetWpl(pixd);
    if (!pixm) {
        for (i = 0; i < h; i++) {
            lines = datas + i * wpls;
            lined = datad + i * wpld;
            for (j = 0; j < w; j++) {
                vals = GET_DATA_BYTE(lines, j);
                vald = lut[vals];
                SET_DATA_BYTE(lined, j, vald);
            }
        }
        LEPT_FREE(lut);
        return pixd;
    }

    datam = pixGetData(pixmr);
    wplm = pixGetWpl(pixmr);
    for (i = 0; i < h; i++) {
        lines = datas + i * wpls;
        linem = datam + i * wplm;
        lined = datad + i * wpld;
        for (j = 0; j < w; j++) {
            if (!GET_DATA_BIT(linem, j))
                continue;
            vals = GET_DATA_BYTE(lines, j);
            vald = lut[vals];
            SET_DATA_BYTE(lined, j, vald);
        }
    }
    pixDestroy(&pixmr);
    LEPT_FREE(lut);
    return pixd;
}


static l_int32
numaFillCmapFromHisto(NUMA      *na,
                      PIXCMAP   *cmap,
                      l_float32  minfract,
                      l_int32    maxsize,
                      l_int32  **plut)
{
l_int32    mincount, index, sum, wtsum, span, istart, i, val, ret;
l_int32   *iahisto, *lut;
l_float32  total;

    PROCNAME("numaFillCmapFromHisto");

    if (!plut)
        return ERROR_INT("&lut not defined", procName, 1);
    *plut = NULL;
    if (!na)
        return ERROR_INT("na not defined", procName, 1);
    if (!cmap)
        return ERROR_INT("cmap not defined", procName, 1);

    numaGetSum(na, &total);
    mincount = (l_int32)(minfract * total);
    iahisto = numaGetIArray(na);
    lut = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32));
    *plut = lut;
    index = pixcmapGetCount(cmap);  /* start with number of colors
                                     * already reserved */

        /* March through, associating colors with sets of adjacent
         * gray levels.  During the process, the LUT that gives
         * the colormap index for each gray level is computed.
         * To complete a color, either the total count must equal
         * or exceed %mincount, or the current span of colors must
         * equal or exceed %maxsize.  An empty span is not converted
         * into a color; it is simply ignored.  When a span is completed for a
         * color, the weighted color in the span is added to the colormap. */
    sum = 0;
    wtsum = 0;
    istart = 0;
    ret = 0;
    for (i = 0; i < 256; i++) {
        lut[i] = index;
        sum += iahisto[i];
        wtsum += i * iahisto[i];
        span = i - istart + 1;
        if (sum < mincount && span < maxsize)
            continue;

        if (sum == 0) {  /* empty span; don't save */
            istart = i + 1;
            continue;
        }

            /* Found new color; sum > 0 */
        val = (l_int32)((l_float32)wtsum / (l_float32)sum + 0.5);
        ret = pixcmapAddColor(cmap, val, val, val);
        istart = i + 1;
        sum = 0;
        wtsum = 0;
        index++;
    }
    if (istart < 256 && sum > 0) {  /* last one */
        span = 256 - istart;
        val = (l_int32)((l_float32)wtsum / (l_float32)sum + 0.5);
        ret = pixcmapAddColor(cmap, val, val, val);
    }

    LEPT_FREE(iahisto);
    return ret;
}


/*----------------------------------------------------------------------*
 *        Color quantize grayscale image using existing colormap        *
 *----------------------------------------------------------------------*/
PIX *
pixGrayQuantFromCmap(PIX      *pixs,
                     PIXCMAP  *cmap,
                     l_int32   mindepth)
{
l_int32    i, j, index, w, h, d, depth, wpls, wpld;
l_int32    hascolor, vals, vald;
l_int32   *tab;
l_uint32  *datas, *datad, *lines, *lined;
PIXCMAP   *cmapd;
PIX       *pixd;

    PROCNAME("pixGrayQuantFromCmap");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (pixGetColormap(pixs) != NULL) {
        L_WARNING("pixs already has a colormap; returning a copy\n", procName);
        return pixCopy(NULL, pixs);
    }
    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 8)
        return (PIX *)ERROR_PTR("pixs not 8 bpp", procName, NULL);
    if (!cmap)
        return (PIX *)ERROR_PTR("cmap not defined", procName, NULL);
    if (mindepth != 2 && mindepth != 4 && mindepth != 8)
        return (PIX *)ERROR_PTR("invalid mindepth", procName, NULL);

        /* Make sure the colormap is gray */
    pixcmapHasColor(cmap, &hascolor);
    if (hascolor) {
        L_WARNING("Converting colormap colors to gray\n", procName);
        cmapd = pixcmapColorToGray(cmap, 0.3, 0.5, 0.2);
    } else {
        cmapd = pixcmapCopy(cmap);
    }

        /* Make LUT into colormap */
    tab = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32));
    for (i = 0; i < 256; i++) {
        pixcmapGetNearestGrayIndex(cmapd, i, &index);
        tab[i] = index;
    }

    pixcmapGetMinDepth(cmap, &depth);
    depth = L_MAX(depth, mindepth);
    pixd = pixCreate(w, h, depth);
    pixSetColormap(pixd, cmapd);
    pixCopyResolution(pixd, pixs);
    pixCopyInputFormat(pixd, pixs);
    datas = pixGetData(pixs);
    datad = pixGetData(pixd);
    wpls = pixGetWpl(pixs);
    wpld = pixGetWpl(pixd);
    for (i = 0; i < h; i++) {
        lines = datas + i * wpls;
        lined = datad + i * wpld;
        for (j = 0; j < w; j++) {
            vals = GET_DATA_BYTE(lines, j);
            vald = tab[vals];
            if (depth == 2)
                SET_DATA_DIBIT(lined, j, vald);
            else if (depth == 4)
                SET_DATA_QBIT(lined, j, vald);
            else  /* depth == 8 */
                SET_DATA_BYTE(lined, j, vald);
        }
    }

    LEPT_FREE(tab);
    return pixd;
}


#if 0   /* Documentation */
/*--------------------------------------------------------------------*
 *       Implementation of binarization by dithering using LUTs       *
 *                        It is archived here.                        *
 *--------------------------------------------------------------------*/
PIX *
pixDitherToBinaryLUT(PIX     *pixs,
                     l_int32  lowerclip,
                     l_int32  upperclip)
{
l_int32    w, h, d, wplt, wpld;
l_int32   *tabval, *tab38, *tab14;
l_uint32  *datat, *datad;
l_uint32  *bufs1, *bufs2;
PIX       *pixt, *pixd;

    PROCNAME("pixDitherToBinaryLUT");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 8)
        return (PIX *)ERROR_PTR("must be 8 bpp for dithering", procName, NULL);
    if (lowerclip < 0)
        lowerclip = DEFAULT_CLIP_LOWER_1;
    if (upperclip < 0)
        upperclip = DEFAULT_CLIP_UPPER_1;

    if ((pixd = pixCreate(w, h, 1)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    pixCopyResolution(pixd, pixs);
    pixCopyInputFormat(pixd, pixs);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);

        /* Remove colormap if it exists */
    pixt = pixRemoveColormap(pixs, REMOVE_CMAP_TO_GRAYSCALE);
    datat = pixGetData(pixt);
    wplt = pixGetWpl(pixt);

        /* Two line buffers, 1 for current line and 2 for next line */
    bufs1 = (l_uint32 *)LEPT_CALLOC(wplt, sizeof(l_uint32));
    bufs2 = (l_uint32 *)LEPT_CALLOC(wplt, sizeof(l_uint32));
    if (!bufs1 || !bufs2) {
        LEPT_FREE(bufs1);
        LEPT_FREE(bufs2);
        pixDestroy(&pixd);
        pixDestroy(&pixt);
        return (PIX *)ERROR_PTR("bufs1, bufs2 not both made", procName, NULL);
    }

        /* 3 lookup tables: 1-bit value, (3/8)excess, and (1/4)excess */
    make8To1DitherTables(&tabval, &tab38, &tab14, lowerclip, upperclip);

    ditherToBinaryLUTLow(datad, w, h, wpld, datat, wplt, bufs1, bufs2,
                         tabval, tab38, tab14);

    LEPT_FREE(bufs1);
    LEPT_FREE(bufs2);
    LEPT_FREE(tabval);
    LEPT_FREE(tab38);
    LEPT_FREE(tab14);
    pixDestroy(&pixt);
    return pixd;
}

void
ditherToBinaryLUTLow(l_uint32  *datad,
                     l_int32    w,
                     l_int32    h,
                     l_int32    wpld,
                     l_uint32  *datas,
                     l_int32    wpls,
                     l_uint32  *bufs1,
                     l_uint32  *bufs2,
                     l_int32   *tabval,
                     l_int32   *tab38,
                     l_int32   *tab14)
{
l_int32      i;
l_uint32    *lined;

        /* do all lines except last line */
    memcpy(bufs2, datas, 4 * wpls);  /* prime the buffer */
    for (i = 0; i < h - 1; i++) {
        memcpy(bufs1, bufs2, 4 * wpls);
        memcpy(bufs2, datas + (i + 1) * wpls, 4 * wpls);
        lined = datad + i * wpld;
        ditherToBinaryLineLUTLow(lined, w, bufs1, bufs2,
                                 tabval, tab38, tab14, 0);
    }

        /* do last line */
    memcpy(bufs1, bufs2, 4 * wpls);
    lined = datad + (h - 1) * wpld;
    ditherToBinaryLineLUTLow(lined, w, bufs1, bufs2, tabval, tab38, tab14,  1);
    return;
}

void
ditherToBinaryLineLUTLow(l_uint32  *lined,
                         l_int32    w,
                         l_uint32  *bufs1,
                         l_uint32  *bufs2,
                         l_int32   *tabval,
                         l_int32   *tab38,
                         l_int32   *tab14,
                         l_int32    lastlineflag)
{
l_int32  j;
l_int32  oval, tab38val, tab14val;
l_uint8  rval, bval, dval;

    if (lastlineflag == 0) {
        for (j = 0; j < w - 1; j++) {
            oval = GET_DATA_BYTE(bufs1, j);
            if (tabval[oval])
                SET_DATA_BIT(lined, j);
            rval = GET_DATA_BYTE(bufs1, j + 1);
            bval = GET_DATA_BYTE(bufs2, j);
            dval = GET_DATA_BYTE(bufs2, j + 1);
            tab38val = tab38[oval];
            if (tab38val == 0)
                continue;
            tab14val = tab14[oval];
            if (tab38val < 0) {
                rval = L_MAX(0, rval + tab38val);
                bval = L_MAX(0, bval + tab38val);
                dval = L_MAX(0, dval + tab14val);
            } else {
                rval = L_MIN(255, rval + tab38val);
                bval = L_MIN(255, bval + tab38val);
                dval = L_MIN(255, dval + tab14val);
            }
            SET_DATA_BYTE(bufs1, j + 1, rval);
            SET_DATA_BYTE(bufs2, j, bval);
            SET_DATA_BYTE(bufs2, j + 1, dval);
        }

            /* do last column: j = w - 1 */
        oval = GET_DATA_BYTE(bufs1, j);
        if (tabval[oval])
            SET_DATA_BIT(lined, j);
        bval = GET_DATA_BYTE(bufs2, j);
        tab38val = tab38[oval];
        if (tab38val < 0) {
            bval = L_MAX(0, bval + tab38val);
            SET_DATA_BYTE(bufs2, j, bval);
        } else if (tab38val > 0 ) {
            bval = L_MIN(255, bval + tab38val);
            SET_DATA_BYTE(bufs2, j, bval);
        }
    } else {   /* lastlineflag == 1 */
        for (j = 0; j < w - 1; j++) {
            oval = GET_DATA_BYTE(bufs1, j);
            if (tabval[oval])
                SET_DATA_BIT(lined, j);
            rval = GET_DATA_BYTE(bufs1, j + 1);
            tab38val = tab38[oval];
            if (tab38val == 0)
                continue;
            if (tab38val < 0)
                rval = L_MAX(0, rval + tab38val);
            else
                rval = L_MIN(255, rval + tab38val);
            SET_DATA_BYTE(bufs1, j + 1, rval);
        }

            /* do last pixel: (i, j) = (h - 1, w - 1) */
        oval = GET_DATA_BYTE(bufs1, j);
        if (tabval[oval])
            SET_DATA_BIT(lined, j);
    }

    return;
}

l_ok
make8To1DitherTables(l_int32 **ptabval,
                     l_int32 **ptab38,
                     l_int32 **ptab14,
                     l_int32   lowerclip,
                     l_int32   upperclip)
{
l_int32   i;
l_int32  *tabval, *tab38, *tab14;

    PROCNAME("make8To1DitherTables");

    if (ptabval) *ptabval = NULL;
    if (ptab38) *ptab38 = NULL;
    if (ptab14) *ptab14 = NULL;
    if (!ptabval || !ptab38 || !ptab14)
        return ERROR_INT("table ptrs not all defined", procName, 1);

        /* 3 lookup tables: 1-bit value, (3/8)excess, and (1/4)excess */
    tabval = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32));
    tab38 = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32));
    tab14 = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32));
    if (!tabval || !tab38 || !tab14)
        return ERROR_INT("calloc failure to make small table", procName, 1);
    *ptabval = tabval;
    *ptab38 = tab38;
    *ptab14 = tab14;

    for (i = 0; i < 256; i++) {
        if (i <= lowerclip) {
            tabval[i] = 1;
            tab38[i] = 0;
            tab14[i] = 0;
        } else if (i < 128) {
            tabval[i] = 1;
            tab38[i] = (3 * i + 4) / 8;
            tab14[i] = (i + 2) / 4;
        } else if (i < 255 - upperclip) {
            tabval[i] = 0;
            tab38[i] = (3 * (i - 255) + 4) / 8;
            tab14[i] = ((i - 255) + 2) / 4;
        } else {  /* i >= 255 - upperclip */
            tabval[i] = 0;
            tab38[i] = 0;
            tab14[i] = 0;
        }
    }

    return 0;
}
#endif   /* Documentation */