diff mupdf-source/thirdparty/leptonica/src/enhance.c @ 2:b50eed0cc0ef upstream

ADD: MuPDF v1.26.7: the MuPDF source as downloaded by a default build of PyMuPDF 1.26.4. The directory name has changed: no version number in the expanded directory now.
author Franz Glasner <fzglas.hg@dom66.de>
date Mon, 15 Sep 2025 11:43:07 +0200
parents
children
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mupdf-source/thirdparty/leptonica/src/enhance.c	Mon Sep 15 11:43:07 2025 +0200
@@ -0,0 +1,2318 @@
+/*====================================================================*
+ -  Copyright (C) 2001 Leptonica.  All rights reserved.
+ -
+ -  Redistribution and use in source and binary forms, with or without
+ -  modification, are permitted provided that the following conditions
+ -  are met:
+ -  1. Redistributions of source code must retain the above copyright
+ -     notice, this list of conditions and the following disclaimer.
+ -  2. Redistributions in binary form must reproduce the above
+ -     copyright notice, this list of conditions and the following
+ -     disclaimer in the documentation and/or other materials
+ -     provided with the distribution.
+ -
+ -  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ -  ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ -  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ -  A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL ANY
+ -  CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ -  EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ -  PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ -  PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
+ -  OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ -  NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ -  SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *====================================================================*/
+
+/*!
+ * \file enhance.c
+ * <pre>
+ *
+ *      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
+ *           l_int32  pixTRCMap()
+ *           l_int32  pixTRCMapGeneral()
+ *
+ *      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()
+ *           PIX     *pixModifyBrightness()
+ *
+ *      Color shifting
+ *           PIX     *pixMosaicColorShiftRGB()
+ *           PIX     *pixColorShiftRGB()
+ *
+ *      Darken gray (unsaturated) pixels
+ *           PIX     *pixDarkenGray()
+ *
+ *      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.
+ * </pre>
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config_auto.h>
+#endif  /* HAVE_CONFIG_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  EnhanceScaleFactor = 5.0;
+
+/*-------------------------------------------------------------*
+ *         Gamma TRC (tone reproduction curve) mapping         *
+ *-------------------------------------------------------------*/
+/*!
+ * \brief   pixGammaTRC()
+ *
+ * \param[in]    pixd     [optional] null or equal to pixs
+ * \param[in]    pixs     8 or 32 bpp; or 2, 4 or 8 bpp with colormap
+ * \param[in]    gamma    gamma correction; must be > 0.0
+ * \param[in]    minval   input value that gives 0 for output; can be < 0
+ * \param[in]    maxval   input value that gives 255 for output; can be > 255
+ * \return  pixd always
+ *
+ * <pre>
+ * 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.
+ *      (12) Use pixTRCMapGeneral() if applying different mappings
+ *           to each channel in an RGB image.
+ * </pre>
+ */
+PIX *
+pixGammaTRC(PIX       *pixd,
+            PIX       *pixs,
+            l_float32  gamma,
+            l_int32    minval,
+            l_int32    maxval)
+{
+l_int32   d;
+NUMA     *nag;
+PIXCMAP  *cmap;
+
+    if (!pixs)
+        return (PIX *)ERROR_PTR("pixs not defined", __func__, pixd);
+    if (pixd && (pixd != pixs))
+        return (PIX *)ERROR_PTR("pixd not null or pixs", __func__, pixd);
+    if (gamma <= 0.0) {
+        L_WARNING("gamma must be > 0.0; setting to 1.0\n", __func__);
+        gamma = 1.0;
+    }
+    if (minval >= maxval)
+        return (PIX *)ERROR_PTR("minval not < maxval", __func__, pixd);
+    cmap = pixGetColormap(pixs);
+    d = pixGetDepth(pixs);
+    if (!cmap && d != 8 && d != 32)
+        return (PIX *)ERROR_PTR("depth not 8 or 32 bpp", __func__, 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", __func__, pixd);
+    pixTRCMap(pixd, NULL, nag);
+    numaDestroy(&nag);
+
+    return pixd;
+}
+
+
+/*!
+ * \brief   pixGammaTRCMasked()
+ *
+ * \param[in]    pixd      [optional] null or equal to pixs
+ * \param[in]    pixs      8 or 32 bpp; not colormapped
+ * \param[in]    pixm      [optional] null or 1 bpp
+ * \param[in]    gamma     gamma correction; must be > 0.0
+ * \param[in]    minval    input value that gives 0 for output; can be < 0
+ * \param[in]    maxval    input value that gives 255 for output; can be > 255
+ * \return  pixd always
+ *
+ * <pre>
+ * 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.
+ * </pre>
+ */
+PIX *
+pixGammaTRCMasked(PIX       *pixd,
+                  PIX       *pixs,
+                  PIX       *pixm,
+                  l_float32  gamma,
+                  l_int32    minval,
+                  l_int32    maxval)
+{
+l_int32  d;
+NUMA    *nag;
+
+    if (!pixm)
+        return pixGammaTRC(pixd, pixs, gamma, minval, maxval);
+
+    if (!pixs)
+        return (PIX *)ERROR_PTR("pixs not defined", __func__, pixd);
+    if (pixGetColormap(pixs))
+        return (PIX *)ERROR_PTR("invalid: pixs has a colormap", __func__, pixd);
+    if (pixd && (pixd != pixs))
+        return (PIX *)ERROR_PTR("pixd not null or pixs", __func__, pixd);
+    d = pixGetDepth(pixs);
+    if (d != 8 && d != 32)
+        return (PIX *)ERROR_PTR("depth not 8 or 32 bpp", __func__, pixd);
+    if (minval >= maxval)
+        return (PIX *)ERROR_PTR("minval not < maxval", __func__, pixd);
+    if (gamma <= 0.0) {
+        L_WARNING("gamma must be > 0.0; setting to 1.0\n", __func__);
+        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", __func__, pixd);
+    pixTRCMap(pixd, pixm, nag);
+    numaDestroy(&nag);
+
+    return pixd;
+}
+
+
+/*!
+ * \brief   pixGammaTRCWithAlpha()
+ *
+ * \param[in]    pixd     [optional] null or equal to pixs
+ * \param[in]    pixs     32 bpp
+ * \param[in]    gamma    gamma correction; must be > 0.0
+ * \param[in]    minval   input value that gives 0 for output; can be < 0
+ * \param[in]    maxval   input value that gives 255 for output; can be > 255
+ * \return  pixd always
+ *
+ * <pre>
+ * 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.
+ * </pre>
+ */
+PIX *
+pixGammaTRCWithAlpha(PIX       *pixd,
+                     PIX       *pixs,
+                     l_float32  gamma,
+                     l_int32    minval,
+                     l_int32    maxval)
+{
+NUMA  *nag;
+PIX   *pixalpha;
+
+    if (!pixs || pixGetDepth(pixs) != 32)
+        return (PIX *)ERROR_PTR("pixs undefined or not 32 bpp", __func__, pixd);
+    if (pixd && (pixd != pixs))
+        return (PIX *)ERROR_PTR("pixd not null or pixs", __func__, pixd);
+    if (gamma <= 0.0) {
+        L_WARNING("gamma must be > 0.0; setting to 1.0\n", __func__);
+        gamma = 1.0;
+    }
+    if (minval >= maxval)
+        return (PIX *)ERROR_PTR("minval not < maxval", __func__, pixd);
+
+    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);
+
+    pixalpha = pixGetRGBComponent(pixs, L_ALPHA_CHANNEL);  /* save */
+    if ((nag = numaGammaTRC(gamma, minval, maxval)) == NULL)
+        return (PIX *)ERROR_PTR("nag not made", __func__, pixd);
+    pixTRCMap(pixd, NULL, nag);
+    pixSetRGBComponent(pixd, pixalpha, L_ALPHA_CHANNEL);  /* restore */
+    pixSetSpp(pixd, 4);
+
+    numaDestroy(&nag);
+    pixDestroy(&pixalpha);
+    return pixd;
+}
+
+
+/*!
+ * \brief   numaGammaTRC()
+ *
+ * \param[in]    gamma   gamma factor; must be > 0.0
+ * \param[in]    minval  input value that gives 0 for output
+ * \param[in]    maxval  input value that gives 255 for output
+ * \return  na, or NULL on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) The map is returned as a numa; values are clipped to [0, 255].
+ *      (2) For a linear mapping, set gamma = 1.0.
+ *      (3) To force all intensities into a range within fraction delta
+ *          of white, use: minval = -256 * (1 - delta) / delta
+ *                         maxval = 255
+ *      (4) To force all intensities into a range within fraction delta
+ *          of black, use: minval = 0
+ *                         maxval = 256 * (1 - delta) / delta
+ * </pre>
+ */
+NUMA *
+numaGammaTRC(l_float32  gamma,
+             l_int32    minval,
+             l_int32    maxval)
+{
+l_int32    i, val;
+l_float32  x, invgamma;
+NUMA      *na;
+
+    if (minval >= maxval)
+        return (NUMA *)ERROR_PTR("minval not < maxval", __func__, NULL);
+    if (gamma <= 0.0) {
+        L_WARNING("gamma must be > 0.0; setting to 1.0\n", __func__);
+        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                   *
+ *-------------------------------------------------------------*/
+/*!
+ * \brief   pixContrastTRC()
+ *
+ * \param[in]    pixd     [optional] null or equal to pixs
+ * \param[in]    pixs     8 or 32 bpp; or 2, 4 or 8 bpp with colormap
+ * \param[in]    factor   0.0 is no enhancement
+ * \return  pixd always
+ *
+ * <pre>
+ * 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.
+ * </pre>
+ */
+PIX *
+pixContrastTRC(PIX       *pixd,
+               PIX       *pixs,
+               l_float32  factor)
+{
+l_int32   d;
+NUMA     *nac;
+PIXCMAP  *cmap;
+
+    if (!pixs)
+        return (PIX *)ERROR_PTR("pixs not defined", __func__, pixd);
+    if (pixd && (pixd != pixs))
+        return (PIX *)ERROR_PTR("pixd not null or pixs", __func__, pixd);
+    if (factor < 0.0) {
+        L_WARNING("factor must be >= 0.0; using 0.0\n", __func__);
+        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", __func__, 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", __func__, pixd);
+    pixTRCMap(pixd, NULL, nac);
+    numaDestroy(&nac);
+
+    return pixd;
+}
+
+
+/*!
+ * \brief   pixContrastTRCMasked()
+ *
+ * \param[in]    pixd     [optional] null or equal to pixs
+ * \param[in]    pixs     8 or 32 bpp; or 2, 4 or 8 bpp with colormap
+ * \param[in]    pixm     [optional] null or 1 bpp
+ * \param[in]    factor   0.0 is no enhancement
+ * \return  pixd always
+ *
+ * <pre>
+ * 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.
+ * </pre>
+ */
+PIX *
+pixContrastTRCMasked(PIX       *pixd,
+                     PIX       *pixs,
+                     PIX       *pixm,
+                     l_float32  factor)
+{
+l_int32  d;
+NUMA    *nac;
+
+    if (!pixm)
+        return pixContrastTRC(pixd, pixs, factor);
+
+    if (!pixs)
+        return (PIX *)ERROR_PTR("pixs not defined", __func__, pixd);
+    if (pixGetColormap(pixs))
+        return (PIX *)ERROR_PTR("invalid: pixs has a colormap", __func__, pixd);
+    if (pixd && (pixd != pixs))
+        return (PIX *)ERROR_PTR("pixd not null or pixs", __func__, pixd);
+    d = pixGetDepth(pixs);
+    if (d != 8 && d != 32)
+        return (PIX *)ERROR_PTR("depth not 8 or 32 bpp", __func__, pixd);
+
+    if (factor < 0.0) {
+        L_WARNING("factor must be >= 0.0; using 0.0\n", __func__);
+        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", __func__, pixd);
+    pixTRCMap(pixd, pixm, nac);
+    numaDestroy(&nac);
+
+    return pixd;
+}
+
+
+/*!
+ * \brief   numaContrastTRC()
+ *
+ * \param[in]    factor   generally between 0.0 [no enhancement]
+ *                        and 1.0, but can be larger than 1.0
+ * \return  na, or NULL on error
+ *
+ * <pre>
+ * 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).
+ * </pre>
+ */
+NUMA *
+numaContrastTRC(l_float32  factor)
+{
+l_int32    i, val;
+l_float64  x, ymax, ymin, dely, scale;
+NUMA      *na;
+
+    if (factor < 0.0) {
+        L_WARNING("factor must be >= 0.0; using 0.0; no enhancement\n",
+                  __func__);
+        factor = 0.0;
+    }
+    if (factor == 0.0)
+        return numaMakeSequence(0, 1, 256);  /* linear map */
+
+    scale = EnhanceScaleFactor;
+    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                  *
+ *-------------------------------------------------------------*/
+/*!
+ * \brief   pixEqualizeTRC()
+ *
+ * \param[in]    pixd     [optional] null or equal to pixs
+ * \param[in]    pixs     8 bpp gray, 32 bpp rgb, or colormapped
+ * \param[in]    fract    fraction of equalization movement of pixel values
+ * \param[in]    factor   subsampling factor; integer >= 1
+ * \return  pixd, or NULL on error
+ *
+ * <pre>
+ * 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).
+ * </pre>
+ */
+PIX *
+pixEqualizeTRC(PIX       *pixd,
+               PIX       *pixs,
+               l_float32  fract,
+               l_int32    factor)
+{
+l_int32   d;
+NUMA     *na;
+PIX      *pixt, *pix8;
+PIXCMAP  *cmap;
+
+    if (!pixs)
+        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
+    if (pixd && (pixd != pixs))
+        return (PIX *)ERROR_PTR("pixd not null or pixs", __func__, pixd);
+    cmap = pixGetColormap(pixs);
+    d = pixGetDepth(pixs);
+    if (d != 8 && d != 32 && !cmap)
+        return (PIX *)ERROR_PTR("pixs not 8/32 bpp or cmapped", __func__, NULL);
+    if (fract < 0.0 || fract > 1.0)
+        return (PIX *)ERROR_PTR("fract not in [0.0 ... 1.0]", __func__, NULL);
+    if (factor < 1)
+        return (PIX *)ERROR_PTR("sampling factor < 1", __func__, 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;
+}
+
+
+/*!
+ * \brief   numaEqualizeTRC()
+ *
+ * \param[in]    pix     8 bpp, no colormap
+ * \param[in]    fract   fraction of equalization movement of pixel values
+ * \param[in]    factor  subsampling factor; integer >= 1
+ * \return  nad, or NULL on error
+ *
+ * <pre>
+ * 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).
+ * </pre>
+ */
+NUMA *
+numaEqualizeTRC(PIX       *pix,
+                l_float32  fract,
+                l_int32    factor)
+{
+l_int32    iin, iout, itarg;
+l_float32  val, sum;
+NUMA      *nah, *nasum, *nad;
+
+    if (!pix)
+        return (NUMA *)ERROR_PTR("pix not defined", __func__, NULL);
+    if (pixGetDepth(pix) != 8)
+        return (NUMA *)ERROR_PTR("pix not 8 bpp", __func__, NULL);
+    if (fract < 0.0 || fract > 1.0)
+        return (NUMA *)ERROR_PTR("fract not in [0.0 ... 1.0]", __func__, NULL);
+    if (factor < 1)
+        return (NUMA *)ERROR_PTR("sampling factor < 1", __func__, NULL);
+
+    if (fract == 0.0)
+        L_WARNING("fract = 0.0; no equalization requested\n", __func__);
+
+    if ((nah = pixGetGrayHistogram(pix, factor)) == NULL)
+        return (NUMA *)ERROR_PTR("histogram not made", __func__, 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                   *
+ *-------------------------------------------------------------*/
+/*!
+ * \brief   pixTRCMap()
+ *
+ * \param[in]    pixs    8 grayscale or 32 bpp rgb; not colormapped
+ * \param[in]    pixm    [optional] 1 bpp mask
+ * \param[in]    na      mapping array
+ * \return  0 if OK, 1 on error
+ *
+ * <pre>
+ * 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.
+ * </pre>
+ */
+l_int32
+pixTRCMap(PIX   *pixs,
+          PIX   *pixm,
+          NUMA  *na)
+{
+l_int32    w, h, d, wm, hm, wpl, wplm, i, j, sval8, dval8;
+l_uint32   sval32, dval32;
+l_uint32  *data, *datam, *line, *linem, *tab;
+
+    if (!pixs)
+        return ERROR_INT("pixs not defined", __func__, 1);
+    if (pixGetColormap(pixs))
+        return ERROR_INT("pixs is colormapped", __func__, 1);
+    if (!na)
+        return ERROR_INT("na not defined", __func__, 1);
+    if (numaGetCount(na) != 256)
+        return ERROR_INT("na not of size 256", __func__, 1);
+    pixGetDimensions(pixs, &w, &h, &d);
+    if (d != 8 && d != 32)
+        return ERROR_INT("pixs not 8 or 32 bpp", __func__, 1);
+    if (pixm) {
+        if (pixGetDepth(pixm) != 1)
+            return ERROR_INT("pixm not 1 bpp", __func__, 1);
+    }
+
+    tab = (l_uint32 *)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;
+                }
+            }
+        }
+    }
+
+    LEPT_FREE(tab);
+    return 0;
+}
+
+
+/*!
+ * \brief   pixTRCMapGeneral()
+ *
+ * \param[in]    pixs             32 bpp rgb; not colormapped
+ * \param[in]    pixm             [optional] 1 bpp mask
+ * \param[in]    nar, nag, nab    mapping arrays
+ * \return  0 if OK, 1 on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) This operation is in-place on %pixs.
+ *      (2) Each of the r,g,b mapping arrays is of size 256. They map the
+ *          input value for that color component into values in the
+ *          range [0, 255].
+ *      (3) In the special case where the r, g and b mapping arrays are
+ *          all the same, call pixTRCMap() instead.
+ *      (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) The alpha channel is not saved.
+ * </pre>
+ */
+l_int32
+pixTRCMapGeneral(PIX   *pixs,
+                 PIX   *pixm,
+                 NUMA  *nar,
+                 NUMA  *nag,
+                 NUMA  *nab)
+{
+l_int32    w, h, wm, hm, wpl, wplm, i, j;
+l_uint32   sval32, dval32;
+l_uint32  *data, *datam, *line, *linem, *tabr, *tabg, *tabb;
+
+    if (!pixs || pixGetDepth(pixs) != 32)
+        return ERROR_INT("pixs not defined or not 32 bpp", __func__, 1);
+    if (pixm && pixGetDepth(pixm) != 1)
+        return ERROR_INT("pixm defined and not 1 bpp", __func__, 1);
+    if (!nar || !nag || !nab)
+        return ERROR_INT("na{r,g,b} not all defined", __func__, 1);
+    if (numaGetCount(nar) != 256 || numaGetCount(nag) != 256 ||
+        numaGetCount(nab) != 256)
+        return ERROR_INT("na{r,g,b} not all of size 256", __func__, 1);
+
+        /* Get the arrays for efficiency */
+    tabr = (l_uint32 *)numaGetIArray(nar);
+    tabg = (l_uint32 *)numaGetIArray(nag);
+    tabb = (l_uint32 *)numaGetIArray(nab);
+    pixGetDimensions(pixs, &w, &h, NULL);
+    wpl = pixGetWpl(pixs);
+    data = pixGetData(pixs);
+    if (!pixm) {
+        for (i = 0; i < h; i++) {
+            line = data + i * wpl;
+            for (j = 0; j < w; j++) {
+                sval32 = *(line + j);
+                dval32 =
+                    tabr[(sval32 >> L_RED_SHIFT) & 0xff] << L_RED_SHIFT |
+                    tabg[(sval32 >> L_GREEN_SHIFT) & 0xff] << L_GREEN_SHIFT |
+                    tabb[(sval32 >> L_BLUE_SHIFT) & 0xff] << L_BLUE_SHIFT;
+                *(line + j) = dval32;
+            }
+        }
+    } else {
+        datam = pixGetData(pixm);
+        wplm = pixGetWpl(pixm);
+        pixGetDimensions(pixm, &wm, &hm, NULL);
+        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 =
+                    tabr[(sval32 >> L_RED_SHIFT) & 0xff] << L_RED_SHIFT |
+                    tabg[(sval32 >> L_GREEN_SHIFT) & 0xff] << L_GREEN_SHIFT |
+                    tabb[(sval32 >> L_BLUE_SHIFT) & 0xff] << L_BLUE_SHIFT;
+                *(line + j) = dval32;
+            }
+        }
+    }
+
+    LEPT_FREE(tabr);
+    LEPT_FREE(tabg);
+    LEPT_FREE(tabb);
+    return 0;
+}
+
+
+
+/*-----------------------------------------------------------------------*
+ *                             Unsharp masking                           *
+ *-----------------------------------------------------------------------*/
+/*!
+ * \brief   pixUnsharpMasking()
+ *
+ * \param[in]    pixs       all depths except 1 bpp; with or without colormaps
+ * \param[in]    halfwidth  "half-width" of smoothing filter
+ * \param[in]    fract      fraction of edge added back into image
+ * \return  pixd, or NULL on error
+ *
+ * <pre>
+ * 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.
+ * </pre>
+ */
+PIX *
+pixUnsharpMasking(PIX       *pixs,
+                  l_int32    halfwidth,
+                  l_float32  fract)
+{
+l_int32  d;
+PIX     *pix1, *pixd, *pixr, *pixrs, *pixg, *pixgs, *pixb, *pixbs;
+
+    if (!pixs || (pixGetDepth(pixs) == 1))
+        return (PIX *)ERROR_PTR("pixs not defined or 1 bpp", __func__, NULL);
+    if (fract <= 0.0 || halfwidth <= 0) {
+        L_WARNING("no sharpening requested; clone returned\n", __func__);
+        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 ((pix1 = pixConvertTo8Or32(pixs, L_CLONE, 0)) == NULL)
+        return (PIX *)ERROR_PTR("pix1 not made", __func__, NULL);
+
+        /* Sharpen */
+    d = pixGetDepth(pix1);
+    if (d == 8) {
+        pixd = pixUnsharpMaskingGray(pix1, halfwidth, fract);
+    } else {  /* d == 32 */
+        pixr = pixGetRGBComponent(pix1, COLOR_RED);
+        pixrs = pixUnsharpMaskingGray(pixr, halfwidth, fract);
+        pixDestroy(&pixr);
+        pixg = pixGetRGBComponent(pix1, COLOR_GREEN);
+        pixgs = pixUnsharpMaskingGray(pixg, halfwidth, fract);
+        pixDestroy(&pixg);
+        pixb = pixGetRGBComponent(pix1, COLOR_BLUE);
+        pixbs = pixUnsharpMaskingGray(pixb, halfwidth, fract);
+        pixDestroy(&pixb);
+        pixd = pixCreateRGBImage(pixrs, pixgs, pixbs);
+        pixDestroy(&pixrs);
+        pixDestroy(&pixgs);
+        pixDestroy(&pixbs);
+        if (pixGetSpp(pixs) == 4)
+            pixCopyRGBComponent(pixd, pixs, L_ALPHA_CHANNEL);
+    }
+
+    pixDestroy(&pix1);
+    return pixd;
+}
+
+
+/*!
+ * \brief   pixUnsharpMaskingGray()
+ *
+ * \param[in]    pixs       8 bpp; no colormap
+ * \param[in]    halfwidth  "half-width" of smoothing filter
+ * \param[in]    fract      fraction of edge added back into image
+ * \return  pixd, or NULL on error
+ *
+ * <pre>
+ * 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.
+ * </pre>
+ */
+PIX *
+pixUnsharpMaskingGray(PIX       *pixs,
+                      l_int32    halfwidth,
+                      l_float32  fract)
+{
+l_int32  w, h, d;
+PIX     *pixc, *pixd;
+PIXACC  *pixacc;
+
+    if (!pixs)
+        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
+    pixGetDimensions(pixs, &w, &h, &d);
+    if (d != 8 || pixGetColormap(pixs) != NULL)
+        return (PIX *)ERROR_PTR("pixs not 8 bpp or has cmap", __func__, NULL);
+    if (fract <= 0.0 || halfwidth <= 0) {
+        L_WARNING("no sharpening requested; clone returned\n", __func__);
+        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", __func__, 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;
+}
+
+
+/*!
+ * \brief   pixUnsharpMaskingFast()
+ *
+ * \param[in]    pixs       all depths except 1 bpp; with or without colormaps
+ * \param[in]    halfwidth  "half-width" of smoothing filter; 1 and 2 only
+ * \param[in]    fract      fraction of high frequency added to image
+ * \param[in]    direction  L_HORIZ, L_VERT, L_BOTH_DIRECTIONS
+ * \return  pixd, or NULL on error
+ *
+ * <pre>
+ * 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.
+ * </pre>
+ */
+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;
+
+    if (!pixs || (pixGetDepth(pixs) == 1))
+        return (PIX *)ERROR_PTR("pixs not defined or 1 bpp", __func__, NULL);
+    if (fract <= 0.0 || halfwidth <= 0) {
+        L_WARNING("no sharpening requested; clone returned\n", __func__);
+        return pixClone(pixs);
+    }
+    if (halfwidth != 1 && halfwidth != 2)
+        return (PIX *)ERROR_PTR("halfwidth must be 1 or 2", __func__, NULL);
+    if (direction != L_HORIZ && direction != L_VERT &&
+        direction != L_BOTH_DIRECTIONS)
+        return (PIX *)ERROR_PTR("invalid direction", __func__, NULL);
+
+        /* Remove colormap; clone if possible; result is either 8 or 32 bpp */
+    if ((pixt = pixConvertTo8Or32(pixs, L_CLONE, 0)) == NULL)
+        return (PIX *)ERROR_PTR("pixt not made", __func__, 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);
+        if (pixGetSpp(pixs) == 4)
+            pixCopyRGBComponent(pixd, pixs, L_ALPHA_CHANNEL);
+        pixDestroy(&pixrs);
+        pixDestroy(&pixgs);
+        pixDestroy(&pixbs);
+    }
+
+    pixDestroy(&pixt);
+    return pixd;
+}
+
+
+
+/*!
+ * \brief   pixUnsharpMaskingGrayFast()
+ *
+ * \param[in]    pixs       8 bpp; no colormap
+ * \param[in]    halfwidth  "half-width" of smoothing filter: 1 or 2
+ * \param[in]    fract      fraction of high frequency added to image
+ * \param[in]    direction  L_HORIZ, L_VERT, L_BOTH_DIRECTIONS
+ * \return  pixd, or NULL on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) For usage and explanation of the algorithm, see notes
+ *          in pixUnsharpMaskingFast().
+ *      (2) Returns a clone if no sharpening is requested.
+ * </pre>
+ */
+PIX *
+pixUnsharpMaskingGrayFast(PIX       *pixs,
+                          l_int32    halfwidth,
+                          l_float32  fract,
+                          l_int32    direction)
+{
+PIX  *pixd;
+
+    if (!pixs)
+        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
+    if (pixGetDepth(pixs) != 8 || pixGetColormap(pixs) != NULL)
+        return (PIX *)ERROR_PTR("pixs not 8 bpp or has cmap", __func__, NULL);
+    if (fract <= 0.0 || halfwidth <= 0) {
+        L_WARNING("no sharpening requested; clone returned\n", __func__);
+        return pixClone(pixs);
+    }
+    if (halfwidth != 1 && halfwidth != 2)
+        return (PIX *)ERROR_PTR("halfwidth must be 1 or 2", __func__, NULL);
+    if (direction != L_HORIZ && direction != L_VERT &&
+        direction != L_BOTH_DIRECTIONS)
+        return (PIX *)ERROR_PTR("invalid direction", __func__, NULL);
+
+    if (direction != L_BOTH_DIRECTIONS)
+        pixd = pixUnsharpMaskingGray1D(pixs, halfwidth, fract, direction);
+    else  /* 2D sharpening */
+        pixd = pixUnsharpMaskingGray2D(pixs, halfwidth, fract);
+
+    return pixd;
+}
+
+
+/*!
+ * \brief   pixUnsharpMaskingGray1D()
+ *
+ * \param[in]    pixs        8 bpp; no colormap
+ * \param[in]    halfwidth   "half-width" of smoothing filter: 1 or 2
+ * \param[in]    fract       fraction of high frequency added to image
+ * \param[in]    direction   filtering direction; use L_HORIZ or L_VERT
+ * \return  pixd, or NULL on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) For usage and explanation of the algorithm, see notes
+ *          in pixUnsharpMaskingFast().
+ *      (2) Returns a clone if no sharpening is requested.
+ * </pre>
+ */
+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;
+
+    if (!pixs)
+        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
+    pixGetDimensions(pixs, &w, &h, &d);
+    if (d != 8 || pixGetColormap(pixs) != NULL)
+        return (PIX *)ERROR_PTR("pixs not 8 bpp or has cmap", __func__, NULL);
+    if (fract <= 0.0 || halfwidth <= 0) {
+        L_WARNING("no sharpening requested; clone returned\n", __func__);
+        return pixClone(pixs);
+    }
+    if (halfwidth != 1 && halfwidth != 2)
+        return (PIX *)ERROR_PTR("halfwidth must be 1 or 2", __func__, 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;
+}
+
+
+/*!
+ * \brief   pixUnsharpMaskingGray2D()
+ *
+ * \param[in]    pixs       8 bpp; no colormap
+ * \param[in]    halfwidth  "half-width" of smoothing filter: 1 or 2
+ * \param[in]    fract      fraction of high frequency added to image
+ * \return  pixd, or NULL on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) This is for %halfwidth == 1, 2.
+ *      (2) The lowpass filter is implemented separably.
+ *      (3) Returns a clone if no sharpening is requested.
+ * </pre>
+ */
+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, *lined;
+l_float32   val, norm;
+l_float32  *dataf, *linef, *linef0, *linef1, *linef2, *linef3, *linef4;
+PIX        *pixd;
+FPIX       *fpix;
+
+    if (!pixs)
+        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
+    pixGetDimensions(pixs, &w, &h, &d);
+    if (d != 8 || pixGetColormap(pixs) != NULL)
+        return (PIX *)ERROR_PTR("pixs not 8 bpp or has cmap", __func__, NULL);
+    if (fract <= 0.0 || halfwidth <= 0) {
+        L_WARNING("no sharpening requested; clone returned\n", __func__);
+        return pixClone(pixs);
+    }
+    if (halfwidth != 1 && halfwidth != 2)
+        return (PIX *)ERROR_PTR("halfwidth must be 1 or 2", __func__, NULL);
+
+    if ((pixd = pixCopyBorder(NULL, pixs, halfwidth, halfwidth,
+                              halfwidth, halfwidth)) == NULL)
+        return (PIX *)ERROR_PTR("pixd not made", __func__, NULL);
+    datad = pixGetData(pixd);
+    wpld = pixGetWpl(pixd);
+    datas = pixGetData(pixs);
+    wpls = pixGetWpl(pixs);
+
+        /* Do the low pass separably.  Store the result of horizontal
+         * smoothing in an intermediate fpix.  */
+    if ((fpix = fpixCreate(w, h)) == NULL) {
+        pixDestroy(&pixd);
+        return (PIX *)ERROR_PTR("fpix not made", __func__, NULL);
+    }
+    dataf = fpixGetData(fpix);
+    wplf = fpixGetWpl(fpix);
+    if (halfwidth == 1) {
+        for (i = 0; i < h; i++) {
+            lines = datas + i * wpls;
+            linef = dataf + i * wplf;
+            for (j = 1; j < w - 1; j++) {
+                val = GET_DATA_BYTE(lines, j - 1) +
+                      GET_DATA_BYTE(lines, j) +
+                      GET_DATA_BYTE(lines, j + 1);
+                linef[j] = val;
+            }
+        }
+    } else {
+        for (i = 0; i < h; 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.                 */
+    if (halfwidth == 1) {
+        for (i = 1; i < h - 1; i++) {
+            linef0 = dataf + (i - 1) * wplf;
+            linef1 = dataf + i * wplf;
+            linef2 = dataf + (i + 1) * wplf;
+            lined = datad + i * wpld;
+            lines = datas + i * wpls;
+            norm = 1.0f / 9.0f;
+            for (j = 1; j < w - 1; j++) {
+                val = norm * (linef0[j] + linef1[j] +
+                              linef2[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);
+            }
+        }
+    } else {
+        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;
+            norm = 1.0f / 25.0f;
+            for (j = 2; j < w - 2; j++) {
+                val = norm * (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                    *
+ *-----------------------------------------------------------------------*/
+/*!
+ * \brief   pixModifyHue()
+ *
+ * \param[in]    pixd      [optional] can be null or equal to pixs
+ * \param[in]    pixs      32 bpp rgb
+ * \param[in]    fract     between -1.0 and 1.0
+ * \return  pixd, or NULL on error
+ *
+ * <pre>
+ * 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) 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.
+ *      (3) 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.
+ *      (4) This leaves saturation and intensity invariant.
+ *      (5) See discussion of color-modification methods, in coloring.c.
+ * </pre>
+ */
+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;
+
+    if (!pixs)
+        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
+    if (pixGetColormap(pixs) != NULL)
+        return (PIX *)ERROR_PTR("pixs colormapped", __func__, NULL);
+    if (pixd && (pixd != pixs))
+        return (PIX *)ERROR_PTR("pixd not null or pixs", __func__, pixd);
+    pixGetDimensions(pixs, &w, &h, &d);
+    if (d != 32)
+        return (PIX *)ERROR_PTR("pixs not 32 bpp", __func__, NULL);
+    if (L_ABS(fract) > 1.0)
+        return (PIX *)ERROR_PTR("fract not in [-1.0 ... 1.0]", __func__, NULL);
+
+    pixd = pixCopy(pixd, pixs);
+
+    delhue = (l_int32)(240 * fract);
+    if (delhue == 0 || delhue == 240 || delhue == -240) {
+        L_WARNING("no change requested in hue\n", __func__);
+        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);
+        }
+    }
+    if (pixGetSpp(pixs) == 4)
+        pixCopyRGBComponent(pixd, pixs, L_ALPHA_CHANNEL);
+
+    return pixd;
+}
+
+
+/*!
+ * \brief   pixModifySaturation()
+ *
+ * \param[in]    pixd     [optional] can be null, existing or equal to pixs
+ * \param[in]    pixs     32 bpp rgb
+ * \param[in]    fract    between -1.0 and 1.0
+ * \return  pixd, or NULL on error
+ *
+ * <pre>
+ * 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.
+ *      (3) This leaves hue and intensity invariant.
+ *      (4) See discussion of color-modification methods, in coloring.c.
+ * </pre>
+ */
+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;
+
+    if (!pixs)
+        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
+    pixGetDimensions(pixs, &w, &h, &d);
+    if (d != 32)
+        return (PIX *)ERROR_PTR("pixs not 32 bpp", __func__, NULL);
+    if (L_ABS(fract) > 1.0)
+        return (PIX *)ERROR_PTR("fract not in [-1.0 ... 1.0]", __func__, NULL);
+
+    pixd = pixCopy(pixd, pixs);
+    if (fract == 0.0) {
+        L_WARNING("no change requested in saturation\n", __func__);
+        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);
+        }
+    }
+    if (pixGetSpp(pixs) == 4)
+        pixCopyRGBComponent(pixd, pixs, L_ALPHA_CHANNEL);
+
+    return pixd;
+}
+
+
+/*!
+ * \brief   pixMeasureSaturation()
+ *
+ * \param[in]    pixs     32 bpp rgb
+ * \param[in]    factor   subsampling factor; integer >= 1
+ * \param[out]   psat     average saturation
+ * \return  0 if OK, 1 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;
+
+    if (!psat)
+        return ERROR_INT("pixs not defined", __func__, 1);
+    *psat = 0.0;
+    if (!pixs)
+        return ERROR_INT("pixs not defined", __func__, 1);
+    pixGetDimensions(pixs, &w, &h, &d);
+    if (d != 32)
+        return ERROR_INT("pixs not 32 bpp", __func__, 1);
+    if (factor < 1)
+        return ERROR_INT("subsampling factor < 1", __func__, 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;
+}
+
+
+/*!
+ * \brief   pixModifyBrightness()
+ *
+ * \param[in]    pixd     [optional] can be null, existing or equal to pixs
+ * \param[in]    pixs     32 bpp rgb
+ * \param[in]    fract    between -1.0 and 1.0
+ * \return  pixd, or NULL on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) If fract > 0.0, it gives the fraction that the v-parameter,
+ *          which is max(r,g,b), is moved from its initial value toward 255.
+ *          If fract < 0.0, it gives the fraction that the v-parameter
+ *          is moved from its initial value toward 0.
+ *          The limiting values for fract = -1.0 (1.0) thus set the
+ *          v-parameter 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.
+ *      (3) This leaves hue and saturation invariant.
+ *      (4) See discussion of color-modification methods, in coloring.c.
+ * </pre>
+ */
+PIX  *
+pixModifyBrightness(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;
+
+    if (!pixs)
+        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
+    pixGetDimensions(pixs, &w, &h, &d);
+    if (d != 32)
+        return (PIX *)ERROR_PTR("pixs not 32 bpp", __func__, NULL);
+    if (L_ABS(fract) > 1.0)
+        return (PIX *)ERROR_PTR("fract not in [-1.0 ... 1.0]", __func__, NULL);
+
+    pixd = pixCopy(pixd, pixs);
+    if (fract == 0.0) {
+        L_WARNING("no change requested in brightness\n", __func__);
+        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)
+                vval = (l_int32)(vval + fract * (255.0 - vval));
+            else
+                vval = (l_int32)(vval * (1.0 + fract));
+            convertHSVToRGB(hval, sval, vval, &rval, &gval, &bval);
+            composeRGBPixel(rval, gval, bval, line + j);
+        }
+    }
+    if (pixGetSpp(pixs) == 4)
+        pixCopyRGBComponent(pixd, pixs, L_ALPHA_CHANNEL);
+
+    return pixd;
+}
+
+
+/*-----------------------------------------------------------------------*
+ *                             Color shifting                            *
+ *-----------------------------------------------------------------------*/
+/*!
+ * \brief   pixMosaicColorShiftRGB()
+ *
+ * \param[in]    pixs     32 bpp rgb
+ * \param[in]    roff   center offset of red component
+ * \param[in]    goff   center offset of green component
+ * \param[in]    boff   center offset of blue component
+ * \param[in]    delta  increments from center offsets [0.0 - 0.1];
+ *                      use 0.0 to get the default (0.04)
+ * \param[in]    nincr  number of increments in each (positive and negative)
+ *                      direction; use 0 to get the default (2).
+ * \return  pix, or NULL on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) This generates a mosaic view of the effect of shifting the RGB
+ *          components.  See pixColorShiftRGB() for details on the shifting.
+ *      (2) The offsets (%roff, %goff, %boff) set the color center point,
+ *          and the deviations from this are shown separately for deltas
+ *          in r, g and b.  For each component, we show 2 * %nincr + 1
+ *          images.
+ *      (3) The pix must have minimum dimensions of 100 and an aspect
+ *          ratio not exceeding 5.0.
+ *      (4) Usage: color prints differ from the original due to three factors:
+ *          illumination, calibration of the camera in acquisition,
+ *          and calibration of the printer.  This function can be used
+ *          to iteratively match a color print to the original.  On each
+ *          iteration, the center offsets are set to the best match so
+ *          far, and the %delta increments are typically reduced.
+ * </pre>
+ */
+PIX *
+pixMosaicColorShiftRGB(PIX       *pixs,
+                       l_float32  roff,
+                       l_float32  goff,
+                       l_float32  boff,
+                       l_float32  delta,
+                       l_int32    nincr)
+{
+char       buf[64];
+l_int32    i, w, h;
+l_float32  del, ratio;
+L_BMF     *bmf;
+PIX       *pix1, *pix2, *pix3;
+PIXA      *pixa;
+
+    if (!pixs  || pixGetDepth(pixs) != 32)
+        return (PIX *)ERROR_PTR("pixs undefined or not rgb", __func__, NULL);
+    if (roff < -1.0 || roff > 1.0)
+        return (PIX *)ERROR_PTR("roff not in [-1.0, 1.0]", __func__, NULL);
+    if (goff < -1.0 || goff > 1.0)
+        return (PIX *)ERROR_PTR("goff not in [-1.0, 1.0]", __func__, NULL);
+    if (boff < -1.0 || boff > 1.0)
+        return (PIX *)ERROR_PTR("boff not in [-1.0, 1.0]", __func__, NULL);
+    if (delta < 0.0 || delta > 0.1)
+        return (PIX *)ERROR_PTR("delta not in [0.0, 0.1]", __func__, NULL);
+    if (delta == 0.0) delta = 0.04f;
+    if (nincr < 0 || nincr > 6)
+        return (PIX *)ERROR_PTR("nincr not in [0, 6]", __func__, NULL);
+    if (nincr == 0) nincr = 2;
+
+        /* Require width and height to be >= 100, and the aspect ratio <= 5.0 */
+    pixGetDimensions(pixs, &w, &h, NULL);
+    if (w < 100 || h < 100)
+        return (PIX *)ERROR_PTR("w and h not both >= 100", __func__, NULL);
+    pixMaxAspectRatio(pixs, &ratio);
+    if (ratio < 1.0 || ratio > 5.0) {
+        L_ERROR("invalid aspect ratio %5.1f\n", __func__, ratio);
+        return NULL;
+    }
+
+    pixa = pixaCreate(3 * (2 * nincr + 1));
+    bmf = bmfCreate(NULL, 8);
+    pix1 = pixScaleToSize(pixs, 400, 0);
+    for (i = 0, del = - nincr * delta; i < 2 * nincr + 1; i++, del += delta) {
+        pix2 = pixColorShiftRGB(pix1, roff + del, goff, boff);
+        snprintf(buf, sizeof(buf), "%4.2f, %4.2f, %4.2f",
+                 roff + del, goff, boff);
+        pix3 = pixAddSingleTextblock(pix2, bmf, buf, 0xff000000,
+                                     L_ADD_BELOW, 0);
+        pixaAddPix(pixa, pix3, L_INSERT);
+        pixDestroy(&pix2);
+    }
+    for (i = 0, del = - nincr * delta; i < 2 * nincr + 1; i++, del += delta) {
+        pix2 = pixColorShiftRGB(pix1, roff, goff + del, boff);
+        snprintf(buf, sizeof(buf), "%4.2f, %4.2f, %4.2f",
+                 roff, goff + del, boff);
+        pix3 = pixAddSingleTextblock(pix2, bmf, buf, 0xff000000,
+                                     L_ADD_BELOW, 0);
+        pixaAddPix(pixa, pix3, L_INSERT);
+        pixDestroy(&pix2);
+    }
+    for (i = 0, del = - nincr * delta; i < 2 * nincr + 1; i++, del += delta) {
+        pix2 = pixColorShiftRGB(pix1, roff, goff, boff + del);
+        snprintf(buf, sizeof(buf), "%4.2f, %4.2f, %4.2f",
+                 roff, goff, boff + del);
+        pix3 = pixAddSingleTextblock(pix2, bmf, buf, 0xff000000,
+                                     L_ADD_BELOW, 0);
+        pixaAddPix(pixa, pix3, L_INSERT);
+        pixDestroy(&pix2);
+    }
+    pixDestroy(&pix1);
+
+    pix1 = pixaDisplayTiledAndScaled(pixa, 32, 300, 2 * nincr + 1, 0, 30, 2);
+    pixaDestroy(&pixa);
+    bmfDestroy(&bmf);
+    return pix1;
+}
+
+
+/*!
+ * \brief   pixColorShiftRGB()
+ *
+ * \param[in]    pixs     32 bpp rgb
+ * \param[in]    rfract   fractional shift in red component
+ * \param[in]    gfract   fractional shift in green component
+ * \param[in]    bfract   fractional shift in blue component
+ * \return  pixd, or NULL on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) This allows independent fractional shifts of the r,g and b
+ *          components.  A positive shift pushes to saturation (255);
+ *          a negative shift pushes toward 0 (black).
+ *      (2) The effect can be imagined using a color wheel that consists
+ *          (for our purposes) of these 6 colors, separated by 60 degrees:
+ *             red, magenta, blue, cyan, green, yellow
+ *      (3) So, for example, a negative shift of the blue component
+ *          (bfract < 0) could be accompanied by positive shifts
+ *          of red and green to make an image more yellow.
+ *      (4) Examples of limiting cases:
+ *            rfract = 1 ==> r = 255
+ *            rfract = -1 ==> r = 0
+ * </pre>
+ */
+PIX *
+pixColorShiftRGB(PIX       *pixs,
+                 l_float32  rfract,
+                 l_float32  gfract,
+                 l_float32  bfract)
+{
+l_int32    w, h, i, j, wpls, wpld, rval, gval, bval;
+l_int32   *rlut, *glut, *blut;
+l_uint32  *datas, *datad, *lines, *lined;
+l_float32  fi;
+PIX       *pixd;
+
+    if (!pixs)
+        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
+    if (pixGetDepth(pixs) != 32)
+        return (PIX *)ERROR_PTR("pixs not 32 bpp", __func__, NULL);
+    if (rfract < -1.0 || rfract > 1.0)
+        return (PIX *)ERROR_PTR("rfract not in [-1.0, 1.0]", __func__, NULL);
+    if (gfract < -1.0 || gfract > 1.0)
+        return (PIX *)ERROR_PTR("gfract not in [-1.0, 1.0]", __func__, NULL);
+    if (bfract < -1.0 || bfract > 1.0)
+        return (PIX *)ERROR_PTR("bfract not in [-1.0, 1.0]", __func__, NULL);
+    if (rfract == 0.0 && gfract == 0.0 && bfract == 0.0)
+        return pixCopy(NULL, pixs);
+
+    rlut = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32));
+    glut = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32));
+    blut = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32));
+    for (i = 0; i < 256; i++) {
+        fi = i;
+        if (rfract >= 0) {
+            rlut[i] = (l_int32)(fi + (255.0 - fi) * rfract);
+        } else {
+            rlut[i] = (l_int32)(fi * (1.0 + rfract));
+        }
+        if (gfract >= 0) {
+            glut[i] = (l_int32)(fi + (255.0 - fi) * gfract);
+        } else {
+            glut[i] = (l_int32)(fi * (1.0 + gfract));
+        }
+        if (bfract >= 0) {
+            blut[i] = (l_int32)(fi + (255.0 - fi) * bfract);
+        } else {
+            blut[i] = (l_int32)(fi * (1.0 + bfract));
+        }
+    }
+
+    pixGetDimensions(pixs, &w, &h, NULL);
+    datas = pixGetData(pixs);
+    wpls = pixGetWpl(pixs);
+    pixd = pixCreate(w, h, 32);
+    datad = pixGetData(pixd);
+    wpld = pixGetWpl(pixd);
+    for (i = 0; i < h; i++) {
+        lines = datas + i * wpls;
+        lined = datad + i * wpld;
+        for (j = 0; j < w; j++) {
+            extractRGBValues(lines[j], &rval, &gval, &bval);
+            composeRGBPixel(rlut[rval], glut[gval], blut[bval], lined + j);
+        }
+    }
+
+    LEPT_FREE(rlut);
+    LEPT_FREE(glut);
+    LEPT_FREE(blut);
+    return pixd;
+}
+
+/*-----------------------------------------------------------------------*
+ *                     Darken gray (unsaturated) pixels
+ *-----------------------------------------------------------------------*/
+/*!
+ * \brief   pixDarkenGray()
+ *
+ * \param[in]    pixd      [optional] can be null or equal to pixs
+ * \param[in]    pixs      32 bpp rgb
+ * \param[in]    thresh    pixels with max component >= %thresh are unchanged
+ * \param[in]    satlimit  pixels with saturation >= %satlimit are unchanged
+ * \return  pixd, or NULL on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) This darkens gray pixels, by a fraction (sat/%satlimit), where
+ *          the saturation, sat, is the component difference (max - min).
+ *          The pixel value is unchanged if sat >= %satlimit.  A typical
+ *          value of %satlimit might be 40; the larger the value, the
+ *          more that pixels with a smaller saturation will be darkened.
+ *      (2) Pixels with max component >= %thresh are unchanged. This can be
+ *          used to prevent bright pixels with low saturation from being
+ *          darkened.  Setting thresh == 0 is a no-op; setting %thresh == 255
+ *          causes the darkening to be applied to all pixels.
+ *      (3) This function is useful to enhance pixels relative to a
+ *          gray background.
+ *      (4) A related function that builds a 1 bpp mask over the gray
+ *          pixels is pixMaskOverGrayPixels().
+ * </pre>
+ */
+PIX *
+pixDarkenGray(PIX     *pixd,
+              PIX     *pixs,
+              l_int32  thresh,
+              l_int32  satlimit)
+{
+l_int32    w, h, i, j, wpls, wpld;
+l_int32    rval, gval, bval, minrg, min, maxrg, max, sat;
+l_uint32  *datas, *datad, *lines, *lined;
+l_float32  ratio;
+
+    if (!pixs || pixGetDepth(pixs) != 32)
+        return (PIX *)ERROR_PTR("pixs undefined or not 32 bpp", __func__, NULL);
+    if (thresh < 0 || thresh > 255)
+        return (PIX *)ERROR_PTR("invalid thresh", __func__, NULL);
+    if (satlimit < 1)
+        return (PIX *)ERROR_PTR("invalid satlimit", __func__, NULL);
+    if (pixd && (pixs != pixd))
+        return (PIX *)ERROR_PTR("not new or in-place", __func__, NULL);
+
+    pixGetDimensions(pixs, &w, &h, NULL);
+    datas = pixGetData(pixs);
+    wpls = pixGetWpl(pixs);
+    if ((pixd = pixCopy(pixd, pixs)) == NULL)
+        return (PIX *)ERROR_PTR("pixd not made", __func__, NULL);
+    datad = pixGetData(pixd);
+    wpld = pixGetWpl(pixd);
+
+    for (i = 0; i < h; i++) {
+        lines = datas + i * wpls;
+        lined = datad + i * wpld;
+        for (j = 0; j < w; j++) {
+            extractRGBValues(lines[j], &rval, &gval, &bval);
+            minrg = L_MIN(rval, gval);
+            min = L_MIN(minrg, bval);
+            maxrg = L_MAX(rval, gval);
+            max = L_MAX(maxrg, bval);
+            sat = max - min;
+            if (max >= thresh || sat >= satlimit)
+                continue;
+            ratio = (l_float32)sat / (l_float32)satlimit;
+            composeRGBPixel((l_int32)(ratio * rval), (l_int32)(ratio * gval),
+                            (l_int32)(ratio * bval), &lined[j]);
+        }
+    }
+    return pixd;
+}
+
+
+/*-----------------------------------------------------------------------*
+ *            General multiplicative constant color transform            *
+ *-----------------------------------------------------------------------*/
+/*!
+ * \brief   pixMultConstantColor()
+ *
+ * \param[in]    pixs     colormapped or rgb
+ * \param[in]    rfact    red multiplicative factor
+ * \param[in]    gfact    green multiplicative factor
+ * \param[in]    bfact    blue multiplicative factor
+ * \return  pixd colormapped or rgb, with colors scaled, or NULL on error
+ *
+ * <pre>
+ * 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().
+ * </pre>
+ */
+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;
+
+    if (!pixs)
+        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
+    pixGetDimensions(pixs, &w, &h, &d);
+    cmap = pixGetColormap(pixs);
+    if (!cmap && d != 32)
+        return (PIX *)ERROR_PTR("pixs not cmapped or 32 bpp", __func__, 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", __func__, 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 = pixCreateTemplate(pixs)) == NULL)
+        return (PIX *)ERROR_PTR("pixd not made", __func__, 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;
+}
+
+
+/*!
+ * \brief   pixMultMatrixColor()
+ *
+ * \param[in]    pixs    colormapped or rgb
+ * \param[in]    kel     kernel 3x3 matrix of floats
+ * \return  pixd colormapped or rgb, or NULL on error
+ *
+ * <pre>
+ * 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].
+ * </pre>
+ */
+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;
+
+    if (!pixs)
+        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
+    if (!kel)
+        return (PIX *)ERROR_PTR("kel not defined", __func__, NULL);
+    kernelGetParameters(kel, &kw, &kh, NULL, NULL);
+    if (kw != 3 || kh != 3)
+        return (PIX *)ERROR_PTR("matrix not 3x3", __func__, NULL);
+    pixGetDimensions(pixs, &w, &h, &d);
+    cmap = pixGetColormap(pixs);
+    if (!cmap && d != 32)
+        return (PIX *)ERROR_PTR("pixs not cmapped or 32 bpp", __func__, 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", __func__, 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 = pixCreateTemplate(pixs)) == NULL)
+        return (PIX *)ERROR_PTR("pixd not made", __func__, 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                    *
+ *-------------------------------------------------------------*/
+/*!
+ * \brief   pixHalfEdgeByBandpass()
+ *
+ * \param[in]    pixs         8 bpp gray or 32 bpp rgb
+ * \param[in]    sm1h, sm1v   "half-widths" of smoothing filter sm1
+ * \param[in]    sm2h, sm2v   "half-widths" of smoothing filter sm2;
+ *                            require sm2 != sm1
+ * \return  pixd, or NULL on error
+ *
+ * <pre>
+ * 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.
+ * </pre>
+ */
+PIX *
+pixHalfEdgeByBandpass(PIX     *pixs,
+                      l_int32  sm1h,
+                      l_int32  sm1v,
+                      l_int32  sm2h,
+                      l_int32  sm2v)
+{
+l_int32  d;
+PIX     *pixg, *pixacc, *pixc1, *pixc2;
+
+    if (!pixs)
+        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
+    if (sm1h == sm2h && sm1v == sm2v)
+        return (PIX *)ERROR_PTR("sm2 = sm1", __func__, NULL);
+    d = pixGetDepth(pixs);
+    if (d != 8 && d != 32)
+        return (PIX *)ERROR_PTR("pixs not 8 or 32 bpp", __func__, 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) {
+        pixDestroy(&pixg);
+        return (PIX *)ERROR_PTR("pixacc not made", __func__, NULL);
+    }
+    if ((pixc1 = pixBlockconvGray(pixg, pixacc, sm1h, sm1v)) == NULL) {
+        pixDestroy(&pixg);
+        pixDestroy(&pixacc);
+        return (PIX *)ERROR_PTR("pixc1 not made", __func__, NULL);
+    }
+    pixc2 = pixBlockconvGray(pixg, pixacc, sm2h, sm2v);
+    pixDestroy(&pixg);
+    pixDestroy(&pixacc);
+    if (!pixc2) {
+        pixDestroy(&pixc1);
+        return (PIX *)ERROR_PTR("pixc2 not made", __func__, NULL);
+    }
+
+        /* Compute the half-edge using pixc1 - pixc2.  */
+    pixSubtractGray(pixc1, pixc1, pixc2);
+    pixDestroy(&pixc2);
+    return pixc1;
+}