numafunc1.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"
 
/*----------------------------------------------------------------------*
 *                Arithmetic and logical ops on Numas                   *
 *----------------------------------------------------------------------*/
NUMA *
numaArithOp(NUMA    *nad,
            NUMA    *na1,
            NUMA    *na2,
            l_int32  op)
{
l_int32    i, n;
l_float32  val1, val2;
 
    PROCNAME("numaArithOp");
 
    if (!na1 || !na2)
        return (NUMA *)ERROR_PTR("na1, na2 not both defined", procName, nad);
    n = numaGetCount(na1);
    if (n != numaGetCount(na2))
        return (NUMA *)ERROR_PTR("na1, na2 sizes differ", procName, nad);
    if (nad && nad != na1)
        return (NUMA *)ERROR_PTR("nad defined but not in-place", procName, nad);
    if (op != L_ARITH_ADD && op != L_ARITH_SUBTRACT &&
        op != L_ARITH_MULTIPLY && op != L_ARITH_DIVIDE)
        return (NUMA *)ERROR_PTR("invalid op", procName, nad);
    if (op == L_ARITH_DIVIDE) {
        for (i = 0; i < n; i++) {
            numaGetFValue(na2, i, &val2);
            if (val2 == 0.0)
                return (NUMA *)ERROR_PTR("na2 has 0 element", procName, nad);
        }
    }
 
        /* If nad is not identical to na1, make it an identical copy */
    if (!nad)
        nad = numaCopy(na1);
 
    for (i = 0; i < n; i++) {
        numaGetFValue(nad, i, &val1);
        numaGetFValue(na2, i, &val2);
        switch (op) {
        case L_ARITH_ADD:
            numaSetValue(nad, i, val1 + val2);
            break;
        case L_ARITH_SUBTRACT:
            numaSetValue(nad, i, val1 - val2);
            break;
        case L_ARITH_MULTIPLY:
            numaSetValue(nad, i, val1 * val2);
            break;
        case L_ARITH_DIVIDE:
            numaSetValue(nad, i, val1 / val2);
            break;
        default:
            lept_stderr(" Unknown arith op: %d\n", op);
            return nad;
        }
    }
 
    return nad;
}
 
 
NUMA *
numaLogicalOp(NUMA    *nad,
              NUMA    *na1,
              NUMA    *na2,
              l_int32  op)
{
l_int32  i, n, val1, val2, val;
 
    PROCNAME("numaLogicalOp");
 
    if (!na1 || !na2)
        return (NUMA *)ERROR_PTR("na1, na2 not both defined", procName, nad);
    n = numaGetCount(na1);
    if (n != numaGetCount(na2))
        return (NUMA *)ERROR_PTR("na1, na2 sizes differ", procName, nad);
    if (nad && nad != na1)
        return (NUMA *)ERROR_PTR("nad defined; not in-place", procName, nad);
    if (op != L_UNION && op != L_INTERSECTION &&
        op != L_SUBTRACTION && op != L_EXCLUSIVE_OR)
        return (NUMA *)ERROR_PTR("invalid op", procName, nad);
 
        /* If nad is not identical to na1, make it an identical copy */
    if (!nad)
        nad = numaCopy(na1);
 
    for (i = 0; i < n; i++) {
        numaGetIValue(nad, i, &val1);
        numaGetIValue(na2, i, &val2);
        val1 = (val1 == 0) ? 0 : 1;
        val2 = (val2 == 0) ? 0 : 1;
        switch (op) {
        case L_UNION:
            val = (val1 || val2) ? 1 : 0;
            numaSetValue(nad, i, val);
            break;
        case L_INTERSECTION:
            val = (val1 && val2) ? 1 : 0;
            numaSetValue(nad, i, val);
            break;
        case L_SUBTRACTION:
            val = (val1 && !val2) ? 1 : 0;
            numaSetValue(nad, i, val);
            break;
        case L_EXCLUSIVE_OR:
            val = (val1 != val2) ? 1 : 0;
            numaSetValue(nad, i, val);
            break;
        default:
            lept_stderr(" Unknown logical op: %d\n", op);
            return nad;
        }
    }
 
    return nad;
}
 
 
NUMA *
numaInvert(NUMA  *nad,
           NUMA  *nas)
{
l_int32  i, n, val;
 
    PROCNAME("numaInvert");
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", procName, nad);
    if (nad && nad != nas)
        return (NUMA *)ERROR_PTR("nad defined; not in-place", procName, nad);
 
    if (!nad)
        nad = numaCopy(nas);
    n = numaGetCount(nad);
    for (i = 0; i < n; i++) {
        numaGetIValue(nad, i, &val);
        if (!val)
            val = 1;
        else
            val = 0;
        numaSetValue(nad, i, val);
    }
    return nad;
}
 
 
l_int32
numaSimilar(NUMA      *na1,
            NUMA      *na2,
            l_float32  maxdiff,
            l_int32   *psimilar)
{
l_int32    i, n;
l_float32  val1, val2;
 
    PROCNAME("numaSimilar");
 
    if (!psimilar)
        return ERROR_INT("&similar not defined", procName, 1);
    *psimilar = 0;
    if (!na1 || !na2)
        return ERROR_INT("na1 and na2 not both defined", procName, 1);
    maxdiff = L_ABS(maxdiff);
 
    n = numaGetCount(na1);
    if (n != numaGetCount(na2)) return 0;
 
    for (i = 0; i < n; i++) {
        numaGetFValue(na1, i, &val1);
        numaGetFValue(na2, i, &val2);
        if (L_ABS(val1 - val2) > maxdiff) return 0;
    }
 
    *psimilar = 1;
    return 0;
}
 
 
l_ok
numaAddToNumber(NUMA      *na,
                l_int32    index,
                l_float32  val)
{
l_int32  n;
 
    PROCNAME("numaAddToNumber");
 
    if (!na)
        return ERROR_INT("na not defined", procName, 1);
    if ((n = numaGetCount(na)) == 0)
        return ERROR_INT("na is empty", procName, 1);
    if (index < 0 || index >= n) {
        L_ERROR("index %d not in [0,...,%d]\n", procName, index, n - 1);
        return 1;
    }
 
    na->array[index] += val;
    return 0;
}
 
 
/*----------------------------------------------------------------------*
 *                         Simple extractions                           *
 *----------------------------------------------------------------------*/
l_ok
numaGetMin(NUMA       *na,
           l_float32  *pminval,
           l_int32    *piminloc)
{
l_int32    i, n, iminloc;
l_float32  val, minval;
 
    PROCNAME("numaGetMin");
 
    if (!pminval && !piminloc)
        return ERROR_INT("nothing to do", procName, 1);
    if (pminval) *pminval = 0.0;
    if (piminloc) *piminloc = 0;
    if (!na)
        return ERROR_INT("na not defined", procName, 1);
    if ((n = numaGetCount(na)) == 0)
        return ERROR_INT("na is empty", procName, 1);
 
    minval = +1000000000.;
    iminloc = 0;
    for (i = 0; i < n; i++) {
        numaGetFValue(na, i, &val);
        if (val < minval) {
            minval = val;
            iminloc = i;
        }
    }
 
    if (pminval) *pminval = minval;
    if (piminloc) *piminloc = iminloc;
    return 0;
}
 
 
l_ok
numaGetMax(NUMA       *na,
           l_float32  *pmaxval,
           l_int32    *pimaxloc)
{
l_int32    i, n, imaxloc;
l_float32  val, maxval;
 
    PROCNAME("numaGetMax");
 
    if (!pmaxval && !pimaxloc)
        return ERROR_INT("nothing to do", procName, 1);
    if (pmaxval) *pmaxval = 0.0;
    if (pimaxloc) *pimaxloc = 0;
    if (!na)
        return ERROR_INT("na not defined", procName, 1);
    if ((n = numaGetCount(na)) == 0)
        return ERROR_INT("na is empty", procName, 1);
 
    maxval = -1000000000.;
    imaxloc = 0;
    for (i = 0; i < n; i++) {
        numaGetFValue(na, i, &val);
        if (val > maxval) {
            maxval = val;
            imaxloc = i;
        }
    }
 
    if (pmaxval) *pmaxval = maxval;
    if (pimaxloc) *pimaxloc = imaxloc;
    return 0;
}
 
 
l_ok
numaGetSum(NUMA       *na,
           l_float32  *psum)
{
l_int32    i, n;
l_float32  val, sum;
 
    PROCNAME("numaGetSum");
 
    if (!psum)
        return ERROR_INT("&sum not defined", procName, 1);
    *psum = 0;
    if (!na)
        return ERROR_INT("na not defined", procName, 1);
 
    if ((n = numaGetCount(na)) == 0)
        return ERROR_INT("na is empty", procName, 1);
    sum = 0.0;
    for (i = 0; i < n; i++) {
        numaGetFValue(na, i, &val);
        sum += val;
    }
    *psum = sum;
    return 0;
}
 
 
NUMA *
numaGetPartialSums(NUMA  *na)
{
l_int32    i, n;
l_float32  val, sum;
NUMA      *nasum;
 
    PROCNAME("numaGetPartialSums");
 
    if (!na)
        return (NUMA *)ERROR_PTR("na not defined", procName, NULL);
 
    if ((n = numaGetCount(na)) == 0)
        L_WARNING("na is empty\n", procName);
    nasum = numaCreate(n);
    sum = 0.0;
    for (i = 0; i < n; i++) {
        numaGetFValue(na, i, &val);
        sum += val;
        numaAddNumber(nasum, sum);
    }
    return nasum;
}
 
 
l_ok
numaGetSumOnInterval(NUMA       *na,
                     l_int32     first,
                     l_int32     last,
                     l_float32  *psum)
{
l_int32    i, n;
l_float32  val, sum;
 
    PROCNAME("numaGetSumOnInterval");
 
    if (!psum)
        return ERROR_INT("&sum not defined", procName, 1);
    *psum = 0.0;
    if (!na)
        return ERROR_INT("na not defined", procName, 1);
 
    sum = 0.0;
    if ((n = numaGetCount(na)) == 0)
        return ERROR_INT("na is empty", procName, 1);
    if (first < 0) first = 0;
    if (first >= n || last < -1)  /* not an error */
        return 0;
    if (last == -1)
        last = n - 1;
    last = L_MIN(last, n - 1);
 
    for (i = first; i <= last; i++) {
        numaGetFValue(na, i, &val);
        sum += val;
    }
    *psum = sum;
    return 0;
}
 
 
l_ok
numaHasOnlyIntegers(NUMA     *na,
                    l_int32  *pallints)
{
l_int32    i, n;
l_float32  val;
 
    PROCNAME("numaHasOnlyIntegers");
 
    if (!pallints)
        return ERROR_INT("&allints not defined", procName, 1);
    *pallints = TRUE;
    if (!na)
        return ERROR_INT("na not defined", procName, 1);
 
    if ((n = numaGetCount(na)) == 0)
        return ERROR_INT("na is empty", procName, 1);
    for (i = 0; i < n; i ++) {
        numaGetFValue(na, i, &val);
        if (val != (l_int32)val) {
            *pallints = FALSE;
            return 0;
        }
    }
    return 0;
}
 
 
l_ok
numaGetMean(NUMA       *na,
            l_float32  *pave)
{
l_int32    n;
l_float32  sum;
 
    PROCNAME("numaGetMean");
 
    if (!pave)
        return ERROR_INT("&ave not defined", procName, 1);
    *pave = 0;
    if (!na)
        return ERROR_INT("na not defined", procName, 1);
    if ((n = numaGetCount(na)) == 0)
        return ERROR_INT("na is empty", procName, 1);
 
    numaGetSum(na, &sum);
    *pave = sum / n;
    return 0;
}
 
l_ok
numaGetMeanAbsval(NUMA       *na,
                  l_float32  *paveabs)
{
l_int32  n;
NUMA    *na1;
 
    PROCNAME("numaGetMeanAbsval");
 
    if (!paveabs)
        return ERROR_INT("&aveabs not defined", procName, 1);
    *paveabs = 0;
    if (!na)
        return ERROR_INT("na not defined", procName, 1);
    if ((n = numaGetCount(na)) == 0)
        return ERROR_INT("na is empty", procName, 1);
 
    na1 = numaMakeAbsval(NULL, na);
    numaGetMean(na1, paveabs);
    numaDestroy(&na1);
    return 0;
}
 
 
NUMA *
numaSubsample(NUMA    *nas,
              l_int32  subfactor)
{
l_int32    i, n;
l_float32  val;
NUMA      *nad;
 
    PROCNAME("numaSubsample");
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", procName, NULL);
    if (subfactor < 1)
        return (NUMA *)ERROR_PTR("subfactor < 1", procName, NULL);
 
    nad = numaCreate(0);
    if ((n = numaGetCount(nas)) == 0)
        L_WARNING("nas is empty\n", procName);
    for (i = 0; i < n; i++) {
        if (i % subfactor != 0) continue;
        numaGetFValue(nas, i, &val);
        numaAddNumber(nad, val);
    }
 
    return nad;
}
 
 
NUMA *
numaMakeDelta(NUMA  *nas)
{
l_int32    i, n;
l_float32  prev, cur;
NUMA      *nad;
 
    PROCNAME("numaMakeDelta");
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", procName, NULL);
    if ((n = numaGetCount(nas)) < 2) {
        L_WARNING("n < 2; returning empty numa\n", procName);
        return numaCreate(1);
    }
 
    nad = numaCreate(n - 1);
    numaGetFValue(nas, 0, &prev);
    for (i = 1; i < n; i++) {
        numaGetFValue(nas, i, &cur);
        numaAddNumber(nad, cur - prev);
        prev = cur;
    }
    return nad;
}
 
 
NUMA *
numaMakeSequence(l_float32  startval,
                 l_float32  increment,
                 l_int32    size)
{
l_int32    i;
l_float32  val;
NUMA      *na;
 
    PROCNAME("numaMakeSequence");
 
    if ((na = numaCreate(size)) == NULL)
        return (NUMA *)ERROR_PTR("na not made", procName, NULL);
 
    for (i = 0; i < size; i++) {
        val = startval + i * increment;
        numaAddNumber(na, val);
    }
    return na;
}
 
 
NUMA *
numaMakeConstant(l_float32  val,
                 l_int32    size)
{
    return numaMakeSequence(val, 0.0, size);
}
 
 
NUMA *
numaMakeAbsval(NUMA  *nad,
               NUMA  *nas)
{
l_int32    i, n;
l_float32  val;
 
    PROCNAME("numaMakeAbsval");
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", procName, NULL);
    if (nad && nad != nas)
        return (NUMA *)ERROR_PTR("nad and not in-place", procName, NULL);
 
    if (!nad)
        nad = numaCopy(nas);
    n = numaGetCount(nad);
    for (i = 0; i < n; i++) {
        val = nad->array[i];
        nad->array[i] = L_ABS(val);
    }
 
    return nad;
}
 
 
NUMA *
numaAddBorder(NUMA      *nas,
              l_int32    left,
              l_int32    right,
              l_float32  val)
{
l_int32     i, n, len;
l_float32   startx, delx;
l_float32  *fas, *fad;
NUMA       *nad;
 
    PROCNAME("numaAddBorder");
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", procName, NULL);
    if (left < 0) left = 0;
    if (right < 0) right = 0;
    if (left == 0 && right == 0)
        return numaCopy(nas);
 
    n = numaGetCount(nas);
    len = n + left + right;
    nad = numaMakeConstant(val, len);
    numaGetParameters(nas, &startx, &delx);
    numaSetParameters(nad, startx - delx * left, delx);
    fas = numaGetFArray(nas, L_NOCOPY);
    fad = numaGetFArray(nad, L_NOCOPY);
    for (i = 0; i < n; i++)
        fad[left + i] = fas[i];
 
    return nad;
}
 
 
NUMA *
numaAddSpecifiedBorder(NUMA    *nas,
                       l_int32  left,
                       l_int32  right,
                       l_int32  type)
{
l_int32     i, n;
l_float32  *fa;
NUMA       *nad;
 
    PROCNAME("numaAddSpecifiedBorder");
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", procName, NULL);
    if (left < 0) left = 0;
    if (right < 0) right = 0;
    if (left == 0 && right == 0)
        return numaCopy(nas);
    if (type != L_CONTINUED_BORDER && type != L_MIRRORED_BORDER)
        return (NUMA *)ERROR_PTR("invalid type", procName, NULL);
    n = numaGetCount(nas);
    if (type == L_MIRRORED_BORDER && (left > n || right > n))
        return (NUMA *)ERROR_PTR("border too large", procName, NULL);
 
    nad = numaAddBorder(nas, left, right, 0);
    n = numaGetCount(nad);
    fa = numaGetFArray(nad, L_NOCOPY);
    if (type == L_CONTINUED_BORDER) {
        for (i = 0; i < left; i++)
            fa[i] = fa[left];
        for (i = n - right; i < n; i++)
            fa[i] = fa[n - right - 1];
    } else {  /* type == L_MIRRORED_BORDER */
        for (i = 0; i < left; i++)
            fa[i] = fa[2 * left - 1 - i];
        for (i = 0; i < right; i++)
            fa[n - right + i] = fa[n - right - i - 1];
    }
 
    return nad;
}
 
 
NUMA *
numaRemoveBorder(NUMA      *nas,
                 l_int32    left,
                 l_int32    right)
{
l_int32     i, n, len;
l_float32   startx, delx;
l_float32  *fas, *fad;
NUMA       *nad;
 
    PROCNAME("numaRemoveBorder");
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", procName, NULL);
    if (left < 0) left = 0;
    if (right < 0) right = 0;
    if (left == 0 && right == 0)
        return numaCopy(nas);
 
    n = numaGetCount(nas);
    if ((len = n - left - right) < 0)
        return (NUMA *)ERROR_PTR("len < 0 after removal", procName, NULL);
    nad = numaMakeConstant(0, len);
    numaGetParameters(nas, &startx, &delx);
    numaSetParameters(nad, startx + delx * left, delx);
    fas = numaGetFArray(nas, L_NOCOPY);
    fad = numaGetFArray(nad, L_NOCOPY);
    for (i = 0; i < len; i++)
        fad[i] = fas[left + i];
 
    return nad;
}
 
 
l_ok
numaCountNonzeroRuns(NUMA     *na,
                     l_int32  *pcount)
{
l_int32  n, i, val, count, inrun;
 
    PROCNAME("numaCountNonzeroRuns");
 
    if (!pcount)
        return ERROR_INT("&count not defined", procName, 1);
    *pcount = 0;
    if (!na)
        return ERROR_INT("na not defined", procName, 1);
    if ((n = numaGetCount(na)) == 0)
        return ERROR_INT("na is empty", procName, 1);
 
    count = 0;
    inrun = FALSE;
    for (i = 0; i < n; i++) {
        numaGetIValue(na, i, &val);
        if (!inrun && val > 0) {
            count++;
            inrun = TRUE;
        } else if (inrun && val == 0) {
            inrun = FALSE;
        }
    }
    *pcount = count;
    return 0;
}
 
 
l_ok
numaGetNonzeroRange(NUMA      *na,
                    l_float32  eps,
                    l_int32   *pfirst,
                    l_int32   *plast)
{
l_int32    n, i, found;
l_float32  val;
 
    PROCNAME("numaGetNonzeroRange");
 
    if (pfirst) *pfirst = 0;
    if (plast) *plast = 0;
    if (!pfirst || !plast)
        return ERROR_INT("pfirst and plast not both defined", procName, 1);
    if (!na)
        return ERROR_INT("na not defined", procName, 1);
    if ((n = numaGetCount(na)) == 0)
        return ERROR_INT("na is empty", procName, 1);
 
    found = FALSE;
    for (i = 0; i < n; i++) {
        numaGetFValue(na, i, &val);
        if (val > eps) {
            found = TRUE;
            break;
        }
    }
    if (!found) {
        *pfirst = n - 1;
        *plast = 0;
        return 1;
    }
 
    *pfirst = i;
    for (i = n - 1; i >= 0; i--) {
        numaGetFValue(na, i, &val);
        if (val > eps)
            break;
    }
    *plast = i;
    return 0;
}
 
 
l_ok
numaGetCountRelativeToZero(NUMA     *na,
                           l_int32   type,
                           l_int32  *pcount)
{
l_int32    n, i, count;
l_float32  val;
 
    PROCNAME("numaGetCountRelativeToZero");
 
    if (!pcount)
        return ERROR_INT("&count not defined", procName, 1);
    *pcount = 0;
    if (!na)
        return ERROR_INT("na not defined", procName, 1);
    if ((n = numaGetCount(na)) == 0)
        return ERROR_INT("na is empty", procName, 1);
 
    for (i = 0, count = 0; i < n; i++) {
        numaGetFValue(na, i, &val);
        if (type == L_LESS_THAN_ZERO && val < 0.0)
            count++;
        else if (type == L_EQUAL_TO_ZERO && val == 0.0)
            count++;
        else if (type == L_GREATER_THAN_ZERO && val > 0.0)
            count++;
    }
 
    *pcount = count;
    return 0;
}
 
 
NUMA *
numaClipToInterval(NUMA    *nas,
                   l_int32  first,
                   l_int32  last)
{
l_int32    n, i;
l_float32  val, startx, delx;
NUMA      *nad;
 
    PROCNAME("numaClipToInterval");
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", procName, NULL);
    if ((n = numaGetCount(nas)) == 0)
        return (NUMA *)ERROR_PTR("nas is empty", procName, NULL);
    if (first < 0 || first > last)
        return (NUMA *)ERROR_PTR("range not valid", procName, NULL);
    if (first >= n)
        return (NUMA *)ERROR_PTR("no elements in range", procName, NULL);
 
    last = L_MIN(last, n - 1);
    if ((nad = numaCreate(last - first + 1)) == NULL)
        return (NUMA *)ERROR_PTR("nad not made", procName, NULL);
    for (i = first; i <= last; i++) {
        numaGetFValue(nas, i, &val);
        numaAddNumber(nad, val);
    }
    numaGetParameters(nas, &startx, &delx);
    numaSetParameters(nad, startx + first * delx, delx);
    return nad;
}
 
 
NUMA *
numaMakeThresholdIndicator(NUMA      *nas,
                           l_float32  thresh,
                           l_int32    type)
{
l_int32    n, i, ival;
l_float32  fval;
NUMA      *nai;
 
    PROCNAME("numaMakeThresholdIndicator");
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", procName, NULL);
    if ((n = numaGetCount(nas)) == 0)
        return (NUMA *)ERROR_PTR("nas is empty", procName, NULL);
 
    nai = numaCreate(n);
    for (i = 0; i < n; i++) {
        numaGetFValue(nas, i, &fval);
        ival = 0;
        switch (type)
        {
        case L_SELECT_IF_LT:
            if (fval < thresh) ival = 1;
            break;
        case L_SELECT_IF_GT:
            if (fval > thresh) ival = 1;
            break;
        case L_SELECT_IF_LTE:
            if (fval <= thresh) ival = 1;
            break;
        case L_SELECT_IF_GTE:
            if (fval >= thresh) ival = 1;
            break;
        default:
            numaDestroy(&nai);
            return (NUMA *)ERROR_PTR("invalid type", procName, NULL);
        }
        numaAddNumber(nai, ival);
    }
 
    return nai;
}
 
 
NUMA *
numaUniformSampling(NUMA    *nas,
                    l_int32  nsamp)
{
l_int32     n, i, j, ileft, iright;
l_float32   left, right, binsize, lfract, rfract, sum, startx, delx;
l_float32  *array;
NUMA       *nad;
 
    PROCNAME("numaUniformSampling");
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", procName, NULL);
    if ((n = numaGetCount(nas)) == 0)
        return (NUMA *)ERROR_PTR("nas is empty", procName, NULL);
    if (nsamp <= 0)
        return (NUMA *)ERROR_PTR("nsamp must be > 0", procName, NULL);
 
    nad = numaCreate(nsamp);
    array = numaGetFArray(nas, L_NOCOPY);
    binsize = (l_float32)n / (l_float32)nsamp;
    numaGetParameters(nas, &startx, &delx);
    numaSetParameters(nad, startx, binsize * delx);
    left = 0.0;
    for (i = 0; i < nsamp; i++) {
        sum = 0.0;
        right = left + binsize;
        ileft = (l_int32)left;
        lfract = 1.0 - left + ileft;
        if (lfract >= 1.0)  /* on left bin boundary */
            lfract = 0.0;
        iright = (l_int32)right;
        rfract = right - iright;
        iright = L_MIN(iright, n - 1);
        if (ileft == iright) {  /* both are within the same original sample */
            sum += (lfract + rfract - 1.0) * array[ileft];
        } else {
            if (lfract > 0.0001)  /* left fraction */
                sum += lfract * array[ileft];
            if (rfract > 0.0001)  /* right fraction */
                sum += rfract * array[iright];
            for (j = ileft + 1; j < iright; j++)  /* entire pixels */
                sum += array[j];
        }
 
        numaAddNumber(nad, sum);
        left = right;
    }
    return nad;
}
 
 
NUMA *
numaReverse(NUMA  *nad,
            NUMA  *nas)
{
l_int32    n, i;
l_float32  val1, val2;
 
    PROCNAME("numaReverse");
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", procName, NULL);
    if (nad && nas != nad)
        return (NUMA *)ERROR_PTR("nad defined but != nas", procName, NULL);
 
    n = numaGetCount(nas);
    if (nad) {  /* in-place */
        for (i = 0; i < n / 2; i++) {
            numaGetFValue(nad, i, &val1);
            numaGetFValue(nad, n - i - 1, &val2);
            numaSetValue(nad, i, val2);
            numaSetValue(nad, n - i - 1, val1);
        }
    } else {
        nad = numaCreate(n);
        for (i = n - 1; i >= 0; i--) {
            numaGetFValue(nas, i, &val1);
            numaAddNumber(nad, val1);
        }
    }
 
        /* Reverse the startx and delx fields */
    nad->startx = nas->startx + (n - 1) * nas->delx;
    nad->delx = -nas->delx;
    return nad;
}
 
 
/*----------------------------------------------------------------------*
 *                       Signal feature extraction                      *
 *----------------------------------------------------------------------*/
NUMA *
numaLowPassIntervals(NUMA      *nas,
                     l_float32  thresh,
                     l_float32  maxn)
{
l_int32    n, i, inrun;
l_float32  maxval, threshval, fval, startx, delx, x0, x1;
NUMA      *nad;
 
    PROCNAME("numaLowPassIntervals");
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", procName, NULL);
    if ((n = numaGetCount(nas)) == 0)
        return (NUMA *)ERROR_PTR("nas is empty", procName, NULL);
    if (thresh < 0.0 || thresh > 1.0)
        return (NUMA *)ERROR_PTR("invalid thresh", procName, NULL);
 
        /* The input threshold is a fraction of the max.
         * The first entry in nad is the value of the max. */
    if (maxn == 0.0)
        numaGetMax(nas, &maxval, NULL);
    else
        maxval = maxn;
    numaGetParameters(nas, &startx, &delx);
    threshval = thresh * maxval;
    nad = numaCreate(0);
    numaAddNumber(nad, maxval);
 
        /* Write pairs of pts (x0, x1) for the intervals */
    inrun = FALSE;
    for (i = 0; i < n; i++) {
        numaGetFValue(nas, i, &fval);
        if (fval < threshval && inrun == FALSE) {  /* start a new run */
            inrun = TRUE;
            x0 = startx + i * delx;
        } else if (fval > threshval && inrun == TRUE) {  /* end the run */
            inrun = FALSE;
            x1 = startx + i * delx;
            numaAddNumber(nad, x0);
            numaAddNumber(nad, x1);
        }
    }
    if (inrun == TRUE) {  /* must end the last run */
        x1 = startx + (n - 1) * delx;
        numaAddNumber(nad, x0);
        numaAddNumber(nad, x1);
    }
 
    return nad;
}
 
 
NUMA *
numaThresholdEdges(NUMA      *nas,
                   l_float32  thresh1,
                   l_float32  thresh2,
                   l_float32  maxn)
{
l_int32    n, i, istart, inband, output, sign;
l_int32    startbelow, below, above, belowlast, abovelast;
l_float32  maxval, threshval1, threshval2, fval, startx, delx, x0, x1;
NUMA      *nad;
 
    PROCNAME("numaThresholdEdges");
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", procName, NULL);
    if ((n = numaGetCount(nas)) == 0)
        return (NUMA *)ERROR_PTR("nas is empty", procName, NULL);
    if (thresh1 < 0.0 || thresh1 > 1.0 || thresh2 < 0.0 || thresh2 > 1.0)
        return (NUMA *)ERROR_PTR("invalid thresholds", procName, NULL);
    if (thresh2 < thresh1)
        return (NUMA *)ERROR_PTR("thresh2 < thresh1", procName, NULL);
 
        /* The input thresholds are fractions of the max.
         * The first entry in nad is the value of the max used
         * here for normalization. */
    if (maxn == 0.0)
        numaGetMax(nas, &maxval, NULL);
    else
        maxval = maxn;
    numaGetMax(nas, &maxval, NULL);
    numaGetParameters(nas, &startx, &delx);
    threshval1 = thresh1 * maxval;
    threshval2 = thresh2 * maxval;
    nad = numaCreate(0);
    numaAddNumber(nad, maxval);
 
        /* Write triplets of pts (x0, x1, sign) for the edges.
         * First make sure we start search from outside the band.
         * Only one of {belowlast, abovelast} is true. */
    for (i = 0; i < n; i++) {
        istart = i;
        numaGetFValue(nas, i, &fval);
        belowlast = (fval < threshval1) ? TRUE : FALSE;
        abovelast = (fval > threshval2) ? TRUE : FALSE;
        if (belowlast == TRUE || abovelast == TRUE)
            break;
    }
    if (istart == n)  /* no intervals found */
        return nad;
 
        /* x0 and x1 can only be set from outside the edge.
         * They are the values just before entering the band,
         * and just after entering the band.  We can jump through
         * the band, in which case they differ by one index in nas. */
    inband = FALSE;
    startbelow = belowlast; /* one of these is true */
    output = FALSE;
    x0 = startx + istart * delx;
    for (i = istart + 1; i < n; i++) {
        numaGetFValue(nas, i, &fval);
        below = (fval < threshval1) ? TRUE : FALSE;
        above = (fval > threshval2) ? TRUE : FALSE;
        if (!inband && belowlast && above) {  /* full jump up */
            x1 = startx + i * delx;
            sign = 1;
            startbelow = FALSE;  /* for the next transition */
            output = TRUE;
        } else if (!inband && abovelast && below) {  /* full jump down */
            x1 = startx + i * delx;
            sign = -1;
            startbelow = TRUE;  /* for the next transition */
            output = TRUE;
        } else if (inband && startbelow && above) {  /* exit rising; success */
            x1 = startx + i * delx;
            sign = 1;
            inband = FALSE;
            startbelow = FALSE;  /* for the next transition */
            output = TRUE;
        } else if (inband && !startbelow && below) {
                /* exit falling; success */
            x1 = startx + i * delx;
            sign = -1;
            inband = FALSE;
            startbelow = TRUE;  /* for the next transition */
            output = TRUE;
        } else if (inband && !startbelow && above) {  /* exit rising; failure */
            x0 = startx + i * delx;
            inband = FALSE;
        } else if (inband && startbelow && below) {  /* exit falling; failure */
            x0 = startx + i * delx;
            inband = FALSE;
        } else if (!inband && !above && !below) {  /* enter */
            inband = TRUE;
            startbelow = belowlast;
        } else if (!inband && (above || below)) {  /* outside and remaining */
            x0 = startx + i * delx;  /* update position */
        }
        belowlast = below;
        abovelast = above;
        if (output) {  /* we have exited; save new x0 */
            numaAddNumber(nad, x0);
            numaAddNumber(nad, x1);
            numaAddNumber(nad, sign);
            output = FALSE;
            x0 = startx + i * delx;
        }
    }
 
    return nad;
}
 
 
l_int32
numaGetSpanValues(NUMA    *na,
                  l_int32  span,
                  l_int32 *pstart,
                  l_int32 *pend)
{
l_int32  n, nspans;
 
    PROCNAME("numaGetSpanValues");
 
    if (!na)
        return ERROR_INT("na not defined", procName, 1);
    if ((n = numaGetCount(na)) == 0)
        return ERROR_INT("na is empty", procName, 1);
    if (n % 2 != 1)
        return ERROR_INT("n is not odd", procName, 1);
    nspans = n / 2;
    if (nspans < 0 || span >= nspans)
        return ERROR_INT("invalid span", procName, 1);
 
    if (pstart) numaGetIValue(na, 2 * span + 1, pstart);
    if (pend) numaGetIValue(na, 2 * span + 2, pend);
    return 0;
}
 
 
l_int32
numaGetEdgeValues(NUMA    *na,
                  l_int32  edge,
                  l_int32 *pstart,
                  l_int32 *pend,
                  l_int32 *psign)
{
l_int32  n, nedges;
 
    PROCNAME("numaGetEdgeValues");
 
    if (!na)
        return ERROR_INT("na not defined", procName, 1);
    if ((n = numaGetCount(na)) == 0)
        return ERROR_INT("na is empty", procName, 1);
    if (n % 3 != 1)
        return ERROR_INT("n % 3 is not 1", procName, 1);
    nedges = (n - 1) / 3;
    if (edge < 0 || edge >= nedges)
        return ERROR_INT("invalid edge", procName, 1);
 
    if (pstart) numaGetIValue(na, 3 * edge + 1, pstart);
    if (pend) numaGetIValue(na, 3 * edge + 2, pend);
    if (psign) numaGetIValue(na, 3 * edge + 3, psign);
    return 0;
}
 
 
/*----------------------------------------------------------------------*
 *                             Interpolation                            *
 *----------------------------------------------------------------------*/
l_ok
numaInterpolateEqxVal(l_float32   startx,
                      l_float32   deltax,
                      NUMA       *nay,
                      l_int32     type,
                      l_float32   xval,
                      l_float32  *pyval)
{
l_int32     i, n, i1, i2, i3;
l_float32   x1, x2, x3, fy1, fy2, fy3, d1, d2, d3, del, fi, maxx;
l_float32  *fa;
 
    PROCNAME("numaInterpolateEqxVal");
 
    if (!pyval)
        return ERROR_INT("&yval not defined", procName, 1);
    *pyval = 0.0;
    if (!nay)
        return ERROR_INT("nay not defined", procName, 1);
    if (deltax <= 0.0)
        return ERROR_INT("deltax not > 0", procName, 1);
    if (type != L_LINEAR_INTERP && type != L_QUADRATIC_INTERP)
        return ERROR_INT("invalid interp type", procName, 1);
    if ((n = numaGetCount(nay)) < 2)
        return ERROR_INT("not enough points", procName, 1);
    if (type == L_QUADRATIC_INTERP && n == 2) {
        type = L_LINEAR_INTERP;
        L_WARNING("only 2 points; using linear interp\n", procName);
    }
    maxx = startx + deltax * (n - 1);
    if (xval < startx || xval > maxx)
        return ERROR_INT("xval is out of bounds", procName, 1);
 
    fa = numaGetFArray(nay, L_NOCOPY);
    fi = (xval - startx) / deltax;
    i = (l_int32)fi;
    del = fi - i;
    if (del == 0.0) {  /* no interpolation required */
        *pyval = fa[i];
        return 0;
    }
 
    if (type == L_LINEAR_INTERP) {
        *pyval = fa[i] + del * (fa[i + 1] - fa[i]);
        return 0;
    }
 
        /* Quadratic interpolation */
    d1 = d3 = 0.5 / (deltax * deltax);
    d2 = -2. * d1;
    if (i == 0) {
        i1 = i;
        i2 = i + 1;
        i3 = i + 2;
    } else {
        i1 = i - 1;
        i2 = i;
        i3 = i + 1;
    }
    x1 = startx + i1 * deltax;
    x2 = startx + i2 * deltax;
    x3 = startx + i3 * deltax;
    fy1 = d1 * fa[i1];
    fy2 = d2 * fa[i2];
    fy3 = d3 * fa[i3];
    *pyval = fy1 * (xval - x2) * (xval - x3) +
             fy2 * (xval - x1) * (xval - x3) +
             fy3 * (xval - x1) * (xval - x2);
    return 0;
}
 
 
l_ok
numaInterpolateArbxVal(NUMA       *nax,
                       NUMA       *nay,
                       l_int32     type,
                       l_float32   xval,
                       l_float32  *pyval)
{
l_int32     i, im, nx, ny, i1, i2, i3;
l_float32   delu, dell, fract, d1, d2, d3;
l_float32   minx, maxx;
l_float32  *fax, *fay;
 
    PROCNAME("numaInterpolateArbxVal");
 
    if (!pyval)
        return ERROR_INT("&yval not defined", procName, 1);
    *pyval = 0.0;
    if (!nax)
        return ERROR_INT("nax not defined", procName, 1);
    if (!nay)
        return ERROR_INT("nay not defined", procName, 1);
    if (type != L_LINEAR_INTERP && type != L_QUADRATIC_INTERP)
        return ERROR_INT("invalid interp type", procName, 1);
    ny = numaGetCount(nay);
    nx = numaGetCount(nax);
    if (nx != ny)
        return ERROR_INT("nax and nay not same size arrays", procName, 1);
    if (ny < 2)
        return ERROR_INT("not enough points", procName, 1);
    if (type == L_QUADRATIC_INTERP && ny == 2) {
        type = L_LINEAR_INTERP;
        L_WARNING("only 2 points; using linear interp\n", procName);
    }
    numaGetFValue(nax, 0, &minx);
    numaGetFValue(nax, nx - 1, &maxx);
    if (xval < minx || xval > maxx)
        return ERROR_INT("xval is out of bounds", procName, 1);
 
    fax = numaGetFArray(nax, L_NOCOPY);
    fay = numaGetFArray(nay, L_NOCOPY);
 
        /* Linear search for interval.  We are guaranteed
         * to either return or break out of the loop.
         * In addition, we are assured that fax[i] - fax[im] > 0.0 */
    if (xval == fax[0]) {
        *pyval = fay[0];
        return 0;
    }
    im = 0;
    dell = 0.0;
    for (i = 1; i < nx; i++) {
        delu = fax[i] - xval;
        if (delu >= 0.0) {  /* we've passed it */
            if (delu == 0.0) {
                *pyval = fay[i];
                return 0;
            }
            im = i - 1;
            dell = xval - fax[im];  /* >= 0 */
            break;
        }
    }
    fract = dell / (fax[i] - fax[im]);
 
    if (type == L_LINEAR_INTERP) {
        *pyval = fay[i] + fract * (fay[i + 1] - fay[i]);
        return 0;
    }
 
        /* Quadratic interpolation */
    if (im == 0) {
        i1 = im;
        i2 = im + 1;
        i3 = im + 2;
    } else {
        i1 = im - 1;
        i2 = im;
        i3 = im + 1;
    }
    d1 = (fax[i1] - fax[i2]) * (fax[i1] - fax[i3]);
    d2 = (fax[i2] - fax[i1]) * (fax[i2] - fax[i3]);
    d3 = (fax[i3] - fax[i1]) * (fax[i3] - fax[i2]);
    *pyval = fay[i1] * (xval - fax[i2]) * (xval - fax[i3]) / d1 +
             fay[i2] * (xval - fax[i1]) * (xval - fax[i3]) / d2 +
             fay[i3] * (xval - fax[i1]) * (xval - fax[i2]) / d3;
    return 0;
}
 
 
l_ok
numaInterpolateEqxInterval(l_float32  startx,
                           l_float32  deltax,
                           NUMA      *nasy,
                           l_int32    type,
                           l_float32  x0,
                           l_float32  x1,
                           l_int32    npts,
                           NUMA     **pnax,
                           NUMA     **pnay)
{
l_int32     i, n;
l_float32   x, yval, maxx, delx;
NUMA       *nax, *nay;
 
    PROCNAME("numaInterpolateEqxInterval");
 
    if (pnax) *pnax = NULL;
    if (!pnay)
        return ERROR_INT("&nay not defined", procName, 1);
    *pnay = NULL;
    if (!nasy)
        return ERROR_INT("nasy not defined", procName, 1);
    if ((n = numaGetCount(nasy)) < 2)
        return ERROR_INT("n < 2", procName, 1);
    if (deltax <= 0.0)
        return ERROR_INT("deltax not > 0", procName, 1);
    if (type != L_LINEAR_INTERP && type != L_QUADRATIC_INTERP)
        return ERROR_INT("invalid interp type", procName, 1);
    if (type == L_QUADRATIC_INTERP && n == 2) {
        type = L_LINEAR_INTERP;
        L_WARNING("only 2 points; using linear interp\n", procName);
    }
    maxx = startx + deltax * (n - 1);
    if (x0 < startx || x1 > maxx || x1 <= x0)
        return ERROR_INT("[x0 ... x1] is not valid", procName, 1);
    if (npts < 3)
        return ERROR_INT("npts < 3", procName, 1);
    delx = (x1 - x0) / (l_float32)(npts - 1);  /* delx is for output nay */
 
    if ((nay = numaCreate(npts)) == NULL)
        return ERROR_INT("nay not made", procName, 1);
    numaSetParameters(nay, x0, delx);
    *pnay = nay;
    if (pnax) {
        nax = numaCreate(npts);
        *pnax = nax;
    }
 
    for (i = 0; i < npts; i++) {
        x = x0 + i * delx;
        if (pnax)
            numaAddNumber(nax, x);
        numaInterpolateEqxVal(startx, deltax, nasy, type, x, &yval);
        numaAddNumber(nay, yval);
    }
 
    return 0;
}
 
 
l_ok
numaInterpolateArbxInterval(NUMA       *nax,
                            NUMA       *nay,
                            l_int32     type,
                            l_float32   x0,
                            l_float32   x1,
                            l_int32     npts,
                            NUMA      **pnadx,
                            NUMA      **pnady)
{
l_int32     i, im, j, nx, ny, i1, i2, i3, sorted;
l_int32    *index;
l_float32   del, xval, yval, excess, fract, minx, maxx, d1, d2, d3;
l_float32  *fax, *fay;
NUMA       *nasx, *nasy, *nadx, *nady;
 
    PROCNAME("numaInterpolateArbxInterval");
 
    if (pnadx) *pnadx = NULL;
    if (!pnady)
        return ERROR_INT("&nady not defined", procName, 1);
    *pnady = NULL;
    if (!nay)
        return ERROR_INT("nay not defined", procName, 1);
    if (!nax)
        return ERROR_INT("nax not defined", procName, 1);
    if (type != L_LINEAR_INTERP && type != L_QUADRATIC_INTERP)
        return ERROR_INT("invalid interp type", procName, 1);
    if (x0 > x1)
        return ERROR_INT("x0 > x1", procName, 1);
    ny = numaGetCount(nay);
    nx = numaGetCount(nax);
    if (nx != ny)
        return ERROR_INT("nax and nay not same size arrays", procName, 1);
    if (ny < 2)
        return ERROR_INT("not enough points", procName, 1);
    if (type == L_QUADRATIC_INTERP && ny == 2) {
        type = L_LINEAR_INTERP;
        L_WARNING("only 2 points; using linear interp\n", procName);
    }
    numaGetMin(nax, &minx, NULL);
    numaGetMax(nax, &maxx, NULL);
    if (x0 < minx || x1 > maxx)
        return ERROR_INT("xval is out of bounds", procName, 1);
 
        /* Make sure that nax is sorted in increasing order */
    numaIsSorted(nax, L_SORT_INCREASING, &sorted);
    if (!sorted) {
        L_WARNING("we are sorting nax in increasing order\n", procName);
        numaSortPair(nax, nay, L_SORT_INCREASING, &nasx, &nasy);
    } else {
        nasx = numaClone(nax);
        nasy = numaClone(nay);
    }
 
    fax = numaGetFArray(nasx, L_NOCOPY);
    fay = numaGetFArray(nasy, L_NOCOPY);
 
        /* Get array of indices into fax for interpolated locations */
    if ((index = (l_int32 *)LEPT_CALLOC(npts, sizeof(l_int32))) == NULL) {
        numaDestroy(&nasx);
        numaDestroy(&nasy);
        return ERROR_INT("ind not made", procName, 1);
    }
    del = (x1 - x0) / (npts - 1.0);
    for (i = 0, j = 0; j < nx && i < npts; i++) {
        xval = x0 + i * del;
        while (j < nx - 1 && xval > fax[j])
            j++;
        if (xval == fax[j])
            index[i] = L_MIN(j, nx - 1);
        else    /* the index of fax[] is just below xval */
            index[i] = L_MAX(j - 1, 0);
    }
 
        /* For each point to be interpolated, get the y value */
    nady = numaCreate(npts);
    *pnady = nady;
    if (pnadx) {
        nadx = numaCreate(npts);
        *pnadx = nadx;
    }
    for (i = 0; i < npts; i++) {
        xval = x0 + i * del;
        if (pnadx)
            numaAddNumber(nadx, xval);
        im = index[i];
        excess = xval - fax[im];
        if (excess == 0.0) {
            numaAddNumber(nady, fay[im]);
            continue;
        }
        fract = excess / (fax[im + 1] - fax[im]);
 
        if (type == L_LINEAR_INTERP) {
            yval = fay[im] + fract * (fay[im + 1] - fay[im]);
            numaAddNumber(nady, yval);
            continue;
        }
 
            /* Quadratic interpolation */
        if (im == 0) {
            i1 = im;
            i2 = im + 1;
            i3 = im + 2;
        } else {
            i1 = im - 1;
            i2 = im;
            i3 = im + 1;
        }
        d1 = (fax[i1] - fax[i2]) * (fax[i1] - fax[i3]);
        d2 = (fax[i2] - fax[i1]) * (fax[i2] - fax[i3]);
        d3 = (fax[i3] - fax[i1]) * (fax[i3] - fax[i2]);
        yval = fay[i1] * (xval - fax[i2]) * (xval - fax[i3]) / d1 +
               fay[i2] * (xval - fax[i1]) * (xval - fax[i3]) / d2 +
               fay[i3] * (xval - fax[i1]) * (xval - fax[i2]) / d3;
        numaAddNumber(nady, yval);
    }
 
    LEPT_FREE(index);
    numaDestroy(&nasx);
    numaDestroy(&nasy);
    return 0;
}
 
 
/*----------------------------------------------------------------------*
 *                     Functions requiring interpolation                *
 *----------------------------------------------------------------------*/
l_ok
numaFitMax(NUMA       *na,
           l_float32  *pmaxval,
           NUMA       *naloc,
           l_float32  *pmaxloc)
{
l_float32  val;
l_float32  smaxval;  /* start value of maximum sample, before interpolating */
l_int32    n, imaxloc;
l_float32  x1, x2, x3, y1, y2, y3, c1, c2, c3, a, b, xmax, ymax;
 
    PROCNAME("numaFitMax");
 
    if (pmaxval) *pmaxval = 0.0;
    if (pmaxloc) *pmaxloc = 0.0;
    if (!na)
        return ERROR_INT("na not defined", procName, 1);
    if ((n = numaGetCount(na)) == 0)
        return ERROR_INT("na is empty", procName, 1);
    if (!pmaxval)
        return ERROR_INT("&maxval not defined", procName, 1);
    if (!pmaxloc)
        return ERROR_INT("&maxloc not defined", procName, 1);
    if (naloc) {
        if (n != numaGetCount(naloc))
            return ERROR_INT("na and naloc of unequal size", procName, 1);
    }
 
    numaGetMax(na, &smaxval, &imaxloc);
 
        /* Simple case: max is at end point */
    if (imaxloc == 0 || imaxloc == n - 1) {
        *pmaxval = smaxval;
        if (naloc) {
            numaGetFValue(naloc, imaxloc, &val);
            *pmaxloc = val;
        } else {
            *pmaxloc = imaxloc;
        }
        return 0;
    }
 
        /* Interior point; use quadratic interpolation */
    y2 = smaxval;
    numaGetFValue(na, imaxloc - 1, &val);
    y1 = val;
    numaGetFValue(na, imaxloc + 1, &val);
    y3 = val;
    if (naloc) {
        numaGetFValue(naloc, imaxloc - 1, &val);
        x1 = val;
        numaGetFValue(naloc, imaxloc, &val);
        x2 = val;
        numaGetFValue(naloc, imaxloc + 1, &val);
        x3 = val;
    } else {
        x1 = imaxloc - 1;
        x2 = imaxloc;
        x3 = imaxloc + 1;
    }
 
        /* Can't interpolate; just use the max val in na
         * and the corresponding one in naloc */
    if (x1 == x2 || x1 == x3 || x2 == x3) {
        *pmaxval = y2;
        *pmaxloc = x2;
        return 0;
    }
 
        /* Use lagrangian interpolation; set dy/dx = 0 */
    c1 = y1 / ((x1 - x2) * (x1 - x3));
    c2 = y2 / ((x2 - x1) * (x2 - x3));
    c3 = y3 / ((x3 - x1) * (x3 - x2));
    a = c1 + c2 + c3;
    b = c1 * (x2 + x3) + c2 * (x1 + x3) + c3 * (x1 + x2);
    xmax = b / (2 * a);
    ymax = c1 * (xmax - x2) * (xmax - x3) +
           c2 * (xmax - x1) * (xmax - x3) +
           c3 * (xmax - x1) * (xmax - x2);
    *pmaxval = ymax;
    *pmaxloc = xmax;
 
    return 0;
}
 
 
l_ok
numaDifferentiateInterval(NUMA       *nax,
                          NUMA       *nay,
                          l_float32   x0,
                          l_float32   x1,
                          l_int32     npts,
                          NUMA      **pnadx,
                          NUMA      **pnady)
{
l_int32     i, nx, ny;
l_float32   minx, maxx, der, invdel;
l_float32  *fay;
NUMA       *nady, *naiy;
 
    PROCNAME("numaDifferentiateInterval");
 
    if (pnadx) *pnadx = NULL;
    if (!pnady)
        return ERROR_INT("&nady not defined", procName, 1);
    *pnady = NULL;
    if (!nay)
        return ERROR_INT("nay not defined", procName, 1);
    if (!nax)
        return ERROR_INT("nax not defined", procName, 1);
    if (x0 > x1)
        return ERROR_INT("x0 > x1", procName, 1);
    ny = numaGetCount(nay);
    nx = numaGetCount(nax);
    if (nx != ny)
        return ERROR_INT("nax and nay not same size arrays", procName, 1);
    if (ny < 2)
        return ERROR_INT("not enough points", procName, 1);
    numaGetMin(nax, &minx, NULL);
    numaGetMax(nax, &maxx, NULL);
    if (x0 < minx || x1 > maxx)
        return ERROR_INT("xval is out of bounds", procName, 1);
    if (npts < 2)
        return ERROR_INT("npts < 2", procName, 1);
 
        /* Generate interpolated array over specified interval */
    if (numaInterpolateArbxInterval(nax, nay, L_LINEAR_INTERP, x0, x1,
                                    npts, pnadx, &naiy))
        return ERROR_INT("interpolation failed", procName, 1);
 
    nady = numaCreate(npts);
    *pnady = nady;
    invdel = 0.5 * ((l_float32)npts - 1.0) / (x1 - x0);
    fay = numaGetFArray(naiy, L_NOCOPY);
 
        /* Compute and save derivatives */
    der = 0.5 * invdel * (fay[1] - fay[0]);
    numaAddNumber(nady, der);
    for (i = 1; i < npts - 1; i++)  {
        der = invdel * (fay[i + 1] - fay[i - 1]);
        numaAddNumber(nady, der);
    }
    der = 0.5 * invdel * (fay[npts - 1] - fay[npts - 2]);
    numaAddNumber(nady, der);
 
    numaDestroy(&naiy);
    return 0;
}
 
 
l_ok
numaIntegrateInterval(NUMA       *nax,
                      NUMA       *nay,
                      l_float32   x0,
                      l_float32   x1,
                      l_int32     npts,
                      l_float32  *psum)
{
l_int32     i, nx, ny;
l_float32   minx, maxx, sum, del;
l_float32  *fay;
NUMA       *naiy;
 
    PROCNAME("numaIntegrateInterval");
 
    if (!psum)
        return ERROR_INT("&sum not defined", procName, 1);
    *psum = 0.0;
    if (!nay)
        return ERROR_INT("nay not defined", procName, 1);
    if (!nax)
        return ERROR_INT("nax not defined", procName, 1);
    if (x0 > x1)
        return ERROR_INT("x0 > x1", procName, 1);
    if (npts < 2)
        return ERROR_INT("npts < 2", procName, 1);
    ny = numaGetCount(nay);
    nx = numaGetCount(nax);
    if (nx != ny)
        return ERROR_INT("nax and nay not same size arrays", procName, 1);
    if (ny < 2)
        return ERROR_INT("not enough points", procName, 1);
    numaGetMin(nax, &minx, NULL);
    numaGetMax(nax, &maxx, NULL);
    if (x0 < minx || x1 > maxx)
        return ERROR_INT("xval is out of bounds", procName, 1);
 
        /* Generate interpolated array over specified interval */
    if (numaInterpolateArbxInterval(nax, nay, L_LINEAR_INTERP, x0, x1,
                                    npts, NULL, &naiy))
        return ERROR_INT("interpolation failed", procName, 1);
 
    del = (x1 - x0) / ((l_float32)npts - 1.0);
    fay = numaGetFArray(naiy, L_NOCOPY);
 
        /* Compute integral (simple trapezoid) */
    sum = 0.5 * (fay[0] + fay[npts - 1]);
    for (i = 1; i < npts - 1; i++)
        sum += fay[i];
    *psum = del * sum;
 
    numaDestroy(&naiy);
    return 0;
}
 
 
/*----------------------------------------------------------------------*
 *                                Sorting                               *
 *----------------------------------------------------------------------*/
l_ok
numaSortGeneral(NUMA    *na,
                NUMA   **pnasort,
                NUMA   **pnaindex,
                NUMA   **pnainvert,
                l_int32  sortorder,
                l_int32  sorttype)
{
l_int32    isize;
l_float32  size;
NUMA      *naindex = NULL;
 
    PROCNAME("numaSortGeneral");
 
    if (pnasort) *pnasort = NULL;
    if (pnaindex) *pnaindex = NULL;
    if (pnainvert) *pnainvert = NULL;
    if (!na)
        return ERROR_INT("na not defined", procName, 1);
    if (sortorder != L_SORT_INCREASING && sortorder != L_SORT_DECREASING)
        return ERROR_INT("invalid sort order", procName, 1);
    if (sorttype != L_SHELL_SORT && sorttype != L_BIN_SORT)
        return ERROR_INT("invalid sort type", procName, 1);
    if (!pnasort && !pnaindex && !pnainvert)
        return ERROR_INT("nothing to do", procName, 1);
 
    if (sorttype == L_BIN_SORT) {
        numaGetMax(na, &size, NULL);
        isize = (l_int32)size;
        if (isize > MaxInitPtraSize - 1) {
            L_WARNING("array too large; using shell sort\n", procName);
            sorttype = L_SHELL_SORT;
        } else {
            naindex = numaGetBinSortIndex(na, sortorder);
        }
    }
 
    if (sorttype == L_SHELL_SORT)
        naindex = numaGetSortIndex(na, sortorder);
 
    if (pnasort)
        *pnasort = numaSortByIndex(na, naindex);
    if (pnainvert)
        *pnainvert = numaInvertMap(naindex);
    if (pnaindex)
        *pnaindex = naindex;
    else
        numaDestroy(&naindex);
    return 0;
}
 
 
NUMA *
numaSortAutoSelect(NUMA    *nas,
                   l_int32  sortorder)
{
l_int32  type;
 
    PROCNAME("numaSortAutoSelect");
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", procName, NULL);
    if (numaGetCount(nas) == 0) {
        L_WARNING("nas is empty; returning copy\n", procName);
        return numaCopy(nas);
    }
    if (sortorder != L_SORT_INCREASING && sortorder != L_SORT_DECREASING)
        return (NUMA *)ERROR_PTR("invalid sort order", procName, NULL);
 
    type = numaChooseSortType(nas);
    if (type != L_SHELL_SORT && type != L_BIN_SORT)
        return (NUMA *)ERROR_PTR("invalid sort type", procName, NULL);
 
    if (type == L_BIN_SORT)
        return numaBinSort(nas, sortorder);
    else  /* shell sort */
        return numaSort(NULL, nas, sortorder);
}
 
 
NUMA *
numaSortIndexAutoSelect(NUMA    *nas,
                        l_int32  sortorder)
{
l_int32  type;
 
    PROCNAME("numaSortIndexAutoSelect");
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", procName, NULL);
    if (numaGetCount(nas) == 0) {
        L_WARNING("nas is empty; returning copy\n", procName);
        return numaCopy(nas);
    }
    if (sortorder != L_SORT_INCREASING && sortorder != L_SORT_DECREASING)
        return (NUMA *)ERROR_PTR("invalid sort order", procName, NULL);
    type = numaChooseSortType(nas);
    if (type != L_SHELL_SORT && type != L_BIN_SORT)
        return (NUMA *)ERROR_PTR("invalid sort type", procName, NULL);
 
    if (type == L_BIN_SORT)
        return numaGetBinSortIndex(nas, sortorder);
    else  /* shell sort */
        return numaGetSortIndex(nas, sortorder);
}
 
 
l_int32
numaChooseSortType(NUMA  *nas)
{
l_int32    n;
l_float32  minval, maxval;
 
    PROCNAME("numaChooseSortType");
 
    if (!nas)
        return ERROR_INT("nas not defined", procName, UNDEF);
 
        /* If small histogram or negative values; use shell sort */
    numaGetMin(nas, &minval, NULL);
    n = numaGetCount(nas);
    if (minval < 0.0 || n < 200)
        return L_SHELL_SORT;
 
        /* If large maxval, use shell sort */
    numaGetMax(nas, &maxval, NULL);
    if (maxval > MaxInitPtraSize - 1)
        return L_SHELL_SORT;
 
        /* Otherwise, need to compare nlog(n) with maxval.
         * The factor of 0.003 was determined by comparing times for
         * different histogram sizes and maxval.  It is very small
         * because binsort is fast and shell sort gets slow for large n. */
    if (n * log((l_float32)n) < 0.003 * maxval)
        return L_SHELL_SORT;
    else
        return L_BIN_SORT;
}
 
 
NUMA *
numaSort(NUMA    *naout,
         NUMA    *nain,
         l_int32  sortorder)
{
l_int32     i, n, gap, j;
l_float32   tmp;
l_float32  *array;
 
    PROCNAME("numaSort");
 
    if (!nain)
        return (NUMA *)ERROR_PTR("nain not defined", procName, NULL);
    if (sortorder != L_SORT_INCREASING && sortorder != L_SORT_DECREASING)
        return (NUMA *)ERROR_PTR("invalid sort order", procName, NULL);
 
        /* Make naout if necessary; otherwise do in-place */
    if (!naout)
        naout = numaCopy(nain);
    else if (nain != naout)
        return (NUMA *)ERROR_PTR("invalid: not in-place", procName, NULL);
    if ((n = numaGetCount(naout)) == 0) {
        L_WARNING("naout is empty\n", procName);
        return naout;
    }
    array = naout->array;  /* operate directly on the array */
    n = numaGetCount(naout);
 
        /* Shell sort */
    for (gap = n/2; gap > 0; gap = gap / 2) {
        for (i = gap; i < n; i++) {
            for (j = i - gap; j >= 0; j -= gap) {
                if ((sortorder == L_SORT_INCREASING &&
                     array[j] > array[j + gap]) ||
                    (sortorder == L_SORT_DECREASING &&
                     array[j] < array[j + gap]))
                {
                    tmp = array[j];
                    array[j] = array[j + gap];
                    array[j + gap] = tmp;
                }
            }
        }
    }
 
    return naout;
}
 
 
NUMA *
numaBinSort(NUMA    *nas,
            l_int32  sortorder)
{
NUMA  *nat, *nad;
 
    PROCNAME("numaBinSort");
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", procName, NULL);
    if (numaGetCount(nas) == 0) {
        L_WARNING("nas is empty; returning copy\n", procName);
        return numaCopy(nas);
    }
    if (sortorder != L_SORT_INCREASING && sortorder != L_SORT_DECREASING)
        return (NUMA *)ERROR_PTR("invalid sort order", procName, NULL);
 
    if ((nat = numaGetBinSortIndex(nas, sortorder)) == NULL)
        return (NUMA *)ERROR_PTR("bin sort failed", procName, NULL);
    nad = numaSortByIndex(nas, nat);
    numaDestroy(&nat);
    return nad;
}
 
 
NUMA *
numaGetSortIndex(NUMA    *na,
                 l_int32  sortorder)
{
l_int32     i, n, gap, j;
l_float32   tmp;
l_float32  *array;   /* copy of input array */
l_float32  *iarray;  /* array of indices */
NUMA       *naisort;
 
    PROCNAME("numaGetSortIndex");
 
    if (!na)
        return (NUMA *)ERROR_PTR("na not defined", procName, NULL);
    if (numaGetCount(na) == 0) {
        L_WARNING("na is empty\n", procName);
        return numaCreate(1);
    }
    if (sortorder != L_SORT_INCREASING && sortorder != L_SORT_DECREASING)
        return (NUMA *)ERROR_PTR("invalid sortorder", procName, NULL);
 
    n = numaGetCount(na);
    if ((array = numaGetFArray(na, L_COPY)) == NULL)
        return (NUMA *)ERROR_PTR("array not made", procName, NULL);
    if ((iarray = (l_float32 *)LEPT_CALLOC(n, sizeof(l_float32))) == NULL) {
        LEPT_FREE(array);
        return (NUMA *)ERROR_PTR("iarray not made", procName, NULL);
    }
    for (i = 0; i < n; i++)
        iarray[i] = i;
 
        /* Shell sort */
    for (gap = n/2; gap > 0; gap = gap / 2) {
        for (i = gap; i < n; i++) {
            for (j = i - gap; j >= 0; j -= gap) {
                if ((sortorder == L_SORT_INCREASING &&
                     array[j] > array[j + gap]) ||
                    (sortorder == L_SORT_DECREASING &&
                     array[j] < array[j + gap]))
                {
                    tmp = array[j];
                    array[j] = array[j + gap];
                    array[j + gap] = tmp;
                    tmp = iarray[j];
                    iarray[j] = iarray[j + gap];
                    iarray[j + gap] = tmp;
                }
            }
        }
    }
 
    naisort = numaCreate(n);
    for (i = 0; i < n; i++)
        numaAddNumber(naisort, iarray[i]);
 
    LEPT_FREE(array);
    LEPT_FREE(iarray);
    return naisort;
}
 
 
NUMA *
numaGetBinSortIndex(NUMA    *nas,
                    l_int32  sortorder)
{
l_int32    i, n, isize, ival, imax;
l_float32  minsize, size;
NUMA      *na, *nai, *nad;
L_PTRA    *paindex;
 
    PROCNAME("numaGetBinSortIndex");
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", procName, NULL);
    if (numaGetCount(nas) == 0) {
        L_WARNING("nas is empty\n", procName);
        return numaCreate(1);
    }
    if (sortorder != L_SORT_INCREASING && sortorder != L_SORT_DECREASING)
        return (NUMA *)ERROR_PTR("invalid sort order", procName, NULL);
    numaGetMin(nas, &minsize, NULL);
    if (minsize < 0)
        return (NUMA *)ERROR_PTR("nas has negative numbers", procName, NULL);
    numaGetMax(nas, &size, NULL);
    isize = (l_int32)size;
    if (isize > MaxInitPtraSize - 1) {
        L_ERROR("array too large: %d elements > max size = %d\n",
                procName, isize, MaxInitPtraSize - 1);
        return NULL;
    }
 
        /* Set up a ptra holding numa at indices for which there
         * are values in nas.  Suppose nas has the value 230 at index
         * 7355.  A numa holding the index 7355 is created and stored
         * at the ptra index 230.  If there is another value of 230
         * in nas, its index is added to the same numa (at index 230
         * in the ptra).  When finished, the ptra can be scanned for numa,
         * and the original indices in the nas can be read out.  In this
         * way, the ptra effectively sorts the input numbers in the nas. */
    paindex = ptraCreate(isize + 1);
    n = numaGetCount(nas);
    for (i = 0; i < n; i++) {
        numaGetIValue(nas, i, &ival);
        nai = (NUMA *)ptraGetPtrToItem(paindex, ival);
        if (!nai) {  /* make it; no shifting will occur */
            nai = numaCreate(1);
            ptraInsert(paindex, ival, nai, L_MIN_DOWNSHIFT);
        }
        numaAddNumber(nai, i);
    }
 
        /* Sort by scanning the ptra, extracting numas and pulling
         * the (index into nas) numbers out of each numa, taken
         * successively in requested order. */
    ptraGetMaxIndex(paindex, &imax);
    nad = numaCreate(0);
    if (sortorder == L_SORT_INCREASING) {
        for (i = 0; i <= imax; i++) {
            na = (NUMA *)ptraRemove(paindex, i, L_NO_COMPACTION);
            if (!na) continue;
            numaJoin(nad, na, 0, -1);
            numaDestroy(&na);
        }
    } else {  /* L_SORT_DECREASING */
        for (i = imax; i >= 0; i--) {
            na = (NUMA *)ptraRemoveLast(paindex);
            if (!na) break;  /* they've all been removed */
            numaJoin(nad, na, 0, -1);
            numaDestroy(&na);
        }
    }
 
    ptraDestroy(&paindex, FALSE, FALSE);
    return nad;
}
 
 
NUMA *
numaSortByIndex(NUMA  *nas,
                NUMA  *naindex)
{
l_int32    i, n, ni, index;
l_float32  val;
NUMA      *nad;
 
    PROCNAME("numaSortByIndex");
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", procName, NULL);
    if (!naindex)
        return (NUMA *)ERROR_PTR("naindex not defined", procName, NULL);
    n = numaGetCount(nas);
    ni = numaGetCount(naindex);
    if (n != ni)
        return (NUMA *)ERROR_PTR("numa sizes differ", procName, NULL);
    if (n == 0) {
        L_WARNING("nas is empty\n", procName);
        return numaCopy(nas);
    }
 
    nad = numaCreate(n);
    for (i = 0; i < n; i++) {
        numaGetIValue(naindex, i, &index);
        numaGetFValue(nas, index, &val);
        numaAddNumber(nad, val);
    }
 
    return nad;
}
 
 
l_int32
numaIsSorted(NUMA     *nas,
             l_int32   sortorder,
             l_int32  *psorted)
{
l_int32    i, n;
l_float32  prevval, val;
 
    PROCNAME("numaIsSorted");
 
    if (!psorted)
        return ERROR_INT("&sorted not defined", procName, 1);
    *psorted = FALSE;
    if (!nas)
        return ERROR_INT("nas not defined", procName, 1);
    if ((n = numaGetCount(nas))== 0) {
        L_WARNING("nas is empty\n", procName);
        *psorted = TRUE;
        return 0;
    }
    if (sortorder != L_SORT_INCREASING && sortorder != L_SORT_DECREASING)
        return ERROR_INT("invalid sortorder", procName, 1);
 
    n = numaGetCount(nas);
    numaGetFValue(nas, 0, &prevval);
    for (i = 1; i < n; i++) {
        numaGetFValue(nas, i, &val);
        if ((sortorder == L_SORT_INCREASING && val < prevval) ||
            (sortorder == L_SORT_DECREASING && val > prevval))
            return 0;
    }
 
    *psorted = TRUE;
    return 0;
}
 
 
l_ok
numaSortPair(NUMA    *nax,
             NUMA    *nay,
             l_int32  sortorder,
             NUMA   **pnasx,
             NUMA   **pnasy)
{
l_int32  sorted;
NUMA    *naindex;
 
    PROCNAME("numaSortPair");
 
    if (pnasx) *pnasx = NULL;
    if (pnasy) *pnasy = NULL;
    if (!pnasx || !pnasy)
        return ERROR_INT("&nasx and/or &nasy not defined", procName, 1);
    if (!nax)
        return ERROR_INT("nax not defined", procName, 1);
    if (!nay)
        return ERROR_INT("nay not defined", procName, 1);
    if (sortorder != L_SORT_INCREASING && sortorder != L_SORT_DECREASING)
        return ERROR_INT("invalid sortorder", procName, 1);
 
    numaIsSorted(nax, sortorder, &sorted);
    if (sorted == TRUE) {
        *pnasx = numaCopy(nax);
        *pnasy = numaCopy(nay);
    } else {
        naindex = numaGetSortIndex(nax, sortorder);
        *pnasx = numaSortByIndex(nax, naindex);
        *pnasy = numaSortByIndex(nay, naindex);
        numaDestroy(&naindex);
    }
 
    return 0;
}
 
 
NUMA *
numaInvertMap(NUMA  *nas)
{
l_int32   i, n, val, error;
l_int32  *test;
NUMA     *nad;
 
    PROCNAME("numaInvertMap");
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", procName, NULL);
    if ((n = numaGetCount(nas)) == 0) {
        L_WARNING("nas is empty\n", procName);
        return numaCopy(nas);
    }
 
    nad = numaMakeConstant(0.0, n);
    test = (l_int32 *)LEPT_CALLOC(n, sizeof(l_int32));
    error = 0;
    for (i = 0; i < n; i++) {
        numaGetIValue(nas, i, &val);
        if (val >= n) {
            error = 1;
            break;
        }
        numaReplaceNumber(nad, val, i);
        if (test[val] == 0) {
            test[val] = 1;
        } else {
            error = 1;
            break;
        }
    }
 
    LEPT_FREE(test);
    if (error) {
        numaDestroy(&nad);
        return (NUMA *)ERROR_PTR("nas not invertible", procName, NULL);
    }
 
    return nad;
}
 
l_ok
numaAddSorted(NUMA      *na,
              l_float32  val)
{
l_int32  index;
 
    PROCNAME("numaAddSorted");
 
    if (!na)
        return ERROR_INT("na not defined", procName, 1);
 
    if (numaFindSortedLoc(na, val, &index) == 1)
        return ERROR_INT("insert failure", procName, 1);
    numaInsertNumber(na, index, val);
    return 0;
}
 
 
l_ok
numaFindSortedLoc(NUMA      *na,
                  l_float32  val,
                  l_int32   *pindex)
{
l_int32    n, increasing, lindex, rindex, midindex;
l_float32  val0, valn, valmid;
 
    PROCNAME("numaFindSortedLoc");
 
    if (!pindex)
        return ERROR_INT("&index not defined", procName, 1);
    *pindex = 0;
    if (!na)
        return ERROR_INT("na not defined", procName, 1);
 
    n = numaGetCount(na);
    if (n == 0) return 0;
    numaGetFValue(na, 0, &val0);
    if (n == 1) {  /* use increasing sort order */
        if (val >= val0)
            *pindex = 1;
        return 0;
    }
 
        /* -----------------  n >= 2 ----------------- */
    numaGetFValue(na, n - 1, &valn);
    increasing = (valn >= val0) ? 1 : 0;  /* sort order */
 
        /* Check if outside bounds of existing array */
    if (increasing) {
        if (val < val0) {
            *pindex = 0;
            return 0;
        } else if (val > valn) {
            *pindex = n;
            return 0;
        }
    } else {  /* decreasing */
        if (val > val0) {
            *pindex = 0;
            return 0;
        } else if (val < valn) {
            *pindex = n;
            return 0;
        }
    }
 
        /* Within bounds of existing array; search */
    lindex = 0;
    rindex = n - 1;
    while (1) {
        midindex = (lindex + rindex) / 2;
        if (midindex == lindex || midindex == rindex) break;
        numaGetFValue(na, midindex, &valmid);
        if (increasing) {
            if (val > valmid)
                lindex = midindex;
            else
                rindex = midindex;
        } else {  /* decreasing */
            if (val > valmid)
                rindex = midindex;
            else
                lindex = midindex;
        }
    }
    *pindex = rindex;
    return 0;
}
 
 
/*----------------------------------------------------------------------*
 *                          Random permutation                          *
 *----------------------------------------------------------------------*/
NUMA *
numaPseudorandomSequence(l_int32  size,
                         l_int32  seed)
{
l_int32   i, index, temp;
l_int32  *array;
NUMA     *na;
 
    PROCNAME("numaPseudorandomSequence");
 
    if (size <= 0)
        return (NUMA *)ERROR_PTR("size <= 0", procName, NULL);
 
    if ((array = (l_int32 *)LEPT_CALLOC(size, sizeof(l_int32))) == NULL)
        return (NUMA *)ERROR_PTR("array not made", procName, NULL);
    for (i = 0; i < size; i++)
        array[i] = i;
    srand(seed);
    for (i = size - 1; i > 0; i--) {
        index = (l_int32)((i + 1) * ((l_float64)rand() / (l_float64)RAND_MAX));
        index = L_MIN(index, i);
        temp = array[i];
        array[i] = array[index];
        array[index] = temp;
    }
 
    na = numaCreateFromIArray(array, size);
    LEPT_FREE(array);
    return na;
}
 
 
NUMA *
numaRandomPermutation(NUMA    *nas,
                      l_int32  seed)
{
l_int32    i, index, size;
l_float32  val;
NUMA      *naindex, *nad;
 
    PROCNAME("numaRandomPermutation");
 
    if (!nas)
        return (NUMA *)ERROR_PTR("nas not defined", procName, NULL);
    if ((size = numaGetCount(nas)) == 0) {
        L_WARNING("nas is empty\n", procName);
        return numaCopy(nas);
    }
 
    naindex = numaPseudorandomSequence(size, seed);
    nad = numaCreate(size);
    for (i = 0; i < size; i++) {
        numaGetIValue(naindex, i, &index);
        numaGetFValue(nas, index, &val);
        numaAddNumber(nad, val);
    }
    numaDestroy(&naindex);
    return nad;
}
 
 
/*----------------------------------------------------------------------*
 *                     Functions requiring sorting                      *
 *----------------------------------------------------------------------*/
l_ok
numaGetRankValue(NUMA       *na,
                 l_float32   fract,
                 NUMA       *nasort,
                 l_int32     usebins,
                 l_float32  *pval)
{
l_int32  n, index;
NUMA    *nas;
 
    PROCNAME("numaGetRankValue");
 
    if (!pval)
        return ERROR_INT("&val not defined", procName, 1);
    *pval = 0.0;  /* init */
    if (!na)
        return ERROR_INT("na not defined", procName, 1);
    if ((n = numaGetCount(na)) == 0)
        return ERROR_INT("na empty", procName, 1);
    if (fract < 0.0 || fract > 1.0)
        return ERROR_INT("fract not in [0.0 ... 1.0]", procName, 1);
 
    if (nasort) {
        nas = nasort;
    } else {
        if (usebins == 0)
            nas = numaSort(NULL, na, L_SORT_INCREASING);
        else
            nas = numaBinSort(na, L_SORT_INCREASING);
        if (!nas)
            return ERROR_INT("nas not made", procName, 1);
    }
    index = (l_int32)(fract * (l_float32)(n - 1) + 0.5);
    numaGetFValue(nas, index, pval);
 
    if (!nasort) numaDestroy(&nas);
    return 0;
}
 
 
l_ok
numaGetMedian(NUMA       *na,
              l_float32  *pval)
{
    PROCNAME("numaGetMedian");
 
    if (!pval)
        return ERROR_INT("&val not defined", procName, 1);
    *pval = 0.0;  /* init */
    if (!na || numaGetCount(na) == 0)
        return ERROR_INT("na not defined or empty", procName, 1);
 
    return numaGetRankValue(na, 0.5, NULL, 0, pval);
}
 
 
l_ok
numaGetBinnedMedian(NUMA     *na,
                    l_int32  *pval)
{
l_int32    ret;
l_float32  fval;
 
    PROCNAME("numaGetBinnedMedian");
 
    if (!pval)
        return ERROR_INT("&val not defined", procName, 1);
    *pval = 0;  /* init */
    if (!na || numaGetCount(na) == 0)
        return ERROR_INT("na not defined or empty", procName, 1);
 
    ret = numaGetRankValue(na, 0.5, NULL, 1, &fval);
    *pval = lept_roundftoi(fval);
    return ret;
}
 
 
l_ok
numaGetMeanDevFromMedian(NUMA       *na,
                         l_float32   med,
                         l_float32  *pdev)
{
l_int32    i, n;
l_float32  val, dev;
 
    PROCNAME("numaGetMeanDevFromMedian");
 
    if (!pdev)
        return ERROR_INT("&dev not defined", procName, 1);
    *pdev = 0.0;  /* init */
    if (!na)
        return ERROR_INT("na not defined", procName, 1);
    if ((n = numaGetCount(na)) == 0)
        return ERROR_INT("na is empty", procName, 1);
 
    dev = 0.0;
    for (i = 0; i < n; i++) {
        numaGetFValue(na, i, &val);
        dev += L_ABS(val - med);
    }
    *pdev = dev / (l_float32)n;
    return 0;
}
 
 
l_ok
numaGetMedianDevFromMedian(NUMA       *na,
                           l_float32  *pmed,
                           l_float32  *pdev)
{
l_int32    n, i;
l_float32  val, med;
NUMA      *nadev;
 
    PROCNAME("numaGetMedianDevFromMedian");
 
    if (pmed) *pmed = 0.0;
    if (!pdev)
        return ERROR_INT("&dev not defined", procName, 1);
    *pdev = 0.0;
    if (!na || numaGetCount(na) == 0)
        return ERROR_INT("na not defined or empty", procName, 1);
 
    numaGetMedian(na, &med);
    if (pmed) *pmed = med;
    n = numaGetCount(na);
    nadev = numaCreate(n);
    for (i = 0; i < n; i++) {
        numaGetFValue(na, i, &val);
        numaAddNumber(nadev, L_ABS(val - med));
    }
    numaGetMedian(nadev, pdev);
 
    numaDestroy(&nadev);
    return 0;
}
 
 
l_ok
numaGetMode(NUMA       *na,
            l_float32  *pval,
            l_int32    *pcount)
{
l_int32     i, n, maxcount, prevcount;
l_float32   val, maxval, prevval;
l_float32  *array;
NUMA       *nasort;
 
    PROCNAME("numaGetMode");
 
    if (pcount) *pcount = 0;
    if (!pval)
        return ERROR_INT("&val not defined", procName, 1);
    *pval = 0.0;
    if (!na)
        return ERROR_INT("na not defined", procName, 1);
    if ((n = numaGetCount(na)) == 0)
        return ERROR_INT("na is empty", procName, 1);
 
    if ((nasort = numaSort(NULL, na, L_SORT_DECREASING)) == NULL)
        return ERROR_INT("nas not made", procName, 1);
    array = numaGetFArray(nasort, L_NOCOPY);
 
        /* Initialize with array[0] */
    prevval = array[0];
    prevcount = 1;
    maxval = prevval;
    maxcount = prevcount;
 
        /* Scan the sorted array, aggregating duplicates */
    for (i = 1; i < n; i++) {
        val = array[i];
        if (val == prevval) {
            prevcount++;
        } else {  /* new value */
            if (prevcount > maxcount) {  /* new max */
                maxcount = prevcount;
                maxval = prevval;
            }
            prevval = val;
            prevcount = 1;
        }
    }
 
        /* Was the mode the last run of elements? */
    if (prevcount > maxcount) {
        maxcount = prevcount;
        maxval = prevval;
    }
 
    *pval = maxval;
    if (pcount)
        *pcount = maxcount;
 
    numaDestroy(&nasort);
    return 0;
}
 
 
/*----------------------------------------------------------------------*
 *                            Rearrangements                            *
 *----------------------------------------------------------------------*/
l_ok
numaJoin(NUMA    *nad,
         NUMA    *nas,
         l_int32  istart,
         l_int32  iend)
{
l_int32    n, i;
l_float32  val;
 
    PROCNAME("numaJoin");
 
    if (!nad)
        return ERROR_INT("nad not defined", procName, 1);
    if (!nas)
        return 0;
 
    if (istart < 0)
        istart = 0;
    n = numaGetCount(nas);
    if (iend < 0 || iend >= n)
        iend = n - 1;
    if (istart > iend)
        return ERROR_INT("istart > iend; nothing to add", procName, 1);
 
    for (i = istart; i <= iend; i++) {
        numaGetFValue(nas, i, &val);
        numaAddNumber(nad, val);
    }
 
    return 0;
}
 
 
l_ok
numaaJoin(NUMAA   *naad,
          NUMAA   *naas,
          l_int32  istart,
          l_int32  iend)
{
l_int32  n, i;
NUMA    *na;
 
    PROCNAME("numaaJoin");
 
    if (!naad)
        return ERROR_INT("naad not defined", procName, 1);
    if (!naas)
        return 0;
 
    if (istart < 0)
        istart = 0;
    n = numaaGetCount(naas);
    if (iend < 0 || iend >= n)
        iend = n - 1;
    if (istart > iend)
        return ERROR_INT("istart > iend; nothing to add", procName, 1);
 
    for (i = istart; i <= iend; i++) {
        na = numaaGetNuma(naas, i, L_CLONE);
        numaaAddNuma(naad, na, L_INSERT);
    }
 
    return 0;
}
 
 
NUMA *
numaaFlattenToNuma(NUMAA  *naa)
{
l_int32  i, nalloc;
NUMA    *na, *nad;
NUMA   **array;
 
    PROCNAME("numaaFlattenToNuma");
 
    if (!naa)
        return (NUMA *)ERROR_PTR("naa not defined", procName, NULL);
 
    nalloc = naa->nalloc;
    array = numaaGetPtrArray(naa);
    nad = numaCreate(0);
    for (i = 0; i < nalloc; i++) {
        na = array[i];
        if (!na) continue;
        numaJoin(nad, na, 0, -1);
    }
 
    return nad;
}