enhance.c

Go to the documentation of this file.

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

/*
 *  enhance.c
 *
 *      Gamma TRC (tone reproduction curve) mapping
 *           PIX     *pixGammaTRC()
 *           PIX     *pixGammaTRCMasked()
 *           PIX     *pixGammaTRCWithAlpha()
 *           NUMA    *numaGammaTRC()
 *
 *      Contrast enhancement
 *           PIX     *pixContrastTRC()
 *           PIX     *pixContrastTRCMasked()
 *           NUMA    *numaContrastTRC()
 *
 *      Histogram equalization
 *           PIX     *pixEqualizeTRC()
 *           NUMA    *numaEqualizeTRC()
 *
 *      Generic TRC mapper
 *           PIX     *pixTRCMap()
 *
 *      Unsharp-masking
 *           PIX     *pixUnsharpMasking()
 *           PIX     *pixUnsharpMaskingGray()
 *           PIX     *pixUnsharpMaskingFast()
 *           PIX     *pixUnsharpMaskingGrayFast()
 *           PIX     *pixUnsharpMaskingGray1D()
 *           PIX     *pixUnsharpMaskingGray2D()
 *
 *      Hue and saturation modification
 *           PIX     *pixModifyHue()
 *           PIX     *pixModifySaturation()
 *           l_int32  pixMeasureSaturation()
 *
 *      General multiplicative constant color transform
 *           PIX     *pixMultConstantColor()
 *           PIX     *pixMultMatrixColor()
 *
 *      Edge by bandpass
 *           PIX     *pixHalfEdgeByBandpass()
 *
 *      Gamma correction, contrast enhancement and histogram equalization
 *      apply a simple mapping function to each pixel (or, for color
 *      images, to each sample (i.e., r,g,b) of the pixel).
 *
 *       - Gamma correction either lightens the image or darkens
 *         it, depending on whether the gamma factor is greater
 *         or less than 1.0, respectively.
 *
 *       - Contrast enhancement darkens the pixels that are already
 *         darker than the middle of the dynamic range (128)
 *         and lightens pixels that are lighter than 128.
 *
 *       - Histogram equalization remaps to have the same number
 *         of image pixels at each of 256 intensity values.  This is
 *         a quick and dirty method of adjusting contrast and brightness
 *         to bring out details in both light and dark regions.
 *
 *      Unsharp masking is a more complicated enhancement.
 *      A "high frequency" image, generated by subtracting
 *      the smoothed ("low frequency") part of the image from
 *      itself, has all the energy at the edges.  This "edge image"
 *      has 0 average value.  A fraction of the edge image is
 *      then added to the original, enhancing the differences
 *      between pixel values at edges.  Because we represent
 *      images as l_uint8 arrays, we preserve dynamic range and
 *      handle negative values by doing all the arithmetic on
 *      shifted l_uint16 arrays; the l_uint8 values are recovered
 *      at the end.
 *
 *      Hue and saturation modification work in HSV space.  Because
 *      this is too large for efficient table lookup, each pixel value
 *      is transformed to HSV, modified, and transformed back.
 *      It's not the fastest way to do this, but the method is
 *      easily understood.
 *
 *      Unsharp masking is never in-place, and returns a clone if no
 *      operation is to be performed.
 */


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

    /* Scales contrast enhancement factor to have a useful range
     * between 0.0 and 1.0 */
static const l_float32  ENHANCE_SCALE_FACTOR = 5.;

    /* Default number of pixels sampled to determine histogram */
static const l_int32  DEFAULT_HISTO_SAMPLES = 100000;


/*-------------------------------------------------------------*
 *         Gamma TRC (tone reproduction curve) mapping         *
 *-------------------------------------------------------------*/
/*!
 *  pixGammaTRC()
 *
 *      Input:  pixd (<optional> null or equal to pixs)
 *              pixs (8 or 32 bpp; or 2, 4 or 8 bpp with colormap)
 *              gamma (gamma correction; must be > 0.0)
 *              minval  (input value that gives 0 for output; can be < 0)
 *              maxval  (input value that gives 255 for output; can be > 255)
 *      Return: pixd always
 *
 *  Notes:
 *      (1) pixd must either be null or equal to pixs.
 *          For in-place operation, set pixd == pixs:
 *             pixGammaTRC(pixs, pixs, ...);
 *          To get a new image, set pixd == null:
 *             pixd = pixGammaTRC(NULL, pixs, ...);
 *      (2) If pixs is colormapped, the colormap is transformed,
 *          either in-place or in a copy of pixs.
 *      (3) We use a gamma mapping between minval and maxval.
 *      (4) If gamma < 1.0, the image will appear darker;
 *          if gamma > 1.0, the image will appear lighter;
 *      (5) If gamma = 1.0 and minval = 0 and maxval = 255, no
 *          enhancement is performed; return a copy unless in-place,
 *          in which case this is a no-op.
 *      (6) For color images that are not colormapped, the mapping
 *          is applied to each component.
 *      (7) minval and maxval are not restricted to the interval [0, 255].
 *          If minval < 0, an input value of 0 is mapped to a
 *          nonzero output.  This will turn black to gray.
 *          If maxval > 255, an input value of 255 is mapped to
 *          an output value less than 255.  This will turn
 *          white (e.g., in the background) to gray.
 *      (8) Increasing minval darkens the image.
 *      (9) Decreasing maxval bleaches the image.
 *      (10) Simultaneously increasing minval and decreasing maxval
 *           will darken the image and make the colors more intense;
 *           e.g., minval = 50, maxval = 200.
 *      (11) See numaGammaTRC() for further examples of use.
 */
PIX *
pixGammaTRC(PIX       *pixd,
            PIX       *pixs,
            l_float32  gamma,
            l_int32    minval,
            l_int32    maxval)
{
l_int32   d;
NUMA     *nag;
PIXCMAP  *cmap;

    PROCNAME("pixGammaTRC");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, pixd);
    if (pixd && (pixd != pixs))
        return (PIX *)ERROR_PTR("pixd not null or pixs", procName, pixd);
    if (gamma <= 0.0) {
        L_WARNING("gamma must be > 0.0; setting to 1.0", procName);
        gamma = 1.0;
    }
    if (minval >= maxval)
        return (PIX *)ERROR_PTR("minval not < maxval", procName, pixd);
    cmap = pixGetColormap(pixs);
    d = pixGetDepth(pixs);
    if (!cmap && d != 8 && d != 32)
        return (PIX *)ERROR_PTR("depth not 8 or 32 bpp", procName, pixd);

    if (gamma == 1.0 && minval == 0 && maxval == 255)  /* no-op */
        return pixCopy(pixd, pixs);

    if (!pixd)  /* start with a copy if not in-place */
        pixd = pixCopy(NULL, pixs);

    if (cmap) {
        pixcmapGammaTRC(pixGetColormap(pixd), gamma, minval, maxval);
        return pixd;
    }

        /* pixd is 8 or 32 bpp */
    if ((nag = numaGammaTRC(gamma, minval, maxval)) == NULL)
        return (PIX *)ERROR_PTR("nag not made", procName, pixd);
    pixTRCMap(pixd, NULL, nag);
    numaDestroy(&nag);

    return pixd;
}


/*!
 *  pixGammaTRCMasked()
 *
 *      Input:  pixd (<optional> null or equal to pixs)
 *              pixs (8 or 32 bpp; not colormapped)
 *              pixm (<optional> null or 1 bpp)
 *              gamma (gamma correction; must be > 0.0)
 *              minval  (input value that gives 0 for output; can be < 0)
 *              maxval  (input value that gives 255 for output; can be > 255)
 *      Return: pixd always
 *
 *  Notes:
 *      (1) Same as pixGammaTRC() except mapping is optionally over
 *          a subset of pixels described by pixm.
 *      (2) Masking does not work for colormapped images.
 *      (3) See pixGammaTRC() for details on how to use the parameters.
 */
PIX *
pixGammaTRCMasked(PIX       *pixd,
                  PIX       *pixs,
                  PIX       *pixm,
                  l_float32  gamma,
                  l_int32    minval,
                  l_int32    maxval)
{
l_int32  d;
NUMA    *nag;

    PROCNAME("pixGammaTRCMasked");

    if (!pixm)
        return pixGammaTRC(pixd, pixs, gamma, minval, maxval);

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, pixd);
    if (pixGetColormap(pixs))
        return (PIX *)ERROR_PTR("invalid: pixs has a colormap", procName, pixd);
    if (pixd && (pixd != pixs))
        return (PIX *)ERROR_PTR("pixd not null or pixs", procName, pixd);
    d = pixGetDepth(pixs);
    if (d != 8 && d != 32)
        return (PIX *)ERROR_PTR("depth not 8 or 32 bpp", procName, pixd);
    if (minval >= maxval)
        return (PIX *)ERROR_PTR("minval not < maxval", procName, pixd);
    if (gamma <= 0.0) {
        L_WARNING("gamma must be > 0.0; setting to 1.0", procName);
        gamma = 1.0;
    }

    if (gamma == 1.0 && minval == 0 && maxval == 255)
        return pixCopy(pixd, pixs);

    if (!pixd)  /* start with a copy if not in-place */
        pixd = pixCopy(NULL, pixs);

    if ((nag = numaGammaTRC(gamma, minval, maxval)) == NULL)
        return (PIX *)ERROR_PTR("nag not made", procName, pixd);
    pixTRCMap(pixd, pixm, nag);
    numaDestroy(&nag);

    return pixd;
}


/*!
 *  pixGammaTRCWithAlpha()
 *
 *      Input:  pixd (<optional> null or equal to pixs)
 *              pixs (32 bpp)
 *              gamma (gamma correction; must be > 0.0)
 *              minval  (input value that gives 0 for output; can be < 0)
 *              maxval  (input value that gives 255 for output; can be > 255)
 *      Return: pixd always
 *
 *  Notes:
 *      (1) See usage notes in pixGammaTRC().
 *      (2) This version saves the alpha channel.  It is only valid
 *          for 32 bpp (no colormap), and is a bit slower.
 */
PIX *
pixGammaTRCWithAlpha(PIX       *pixd,
                     PIX       *pixs,
                     l_float32  gamma,
                     l_int32    minval,
                     l_int32    maxval)
{
NUMA  *nag;
PIX   *pixalpha;

    PROCNAME("pixGammaTRCWithAlpha");

    if (!pixs || pixGetDepth(pixs) != 32)
        return (PIX *)ERROR_PTR("pixs undefined or not 32 bpp", procName, pixd);
    if (pixd && (pixd != pixs))
        return (PIX *)ERROR_PTR("pixd not null or pixs", procName, pixd);
    if (gamma <= 0.0) {
        L_WARNING("gamma must be > 0.0; setting to 1.0", procName);
        gamma = 1.0;
    }
    if (minval >= maxval)
        return (PIX *)ERROR_PTR("minval not < maxval", procName, pixd);

    if (gamma == 1.0 && minval == 0 && maxval == 255)
        return pixCopy(pixd, pixs);

    pixalpha = pixGetRGBComponent(pixs, L_ALPHA_CHANNEL);  /* save */
    if (!pixd)  /* start with a copy if not in-place */
        pixd = pixCopy(NULL, pixs);

    if ((nag = numaGammaTRC(gamma, minval, maxval)) == NULL)
        return (PIX *)ERROR_PTR("nag not made", procName, pixd);
    pixTRCMap(pixd, NULL, nag);
    pixSetRGBComponent(pixd, pixalpha, L_ALPHA_CHANNEL);  /* restore */

    numaDestroy(&nag);
    pixDestroy(&pixalpha);
    return pixd;
}


/*!
 *  numaGammaTRC()
 *
 *      Input:  gamma   (gamma factor; must be > 0.0)
 *              minval  (input value that gives 0 for output)
 *              maxval  (input value that gives 255 for output)
 *      Return: na, or null on error
 *
 *  Notes:
 *      (1) The map is returned as a numa; values are clipped to [0, 255].
 *      (2) To force all intensities into a range within fraction delta
 *          of white, use: minval = -256 * (1 - delta) / delta
 *                         maxval = 255
 *      (3) To force all intensities into a range within fraction delta
 *          of black, use: minval = 0
 *                         maxval = 256 * (1 - delta) / delta
 */
NUMA *
numaGammaTRC(l_float32  gamma,
             l_int32    minval,
             l_int32    maxval)
{
l_int32    i, val;
l_float32  x, invgamma;
NUMA      *na;

    PROCNAME("numaGammaTRC");

    if (minval >= maxval)
        return (NUMA *)ERROR_PTR("minval not < maxval", procName, NULL);
    if (gamma <= 0.0) {
        L_WARNING("gamma must be > 0.0; setting to 1.0", procName);
        gamma = 1.0;
    }

    invgamma = 1. / gamma;
    na = numaCreate(256);
    for (i = 0; i < minval; i++)
        numaAddNumber(na, 0);
    for (i = minval; i <= maxval; i++) {
        if (i < 0) continue;
        if (i > 255) continue;
        x = (l_float32)(i - minval) / (l_float32)(maxval - minval);
        val = (l_int32)(255. * powf(x, invgamma) + 0.5);
        val = L_MAX(val, 0);
        val = L_MIN(val, 255);
        numaAddNumber(na, val);
    }
    for (i = maxval + 1; i < 256; i++)
        numaAddNumber(na, 255);

    return na;
}


/*-------------------------------------------------------------*
 *                      Contrast enhancement                   *
 *-------------------------------------------------------------*/
/*!
 *  pixContrastTRC()
 *
 *      Input:  pixd (<optional> null or equal to pixs)
 *              pixs (8 or 32 bpp; or 2, 4 or 8 bpp with colormap)
 *              factor  (0.0 is no enhancement)
 *      Return: pixd always
 *
 *  Notes:
 *      (1) pixd must either be null or equal to pixs.
 *          For in-place operation, set pixd == pixs:
 *             pixContrastTRC(pixs, pixs, ...);
 *          To get a new image, set pixd == null:
 *             pixd = pixContrastTRC(NULL, pixs, ...);
 *      (2) If pixs is colormapped, the colormap is transformed,
 *          either in-place or in a copy of pixs.
 *      (3) Contrast is enhanced by mapping each color component
 *          using an atan function with maximum slope at 127.
 *          Pixels below 127 are lowered in intensity and pixels
 *          above 127 are increased.
 *      (4) The useful range for the contrast factor is scaled to
 *          be in (0.0 to 1.0), but larger values can also be used.
 *      (5) If factor == 0.0, no enhancement is performed; return a copy
 *          unless in-place, in which case this is a no-op.
 *      (6) For color images that are not colormapped, the mapping
 *          is applied to each component.
 */
PIX *
pixContrastTRC(PIX       *pixd,
               PIX       *pixs,
               l_float32  factor)
{
l_int32   d;
NUMA     *nac;
PIXCMAP  *cmap;

    PROCNAME("pixContrastTRC");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, pixd);
    if (pixd && (pixd != pixs))
        return (PIX *)ERROR_PTR("pixd not null or pixs", procName, pixd);
    if (factor < 0.0) {
        L_WARNING("factor must be >= 0.0; using 0.0", procName);
        factor = 0.0;
    }
    if (factor == 0.0)
        return pixCopy(pixd, pixs);

    cmap = pixGetColormap(pixs);
    d = pixGetDepth(pixs);
    if (!cmap && d != 8 && d != 32)
        return (PIX *)ERROR_PTR("depth not 8 or 32 bpp", procName, pixd);

    if (!pixd)  /* start with a copy if not in-place */
        pixd = pixCopy(NULL, pixs);

    if (cmap) {
        pixcmapContrastTRC(pixGetColormap(pixd), factor);
        return pixd;
    }

        /* pixd is 8 or 32 bpp */
    if ((nac = numaContrastTRC(factor)) == NULL)
        return (PIX *)ERROR_PTR("nac not made", procName, pixd);
    pixTRCMap(pixd, NULL, nac);
    numaDestroy(&nac);

    return pixd;
}


/*!
 *  pixContrastTRCMasked()
 *
 *      Input:  pixd (<optional> null or equal to pixs)
 *              pixs (8 or 32 bpp; or 2, 4 or 8 bpp with colormap)
 *              pixm (<optional> null or 1 bpp)
 *              factor  (0.0 is no enhancement)
 *      Return: pixd always
 *
 *  Notes:
 *      (1) Same as pixContrastTRC() except mapping is optionally over
 *          a subset of pixels described by pixm.
 *      (2) Masking does not work for colormapped images.
 *      (3) See pixContrastTRC() for details on how to use the parameters.
 */
PIX *
pixContrastTRCMasked(PIX       *pixd,
                     PIX       *pixs,
                     PIX       *pixm,
                     l_float32  factor)
{
l_int32  d;
NUMA    *nac;

    PROCNAME("pixContrastTRCMasked");

    if (!pixm)
        return pixContrastTRC(pixd, pixs, factor);

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, pixd);
    if (pixGetColormap(pixs))
        return (PIX *)ERROR_PTR("invalid: pixs has a colormap", procName, pixd);
    if (pixd && (pixd != pixs))
        return (PIX *)ERROR_PTR("pixd not null or pixs", procName, pixd);
    d = pixGetDepth(pixs);
    if (d != 8 && d != 32)
        return (PIX *)ERROR_PTR("depth not 8 or 32 bpp", procName, pixd);

    if (factor < 0.0) {
        L_WARNING("factor must be >= 0.0; using 0.0", procName);
        factor = 0.0;
    }
    if (factor == 0.0)
        return pixCopy(pixd, pixs);

    if (!pixd)  /* start with a copy if not in-place */
        pixd = pixCopy(NULL, pixs);

    if ((nac = numaContrastTRC(factor)) == NULL)
        return (PIX *)ERROR_PTR("nac not made", procName, pixd);
    pixTRCMap(pixd, pixm, nac);
    numaDestroy(&nac);

    return pixd;
}


/*!
 *  numaContrastTRC()
 *
 *      Input:  factor (generally between 0.0 (no enhancement)
 *              and 1.0, but can be larger than 1.0)
 *      Return: na, or null on error
 *
 *  Notes:
 *      (1) The mapping is monotonic increasing, where 0 is mapped
 *          to 0 and 255 is mapped to 255. 
 *      (2) As 'factor' is increased from 0.0 (where the mapping is linear),
 *          the map gets closer to its limit as a step function that
 *          jumps from 0 to 255 at the center (input value = 127).
 */
NUMA *
numaContrastTRC(l_float32  factor)
{
l_int32    i, val;
l_float64  x, ymax, ymin, dely, scale;
NUMA      *na;

    PROCNAME("numaContrastTRC");

    if (factor < 0.0) {
        L_WARNING("factor must be >= 0.0; using 0.0; no enhancement", procName);
        factor = 0.0;
    }
    if (factor == 0.0)
        return numaMakeSequence(0, 1, 256);  /* linear map */

    scale = ENHANCE_SCALE_FACTOR;
    ymax = atan((l_float64)(1.0 * factor * scale));
    ymin = atan((l_float64)(-127. * factor * scale / 128.));
    dely = ymax - ymin;
    na = numaCreate(256);
    for (i = 0; i < 256; i++) {
        x = (l_float64)i;
        val = (l_int32)((255. / dely) *
             (-ymin + atan((l_float64)(factor * scale * (x - 127.) / 128.))) +
                 0.5);
        numaAddNumber(na, val);
    }

    return na;
}


/*-------------------------------------------------------------*
 *                     Histogram equalization                  *
 *-------------------------------------------------------------*/
/*!
 *  pixEqualizeTRC()
 *
 *      Input:  pixd (<optional> null or equal to pixs)
 *              pixs (8 bpp gray, 32 bpp rgb, or colormapped)
 *              fract (fraction of equalization movement of pixel values)
 *              factor (subsampling factor; integer >= 1)
 *      Return: pixd, or null on error
 *
 *  Notes:
 *      (1) pixd must either be null or equal to pixs.
 *          For in-place operation, set pixd == pixs:
 *             pixEqualizeTRC(pixs, pixs, ...);
 *          To get a new image, set pixd == null:
 *             pixd = pixEqualizeTRC(NULL, pixs, ...);
 *      (2) In histogram equalization, a tone reproduction curve
 *          mapping is used to make the number of pixels at each
 *          intensity equal.
 *      (3) If fract == 0.0, no equalization is performed; return a copy
 *          unless in-place, in which case this is a no-op.
 *          If fract == 1.0, equalization is complete.
 *      (4) Set the subsampling factor > 1 to reduce the amount of computation.
 *      (5) If pixs is colormapped, the colormap is removed and
 *          converted to rgb or grayscale.
 *      (6) If pixs has color, equalization is done in each channel
 *          separately.
 *      (7) Note that even if there is a colormap, we can get an
 *          in-place operation because the intermediate image pixt
 *          is copied back to pixs (which for in-place is the same
 *          as pixd).
 */
PIX *
pixEqualizeTRC(PIX       *pixd,
               PIX       *pixs,
               l_float32  fract,
               l_int32    factor)
{
l_int32   d;
NUMA     *na;
PIX      *pixt, *pix8;
PIXCMAP  *cmap;

    PROCNAME("pixEqualizeTRC");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (pixd && (pixd != pixs))
        return (PIX *)ERROR_PTR("pixd not null or pixs", procName, pixd);
    cmap = pixGetColormap(pixs);
    d = pixGetDepth(pixs);
    if (d != 8 && d != 32 && !cmap)
        return (PIX *)ERROR_PTR("pixs not 8/32 bpp or cmapped", procName, NULL);
    if (fract < 0.0 || fract > 1.0)
        return (PIX *)ERROR_PTR("fract not in [0.0 ... 1.0]", procName, NULL);
    if (factor < 1)
        return (PIX *)ERROR_PTR("sampling factor < 1", procName, NULL);

    if (fract == 0.0)
        return pixCopy(pixd, pixs);
 
        /* If there is a colormap, remove it. */
    if (cmap)
        pixt = pixRemoveColormap(pixs, REMOVE_CMAP_BASED_ON_SRC);
    else
        pixt = pixClone(pixs);

        /* Make a copy if necessary */
    pixd = pixCopy(pixd, pixt);
    pixDestroy(&pixt);

    d = pixGetDepth(pixd);
    if (d == 8) {
        na = numaEqualizeTRC(pixd, fract, factor);
        pixTRCMap(pixd, NULL, na);
        numaDestroy(&na);
    }
    else {  /* 32 bpp */
        pix8 = pixGetRGBComponent(pixd, COLOR_RED);
        na = numaEqualizeTRC(pix8, fract, factor);
        pixTRCMap(pix8, NULL, na);
        pixSetRGBComponent(pixd, pix8, COLOR_RED);
        numaDestroy(&na);
        pixDestroy(&pix8);
        pix8 = pixGetRGBComponent(pixd, COLOR_GREEN);
        na = numaEqualizeTRC(pix8, fract, factor);
        pixTRCMap(pix8, NULL, na);
        pixSetRGBComponent(pixd, pix8, COLOR_GREEN);
        numaDestroy(&na);
        pixDestroy(&pix8);
        pix8 = pixGetRGBComponent(pixd, COLOR_BLUE);
        na = numaEqualizeTRC(pix8, fract, factor);
        pixTRCMap(pix8, NULL, na);
        pixSetRGBComponent(pixd, pix8, COLOR_BLUE);
        numaDestroy(&na);
        pixDestroy(&pix8);
    }

    return pixd;
}


/*!
 *  numaEqualizeTRC()
 *
 *      Input:  pix (8 bpp, no colormap)
 *              fract (fraction of equalization movement of pixel values)
 *              factor (subsampling factor; integer >= 1)
 *      Return: nad, or null on error
 *
 *  Notes:
 *      (1) If fract == 0.0, no equalization will be performed.
 *          If fract == 1.0, equalization is complete.
 *      (2) Set the subsampling factor > 1 to reduce the amount of computation.
 *      (3) The map is returned as a numa with 256 values, specifying
 *          the equalized value (array value) for every input value
 *          (the array index).
 */
NUMA *
numaEqualizeTRC(PIX       *pix,
                l_float32  fract,
                l_int32    factor)
{
l_int32    iin, iout, itarg;
l_float32  val, sum;
NUMA      *nah, *nasum, *nad;

    PROCNAME("numaEqualizeTRC");

    if (!pix)
        return (NUMA *)ERROR_PTR("pix not defined", procName, NULL);
    if (pixGetDepth(pix) != 8)
        return (NUMA *)ERROR_PTR("pix not 8 bpp", procName, NULL);
    if (fract < 0.0 || fract > 1.0)
        return (NUMA *)ERROR_PTR("fract not in [0.0 ... 1.0]", procName, NULL);
    if (factor < 1)
        return (NUMA *)ERROR_PTR("sampling factor < 1", procName, NULL);

    if (fract == 0.0)
        L_WARNING("fract = 0.0; no equalization requested", procName);

    if ((nah = pixGetGrayHistogram(pix, factor)) == NULL)
        return (NUMA *)ERROR_PTR("histogram not made", procName, NULL);
    numaGetSum(nah, &sum);
    nasum = numaGetPartialSums(nah);

    nad = numaCreate(256);
    for (iin = 0; iin < 256; iin++) {
        numaGetFValue(nasum, iin, &val);
        itarg = (l_int32)(255. * val / sum + 0.5);
        iout = iin + (l_int32)(fract * (itarg - iin));
        iout = L_MIN(iout, 255);  /* to be safe */
        numaAddNumber(nad, iout);
    }

    numaDestroy(&nah);
    numaDestroy(&nasum);
    return nad;
}


/*-------------------------------------------------------------*
 *                       Generic TRC mapping                   *
 *-------------------------------------------------------------*/
/*!
 *  pixTRCMap()
 *
 *      Input:  pixs (8 grayscale or 32 bpp rgb; not colormapped)
 *              pixm (<optional> 1 bpp mask)
 *              na (mapping array)
 *      Return: pixd, or null on error
 *
 *  Notes:
 *      (1) This operation is in-place on pixs.
 *      (2) For 32 bpp, this applies the same map to each of the r,g,b
 *          components.
 *      (3) The mapping array is of size 256, and it maps the input
 *          index into values in the range [0, 255].
 *      (4) If defined, the optional 1 bpp mask pixm has its origin
 *          aligned with pixs, and the map function is applied only
 *          to pixels in pixs under the fg of pixm.
 *      (5) For 32 bpp, this does not save the alpha channel.
 */
l_int32
pixTRCMap(PIX   *pixs,
          PIX   *pixm,
          NUMA  *na)
{
l_int32    w, h, d, wm, hm, wpl, wplm, i, j, sval8, dval8;
l_int32   *tab;
l_uint32   sval32, dval32;
l_uint32  *data, *datam, *line, *linem;

    PROCNAME("pixTRCMap");

    if (!pixs)
        return ERROR_INT("pixs not defined", procName, 1);
    if (pixGetColormap(pixs))
        return ERROR_INT("pixs is colormapped", procName, 1);
    if (!na)
        return ERROR_INT("na not defined", procName, 1);
    if (numaGetCount(na) != 256)
        return ERROR_INT("na not of size 256", procName, 1);
    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 8 && d != 32)
        return ERROR_INT("pixs not 8 or 32 bpp", procName, 1);
    if (pixm) {
        if (pixGetDepth(pixm) != 1)
            return ERROR_INT("pixm not 1 bpp", procName, 1);
    }

    tab = numaGetIArray(na);  /* get the array for efficiency */
    wpl = pixGetWpl(pixs);
    data = pixGetData(pixs);
    if (!pixm) {
        if (d == 8) {
            for (i = 0; i < h; i++) {
                line = data + i * wpl;
                for (j = 0; j < w; j++) {
                    sval8 = GET_DATA_BYTE(line, j);
                    dval8 = tab[sval8];
                    SET_DATA_BYTE(line, j, dval8);
                }
            }
        }
        else {  /* d == 32 */
            for (i = 0; i < h; i++) {
                line = data + i * wpl;
                for (j = 0; j < w; j++) {
                    sval32 = *(line + j);
                    dval32 =
                        tab[(sval32 >> L_RED_SHIFT) & 0xff] << L_RED_SHIFT |
                        tab[(sval32 >> L_GREEN_SHIFT) & 0xff] << L_GREEN_SHIFT |
                        tab[(sval32 >> L_BLUE_SHIFT) & 0xff] << L_BLUE_SHIFT;
                    *(line + j) = dval32;
                }
            }
        }
    }
    else {
        datam = pixGetData(pixm);
        wplm = pixGetWpl(pixm);
        pixGetDimensions(pixm, &wm, &hm, NULL);
        if (d == 8) {
            for (i = 0; i < h; i++) {
                if (i >= hm)
                    break;
                line = data + i * wpl;
                linem = datam + i * wplm;
                for (j = 0; j < w; j++) {
                    if (j >= wm)
                        break;
                    if (GET_DATA_BIT(linem, j) == 0)
                        continue;
                    sval8 = GET_DATA_BYTE(line, j);
                    dval8 = tab[sval8];
                    SET_DATA_BYTE(line, j, dval8);
                }
            }
        }
        else {  /* d == 32 */
            for (i = 0; i < h; i++) {
                if (i >= hm)
                    break;
                line = data + i * wpl;
                linem = datam + i * wplm;
                for (j = 0; j < w; j++) {
                    if (j >= wm)
                        break;
                    if (GET_DATA_BIT(linem, j) == 0)
                        continue;
                    sval32 = *(line + j);
                    dval32 =
                        tab[(sval32 >> L_RED_SHIFT) & 0xff] << L_RED_SHIFT |
                        tab[(sval32 >> L_GREEN_SHIFT) & 0xff] << L_GREEN_SHIFT |
                        tab[(sval32 >> L_BLUE_SHIFT) & 0xff] << L_BLUE_SHIFT;
                    *(line + j) = dval32;
                }
            }
        }
    }

    FREE(tab);
    return 0;
}



/*-----------------------------------------------------------------------*
 *                             Unsharp masking                           *
 *-----------------------------------------------------------------------*/
/*!
 *  pixUnsharpMasking()
 *
 *      Input:  pixs (all depths except 1 bpp; with or without colormaps)
 *              halfwidth  ("half-width" of smoothing filter)
 *              fract  (fraction of edge added back into image)
 *      Return: pixd, or null on error
 *
 *  Notes:
 *      (1) We use symmetric smoothing filters of odd dimension,
 *          typically use sizes of 3, 5, 7, etc.  The @halfwidth parameter
 *          for these is (size - 1)/2; i.e., 1, 2, 3, etc.
 *      (2) The fract parameter is typically taken in the
 *          range:  0.2 < fract < 0.7
 *      (3) Returns a clone if no sharpening is requested.
 */
PIX *
pixUnsharpMasking(PIX       *pixs,
                  l_int32    halfwidth,
                  l_float32  fract)
{
l_int32  d;
PIX     *pixt, *pixd, *pixr, *pixrs, *pixg, *pixgs, *pixb, *pixbs;

    PROCNAME("pixUnsharpMasking");

    if (!pixs || (pixGetDepth(pixs) == 1))
        return (PIX *)ERROR_PTR("pixs not defined or 1 bpp", procName, NULL);
    if (fract <= 0.0 || halfwidth <= 0) {
        L_WARNING("no sharpening requested; clone returned", procName);
        return pixClone(pixs);
    }

    if (halfwidth == 1 || halfwidth == 2)
        return pixUnsharpMaskingFast(pixs, halfwidth, fract, L_BOTH_DIRECTIONS);

        /* Remove colormap; clone if possible; result is either 8 or 32 bpp */
    if ((pixt = pixConvertTo8Or32(pixs, 0, 1)) == NULL)
        return (PIX *)ERROR_PTR("pixt not made", procName, NULL);

        /* Sharpen */
    d = pixGetDepth(pixt);
    if (d == 8)
        pixd = pixUnsharpMaskingGray(pixt, halfwidth, fract);
    else {  /* d == 32 */
        pixr = pixGetRGBComponent(pixs, COLOR_RED);
        pixrs = pixUnsharpMaskingGray(pixr, halfwidth, fract);
        pixDestroy(&pixr);
        pixg = pixGetRGBComponent(pixs, COLOR_GREEN);
        pixgs = pixUnsharpMaskingGray(pixg, halfwidth, fract);
        pixDestroy(&pixg);
        pixb = pixGetRGBComponent(pixs, COLOR_BLUE);
        pixbs = pixUnsharpMaskingGray(pixb, halfwidth, fract);
        pixDestroy(&pixb);
        pixd = pixCreateRGBImage(pixrs, pixgs, pixbs);
        pixDestroy(&pixrs);
        pixDestroy(&pixgs);
        pixDestroy(&pixbs);
    }

    pixDestroy(&pixt);
    return pixd;
}


/*!
 *  pixUnsharpMaskingGray()
 *
 *      Input:  pixs (8 bpp; no colormap)
 *              halfwidth  ("half-width" of smoothing filter)
 *              fract  (fraction of edge added back into image)
 *      Return: pixd, or null on error
 *
 *  Notes:
 *      (1) We use symmetric smoothing filters of odd dimension,
 *          typically use sizes of 3, 5, 7, etc.  The @halfwidth parameter
 *          for these is (size - 1)/2; i.e., 1, 2, 3, etc.
 *      (2) The fract parameter is typically taken in the range:
 *          0.2 < fract < 0.7
 *      (3) Returns a clone if no sharpening is requested.
 */
PIX *
pixUnsharpMaskingGray(PIX       *pixs,
                      l_int32    halfwidth,
                      l_float32  fract)
{
l_int32  w, h, d;
PIX     *pixc, *pixd;
PIXACC  *pixacc;

    PROCNAME("pixUnsharpMaskingGray");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 8 || pixGetColormap(pixs) != NULL)
        return (PIX *)ERROR_PTR("pixs not 8 bpp or has cmap", procName, NULL);
    if (fract <= 0.0 || halfwidth <= 0) {
        L_WARNING("no sharpening requested; clone returned", procName);
        return pixClone(pixs);
    }
    if (halfwidth == 1 || halfwidth == 2)
        return pixUnsharpMaskingGrayFast(pixs, halfwidth, fract,
                                         L_BOTH_DIRECTIONS);

    if ((pixc = pixBlockconvGray(pixs, NULL, halfwidth, halfwidth)) == NULL)
        return (PIX *)ERROR_PTR("pixc not made", procName, NULL);

        /* Steps:
         *    (1) edge image is pixs - pixc  (this is highpass part)
         *    (2) multiply edge image by fract
         *    (3) add fraction of edge to pixs
         *
         * To show how this is done with both interfaces to arithmetic
         * on integer Pix, here is the implementation in the lower-level
         * function calls:
         *    pixt = pixInitAccumulate(w, h, 0x10000000)) == NULL)
         *    pixAccumulate(pixt, pixs, L_ARITH_ADD);
         *    pixAccumulate(pixt, pixc, L_ARITH_SUBTRACT);
         *    pixMultConstAccumulate(pixt, fract, 0x10000000);
         *    pixAccumulate(pixt, pixs, L_ARITH_ADD);
         *    pixd = pixFinalAccumulate(pixt, 0x10000000, 8)) == NULL)
         *    pixDestroy(&pixt);
         *
         * The code below does the same thing using the Pixacc accumulator,
         * hiding the details of the offset that is needed for subtraction.
         */
    pixacc = pixaccCreate(w, h, 1);
    pixaccAdd(pixacc, pixs);
    pixaccSubtract(pixacc, pixc);
    pixaccMultConst(pixacc, fract);
    pixaccAdd(pixacc, pixs);
    pixd = pixaccFinal(pixacc, 8);
    pixaccDestroy(&pixacc);

    pixDestroy(&pixc);
    return pixd;
}


/*!
 *  pixUnsharpMaskingFast()
 *
 *      Input:  pixs (all depths except 1 bpp; with or without colormaps)
 *              halfwidth  ("half-width" of smoothing filter; 1 and 2 only)
 *              fract  (fraction of high frequency added to image)
 *              direction (L_HORIZ, L_VERT, L_BOTH_DIRECTIONS)
 *      Return: pixd, or null on error
 *
 *  Notes:
 *      (1) The fast version uses separable 1-D filters directly on
 *          the input image.  The halfwidth is either 1 (full width = 3)
 *          or 2 (full width = 5).
 *      (2) The fract parameter is typically taken in the
 *            range:  0.2 < fract < 0.7
 *      (3) To skip horizontal sharpening, use @fracth = 0.0; ditto for @fractv
 *      (4) For one dimensional filtering (as an example):
 *          For @halfwidth = 1, the low-pass filter is
 *              L:    1/3    1/3   1/3
 *          and the high-pass filter is
 *              H = I - L:   -1/3   2/3   -1/3
 *          For @halfwidth = 2, the low-pass filter is
 *              L:    1/5    1/5   1/5    1/5    1/5
 *          and the high-pass filter is
 *              H = I - L:   -1/5  -1/5   4/5  -1/5   -1/5
 *          The new sharpened pixel value is found by adding some fraction
 *          of the high-pass filter value (which sums to 0) to the
 *          initial pixel value:
 *              N = I + fract * H
 *      (5) For 2D, the sharpening filter is not separable, because the
 *          vertical filter depends on the horizontal location relative
 *          to the filter origin, and v.v.   So we either do the full
 *          2D filter (for @halfwidth == 1) or do the low-pass
 *          convolution separably and then compose with the original pix.
 *      (6) Returns a clone if no sharpening is requested.
 */
PIX *
pixUnsharpMaskingFast(PIX       *pixs,
                      l_int32    halfwidth,
                      l_float32  fract,
                      l_int32    direction)
{
l_int32  d;
PIX     *pixt, *pixd, *pixr, *pixrs, *pixg, *pixgs, *pixb, *pixbs;

    PROCNAME("pixUnsharpMaskingFast");

    if (!pixs || (pixGetDepth(pixs) == 1))
        return (PIX *)ERROR_PTR("pixs not defined or 1 bpp", procName, NULL);
    if (fract <= 0.0 || halfwidth <= 0) {
        L_WARNING("no sharpening requested; clone returned", procName);
        return pixClone(pixs);
    }
    if (halfwidth != 1 && halfwidth != 2)
        return (PIX *)ERROR_PTR("halfwidth must be 1 or 2", procName, NULL);
    if (direction != L_HORIZ && direction != L_VERT &&
        direction != L_BOTH_DIRECTIONS)
        return (PIX *)ERROR_PTR("invalid direction", procName, NULL);

        /* Remove colormap; clone if possible; result is either 8 or 32 bpp */
    if ((pixt = pixConvertTo8Or32(pixs, 0, 1)) == NULL)
        return (PIX *)ERROR_PTR("pixt not made", procName, NULL);

        /* Sharpen */
    d = pixGetDepth(pixt);
    if (d == 8)
        pixd = pixUnsharpMaskingGrayFast(pixt, halfwidth, fract, direction);
    else {  /* d == 32 */
        pixr = pixGetRGBComponent(pixs, COLOR_RED);
        pixrs = pixUnsharpMaskingGrayFast(pixr, halfwidth, fract, direction);
        pixDestroy(&pixr);
        pixg = pixGetRGBComponent(pixs, COLOR_GREEN);
        pixgs = pixUnsharpMaskingGrayFast(pixg, halfwidth, fract, direction);
        pixDestroy(&pixg);
        pixb = pixGetRGBComponent(pixs, COLOR_BLUE);
        pixbs = pixUnsharpMaskingGrayFast(pixb, halfwidth, fract, direction);
        pixDestroy(&pixb);
        pixd = pixCreateRGBImage(pixrs, pixgs, pixbs);
        pixDestroy(&pixrs);
        pixDestroy(&pixgs);
        pixDestroy(&pixbs);
    }

    pixDestroy(&pixt);
    return pixd;
}



/*!
 *  pixUnsharpMaskingGrayFast()
 *
 *      Input:  pixs (8 bpp; no colormap)
 *              halfwidth  ("half-width" of smoothing filter: 1 or 2)
 *              fract  (fraction of high frequency added to image)
 *              direction (L_HORIZ, L_VERT, L_BOTH_DIRECTIONS)
 *      Return: pixd, or null on error
 *
 *  Notes:
 *      (1) For usage and explanation of the algorithm, see notes
 *          in pixUnsharpMaskingFast().
 *      (2) Returns a clone if no sharpening is requested.
 */
PIX *
pixUnsharpMaskingGrayFast(PIX       *pixs,
                          l_int32    halfwidth,
                          l_float32  fract,
                          l_int32    direction)
{
PIX  *pixd;

    PROCNAME("pixUnsharpMaskingGrayFast");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (pixGetDepth(pixs) != 8 || pixGetColormap(pixs) != NULL)
        return (PIX *)ERROR_PTR("pixs not 8 bpp or has cmap", procName, NULL);
    if (fract <= 0.0 || halfwidth <= 0) {
        L_WARNING("no sharpening requested; clone returned", procName);
        return pixClone(pixs);
    }
    if (halfwidth != 1 && halfwidth != 2)
        return (PIX *)ERROR_PTR("halfwidth must be 1 or 2", procName, NULL);
    if (direction != L_HORIZ && direction != L_VERT &&
        direction != L_BOTH_DIRECTIONS)
        return (PIX *)ERROR_PTR("invalid direction", procName, NULL);

    if (direction != L_BOTH_DIRECTIONS)
        pixd = pixUnsharpMaskingGray1D(pixs, halfwidth, fract, direction);
    else  /* 2D sharpening */
        pixd = pixUnsharpMaskingGray2D(pixs, halfwidth, fract);

    return pixd;
}


/*!
 *  pixUnsharpMaskingGray1D()
 *
 *      Input:  pixs (8 bpp; no colormap)
 *              halfwidth  ("half-width" of smoothing filter: 1 or 2)
 *              fract  (fraction of high frequency added to image)
 *              direction (of filtering; use L_HORIZ or L_VERT)
 *      Return: pixd, or null on error
 *
 *  Notes:
 *      (1) For usage and explanation of the algorithm, see notes
 *          in pixUnsharpMaskingFast().
 *      (2) Returns a clone if no sharpening is requested.
 */
PIX *
pixUnsharpMaskingGray1D(PIX       *pixs,
                        l_int32    halfwidth,
                        l_float32  fract,
                        l_int32    direction)
{
l_int32    w, h, d, wpls, wpld, i, j, ival;
l_uint32  *datas, *datad;
l_uint32  *lines, *lines0, *lines1, *lines2, *lines3, *lines4, *lined;
l_float32  val, a[5];
PIX       *pixd;

    PROCNAME("pixUnsharpMaskingGray1D");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 8 || pixGetColormap(pixs) != NULL)
        return (PIX *)ERROR_PTR("pixs not 8 bpp or has cmap", procName, NULL);
    if (fract <= 0.0 || halfwidth <= 0) {
        L_WARNING("no sharpening requested; clone returned", procName);
        return pixClone(pixs);
    }
    if (halfwidth != 1 && halfwidth != 2)
        return (PIX *)ERROR_PTR("halfwidth must be 1 or 2", procName, NULL);

        /* Initialize pixd with pixels from pixs that will not be
         * set when computing the sharpened values. */
    pixd = pixCopyBorder(NULL, pixs, halfwidth, halfwidth,
                         halfwidth, halfwidth);
    datas = pixGetData(pixs);
    datad = pixGetData(pixd);
    wpls = pixGetWpl(pixs);
    wpld = pixGetWpl(pixd);

    if (halfwidth == 1) {
        a[0] = -fract / 3.0;
        a[1] = 1.0 + fract * 2.0 / 3.0;
        a[2] = a[0];
    }
    else {  /* halfwidth == 2 */
        a[0] = -fract / 5.0;
        a[1] = a[0];
        a[2] = 1.0 + fract * 4.0 / 5.0;
        a[3] = a[0];
        a[4] = a[0];
    }

    if (direction == L_HORIZ) {
        for (i = 0; i < h; i++) {
            lines = datas + i * wpls;
            lined = datad + i * wpld;
            if (halfwidth == 1) {
                for (j = 1; j < w - 1; j++) {
                    val = a[0] * GET_DATA_BYTE(lines, j - 1) +
                          a[1] * GET_DATA_BYTE(lines, j) +
                          a[2] * GET_DATA_BYTE(lines, j + 1);
                    ival = (l_int32)val;
                    ival = L_MAX(0, ival);
                    ival = L_MIN(255, ival);
                    SET_DATA_BYTE(lined, j, ival);
                }
            }
            else {  /* halfwidth == 2 */
                for (j = 2; j < w - 2; j++) {
                    val = a[0] * GET_DATA_BYTE(lines, j - 2) +
                          a[1] * GET_DATA_BYTE(lines, j - 1) +
                          a[2] * GET_DATA_BYTE(lines, j) +
                          a[3] * GET_DATA_BYTE(lines, j + 1) +
                          a[4] * GET_DATA_BYTE(lines, j + 2);
                    ival = (l_int32)val;
                    ival = L_MAX(0, ival);
                    ival = L_MIN(255, ival);
                    SET_DATA_BYTE(lined, j, ival);
                }
            }
        }
    }
    else {  /* direction == L_VERT */
        if (halfwidth == 1) {
            for (i = 1; i < h - 1; i++) {
                lines0 = datas + (i - 1) * wpls;
                lines1 = datas + i * wpls;
                lines2 = datas + (i + 1) * wpls;
                lined = datad + i * wpld;
                for (j = 0; j < w; j++) {
                    val = a[0] * GET_DATA_BYTE(lines0, j) +
                          a[1] * GET_DATA_BYTE(lines1, j) +
                          a[2] * GET_DATA_BYTE(lines2, j);
                    ival = (l_int32)val;
                    ival = L_MAX(0, ival);
                    ival = L_MIN(255, ival);
                    SET_DATA_BYTE(lined, j, ival);
                }
            }
        }
        else {  /* halfwidth == 2 */
            for (i = 2; i < h - 2; i++) {
                lines0 = datas + (i - 2) * wpls;
                lines1 = datas + (i - 1) * wpls;
                lines2 = datas + i * wpls;
                lines3 = datas + (i + 1) * wpls;
                lines4 = datas + (i + 2) * wpls;
                lined = datad + i * wpld;
                for (j = 0; j < w; j++) {
                    val = a[0] * GET_DATA_BYTE(lines0, j) +
                          a[1] * GET_DATA_BYTE(lines1, j) +
                          a[2] * GET_DATA_BYTE(lines2, j) +
                          a[3] * GET_DATA_BYTE(lines3, j) +
                          a[4] * GET_DATA_BYTE(lines4, j);
                    ival = (l_int32)val;
                    ival = L_MAX(0, ival);
                    ival = L_MIN(255, ival);
                    SET_DATA_BYTE(lined, j, ival);
                }
            }
        }
    }

    return pixd;
}


/*!
 *  pixUnsharpMaskingGray2D()
 *
 *      Input:  pixs (8 bpp; no colormap)
 *              halfwidth  ("half-width" of smoothing filter: 1 or 2)
 *              fract  (fraction of high frequency added to image)
 *      Return: pixd, or null on error
 *
 *  Notes:
 *      (1) For halfwidth == 1, we implement the full sharpening filter
 *          directly.  For halfwidth == 2, we implement the the lowpass
 *          filter separably and then compute the sharpening result locally.
 *      (2) Returns a clone if no sharpening is requested.
 */
PIX *
pixUnsharpMaskingGray2D(PIX       *pixs,
                        l_int32    halfwidth,
                        l_float32  fract)
{
l_int32     w, h, d, wpls, wpld, wplf, i, j, ival, sval;
l_uint32   *datas, *datad, *lines, *lines0, *lines1, *lines2, *lined;
l_float32   val, a[9];
l_float32  *dataf, *linef, *linef0, *linef1, *linef2, *linef3, *linef4;
PIX        *pixd;
FPIX       *fpix;

    PROCNAME("pixUnsharpMaskingGray2D");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 8 || pixGetColormap(pixs) != NULL)
        return (PIX *)ERROR_PTR("pixs not 8 bpp or has cmap", procName, NULL);
    if (fract <= 0.0 || halfwidth <= 0) {
        L_WARNING("no sharpening requested; clone returned", procName);
        return pixClone(pixs);
    }
    if (halfwidth != 1 && halfwidth != 2)
        return (PIX *)ERROR_PTR("halfwidth must be 1 or 2", procName, NULL);

    pixd = pixCopyBorder(NULL, pixs, halfwidth, halfwidth,
                         halfwidth, halfwidth);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);
    datas = pixGetData(pixs);
    wpls = pixGetWpl(pixs);

    if (halfwidth == 1) {
        for (i = 0; i < 9; i++)
            a[i] = -fract / 9.0;
        a[4] = 1.0 + fract * 8.0 / 9.0;
        for (i = 1; i < h - 1; i++) {
            lines0 = datas + (i - 1) * wpls;
            lines1 = datas + i * wpls;
            lines2 = datas + (i + 1) * wpls;
            lined = datad + i * wpld;
            for (j = 1; j < w - 1; j++) {
                val = a[0] * GET_DATA_BYTE(lines0, j - 1) +
                      a[1] * GET_DATA_BYTE(lines0, j) +
                      a[2] * GET_DATA_BYTE(lines0, j + 1) +
                      a[3] * GET_DATA_BYTE(lines1, j - 1) +
                      a[4] * GET_DATA_BYTE(lines1, j) +
                      a[5] * GET_DATA_BYTE(lines1, j + 1) +
                      a[6] * GET_DATA_BYTE(lines2, j - 1) +
                      a[7] * GET_DATA_BYTE(lines2, j) +
                      a[8] * GET_DATA_BYTE(lines2, j + 1);
                ival = (l_int32)(val + 0.5);
                ival = L_MAX(0, ival);
                ival = L_MIN(255, ival);
                SET_DATA_BYTE(lined, j, ival);
            }
        }

        return pixd;
    }

        /* For halfwidth == 2, do the low pass separably.  Store
         * the result of horizontal smoothing in an intermediate fpix. */
    fpix = fpixCreate(w, h);
    dataf = fpixGetData(fpix);
    wplf = fpixGetWpl(fpix);
    for (i = 2; i < h - 2; i++) {
        lines = datas + i * wpls;
        linef = dataf + i * wplf;
        for (j = 2; j < w - 2; j++) {
            val = GET_DATA_BYTE(lines, j - 2) +
                  GET_DATA_BYTE(lines, j - 1) +
                  GET_DATA_BYTE(lines, j) +
                  GET_DATA_BYTE(lines, j + 1) +
                  GET_DATA_BYTE(lines, j + 2);
            linef[j] = val;
        }
    }

        /* Do vertical smoothing to finish the low-pass filter.
         * At each pixel, if L is the lowpass value, I is the
         * src pixel value and f is the fraction of highpass to
         * be added to I, then the highpass filter value is
         *     H = I - L
         * and the new sharpened value is
         *     N = I + f * H.
         */
    for (i = 2; i < h - 2; i++) {
        linef0 = dataf + (i - 2) * wplf;
        linef1 = dataf + (i - 1) * wplf;
        linef2 = dataf + i * wplf;
        linef3 = dataf + (i + 1) * wplf;
        linef4 = dataf + (i + 2) * wplf;
        lined = datad + i * wpld;
        lines = datas + i * wpls;
        for (j = 2; j < w - 2; j++) {
            val = 0.04 * (linef0[j] + linef1[j] + linef2[j] +
                          linef3[j] + linef4[j]);  /* L: lowpass filter value */
            sval = GET_DATA_BYTE(lines, j);   /* I: source pixel */
            ival = (l_int32)(sval + fract * (sval - val) + 0.5);
            ival = L_MAX(0, ival);
            ival = L_MIN(255, ival);
            SET_DATA_BYTE(lined, j, ival);
        }
    }

    fpixDestroy(&fpix);
    return pixd;
}



/*-----------------------------------------------------------------------*
 *                    Hue and saturation modification                    *
 *-----------------------------------------------------------------------*/
/*!
 *  pixModifyHue()
 *
 *      Input:  pixd (<optional> can be null or equal to pixs)
 *              pixs (32 bpp rgb)
 *              fract (between -1.0 and 1.0)
 *      Return: pixd, or null on error
 *
 *  Notes:
 *      (1) pixd must either be null or equal to pixs.
 *          For in-place operation, set pixd == pixs:
 *             pixEqualizeTRC(pixs, pixs, ...);
 *          To get a new image, set pixd == null:
 *             pixd = pixEqualizeTRC(NULL, pixs, ...);
 *      (1) Use fract > 0.0 to increase hue value; < 0.0 to decrease it.
 *          1.0 (or -1.0) represents a 360 degree rotation; i.e., no change.
 *      (2) If no modification is requested (fract = -1.0 or 0 or 1.0),
 *          return a copy unless in-place, in which case this is a no-op.
 */
PIX  *
pixModifyHue(PIX       *pixd,
             PIX       *pixs,
             l_float32  fract)
{
l_int32    w, h, d, i, j, wpl, delhue;
l_int32    rval, gval, bval, hval, sval, vval;
l_uint32  *data, *line;

    PROCNAME("pixModifyHue");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (pixGetColormap(pixs) != NULL)
        return (PIX *)ERROR_PTR("pixs colormapped", procName, NULL);
    if (pixd && (pixd != pixs))
        return (PIX *)ERROR_PTR("pixd not null or pixs", procName, pixd);
    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 32) 
        return (PIX *)ERROR_PTR("pixs not 32 bpp", procName, NULL);
    if (L_ABS(fract) > 1.0)
        return (PIX *)ERROR_PTR("fract not in [-1.0 ... 1.0]", procName, NULL);

    pixd = pixCopy(pixd, pixs);

    delhue = (l_int32)(240 * fract);
    if (delhue == 0 || delhue == 240 || delhue == -240) {
        L_WARNING("no change requested in hue", procName);
        return pixd;
    }
    if (delhue < 0)
        delhue += 240;

    data = pixGetData(pixd);
    wpl = pixGetWpl(pixd);
    for (i = 0; i < h; i++) {
        line = data + i * wpl;
        for (j = 0; j < w; j++) {
            extractRGBValues(line[j], &rval, &gval, &bval);
            convertRGBToHSV(rval, gval, bval, &hval, &sval, &vval);
            hval = (hval + delhue) % 240;
            convertHSVToRGB(hval, sval, vval, &rval, &gval, &bval);
            composeRGBPixel(rval, gval, bval, line + j);
        }
    }

    return pixd;
}


/*!
 *  pixModifySaturation()
 *
 *      Input:  pixd (<optional> can be null, existing or equal to pixs)
 *              pixs (32 bpp rgb)
 *              fract (between -1.0 and 1.0)
 *      Return: pixd, or null on error
 *
 *  Notes:
 *      (1) If fract > 0.0, it gives the fraction that the pixel
 *          saturation is moved from its initial value toward 255.
 *          If fract < 0.0, it gives the fraction that the pixel
 *          saturation is moved from its initial value toward 0.
 *          The limiting values for fract = -1.0 (1.0) thus set the
 *          saturation to 0 (255).
 *      (2) If fract = 0, no modification is requested; return a copy
 *          unless in-place, in which case this is a no-op.
 */
PIX  *
pixModifySaturation(PIX       *pixd,
                    PIX       *pixs,
                    l_float32  fract)
{
l_int32    w, h, d, i, j, wpl;
l_int32    rval, gval, bval, hval, sval, vval;
l_uint32  *data, *line;

    PROCNAME("pixModifySaturation");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 32) 
        return (PIX *)ERROR_PTR("pixs not 32 bpp", procName, NULL);
    if (L_ABS(fract) > 1.0)
        return (PIX *)ERROR_PTR("fract not in [-1.0 ... 1.0]", procName, NULL);

    pixd = pixCopy(pixd, pixs);
    if (fract == 0.0) {
        L_WARNING("no change requested in saturation", procName);
        return pixd;
    }

    data = pixGetData(pixd);
    wpl = pixGetWpl(pixd);
    for (i = 0; i < h; i++) {
        line = data + i * wpl;
        for (j = 0; j < w; j++) {
            extractRGBValues(line[j], &rval, &gval, &bval);
            convertRGBToHSV(rval, gval, bval, &hval, &sval, &vval);
            if (fract < 0.0)
                sval = (l_int32)(sval * (1.0 + fract));
            else
                sval = (l_int32)(sval + fract * (255 - sval));
            convertHSVToRGB(hval, sval, vval, &rval, &gval, &bval);
            composeRGBPixel(rval, gval, bval, line + j);
        }
    }

    return pixd;
}


/*!
 *  pixMeasureSaturation()
 *
 *      Input:  pixs (32 bpp rgb)
 *              factor (subsampling factor; integer >= 1)
 *              &sat (<return> average saturation)
 *      Return: pixd, or null on error
 */
l_int32
pixMeasureSaturation(PIX        *pixs,
                     l_int32     factor,
                     l_float32  *psat)
{
l_int32    w, h, d, i, j, wpl, sum, count;
l_int32    rval, gval, bval, hval, sval, vval;
l_uint32  *data, *line;

    PROCNAME("pixMeasureSaturation");

    if (!psat)
        return ERROR_INT("pixs not defined", procName, 1);
    *psat = 0.0;
    if (!pixs)
        return ERROR_INT("pixs not defined", procName, 1);
    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 32) 
        return ERROR_INT("pixs not 32 bpp", procName, 1);
    if (factor < 1)
        return ERROR_INT("subsampling factor < 1", procName, 1);

    data = pixGetData(pixs);
    wpl = pixGetWpl(pixs);
    for (i = 0, sum = 0, count = 0; i < h; i += factor) {
        line = data + i * wpl;
        for (j = 0; j < w; j += factor) {
            extractRGBValues(line[j], &rval, &gval, &bval);
            convertRGBToHSV(rval, gval, bval, &hval, &sval, &vval);
            sum += sval;
            count++;
        }
    }

    if (count > 0)
        *psat = (l_float32)sum / (l_float32)count;
    return 0;
}


/*-----------------------------------------------------------------------*
 *            General multiplicative constant color transform            *
 *-----------------------------------------------------------------------*/
/*
 *  pixMultConstantColor()
 *
 *      Input:  pixs (colormapped or rgb)
 *              rfact, gfact, bfact (multiplicative factors on each component)
 *      Return: pixd (colormapped or rgb, with colors scaled), or null on error
 *
 *  Notes:
 *      (1) rfact, gfact and bfact can only have non-negative values.
 *          They can be greater than 1.0.  All transformed component
 *          values are clipped to the interval [0, 255].
 *      (2) For multiplication with a general 3x3 matrix of constants,
 *          use pixMultMatrixColor().
 */
PIX *
pixMultConstantColor(PIX       *pixs,
                     l_float32  rfact,
                     l_float32  gfact,
                     l_float32  bfact)
{
l_int32    i, j, w, h, d, wpls, wpld;
l_int32    ncolors, rval, gval, bval, nrval, ngval, nbval;
l_uint32   nval;
l_uint32  *datas, *datad, *lines, *lined;
PIX       *pixd;
PIXCMAP   *cmap;

    PROCNAME("pixMultConstantColor");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    pixGetDimensions(pixs, &w, &h, &d);
    cmap = pixGetColormap(pixs);
    if (!cmap && d != 32)
        return (PIX *)ERROR_PTR("pixs not cmapped or 32 bpp", procName, NULL);
    rfact = L_MAX(0.0, rfact);
    gfact = L_MAX(0.0, gfact);
    bfact = L_MAX(0.0, bfact);

    if (cmap) {
        if ((pixd = pixCopy(NULL, pixs)) == NULL)
            return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
        cmap = pixGetColormap(pixd);
        ncolors = pixcmapGetCount(cmap);
        for (i = 0; i < ncolors; i++) {
            pixcmapGetColor(cmap, i, &rval, &gval, &bval);
            nrval = (l_int32)(rfact * rval);
            ngval = (l_int32)(gfact * gval);
            nbval = (l_int32)(bfact * bval);
            nrval = L_MIN(255, nrval);
            ngval = L_MIN(255, ngval);
            nbval = L_MIN(255, nbval);
            pixcmapResetColor(cmap, i, nrval, ngval, nbval);
        }
        return pixd;
    }

    if ((pixd = pixCreateTemplateNoInit(pixs)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    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++) {
            extractRGBValues(lines[j], &rval, &gval, &bval);
            nrval = (l_int32)(rfact * rval);
            ngval = (l_int32)(gfact * gval);
            nbval = (l_int32)(bfact * bval);
            nrval = L_MIN(255, nrval);
            ngval = L_MIN(255, ngval);
            nbval = L_MIN(255, nbval);
            composeRGBPixel(nrval, ngval, nbval, &nval);
            *(lined + j) = nval;
        }
    }

    return pixd;
}


/*
 *  pixMultMatrixColor()
 *
 *      Input:  pixs (colormapped or rgb)
 *              kernel (3x3 matrix of floats)
 *      Return: pixd (colormapped or rgb), or null on error
 *
 *  Notes:
 *      (1) The kernel is a data structure used mostly for floating point
 *          convolution.  Here it is a 3x3 matrix of floats that are used
 *          to transform the pixel values by matrix multiplication:
 *            nrval = a[0,0] * rval + a[0,1] * gval + a[0,2] * bval
 *            ngval = a[1,0] * rval + a[1,1] * gval + a[1,2] * bval
 *            nbval = a[2,0] * rval + a[2,1] * gval + a[2,2] * bval
 *      (2) The matrix can be generated in several ways.
 *          See kernel.c for details.  Here are two of them:
 *            (a) kel = kernelCreate(3, 3);
 *                kernelSetElement(kel, 0, 0, val00);
 *                kernelSetElement(kel, 0, 1, val01);
 *                ...
 *            (b) from a static string; e.g.,:
 *                const char *kdata = " 0.6  0.3 -0.2 "
 *                                    " 0.1  1.2  0.4 "
 *                                    " -0.4 0.2  0.9 ";
 *                kel = kernelCreateFromString(3, 3, 0, 0, kdata);
 *      (3) For the special case where the matrix is diagonal, it is easier
 *          to use pixMultConstantColor().
 *      (4) Matrix entries can have positive and negative values, and can
 *          be larger than 1.0.  All transformed component values
 *          are clipped to [0, 255].
 */
PIX *
pixMultMatrixColor(PIX       *pixs,
                   L_KERNEL  *kel)
{
l_int32    i, j, index, kw, kh, w, h, d, wpls, wpld;
l_int32    ncolors, rval, gval, bval, nrval, ngval, nbval;
l_uint32   nval;
l_uint32  *datas, *datad, *lines, *lined;
l_float32  v[9];  /* use linear array for convenience */
PIX       *pixd;
PIXCMAP   *cmap;

    PROCNAME("pixMultMatrixColor");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (!kel)
        return (PIX *)ERROR_PTR("kel not defined", procName, NULL);
    kernelGetParameters(kel, &kw, &kh, NULL, NULL);
    if (kw != 3 || kh != 3)
        return (PIX *)ERROR_PTR("matrix not 3x3", procName, NULL);
    pixGetDimensions(pixs, &w, &h, &d);
    cmap = pixGetColormap(pixs);
    if (!cmap && d != 32)
        return (PIX *)ERROR_PTR("pixs not cmapped or 32 bpp", procName, NULL);

    for (i = 0, index = 0; i < 3; i++)
        for (j = 0; j < 3; j++, index++)
            kernelGetElement(kel, i, j, v + index);

    if (cmap) {
        if ((pixd = pixCopy(NULL, pixs)) == NULL)
            return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
        cmap = pixGetColormap(pixd);
        ncolors = pixcmapGetCount(cmap);
        for (i = 0; i < ncolors; i++) {
            pixcmapGetColor(cmap, i, &rval, &gval, &bval);
            nrval = (l_int32)(v[0] * rval + v[1] * gval + v[2] * bval);
            ngval = (l_int32)(v[3] * rval + v[4] * gval + v[5] * bval);
            nbval = (l_int32)(v[6] * rval + v[7] * gval + v[8] * bval);
            nrval = L_MAX(0, L_MIN(255, nrval));
            ngval = L_MAX(0, L_MIN(255, ngval));
            nbval = L_MAX(0, L_MIN(255, nbval));
            pixcmapResetColor(cmap, i, nrval, ngval, nbval);
        }
        return pixd;
    }

    if ((pixd = pixCreateTemplateNoInit(pixs)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    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++) {
            extractRGBValues(lines[j], &rval, &gval, &bval);
            nrval = (l_int32)(v[0] * rval + v[1] * gval + v[2] * bval);
            ngval = (l_int32)(v[3] * rval + v[4] * gval + v[5] * bval);
            nbval = (l_int32)(v[6] * rval + v[7] * gval + v[8] * bval);
            nrval = L_MAX(0, L_MIN(255, nrval));
            ngval = L_MAX(0, L_MIN(255, ngval));
            nbval = L_MAX(0, L_MIN(255, nbval));
            composeRGBPixel(nrval, ngval, nbval, &nval);
            *(lined + j) = nval;
        }
    }

    return pixd;
}


/*-------------------------------------------------------------*
 *                    Half-edge by bandpass                    *
 *-------------------------------------------------------------*/
/*!
 *  pixHalfEdgeByBandpass()
 *
 *      Input:  pixs (8 bpp gray or 32 bpp rgb)
 *              sm1h, sm1v ("half-widths" of smoothing filter sm1)
 *              sm2h, sm2v ("half-widths" of smoothing filter sm2)
 *                      (require sm2 != sm1)
 *      Return: pixd, or null on error
 *
 *  Notes:
 *      (1) We use symmetric smoothing filters of odd dimension,
 *          typically use 3, 5, 7, etc.  The smoothing parameters
 *          for these are 1, 2, 3, etc.  The filter size is related
 *          to the smoothing parameter by
 *               size = 2 * smoothing + 1
 *      (2) Because we take the difference of two lowpass filters,
 *          this is actually a bandpass filter.
 *      (3) We allow both filters to be anisotropic.
 *      (4) Consider either the h or v component of the 2 filters.
 *          Depending on whether sm1 > sm2 or sm2 > sm1, we get
 *          different halves of the smoothed gradients (or "edges").
 *          This difference of smoothed signals looks more like
 *          a second derivative of a transition, which we rectify
 *          by not allowing the signal to go below zero.  If sm1 < sm2,
 *          the sm2 transition is broader, so the difference between
 *          sm1 and sm2 signals is positive on the upper half of
 *          the transition.  Likewise, if sm1 > sm2, the sm1 - sm2
 *          signal difference is positive on the lower half of
 *          the transition.
 */
PIX *
pixHalfEdgeByBandpass(PIX     *pixs,
                      l_int32  sm1h,
                      l_int32  sm1v,
                      l_int32  sm2h,
                      l_int32  sm2v)
{
l_int32  d;
PIX     *pixg, *pixacc, *pixc1, *pixc2;

    PROCNAME("pixHalfEdgeByBandpass");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (sm1h == sm2h && sm1v == sm2v)
        return (PIX *)ERROR_PTR("sm2 = sm1", procName, NULL);
    d = pixGetDepth(pixs);
    if (d != 8 && d != 32)
        return (PIX *)ERROR_PTR("pixs not 8 or 32 bpp", procName, NULL);
    if (d == 32)
        pixg = pixConvertRGBToLuminance(pixs);
    else   /* d == 8 */
        pixg = pixClone(pixs);

        /* Make a convolution accumulator and use it twice */
    if ((pixacc = pixBlockconvAccum(pixg)) == NULL)
        return (PIX *)ERROR_PTR("pixacc not made", procName, NULL);
    if ((pixc1 = pixBlockconvGray(pixg, pixacc, sm1h, sm1v)) == NULL)
        return (PIX *)ERROR_PTR("pixc1 not made", procName, NULL);
    if ((pixc2 = pixBlockconvGray(pixg, pixacc, sm2h, sm2v)) == NULL)
        return (PIX *)ERROR_PTR("pixc2 not made", procName, NULL);
    pixDestroy(&pixacc);

        /* Compute the half-edge using pixc1 - pixc2.  */
    pixSubtractGray(pixc1, pixc1, pixc2);

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