/*!
 *  pixAddRGB()
 *
 *      Input:  pixs1, pixs2  (32 bpp RGB, or colormapped)
 *      Return: pixd, or null on error
 *
 *  Notes:
 *      (1) Clips computation to the minimum size, aligning the UL corners.
 *      (2) Removes any colormap to RGB, and ignores the LSB of each
 *          pixel word (the alpha channel).
 *      (3) Adds each component value, pixelwise, clipping to 255.
 *      (4) This is useful to combine two images where most of the
 *          pixels are essentially black, such as in pixPerceptualDiff().
 */
PIX *
pixAddRGB(PIX  *pixs1,
          PIX  *pixs2)
{
l_int32    i, j, w, h, d, w2, h2, d2, wplc1, wplc2, wpld;
l_int32    rval1, gval1, bval1, rval2, gval2, bval2, rval, gval, bval;
l_uint32  *datac1, *datac2, *datad, *linec1, *linec2, *lined;
PIX       *pixc1, *pixc2, *pixd;

    PROCNAME("pixAddRGB");

    if (!pixs1)
        return (PIX *)ERROR_PTR("pixs1 not defined", procName, NULL);
    if (!pixs2)
        return (PIX *)ERROR_PTR("pixs2 not defined", procName, NULL);
    pixGetDimensions(pixs1, &w, &h, &d);
    pixGetDimensions(pixs2, &w2, &h2, &d2);
    if (!pixGetColormap(pixs1) && d != 32)
        return (PIX *)ERROR_PTR("pixs1 not cmapped or rgb", procName, NULL);
    if (!pixGetColormap(pixs2) && d2 != 32)
        return (PIX *)ERROR_PTR("pixs2 not cmapped or rgb", procName, NULL);
    if (pixGetColormap(pixs1))
        pixc1 = pixRemoveColormap(pixs1, REMOVE_CMAP_TO_FULL_COLOR);
    else
        pixc1 = pixClone(pixs1);
    if (pixGetColormap(pixs2))
        pixc2 = pixRemoveColormap(pixs2, REMOVE_CMAP_TO_FULL_COLOR);
    else
        pixc2 = pixClone(pixs2);

    w = L_MIN(w, w2);
    h = L_MIN(h, h2);
    pixd = pixCreate(w, h, 32);
    pixCopyResolution(pixd, pixs1);
    datac1 = pixGetData(pixc1);
    datac2 = pixGetData(pixc2);
    datad = pixGetData(pixd);
    wplc1 = pixGetWpl(pixc1);
    wplc2 = pixGetWpl(pixc2);
    wpld = pixGetWpl(pixd);
    for (i = 0; i < h; i++) {
        linec1 = datac1 + i * wplc1;
        linec2 = datac2 + i * wplc2;
        lined = datad + i * wpld;
        for (j = 0; j < w; j++) {
            extractRGBValues(linec1[j], &rval1, &gval1, &bval1);
            extractRGBValues(linec2[j], &rval2, &gval2, &bval2);
            rval = L_MIN(255, rval1 + rval2);
            gval = L_MIN(255, gval1 + gval2);
            bval = L_MIN(255, bval1 + bval2);
            composeRGBPixel(rval, gval, bval, lined + j);
        }
    }

    pixDestroy(&pixc1);
    pixDestroy(&pixc2);
    return pixd;
}

/*!
 *  pixRotateShearIP()
 *
 *      Input:  pixs (any depth; not colormapped)
 *              xcen, ycen (center of rotation)
 *              angle (radians)
 *              incolor (L_BRING_IN_WHITE, L_BRING_IN_BLACK)
 *      Return: 0 if OK; 1 on error
 *
 *  Notes:
 *      (1) This does an in-place rotation of the image
 *          about the image center, using the 3-shear method.
 *      (2) A positive angle gives a clockwise rotation.
 *      (3) 3-shear rotation by a specified angle is equivalent
 *          to the sequential transformations
 *            y' = y + tan(angle/2) * (x - xcen)      for first y-shear
 *            x' = x + sin(angle) * (y - ycen)        for x-shear
 *            y' = y + tan(angle/2) * (x - xcen)      for second y-shear
 *      (4) Computation of tan(angle) is performed in the shear operations.
 *      (5) This brings in 'incolor' pixels from outside the image.
 *      (6) The pix cannot be colormapped, because the in-place operation
 *          only blits in 0 or 1 bits, not an arbitrary colormap index.
 */
l_int32
pixRotateShearIP(PIX       *pixs,
                 l_int32    xcen,
                 l_int32    ycen,
                 l_float32  angle,
                 l_int32    incolor)
{
l_float32  hangle;

    PROCNAME("pixRotateShearIP");

    if (!pixs)
        return ERROR_INT("pixs not defined", procName, 1);
    if (incolor != L_BRING_IN_WHITE && incolor != L_BRING_IN_BLACK)
        return ERROR_INT("invalid value for incolor", procName, 1);
    if (pixGetColormap(pixs) != NULL)
        return ERROR_INT("pixs is colormapped", procName, 1);

    if (angle == 0.0)
        return 0;

    hangle = atan(sin(angle));
    pixHShearIP(pixs, ycen, angle / 2., incolor);
    pixVShearIP(pixs, xcen, hangle, incolor);
    pixHShearIP(pixs, ycen, angle / 2., incolor);

    return 0;
}

static l_int32
test_8bpp_trans(L_REGPARAMS  *rp)
{
l_int32   same, transp;
FILE     *fp;
PIX      *pix1, *pix2, *pix3;
PIXCMAP  *cmap;

    pix1 = pixRead("wyom.jpg");
    pix2 = pixColorSegment(pix1, 75, 10, 8, 7);
    cmap = pixGetColormap(pix2);
    pixcmapSetAlpha(cmap, 0, 0);  /* set blueish sky color to transparent */
    pixWrite("/tmp/regout/8bpp-trans.png", pix2, IFF_PNG);
    pix3 = pixRead("/tmp/regout/8bpp-trans.png");
    pixEqual(pix2, pix3, &same);
    if (same)
        fprintf(stderr, "8bpp_trans: success\n");
    else
        fprintf(stderr, "8bpp_trans: bad output\n");
    pixDisplayWithTitle(pix3, 700, 0, NULL, rp->display);
    pixDestroy(&pix1);
    pixDestroy(&pix2);
    pixDestroy(&pix3);
    fp = fopenReadStream("/tmp/regout/8bpp-trans.png");
    fgetPngColormapInfo(fp, &cmap, &transp);
    fclose(fp);
    if (transp)
        fprintf(stderr, "8bpp_trans: correct -- transparency found\n");
    else
        fprintf(stderr, "8bpp_trans: error -- no transparency found!\n");
    if (rp->display) pixcmapWriteStream(stderr, cmap);
    pixcmapDestroy(&cmap);
    return same;
}

/*!
 * \brief   pixRotateShearIP()
 *
 * \param[in]    pixs any depth; not colormapped
 * \param[in]    xcen, ycen center of rotation
 * \param[in]    angle radians
 * \param[in]    incolor L_BRING_IN_WHITE, L_BRING_IN_BLACK
 * \return  0 if OK; 1 on error
 *
 * <pre>
 * Notes:
 *      (1) This does an in-place rotation of the image about the
 *          specified point, using the 3-shear method.  It should only
 *          be used for angles smaller than LIMIT_SHEAR_ANGLE.
 *          For larger angles, a warning is issued.
 *      (2) A positive angle gives a clockwise rotation.
 *      (3) 3-shear rotation by a specified angle is equivalent
 *          to the sequential transformations
 *            y' = y + tan(angle/2) * (x - xcen)      for first y-shear
 *            x' = x + sin(angle) * (y - ycen)        for x-shear
 *            y' = y + tan(angle/2) * (x - xcen)      for second y-shear
 *      (4) Computation of tan(angle) is performed in the shear operations.
 *      (5) This brings in 'incolor' pixels from outside the image.
 *      (6) The pix cannot be colormapped, because the in-place operation
 *          only blits in 0 or 1 bits, not an arbitrary colormap index.
 * </pre>
 */
l_int32
pixRotateShearIP(PIX       *pixs,
                 l_int32    xcen,
                 l_int32    ycen,
                 l_float32  angle,
                 l_int32    incolor)
{
l_float32  hangle;

    PROCNAME("pixRotateShearIP");

    if (!pixs)
        return ERROR_INT("pixs not defined", procName, 1);
    if (incolor != L_BRING_IN_WHITE && incolor != L_BRING_IN_BLACK)
        return ERROR_INT("invalid value for incolor", procName, 1);
    if (pixGetColormap(pixs) != NULL)
        return ERROR_INT("pixs is colormapped", procName, 1);

    if (angle == 0.0)
        return 0;
    if (L_ABS(angle) > LIMIT_SHEAR_ANGLE) {
        L_WARNING("%6.2f radians; large angle for in-place 3-shear rotation\n",
                  procName, L_ABS(angle));
    }

    hangle = atan(sin(angle));
    pixHShearIP(pixs, ycen, angle / 2., incolor);
    pixVShearIP(pixs, xcen, hangle, incolor);
    pixHShearIP(pixs, ycen, angle / 2., incolor);
    return 0;
}

/*!
 *  pixRankFilter()
 *
 *      Input:  pixs (8 or 32 bpp; no colormap)
 *              wf, hf  (width and height of filter; each is >= 1)
 *              rank (in [0.0 ... 1.0])
 *      Return: pixd (of rank values), or null on error
 *
 *  Notes:
 *      (1) This defines, for each pixel in pixs, a neighborhood of
 *          pixels given by a rectangle "centered" on the pixel.
 *          This set of wf*hf pixels has a distribution of values.
 *          For each component, if the values are sorted in increasing
 *          order, we choose the component such that rank*(wf*hf-1)
 *          pixels have a lower or equal value and
 *          (1-rank)*(wf*hf-1) pixels have an equal or greater value.
 *      (2) See notes in pixRankFilterGray() for further details.
 */
PIX  *
pixRankFilter(PIX       *pixs,
              l_int32    wf,
              l_int32    hf,
              l_float32  rank)
{
l_int32  d;

    PROCNAME("pixRankFilter");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (pixGetColormap(pixs) != NULL)
        return (PIX *)ERROR_PTR("pixs has colormap", procName, NULL);
    d = pixGetDepth(pixs);
    if (d != 8 && d != 32)
        return (PIX *)ERROR_PTR("pixs not 8 or 32 bpp", procName, NULL);
    if (wf < 1 || hf < 1)
        return (PIX *)ERROR_PTR("wf < 1 || hf < 1", procName, NULL);
    if (rank < 0.0 || rank > 1.0)
        return (PIX *)ERROR_PTR("rank must be in [0.0, 1.0]", procName, NULL);
    if (wf == 1 && hf == 1)   /* no-op */
        return pixCopy(NULL, pixs);

    if (d == 8)
        return pixRankFilterGray(pixs, wf, hf, rank);
    else  /* d == 32 */
        return pixRankFilterRGB(pixs, wf, hf, rank);
}

/*!
 *  pixaAnyColormaps()
 *
 *      Input:  pixa
 *              &hascmap (<return> 1 if any pix has a colormap; 0 otherwise)
 *      Return: 0 if OK; 1 on error
 */
l_int32
pixaAnyColormaps(PIXA     *pixa,
                 l_int32  *phascmap)
{
l_int32   i, n;
PIX      *pix;
PIXCMAP  *cmap;

    PROCNAME("pixaAnyColormaps");

    if (!phascmap)
        return ERROR_INT("&hascmap not defined", procName, 1);
    *phascmap = 0;
    if (!pixa)
        return ERROR_INT("pixa not defined", procName, 1);

    n = pixaGetCount(pixa);
    for (i = 0; i < n; i++) {
        pix = pixaGetPix(pixa, i, L_CLONE);
        cmap = pixGetColormap(pix);
        pixDestroy(&pix);
        if (cmap) {
            *phascmap = 1;
            return 0;
        }
    }

    return 0;
}

/*!
 * \brief   pixGetAutoFormat()
 *
 * \param[in]    pix
 * \param[in]    &format
 * \return  0 if OK, 1 on error
 *
 * <pre>
 * Notes:
 *      (1) The output formats are restricted to tiff, jpeg and png
 *          because these are the most commonly used image formats and
 *          the ones that are typically installed with leptonica.
 *      (2) This decides what compression to use based on the pix.
 *          It chooses tiff-g4 if 1 bpp without a colormap, jpeg with
 *          quality 75 if grayscale, rgb or rgba (where it loses
 *          the alpha layer), and lossless png for all other situations.
 * </pre>
 */
l_int32
pixGetAutoFormat(PIX      *pix,
                 l_int32  *pformat)
{
l_int32   d;
PIXCMAP  *cmap;

    PROCNAME("pixGetAutoFormat");

    if (!pformat)
        return ERROR_INT("&format not defined", procName, 0);
    *pformat = IFF_UNKNOWN;
    if (!pix)
        return ERROR_INT("pix not defined", procName, 0);

    d = pixGetDepth(pix);
    cmap = pixGetColormap(pix);
    if (d == 1 && !cmap) {
        *pformat = IFF_TIFF_G4;
    } else if ((d == 8 && !cmap) || d == 24 || d == 32) {
        *pformat = IFF_JFIF_JPEG;
    } else {
        *pformat = IFF_PNG;
    }

    return 0;
}

/*!
 *  pixPrintStreamInfo()
 *
 *      Input:  fp (file stream)
 *              pix
 *              text (<optional> identifying string; can be null)
 *      Return: 0 if OK, 1 on error
 */
l_int32
pixPrintStreamInfo(FILE        *fp,
                   PIX         *pix,
                   const char  *text)
{
PIXCMAP  *cmap;

    PROCNAME("pixPrintStreamInfo");

    if (!fp)
        return ERROR_INT("fp not defined", procName, 1);
    if (!pix)
        return ERROR_INT("pix not defined", procName, 1);

    if (text)
        fprintf(fp, "  Pix Info for %s:\n", text);
    fprintf(fp, "    width = %d, height = %d, depth = %d\n",
               pixGetWidth(pix), pixGetHeight(pix), pixGetDepth(pix));
    fprintf(fp, "    wpl = %d, data = %p, refcount = %d\n",
               pixGetWpl(pix), pixGetData(pix), pixGetRefcount(pix));
    if ((cmap = pixGetColormap(pix)) != NULL)
        pixcmapWriteStream(fp, cmap);
    else
        fprintf(fp, "    no colormap\n");

    return 0;
}

/*!
 *  pixBilateralGray()
 *
 *      Input:  pixs (8 bpp gray)
 *              spatial_stdev  (of gaussian kernel; in pixels, > 0.5)
 *              range_stdev  (of gaussian range kernel; > 5.0; typ. 50.0)
 *              ncomps (number of intermediate sums J(k,x); in [4 ... 30])
 *              reduction  (1, 2 or 4)
 *      Return: pixd (8 bpp bilateral filtered image), or null on error
 *
 *  Notes:
 *      (1) See pixBilateral() for constraints on the input parameters.
 *      (2) See pixBilateral() for algorithm details.
 */
PIX *
pixBilateralGray(PIX *pixs,
                 l_float32 spatial_stdev,
                 l_float32 range_stdev,
                 l_int32 ncomps,
                 l_int32 reduction) {
    l_float32 sstdev;  /* scaled spatial stdev */
    PIX *pixd;
    L_BILATERAL *bil;

    PROCNAME("pixBilateralGray");

    if (!pixs || pixGetColormap(pixs))
        return (PIX *) ERROR_PTR("pixs not defined or cmapped", procName, NULL);
    if (pixGetDepth(pixs) != 8)
        return (PIX *) ERROR_PTR("pixs not 8 bpp gray", procName, NULL);
    if (reduction != 1 && reduction != 2 && reduction != 4)
        return (PIX *) ERROR_PTR("reduction invalid", procName, NULL);
    sstdev = spatial_stdev / (l_float32) reduction;  /* reduced spat. stdev */
    if (sstdev < 0.5)
        return (PIX *) ERROR_PTR("sstdev < 0.5", procName, NULL);
    if (range_stdev <= 5.0)
        return (PIX *) ERROR_PTR("range_stdev <= 5.0", procName, NULL);
    if (ncomps < 4 || ncomps > 30)
        return (PIX *) ERROR_PTR("ncomps not in [4 ... 30]", procName, NULL);
    if (ncomps * range_stdev < 100.0)
        return (PIX *) ERROR_PTR("ncomps * range_stdev < 100.0", procName, NULL);

    bil = bilateralCreate(pixs, spatial_stdev, range_stdev, ncomps, reduction);
    if (!bil) return (PIX *) ERROR_PTR("bil not made", procName, NULL);
    pixd = bilateralApply(bil);
    bilateralDestroy(&bil);
    return pixd;
}

/*!
 *  pixBlockBilateralExact()
 *
 *      Input:  pixs (8 bpp gray or 32 bpp rgb)
 *              spatial_stdev (> 0.0)
 *              range_stdev (> 0.0)
 *      Return: pixd (8 bpp or 32 bpp bilateral filtered image)
 *
 *  Notes:
 *      (1) See pixBilateralExact().  This provides an interface using
 *          the standard deviations of the spatial and range filters.
 *      (2) The convolution window halfwidth is 2 * spatial_stdev,
 *          and the square filter size is 4 * spatial_stdev + 1.
 *          The kernel captures 95% of total energy.  This is compensated
 *          by normalization.
 *      (3) The range_stdev is analogous to spatial_halfwidth in the
 *          grayscale domain [0...255], and determines how much damping of the
 *          smoothing operation is applied across edges.  The larger this
 *          value is, the smaller the damping.  The smaller the value, the
 *          more edge details are preserved.  These approximations are useful
 *          for deciding the appropriate cutoff.
 *              kernel[1 * stdev] ~= 0.6  * kernel[0]
 *              kernel[2 * stdev] ~= 0.14 * kernel[0]
 *              kernel[3 * stdev] ~= 0.01 * kernel[0]
 *          If range_stdev is infinite there is no damping, and this
 *          becomes a conventional gaussian smoothing.
 *          This value does not affect the run time.
 *      (4) If range_stdev is negative or zero, the range kernel is
 *          ignored and this degenerates to a straight gaussian convolution.
 *      (5) This is very slow for large spatial filters.  The time
 *          on a 3GHz pentium is roughly
 *             T = 1.2 * 10^-8 * (A * sh^2)  sec
 *          where A = # of pixels, sh = spatial halfwidth of filter.
 */
PIX *
pixBlockBilateralExact(PIX *pixs,
                       l_float32 spatial_stdev,
                       l_float32 range_stdev) {
    l_int32 d, halfwidth;
    L_KERNEL *spatial_kel, *range_kel;
    PIX *pixd;

    PROCNAME("pixBlockBilateralExact");

    if (!pixs)
        return (PIX *) ERROR_PTR("pixs not defined", procName, NULL);
    d = pixGetDepth(pixs);
    if (d != 8 && d != 32)
        return (PIX *) ERROR_PTR("pixs not 8 or 32 bpp", procName, NULL);
    if (pixGetColormap(pixs) != NULL)
        return (PIX *) ERROR_PTR("pixs is cmapped", procName, NULL);
    if (spatial_stdev <= 0.0)
        return (PIX *) ERROR_PTR("invalid spatial stdev", procName, NULL);
    if (range_stdev <= 0.0)
        return (PIX *) ERROR_PTR("invalid range stdev", procName, NULL);

    halfwidth = 2 * spatial_stdev;
    spatial_kel = makeGaussianKernel(halfwidth, halfwidth, spatial_stdev, 1.0);
    range_kel = makeRangeKernel(range_stdev);
    pixd = pixBilateralExact(pixs, spatial_kel, range_kel);
    kernelDestroy(&spatial_kel);
    kernelDestroy(&range_kel);
    return pixd;
}

/*!
 *  pixSerializeToMemory()
 *
 *      Input:  pixs (all depths, colormap OK)
 *              &data (<return> serialized data in memory)
 *              &nbytes (<return> number of bytes in data string)
 *      Return: 0 if OK, 1 on error
 *
 *  Notes:
 *      (1) This does a fast serialization of the principal elements
 *          of the pix, as follows:
 *            "spix"    (4 bytes) -- ID for file type
 *            w         (4 bytes)
 *            h         (4 bytes)
 *            d         (4 bytes)
 *            wpl       (4 bytes)
 *            ncolors   (4 bytes) -- in colormap; 0 if there is no colormap
 *            cdata     (4 * ncolors)  -- size of serialized colormap array
 *            rdatasize (4 bytes) -- size of serialized raster data
 *                                   = 4 * wpl * h
 *            rdata     (rdatasize)
 */
l_int32
pixSerializeToMemory(PIX        *pixs,
                     l_uint32  **pdata,
                     size_t     *pnbytes)
{
char      *id;
l_int32    w, h, d, wpl, rdatasize, ncolors, nbytes, index;
l_uint8   *cdata;  /* data in colormap array (4 bytes/color table entry) */
l_uint32  *data;
l_uint32  *rdata;  /* data in pix raster */
PIXCMAP   *cmap;

    PROCNAME("pixSerializeToMemory");

    if (!pdata || !pnbytes)
        return ERROR_INT("&data and &nbytes not both defined", procName, 1);
    *pdata = NULL;
    *pnbytes = 0;
    if (!pixs)
        return ERROR_INT("pixs not defined", procName, 1);

    pixGetDimensions(pixs, &w, &h, &d);
    wpl = pixGetWpl(pixs);
    rdata = pixGetData(pixs);
    rdatasize = 4 * wpl * h;
    ncolors = 0;
    cdata = NULL;
    if ((cmap = pixGetColormap(pixs)) != NULL)
        pixcmapSerializeToMemory(cmap, 4, &ncolors, &cdata);

    nbytes = 24 + 4 * ncolors + 4 + rdatasize;
    if ((data = (l_uint32 *)CALLOC(nbytes / 4, sizeof(l_uint32))) == NULL)
        return ERROR_INT("data not made", procName, 1);
    *pdata = data;
    *pnbytes = nbytes;
    id = (char *)data;
    id[0] = 's';
    id[1] = 'p';
    id[2] = 'i';
    id[3] = 'x';
    data[1] = w;
    data[2] = h;
    data[3] = d;
    data[4] = wpl;
    data[5] = ncolors;
    if (ncolors > 0)
        memcpy((char *)(data + 6), (char *)cdata, 4 * ncolors);
    index = 6 + ncolors;
    data[index] = rdatasize;
    memcpy((char *)(data + index + 1), (char *)rdata, rdatasize);

#if  DEBUG_SERIALIZE
    fprintf(stderr, "Serialize:   "
            "raster size = %d, ncolors in cmap = %d, total bytes = %d\n",
            rdatasize, ncolors, nbytes);
#endif  /* DEBUG_SERIALIZE */

    FREE(cdata);
    return 0;
}

/*!
 *  pixDisplayWriteFormat()
 *
 *      Input:  pix (1, 2, 4, 8, 16, 32 bpp)
 *              reduction (-1 to reset/erase; 0 to disable;
 *                         otherwise this is a reduction factor)
 *              format (IFF_PNG or IFF_JFIF_JPEG)
 *      Return: 0 if OK; 1 on error
 *
 *  Notes:
 *      (1) This writes files if reduction > 0.  These can be displayed using
 *            pixDisplayMultiple("/tmp/junk_write_display*");
 *      (2) All previously written files can be erased by calling with
 *          reduction < 0; the value of pixs is ignored.
 *      (3) If reduction > 1 and depth == 1, this does a scale-to-gray
 *          reduction.
 *      (4) This function uses a static internal variable to number
 *          output files written by a single process.  Behavior
 *          with a shared library may be unpredictable.
 *      (5) Output file format is as follows:
 *            format == IFF_JFIF_JPEG:
 *                png if d < 8 or d == 16 or if the output pix
 *                has a colormap.   Otherwise, output is jpg.
 *            format == IFF_PNG:
 *                png (lossless) on all images.
 *      (6) For 16 bpp, the choice of full dynamic range with log scale
 *          is the best for displaying these images.  Alternative outputs are
 *             pix8 = pixMaxDynamicRange(pixt, L_LINEAR_SCALE);
 *             pix8 = pixConvert16To8(pixt, 0);  // low order byte
 *             pix8 = pixConvert16To8(pixt, 1);  // high order byte
 */
l_int32
pixDisplayWriteFormat(PIX     *pixs,
                      l_int32  reduction,
                      l_int32  format)
{
char            buffer[L_BUF_SIZE];
l_float32       scale;
PIX            *pixt, *pix8;
static l_int32  index = 0;  /* caution: not .so or thread safe */

    PROCNAME("pixDisplayWriteFormat");

    if (reduction == 0) return 0;

    if (reduction < 0) {
        index = 0;  /* reset; this will cause erasure at next call to write */
        return 0;
    }

    if (format != IFF_JFIF_JPEG && format != IFF_PNG)
        return ERROR_INT("invalid format", procName, 1);
    if (!pixs)
        return ERROR_INT("pixs not defined", procName, 1);

    if (index == 0) {
        snprintf(buffer, L_BUF_SIZE,
           "rm -f /tmp/junk_write_display.*.png /tmp/junk_write_display.*.jpg");
        system(buffer);
    }
    index++;

    if (reduction == 1)
        pixt = pixClone(pixs);
    else {
        scale = 1. / (l_float32)reduction;
        if (pixGetDepth(pixs) == 1)
            pixt = pixScaleToGray(pixs, scale);
        else
            pixt = pixScale(pixs, scale, scale);
    }

    if (pixGetDepth(pixt) == 16) {
        pix8 = pixMaxDynamicRange(pixt, L_LOG_SCALE);
        snprintf(buffer, L_BUF_SIZE, "/tmp/junk_write_display.%03d.png", index);
        pixWrite(buffer, pix8, IFF_PNG);
        pixDestroy(&pix8);
    }
    else if (pixGetDepth(pixt) < 8 || pixGetColormap(pixt) ||
             format == IFF_PNG) {
        snprintf(buffer, L_BUF_SIZE, "/tmp/junk_write_display.%03d.png", index);
        pixWrite(buffer, pixt, IFF_PNG);
    }
    else {
        snprintf(buffer, L_BUF_SIZE, "/tmp/junk_write_display.%03d.jpg", index);
        pixWrite(buffer, pixt, format);
    }
    pixDestroy(&pixt);

    return 0;
}

/*!
 *  pixCopyColormap()
 *
 *      Input:  src and dest Pix
 *      Return: 0 if OK, 1 on error
 *
 *  Notes:
 *      (1) This always destroys any colormap in pixd (except if
 *          the operation is a no-op.
 */
l_int32
pixCopyColormap(PIX  *pixd,
                PIX  *pixs)
{
PIXCMAP  *cmaps, *cmapd;

    PROCNAME("pixCopyColormap");

    if (!pixs)
        return ERROR_INT("pixs not defined", procName, 1);
    if (!pixd)
        return ERROR_INT("pixd not defined", procName, 1);
    if (pixs == pixd)
        return 0;   /* no-op */

    pixDestroyColormap(pixd);
    if ((cmaps = pixGetColormap(pixs)) == NULL)  /* not an error */
        return 0;

    if ((cmapd = pixcmapCopy(cmaps)) == NULL)
        return ERROR_INT("cmapd not made", procName, 1);
    pixSetColormap(pixd, cmapd);

    return 0;
}

/*!
 *  pixRasteropHip()
 *
 *      Input:  pixd (in-place operation)
 *              by  (top of horizontal band)
 *              bh  (height of horizontal band)
 *              hshift (horizontal shift of band; hshift > 0 is to right)
 *              incolor (L_BRING_IN_WHITE, L_BRING_IN_BLACK)
 *      Return: 0 if OK; 1 on error
 *
 *  Notes:
 *      (1) This rasterop translates a horizontal band of the
 *          image either left or right, bringing in either white
 *          or black pixels from outside the image.
 *      (2) The horizontal band extends the full width of pixd.
 *      (3) If a colormap exists, the nearest color to white or black
 *          is brought in.
 */
l_int32
pixRasteropHip(PIX     *pixd,
               l_int32  by,
               l_int32  bh,
               l_int32  hshift,
               l_int32  incolor)
{
l_int32   w, h, d, index, op;
PIX      *pixt;
PIXCMAP  *cmap;

    PROCNAME("pixRasteropHip");

    if (!pixd)
        return ERROR_INT("pixd not defined", procName, 1);
    if (incolor != L_BRING_IN_WHITE && incolor != L_BRING_IN_BLACK)
        return ERROR_INT("invalid value for incolor", procName, 1);
    if (bh <= 0)
        return ERROR_INT("bh must be > 0", procName, 1);

    if (hshift == 0)
        return 0;

    pixGetDimensions(pixd, &w, &h, &d);
    rasteropHipLow(pixGetData(pixd), h, d, pixGetWpl(pixd), by, bh, hshift);

    cmap = pixGetColormap(pixd);
    if (!cmap) {
        if ((d == 1 && incolor == L_BRING_IN_BLACK) ||
            (d > 1 && incolor == L_BRING_IN_WHITE))
            op = PIX_SET;
        else
            op = PIX_CLR;

            /* Set the pixels brought in at left or right */
        if (hshift > 0)
            pixRasterop(pixd, 0, by, hshift, bh, op, NULL, 0, 0);
        else  /* hshift < 0 */
            pixRasterop(pixd, w + hshift, by, -hshift, bh, op, NULL, 0, 0);
        return 0;
    }

        /* Get the nearest index and fill with that */
    if (incolor == L_BRING_IN_BLACK)
        pixcmapGetRankIntensity(cmap, 0.0, &index);
    else  /* white */
        pixcmapGetRankIntensity(cmap, 1.0, &index);
    pixt = pixCreate(L_ABS(hshift), bh, d);
    pixSetAllArbitrary(pixt, index);
    if (hshift > 0)
        pixRasterop(pixd, 0, by, hshift, bh, PIX_SRC, pixt, 0, 0);
    else  /* hshift < 0 */
        pixRasterop(pixd, w + hshift, by, -hshift, bh, PIX_SRC, pixt, 0, 0);
    pixDestroy(&pixt);
    return 0;
}

l_int32 main(int    argc,
             char **argv)
{
l_uint32     *colors;
l_int32       ncolors;
PIX          *pix1, *pix2, *pix3;
L_REGPARAMS  *rp;

    if (regTestSetup(argc, argv, &rp))
        return 1;

        /* Find the most populated colors */
    pix1 = pixRead("fish24.jpg");
    pixGetMostPopulatedColors(pix1, 2, 3, 10, &colors, NULL);
    pix2 = pixDisplayColorArray(colors, 10, 190, 5, 1);
    pixDisplayWithTitle(pix2, 0, 0, NULL, rp->display);
    regTestWritePixAndCheck(rp, pix2, IFF_PNG);  /* 0 */
    lept_free(colors);
    pixDestroy(&pix2);

        /* Do a simple color quantization with sigbits = 2 */
    pix2 = pixSimpleColorQuantize(pix1, 2, 3, 10);
    pixDisplayWithTitle(pix2, 0, 400, NULL, rp->display);
    regTestWritePixAndCheck(rp, pix2, IFF_PNG);  /* 1 */
    pix3 = pixRemoveColormap(pix2, REMOVE_CMAP_TO_FULL_COLOR);
    regTestComparePix(rp, pix2, pix3);  /* 2 */
    pixNumColors(pix3, 1, &ncolors);
    regTestCompareValues(rp, ncolors, 10, 0.0);  /* 3 */
    pixDestroy(&pix1);
    pixDestroy(&pix2);
    pixDestroy(&pix3);

        /* Do a simple color quantization with sigbits = 3 */
    pix1 = pixRead("wyom.jpg");
    pixNumColors(pix1, 1, &ncolors);  /* >255, so should give 0 */
    regTestCompareValues(rp, ncolors, 0, 0.0);  /* 4 */
    pix2 = pixSimpleColorQuantize(pix1, 3, 3, 20);
    pixDisplayWithTitle(pix2, 1000, 0, NULL, rp->display);
    regTestWritePixAndCheck(rp, pix2, IFF_PNG);  /* 5 */
    ncolors = pixcmapGetCount(pixGetColormap(pix2));
    regTestCompareValues(rp, ncolors, 20, 0.0);  /* 6 */
    pixDestroy(&pix1);
    pixDestroy(&pix2);

        /* Find the number of perceptually significant gray intensities */
    pix1 = pixRead("marge.jpg");
    pix2 = pixConvertTo8(pix1, 0);
    pixNumSignificantGrayColors(pix2, 20, 236, 0.0001, 1, &ncolors);
    regTestCompareValues(rp, ncolors, 219, 0.0);  /* 7 */
    pixDestroy(&pix1);
    pixDestroy(&pix2);

    return regTestCleanup(rp);
}

/*!
 *  pixOtsuThreshOnBackgroundNorm()
 *
 *      Input:  pixs (8 bpp grayscale; not colormapped)
 *              pixim (<optional> 1 bpp 'image' mask; can be null)
 *              sx, sy (tile size in pixels)
 *              thresh (threshold for determining foreground)
 *              mincount (min threshold on counts in a tile)
 *              bgval (target bg val; typ. > 128)
 *              smoothx (half-width of block convolution kernel width)
 *              smoothy (half-width of block convolution kernel height)
 *              scorefract (fraction of the max Otsu score; typ. 0.1)
 *              &thresh (<optional return> threshold value that was
 *                       used on the normalized image)
 *      Return: pixd (1 bpp thresholded image), or null on error
 *
 *  Notes:
 *      (1) This does background normalization followed by Otsu
 *          thresholding.  Otsu binarization attempts to split the
 *          image into two roughly equal sets of pixels, and it does
 *          a very poor job when there are large amounts of dark
 *          background.  By doing a background normalization first,
 *          to get the background near 255, we remove this problem.
 *          Then we use a modified Otsu to estimate the best global
 *          threshold on the normalized image.
 *      (2) See pixBackgroundNorm() for meaning and typical values
 *          of input parameters.  For a start, you can try:
 *            sx, sy = 10, 15
 *            thresh = 100
 *            mincount = 50
 *            bgval = 255
 *            smoothx, smoothy = 2
 */
PIX *
pixOtsuThreshOnBackgroundNorm(PIX       *pixs,
                              PIX       *pixim,
                              l_int32    sx,
                              l_int32    sy,
                              l_int32    thresh,
                              l_int32    mincount,
                              l_int32    bgval,
                              l_int32    smoothx,
                              l_int32    smoothy,
                              l_float32  scorefract,
                              l_int32   *pthresh)
{
l_int32   w, h;
l_uint32  val;
PIX      *pixn, *pixt, *pixd;

    PROCNAME("pixOtsuThreshOnBackgroundNorm");

    if (pthresh) *pthresh = 0;
    if (!pixs || pixGetDepth(pixs) != 8)
        return (PIX *)ERROR_PTR("pixs undefined or not 8 bpp", procName, NULL);
    if (pixGetColormap(pixs))
        return (PIX *)ERROR_PTR("pixs is colormapped", procName, NULL);
    if (sx < 4 || sy < 4)
        return (PIX *)ERROR_PTR("sx and sy must be >= 4", procName, NULL);
    if (mincount > sx * sy) {
        L_WARNING("mincount too large for tile size\n", procName);
        mincount = (sx * sy) / 3;
    }

    pixn = pixBackgroundNorm(pixs, pixim, NULL, sx, sy, thresh,
                             mincount, bgval, smoothx, smoothy);
    if (!pixn)
        return (PIX *)ERROR_PTR("pixn not made", procName, NULL);

        /* Just use 1 tile for a global threshold, which is stored
         * as a single pixel in pixt. */
    pixGetDimensions(pixn, &w, &h, NULL);
    pixOtsuAdaptiveThreshold(pixn, w, h, 0, 0, scorefract, &pixt, &pixd);
    pixDestroy(&pixn);

    if (pixt && pthresh) {
        pixGetPixel(pixt, 0, 0, &val);
        *pthresh = val;
    }
    pixDestroy(&pixt);

    if (!pixd)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    else
        return pixd;
}

/*!
 *  pixCloseGray3()
 *
 *      Input:  pixs (8 bpp, not cmapped)
 *              hsize  (1 or 3)
 *              vsize  (1 or 3)
 *      Return: pixd, or null on error
 *
 *  Notes:
 *      (1) Special case for 1x3, 3x1 or 3x3 brick sel (all hits)
 *      (2) If hsize = vsize = 1, just returns a copy.
 */
PIX *
pixCloseGray3(PIX     *pixs,
              l_int32  hsize,
              l_int32  vsize)
{
PIX  *pixt, *pixb, *pixbd, *pixd;

    PROCNAME("pixCloseGray3");

    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 (pixGetColormap(pixs))
        return (PIX *)ERROR_PTR("pix has colormap", procName, NULL);
    if ((hsize != 1 && hsize != 3) ||
        (vsize != 1 && vsize != 3))
        return (PIX *)ERROR_PTR("invalid size: must be 1 or 3", procName, NULL);

    if (hsize == 1 && vsize == 1)
        return pixCopy(NULL, pixs);

    pixb = pixAddBorderGeneral(pixs, 4, 8, 2, 8, 0);  /* set to min */

    if (vsize == 1) {
        pixt = pixDilateGray3h(pixb);
        pixSetBorderVal(pixt, 4, 8, 2, 8, 255);  /* set to max */
        pixbd = pixErodeGray3h(pixt);
        pixDestroy(&pixt);
    }
    else if (hsize == 1) {
        pixt = pixDilateGray3v(pixb);
        pixSetBorderVal(pixt, 4, 8, 2, 8, 255);
        pixbd = pixErodeGray3v(pixt);
        pixDestroy(&pixt);
    }
    else {  /* vize == hsize == 3 */
        pixt = pixDilateGray3h(pixb);
        pixbd = pixDilateGray3v(pixt);
        pixDestroy(&pixt);
        pixSetBorderVal(pixbd, 4, 8, 2, 8, 255);
        pixt = pixErodeGray3h(pixbd);
        pixDestroy(&pixbd);
        pixbd = pixErodeGray3v(pixt);
        pixDestroy(&pixt);
    }

    pixd = pixRemoveBorderGeneral(pixbd, 4, 8, 2, 8);
    pixDestroy(&pixb);
    pixDestroy(&pixbd);
    return pixd;
}

/*!
 *  pixBilateral()
 *
 *      Input:  pixs (8 bpp gray or 32 bpp rgb, no colormap)
 *              spatial_stdev  (of gaussian kernel; in pixels, > 0.5)
 *              range_stdev  (of gaussian range kernel; > 5.0; typ. 50.0)
 *              ncomps (number of intermediate sums J(k,x); in [4 ... 30])
 *              reduction  (1, 2 or 4)
 *      Return: pixd (bilateral filtered image), or null on error
 *
 *  Notes:
 *      (1) This performs a relatively fast, separable bilateral
 *          filtering operation.  The time is proportional to ncomps
 *          and varies inversely approximately as the cube of the
 *          reduction factor.  See bilateral.h for algorithm details.
 *      (2) We impose minimum values for range_stdev and ncomps to
 *          avoid nasty artifacts when either are too small.  We also
 *          impose a constraint on their product:
 *               ncomps * range_stdev >= 100.
 *          So for values of range_stdev >= 25, ncomps can be as small as 4.
 *          Here is a qualitative, intuitive explanation for this constraint.
 *          Call the difference in k values between the J(k) == 'delta', where
 *              'delta' ~ 200 / ncomps
 *          Then this constraint is roughly equivalent to the condition:
 *              'delta' < 2 * range_stdev
 *          Note that at an intensity difference of (2 * range_stdev), the
 *          range part of the kernel reduces the effect by the factor 0.14.
 *          This constraint requires that we have a sufficient number of
 *          PCBs (i.e, a small enough 'delta'), so that for any value of
 *          image intensity I, there exists a k (and a PCB, J(k), such that
 *              |I - k| < range_stdev
 *          Any fewer PCBs and we don't have enough to support this condition.
 *      (3) The upper limit of 30 on ncomps is imposed because the
 *          gain in accuracy is not worth the extra computation.
 *      (4) The size of the gaussian kernel is twice the spatial_stdev
 *          on each side of the origin.  The minimum value of
 *          spatial_stdev, 0.5, is required to have a finite sized
 *          spatial kernel.  In practice, a much larger value is used.
 *      (5) Computation of the intermediate images goes inversely
 *          as the cube of the reduction factor.  If you can use a
 *          reduction of 2 or 4, it is well-advised.
 *      (6) The range kernel is defined over the absolute value of pixel
 *          grayscale differences, and hence must have size 256 x 1.
 *          Values in the array represent the multiplying weight
 *          depending on the absolute gray value difference between
 *          the source pixel and the neighboring pixel, and should
 *          be monotonically decreasing.
 *      (7) Interesting observation.  Run this on prog/fish24.jpg, with
 *          range_stdev = 60, ncomps = 6, and spatial_dev = {10, 30, 50}.
 *          As spatial_dev gets larger, we get the counter-intuitive
 *          result that the body of the red fish becomes less blurry.
 */
PIX *
pixBilateral(PIX       *pixs,
             l_float32  spatial_stdev,
             l_float32  range_stdev,
             l_int32    ncomps,
             l_int32    reduction)
{
l_int32       d;
l_float32     sstdev;  /* scaled spatial stdev */
PIX          *pixt, *pixr, *pixg, *pixb, *pixd;

    PROCNAME("pixBilateral");

    if (!pixs || pixGetColormap(pixs))
        return (PIX *)ERROR_PTR("pixs not defined or cmapped", procName, NULL);
    d = pixGetDepth(pixs);
    if (d != 8 && d != 32)
        return (PIX *)ERROR_PTR("pixs not 8 or 32 bpp", procName, NULL);
    if (reduction != 1 && reduction != 2 && reduction != 4)
        return (PIX *)ERROR_PTR("reduction invalid", procName, NULL);
    sstdev = spatial_stdev / (l_float32)reduction;  /* reduced spat. stdev */
    if (sstdev < 0.5)
        return (PIX *)ERROR_PTR("sstdev < 0.5", procName, NULL);
    if (range_stdev <= 5.0)
        return (PIX *)ERROR_PTR("range_stdev <= 5.0", procName, NULL);
    if (ncomps < 4 || ncomps > 30)
        return (PIX *)ERROR_PTR("ncomps not in [4 ... 30]", procName, NULL);
    if (ncomps *