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

 
#ifdef HAVE_CONFIG_H
#include <config_auto.h>
#endif  /* HAVE_CONFIG_H */
 
#include <math.h>
#include "allheaders.h"
 
    /* bin sizes in numaMakeHistogram() */
static const l_int32 BinSizeArray[] = {2, 5, 10, 20, 50, 100, 200, 500, 1000,\
                      2000, 5000, 10000, 20000, 50000, 100000, 200000,\
                      500000, 1000000, 2000000, 5000000, 10000000,\
                      200000000, 50000000, 100000000};
static const l_int32 NBinSizes = 24;
 
 
#ifndef  NO_CONSOLE_IO
#define  DEBUG_HISTO        0
#define  DEBUG_CROSSINGS    0
#define  DEBUG_FREQUENCY    0
#endif  /* ~NO_CONSOLE_IO */
 
/*----------------------------------------------------------------------*
 *                     Morphological operations                         *
 *----------------------------------------------------------------------*/
NUMA *
numaErode(NUMA    *nas,
          l_int32  size)
{
l_int32     i, j, n, hsize, len;
l_float32   minval;
l_float32  *fa, *fas, *fad;
NUMA       *nad;
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
    if (size <= 0)
        return (NUMA *)ERROR_PTR("size must be > 0", __func__, NULL);
    if ((size & 1) == 0 ) {
        L_WARNING("sel size must be odd; increasing by 1\n", __func__);
        size++;
    }
 
    if (size == 1)
        return numaCopy(nas);
 
        /* Make a source fa (fas) that has an added (size / 2) boundary
         * on left and right, contains a copy of nas in the interior region
         * (between 'size' and 'size + n', and has large values
         * inserted in the boundary (because it is an erosion). */
    n = numaGetCount(nas);
    hsize = size / 2;
    len = n + 2 * hsize;
    if ((fas = (l_float32 *)LEPT_CALLOC(len, sizeof(l_float32))) == NULL)
        return (NUMA *)ERROR_PTR("fas not made", __func__, NULL);
    for (i = 0; i < hsize; i++)
         fas[i] = 1.0e37;
    for (i = hsize + n; i < len; i++)
         fas[i] = 1.0e37;
    fa = numaGetFArray(nas, L_NOCOPY);
    for (i = 0; i < n; i++)
         fas[hsize + i] = fa[i];
 
    nad = numaMakeConstant(0, n);
    numaCopyParameters(nad, nas);
    fad = numaGetFArray(nad, L_NOCOPY);
    for (i = 0; i < n; i++) {
        minval = 1.0e37;  /* start big */
        for (j = 0; j < size; j++)
            minval = L_MIN(minval, fas[i + j]);
        fad[i] = minval;
    }
 
    LEPT_FREE(fas);
    return nad;
}
 

NUMA *
numaDilate(NUMA    *nas,
           l_int32  size)
{
l_int32     i, j, n, hsize, len;
l_float32   maxval;
l_float32  *fa, *fas, *fad;
NUMA       *nad;
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
    if (size <= 0)
        return (NUMA *)ERROR_PTR("size must be > 0", __func__, NULL);
    if ((size & 1) == 0 ) {
        L_WARNING("sel size must be odd; increasing by 1\n", __func__);
        size++;
    }
 
    if (size == 1)
        return numaCopy(nas);
 
        /* Make a source fa (fas) that has an added (size / 2) boundary
         * on left and right, contains a copy of nas in the interior region
         * (between 'size' and 'size + n', and has small values
         * inserted in the boundary (because it is a dilation). */
    n = numaGetCount(nas);
    hsize = size / 2;
    len = n + 2 * hsize;
    if ((fas = (l_float32 *)LEPT_CALLOC(len, sizeof(l_float32))) == NULL)
        return (NUMA *)ERROR_PTR("fas not made", __func__, NULL);
    for (i = 0; i < hsize; i++)
         fas[i] = -1.0e37;
    for (i = hsize + n; i < len; i++)
         fas[i] = -1.0e37;
    fa = numaGetFArray(nas, L_NOCOPY);
    for (i = 0; i < n; i++)
         fas[hsize + i] = fa[i];
 
    nad = numaMakeConstant(0, n);
    numaCopyParameters(nad, nas);
    fad = numaGetFArray(nad, L_NOCOPY);
    for (i = 0; i < n; i++) {
        maxval = -1.0e37;  /* start small */
        for (j = 0; j < size; j++)
            maxval = L_MAX(maxval, fas[i + j]);
        fad[i] = maxval;
    }
 
    LEPT_FREE(fas);
    return nad;
}
 

NUMA *
numaOpen(NUMA    *nas,
         l_int32  size)
{
NUMA  *nat, *nad;
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
    if (size <= 0)
        return (NUMA *)ERROR_PTR("size must be > 0", __func__, NULL);
    if ((size & 1) == 0 ) {
        L_WARNING("sel size must be odd; increasing by 1\n", __func__);
        size++;
    }
 
    if (size == 1)
        return numaCopy(nas);
 
    nat = numaErode(nas, size);
    nad = numaDilate(nat, size);
    numaDestroy(&nat);
    return nad;
}
 

NUMA *
numaClose(NUMA    *nas,
          l_int32  size)
{
NUMA  *nab, *nat1, *nat2, *nad;
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
    if (size <= 0)
        return (NUMA *)ERROR_PTR("size must be > 0", __func__, NULL);
    if ((size & 1) == 0 ) {
        L_WARNING("sel size must be odd; increasing by 1\n", __func__);
        size++;
    }
 
    if (size == 1)
        return numaCopy(nas);
 
    nab = numaAddBorder(nas, size, size, 0);  /* to preserve extensivity */
    nat1 = numaDilate(nab, size);
    nat2 = numaErode(nat1, size);
    nad = numaRemoveBorder(nat2, size, size);
    numaDestroy(&nab);
    numaDestroy(&nat1);
    numaDestroy(&nat2);
    return nad;
}
 
 
/*----------------------------------------------------------------------*
 *                            Other transforms                          *
 *----------------------------------------------------------------------*/
NUMA *
numaTransform(NUMA      *nas,
              l_float32  shift,
              l_float32  scale)
{
l_int32    i, n;
l_float32  val;
NUMA      *nad;
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
    n = numaGetCount(nas);
    if ((nad = numaCreate(n)) == NULL)
        return (NUMA *)ERROR_PTR("nad not made", __func__, NULL);
    numaCopyParameters(nad, nas);
    for (i = 0; i < n; i++) {
        numaGetFValue(nas, i, &val);
        val = scale * (val + shift);
        numaAddNumber(nad, val);
    }
    return nad;
}
 

l_ok
numaSimpleStats(NUMA       *na,
                l_int32     first,
                l_int32     last,
                l_float32  *pmean,
                l_float32  *pvar,
                l_float32  *prvar)
{
l_int32    i, n, ni;
l_float32  sum, sumsq, val, mean, var;
 
    if (pmean) *pmean = 0.0;
    if (pvar) *pvar = 0.0;
    if (prvar) *prvar = 0.0;
    if (!pmean && !pvar && !prvar)
        return ERROR_INT("nothing requested", __func__, 1);
    if (!na)
        return ERROR_INT("na not defined", __func__, 1);
    if ((n = numaGetCount(na)) == 0)
        return ERROR_INT("na is empty", __func__, 1);
    first = L_MAX(0, first);
    if (last < 0) last = n - 1;
    if (first >= n)
        return ERROR_INT("invalid first", __func__, 1);
    if (last >= n) {
        L_WARNING("last = %d is beyond max index = %d; adjusting\n",
                  __func__, last, n - 1);
        last = n - 1;
    }
    if (first > last)
        return ERROR_INT("first > last\n", __func__, 1);
    ni = last - first + 1;
    sum = sumsq = 0.0;
    for (i = first; i <= last; i++) {
        numaGetFValue(na, i, &val);
        sum += val;
        sumsq += val * val;
    }
 
    mean = sum / ni;
    if (pmean)
        *pmean = mean;
    if (pvar || prvar) {
        var = sumsq / ni - mean * mean;
        if (pvar) *pvar = var;
        if (prvar) *prvar = sqrtf(var);
    }
 
    return 0;
}
 

l_ok
numaWindowedStats(NUMA    *nas,
                  l_int32  wc,
                  NUMA   **pnam,
                  NUMA   **pnams,
                  NUMA   **pnav,
                  NUMA   **pnarv)
{
NUMA  *nam, *nams;
 
    if (!nas)
        return ERROR_INT("nas not defined", __func__, 1);
    if (2 * wc + 1 > numaGetCount(nas))
        L_WARNING("filter wider than input array!\n", __func__);
 
    if (!pnav && !pnarv) {
        if (pnam) *pnam = numaWindowedMean(nas, wc);
        if (pnams) *pnams = numaWindowedMeanSquare(nas, wc);
        return 0;
    }
 
    nam = numaWindowedMean(nas, wc);
    nams = numaWindowedMeanSquare(nas, wc);
    numaWindowedVariance(nam, nams, pnav, pnarv);
    if (pnam)
        *pnam = nam;
    else
        numaDestroy(&nam);
    if (pnams)
        *pnams = nams;
    else
        numaDestroy(&nams);
    return 0;
}
 

NUMA *
numaWindowedMean(NUMA    *nas,
                 l_int32  wc)
{
l_int32     i, n, n1, width;
l_float32   sum, norm;
l_float32  *fa1, *fad, *suma;
NUMA       *na1, *nad;
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
    n = numaGetCount(nas);
    width = 2 * wc + 1;  /* filter width */
    if (width > n)
        L_WARNING("filter wider than input array!\n", __func__);
 
    na1 = numaAddSpecifiedBorder(nas, wc, wc, L_MIRRORED_BORDER);
    n1 = n + 2 * wc;
    fa1 = numaGetFArray(na1, L_NOCOPY);
    nad = numaMakeConstant(0, n);
    fad = numaGetFArray(nad, L_NOCOPY);
 
        /* Make sum array; note the indexing */
    if ((suma = (l_float32 *)LEPT_CALLOC(n1 + 1, sizeof(l_float32))) == NULL) {
        numaDestroy(&na1);
        numaDestroy(&nad);
        return (NUMA *)ERROR_PTR("suma not made", __func__, NULL);
    }
    sum = 0.0;
    suma[0] = 0.0;
    for (i = 0; i < n1; i++) {
        sum += fa1[i];
        suma[i + 1] = sum;
    }
 
    norm = 1. / (2 * wc + 1);
    for (i = 0; i < n; i++)
        fad[i] = norm * (suma[width + i] - suma[i]);
 
    LEPT_FREE(suma);
    numaDestroy(&na1);
    return nad;
}
 

NUMA *
numaWindowedMeanSquare(NUMA    *nas,
                       l_int32  wc)
{
l_int32     i, n, n1, width;
l_float32   sum, norm;
l_float32  *fa1, *fad, *suma;
NUMA       *na1, *nad;
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
    n = numaGetCount(nas);
    width = 2 * wc + 1;  /* filter width */
    if (width > n)
        L_WARNING("filter wider than input array!\n", __func__);
 
    na1 = numaAddSpecifiedBorder(nas, wc, wc, L_MIRRORED_BORDER);
    n1 = n + 2 * wc;
    fa1 = numaGetFArray(na1, L_NOCOPY);
    nad = numaMakeConstant(0, n);
    fad = numaGetFArray(nad, L_NOCOPY);
 
        /* Make sum array; note the indexing */
    if ((suma = (l_float32 *)LEPT_CALLOC(n1 + 1, sizeof(l_float32))) == NULL) {
        numaDestroy(&na1);
        numaDestroy(&nad);
        return (NUMA *)ERROR_PTR("suma not made", __func__, NULL);
    }
    sum = 0.0;
    suma[0] = 0.0;
    for (i = 0; i < n1; i++) {
        sum += fa1[i] * fa1[i];
        suma[i + 1] = sum;
    }
 
    norm = 1. / (2 * wc + 1);
    for (i = 0; i < n; i++)
        fad[i] = norm * (suma[width + i] - suma[i]);
 
    LEPT_FREE(suma);
    numaDestroy(&na1);
    return nad;
}
 

l_ok
numaWindowedVariance(NUMA   *nam,
                     NUMA   *nams,
                     NUMA  **pnav,
                     NUMA  **pnarv)
{
l_int32     i, nm, nms;
l_float32   var;
l_float32  *fam, *fams, *fav, *farv;
NUMA       *nav, *narv;  /* variance and square root of variance */
 
    if (pnav) *pnav = NULL;
    if (pnarv) *pnarv = NULL;
    if (!pnav && !pnarv)
        return ERROR_INT("neither &nav nor &narv are defined", __func__, 1);
    if (!nam)
        return ERROR_INT("nam not defined", __func__, 1);
    if (!nams)
        return ERROR_INT("nams not defined", __func__, 1);
    nm = numaGetCount(nam);
    nms = numaGetCount(nams);
    if (nm != nms)
        return ERROR_INT("sizes of nam and nams differ", __func__, 1);
 
    if (pnav) {
        nav = numaMakeConstant(0, nm);
        *pnav = nav;
        fav = numaGetFArray(nav, L_NOCOPY);
    }
    if (pnarv) {
        narv = numaMakeConstant(0, nm);
        *pnarv = narv;
        farv = numaGetFArray(narv, L_NOCOPY);
    }
    fam = numaGetFArray(nam, L_NOCOPY);
    fams = numaGetFArray(nams, L_NOCOPY);
 
    for (i = 0; i < nm; i++) {
        var = fams[i] - fam[i] * fam[i];
        if (pnav)
            fav[i] = var;
        if (pnarv)
            farv[i] = sqrtf(var);
    }
 
    return 0;
}
 

NUMA *
numaWindowedMedian(NUMA    *nas,
                   l_int32  halfwin)
{
l_int32    i, n;
l_float32  medval;
NUMA      *na1, *na2, *nad;
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
    if ((n = numaGetCount(nas)) < 3)
        return numaCopy(nas);
    if (halfwin <= 0) {
        L_ERROR("filter too small; returning a copy\n", __func__);
        return numaCopy(nas);
    }
 
    if (halfwin > (n - 1) / 2) {
        halfwin = (n - 1) / 2;
        L_INFO("reducing filter to halfwin = %d\n", __func__, halfwin);
    }
 
        /* Add a border to both ends */
    na1 = numaAddSpecifiedBorder(nas, halfwin, halfwin, L_MIRRORED_BORDER);
 
        /* Get the median value at the center of each window, corresponding
         * to locations in the input nas. */
    nad = numaCreate(n);
    for (i = 0; i < n; i++) {
        na2 = numaClipToInterval(na1, i, i + 2 * halfwin);
        numaGetMedian(na2, &medval);
        numaAddNumber(nad, medval);
        numaDestroy(&na2);
    }
 
    numaDestroy(&na1);
    return nad;
}
 

NUMA *
numaConvertToInt(NUMA  *nas)
{
l_int32  i, n, ival;
NUMA    *nad;
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
 
    n = numaGetCount(nas);
    if ((nad = numaCreate(n)) == NULL)
        return (NUMA *)ERROR_PTR("nad not made", __func__, NULL);
    numaCopyParameters(nad, nas);
    for (i = 0; i < n; i++) {
        numaGetIValue(nas, i, &ival);
        numaAddNumber(nad, ival);
    }
    return nad;
}
 
 
/*----------------------------------------------------------------------*
 *                 Histogram generation and statistics                  *
 *----------------------------------------------------------------------*/
NUMA *
numaMakeHistogram(NUMA     *na,
                  l_int32   maxbins,
                  l_int32  *pbinsize,
                  l_int32  *pbinstart)
{
l_int32    i, n, ival, hval;
l_int32    iminval, imaxval, range, binsize, nbins, ibin;
l_float32  val, ratio;
NUMA      *nai, *nahist;
 
    if (pbinsize) *pbinsize = 0;
    if (pbinstart) *pbinstart = 0;
    if (!na)
        return (NUMA *)ERROR_PTR("na not defined", __func__, NULL);
    if (maxbins < 1)
        return (NUMA *)ERROR_PTR("maxbins < 1", __func__, NULL);
 
        /* Determine input range */
    numaGetMin(na, &val, NULL);
    iminval = (l_int32)(val + 0.5);
    numaGetMax(na, &val, NULL);
    imaxval = (l_int32)(val + 0.5);
    if (pbinstart == NULL) {  /* clip negative vals; start from 0 */
        iminval = 0;
        if (imaxval < 0)
            return (NUMA *)ERROR_PTR("all values < 0", __func__, NULL);
    }
 
        /* Determine binsize */
    range = imaxval - iminval + 1;
    if (range > maxbins - 1) {
        ratio = (l_float64)range / (l_float64)maxbins;
        binsize = 0;
        for (i = 0; i < NBinSizes; i++) {
            if (ratio < BinSizeArray[i]) {
                binsize = BinSizeArray[i];
                break;
            }
        }
        if (binsize == 0)
            return (NUMA *)ERROR_PTR("numbers too large", __func__, NULL);
    } else {
        binsize = 1;
    }
    if (pbinsize) *pbinsize = binsize;
    nbins = 1 + range / binsize;  /* +1 seems to be sufficient */
 
        /* Redetermine iminval */
    if (pbinstart && binsize > 1) {
        if (iminval >= 0)
            iminval = binsize * (iminval / binsize);
        else
            iminval = binsize * ((iminval - binsize + 1) / binsize);
    }
    if (pbinstart) *pbinstart = iminval;
 
#if  DEBUG_HISTO
    lept_stderr(" imaxval = %d, range = %d, nbins = %d\n",
                imaxval, range, nbins);
#endif  /* DEBUG_HISTO */
 
        /* Use integerized data for input */
    if ((nai = numaConvertToInt(na)) == NULL)
        return (NUMA *)ERROR_PTR("nai not made", __func__, NULL);
    n = numaGetCount(nai);
 
        /* Make histogram, converting value in input array
         * into a bin number for this histogram array. */
    if ((nahist = numaCreate(nbins)) == NULL) {
        numaDestroy(&nai);
        return (NUMA *)ERROR_PTR("nahist not made", __func__, NULL);
    }
    numaSetCount(nahist, nbins);
    numaSetParameters(nahist, iminval, binsize);
    for (i = 0; i < n; i++) {
        numaGetIValue(nai, i, &ival);
        ibin = (ival - iminval) / binsize;
        if (ibin >= 0 && ibin < nbins) {
            numaGetIValue(nahist, ibin, &hval);
            numaSetValue(nahist, ibin, hval + 1.0);
        }
    }
 
    numaDestroy(&nai);
    return nahist;
}
 

NUMA *
numaMakeHistogramAuto(NUMA    *na,
                      l_int32  maxbins)
{
l_int32    i, n, imin, imax, irange, ibin, ival, allints;
l_float32  minval, maxval, range, binsize, fval;
NUMA      *nah;
 
    if (!na)
        return (NUMA *)ERROR_PTR("na not defined", __func__, NULL);
    maxbins = L_MAX(1, maxbins);
 
        /* Determine input range */
    numaGetMin(na, &minval, NULL);
    numaGetMax(na, &maxval, NULL);
 
        /* Determine if values are all integers */
    n = numaGetCount(na);
    numaHasOnlyIntegers(na, &allints);
 
        /* Do simple integer binning if possible */
    if (allints && (maxval - minval < maxbins)) {
        imin = (l_int32)minval;
        imax = (l_int32)maxval;
        irange = imax - imin + 1;
        nah = numaCreate(irange);
        numaSetCount(nah, irange);  /* init */
        numaSetParameters(nah, minval, 1.0);
        for (i = 0; i < n; i++) {
            numaGetIValue(na, i, &ival);
            ibin = ival - imin;
            numaGetIValue(nah, ibin, &ival);
            numaSetValue(nah, ibin, ival + 1.0);
        }
 
        return nah;
    }
 
        /* Do float binning, even if the data is integers. */
    range = maxval - minval;
    binsize = range / (l_float32)maxbins;
    if (range == 0.0) {
        nah = numaCreate(1);
        numaSetParameters(nah, minval, binsize);
        numaAddNumber(nah, n);
        return nah;
    }
    nah = numaCreate(maxbins);
    numaSetCount(nah, maxbins);
    numaSetParameters(nah, minval, binsize);
    for (i = 0; i < n; i++) {
        numaGetFValue(na, i, &fval);
        ibin = (l_int32)((fval - minval) / binsize);
        ibin = L_MIN(ibin, maxbins - 1);  /* "edge" case; stay in bounds */
        numaGetIValue(nah, ibin, &ival);
        numaSetValue(nah, ibin, ival + 1.0);
    }
 
    return nah;
}
 

NUMA *
numaMakeHistogramClipped(NUMA      *na,
                         l_float32  binsize,
                         l_float32  maxsize)
{
l_int32    i, n, nbins, ival, ibin;
l_float32  val, maxval;
NUMA      *nad;
 
    if (!na)
        return (NUMA *)ERROR_PTR("na not defined", __func__, NULL);
    if (binsize <= 0.0)
        return (NUMA *)ERROR_PTR("binsize must be > 0.0", __func__, NULL);
    if (binsize > maxsize)
        binsize = maxsize;  /* just one bin */
 
    numaGetMax(na, &maxval, NULL);
    n = numaGetCount(na);
    maxsize = L_MIN(maxsize, maxval);
    nbins = (l_int32)(maxsize / binsize) + 1;
 
/*    lept_stderr("maxsize = %7.3f, nbins = %d\n", maxsize, nbins); */
 
    if ((nad = numaCreate(nbins)) == NULL)
        return (NUMA *)ERROR_PTR("nad not made", __func__, NULL);
    numaSetParameters(nad, 0.0, binsize);
    numaSetCount(nad, nbins);  /* interpret zeroes in bins as data */
    for (i = 0; i < n; i++) {
        numaGetFValue(na, i, &val);
        ibin = (l_int32)(val / binsize);
        if (ibin >= 0 && ibin < nbins) {
            numaGetIValue(nad, ibin, &ival);
            numaSetValue(nad, ibin, ival + 1.0);
        }
    }
 
    return nad;
}
 

NUMA *
numaRebinHistogram(NUMA    *nas,
                   l_int32  newsize)
{
l_int32    i, j, ns, nd, index, count, val;
l_float32  start, oldsize;
NUMA      *nad;
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
    if (newsize <= 1)
        return (NUMA *)ERROR_PTR("newsize must be > 1", __func__, NULL);
    if ((ns = numaGetCount(nas)) == 0)
        return (NUMA *)ERROR_PTR("no bins in nas", __func__, NULL);
 
    nd = (ns + newsize - 1) / newsize;
    if ((nad = numaCreate(nd)) == NULL)
        return (NUMA *)ERROR_PTR("nad not made", __func__, NULL);
    numaGetParameters(nad, &start, &oldsize);
    numaSetParameters(nad, start, oldsize * newsize);
 
    for (i = 0; i < nd; i++) {  /* new bins */
        count = 0;
        index = i * newsize;
        for (j = 0; j < newsize; j++) {
            if (index < ns) {
                numaGetIValue(nas, index, &val);
                count += val;
                index++;
            }
        }
        numaAddNumber(nad, count);
    }
 
    return nad;
}
 

NUMA *
numaNormalizeHistogram(NUMA      *nas,
                       l_float32  tsum)
{
l_int32    i, ns;
l_float32  sum, factor, fval;
NUMA      *nad;
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
    if (tsum <= 0.0)
        return (NUMA *)ERROR_PTR("tsum must be > 0.0", __func__, NULL);
    if ((ns = numaGetCount(nas)) == 0)
        return (NUMA *)ERROR_PTR("no bins in nas", __func__, NULL);
 
    numaGetSum(nas, &sum);
    factor = tsum / sum;
 
    if ((nad = numaCreate(ns)) == NULL)
        return (NUMA *)ERROR_PTR("nad not made", __func__, NULL);
    numaCopyParameters(nad, nas);
 
    for (i = 0; i < ns; i++) {
        numaGetFValue(nas, i, &fval);
        fval *= factor;
        numaAddNumber(nad, fval);
    }
 
    return nad;
}
 

l_ok
numaGetStatsUsingHistogram(NUMA       *na,
                           l_int32     maxbins,
                           l_float32  *pmin,
                           l_float32  *pmax,
                           l_float32  *pmean,
                           l_float32  *pvariance,
                           l_float32  *pmedian,
                           l_float32   rank,
                           l_float32  *prval,
                           NUMA      **phisto)
{
l_int32    i, n;
l_float32  minval, maxval, fval, mean, sum;
NUMA      *nah;
 
    if (pmin) *pmin = 0.0;
    if (pmax) *pmax = 0.0;
    if (pmean) *pmean = 0.0;
    if (pvariance) *pvariance = 0.0;
    if (pmedian) *pmedian = 0.0;
    if (prval) *prval = 0.0;
    if (phisto) *phisto = NULL;
    if (!na)
        return ERROR_INT("na not defined", __func__, 1);
    if ((n = numaGetCount(na)) == 0)
        return ERROR_INT("numa is empty", __func__, 1);
 
    numaGetMin(na, &minval, NULL);
    numaGetMax(na, &maxval, NULL);
    if (pmin) *pmin = minval;
    if (pmax) *pmax = maxval;
    if (pmean || pvariance) {
        sum = 0.0;
        for (i = 0; i < n; i++) {
            numaGetFValue(na, i, &fval);
            sum += fval;
        }
        mean = sum / (l_float32)n;
        if (pmean) *pmean = mean;
    }
    if (pvariance) {
        sum = 0.0;
        for (i = 0; i < n; i++) {
            numaGetFValue(na, i, &fval);
            sum += fval * fval;
        }
        *pvariance = sum / (l_float32)n - mean * mean;
    }
 
    if (!pmedian && !prval && !phisto)
        return 0;
 
    nah = numaMakeHistogramAuto(na, maxbins);
    if (pmedian)
        numaHistogramGetValFromRank(nah, 0.5, pmedian);
    if (prval)
        numaHistogramGetValFromRank(nah, rank, prval);
    if (phisto)
        *phisto = nah;
    else
        numaDestroy(&nah);
    return 0;
}
 

l_ok
numaGetHistogramStats(NUMA       *nahisto,
                      l_float32   startx,
                      l_float32   deltax,
                      l_float32  *pxmean,
                      l_float32  *pxmedian,
                      l_float32  *pxmode,
                      l_float32  *pxvariance)
{
    if (pxmean) *pxmean = 0.0;
    if (pxmedian) *pxmedian = 0.0;
    if (pxmode) *pxmode = 0.0;
    if (pxvariance) *pxvariance = 0.0;
    if (!nahisto)
        return ERROR_INT("nahisto not defined", __func__, 1);
 
    return numaGetHistogramStatsOnInterval(nahisto, startx, deltax, 0, -1,
                                           pxmean, pxmedian, pxmode,
                                           pxvariance);
}
 

l_ok
numaGetHistogramStatsOnInterval(NUMA       *nahisto,
                                l_float32   startx,
                                l_float32   deltax,
                                l_int32     ifirst,
                                l_int32     ilast,
                                l_float32  *pxmean,
                                l_float32  *pxmedian,
                                l_float32  *pxmode,
                                l_float32  *pxvariance)
{
l_int32    i, n, imax;
l_float32  sum, sumval, halfsum, moment, var, x, y, ymax;
 
    if (pxmean) *pxmean = 0.0;
    if (pxmedian) *pxmedian = 0.0;
    if (pxmode) *pxmode = 0.0;
    if (pxvariance) *pxvariance = 0.0;
    if (!nahisto)
        return ERROR_INT("nahisto not defined", __func__, 1);
    if (!pxmean && !pxmedian && !pxmode && !pxvariance)
        return ERROR_INT("nothing to compute", __func__, 1);
 
    n = numaGetCount(nahisto);
    ifirst = L_MAX(0, ifirst);
    if (ilast < 0) ilast = n - 1;
    if (ifirst >= n)
        return ERROR_INT("invalid ifirst", __func__, 1);
    if (ilast >= n) {
        L_WARNING("ilast = %d is beyond max index = %d; adjusting\n",
                  __func__, ilast, n - 1);
        ilast = n - 1;
    }
    if (ifirst > ilast)
        return ERROR_INT("ifirst > ilast", __func__, 1);
    for (sum = 0.0, moment = 0.0, var = 0.0, i = ifirst; i <= ilast ; i++) {
        x = startx + i * deltax;
        numaGetFValue(nahisto, i, &y);
        sum += y;
        moment += x * y;
        var += x * x * y;
    }
    if (sum == 0.0) {
        L_INFO("sum is 0\n", __func__);
        return 0;
    }
 
    if (pxmean)
        *pxmean = moment / sum;
    if (pxvariance)
        *pxvariance = var / sum - moment * moment / (sum * sum);
 
    if (pxmedian) {
        halfsum = sum / 2.0;
        for (sumval = 0.0, i = ifirst; i <= ilast; i++) {
            numaGetFValue(nahisto, i, &y);
            sumval += y;
            if (sumval >= halfsum) {
                *pxmedian = startx + i * deltax;
                break;
            }
        }
    }
 
    if (pxmode) {
        imax = -1;
        ymax = -1.0e10;
        for (i = ifirst; i <= ilast; i++) {
            numaGetFValue(nahisto, i, &y);
            if (y > ymax) {
                ymax = y;
                imax = i;
            }
        }
        *pxmode = startx + imax * deltax;
    }
 
    return 0;
}
 

l_ok
numaMakeRankFromHistogram(l_float32  startx,
                          l_float32  deltax,
                          NUMA      *nasy,
                          l_int32    npts,
                          NUMA     **pnax,
                          NUMA     **pnay)
{
l_int32    i, n;
l_float32  sum, fval;
NUMA      *nan, *nar;
 
    if (pnax) *pnax = NULL;
    if (!pnay)
        return ERROR_INT("&nay not defined", __func__, 1);
    *pnay = NULL;
    if (!nasy)
        return ERROR_INT("nasy not defined", __func__, 1);
    if ((n = numaGetCount(nasy)) == 0)
        return ERROR_INT("no bins in nas", __func__, 1);
 
        /* Normalize and generate the rank array corresponding to
         * the binned histogram. */
    nan = numaNormalizeHistogram(nasy, 1.0);
    nar = numaCreate(n + 1);  /* rank numa corresponding to nan */
    sum = 0.0;
    numaAddNumber(nar, sum);  /* first element is 0.0 */
    for (i = 0; i < n; i++) {
        numaGetFValue(nan, i, &fval);
        sum += fval;
        numaAddNumber(nar, sum);
    }
 
        /* Compute rank array on full range with specified
         * number of points and correspondence to x-values. */
    numaInterpolateEqxInterval(startx, deltax, nar, L_LINEAR_INTERP,
                               startx, startx + n * deltax, npts,
                               pnax, pnay);
    numaDestroy(&nan);
    numaDestroy(&nar);
    return 0;
}
 

l_ok
numaHistogramGetRankFromVal(NUMA       *na,
                            l_float32   rval,
                            l_float32  *prank)
{
l_int32    i, ibinval, n;
l_float32  startval, binsize, binval, maxval, fractval, total, sum, val;
 
    if (!prank)
        return ERROR_INT("prank not defined", __func__, 1);
    *prank = 0.0;
    if (!na)
        return ERROR_INT("na not defined", __func__, 1);
    numaGetParameters(na, &startval, &binsize);
    n = numaGetCount(na);
    if (rval < startval)
        return 0;
    maxval = startval + n * binsize;
    if (rval > maxval) {
        *prank = 1.0;
        return 0;
    }
 
    binval = (rval - startval) / binsize;
    ibinval = (l_int32)binval;
    if (ibinval >= n) {
        *prank = 1.0;
        return 0;
    }
    fractval = binval - (l_float32)ibinval;
 
    sum = 0.0;
    for (i = 0; i < ibinval; i++) {
        numaGetFValue(na, i, &val);
        sum += val;
    }
    numaGetFValue(na, ibinval, &val);
    sum += fractval * val;
    numaGetSum(na, &total);
    *prank = sum / total;
 
/*    lept_stderr("binval = %7.3f, rank = %7.3f\n", binval, *prank); */
 
    return 0;
}
 

l_ok
numaHistogramGetValFromRank(NUMA       *na,
                            l_float32   rank,
                            l_float32  *prval)
{
l_int32    i, n;
l_float32  startval, binsize, rankcount, total, sum, fract, val;
 
    if (!prval)
        return ERROR_INT("prval not defined", __func__, 1);
    *prval = 0.0;
    if (!na)
        return ERROR_INT("na not defined", __func__, 1);
    if (rank < 0.0) {
        L_WARNING("rank < 0; setting to 0.0\n", __func__);
        rank = 0.0;
    }
    if (rank > 1.0) {
        L_WARNING("rank > 1.0; setting to 1.0\n", __func__);
        rank = 1.0;
    }
 
    n = numaGetCount(na);
    numaGetParameters(na, &startval, &binsize);
    numaGetSum(na, &total);
    rankcount = rank * total;  /* count that corresponds to rank */
    sum = 0.0;
    for (i = 0; i < n; i++) {
        numaGetFValue(na, i, &val);
        if (sum + val >= rankcount)
            break;
        sum += val;
    }
    if (val <= 0.0)  /* can == 0 if rank == 0.0 */
        fract = 0.0;
    else  /* sum + fract * val = rankcount */
        fract = (rankcount - sum) / val;
 
    /* The use of the fraction of a bin allows a simple calculation
     * for the histogram value at the given rank. */
    *prval = startval + binsize * ((l_float32)i + fract);
 
/*    lept_stderr("rank = %7.3f, val = %7.3f\n", rank, *prval); */
 
    return 0;
}
 

l_ok
numaDiscretizeSortedInBins(NUMA    *na,
                           l_int32  nbins,
                           NUMA   **pnabinval)
{
NUMA      *nabinval;  /* average gray value in the bins */
NUMA      *naeach;
l_int32    i, ntot, bincount, binindex, binsize;
l_float32  sum, val, ave;
 
    if (!pnabinval)
        return ERROR_INT("&nabinval not defined", __func__, 1);
    *pnabinval = NULL;
    if (!na)
        return ERROR_INT("na not defined", __func__, 1);
    if (nbins < 2)
        return ERROR_INT("nbins must be > 1", __func__, 1);
 
        /* Get the number of items in each bin */
    ntot = numaGetCount(na);
    if ((naeach = numaGetUniformBinSizes(ntot, nbins)) == NULL)
        return ERROR_INT("naeach not made", __func__, 1);
 
        /* Get the average value in each bin */
    sum = 0.0;
    bincount = 0;
    binindex = 0;
    numaGetIValue(naeach, 0, &binsize);
    nabinval = numaCreate(nbins);
    for (i = 0; i < ntot; i++) {
        numaGetFValue(na, i, &val);
        bincount++;
        sum += val;
        if (bincount == binsize) {  /* add bin entry */
            ave = sum / binsize;
            numaAddNumber(nabinval, ave);
            sum = 0.0;
            bincount = 0;
            binindex++;
            if (binindex == nbins) break;
            numaGetIValue(naeach, binindex, &binsize);
        }
    }
    *pnabinval = nabinval;
 
    numaDestroy(&naeach);
    return 0;
}
 

l_ok
numaDiscretizeHistoInBins(NUMA    *na,
                          l_int32  nbins,
                          NUMA   **pnabinval,
                          NUMA   **pnarank)
{
NUMA      *nabinval;  /* average gray value in the bins */
NUMA      *naeach, *nan;
l_int32    i, j, nxvals, occup, count, bincount, binindex, binsize;
l_float32  sum, ave, ntot;
 
    if (pnarank) *pnarank = NULL;
    if (!pnabinval)
        return ERROR_INT("&nabinval not defined", __func__, 1);
    *pnabinval = NULL;
    if (!na)
        return ERROR_INT("na not defined", __func__, 1);
    if (nbins < 2)
        return ERROR_INT("nbins must be > 1", __func__, 1);
 
    nxvals = numaGetCount(na);
    numaGetSum(na, &ntot);
    occup = ntot / nxvals;
    if (occup < 1) L_INFO("average occupancy %d < 1\n", __func__, occup);
 
        /* Get the number of items in each bin */
    if ((naeach = numaGetUniformBinSizes(ntot, nbins)) == NULL)
        return ERROR_INT("naeach not made", __func__, 1);
 
        /* Get the average value in each bin */
    sum = 0.0;
    bincount = 0;
    binindex = 0;
    numaGetIValue(naeach, 0, &binsize);
    nabinval = numaCreate(nbins);
    for (i = 0; i < nxvals; i++) {
        numaGetIValue(na, i, &count);
        for (j = 0; j < count; j++) {
            bincount++;
            sum += i;
            if (bincount == binsize) {  /* add bin entry */
                ave = sum / binsize;
                numaAddNumber(nabinval, ave);
                sum = 0.0;
                bincount = 0;
                binindex++;
                if (binindex == nbins) break;
                numaGetIValue(naeach, binindex, &binsize);
            }
        }
        if (binindex == nbins) break;
    }
    *pnabinval = nabinval;
    if (binindex != nbins)
        L_ERROR("binindex = %d != nbins = %d\n", __func__, binindex, nbins);
 
        /* Get cumulative normalized histogram (rank[gray value]).
         * This is the partial sum operating on the normalized histogram. */
    if (pnarank) {
        nan = numaNormalizeHistogram(na, 1.0);
        *pnarank = numaGetPartialSums(nan);
        numaDestroy(&nan);
    }
    numaDestroy(&naeach);
    return 0;
}
 

l_ok
numaGetRankBinValues(NUMA    *na,
                     l_int32  nbins,
                     NUMA   **pnam)
{
NUMA      *na1;
l_int32    maxbins, type;
l_float32  maxval, delx;
 
    if (!pnam)
        return ERROR_INT("&pnam not defined", __func__, 1);
    *pnam = NULL;
    if (!na)
        return ERROR_INT("na not defined", __func__, 1);
    if (numaGetCount(na) == 0)
        return ERROR_INT("na is empty", __func__, 1);
    if (nbins < 2)
        return ERROR_INT("nbins must be > 1", __func__, 1);
 
        /* Choose between sorted array and a histogram.
         * If the input array is has a small number of numbers with
         * a large maximum, we will sort it.  At the other extreme, if
         * the array has many numbers with a small maximum, such as the
         * values of pixels in an 8 bpp grayscale image, generate a histogram.
         * If type comes back as L_BIN_SORT, use a histogram. */
    type = numaChooseSortType(na);
    if (type == L_SHELL_SORT) {  /* sort the array */
        L_INFO("sort the array: input size = %d\n", __func__, numaGetCount(na));
        na1 = numaSort(NULL, na, L_SORT_INCREASING);
        numaDiscretizeSortedInBins(na1, nbins, pnam);
        numaDestroy(&na1);
        return 0;
    }
 
        /* Make the histogram.  Assuming there are no negative values
         * in the array, if the max value in the array does not exceed
         * about 100000, the bin size for generating the histogram will
         * be 1; maxbins refers to the number of entries in the histogram. */
    L_INFO("use a histogram: input size = %d\n", __func__, numaGetCount(na));
    numaGetMax(na, &maxval, NULL);
    maxbins = L_MIN(100002, (l_int32)maxval + 2);
    na1 = numaMakeHistogram(na, maxbins, NULL, NULL);
 
        /* Warn if there is a scale change.  This shouldn't happen
         * unless the max value is above 100000.  */
    numaGetParameters(na1, NULL, &delx);
    if (delx > 1.0)
        L_WARNING("scale change: delx = %6.2f\n", __func__, delx);
 
        /* Rank bin the results */
    numaDiscretizeHistoInBins(na1, nbins, pnam, NULL);
    numaDestroy(&na1);
    return 0;
}
 

NUMA *
numaGetUniformBinSizes(l_int32  ntotal,
                       l_int32  nbins)
{
l_int32  i, start, end;
NUMA    *naeach;
 
    if (ntotal <= 0)
        return (NUMA *)ERROR_PTR("ntotal <= 0", __func__, NULL);
    if (nbins <= 0)
        return (NUMA *)ERROR_PTR("nbins <= 0", __func__, NULL);
 
    if ((naeach = numaCreate(nbins)) == NULL)
        return (NUMA *)ERROR_PTR("naeach not made", __func__, NULL);
    start = 0;
    for (i = 0; i < nbins; i++) {
        end = ntotal * (i + 1) / nbins;
        numaAddNumber(naeach, end - start);
        start = end;
    }
    return naeach;
}
 
 
/*----------------------------------------------------------------------*
 *                      Splitting a distribution                        *
 *----------------------------------------------------------------------*/
l_ok
numaSplitDistribution(NUMA       *na,
                      l_float32   scorefract,
                      l_int32    *psplitindex,
                      l_float32  *pave1,
                      l_float32  *pave2,
                      l_float32  *pnum1,
                      l_float32  *pnum2,
                      NUMA      **pnascore)
{
l_int32    i, n, bestsplit, minrange, maxrange, maxindex;
l_float32  ave1, ave2, ave1prev, ave2prev;
l_float32  num1, num2, num1prev, num2prev;
l_float32  val, minval, sum, fract1;
l_float32  norm, score, minscore, maxscore;
NUMA      *nascore, *naave1, *naave2, *nanum1, *nanum2;
 
    if (psplitindex) *psplitindex = 0;
    if (pave1) *pave1 = 0.0;
    if (pave2) *pave2 = 0.0;
    if (pnum1) *pnum1 = 0.0;
    if (pnum2) *pnum2 = 0.0;
    if (pnascore) *pnascore = NULL;
    if (!na)
        return ERROR_INT("na not defined", __func__, 1);
 
    n = numaGetCount(na);
    if (n <= 1)
        return ERROR_INT("n = 1 in histogram", __func__, 1);
    numaGetSum(na, &sum);
    if (sum <= 0.0)
        return ERROR_INT("sum <= 0.0", __func__, 1);
    norm = 4.0 / ((l_float32)(n - 1) * (n - 1));
    ave1prev = 0.0;
    numaGetHistogramStats(na, 0.0, 1.0, &ave2prev, NULL, NULL, NULL);
    num1prev = 0.0;
    num2prev = sum;
    maxindex = n / 2;  /* initialize with something */
 
        /* Split the histogram with [0 ... i] in the lower part
         * and [i+1 ... n-1] in upper part.  First, compute an otsu
         * score for each possible splitting.  */
    if ((nascore = numaCreate(n)) == NULL)
        return ERROR_INT("nascore not made", __func__, 1);
    naave1 = (pave1) ? numaCreate(n) : NULL;
    naave2 = (pave2) ? numaCreate(n) : NULL;
    nanum1 = (pnum1) ? numaCreate(n) : NULL;
    nanum2 = (pnum2) ? numaCreate(n) : NULL;
    maxscore = 0.0;
    for (i = 0; i < n; i++) {
        numaGetFValue(na, i, &val);
        num1 = num1prev + val;
        if (num1 == 0)
            ave1 = ave1prev;
        else
            ave1 = (num1prev * ave1prev + i * val) / num1;
        num2 = num2prev - val;
        if (num2 == 0)
            ave2 = ave2prev;
        else
            ave2 = (num2prev * ave2prev - i * val) / num2;
        fract1 = num1 / sum;
        score = norm * (fract1 * (1 - fract1)) * (ave2 - ave1) * (ave2 - ave1);
        numaAddNumber(nascore, score);
        if (pave1) numaAddNumber(naave1, ave1);
        if (pave2) numaAddNumber(naave2, ave2);
        if (pnum1) numaAddNumber(nanum1, num1);
        if (pnum2) numaAddNumber(nanum2, num2);
        if (score > maxscore) {
            maxscore = score;
            maxindex = i;
        }
        num1prev = num1;
        num2prev = num2;
        ave1prev = ave1;
        ave2prev = ave2;
    }
 
        /* Next, for all contiguous scores within a specified fraction
         * of the max, choose the split point as the value with the
         * minimum in the histogram. */
    minscore = (1. - scorefract) * maxscore;
    for (i = maxindex - 1; i >= 0; i--) {
        numaGetFValue(nascore, i, &val);
        if (val < minscore)
            break;
    }
    minrange = i + 1;
    for (i = maxindex + 1; i < n; i++) {
        numaGetFValue(nascore, i, &val);
        if (val < minscore)
            break;
    }
    maxrange = i - 1;
    numaGetFValue(na, minrange, &minval);
    bestsplit = minrange;
    for (i = minrange + 1; i <= maxrange; i++) {
        numaGetFValue(na, i, &val);
        if (val < minval) {
            minval = val;
            bestsplit = i;
        }
    }
 
        /* Add one to the bestsplit value to get the threshold value,
         * because when we take a threshold, as in pixThresholdToBinary(),
         * we always choose the set with values below the threshold. */
    bestsplit = L_MIN(255, bestsplit + 1);
 
    if (psplitindex) *psplitindex = bestsplit;
    if (pave1) numaGetFValue(naave1, bestsplit, pave1);
    if (pave2) numaGetFValue(naave2, bestsplit, pave2);
    if (pnum1) numaGetFValue(nanum1, bestsplit, pnum1);
    if (pnum2) numaGetFValue(nanum2, bestsplit, pnum2);
 
    if (pnascore) {  /* debug mode */
        lept_stderr("minrange = %d, maxrange = %d\n", minrange, maxrange);
        lept_stderr("minval = %10.0f\n", minval);
        gplotSimple1(nascore, GPLOT_PNG, "/tmp/lept/nascore",
                     "Score for split distribution");
        *pnascore = nascore;
    } else {
        numaDestroy(&nascore);
    }
 
    if (pave1) numaDestroy(&naave1);
    if (pave2) numaDestroy(&naave2);
    if (pnum1) numaDestroy(&nanum1);
    if (pnum2) numaDestroy(&nanum2);
    return 0;
}
 
 
/*----------------------------------------------------------------------*
 *                         Comparing histograms                         *
 *----------------------------------------------------------------------*/
l_ok
grayHistogramsToEMD(NUMAA  *naa1,
                    NUMAA  *naa2,
                    NUMA  **pnad)
{
l_int32     i, n, nt;
l_float32   dist;
NUMA       *na1, *na2, *nad;
 
    if (!pnad)
        return ERROR_INT("&nad not defined", __func__, 1);
    *pnad = NULL;
    if (!naa1 || !naa2)
        return ERROR_INT("na1 and na2 not both defined", __func__, 1);
    n = numaaGetCount(naa1);
    if (n != numaaGetCount(naa2))
        return ERROR_INT("naa1 and naa2 numa counts differ", __func__, 1);
    nt = numaaGetNumberCount(naa1);
    if (nt != numaaGetNumberCount(naa2))
        return ERROR_INT("naa1 and naa2 number counts differ", __func__, 1);
    if (256 * n != nt)  /* good enough check */
        return ERROR_INT("na sizes must be 256", __func__, 1);
 
    nad = numaCreate(n);
    *pnad = nad;
    for (i = 0; i < n; i++) {
        na1 = numaaGetNuma(naa1, i, L_CLONE);
        na2 = numaaGetNuma(naa2, i, L_CLONE);
        numaEarthMoverDistance(na1, na2, &dist);
        numaAddNumber(nad, dist / 255.);  /* normalize to [0.0 - 1.0] */
        numaDestroy(&na1);
        numaDestroy(&na2);
    }
    return 0;
}
 

l_ok
numaEarthMoverDistance(NUMA       *na1,
                       NUMA       *na2,
                       l_float32  *pdist)
{
l_int32     n, norm, i;
l_float32   sum1, sum2, diff, total;
l_float32  *array1, *array3;
NUMA       *na3;
 
    if (!pdist)
        return ERROR_INT("&dist not defined", __func__, 1);
    *pdist = 0.0;
    if (!na1 || !na2)
        return ERROR_INT("na1 and na2 not both defined", __func__, 1);
    n = numaGetCount(na1);
    if (n != numaGetCount(na2))
        return ERROR_INT("na1 and na2 have different size", __func__, 1);
 
        /* Generate na3; normalize to na1 if necessary */
    numaGetSum(na1, &sum1);
    numaGetSum(na2, &sum2);
    norm = (L_ABS(sum1 - sum2) < 0.00001 * L_ABS(sum1)) ? 1 : 0;
    if (!norm)
        na3 = numaTransform(na2, 0, sum1 / sum2);
    else
        na3 = numaCopy(na2);
    array1 = numaGetFArray(na1, L_NOCOPY);
    array3 = numaGetFArray(na3, L_NOCOPY);
 
        /* Move earth in n3 from array elements, to match n1 */
    total = 0;
    for (i = 1; i < n; i++) {
        diff = array1[i - 1] - array3[i - 1];
        array3[i] -= diff;
        total += L_ABS(diff);
    }
    *pdist = total / sum1;
 
    numaDestroy(&na3);
    return 0;
}
 

l_ok
grayInterHistogramStats(NUMAA   *naa,
                        l_int32  wc,
                        NUMA   **pnam,
                        NUMA   **pnams,
                        NUMA   **pnav,
                        NUMA   **pnarv)
{
l_int32      i, j, n, nn;
l_float32  **arrays;
l_float32    mean, var, rvar;
NUMA        *na1, *na2, *na3, *na4;
 
    if (pnam) *pnam = NULL;
    if (pnams) *pnams = NULL;
    if (pnav) *pnav = NULL;
    if (pnarv) *pnarv = NULL;
    if (!pnam && !pnams && !pnav && !pnarv)
        return ERROR_INT("nothing requested", __func__, 1);
    if (!naa)
        return ERROR_INT("naa not defined", __func__, 1);
    n = numaaGetCount(naa);
    for (i = 0; i < n; i++) {
        nn = numaaGetNumaCount(naa, i);
        if (nn != 256) {
            L_ERROR("%d numbers in numa[%d]\n", __func__, nn, i);
            return 1;
        }
    }
 
    if (pnam) *pnam = numaCreate(256);
    if (pnams) *pnams = numaCreate(256);
    if (pnav) *pnav = numaCreate(256);
    if (pnarv) *pnarv = numaCreate(256);
 
        /* First, use mean smoothing, normalize each histogram,
         * and save all results in a 2D matrix. */
    arrays = (l_float32 **)LEPT_CALLOC(n, sizeof(l_float32 *));
    for (i = 0; i < n; i++) {
        na1 = numaaGetNuma(naa, i, L_CLONE);
        na2 = numaWindowedMean(na1, wc);
        na3 = numaNormalizeHistogram(na2, 10000.);
        arrays[i] = numaGetFArray(na3, L_COPY);
        numaDestroy(&na1);
        numaDestroy(&na2);
        numaDestroy(&na3);
    }
 
        /* Get stats between histograms */
    for (j = 0; j < 256; j++) {
        na4 = numaCreate(n);
        for (i = 0; i < n; i++) {
            numaAddNumber(na4, arrays[i][j]);
        }
        numaSimpleStats(na4, 0, -1, &mean, &var, &rvar);
        if (pnam) numaAddNumber(*pnam, mean);
        if (pnams) numaAddNumber(*pnams, mean * mean);
        if (pnav) numaAddNumber(*pnav, var);
        if (pnarv) numaAddNumber(*pnarv, rvar);
        numaDestroy(&na4);
    }
 
    for (i = 0; i < n; i++)
        LEPT_FREE(arrays[i]);
    LEPT_FREE(arrays);
    return 0;
}
 
 
/*----------------------------------------------------------------------*
 *                             Extrema finding                          *
 *----------------------------------------------------------------------*/
NUMA *
numaFindPeaks(NUMA      *nas,
              l_int32    nmax,
              l_float32  fract1,
              l_float32  fract2)
{
l_int32    i, k, n, maxloc, lloc, rloc;
l_float32  fmaxval, sum, total, newtotal, val, lastval;
l_float32  peakfract;
NUMA      *na, *napeak;
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
    n = numaGetCount(nas);
    numaGetSum(nas, &total);
 
        /* We munge this copy */
    if ((na = numaCopy(nas)) == NULL)
        return (NUMA *)ERROR_PTR("na not made", __func__, NULL);
    if ((napeak = numaCreate(4 * nmax)) == NULL) {
        numaDestroy(&na);
        return (NUMA *)ERROR_PTR("napeak not made", __func__, NULL);
    }
 
    for (k = 0; k < nmax; k++) {
        numaGetSum(na, &newtotal);
        if (newtotal == 0.0)   /* sanity check */
            break;
        numaGetMax(na, &fmaxval, &maxloc);
        sum = fmaxval;
        lastval = fmaxval;
        lloc = 0;
        for (i = maxloc - 1; i >= 0; --i) {
            numaGetFValue(na, i, &val);
            if (val == 0.0) {
                lloc = i + 1;
                break;
            }
            if (val > fract1 * fmaxval) {
                sum += val;
                lastval = val;
                continue;
            }
            if (lastval - val > fract2 * lastval) {
                sum += val;
                lastval = val;
                continue;
            }
            lloc = i;
            break;
        }
        lastval = fmaxval;
        rloc = n - 1;
        for (i = maxloc + 1; i < n; ++i) {
            numaGetFValue(na, i, &val);
            if (val == 0.0) {
                rloc = i - 1;
                break;
            }
            if (val > fract1 * fmaxval) {
                sum += val;
                lastval = val;
                continue;
            }
            if (lastval - val > fract2 * lastval) {
                sum += val;
                lastval = val;
                continue;
            }
            rloc = i;
            break;
        }
        peakfract = sum / total;
        numaAddNumber(napeak, lloc);
        numaAddNumber(napeak, maxloc);
        numaAddNumber(napeak, rloc);
        numaAddNumber(napeak, peakfract);
 
        for (i = lloc; i <= rloc; i++)
            numaSetValue(na, i, 0.0);
    }
 
    numaDestroy(&na);
    return napeak;
}
 

NUMA *
numaFindExtrema(NUMA      *nas,
                l_float32  delta,
                NUMA     **pnav)
{
l_int32    i, n, found, loc, direction;
l_float32  startval, val, maxval, minval;
NUMA      *nav, *nad;
 
    if (pnav) *pnav = NULL;
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
    if (delta < 0.0)
        return (NUMA *)ERROR_PTR("delta < 0", __func__, NULL);
 
    n = numaGetCount(nas);
    nad = numaCreate(0);
    nav =  NULL;
    if (pnav) {
        nav = numaCreate(0);
        *pnav = nav;
    }
 
        /* We don't know if we'll find a peak or valley first,
         * but use the first element of nas as the reference point.
         * Break when we deviate by 'delta' from the first point. */
    numaGetFValue(nas, 0, &startval);
    found = FALSE;
    for (i = 1; i < n; i++) {
        numaGetFValue(nas, i, &val);
        if (L_ABS(val - startval) >= delta) {
            found = TRUE;
            break;
        }
    }
 
    if (!found)
        return nad;  /* it's empty */
 
        /* Are we looking for a peak or a valley? */
    if (val > startval) {  /* peak */
        direction = 1;
        maxval = val;
    } else {
        direction = -1;
        minval = val;
    }
    loc = i;
 
        /* Sweep through the rest of the array, recording alternating
         * peak/valley extrema. */
    for (i = i + 1; i < n; i++) {
        numaGetFValue(nas, i, &val);
        if (direction == 1 && val > maxval ) {  /* new local max */
            maxval = val;
            loc = i;
        } else if (direction == -1 && val < minval ) {  /* new local min */
            minval = val;
            loc = i;
        } else if (direction == 1 && (maxval - val >= delta)) {
            numaAddNumber(nad, loc);  /* save the current max location */
            if (nav) numaAddNumber(nav, maxval);
            direction = -1;  /* reverse: start looking for a min */
            minval = val;
            loc = i;  /* current min location */
        } else if (direction == -1 && (val - minval >= delta)) {
            numaAddNumber(nad, loc);  /* save the current min location */
            if (nav) numaAddNumber(nav, minval);
            direction = 1;  /* reverse: start looking for a max */
            maxval = val;
            loc = i;  /* current max location */
        }
    }
 
        /* Save the final extremum */
/*    numaAddNumber(nad, loc); */
    return nad;
}
 

l_ok
numaFindLocForThreshold(NUMA       *na,
                        l_int32     skip,
                        l_int32    *pthresh,
                        l_float32  *pfract)
{
l_int32     i, n, start, index, minloc, found;
l_float32   val, pval, jval, minval, maxval, sum, partsum;
l_float32  *fa;
 
    if (pfract) *pfract = 0.0;
    if (!pthresh)
        return ERROR_INT("&thresh not defined", __func__, 1);
    *pthresh = 0;
    if (!na)
        return ERROR_INT("na not defined", __func__, 1);
    if (skip <= 0) skip = 20;
 
        /* Test for constant value */
    numaGetMin(na, &minval, NULL);
    numaGetMax(na, &maxval, NULL);
    if (minval == maxval)
       return ERROR_INT("all array values are the same", __func__, 1);
 
        /* Look for the top of the first peak */
    n = numaGetCount(na);
    if (n < 256)
        L_WARNING("array size %d < 256\n", __func__, n);
    fa = numaGetFArray(na, L_NOCOPY);
    pval = fa[0];
    for (i = 1; i < n; i++) {
        val = fa[i];
        index = L_MIN(i + skip, n - 1);
        jval = fa[index];
        if (val < pval && jval < pval)  /* near the top if not there */
            break;
        pval = val;
    }
 
    if (i > n - 5)  /* just an increasing function */
       return ERROR_INT("top of first peak not found", __func__, 1);
 
        /* Look for the low point in the valley */
    found = FALSE;
    start = i;
    pval = fa[start];
    for (i = start + 1; i < n; i++) {
        val = fa[i];
        if (val <= pval) {  /* appears to be going down */
            pval = val;
        } else {  /* appears to be going up */
            index = L_MIN(i + skip, n - 1);
            jval = fa[index];  /* junp ahead by 'skip' */
            if (val > jval) {  /* still going down; jump ahead */
                pval = jval;
                i = index;
            } else {  /* really going up; passed the min */
                found = TRUE;
                break;
            }
        }
    }
    if (!found)
       return ERROR_INT("no minimum found", __func__, 1);
 
        /* Find the location of the minimum in the interval */
    minloc = index;  /* likely passed the min; look backward */
    minval = fa[index];
    for (i = index - 1; i > index - skip; i--) {
        if (fa[i] < minval) {
            minval = fa[i];
            minloc = i;
        }
    }
 
        /* Is the minimum very near the end of the array? */
    if (minloc > n - 10)
       return ERROR_INT("minimum at end of array; invalid", __func__, 1);
    *pthresh = minloc;
 
        /* Find the fraction under the first peak */
    if (pfract) {
        numaGetSumOnInterval(na, 0, minloc, &partsum);
        numaGetSum(na, &sum);
        if (sum > 0.0)
           *pfract = partsum / sum;
    }
    return 0;
}
 

l_ok
numaCountReversals(NUMA       *nas,
                   l_float32   minreversal,
                   l_int32    *pnr,
                   l_float32  *prd)
{
l_int32    i, n, nr, ival, binvals;
l_int32   *ia;
l_float32  fval, delx, len;
NUMA      *nat;
 
    if (pnr) *pnr = 0;
    if (prd) *prd = 0.0;
    if (!pnr && !prd)
        return ERROR_INT("neither &nr nor &rd are defined", __func__, 1);
    if (!nas)
        return ERROR_INT("nas not defined", __func__, 1);
    if ((n = numaGetCount(nas)) == 0) {
        L_INFO("nas is empty\n", __func__);
        return 0;
    }
    if (minreversal < 0.0)
        return ERROR_INT("minreversal < 0", __func__, 1);
 
        /* Decide if the only values are 0 and 1 */
    binvals = TRUE;
    for (i = 0; i < n; i++) {
        numaGetFValue(nas, i, &fval);
        if (fval != 0.0 && fval != 1.0) {
            binvals = FALSE;
            break;
        }
    }
 
    nr = 0;
    if (binvals) {
        if (minreversal > 1.0) {
            L_WARNING("binary values but minreversal > 1\n", __func__);
        } else {
            ia = numaGetIArray(nas);
            ival = ia[0];
            for (i = 1; i < n; i++) {
                if (ia[i] != ival) {
                    nr++;
                    ival = ia[i];
                }
            }
            LEPT_FREE(ia);
        }
    } else {
        nat = numaFindExtrema(nas, minreversal, NULL);
        nr = numaGetCount(nat);
        numaDestroy(&nat);
    }
    if (pnr) *pnr = nr;
    if (prd) {
        numaGetParameters(nas, NULL, &delx);
        len = delx * n;
        *prd = (l_float32)nr / len;
    }
 
    return 0;
}
 
 
/*----------------------------------------------------------------------*
 *                Threshold crossings and frequency analysis            *
 *----------------------------------------------------------------------*/
l_ok
numaSelectCrossingThreshold(NUMA       *nax,
                            NUMA       *nay,
                            l_float32   estthresh,
                            l_float32  *pbestthresh)
{
l_int32    i, inrun, istart, iend, maxstart, maxend, runlen, maxrunlen;
l_int32    val, maxval, nmax, count;
l_float32  thresh, fmaxval, fmodeval;
NUMA      *nat, *nac;
 
    if (!pbestthresh)
        return ERROR_INT("&bestthresh not defined", __func__, 1);
    *pbestthresh = 0.0;
    if (!nay)
        return ERROR_INT("nay not defined", __func__, 1);
    if (numaGetCount(nay) < 2) {
        L_WARNING("nay count < 2; no threshold crossing\n", __func__);
        return 1;
    }
 
        /* Compute the number of crossings for different thresholds */
    nat = numaCreate(41);
    for (i = 0; i < 41; i++) {
        thresh = estthresh - 80.0 + 4.0 * i;
        nac = numaCrossingsByThreshold(nax, nay, thresh);
        numaAddNumber(nat, numaGetCount(nac));
        numaDestroy(&nac);
    }
 
        /* Find the center of the plateau of max crossings, which
         * extends from thresh[istart] to thresh[iend]. */
    numaGetMax(nat, &fmaxval, NULL);
    maxval = (l_int32)fmaxval;
    nmax = 0;
    for (i = 0; i < 41; i++) {
        numaGetIValue(nat, i, &val);
        if (val == maxval)
            nmax++;
    }
    if (nmax < 3) {  /* likely accidental max; try the mode */
        numaGetMode(nat, &fmodeval, &count);
        if (count > nmax && fmodeval > 0.5 * fmaxval)
            maxval = (l_int32)fmodeval;  /* use the mode */
    }
 
    inrun = FALSE;
    iend = 40;
    maxrunlen = 0, maxstart = 0, maxend = 0;
    for (i = 0; i < 41; i++) {
        numaGetIValue(nat, i, &val);
        if (val == maxval) {
            if (!inrun) {
                istart = i;
                inrun = TRUE;
            }
            continue;
        }
        if (inrun && (val != maxval)) {
            iend = i - 1;
            runlen = iend - istart + 1;
            inrun = FALSE;
            if (runlen > maxrunlen) {
                maxstart = istart;
                maxend = iend;
                maxrunlen = runlen;
            }
        }
    }
    if (inrun) {
        runlen = i - istart;
        if (runlen > maxrunlen) {
            maxstart = istart;
            maxend = i - 1;
            maxrunlen = runlen;
        }
    }
 
    *pbestthresh = estthresh - 80.0 + 2.0 * (l_float32)(maxstart + maxend);
 
#if  DEBUG_CROSSINGS
    lept_stderr("\nCrossings attain a maximum at %d thresholds, between:\n"
                "  thresh[%d] = %5.1f and thresh[%d] = %5.1f\n",
                nmax, maxstart, estthresh - 80.0 + 4.0 * maxstart,
                maxend, estthresh - 80.0 + 4.0 * maxend);
    lept_stderr("The best choice: %5.1f\n", *pbestthresh);
    lept_stderr("Number of crossings at the 41 thresholds:");
    numaWriteStderr(nat);
#endif  /* DEBUG_CROSSINGS */
 
    numaDestroy(&nat);
    return 0;
}
 

NUMA *
numaCrossingsByThreshold(NUMA      *nax,
                         NUMA      *nay,
                         l_float32  thresh)
{
l_int32    i, n;
l_float32  startx, delx;
l_float32  xval1, xval2, yval1, yval2, delta1, delta2, crossval, fract;
NUMA      *nad;
 
    if (!nay)
        return (NUMA *)ERROR_PTR("nay not defined", __func__, NULL);
    n = numaGetCount(nay);
 
    if (nax && (numaGetCount(nax) != n))
        return (NUMA *)ERROR_PTR("nax and nay sizes differ", __func__, NULL);
 
    nad = numaCreate(0);
    if (n < 2) return nad;
    numaGetFValue(nay, 0, &yval1);
    numaGetParameters(nay, &startx, &delx);
    if (nax)
        numaGetFValue(nax, 0, &xval1);
    else
        xval1 = startx;
    for (i = 1; i < n; i++) {
        numaGetFValue(nay, i, &yval2);
        if (nax)
            numaGetFValue(nax, i, &xval2);
        else
            xval2 = startx + i * delx;
        delta1 = yval1 - thresh;
        delta2 = yval2 - thresh;
        if (delta1 == 0.0) {
            numaAddNumber(nad, xval1);
        } else if (delta2 == 0.0) {
            numaAddNumber(nad, xval2);
        } else if (delta1 * delta2 < 0.0) {  /* crossing */
            fract = L_ABS(delta1) / L_ABS(yval1 - yval2);
            crossval = xval1 + fract * (xval2 - xval1);
            numaAddNumber(nad, crossval);
        }
        xval1 = xval2;
        yval1 = yval2;
    }
 
    return nad;
}
 

NUMA *
numaCrossingsByPeaks(NUMA      *nax,
                     NUMA      *nay,
                     l_float32  delta)
{
l_int32    i, j, n, np, previndex, curindex;
l_float32  startx, delx;
l_float32  xval1, xval2, yval1, yval2, delta1, delta2;
l_float32  prevval, curval, thresh, crossval, fract;
NUMA      *nap, *nad;
 
    if (!nay)
        return (NUMA *)ERROR_PTR("nay not defined", __func__, NULL);
 
    n = numaGetCount(nay);
    if (nax && (numaGetCount(nax) != n))
        return (NUMA *)ERROR_PTR("nax and nay sizes differ", __func__, NULL);
 
        /* Find the extrema.  Also add last point in nay to get
         * the last transition (from the last peak to the end).
         * The number of crossings is 1 more than the number of extrema. */
    nap = numaFindExtrema(nay, delta, NULL);
    numaAddNumber(nap, n - 1);
    np = numaGetCount(nap);
    L_INFO("Number of crossings: %d\n", __func__, np);
 
        /* Do all computation in index units of nax or the delx of nay */
    nad = numaCreate(np);  /* output crossing locations, in nax units */
    previndex = 0;  /* prime the search with 1st point */
    numaGetFValue(nay, 0, &prevval);  /* prime the search with 1st point */
    numaGetParameters(nay, &startx, &delx);
    for (i = 0; i < np; i++) {
        numaGetIValue(nap, i, &curindex);
        numaGetFValue(nay, curindex, &curval);
        thresh = (prevval + curval) / 2.0;
        if (nax)
            numaGetFValue(nax, previndex, &xval1);
        else
            xval1 = startx + previndex * delx;
        numaGetFValue(nay, previndex, &yval1);
        for (j = previndex + 1; j <= curindex; j++) {
            if (nax)
                numaGetFValue(nax, j, &xval2);
            else
                xval2 = startx + j * delx;
            numaGetFValue(nay, j, &yval2);
            delta1 = yval1 - thresh;
            delta2 = yval2 - thresh;
            if (delta1 == 0.0) {
                numaAddNumber(nad, xval1);
                break;
            } else if (delta2 == 0.0) {
                numaAddNumber(nad, xval2);
                break;
            } else if (delta1 * delta2 < 0.0) {  /* crossing */
                fract = L_ABS(delta1) / L_ABS(yval1 - yval2);
                crossval = xval1 + fract * (xval2 - xval1);
                numaAddNumber(nad, crossval);
                break;
            }
            xval1 = xval2;
            yval1 = yval2;
        }
        previndex = curindex;
        prevval = curval;
    }
 
    numaDestroy(&nap);
    return nad;
}
 

l_ok
numaEvalBestHaarParameters(NUMA       *nas,
                           l_float32   relweight,
                           l_int32     nwidth,
                           l_int32     nshift,
                           l_float32   minwidth,
                           l_float32   maxwidth,
                           l_float32  *pbestwidth,
                           l_float32  *pbestshift,
                           l_float32  *pbestscore)
{
l_int32    i, j;
l_float32  delwidth, delshift, width, shift, score;
l_float32  bestwidth, bestshift, bestscore;
 
    if (pbestscore) *pbestscore = 0.0;
    if (pbestwidth) *pbestwidth = 0.0;
    if (pbestshift) *pbestshift = 0.0;
    if (!pbestwidth || !pbestshift)
        return ERROR_INT("&bestwidth and &bestshift not defined", __func__, 1);
    if (!nas)
        return ERROR_INT("nas not defined", __func__, 1);
 
    bestscore = bestwidth = bestshift = 0.0;
    delwidth = (maxwidth - minwidth) / (nwidth - 1.0);
    for (i = 0; i < nwidth; i++) {
        width = minwidth + delwidth * i;
        delshift = width / (l_float32)(nshift);
        for (j = 0; j < nshift; j++) {
            shift = j * delshift;
            numaEvalHaarSum(nas, width, shift, relweight, &score);
            if (score > bestscore) {
                bestscore = score;
                bestwidth = width;
                bestshift = shift;
#if  DEBUG_FREQUENCY
                lept_stderr("width = %7.3f, shift = %7.3f, score = %7.3f\n",
                            width, shift, score);
#endif  /* DEBUG_FREQUENCY */
            }
        }
    }
 
    *pbestwidth = bestwidth;
    *pbestshift = bestshift;
    if (pbestscore)
        *pbestscore = bestscore;
    return 0;
}
 

l_ok
numaEvalHaarSum(NUMA       *nas,
                l_float32   width,
                l_float32   shift,
                l_float32   relweight,
                l_float32  *pscore)
{
l_int32    i, n, nsamp, index;
l_float32  score, weight, val;
 
    if (!pscore)
        return ERROR_INT("&score not defined", __func__, 1);
    *pscore = 0.0;
    if (!nas)
        return ERROR_INT("nas not defined", __func__, 1);
    if ((n = numaGetCount(nas)) < 2 * width)
        return ERROR_INT("nas size too small", __func__, 1);
 
    score = 0.0;
    nsamp = (l_int32)((n - shift) / width);
    for (i = 0; i < nsamp; i++) {
        index = (l_int32)(shift + i * width);
        weight = (i % 2) ? 1.0 : -1.0 * relweight;
        numaGetFValue(nas, index, &val);
        score += weight * val;
    }
 
    *pscore = 2.0 * width * score / (l_float32)n;
    return 0;
}
 
 
/*----------------------------------------------------------------------*
 *            Generating numbers in a range under constraints           *
 *----------------------------------------------------------------------*/
NUMA *
genConstrainedNumaInRange(l_int32  first,
                          l_int32  last,
                          l_int32  nmax,
                          l_int32  use_pairs)
{
l_int32    i, nsets, val;
l_float32  delta;
NUMA      *na;
 
    first = L_MAX(0, first);
    if (last < first)
        return (NUMA *)ERROR_PTR("last < first!", __func__, NULL);
    if (nmax < 1)
        return (NUMA *)ERROR_PTR("nmax < 1!", __func__, NULL);
 
    nsets = L_MIN(nmax, last - first + 1);
    if (use_pairs == 1)
        nsets = nsets / 2;
    if (nsets == 0)
        return (NUMA *)ERROR_PTR("nsets == 0", __func__, NULL);
 
        /* Select delta so that selection covers the full range if possible */
    if (nsets == 1) {
        delta = 0.0;
    } else {
        if (use_pairs == 0)
            delta = (l_float32)(last - first) / (nsets - 1);
        else
            delta = (l_float32)(last - first - 1) / (nsets - 1);
    }
 
    na = numaCreate(nsets);
    for (i = 0; i < nsets; i++) {
        val = (l_int32)(first + i * delta + 0.5);
        numaAddNumber(na, val);
        if (use_pairs == 1)
            numaAddNumber(na, val + 1);
    }
 
    return na;
}