--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mupdf-source/thirdparty/leptonica/src/kernel.c	Mon Sep 15 11:43:07 2025 +0200
@@ -0,0 +1,1239 @@
+/*====================================================================*
+ -  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 kernel.c
+ * <pre>
+ *
+ *      Basic operations on kernels for image convolution
+ *
+ *         Create/destroy/copy
+ *            L_KERNEL   *kernelCreate()
+ *            void        kernelDestroy()
+ *            L_KERNEL   *kernelCopy()
+ *
+ *         Accessors:
+ *            l_int32     kernelGetElement()
+ *            l_int32     kernelSetElement()
+ *            l_int32     kernelGetParameters()
+ *            l_int32     kernelSetOrigin()
+ *            l_int32     kernelGetSum()
+ *            l_int32     kernelGetMinMax()
+ *
+ *         Normalize/invert
+ *            L_KERNEL   *kernelNormalize()
+ *            L_KERNEL   *kernelInvert()
+ *
+ *         Helper function
+ *            l_float32 **create2dFloatArray()
+ *
+ *         Serialized I/O
+ *            L_KERNEL   *kernelRead()
+ *            L_KERNEL   *kernelReadStream()
+ *            l_int32     kernelWrite()
+ *            l_int32     kernelWriteStream()
+ *
+ *         Making a kernel from a compiled string
+ *            L_KERNEL   *kernelCreateFromString()
+ *
+ *         Making a kernel from a simple file format
+ *            L_KERNEL   *kernelCreateFromFile()
+ *
+ *         Making a kernel from a Pix
+ *            L_KERNEL   *kernelCreateFromPix()
+ *
+ *         Display a kernel in a pix
+ *            PIX        *kernelDisplayInPix()
+ *
+ *         Parse string to extract numbers
+ *            NUMA       *parseStringForNumbers()
+ *
+ *      Simple parametric kernels
+ *            L_KERNEL   *makeFlatKernel()
+ *            L_KERNEL   *makeGaussianKernel()
+ *            L_KERNEL   *makeGaussianKernelSep()
+ *            L_KERNEL   *makeDoGKernel()
+ * </pre>
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config_auto.h>
+#endif  /* HAVE_CONFIG_H */
+
+#include <string.h>
+#include <math.h>
+#include "allheaders.h"
+
+    /* Array size must be > 0 and not larger than this */
+static const l_uint32  MaxArraySize = 100000;
+
+/*------------------------------------------------------------------------*
+ *                           Create / Destroy                             *
+ *------------------------------------------------------------------------*/
+/*!
+ * \brief   kernelCreate()
+ *
+ * \param[in]    height, width
+ * \return  kernel, or NULL on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) kernelCreate() initializes all values to 0.
+ *      (2) After this call, (cy,cx) and nonzero data values must be
+ *          assigned.
+ *      (2) The number of kernel elements must be less than 2^29.
+ * </pre>
+ */
+L_KERNEL *
+kernelCreate(l_int32  height,
+             l_int32  width)
+{
+l_uint64   size64;
+L_KERNEL  *kel;
+
+    if (width <= 0)
+        return (L_KERNEL *)ERROR_PTR("width must be > 0", __func__, NULL);
+    if (height <= 0)
+        return (L_KERNEL *)ERROR_PTR("height must be > 0", __func__, NULL);
+
+        /* Avoid overflow in malloc arg */
+    size64 = (l_uint64)width * (l_uint64)height;
+    if (size64 >= (1LL << 29)) {
+        L_ERROR("requested width = %d, height = %d\n", __func__, width, height);
+        return (L_KERNEL *)ERROR_PTR("size >= 2^29", __func__, NULL);
+    }
+
+    kel = (L_KERNEL *)LEPT_CALLOC(1, sizeof(L_KERNEL));
+    kel->sy = height;
+    kel->sx = width;
+    if ((kel->data = create2dFloatArray(height, width)) == NULL) {
+        LEPT_FREE(kel);
+        return (L_KERNEL *)ERROR_PTR("data not allocated", __func__, NULL);
+    }
+    return kel;
+}
+
+
+/*!
+ * \brief   kernelDestroy()
+ *
+ * \param[in,out]   pkel    will be set to null before returning
+ * \return  void
+ */
+void
+kernelDestroy(L_KERNEL  **pkel)
+{
+l_int32    i;
+L_KERNEL  *kel;
+
+    if (pkel == NULL)  {
+        L_WARNING("ptr address is NULL!\n", __func__);
+        return;
+    }
+    if ((kel = *pkel) == NULL)
+        return;
+
+    for (i = 0; i < kel->sy; i++)
+        LEPT_FREE(kel->data[i]);
+    LEPT_FREE(kel->data);
+    LEPT_FREE(kel);
+    *pkel = NULL;
+}
+
+
+/*!
+ * \brief   kernelCopy()
+ *
+ * \param[in]    kels    source kernel
+ * \return  keld   copy of kels, or NULL on error
+ */
+L_KERNEL *
+kernelCopy(L_KERNEL  *kels)
+{
+l_int32    i, j, sx, sy, cx, cy;
+L_KERNEL  *keld;
+
+    if (!kels)
+        return (L_KERNEL *)ERROR_PTR("kels not defined", __func__, NULL);
+
+    kernelGetParameters(kels, &sy, &sx, &cy, &cx);
+    if ((keld = kernelCreate(sy, sx)) == NULL)
+        return (L_KERNEL *)ERROR_PTR("keld not made", __func__, NULL);
+    keld->cy = cy;
+    keld->cx = cx;
+    for (i = 0; i < sy; i++)
+        for (j = 0; j < sx; j++)
+            keld->data[i][j] = kels->data[i][j];
+
+    return keld;
+}
+
+
+/*----------------------------------------------------------------------*
+ *                               Accessors                              *
+ *----------------------------------------------------------------------*/
+/*!
+ * \brief   kernelGetElement()
+ *
+ * \param[in]    kel
+ * \param[in]    row
+ * \param[in]    col
+ * \param[out]   pval
+ * \return  0 if OK; 1 on error
+ */
+l_ok
+kernelGetElement(L_KERNEL   *kel,
+                 l_int32     row,
+                 l_int32     col,
+                 l_float32  *pval)
+{
+    if (!pval)
+        return ERROR_INT("&val not defined", __func__, 1);
+    *pval = 0;
+    if (!kel)
+        return ERROR_INT("kernel not defined", __func__, 1);
+    if (row < 0 || row >= kel->sy)
+        return ERROR_INT("kernel row out of bounds", __func__, 1);
+    if (col < 0 || col >= kel->sx)
+        return ERROR_INT("kernel col out of bounds", __func__, 1);
+
+    *pval = kel->data[row][col];
+    return 0;
+}
+
+
+/*!
+ * \brief   kernelSetElement()
+ *
+ * \param[in]    kel kernel
+ * \param[in]    row
+ * \param[in]    col
+ * \param[in]    val
+ * \return  0 if OK; 1 on error
+ */
+l_ok
+kernelSetElement(L_KERNEL  *kel,
+                 l_int32    row,
+                 l_int32    col,
+                 l_float32  val)
+{
+    if (!kel)
+        return ERROR_INT("kel not defined", __func__, 1);
+    if (row < 0 || row >= kel->sy)
+        return ERROR_INT("kernel row out of bounds", __func__, 1);
+    if (col < 0 || col >= kel->sx)
+        return ERROR_INT("kernel col out of bounds", __func__, 1);
+
+    kel->data[row][col] = val;
+    return 0;
+}
+
+
+/*!
+ * \brief   kernelGetParameters()
+ *
+ * \param[in]    kel                  kernel
+ * \param[out]   psy, psx, pcy, pcx   [optional] each can be null
+ * \return  0 if OK, 1 on error
+ */
+l_ok
+kernelGetParameters(L_KERNEL  *kel,
+                    l_int32   *psy,
+                    l_int32   *psx,
+                    l_int32   *pcy,
+                    l_int32   *pcx)
+{
+    if (psy) *psy = 0;
+    if (psx) *psx = 0;
+    if (pcy) *pcy = 0;
+    if (pcx) *pcx = 0;
+    if (!kel)
+        return ERROR_INT("kernel not defined", __func__, 1);
+    if (psy) *psy = kel->sy;
+    if (psx) *psx = kel->sx;
+    if (pcy) *pcy = kel->cy;
+    if (pcx) *pcx = kel->cx;
+    return 0;
+}
+
+
+/*!
+ * \brief   kernelSetOrigin()
+ *
+ * \param[in]    kel       kernel
+ * \param[in]    cy, cx
+ * \return  0 if OK; 1 on error
+ */
+l_ok
+kernelSetOrigin(L_KERNEL  *kel,
+                l_int32    cy,
+                l_int32    cx)
+{
+    if (!kel)
+        return ERROR_INT("kel not defined", __func__, 1);
+    kel->cy = cy;
+    kel->cx = cx;
+    return 0;
+}
+
+
+/*!
+ * \brief   kernelGetSum()
+ *
+ * \param[in]    kel      kernel
+ * \param[out]   psum     sum of all kernel values
+ * \return  0 if OK, 1 on error
+ */
+l_ok
+kernelGetSum(L_KERNEL   *kel,
+             l_float32  *psum)
+{
+l_int32    sx, sy, i, j;
+
+    if (!psum)
+        return ERROR_INT("&sum not defined", __func__, 1);
+    *psum = 0.0;
+    if (!kel)
+        return ERROR_INT("kernel not defined", __func__, 1);
+
+    kernelGetParameters(kel, &sy, &sx, NULL, NULL);
+    for (i = 0; i < sy; i++) {
+        for (j = 0; j < sx; j++) {
+            *psum += kel->data[i][j];
+        }
+    }
+    return 0;
+}
+
+
+/*!
+ * \brief   kernelGetMinMax()
+ *
+ * \param[in]    kel      kernel
+ * \param[out]   pmin     [optional] minimum value
+ * \param[out]   pmax     [optional] maximum value
+ * \return  0 if OK, 1 on error
+ */
+l_ok
+kernelGetMinMax(L_KERNEL   *kel,
+                l_float32  *pmin,
+                l_float32  *pmax)
+{
+l_int32    sx, sy, i, j;
+l_float32  val, minval, maxval;
+
+    if (!pmin && !pmax)
+        return ERROR_INT("neither &min nor &max defined", __func__, 1);
+    if (pmin) *pmin = 0.0;
+    if (pmax) *pmax = 0.0;
+    if (!kel)
+        return ERROR_INT("kernel not defined", __func__, 1);
+
+    kernelGetParameters(kel, &sy, &sx, NULL, NULL);
+    minval = 10000000.0;
+    maxval = -10000000.0;
+    for (i = 0; i < sy; i++) {
+        for (j = 0; j < sx; j++) {
+            val = kel->data[i][j];
+            if (val < minval)
+                minval = val;
+            if (val > maxval)
+                maxval = val;
+        }
+    }
+    if (pmin)
+        *pmin = minval;
+    if (pmax)
+        *pmax = maxval;
+
+    return 0;
+}
+
+
+/*----------------------------------------------------------------------*
+ *                          Normalize/Invert                            *
+ *----------------------------------------------------------------------*/
+/*!
+ * \brief   kernelNormalize()
+ *
+ * \param[in]    kels      source kel, to be normalized
+ * \param[in]    normsum   desired sum of elements in keld
+ * \return  keld   normalized version of kels, or NULL on error
+ *                 or if sum of elements is very close to 0)
+ *
+ * <pre>
+ * Notes:
+ *      (1) If the sum of kernel elements is close to 0, do not
+ *          try to calculate the normalized kernel.  Instead,
+ *          return a copy of the input kernel, with a warning.
+ * </pre>
+ */
+L_KERNEL *
+kernelNormalize(L_KERNEL  *kels,
+                l_float32  normsum)
+{
+l_int32    i, j, sx, sy, cx, cy;
+l_float32  sum, factor;
+L_KERNEL  *keld;
+
+    if (!kels)
+        return (L_KERNEL *)ERROR_PTR("kels not defined", __func__, NULL);
+
+    kernelGetSum(kels, &sum);
+    if (L_ABS(sum) < 0.00001) {
+        L_WARNING("null sum; not normalizing; returning a copy\n", __func__);
+        return kernelCopy(kels);
+    }
+
+    kernelGetParameters(kels, &sy, &sx, &cy, &cx);
+    if ((keld = kernelCreate(sy, sx)) == NULL)
+        return (L_KERNEL *)ERROR_PTR("keld not made", __func__, NULL);
+    keld->cy = cy;
+    keld->cx = cx;
+
+    factor = normsum / sum;
+    for (i = 0; i < sy; i++)
+        for (j = 0; j < sx; j++)
+            keld->data[i][j] = factor * kels->data[i][j];
+
+    return keld;
+}
+
+
+/*!
+ * \brief   kernelInvert()
+ *
+ * \param[in]    kels   source kel, to be inverted
+ * \return  keld   spatially inverted, about the origin, or NULL on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) For convolution, the kernel is spatially inverted before
+ *          a "correlation" operation is done between the kernel and the image.
+ * </pre>
+ */
+L_KERNEL *
+kernelInvert(L_KERNEL  *kels)
+{
+l_int32    i, j, sx, sy, cx, cy;
+L_KERNEL  *keld;
+
+    if (!kels)
+        return (L_KERNEL *)ERROR_PTR("kels not defined", __func__, NULL);
+
+    kernelGetParameters(kels, &sy, &sx, &cy, &cx);
+    if ((keld = kernelCreate(sy, sx)) == NULL)
+        return (L_KERNEL *)ERROR_PTR("keld not made", __func__, NULL);
+    keld->cy = sy - 1 - cy;
+    keld->cx = sx - 1 - cx;
+
+    for (i = 0; i < sy; i++)
+        for (j = 0; j < sx; j++)
+            keld->data[i][j] = kels->data[sy - 1 - i][sx - 1 - j];
+
+    return keld;
+}
+
+
+/*----------------------------------------------------------------------*
+ *                            Helper function                           *
+ *----------------------------------------------------------------------*/
+/*!
+ * \brief   create2dFloatArray()
+ *
+ * \param[in]    sy   rows == height
+ * \param[in]    sx   columns == width
+ * \return  doubly indexed array i.e., an array of sy row pointers,
+ *              each of which points to an array of sx floats
+ *
+ * <pre>
+ * Notes:
+ *      (1) The array[%sy][%sx] is indexed in standard "matrix notation",
+ *          with the row index first.
+ *      (2) The caller kernelCreate() limits the size to < 2^29 pixels.
+ * </pre>
+ */
+l_float32 **
+create2dFloatArray(l_int32  sy,
+                   l_int32  sx)
+{
+l_int32      i;
+l_float32  **array;
+
+    if (sx <= 0 || sx > MaxArraySize)
+        return (l_float32 **)ERROR_PTR("sx out of bounds", __func__, NULL);
+    if (sy <= 0 || sy > MaxArraySize)
+        return (l_float32 **)ERROR_PTR("sy out of bounds", __func__, NULL);
+
+    array = (l_float32 **)LEPT_CALLOC(sy, sizeof(l_float32 *));
+    for (i = 0; i < sy; i++)
+        array[i] = (l_float32 *)LEPT_CALLOC(sx, sizeof(l_float32));
+    return array;
+}
+
+
+/*----------------------------------------------------------------------*
+ *                            Kernel serialized I/O                     *
+ *----------------------------------------------------------------------*/
+/*!
+ * \brief   kernelRead()
+ *
+ * \param[in]    fname    filename
+ * \return  kernel, or NULL on error
+ */
+L_KERNEL *
+kernelRead(const char  *fname)
+{
+FILE      *fp;
+L_KERNEL  *kel;
+
+    if (!fname)
+        return (L_KERNEL *)ERROR_PTR("fname not defined", __func__, NULL);
+
+    if ((fp = fopenReadStream(fname)) == NULL)
+        return (L_KERNEL *)ERROR_PTR_1("stream not opened",
+                                       fname, __func__, NULL);
+    if ((kel = kernelReadStream(fp)) == NULL) {
+        fclose(fp);
+        return (L_KERNEL *)ERROR_PTR_1("kel not returned",
+                                       fname, __func__, NULL);
+    }
+    fclose(fp);
+
+    return kel;
+}
+
+
+/*!
+ * \brief   kernelReadStream()
+ *
+ * \param[in]    fp    file stream
+ * \return  kernel, or NULL on error
+ */
+L_KERNEL *
+kernelReadStream(FILE  *fp)
+{
+l_int32    sy, sx, cy, cx, i, j, ret, version, ignore;
+L_KERNEL  *kel;
+
+    if (!fp)
+        return (L_KERNEL *)ERROR_PTR("stream not defined", __func__, NULL);
+
+    ret = fscanf(fp, "  Kernel Version %d\n", &version);
+    if (ret != 1)
+        return (L_KERNEL *)ERROR_PTR("not a kernel file", __func__, NULL);
+    if (version != KERNEL_VERSION_NUMBER)
+        return (L_KERNEL *)ERROR_PTR("invalid kernel version", __func__, NULL);
+
+    if (fscanf(fp, "  sy = %d, sx = %d, cy = %d, cx = %d\n",
+            &sy, &sx, &cy, &cx) != 4)
+        return (L_KERNEL *)ERROR_PTR("dimensions not read", __func__, NULL);
+    if (sx > MaxArraySize || sy > MaxArraySize) {
+        L_ERROR("sx = %d or sy = %d > %d\n", __func__, sx, sy, MaxArraySize);
+        return NULL;
+    }
+    if ((kel = kernelCreate(sy, sx)) == NULL)
+        return (L_KERNEL *)ERROR_PTR("kel not made", __func__, NULL);
+    kernelSetOrigin(kel, cy, cx);
+
+    for (i = 0; i < sy; i++) {
+        for (j = 0; j < sx; j++)
+            ignore = fscanf(fp, "%15f", &kel->data[i][j]);
+        ignore = fscanf(fp, "\n");
+    }
+    ignore = fscanf(fp, "\n");
+
+    return kel;
+}
+
+
+/*!
+ * \brief   kernelWrite()
+ *
+ * \param[in]    fname    output file
+ * \param[in]    kel      kernel
+ * \return  0 if OK, 1 on error
+ */
+l_ok
+kernelWrite(const char  *fname,
+            L_KERNEL    *kel)
+{
+FILE  *fp;
+
+    if (!fname)
+        return ERROR_INT("fname not defined", __func__, 1);
+    if (!kel)
+        return ERROR_INT("kel not defined", __func__, 1);
+
+    if ((fp = fopenWriteStream(fname, "wb")) == NULL)
+        return ERROR_INT_1("stream not opened", fname, __func__, 1);
+    kernelWriteStream(fp, kel);
+    fclose(fp);
+
+    return 0;
+}
+
+
+/*!
+ * \brief   kernelWriteStream()
+ *
+ * \param[in]    fp    file stream
+ * \param[in]    kel
+ * \return  0 if OK, 1 on error
+ */
+l_ok
+kernelWriteStream(FILE      *fp,
+                  L_KERNEL  *kel)
+{
+l_int32  sx, sy, cx, cy, i, j;
+
+    if (!fp)
+        return ERROR_INT("stream not defined", __func__, 1);
+    if (!kel)
+        return ERROR_INT("kel not defined", __func__, 1);
+    kernelGetParameters(kel, &sy, &sx, &cy, &cx);
+
+    fprintf(fp, "  Kernel Version %d\n", KERNEL_VERSION_NUMBER);
+    fprintf(fp, "  sy = %d, sx = %d, cy = %d, cx = %d\n", sy, sx, cy, cx);
+    for (i = 0; i < sy; i++) {
+        for (j = 0; j < sx; j++)
+            fprintf(fp, "%15.4f", kel->data[i][j]);
+        fprintf(fp, "\n");
+    }
+    fprintf(fp, "\n");
+
+    return 0;
+}
+
+
+/*----------------------------------------------------------------------*
+ *                 Making a kernel from a compiled string               *
+ *----------------------------------------------------------------------*/
+/*!
+ * \brief   kernelCreateFromString()
+ *
+ * \param[in]    h, w     height, width
+ * \param[in]    cy, cx   origin
+ * \param[in]    kdata
+ * \return  kernel of the given size, or NULL on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) The data is an array of chars, in row-major order, giving
+ *          space separated integers in the range [-255 ... 255].
+ *      (2) The only other formatting limitation is that you must
+ *          leave space between the last number in each row and
+ *          the double-quote.  If possible, it's also nice to have each
+ *          line in the string represent a line in the kernel; e.g.,
+ *              static const char *kdata =
+ *                  " 20   50   20 "
+ *                  " 70  140   70 "
+ *                  " 20   50   20 ";
+ * </pre>
+ */
+L_KERNEL *
+kernelCreateFromString(l_int32      h,
+                       l_int32      w,
+                       l_int32      cy,
+                       l_int32      cx,
+                       const char  *kdata)
+{
+l_int32    n, i, j, index;
+l_float32  val;
+L_KERNEL  *kel;
+NUMA      *na;
+
+    if (h < 1)
+        return (L_KERNEL *)ERROR_PTR("height must be > 0", __func__, NULL);
+    if (w < 1)
+        return (L_KERNEL *)ERROR_PTR("width must be > 0", __func__, NULL);
+    if (cy < 0 || cy >= h)
+        return (L_KERNEL *)ERROR_PTR("cy invalid", __func__, NULL);
+    if (cx < 0 || cx >= w)
+        return (L_KERNEL *)ERROR_PTR("cx invalid", __func__, NULL);
+
+    kel = kernelCreate(h, w);
+    kernelSetOrigin(kel, cy, cx);
+    na = parseStringForNumbers(kdata, " \t\n");
+    n = numaGetCount(na);
+    if (n != w * h) {
+        kernelDestroy(&kel);
+        numaDestroy(&na);
+        lept_stderr("w = %d, h = %d, num ints = %d\n", w, h, n);
+        return (L_KERNEL *)ERROR_PTR("invalid integer data", __func__, NULL);
+    }
+
+    index = 0;
+    for (i = 0; i < h; i++) {
+        for (j = 0; j < w; j++) {
+            numaGetFValue(na, index, &val);
+            kernelSetElement(kel, i, j, val);
+            index++;
+        }
+    }
+
+    numaDestroy(&na);
+    return kel;
+}
+
+
+/*----------------------------------------------------------------------*
+ *                Making a kernel from a simple file format             *
+ *----------------------------------------------------------------------*/
+/*!
+ * \brief   kernelCreateFromFile()
+ *
+ * \param[in]    filename
+ * \return  kernel, or NULL on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) The file contains, in the following order:
+ *           ~ Any number of comment lines starting with '#' are ignored
+ *           ~ The height and width of the kernel
+ *           ~ The y and x values of the kernel origin
+ *           ~ The kernel data, formatted as lines of numbers (integers
+ *             or floats) for the kernel values in row-major order,
+ *             and with no other punctuation.
+ *             (Note: this differs from kernelCreateFromString(),
+ *             where each line must begin and end with a double-quote
+ *             to tell the compiler it's part of a string.)
+ *           ~ The kernel specification ends when a blank line,
+ *             a comment line, or the end of file is reached.
+ *      (2) All lines must be left-justified.
+ *      (3) See kernelCreateFromString() for a description of the string
+ *          format for the kernel data.  As an example, here are the lines
+ *          of a valid kernel description file  In the file, all lines
+ *          are left-justified:
+ * \code
+ *                    # small 3x3 kernel
+ *                    3 3
+ *                    1 1
+ *                    25.5   51    24.3
+ *                    70.2  146.3  73.4
+ *                    20     50.9  18.4
+ * \endcode
+ * </pre>
+ */
+L_KERNEL *
+kernelCreateFromFile(const char  *filename)
+{
+char      *filestr, *line;
+l_int32    nlines, i, j, first, index, w, h, cx, cy, n;
+l_float32  val;
+size_t     size;
+NUMA      *na, *nat;
+SARRAY    *sa;
+L_KERNEL  *kel;
+
+    if (!filename)
+        return (L_KERNEL *)ERROR_PTR("filename not defined", __func__, NULL);
+
+    if ((filestr = (char *)l_binaryRead(filename, &size)) == NULL)
+        return (L_KERNEL *)ERROR_PTR_1("file not found",
+                                       filename, __func__, NULL);
+    if (size == 0) {
+        LEPT_FREE(filestr);
+        return (L_KERNEL *)ERROR_PTR_1("file is empty",
+                                       filename, __func__, NULL);
+    }
+
+    sa = sarrayCreateLinesFromString(filestr, 1);
+    LEPT_FREE(filestr);
+    nlines = sarrayGetCount(sa);
+
+        /* Find the first data line. */
+    for (i = 0, first = 0; i < nlines; i++) {
+        line = sarrayGetString(sa, i, L_NOCOPY);
+        if (line[0] != '#') {
+            first = i;
+            break;
+        }
+    }
+
+        /* Find the kernel dimensions and origin location. */
+    line = sarrayGetString(sa, first, L_NOCOPY);
+    if (sscanf(line, "%d %d", &h, &w) != 2) {
+        sarrayDestroy(&sa);
+        return (L_KERNEL *)ERROR_PTR("error reading h,w", __func__, NULL);
+    }
+    if (h > MaxArraySize || w > MaxArraySize) {
+        L_ERROR("h = %d or w = %d > %d\n", __func__, h, w, MaxArraySize);
+        sarrayDestroy(&sa);
+        return NULL;
+    }
+    line = sarrayGetString(sa, first + 1, L_NOCOPY);
+    if (sscanf(line, "%d %d", &cy, &cx) != 2) {
+        sarrayDestroy(&sa);
+        return (L_KERNEL *)ERROR_PTR("error reading cy,cx", __func__, NULL);
+    }
+
+        /* Extract the data.  This ends when we reach eof, or when we
+         * encounter a line of data that is either a null string or
+         * contains just a newline. */
+    na = numaCreate(0);
+    for (i = first + 2; i < nlines; i++) {
+        line = sarrayGetString(sa, i, L_NOCOPY);
+        if (line[0] == '\0' || line[0] == '\n' || line[0] == '#')
+            break;
+        nat = parseStringForNumbers(line, " \t\n");
+        numaJoin(na, nat, 0, -1);
+        numaDestroy(&nat);
+    }
+    sarrayDestroy(&sa);
+
+    n = numaGetCount(na);
+    if (n != w * h) {
+        numaDestroy(&na);
+        lept_stderr("w = %d, h = %d, num ints = %d\n", w, h, n);
+        return (L_KERNEL *)ERROR_PTR("invalid integer data", __func__, NULL);
+    }
+
+    kel = kernelCreate(h, w);
+    kernelSetOrigin(kel, cy, cx);
+    index = 0;
+    for (i = 0; i < h; i++) {
+        for (j = 0; j < w; j++) {
+            numaGetFValue(na, index, &val);
+            kernelSetElement(kel, i, j, val);
+            index++;
+        }
+    }
+
+    numaDestroy(&na);
+    return kel;
+}
+
+
+/*----------------------------------------------------------------------*
+ *                       Making a kernel from a Pix                     *
+ *----------------------------------------------------------------------*/
+/*!
+ * \brief   kernelCreateFromPix()
+ *
+ * \param[in]    pix
+ * \param[in]    cy, cx    origin of kernel
+ * \return  kernel, or NULL on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) The origin must be positive and within the dimensions of the pix.
+ * </pre>
+ */
+L_KERNEL *
+kernelCreateFromPix(PIX         *pix,
+                    l_int32      cy,
+                    l_int32      cx)
+{
+l_int32    i, j, w, h, d;
+l_uint32   val;
+L_KERNEL  *kel;
+
+    if (!pix)
+        return (L_KERNEL *)ERROR_PTR("pix not defined", __func__, NULL);
+    pixGetDimensions(pix, &w, &h, &d);
+    if (d != 8)
+        return (L_KERNEL *)ERROR_PTR("pix not 8 bpp", __func__, NULL);
+    if (cy < 0 || cx < 0 || cy >= h || cx >= w)
+        return (L_KERNEL *)ERROR_PTR("(cy, cx) invalid", __func__, NULL);
+
+    kel = kernelCreate(h, w);
+    kernelSetOrigin(kel, cy, cx);
+    for (i = 0; i < h; i++) {
+        for (j = 0; j < w; j++) {
+            pixGetPixel(pix, j, i, &val);
+            kernelSetElement(kel, i, j, (l_float32)val);
+        }
+    }
+
+    return kel;
+}
+
+
+/*----------------------------------------------------------------------*
+ *                     Display a kernel in a pix                        *
+ *----------------------------------------------------------------------*/
+/*!
+ * \brief   kernelDisplayInPix()
+ *
+ * \param[in]    kel       kernel
+ * \param[in]    size      of grid interiors; odd; either 1 or a minimum size
+ *                         of 17 is enforced
+ * \param[in]    gthick    grid thickness; either 0 or a minimum size of 2
+ *                         is enforced
+ * \return  pix   display of kernel, or NULL on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) This gives a visual representation of a kernel.
+ *      (2) There are two modes of display:
+ *          (a) Grid lines of minimum width 2, surrounding regions
+ *              representing kernel elements of minimum size 17,
+ *              with a "plus" mark at the kernel origin, or
+ *          (b) A pix without grid lines and using 1 pixel per kernel element.
+ *      (3) For both cases, the kernel absolute value is displayed,
+ *          normalized such that the maximum absolute value is 255.
+ *      (4) Large 2D separable kernels should be used for convolution
+ *          with two 1D kernels.  However, for the bilateral filter,
+ *          the computation time is independent of the size of the
+ *          2D content kernel.
+ * </pre>
+ */
+PIX *
+kernelDisplayInPix(L_KERNEL     *kel,
+                   l_int32       size,
+                   l_int32       gthick)
+{
+l_int32    i, j, w, h, sx, sy, cx, cy, width, x0, y0;
+l_int32    normval;
+l_float32  minval, maxval, max, val, norm;
+PIX       *pixd, *pixt0, *pixt1;
+
+    if (!kel)
+        return (PIX *)ERROR_PTR("kernel not defined", __func__, NULL);
+
+        /* Normalize the max value to be 255 for display */
+    kernelGetParameters(kel, &sy, &sx, &cy, &cx);
+    kernelGetMinMax(kel, &minval, &maxval);
+    max = L_MAX(maxval, -minval);
+    if (max == 0.0)
+        return (PIX *)ERROR_PTR("kernel elements all 0.0", __func__, NULL);
+    norm = 255.f / (l_float32)max;
+
+        /* Handle the 1 element/pixel case; typically with large kernels */
+    if (size == 1 && gthick == 0) {
+        pixd = pixCreate(sx, sy, 8);
+        for (i = 0; i < sy; i++) {
+            for (j = 0; j < sx; j++) {
+                kernelGetElement(kel, i, j, &val);
+                normval = (l_int32)(norm * L_ABS(val));
+                pixSetPixel(pixd, j, i, normval);
+            }
+        }
+        return pixd;
+    }
+
+        /* Enforce the constraints for the grid line version */
+    if (size < 17) {
+        L_WARNING("size < 17; setting to 17\n", __func__);
+        size = 17;
+    }
+    if (size % 2 == 0)
+        size++;
+    if (gthick < 2) {
+        L_WARNING("grid thickness < 2; setting to 2\n", __func__);
+        gthick = 2;
+    }
+
+    w = size * sx + gthick * (sx + 1);
+    h = size * sy + gthick * (sy + 1);
+    pixd = pixCreate(w, h, 8);
+
+        /* Generate grid lines */
+    for (i = 0; i <= sy; i++)
+        pixRenderLine(pixd, 0, gthick / 2 + i * (size + gthick),
+                      w - 1, gthick / 2 + i * (size + gthick),
+                      gthick, L_SET_PIXELS);
+    for (j = 0; j <= sx; j++)
+        pixRenderLine(pixd, gthick / 2 + j * (size + gthick), 0,
+                      gthick / 2 + j * (size + gthick), h - 1,
+                      gthick, L_SET_PIXELS);
+
+        /* Generate mask for each element */
+    pixt0 = pixCreate(size, size, 1);
+    pixSetAll(pixt0);
+
+        /* Generate crossed lines for origin pattern */
+    pixt1 = pixCreate(size, size, 1);
+    width = size / 8;
+    pixRenderLine(pixt1, size / 2, (l_int32)(0.12 * size),
+                           size / 2, (l_int32)(0.88 * size),
+                           width, L_SET_PIXELS);
+    pixRenderLine(pixt1, (l_int32)(0.15 * size), size / 2,
+                           (l_int32)(0.85 * size), size / 2,
+                           width, L_FLIP_PIXELS);
+    pixRasterop(pixt1, size / 2 - width, size / 2 - width,
+                2 * width, 2 * width, PIX_NOT(PIX_DST), NULL, 0, 0);
+
+        /* Paste the patterns in */
+    y0 = gthick;
+    for (i = 0; i < sy; i++) {
+        x0 = gthick;
+        for (j = 0; j < sx; j++) {
+            kernelGetElement(kel, i, j, &val);
+            normval = (l_int32)(norm * L_ABS(val));
+            pixSetMaskedGeneral(pixd, pixt0, normval, x0, y0);
+            if (i == cy && j == cx)
+                pixPaintThroughMask(pixd, pixt1, x0, y0, 255 - normval);
+            x0 += size + gthick;
+        }
+        y0 += size + gthick;
+    }
+
+    pixDestroy(&pixt0);
+    pixDestroy(&pixt1);
+    return pixd;
+}
+
+
+/*------------------------------------------------------------------------*
+ *                     Parse string to extract numbers                    *
+ *------------------------------------------------------------------------*/
+/*!
+ * \brief   parseStringForNumbers()
+ *
+ * \param[in]    str     string containing numbers; not changed
+ * \param[in]    seps    string of characters that can be used between ints
+ * \return  numa   of numbers found, or NULL on error
+ *
+ * <pre>
+ * Notes:
+ *     (1) The numbers can be ints or floats.
+ * </pre>
+ */
+NUMA *
+parseStringForNumbers(const char  *str,
+                      const char  *seps)
+{
+char      *newstr, *head;
+char      *tail = NULL;
+l_float32  val;
+NUMA      *na;
+
+    if (!str)
+        return (NUMA *)ERROR_PTR("str not defined", __func__, NULL);
+
+    newstr = stringNew(str);  /* to enforce const-ness of str */
+    na = numaCreate(0);
+    head = strtokSafe(newstr, seps, &tail);
+    val = atof(head);
+    numaAddNumber(na, val);
+    LEPT_FREE(head);
+    while ((head = strtokSafe(NULL, seps, &tail)) != NULL) {
+        val = atof(head);
+        numaAddNumber(na, val);
+        LEPT_FREE(head);
+    }
+
+    LEPT_FREE(newstr);
+    return na;
+}
+
+
+/*------------------------------------------------------------------------*
+ *                        Simple parametric kernels                       *
+ *------------------------------------------------------------------------*/
+/*!
+ * \brief   makeFlatKernel()
+ *
+ * \param[in]    height, width
+ * \param[in]    cy, cx          origin of kernel
+ * \return  kernel, or NULL on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) This is the same low-pass filtering kernel that is used
+ *          in the block convolution functions.
+ *      (2) The kernel origin (%cy, %cx) is typically placed as near
+ *          the center of the kernel as possible.  If height and
+ *          width are odd, then using %cy = height / 2 and
+ *          %cx = width / 2 places the origin at the exact center.
+ *      (3) This returns a normalized kernel.
+ * </pre>
+ */
+L_KERNEL *
+makeFlatKernel(l_int32  height,
+               l_int32  width,
+               l_int32  cy,
+               l_int32  cx)
+{
+l_int32    i, j;
+l_float32  normval;
+L_KERNEL  *kel;
+
+    if ((kel = kernelCreate(height, width)) == NULL)
+        return (L_KERNEL *)ERROR_PTR("kel not made", __func__, NULL);
+    kernelSetOrigin(kel, cy, cx);
+    normval = 1.0f / (l_float32)(height * width);
+    for (i = 0; i < height; i++) {
+        for (j = 0; j < width; j++) {
+            kernelSetElement(kel, i, j, normval);
+        }
+    }
+
+    return kel;
+}
+
+
+/*!
+ * \brief   makeGaussianKernel()
+ *
+ * \param[in]    halfh     sy = 2 * halfh + 1
+ * \param[in]    halfw     sx = 2 * halfw + 1
+ * \param[in]    stdev     standard deviation
+ * \param[in]    max       value at (cx,cy)
+ * \return  kernel, or NULL on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) The kernel size (sx, sy) = (2 * %halfw + 1, 2 * %halfh + 1)
+ *      (2) The kernel center (cx, cy) = (%halfw, %halfh).
+ *      (3) %halfw and %halfh are typically equal, and
+ *          are typically several times larger than the standard deviation.
+ *      (4) If pixConvolve() is invoked with normalization (the sum of
+ *          kernel elements = 1.0), use 1.0 for max (or any number that's
+ *          not too small or too large).
+ * </pre>
+ */
+L_KERNEL *
+makeGaussianKernel(l_int32    halfh,
+                   l_int32    halfw,
+                   l_float32  stdev,
+                   l_float32  max)
+{
+l_int32    sx, sy, i, j;
+l_float32  val;
+L_KERNEL  *kel;
+
+    sx = 2 * halfw + 1;
+    sy = 2 * halfh + 1;
+    if ((kel = kernelCreate(sy, sx)) == NULL)
+        return (L_KERNEL *)ERROR_PTR("kel not made", __func__, NULL);
+    kernelSetOrigin(kel, halfh, halfw);
+    for (i = 0; i < sy; i++) {
+        for (j = 0; j < sx; j++) {
+            val = expf(-(l_float32)((i - halfh) * (i - halfh) +
+                                    (j - halfw) * (j - halfw)) /
+                        (2. * stdev * stdev));
+            kernelSetElement(kel, i, j, max * val);
+        }
+    }
+
+    return kel;
+}
+
+
+/*!
+ * \brief   makeGaussianKernelSep()
+ *
+ * \param[in]    halfh     sy = 2 * halfh + 1
+ * \param[in]    halfw     sx = 2 * halfw + 1
+ * \param[in]    stdev     standard deviation
+ * \param[in]    max       value at (cx,cy)
+ * \param[out]   pkelx     x part of kernel
+ * \param[out]   pkely     y part of kernel
+ * \return  0 if OK, 1 on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) See makeGaussianKernel() for description of input parameters.
+ *      (2) These kernels are constructed so that the result of both
+ *          normalized and un-normalized convolution will be the same
+ *          as when convolving with pixConvolve() using the full kernel.
+ *      (3) The trick for the un-normalized convolution is to have the
+ *          product of the two kernel elements at (cx,cy) be equal to %max,
+ *          not max**2.  That's why %max for kely is 1.0.  If instead
+ *          we use sqrt(%max) for both, the results are slightly less
+ *          accurate, when compared to using the full kernel in
+ *          makeGaussianKernel().
+ * </pre>
+ */
+l_ok
+makeGaussianKernelSep(l_int32    halfh,
+                      l_int32    halfw,
+                      l_float32  stdev,
+                      l_float32  max,
+                      L_KERNEL **pkelx,
+                      L_KERNEL **pkely)
+{
+    if (!pkelx || !pkely)
+        return ERROR_INT("&kelx and &kely not defined", __func__, 1);
+
+    *pkelx = makeGaussianKernel(0, halfw, stdev, max);
+    *pkely = makeGaussianKernel(halfh, 0, stdev, 1.0);
+    return 0;
+}
+
+
+/*!
+ * \brief   makeDoGKernel()
+ *
+ * \param[in]    halfh     sy = 2 * halfh + 1
+ * \param[in]    halfw     sx = 2 * halfw + 1
+ * \param[in]    stdev     standard deviation of narrower gaussian
+ * \param[in]    ratio     of stdev for wide filter to stdev for narrow one
+ * \return  kernel, or NULL on error
+ *
+ * <pre>
+ * Notes:
+ *      (1) The DoG (difference of gaussians) is a wavelet mother
+ *          function with null total sum.  By subtracting two blurred
+ *          versions of the image, it acts as a bandpass filter for
+ *          frequencies passed by the narrow gaussian but stopped
+ *          by the wide one.See:
+ *               http://en.wikipedia.org/wiki/Difference_of_Gaussians
+ *      (2) The kernel size (sx, sy) = (2 * halfw + 1, 2 * halfh + 1).
+ *      (3) The kernel center (cx, cy) = (halfw, halfh).
+ *      (4) %halfw and %halfh are typically equal, and are typically
+ *          several times larger than the standard deviation.
+ *      (5) %ratio is the ratio of standard deviations of the wide
+ *          to narrow gaussian.  It must be >= 1.0; 1.0 is a no-op.
+ *      (6) Because the kernel is a null sum, it must be invoked without
+ *          normalization in pixConvolve().
+ * </pre>
+ */
+L_KERNEL *
+makeDoGKernel(l_int32    halfh,
+              l_int32    halfw,
+              l_float32  stdev,
+              l_float32  ratio)
+{
+l_int32    sx, sy, i, j;
+l_float32  pi, squaredist, highnorm, lownorm, val;
+L_KERNEL  *kel;
+
+    sx = 2 * halfw + 1;
+    sy = 2 * halfh + 1;
+    if ((kel = kernelCreate(sy, sx)) == NULL)
+        return (L_KERNEL *)ERROR_PTR("kel not made", __func__, NULL);
+    kernelSetOrigin(kel, halfh, halfw);
+
+    pi = 3.1415926535f;
+    for (i = 0; i < sy; i++) {
+        for (j = 0; j < sx; j++) {
+            squaredist = (l_float32)((i - halfh) * (i - halfh) +
+                                     (j - halfw) * (j - halfw));
+            highnorm = 1.f / (2 * stdev * stdev);
+            lownorm = highnorm / (ratio * ratio);
+            val = (highnorm / pi) * expf(-(highnorm * squaredist))
+                  - (lownorm / pi) * expf(-(lownorm * squaredist));
+            kernelSetElement(kel, i, j, val);
+        }
+    }
+
+    return kel;
+}