kernel.c

Go to the documentation of this file.

/*====================================================================*
 -  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.
 *====================================================================*/


#include <string.h>
#include <math.h>
#include "allheaders.h"


/*------------------------------------------------------------------------*
 *                           Create / Destroy                             *
 *------------------------------------------------------------------------*/
L_KERNEL *
kernelCreate(l_int32  height,
             l_int32  width)
{
l_uint64   size64;
L_KERNEL  *kel;

    PROCNAME("kernelCreate");

    if (width <= 0)
        return (L_KERNEL *)ERROR_PTR("width must be > 0", procName, NULL);
    if (height <= 0)
        return (L_KERNEL *)ERROR_PTR("height must be > 0", procName, 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", procName, width, height);
        return (L_KERNEL *)ERROR_PTR("size >= 2^29", procName, 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", procName, NULL);
    }
    return kel;
}


void
kernelDestroy(L_KERNEL  **pkel)
{
l_int32    i;
L_KERNEL  *kel;

    PROCNAME("kernelDestroy");

    if (pkel == NULL)  {
        L_WARNING("ptr address is NULL!\n", procName);
        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;
    return;
}


L_KERNEL *
kernelCopy(L_KERNEL  *kels)
{
l_int32    i, j, sx, sy, cx, cy;
L_KERNEL  *keld;

    PROCNAME("kernelCopy");

    if (!kels)
        return (L_KERNEL *)ERROR_PTR("kels not defined", procName, NULL);

    kernelGetParameters(kels, &sy, &sx, &cy, &cx);
    if ((keld = kernelCreate(sy, sx)) == NULL)
        return (L_KERNEL *)ERROR_PTR("keld not made", procName, 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                              *
 *----------------------------------------------------------------------*/
l_ok
kernelGetElement(L_KERNEL   *kel,
                 l_int32     row,
                 l_int32     col,
                 l_float32  *pval)
{
    PROCNAME("kernelGetElement");

    if (!pval)
        return ERROR_INT("&val not defined", procName, 1);
    *pval = 0;
    if (!kel)
        return ERROR_INT("kernel not defined", procName, 1);
    if (row < 0 || row >= kel->sy)
        return ERROR_INT("kernel row out of bounds", procName, 1);
    if (col < 0 || col >= kel->sx)
        return ERROR_INT("kernel col out of bounds", procName, 1);

    *pval = kel->data[row][col];
    return 0;
}


l_ok
kernelSetElement(L_KERNEL  *kel,
                 l_int32    row,
                 l_int32    col,
                 l_float32  val)
{
    PROCNAME("kernelSetElement");

    if (!kel)
        return ERROR_INT("kel not defined", procName, 1);
    if (row < 0 || row >= kel->sy)
        return ERROR_INT("kernel row out of bounds", procName, 1);
    if (col < 0 || col >= kel->sx)
        return ERROR_INT("kernel col out of bounds", procName, 1);

    kel->data[row][col] = val;
    return 0;
}


l_ok
kernelGetParameters(L_KERNEL  *kel,
                    l_int32   *psy,
                    l_int32   *psx,
                    l_int32   *pcy,
                    l_int32   *pcx)
{
    PROCNAME("kernelGetParameters");

    if (psy) *psy = 0;
    if (psx) *psx = 0;
    if (pcy) *pcy = 0;
    if (pcx) *pcx = 0;
    if (!kel)
        return ERROR_INT("kernel not defined", procName, 1);
    if (psy) *psy = kel->sy;
    if (psx) *psx = kel->sx;
    if (pcy) *pcy = kel->cy;
    if (pcx) *pcx = kel->cx;
    return 0;
}


l_ok
kernelSetOrigin(L_KERNEL  *kel,
                l_int32    cy,
                l_int32    cx)
{
    PROCNAME("kernelSetOrigin");

    if (!kel)
        return ERROR_INT("kel not defined", procName, 1);
    kel->cy = cy;
    kel->cx = cx;
    return 0;
}


l_ok
kernelGetSum(L_KERNEL   *kel,
             l_float32  *psum)
{
l_int32    sx, sy, i, j;

    PROCNAME("kernelGetSum");

    if (!psum)
        return ERROR_INT("&sum not defined", procName, 1);
    *psum = 0.0;
    if (!kel)
        return ERROR_INT("kernel not defined", procName, 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;
}


l_ok
kernelGetMinMax(L_KERNEL   *kel,
                l_float32  *pmin,
                l_float32  *pmax)
{
l_int32    sx, sy, i, j;
l_float32  val, minval, maxval;

    PROCNAME("kernelGetMinmax");

    if (!pmin && !pmax)
        return ERROR_INT("neither &min nor &max defined", procName, 1);
    if (pmin) *pmin = 0.0;
    if (pmax) *pmax = 0.0;
    if (!kel)
        return ERROR_INT("kernel not defined", procName, 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                            *
 *----------------------------------------------------------------------*/
L_KERNEL *
kernelNormalize(L_KERNEL  *kels,
                l_float32  normsum)
{
l_int32    i, j, sx, sy, cx, cy;
l_float32  sum, factor;
L_KERNEL  *keld;

    PROCNAME("kernelNormalize");

    if (!kels)
        return (L_KERNEL *)ERROR_PTR("kels not defined", procName, NULL);

    kernelGetSum(kels, &sum);
    if (L_ABS(sum) < 0.00001) {
        L_WARNING("null sum; not normalizing; returning a copy\n", procName);
        return kernelCopy(kels);
    }

    kernelGetParameters(kels, &sy, &sx, &cy, &cx);
    if ((keld = kernelCreate(sy, sx)) == NULL)
        return (L_KERNEL *)ERROR_PTR("keld not made", procName, 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;
}


L_KERNEL *
kernelInvert(L_KERNEL  *kels)
{
l_int32    i, j, sx, sy, cx, cy;
L_KERNEL  *keld;

    PROCNAME("kernelInvert");

    if (!kels)
        return (L_KERNEL *)ERROR_PTR("kels not defined", procName, NULL);

    kernelGetParameters(kels, &sy, &sx, &cy, &cx);
    if ((keld = kernelCreate(sy, sx)) == NULL)
        return (L_KERNEL *)ERROR_PTR("keld not made", procName, 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                           *
 *----------------------------------------------------------------------*/
l_float32 **
create2dFloatArray(l_int32  sy,
                   l_int32  sx)
{
l_int32      i;
l_float32  **array;

    PROCNAME("create2dFloatArray");

    if ((array = (l_float32 **)LEPT_CALLOC(sy, sizeof(l_float32 *))) == NULL)
        return (l_float32 **)ERROR_PTR("ptr array not made", procName, NULL);

    for (i = 0; i < sy; i++)
        array[i] = (l_float32 *)LEPT_CALLOC(sx, sizeof(l_float32));
    return array;
}


/*----------------------------------------------------------------------*
 *                            Kernel serialized I/O                     *
 *----------------------------------------------------------------------*/
L_KERNEL *
kernelRead(const char  *fname)
{
FILE      *fp;
L_KERNEL  *kel;

    PROCNAME("kernelRead");

    if (!fname)
        return (L_KERNEL *)ERROR_PTR("fname not defined", procName, NULL);

    if ((fp = fopenReadStream(fname)) == NULL)
        return (L_KERNEL *)ERROR_PTR("stream not opened", procName, NULL);
    if ((kel = kernelReadStream(fp)) == NULL) {
        fclose(fp);
        return (L_KERNEL *)ERROR_PTR("kel not returned", procName, NULL);
    }
    fclose(fp);

    return kel;
}


L_KERNEL *
kernelReadStream(FILE  *fp)
{
l_int32    sy, sx, cy, cx, i, j, ret, version, ignore;
L_KERNEL  *kel;

    PROCNAME("kernelReadStream");

    if (!fp)
        return (L_KERNEL *)ERROR_PTR("stream not defined", procName, NULL);

    ret = fscanf(fp, "  Kernel Version %d\n", &version);
    if (ret != 1)
        return (L_KERNEL *)ERROR_PTR("not a kernel file", procName, NULL);
    if (version != KERNEL_VERSION_NUMBER)
        return (L_KERNEL *)ERROR_PTR("invalid kernel version", procName, 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", procName, NULL);

    if ((kel = kernelCreate(sy, sx)) == NULL)
        return (L_KERNEL *)ERROR_PTR("kel not made", procName, 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;
}


l_ok
kernelWrite(const char  *fname,
            L_KERNEL    *kel)
{
FILE  *fp;

    PROCNAME("kernelWrite");

    if (!fname)
        return ERROR_INT("fname not defined", procName, 1);
    if (!kel)
        return ERROR_INT("kel not defined", procName, 1);

    if ((fp = fopenWriteStream(fname, "wb")) == NULL)
        return ERROR_INT("stream not opened", procName, 1);
    kernelWriteStream(fp, kel);
    fclose(fp);

    return 0;
}


l_ok
kernelWriteStream(FILE      *fp,
                  L_KERNEL  *kel)
{
l_int32  sx, sy, cx, cy, i, j;

    PROCNAME("kernelWriteStream");

    if (!fp)
        return ERROR_INT("stream not defined", procName, 1);
    if (!kel)
        return ERROR_INT("kel not defined", procName, 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               *
 *----------------------------------------------------------------------*/
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;

    PROCNAME("kernelCreateFromString");

    if (h < 1)
        return (L_KERNEL *)ERROR_PTR("height must be > 0", procName, NULL);
    if (w < 1)
        return (L_KERNEL *)ERROR_PTR("width must be > 0", procName, NULL);
    if (cy < 0 || cy >= h)
        return (L_KERNEL *)ERROR_PTR("cy invalid", procName, NULL);
    if (cx < 0 || cx >= w)
        return (L_KERNEL *)ERROR_PTR("cx invalid", procName, 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);
        fprintf(stderr, "w = %d, h = %d, num ints = %d\n", w, h, n);
        return (L_KERNEL *)ERROR_PTR("invalid integer data", procName, 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             *
 *----------------------------------------------------------------------*/
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;

    PROCNAME("kernelCreateFromFile");

    if (!filename)
        return (L_KERNEL *)ERROR_PTR("filename not defined", procName, NULL);

    if ((filestr = (char *)l_binaryRead(filename, &size)) == NULL)
        return (L_KERNEL *)ERROR_PTR("file not found", procName, NULL);
    if (size == 0) {
        LEPT_FREE(filestr);
        return (L_KERNEL *)ERROR_PTR("file is empty", procName, 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", procName, 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", procName, 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);
        fprintf(stderr, "w = %d, h = %d, num ints = %d\n", w, h, n);
        return (L_KERNEL *)ERROR_PTR("invalid integer data", procName, 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                     *
 *----------------------------------------------------------------------*/
L_KERNEL *
kernelCreateFromPix(PIX         *pix,
                    l_int32      cy,
                    l_int32      cx)
{
l_int32    i, j, w, h, d;
l_uint32   val;
L_KERNEL  *kel;

    PROCNAME("kernelCreateFromPix");

    if (!pix)
        return (L_KERNEL *)ERROR_PTR("pix not defined", procName, NULL);
    pixGetDimensions(pix, &w, &h, &d);
    if (d != 8)
        return (L_KERNEL *)ERROR_PTR("pix not 8 bpp", procName, NULL);
    if (cy < 0 || cx < 0 || cy >= h || cx >= w)
        return (L_KERNEL *)ERROR_PTR("(cy, cx) invalid", procName, 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                        *
 *----------------------------------------------------------------------*/
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;

    PROCNAME("kernelDisplayInPix");

    if (!kel)
        return (PIX *)ERROR_PTR("kernel not defined", procName, 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", procName, NULL);
    norm = 255. / (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", procName);
        size = 17;
    }
    if (size % 2 == 0)
        size++;
    if (gthick < 2) {
        L_WARNING("grid thickness < 2; setting to 2\n", procName);
        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                    *
 *------------------------------------------------------------------------*/
NUMA *
parseStringForNumbers(const char  *str,
                      const char  *seps)
{
char      *newstr, *head;
char      *tail = NULL;
l_float32  val;
NUMA      *na;

    PROCNAME("parseStringForNumbers");

    if (!str)
        return (NUMA *)ERROR_PTR("str not defined", procName, 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                       *
 *------------------------------------------------------------------------*/
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;

    PROCNAME("makeFlatKernel");

    if ((kel = kernelCreate(height, width)) == NULL)
        return (L_KERNEL *)ERROR_PTR("kel not made", procName, NULL);
    kernelSetOrigin(kel, cy, cx);
    normval = 1.0 / (l_float32)(height * width);
    for (i = 0; i < height; i++) {
        for (j = 0; j < width; j++) {
            kernelSetElement(kel, i, j, normval);
        }
    }

    return kel;
}


L_KERNEL *
makeGaussianKernel(l_int32    halfheight,
                   l_int32    halfwidth,
                   l_float32  stdev,
                   l_float32  max)
{
l_int32    sx, sy, i, j;
l_float32  val;
L_KERNEL  *kel;

    PROCNAME("makeGaussianKernel");

    sx = 2 * halfwidth + 1;
    sy = 2 * halfheight + 1;
    if ((kel = kernelCreate(sy, sx)) == NULL)
        return (L_KERNEL *)ERROR_PTR("kel not made", procName, NULL);
    kernelSetOrigin(kel, halfheight, halfwidth);
    for (i = 0; i < sy; i++) {
        for (j = 0; j < sx; j++) {
            val = expf(-(l_float32)((i - halfheight) * (i - halfheight) +
                                    (j - halfwidth) * (j - halfwidth)) /
                        (2. * stdev * stdev));
            kernelSetElement(kel, i, j, max * val);
        }
    }

    return kel;
}


l_ok
makeGaussianKernelSep(l_int32    halfheight,
                      l_int32    halfwidth,
                      l_float32  stdev,
                      l_float32  max,
                      L_KERNEL **pkelx,
                      L_KERNEL **pkely)
{
    PROCNAME("makeGaussianKernelSep");

    if (!pkelx || !pkely)
        return ERROR_INT("&kelx and &kely not defined", procName, 1);

    *pkelx = makeGaussianKernel(0, halfwidth, stdev, max);
    *pkely = makeGaussianKernel(halfheight, 0, stdev, 1.0);
    return 0;
}


L_KERNEL *
makeDoGKernel(l_int32    halfheight,
              l_int32    halfwidth,
              l_float32  stdev,
              l_float32  ratio)
{
l_int32    sx, sy, i, j;
l_float32  pi, squaredist, highnorm, lownorm, val;
L_KERNEL  *kel;

    PROCNAME("makeDoGKernel");

    sx = 2 * halfwidth + 1;
    sy = 2 * halfheight + 1;
    if ((kel = kernelCreate(sy, sx)) == NULL)
        return (L_KERNEL *)ERROR_PTR("kel not made", procName, NULL);
    kernelSetOrigin(kel, halfheight, halfwidth);

    pi = 3.1415926535;
    for (i = 0; i < sy; i++) {
        for (j = 0; j < sx; j++) {
            squaredist = (l_float32)((i - halfheight) * (i - halfheight) +
                                     (j - halfwidth) * (j - halfwidth));
            highnorm = 1. / (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;
}