--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mupdf-source/thirdparty/leptonica/src/pixarith.c	Mon Sep 15 11:43:07 2025 +0200
@@ -0,0 +1,1597 @@
+/*====================================================================*
+ -  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 pixarith.c
+ * <pre>
+ *
+ *      One-image grayscale arithmetic operations (8, 16, 32 bpp)
+ *           l_int32     pixAddConstantGray()
+ *           l_int32     pixMultConstantGray()
+ *
+ *      Two-image grayscale arithmetic operations (8, 16, 32 bpp)
+ *           PIX        *pixAddGray()
+ *           PIX        *pixSubtractGray()
+ *           PIX        *pixMultiplyGray()
+ *
+ *      Grayscale threshold operation (8, 16, 32 bpp)
+ *           PIX        *pixThresholdToValue()
+ *
+ *      Image accumulator arithmetic operations
+ *           PIX        *pixInitAccumulate()
+ *           PIX        *pixFinalAccumulate()
+ *           PIX        *pixFinalAccumulateThreshold()
+ *           l_int32     pixAccumulate()
+ *           l_int32     pixMultConstAccumulate()
+ *
+ *      Absolute value of difference
+ *           PIX        *pixAbsDifference()
+ *
+ *      Sum of color images
+ *           PIX        *pixAddRGB()
+ *
+ *      Two-image min and max operations (8 and 16 bpp)
+ *           PIX        *pixMinOrMax()
+ *
+ *      Scale pix for maximum dynamic range
+ *           PIX        *pixMaxDynamicRange()
+ *           PIX        *pixMaxDynamicRangeRGB()
+ *
+ *      RGB pixel value scaling
+ *           l_uint32    linearScaleRGBVal()
+ *           l_uint32    logScaleRGBVal()
+ *
+ *      Log base2 lookup
+ *           l_float32  *makeLogBase2Tab()
+ *           l_float32   getLogBase2()
+ *
+ *      The image accumulator operations are used when you expect
+ *      overflow from 8 bits on intermediate results.  For example,
+ *      you might want a tophat contrast operator which is
+ *         3*I - opening(I,S) - closing(I,S)
+ *      To use these operations, first use the init to generate
+ *      a 16 bpp image, use the accumulate to add or subtract 8 bpp
+ *      images from that, or the multiply constant to multiply
+ *      by a small constant (much less than 256 -- we don't want
+ *      overflow from the 16 bit images!), and when you're finished
+ *      use final to bring the result back to 8 bpp, clipped
+ *      if necessary.  There is also a divide function, which
+ *      can be used to divide one image by another, scaling the
+ *      result for maximum dynamic range, and giving back the
+ *      8 bpp result.
+ *
+ *      A simpler interface to the arithmetic operations is
+ *      provided in pixacc.c.
+ * </pre>
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config_auto.h>
+#endif  /* HAVE_CONFIG_H */
+
+#include <string.h>
+#include <math.h>
+#include "allheaders.h"
+
+/*-------------------------------------------------------------*
+ *          One-image grayscale arithmetic operations          *
+ *-------------------------------------------------------------*/
+/*!
+ * \brief   pixAddConstantGray()
+ *
+ * \param[in]    pixs   8, 16 or 32 bpp
+ * \param[in]    val    amount to add to each pixel
+ * \return  0 if OK, 1 on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) In-place operation.
+ *      (2) No clipping for 32 bpp.
+ *      (3) For 8 and 16 bpp, if val > 0 the result is clipped
+ *          to 0xff and 0xffff, rsp.
+ *      (4) For 8 and 16 bpp, if val < 0 the result is clipped to 0.
+ * </pre>
+ */
+l_ok
+pixAddConstantGray(PIX      *pixs,
+                   l_int32   val)
+{
+l_int32    i, j, w, h, d, wpl, pval;
+l_uint32  *data, *line;
+
+    if (!pixs)
+        return ERROR_INT("pixs not defined", __func__, 1);
+    pixGetDimensions(pixs, &w, &h, &d);
+    if (d != 8 && d != 16 && d != 32)
+        return ERROR_INT("pixs not 8, 16 or 32 bpp", __func__, 1);
+
+    data = pixGetData(pixs);
+    wpl = pixGetWpl(pixs);
+    for (i = 0; i < h; i++) {
+        line = data + i * wpl;
+        if (d == 8) {
+            if (val < 0) {
+                for (j = 0; j < w; j++) {
+                    pval = GET_DATA_BYTE(line, j);
+                    pval = L_MAX(0, pval + val);
+                    SET_DATA_BYTE(line, j, pval);
+                }
+            } else {  /* val >= 0 */
+                for (j = 0; j < w; j++) {
+                    pval = GET_DATA_BYTE(line, j);
+                    pval = L_MIN(255, pval + val);
+                    SET_DATA_BYTE(line, j, pval);
+                }
+            }
+        } else if (d == 16) {
+            if (val < 0) {
+                for (j = 0; j < w; j++) {
+                    pval = GET_DATA_TWO_BYTES(line, j);
+                    pval = L_MAX(0, pval + val);
+                    SET_DATA_TWO_BYTES(line, j, pval);
+                }
+            } else {  /* val >= 0 */
+                for (j = 0; j < w; j++) {
+                    pval = GET_DATA_TWO_BYTES(line, j);
+                    pval = L_MIN(0xffff, pval + val);
+                    SET_DATA_TWO_BYTES(line, j, pval);
+                }
+            }
+        } else {  /* d == 32; no check for overflow (< 0 or > 0xffffffff) */
+            for (j = 0; j < w; j++)
+                *(line + j) += val;
+        }
+    }
+
+    return 0;
+}
+
+
+/*!
+ * \brief   pixMultConstantGray()
+ *
+ * \param[in]    pixs   8, 16 or 32 bpp
+ * \param[in]    val    >= 0.0; amount to multiply by each pixel
+ * \return  0 if OK, 1 on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) In-place operation; val must be >= 0.
+ *      (2) No clipping for 32 bpp.
+ *      (3) For 8 and 16 bpp, the result is clipped to 0xff and 0xffff, rsp.
+ * </pre>
+ */
+l_ok
+pixMultConstantGray(PIX       *pixs,
+                    l_float32  val)
+{
+l_int32    i, j, w, h, d, wpl, pval;
+l_uint32   upval;
+l_uint32  *data, *line;
+
+    if (!pixs)
+        return ERROR_INT("pixs not defined", __func__, 1);
+    pixGetDimensions(pixs, &w, &h, &d);
+    if (d != 8 && d != 16 && d != 32)
+        return ERROR_INT("pixs not 8, 16 or 32 bpp", __func__, 1);
+    if (val < 0.0)
+        return ERROR_INT("val < 0.0", __func__, 1);
+
+    data = pixGetData(pixs);
+    wpl = pixGetWpl(pixs);
+    for (i = 0; i < h; i++) {
+        line = data + i * wpl;
+        if (d == 8) {
+            for (j = 0; j < w; j++) {
+                pval = GET_DATA_BYTE(line, j);
+                pval = (l_int32)(val * pval);
+                pval = L_MIN(255, pval);
+                SET_DATA_BYTE(line, j, pval);
+            }
+        } else if (d == 16) {
+            for (j = 0; j < w; j++) {
+                pval = GET_DATA_TWO_BYTES(line, j);
+                pval = (l_int32)(val * pval);
+                pval = L_MIN(0xffff, pval);
+                SET_DATA_TWO_BYTES(line, j, pval);
+            }
+        } else {  /* d == 32; no clipping */
+            for (j = 0; j < w; j++) {
+                upval = *(line + j);
+                upval = (l_uint32)(val * upval);
+                *(line + j) = upval;
+            }
+        }
+    }
+
+    return 0;
+}
+
+
+/*-------------------------------------------------------------*
+ *             Two-image grayscale arithmetic ops              *
+ *-------------------------------------------------------------*/
+/*!
+ * \brief   pixAddGray()
+ *
+ * \param[in]    pixd    [optional]; this can be null, equal to pixs1, or
+ *                       different from pixs1
+ * \param[in]    pixs1   can be equal to pixd
+ * \param[in]    pixs2
+ * \return  pixd always
+ *
+ * <pre>
+ * Notes:
+ *      (1) Arithmetic addition of two 8, 16 or 32 bpp images.
+ *      (2) For 8 and 16 bpp, we do explicit clipping to 0xff and 0xffff,
+ *          respectively.
+ *      (3) Alignment is to UL corner.
+ *      (4) There are 3 cases.  The result can go to a new dest,
+ *          in-place to pixs1, or to an existing input dest:
+ *          * pixd == null:   (src1 + src2) --> new pixd
+ *          * pixd == pixs1:  (src1 + src2) --> src1  (in-place)
+ *          * pixd != pixs1:  (src1 + src2) --> input pixd
+ *      (5) pixs2 must be different from both pixd and pixs1.
+ * </pre>
+ */
+PIX *
+pixAddGray(PIX  *pixd,
+           PIX  *pixs1,
+           PIX  *pixs2)
+{
+l_int32    i, j, d, ws, hs, w, h, wpls, wpld, val, sum;
+l_uint32  *datas, *datad, *lines, *lined;
+
+    if (!pixs1)
+        return (PIX *)ERROR_PTR("pixs1 not defined", __func__, pixd);
+    if (!pixs2)
+        return (PIX *)ERROR_PTR("pixs2 not defined", __func__, pixd);
+    if (pixs2 == pixs1)
+        return (PIX *)ERROR_PTR("pixs2 and pixs1 must differ", __func__, pixd);
+    if (pixs2 == pixd)
+        return (PIX *)ERROR_PTR("pixs2 and pixd must differ", __func__, pixd);
+    d = pixGetDepth(pixs1);
+    if (d != 8 && d != 16 && d != 32)
+        return (PIX *)ERROR_PTR("pix are not 8, 16 or 32 bpp", __func__, pixd);
+    if (pixGetDepth(pixs2) != d)
+        return (PIX *)ERROR_PTR("depths differ (pixs1, pixs2)", __func__, pixd);
+    if (pixd && (pixGetDepth(pixd) != d))
+        return (PIX *)ERROR_PTR("depths differ (pixs1, pixd)", __func__, pixd);
+
+    if (!pixSizesEqual(pixs1, pixs2))
+        L_WARNING("pixs1 and pixs2 not equal in size\n", __func__);
+    if (pixd && !pixSizesEqual(pixs1, pixd))
+        L_WARNING("pixs1 and pixd not equal in size\n", __func__);
+
+    if (pixs1 != pixd)
+        pixd = pixCopy(pixd, pixs1);
+
+        /* pixd + pixs2 ==> pixd  */
+    datas = pixGetData(pixs2);
+    datad = pixGetData(pixd);
+    wpls = pixGetWpl(pixs2);
+    wpld = pixGetWpl(pixd);
+    pixGetDimensions(pixs2, &ws, &hs, NULL);
+    pixGetDimensions(pixd, &w, &h, NULL);
+    w = L_MIN(ws, w);
+    h = L_MIN(hs, h);
+    for (i = 0; i < h; i++) {
+        lined = datad + i * wpld;
+        lines = datas + i * wpls;
+        if (d == 8) {
+            for (j = 0; j < w; j++) {
+                sum = GET_DATA_BYTE(lines, j) + GET_DATA_BYTE(lined, j);
+                val = L_MIN(sum, 255);
+                SET_DATA_BYTE(lined, j, val);
+            }
+        } else if (d == 16) {
+            for (j = 0; j < w; j++) {
+                sum = GET_DATA_TWO_BYTES(lines, j)
+                    + GET_DATA_TWO_BYTES(lined, j);
+                val = L_MIN(sum, 0xffff);
+                SET_DATA_TWO_BYTES(lined, j, val);
+            }
+        } else {   /* d == 32; no clipping */
+            for (j = 0; j < w; j++)
+                *(lined + j) += *(lines + j);
+        }
+    }
+
+    return pixd;
+}
+
+
+/*!
+ * \brief   pixSubtractGray()
+ *
+ * \param[in]    pixd     [optional]; this can be null, equal to pixs1, or
+ *                        different from pixs1
+ * \param[in]    pixs1    can be equal to pixd
+ * \param[in]    pixs2
+ * \return  pixd always
+ *
+ * <pre>
+ * Notes:
+ *      (1) Arithmetic subtraction of two 8, 16 or 32 bpp images.
+ *      (2) Source pixs2 is always subtracted from source pixs1.
+ *      (3) Do explicit clipping to 0.
+ *      (4) Alignment is to UL corner.
+ *      (5) There are 3 cases.  The result can go to a new dest,
+ *          in-place to pixs1, or to an existing input dest:
+ *          (a) pixd == null   (src1 - src2) --> new pixd
+ *          (b) pixd == pixs1  (src1 - src2) --> src1  (in-place)
+ *          (d) pixd != pixs1  (src1 - src2) --> input pixd
+ *      (6) pixs2 must be different from both pixd and pixs1.
+ * </pre>
+ */
+PIX *
+pixSubtractGray(PIX  *pixd,
+                PIX  *pixs1,
+                PIX  *pixs2)
+{
+l_int32    i, j, w, h, ws, hs, d, wpls, wpld, val, diff;
+l_uint32  *datas, *datad, *lines, *lined;
+
+    if (!pixs1)
+        return (PIX *)ERROR_PTR("pixs1 not defined", __func__, pixd);
+    if (!pixs2)
+        return (PIX *)ERROR_PTR("pixs2 not defined", __func__, pixd);
+    if (pixs2 == pixs1)
+        return (PIX *)ERROR_PTR("pixs2 and pixs1 must differ", __func__, pixd);
+    if (pixs2 == pixd)
+        return (PIX *)ERROR_PTR("pixs2 and pixd must differ", __func__, pixd);
+    d = pixGetDepth(pixs1);
+    if (d != 8 && d != 16 && d != 32)
+        return (PIX *)ERROR_PTR("pix are not 8, 16 or 32 bpp", __func__, pixd);
+    if (pixGetDepth(pixs2) != d)
+        return (PIX *)ERROR_PTR("depths differ (pixs1, pixs2)", __func__, pixd);
+    if (pixd && (pixGetDepth(pixd) != d))
+        return (PIX *)ERROR_PTR("depths differ (pixs1, pixd)", __func__, pixd);
+
+    if (!pixSizesEqual(pixs1, pixs2))
+        L_WARNING("pixs1 and pixs2 not equal in size\n", __func__);
+    if (pixd && !pixSizesEqual(pixs1, pixd))
+        L_WARNING("pixs1 and pixd not equal in size\n", __func__);
+
+    if (pixs1 != pixd)
+        pixd = pixCopy(pixd, pixs1);
+
+        /* pixd - pixs2 ==> pixd  */
+    datas = pixGetData(pixs2);
+    datad = pixGetData(pixd);
+    wpls = pixGetWpl(pixs2);
+    wpld = pixGetWpl(pixd);
+    pixGetDimensions(pixs2, &ws, &hs, NULL);
+    pixGetDimensions(pixd, &w, &h, NULL);
+    w = L_MIN(ws, w);
+    h = L_MIN(hs, h);
+    for (i = 0; i < h; i++) {
+        lined = datad + i * wpld;
+        lines = datas + i * wpls;
+        if (d == 8) {
+            for (j = 0; j < w; j++) {
+                diff = GET_DATA_BYTE(lined, j) - GET_DATA_BYTE(lines, j);
+                val = L_MAX(diff, 0);
+                SET_DATA_BYTE(lined, j, val);
+            }
+        } else if (d == 16) {
+            for (j = 0; j < w; j++) {
+                diff = GET_DATA_TWO_BYTES(lined, j)
+                       - GET_DATA_TWO_BYTES(lines, j);
+                val = L_MAX(diff, 0);
+                SET_DATA_TWO_BYTES(lined, j, val);
+            }
+        } else {  /* d == 32; no clipping */
+            for (j = 0; j < w; j++)
+                *(lined + j) -= *(lines + j);
+        }
+    }
+
+    return pixd;
+}
+
+
+/*!
+ * \brief   pixMultiplyGray()
+ *
+ * \param[in]    pixs    32 bpp rgb or 8 bpp gray
+ * \param[in]    pixg    8 bpp gray
+ * \param[in]    norm    multiplicative factor to avoid overflow; 0 for default
+ * \return  pixd, or null on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) This function can be used for correcting a scanned image
+ *          under non-uniform illumination.  For that application,
+ *          %pixs is the scanned image, %pixg is an image whose values
+ *          are inversely related to light from a uniform (say, white)
+ *          target, and %norm is typically the inverse of the maximum
+ *          pixel value in %pixg.
+ *      (2) Set norm = 0 to get the default value, which is the inverse
+ *          of the max value in %pixg.  This avoids overflow in the product.
+ *      (3) For 32 bpp %pixs, all 3 components are multiplied by the
+ *          same number.
+ *      (4) Alignment is to UL corner.
+ * </pre>
+ */
+PIX *
+pixMultiplyGray(PIX        *pixs,
+                PIX        *pixg,
+                l_float32   norm)
+{
+l_int32    i, j, w, h, d, ws, hs, ds, wpls, wplg, wpld;
+l_int32    rval, gval, bval, rval2, gval2, bval2, vals, valg, val, maxgray;
+l_uint32   val32;
+l_uint32  *datas, *datag, *datad, *lines, *lineg, *lined;
+PIX       *pixd;
+
+    if (!pixs)
+        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
+    pixGetDimensions(pixs, &ws, &hs, &ds);
+    if (ds != 8 && ds != 32)
+        return (PIX *)ERROR_PTR("pixs not 8 or 32 bpp", __func__, NULL);
+    if (!pixg)
+        return (PIX *)ERROR_PTR("pixg not defined", __func__, NULL);
+    pixGetDimensions(pixg, &w, &h, &d);
+    if (d != 8)
+        return (PIX *)ERROR_PTR("pixg not 8 bpp", __func__, NULL);
+
+    if (norm <= 0.0) {
+        pixGetExtremeValue(pixg, 1, L_SELECT_MAX, NULL, NULL, NULL, &maxgray);
+        norm = (maxgray > 0) ? 1.0f / (l_float32)maxgray : 1.0f;
+    }
+
+    if ((pixd = pixCreateTemplate(pixs)) == NULL)
+        return (PIX *)ERROR_PTR("pixd not made", __func__, NULL);
+    datas = pixGetData(pixs);
+    datag = pixGetData(pixg);
+    datad = pixGetData(pixd);
+    wpls = pixGetWpl(pixs);
+    wplg = pixGetWpl(pixg);
+    wpld = pixGetWpl(pixd);
+    w = L_MIN(ws, w);
+    h = L_MIN(hs, h);
+    for (i = 0; i < h; i++) {
+        lines = datas + i * wpls;
+        lineg = datag + i * wplg;
+        lined = datad + i * wpld;
+        if (ds == 8) {
+            for (j = 0; j < w; j++) {
+                vals = GET_DATA_BYTE(lines, j);
+                valg = GET_DATA_BYTE(lineg, j);
+                val = (l_int32)(vals * valg * norm + 0.5);
+                val = L_MIN(255, val);
+                SET_DATA_BYTE(lined, j, val);
+            }
+        } else {  /* ds == 32 */
+            for (j = 0; j < w; j++) {
+                val32 = *(lines + j);
+                extractRGBValues(val32, &rval, &gval, &bval);
+                valg = GET_DATA_BYTE(lineg, j);
+                rval2 = (l_int32)(rval * valg * norm + 0.5);
+                rval2 = L_MIN(255, rval2);
+                gval2 = (l_int32)(gval * valg * norm + 0.5);
+                gval2 = L_MIN(255, gval2);
+                bval2 = (l_int32)(bval * valg * norm + 0.5);
+                bval2 = L_MIN(255, bval2);
+                composeRGBPixel(rval2, gval2, bval2, lined + j);
+            }
+        }
+    }
+
+    return pixd;
+}
+
+
+/*-------------------------------------------------------------*
+ *                Grayscale threshold operation                *
+ *-------------------------------------------------------------*/
+/*!
+ * \brief   pixThresholdToValue()
+ *
+ * \param[in]    pixd       [optional]; if not null, must be equal to pixs
+ * \param[in]    pixs       8, 16, 32 bpp
+ * \param[in]    threshval
+ * \param[in]    setval
+ * \return  pixd always
+ *
+ * <pre>
+ * Notes:
+ *    ~ operation can be in-place (pixs == pixd) or to a new pixd
+ *    ~ if %setval > %threshval, sets pixels with a value >= threshval to setval
+ *    ~ if %setval < %threshval, sets pixels with a value <= threshval to setval
+ *    ~ if %setval == %threshval, no-op
+ * </pre>
+ */
+PIX *
+pixThresholdToValue(PIX      *pixd,
+                    PIX      *pixs,
+                    l_int32   threshval,
+                    l_int32   setval)
+{
+l_int32    i, j, w, h, d, wpld, setabove;
+l_uint32  *datad, *lined;
+
+    if (!pixs)
+        return (PIX *)ERROR_PTR("pixs not defined", __func__, pixd);
+    d = pixGetDepth(pixs);
+    if (d != 8 && d != 16 && d != 32)
+        return (PIX *)ERROR_PTR("pixs not 8, 16 or 32 bpp", __func__, pixd);
+    if (pixd && (pixs != pixd))
+        return (PIX *)ERROR_PTR("pixd exists and is not pixs", __func__, pixd);
+    if (threshval < 0 || setval < 0)
+        return (PIX *)ERROR_PTR("threshval & setval not < 0", __func__, pixd);
+    if (d == 8 && setval > 255)
+        return (PIX *)ERROR_PTR("setval > 255 for 8 bpp", __func__, pixd);
+    if (d == 16 && setval > 0xffff)
+        return (PIX *)ERROR_PTR("setval > 0xffff for 16 bpp", __func__, pixd);
+
+    if (!pixd)
+        pixd = pixCopy(NULL, pixs);
+    if (setval == threshval) {
+        L_WARNING("setval == threshval; no operation\n", __func__);
+        return pixd;
+    }
+
+    datad = pixGetData(pixd);
+    pixGetDimensions(pixd, &w, &h, NULL);
+    wpld = pixGetWpl(pixd);
+    if (setval > threshval)
+        setabove = TRUE;
+    else
+        setabove = FALSE;
+
+    for (i = 0; i < h; i++) {
+        lined = datad + i * wpld;
+        if (setabove == TRUE) {
+            if (d == 8) {
+                for (j = 0; j < w; j++) {
+                    if (GET_DATA_BYTE(lined, j) - threshval >= 0)
+                        SET_DATA_BYTE(lined, j, setval);
+                }
+            } else if (d == 16) {
+                for (j = 0; j < w; j++) {
+                    if (GET_DATA_TWO_BYTES(lined, j) - threshval >= 0)
+                        SET_DATA_TWO_BYTES(lined, j, setval);
+                }
+            } else {  /* d == 32 */
+                for (j = 0; j < w; j++) {
+                    if (*(lined + j) >= threshval)
+                        *(lined + j) = setval;
+                }
+            }
+        } else { /* set if below or at threshold */
+            if (d == 8) {
+                for (j = 0; j < w; j++) {
+                    if (GET_DATA_BYTE(lined, j) - threshval <= 0)
+                        SET_DATA_BYTE(lined, j, setval);
+                }
+            } else if (d == 16) {
+                for (j = 0; j < w; j++) {
+                    if (GET_DATA_TWO_BYTES(lined, j) - threshval <= 0)
+                        SET_DATA_TWO_BYTES(lined, j, setval);
+                }
+            } else {  /* d == 32 */
+                for (j = 0; j < w; j++) {
+                    if (*(lined + j) <= threshval)
+                        *(lined + j) = setval;
+                }
+            }
+        }
+    }
+
+    return pixd;
+}
+
+
+/*-------------------------------------------------------------*
+ *            Image accumulator arithmetic operations          *
+ *-------------------------------------------------------------*/
+/*!
+ * \brief   pixInitAccumulate()
+ *
+ * \param[in]    w, h      of accumulate array
+ * \param[in]    offset    initialize the 32 bpp to have this
+ *                         value; not more than 0x40000000
+ * \return  pixd   32 bpp, or NULL on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) %offset must be >= 0.
+ *      (2) %offset is used so that we can do arithmetic
+ *          with negative number results on l_uint32 data; it
+ *          prevents the l_uint32 data from going negative.
+ *      (3) Because we use l_int32 intermediate data results,
+ *          these should never exceed the max of l_int32 (0x7fffffff).
+ *          We do not permit the offset to be above 0x40000000,
+ *          which is half way between 0 and the max of l_int32.
+ *      (4) The same offset should be used for initialization,
+ *          multiplication by a constant, and final extraction!
+ *      (5) If you're only adding positive values, %offset can be 0.
+ * </pre>
+ */
+PIX *
+pixInitAccumulate(l_int32   w,
+                  l_int32   h,
+                  l_uint32  offset)
+{
+PIX  *pixd;
+
+    if ((pixd = pixCreate(w, h, 32)) == NULL)
+        return (PIX *)ERROR_PTR("pixd not made", __func__, NULL);
+    if (offset > 0x40000000)
+        offset = 0x40000000;
+    pixSetAllArbitrary(pixd, offset);
+    return pixd;
+}
+
+
+/*!
+ * \brief   pixFinalAccumulate()
+ *
+ * \param[in]    pixs     32 bpp
+ * \param[in]    offset   same as used for initialization
+ * \param[in]    depth    8, 16 or 32 bpp, of destination
+ * \return  pixd   8, 16 or 32 bpp, or NULL on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) %offset must be >= 0 and should not exceed 0x40000000.
+ *      (2) %offset is subtracted from the src 32 bpp image
+ *      (3) For 8 bpp dest, the result is clipped to [0, 0xff]
+ *      (4) For 16 bpp dest, the result is clipped to [0, 0xffff]
+ * </pre>
+ */
+PIX *
+pixFinalAccumulate(PIX      *pixs,
+                   l_uint32  offset,
+                   l_int32   depth)
+{
+l_int32    i, j, w, h, wpls, wpld, val;
+l_uint32  *datas, *datad, *lines, *lined;
+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 (depth != 8 && depth != 16 && depth != 32)
+        return (PIX *)ERROR_PTR("dest depth not 8, 16, 32 bpp", __func__, NULL);
+    if (offset > 0x40000000)
+        offset = 0x40000000;
+
+    pixGetDimensions(pixs, &w, &h, NULL);
+    if ((pixd = pixCreate(w, h, depth)) == NULL)
+        return (PIX *)ERROR_PTR("pixd not made", __func__, NULL);
+    pixCopyResolution(pixd, pixs);  /* but how did pixs get it initially? */
+    datas = pixGetData(pixs);
+    datad = pixGetData(pixd);
+    wpls = pixGetWpl(pixs);
+    wpld = pixGetWpl(pixd);
+    if (depth == 8) {
+        for (i = 0; i < h; i++) {
+            lines = datas + i * wpls;
+            lined = datad + i * wpld;
+            for (j = 0; j < w; j++) {
+                val = lines[j] - offset;
+                val = L_MAX(0, val);
+                val = L_MIN(255, val);
+                SET_DATA_BYTE(lined, j, (l_uint8)val);
+            }
+        }
+    } else if (depth == 16) {
+        for (i = 0; i < h; i++) {
+            lines = datas + i * wpls;
+            lined = datad + i * wpld;
+            for (j = 0; j < w; j++) {
+                val = lines[j] - offset;
+                val = L_MAX(0, val);
+                val = L_MIN(0xffff, val);
+                SET_DATA_TWO_BYTES(lined, j, (l_uint16)val);
+            }
+        }
+    } else {  /* depth == 32 */
+        for (i = 0; i < h; i++) {
+            lines = datas + i * wpls;
+            lined = datad + i * wpld;
+            for (j = 0; j < w; j++)
+                lined[j] = lines[j] - offset;
+        }
+    }
+
+    return pixd;
+}
+
+
+/*!
+ * \brief   pixFinalAccumulateThreshold()
+ *
+ * \param[in]    pixs        32 bpp
+ * \param[in]    offset      same as used for initialization
+ * \param[in]    threshold   values less than this are set in the destination
+ * \return  pixd   1 bpp, or NULL on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) %offset must be >= 0 and should not exceed 0x40000000.
+ *      (2) %offset is subtracted from the src 32 bpp image
+ * </pre>
+ */
+PIX *
+pixFinalAccumulateThreshold(PIX      *pixs,
+                            l_uint32  offset,
+                            l_uint32  threshold)
+{
+l_int32    i, j, w, h, wpls, wpld, val;
+l_uint32  *datas, *datad, *lines, *lined;
+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 (offset > 0x40000000)
+        offset = 0x40000000;
+
+    pixGetDimensions(pixs, &w, &h, NULL);
+    if ((pixd = pixCreate(w, h, 1)) == NULL)
+        return (PIX *)ERROR_PTR("pixd not made", __func__, NULL);
+    pixCopyResolution(pixd, pixs);  /* but how did pixs get it initially? */
+    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++) {
+            val = lines[j] - offset;
+            if (val >= threshold) {
+                SET_DATA_BIT(lined, j);
+            }
+        }
+    }
+
+    return pixd;
+}
+
+
+/*!
+ * \brief   pixAccumulate()
+ *
+ * \param[in]    pixd    32 bpp
+ * \param[in]    pixs    1, 8, 16 or 32 bpp
+ * \param[in]    op      L_ARITH_ADD or L_ARITH_SUBTRACT
+ * \return  0 if OK; 1 on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) This adds or subtracts each pixs value from pixd.
+ *      (2) This clips to the minimum of pixs and pixd, so they
+ *          do not need to be the same size.
+ *      (3) The alignment is to the origin [UL corner] of pixs & pixd.
+ * </pre>
+ */
+l_ok
+pixAccumulate(PIX     *pixd,
+              PIX     *pixs,
+              l_int32  op)
+{
+l_int32    i, j, w, h, d, wd, hd, wpls, wpld;
+l_uint32  *datas, *datad, *lines, *lined;
+
+
+    if (!pixd || (pixGetDepth(pixd) != 32))
+        return ERROR_INT("pixd not defined or not 32 bpp", __func__, 1);
+    if (!pixs)
+        return ERROR_INT("pixs not defined", __func__, 1);
+    d = pixGetDepth(pixs);
+    if (d != 1 && d != 8 && d != 16 && d != 32)
+        return ERROR_INT("pixs not 1, 8, 16 or 32 bpp", __func__, 1);
+    if (op != L_ARITH_ADD && op != L_ARITH_SUBTRACT)
+        return ERROR_INT("op must be in {L_ARITH_ADD, L_ARITH_SUBTRACT}",
+                         __func__, 1);
+
+    datas = pixGetData(pixs);
+    datad = pixGetData(pixd);
+    wpls = pixGetWpl(pixs);
+    wpld = pixGetWpl(pixd);
+    pixGetDimensions(pixs, &w, &h, NULL);
+    pixGetDimensions(pixd, &wd, &hd, NULL);
+    w = L_MIN(w, wd);
+    h = L_MIN(h, hd);
+    if (d == 1) {
+        for (i = 0; i < h; i++) {
+            lines = datas + i * wpls;
+            lined = datad + i * wpld;
+            if (op == L_ARITH_ADD) {
+                for (j = 0; j < w; j++)
+                    lined[j] += GET_DATA_BIT(lines, j);
+            } else {  /* op == L_ARITH_SUBTRACT */
+                for (j = 0; j < w; j++)
+                    lined[j] -= GET_DATA_BIT(lines, j);
+            }
+        }
+    } else if (d == 8) {
+        for (i = 0; i < h; i++) {
+            lines = datas + i * wpls;
+            lined = datad + i * wpld;
+            if (op == L_ARITH_ADD) {
+                for (j = 0; j < w; j++)
+                    lined[j] += GET_DATA_BYTE(lines, j);
+            } else {  /* op == L_ARITH_SUBTRACT */
+                for (j = 0; j < w; j++)
+                    lined[j] -= GET_DATA_BYTE(lines, j);
+            }
+        }
+    } else if (d == 16) {
+        for (i = 0; i < h; i++) {
+            lines = datas + i * wpls;
+            lined = datad + i * wpld;
+            if (op == L_ARITH_ADD) {
+                for (j = 0; j < w; j++)
+                    lined[j] += GET_DATA_TWO_BYTES(lines, j);
+            } else {  /* op == L_ARITH_SUBTRACT */
+                for (j = 0; j < w; j++)
+                    lined[j] -= GET_DATA_TWO_BYTES(lines, j);
+            }
+        }
+    } else {  /* d == 32 */
+        for (i = 0; i < h; i++) {
+            lines = datas + i * wpls;
+            lined = datad + i * wpld;
+            if (op == L_ARITH_ADD) {
+                for (j = 0; j < w; j++)
+                    lined[j] += lines[j];
+            } else {  /* op == L_ARITH_SUBTRACT */
+                for (j = 0; j < w; j++)
+                    lined[j] -= lines[j];
+            }
+        }
+    }
+
+    return 0;
+}
+
+
+/*!
+ * \brief   pixMultConstAccumulate()
+ *
+ * \param[in]    pixs      32 bpp
+ * \param[in]    factor
+ * \param[in]    offset    same as used for initialization
+ * \return  0 if OK; 1 on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) %offset must be >= 0 and should not exceed 0x40000000.
+ *      (2) This multiplies each pixel, relative to offset, by %factor.
+ *      (3) The result is returned with %offset back in place.
+ * </pre>
+ */
+l_ok
+pixMultConstAccumulate(PIX       *pixs,
+                       l_float32  factor,
+                       l_uint32   offset)
+{
+l_int32    i, j, w, h, wpl, val;
+l_uint32  *data, *line;
+
+    if (!pixs)
+        return ERROR_INT("pixs not defined", __func__, 1);
+    if (pixGetDepth(pixs) != 32)
+        return ERROR_INT("pixs not 32 bpp", __func__, 1);
+    if (offset > 0x40000000)
+        offset = 0x40000000;
+
+    pixGetDimensions(pixs, &w, &h, NULL);
+    data = pixGetData(pixs);
+    wpl = pixGetWpl(pixs);
+    for (i = 0; i < h; i++) {
+        line = data + i * wpl;
+        for (j = 0; j < w; j++) {
+            val = line[j] - offset;
+            val = (l_int32)(val * factor);
+            val += offset;
+            line[j] = (l_uint32)val;
+        }
+    }
+
+    return 0;
+}
+
+
+/*-----------------------------------------------------------------------*
+ *                      Absolute value of difference                     *
+ *-----------------------------------------------------------------------*/
+/*!
+ * \brief   pixAbsDifference()
+ *
+ * \param[in]    pixs1, pixs2    both either 8 or 16 bpp gray, or 32 bpp RGB
+ * \return  pixd, or NULL on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) The depth of pixs1 and pixs2 must be equal.
+ *      (2) Clips computation to the min size, aligning the UL corners
+ *      (3) For 8 and 16 bpp, assumes one gray component.
+ *      (4) For 32 bpp, assumes 3 color components, and ignores the
+ *          LSB of each word (the alpha channel)
+ *      (5) Computes the absolute value of the difference between
+ *          each component value.
+ * </pre>
+ */
+PIX *
+pixAbsDifference(PIX  *pixs1,
+                 PIX  *pixs2)
+{
+l_int32    i, j, w, h, w2, h2, d, wpls1, wpls2, wpld, val1, val2, diff;
+l_int32    rval1, gval1, bval1, rval2, gval2, bval2, rdiff, gdiff, bdiff;
+l_uint32  *datas1, *datas2, *datad, *lines1, *lines2, *lined;
+PIX       *pixd;
+
+    if (!pixs1)
+        return (PIX *)ERROR_PTR("pixs1 not defined", __func__, NULL);
+    if (!pixs2)
+        return (PIX *)ERROR_PTR("pixs2 not defined", __func__, NULL);
+    d = pixGetDepth(pixs1);
+    if (d != pixGetDepth(pixs2))
+        return (PIX *)ERROR_PTR("src1 and src2 depths unequal", __func__, NULL);
+    if (d != 8 && d != 16 && d != 32)
+        return (PIX *)ERROR_PTR("depths not in {8, 16, 32}", __func__, NULL);
+
+    pixGetDimensions(pixs1, &w, &h, NULL);
+    pixGetDimensions(pixs2, &w2, &h2, NULL);
+    w = L_MIN(w, w2);
+    h = L_MIN(h, h2);
+    if ((pixd = pixCreate(w, h, d)) == NULL)
+        return (PIX *)ERROR_PTR("pixd not made", __func__, NULL);
+    pixCopyResolution(pixd, pixs1);
+    datas1 = pixGetData(pixs1);
+    datas2 = pixGetData(pixs2);
+    datad = pixGetData(pixd);
+    wpls1 = pixGetWpl(pixs1);
+    wpls2 = pixGetWpl(pixs2);
+    wpld = pixGetWpl(pixd);
+    if (d == 8) {
+        for (i = 0; i < h; i++) {
+            lines1 = datas1 + i * wpls1;
+            lines2 = datas2 + i * wpls2;
+            lined = datad + i * wpld;
+            for (j = 0; j < w; j++) {
+                val1 = GET_DATA_BYTE(lines1, j);
+                val2 = GET_DATA_BYTE(lines2, j);
+                diff = L_ABS(val1 - val2);
+                SET_DATA_BYTE(lined, j, diff);
+            }
+        }
+    } else if (d == 16) {
+        for (i = 0; i < h; i++) {
+            lines1 = datas1 + i * wpls1;
+            lines2 = datas2 + i * wpls2;
+            lined = datad + i * wpld;
+            for (j = 0; j < w; j++) {
+                val1 = GET_DATA_TWO_BYTES(lines1, j);
+                val2 = GET_DATA_TWO_BYTES(lines2, j);
+                diff = L_ABS(val1 - val2);
+                SET_DATA_TWO_BYTES(lined, j, diff);
+            }
+        }
+    } else {  /* d == 32 */
+        for (i = 0; i < h; i++) {
+            lines1 = datas1 + i * wpls1;
+            lines2 = datas2 + i * wpls2;
+            lined = datad + i * wpld;
+            for (j = 0; j < w; j++) {
+                extractRGBValues(lines1[j], &rval1, &gval1, &bval1);
+                extractRGBValues(lines2[j], &rval2, &gval2, &bval2);
+                rdiff = L_ABS(rval1 - rval2);
+                gdiff = L_ABS(gval1 - gval2);
+                bdiff = L_ABS(bval1 - bval2);
+                composeRGBPixel(rdiff, gdiff, bdiff, lined + j);
+            }
+        }
+    }
+
+    return pixd;
+}
+
+
+/*-----------------------------------------------------------------------*
+ *                           Sum of color images                         *
+ *-----------------------------------------------------------------------*/
+/*!
+ * \brief   pixAddRGB()
+ *
+ * \param[in]    pixs1, pixs2    32 bpp RGB, or colormapped
+ * \return  pixd, or NULL on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) Clips computation to the minimum size, aligning the UL corners.
+ *      (2) Removes any colormap to RGB, and ignores the LSB of each
+ *          pixel word (the alpha channel).
+ *      (3) Adds each component value, pixelwise, clipping to 255.
+ *      (4) This is useful to combine two images where most of the
+ *          pixels are essentially black, such as in pixPerceptualDiff().
+ * </pre>
+ */
+PIX *
+pixAddRGB(PIX  *pixs1,
+          PIX  *pixs2)
+{
+l_int32    i, j, w, h, d, w2, h2, d2, wplc1, wplc2, wpld;
+l_int32    rval1, gval1, bval1, rval2, gval2, bval2, rval, gval, bval;
+l_uint32  *datac1, *datac2, *datad, *linec1, *linec2, *lined;
+PIX       *pixc1, *pixc2, *pixd;
+
+    if (!pixs1)
+        return (PIX *)ERROR_PTR("pixs1 not defined", __func__, NULL);
+    if (!pixs2)
+        return (PIX *)ERROR_PTR("pixs2 not defined", __func__, NULL);
+    pixGetDimensions(pixs1, &w, &h, &d);
+    pixGetDimensions(pixs2, &w2, &h2, &d2);
+    if (!pixGetColormap(pixs1) && d != 32)
+        return (PIX *)ERROR_PTR("pixs1 not cmapped or rgb", __func__, NULL);
+    if (!pixGetColormap(pixs2) && d2 != 32)
+        return (PIX *)ERROR_PTR("pixs2 not cmapped or rgb", __func__, NULL);
+    if (pixGetColormap(pixs1))
+        pixc1 = pixRemoveColormap(pixs1, REMOVE_CMAP_TO_FULL_COLOR);
+    else
+        pixc1 = pixClone(pixs1);
+    if (pixGetColormap(pixs2))
+        pixc2 = pixRemoveColormap(pixs2, REMOVE_CMAP_TO_FULL_COLOR);
+    else
+        pixc2 = pixClone(pixs2);
+
+    w = L_MIN(w, w2);
+    h = L_MIN(h, h2);
+    pixd = pixCreate(w, h, 32);
+    pixCopyResolution(pixd, pixs1);
+    datac1 = pixGetData(pixc1);
+    datac2 = pixGetData(pixc2);
+    datad = pixGetData(pixd);
+    wplc1 = pixGetWpl(pixc1);
+    wplc2 = pixGetWpl(pixc2);
+    wpld = pixGetWpl(pixd);
+    for (i = 0; i < h; i++) {
+        linec1 = datac1 + i * wplc1;
+        linec2 = datac2 + i * wplc2;
+        lined = datad + i * wpld;
+        for (j = 0; j < w; j++) {
+            extractRGBValues(linec1[j], &rval1, &gval1, &bval1);
+            extractRGBValues(linec2[j], &rval2, &gval2, &bval2);
+            rval = L_MIN(255, rval1 + rval2);
+            gval = L_MIN(255, gval1 + gval2);
+            bval = L_MIN(255, bval1 + bval2);
+            composeRGBPixel(rval, gval, bval, lined + j);
+        }
+    }
+
+    pixDestroy(&pixc1);
+    pixDestroy(&pixc2);
+    return pixd;
+}
+
+
+/*-----------------------------------------------------------------------*
+ *             Two-image min and max operations (8 and 16 bpp)           *
+ *-----------------------------------------------------------------------*/
+/*!
+ * \brief   pixMinOrMax()
+ *
+ * \param[in]    pixd     [optional] destination: this can be null,
+ *                        equal to pixs1, or different from pixs1
+ * \param[in]    pixs1    can be equal to pixd
+ * \param[in]    pixs2
+ * \param[in]    type     L_CHOOSE_MIN, L_CHOOSE_MAX
+ * \return  pixd always
+ *
+ * <pre>
+ * Notes:
+ *      (1) This gives the min or max of two images, component-wise.
+ *      (2) The depth can be 8 or 16 bpp for 1 component, and 32 bpp
+ *          for a 3 component image.  For 32 bpp, ignore the LSB
+ *          of each word (the alpha channel)
+ *      (3) There are 3 cases:
+ *          ~  if pixd == null,   MinOrMax(src1, src2) --> new pixd
+ *          ~  if pixd == pixs1,  MinOrMax(src1, src2) --> src1  (in-place)
+ *          ~  if pixd != pixs1,  MinOrMax(src1, src2) --> input pixd
+ * </pre>
+ */
+PIX *
+pixMinOrMax(PIX     *pixd,
+            PIX     *pixs1,
+            PIX     *pixs2,
+            l_int32  type)
+{
+l_int32    d, ws, hs, w, h, wpls, wpld, i, j, vals, vald, val;
+l_int32    rval1, gval1, bval1, rval2, gval2, bval2, rval, gval, bval;
+l_uint32  *datas, *datad, *lines, *lined;
+
+    if (!pixs1)
+        return (PIX *)ERROR_PTR("pixs1 not defined", __func__, pixd);
+    if (!pixs2)
+        return (PIX *)ERROR_PTR("pixs2 not defined", __func__, pixd);
+    if (pixs1 == pixs2)
+        return (PIX *)ERROR_PTR("pixs1 and pixs2 must differ", __func__, pixd);
+    if (type != L_CHOOSE_MIN && type != L_CHOOSE_MAX)
+        return (PIX *)ERROR_PTR("invalid type", __func__, pixd);
+    d = pixGetDepth(pixs1);
+    if (pixGetDepth(pixs2) != d)
+        return (PIX *)ERROR_PTR("depths unequal", __func__, pixd);
+    if (d != 8 && d != 16 && d != 32)
+        return (PIX *)ERROR_PTR("depth not 8, 16 or 32 bpp", __func__, pixd);
+
+    if (pixs1 != pixd)
+        pixd = pixCopy(pixd, pixs1);
+
+    pixGetDimensions(pixs2, &ws, &hs, NULL);
+    pixGetDimensions(pixd, &w, &h, NULL);
+    w = L_MIN(w, ws);
+    h = L_MIN(h, hs);
+    datas = pixGetData(pixs2);
+    datad = pixGetData(pixd);
+    wpls = pixGetWpl(pixs2);
+    wpld = pixGetWpl(pixd);
+    for (i = 0; i < h; i++) {
+        lines = datas + i * wpls;
+        lined = datad + i * wpld;
+        if (d == 8) {
+            for (j = 0; j < w; j++) {
+                vals = GET_DATA_BYTE(lines, j);
+                vald = GET_DATA_BYTE(lined, j);
+                if (type == L_CHOOSE_MIN)
+                    val = L_MIN(vals, vald);
+                else  /* type == L_CHOOSE_MAX */
+                    val = L_MAX(vals, vald);
+                SET_DATA_BYTE(lined, j, val);
+            }
+        } else if (d == 16) {
+            for (j = 0; j < w; j++) {
+                vals = GET_DATA_TWO_BYTES(lines, j);
+                vald = GET_DATA_TWO_BYTES(lined, j);
+                if (type == L_CHOOSE_MIN)
+                    val = L_MIN(vals, vald);
+                else  /* type == L_CHOOSE_MAX */
+                    val = L_MAX(vals, vald);
+                SET_DATA_TWO_BYTES(lined, j, val);
+            }
+        } else {  /* d == 32 */
+            for (j = 0; j < w; j++) {
+                extractRGBValues(lines[j], &rval1, &gval1, &bval1);
+                extractRGBValues(lined[j], &rval2, &gval2, &bval2);
+                if (type == L_CHOOSE_MIN) {
+                    rval = L_MIN(rval1, rval2);
+                    gval = L_MIN(gval1, gval2);
+                    bval = L_MIN(bval1, bval2);
+                } else {  /* type == L_CHOOSE_MAX */
+                    rval = L_MAX(rval1, rval2);
+                    gval = L_MAX(gval1, gval2);
+                    bval = L_MAX(bval1, bval2);
+                }
+                composeRGBPixel(rval, gval, bval, lined + j);
+            }
+        }
+    }
+
+    return pixd;
+}
+
+
+/*-----------------------------------------------------------------------*
+ *                    Scale for maximum dynamic range                    *
+ *-----------------------------------------------------------------------*/
+/*!
+ * \brief   pixMaxDynamicRange()
+ *
+ * \param[in]    pixs    4, 8, 16 or 32 bpp source
+ * \param[in]    type    L_LINEAR_SCALE or L_LOG_SCALE
+ * \return  pixd    8 bpp, or NULL on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) Scales pixel values to fit maximally within the dest 8 bpp pixd
+ *      (2) Assumes the source 'pixels' are a 1-component scalar.  For
+ *          a 32 bpp source, each pixel is treated as a single number --
+ *          not as a 3-component rgb pixel value.
+ *      (3) Uses a LUT for log scaling.
+ * </pre>
+ */
+PIX *
+pixMaxDynamicRange(PIX     *pixs,
+                   l_int32  type)
+{
+l_uint8     dval;
+l_int32     i, j, w, h, d, wpls, wpld, max;
+l_uint32   *datas, *datad;
+l_uint32    word, sval;
+l_uint32   *lines, *lined;
+l_float32   factor;
+l_float32  *tab;
+PIX        *pixd;
+
+    if (!pixs)
+        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
+    pixGetDimensions(pixs, &w, &h, &d);
+    if (d != 4 && d != 8 && d != 16 && d != 32)
+        return (PIX *)ERROR_PTR("pixs not in {4,8,16,32} bpp", __func__, NULL);
+    if (type != L_LINEAR_SCALE && type != L_LOG_SCALE)
+        return (PIX *)ERROR_PTR("invalid type", __func__, NULL);
+
+    if ((pixd = pixCreate(w, h, 8)) == NULL)
+        return (PIX *)ERROR_PTR("pixd not made", __func__, NULL);
+    pixCopyResolution(pixd, pixs);
+    datas = pixGetData(pixs);
+    datad = pixGetData(pixd);
+    wpls = pixGetWpl(pixs);
+    wpld = pixGetWpl(pixd);
+
+        /* Get max */
+    max = 0;
+    for (i = 0; i < h; i++) {
+        lines = datas + i * wpls;
+        for (j = 0; j < wpls; j++) {
+            word = *(lines + j);
+            if (d == 4) {
+                max = L_MAX(max, word >> 28);
+                max = L_MAX(max, (word >> 24) & 0xf);
+                max = L_MAX(max, (word >> 20) & 0xf);
+                max = L_MAX(max, (word >> 16) & 0xf);
+                max = L_MAX(max, (word >> 12) & 0xf);
+                max = L_MAX(max, (word >> 8) & 0xf);
+                max = L_MAX(max, (word >> 4) & 0xf);
+                max = L_MAX(max, word & 0xf);
+            } else if (d == 8) {
+                max = L_MAX(max, word >> 24);
+                max = L_MAX(max, (word >> 16) & 0xff);
+                max = L_MAX(max, (word >> 8) & 0xff);
+                max = L_MAX(max, word & 0xff);
+            } else if (d == 16) {
+                max = L_MAX(max, word >> 16);
+                max = L_MAX(max, word & 0xffff);
+            } else {  /* d == 32 (rgb) */
+                max = L_MAX(max, word);
+            }
+        }
+    }
+
+        /* Map to the full dynamic range */
+    if (d == 4) {
+        if (type == L_LINEAR_SCALE) {
+            factor = 255.f / (l_float32)max;
+            for (i = 0; i < h; i++) {
+                lines = datas + i * wpls;
+                lined = datad + i * wpld;
+                for (j = 0; j < w; j++) {
+                    sval = GET_DATA_QBIT(lines, j);
+                    dval = (l_uint8)(factor * (l_float32)sval + 0.5);
+                    SET_DATA_QBIT(lined, j, dval);
+                }
+            }
+        } else {  /* type == L_LOG_SCALE) */
+            tab = makeLogBase2Tab();
+            factor = 255.f / getLogBase2(max, tab);
+            for (i = 0; i < h; i++) {
+                lines = datas + i * wpls;
+                lined = datad + i * wpld;
+                for (j = 0; j < w; j++) {
+                    sval = GET_DATA_QBIT(lines, j);
+                    dval = (l_uint8)(factor * getLogBase2(sval, tab) + 0.5);
+                    SET_DATA_BYTE(lined, j, dval);
+                }
+            }
+            LEPT_FREE(tab);
+        }
+    } else if (d == 8) {
+        if (type == L_LINEAR_SCALE) {
+            factor = 255.f / (l_float32)max;
+            for (i = 0; i < h; i++) {
+                lines = datas + i * wpls;
+                lined = datad + i * wpld;
+                for (j = 0; j < w; j++) {
+                    sval = GET_DATA_BYTE(lines, j);
+                    dval = (l_uint8)(factor * (l_float32)sval + 0.5);
+                    SET_DATA_BYTE(lined, j, dval);
+                }
+            }
+        } else {  /* type == L_LOG_SCALE) */
+            tab = makeLogBase2Tab();
+            factor = 255.f / getLogBase2(max, tab);
+            for (i = 0; i < h; i++) {
+                lines = datas + i * wpls;
+                lined = datad + i * wpld;
+                for (j = 0; j < w; j++) {
+                    sval = GET_DATA_BYTE(lines, j);
+                    dval = (l_uint8)(factor * getLogBase2(sval, tab) + 0.5);
+                    SET_DATA_BYTE(lined, j, dval);
+                }
+            }
+            LEPT_FREE(tab);
+        }
+    } else if (d == 16) {
+        if (type == L_LINEAR_SCALE) {
+            factor = 255.f / (l_float32)max;
+            for (i = 0; i < h; i++) {
+                lines = datas + i * wpls;
+                lined = datad + i * wpld;
+                for (j = 0; j < w; j++) {
+                    sval = GET_DATA_TWO_BYTES(lines, j);
+                    dval = (l_uint8)(factor * (l_float32)sval + 0.5);
+                    SET_DATA_BYTE(lined, j, dval);
+                }
+            }
+        } else {  /* type == L_LOG_SCALE) */
+            tab = makeLogBase2Tab();
+            factor = 255.f / getLogBase2(max, tab);
+            for (i = 0; i < h; i++) {
+                lines = datas + i * wpls;
+                lined = datad + i * wpld;
+                for (j = 0; j < w; j++) {
+                    sval = GET_DATA_TWO_BYTES(lines, j);
+                    dval = (l_uint8)(factor * getLogBase2(sval, tab) + 0.5);
+                    SET_DATA_BYTE(lined, j, dval);
+                }
+            }
+            LEPT_FREE(tab);
+        }
+    } else {  /* d == 32 */
+        if (type == L_LINEAR_SCALE) {
+            factor = 255.f / (l_float32)max;
+            for (i = 0; i < h; i++) {
+                lines = datas + i * wpls;
+                lined = datad + i * wpld;
+                for (j = 0; j < w; j++) {
+                    sval = lines[j];
+                    dval = (l_uint8)(factor * (l_float32)sval + 0.5);
+                    SET_DATA_BYTE(lined, j, dval);
+                }
+            }
+        } else {  /* type == L_LOG_SCALE) */
+            tab = makeLogBase2Tab();
+            factor = 255.f / getLogBase2(max, tab);
+            for (i = 0; i < h; i++) {
+                lines = datas + i * wpls;
+                lined = datad + i * wpld;
+                for (j = 0; j < w; j++) {
+                    sval = lines[j];
+                    dval = (l_uint8)(factor * getLogBase2(sval, tab) + 0.5);
+                    SET_DATA_BYTE(lined, j, dval);
+                }
+            }
+            LEPT_FREE(tab);
+        }
+    }
+
+    return pixd;
+}
+
+
+/*!
+ * \brief   pixMaxDynamicRangeRGB()
+ *
+ * \param[in]    pixs    32 bpp rgb source
+ * \param[in]    type    L_LINEAR_SCALE or L_LOG_SCALE
+ * \return  pixd   32 bpp, or NULL on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) Scales pixel values to fit maximally within a 32 bpp dest pixd
+ *      (2) All color components are scaled with the same factor, based
+ *          on the maximum r, g or b component in the image.  This should
+ *          not be used if the 32-bit value is a single number (e.g., a
+ *          count in a histogram generated by pixMakeHistoHS()).
+ *      (3) Uses a LUT for log scaling.
+ * </pre>
+ */
+PIX *
+pixMaxDynamicRangeRGB(PIX     *pixs,
+                      l_int32  type)
+{
+l_int32     i, j, w, h, wpls, wpld, max;
+l_uint32    sval, dval, word;
+l_uint32   *datas, *datad;
+l_uint32   *lines, *lined;
+l_float32   factor;
+l_float32  *tab;
+PIX        *pixd;
+
+    if (!pixs || pixGetDepth(pixs) != 32)
+        return (PIX *)ERROR_PTR("pixs undefined or not 32 bpp", __func__, NULL);
+    if (type != L_LINEAR_SCALE && type != L_LOG_SCALE)
+        return (PIX *)ERROR_PTR("invalid type", __func__, NULL);
+
+        /* Get max */
+    pixd = pixCreateTemplate(pixs);
+    datas = pixGetData(pixs);
+    datad = pixGetData(pixd);
+    wpls = pixGetWpl(pixs);
+    wpld = pixGetWpl(pixd);
+    pixGetDimensions(pixs, &w, &h, NULL);
+    max = 0;
+    for (i = 0; i < h; i++) {
+        lines = datas + i * wpls;
+        for (j = 0; j < wpls; j++) {
+            word = lines[j];
+            max = L_MAX(max, word >> 24);
+            max = L_MAX(max, (word >> 16) & 0xff);
+            max = L_MAX(max, (word >> 8) & 0xff);
+        }
+    }
+    if (max == 0) {
+        L_WARNING("max = 0; setting to 1\n", __func__);
+        max = 1;
+    }
+
+        /* Map to the full dynamic range */
+    if (type == L_LINEAR_SCALE) {
+        factor = 255.f / (l_float32)max;
+        for (i = 0; i < h; i++) {
+            lines = datas + i * wpls;
+            lined = datad + i * wpld;
+            for (j = 0; j < w; j++) {
+                sval = lines[j];
+                dval = linearScaleRGBVal(sval, factor);
+                lined[j] = dval;
+            }
+        }
+    } else {  /* type == L_LOG_SCALE) */
+        tab = makeLogBase2Tab();
+        factor = 255.f / getLogBase2(max, tab);
+        for (i = 0; i < h; i++) {
+            lines = datas + i * wpls;
+            lined = datad + i * wpld;
+            for (j = 0; j < w; j++) {
+                sval = lines[j];
+                dval = logScaleRGBVal(sval, tab, factor);
+                lined[j] = dval;
+            }
+        }
+        LEPT_FREE(tab);
+    }
+
+    return pixd;
+}
+
+
+/*-----------------------------------------------------------------------*
+ *                         RGB pixel value scaling                       *
+ *-----------------------------------------------------------------------*/
+/*!
+ * \brief   linearScaleRGBVal()
+ *
+ * \param[in]    sval     32-bit rgb pixel value
+ * \param[in]    factor   multiplication factor on each component
+ * \return  dval  linearly scaled version of %sval
+ *
+ * <pre>
+ * Notes:
+ *      (1) %factor must be chosen to be not greater than (255 / maxcomp),
+ *          where maxcomp is the maximum value of the pixel components.
+ *          Otherwise, the product will overflow a uint8.  In use, factor
+ *          is the same for all pixels in a pix.
+ *      (2) No scaling is performed on the transparency ("A") component.
+ * </pre>
+ */
+l_uint32
+linearScaleRGBVal(l_uint32   sval,
+                  l_float32  factor)
+{
+l_uint32  dval;
+
+    dval = ((l_uint8)(factor * (sval >> 24) + 0.5f) << 24) |
+           ((l_uint8)(factor * ((sval >> 16) & 0xff) + 0.5f) << 16) |
+           ((l_uint8)(factor * ((sval >> 8) & 0xff) + 0.5f) << 8) |
+           (sval & 0xff);
+    return dval;
+}
+
+
+/*!
+ * \brief   logScaleRGBVal()
+ *
+ * \param[in]    sval     32-bit rgb pixel value
+ * \param[in]    tab      256 entry log-base-2 table
+ * \param[in]    factor   multiplication factor on each component
+ * \return  dval  log scaled version of %sval
+ *
+ * <pre>
+ * Notes:
+ *      (1) %tab is made with makeLogBase2Tab().
+ *      (2) %factor must be chosen to be not greater than
+ *          255.0 / log[base2](maxcomp), where maxcomp is the maximum
+ *          value of the pixel components.  Otherwise, the product
+ *          will overflow a uint8.  In use, factor is the same for
+ *          all pixels in a pix.
+ *      (3) No scaling is performed on the transparency ("A") component.
+ * </pre>
+ */
+l_uint32
+logScaleRGBVal(l_uint32    sval,
+               l_float32  *tab,
+               l_float32   factor)
+{
+l_uint32  dval;
+
+    dval = ((l_uint8)(factor * getLogBase2(sval >> 24, tab) + 0.5f) << 24) |
+           ((l_uint8)(factor * getLogBase2(((sval >> 16) & 0xff), tab) + 0.5f)
+                     << 16) |
+           ((l_uint8)(factor * getLogBase2(((sval >> 8) & 0xff), tab) + 0.5f)
+                     << 8) |
+           (sval & 0xff);
+    return dval;
+}
+
+
+/*-----------------------------------------------------------------------*
+ *                            Log base2 lookup                           *
+ *-----------------------------------------------------------------------*/
+/*
+ * \brief   makeLogBase2Tab()
+ *
+ * \return  tab   table giving the log[base2] of values from 1 to 255
+ */
+l_float32 *
+makeLogBase2Tab(void)
+{
+l_int32     i;
+l_float32   log2;
+l_float32  *tab;
+
+    if ((tab = (l_float32 *)LEPT_CALLOC(256, sizeof(l_float32))) == NULL)
+        return (l_float32 *)ERROR_PTR("tab not made", __func__, NULL);
+
+    log2 = (l_float32)log((l_float32)2);
+    for (i = 0; i < 256; i++)
+        tab[i] = (l_float32)log((l_float32)i) / log2;
+
+    return tab;
+}
+
+
+/*
+ * \brief   getLogBase2()
+ *
+ * \param[in]    val       in range [0 ... 255]
+ * \param[in]    logtab    256-entry table of logs
+ * \return       logval    log[base2] of %val, or 0 on error
+ */
+l_float32
+getLogBase2(l_int32     val,
+            l_float32  *logtab)
+{
+    if (!logtab)
+        return ERROR_INT("logtab not defined", __func__, 0);
+
+    if (val < 0x100)
+        return logtab[val];
+    else if (val < 0x10000)
+        return 8.0f + logtab[val >> 8];
+    else if (val < 0x1000000)
+        return 16.0f + logtab[val >> 16];
+    else
+        return 24.0f + logtab[val >> 24];
+}