diff options
Diffstat (limited to 'leptonica/src/numafunc1.c')
| -rw-r--r-- | leptonica/src/numafunc1.c | 3697 |
1 files changed, 3697 insertions, 0 deletions
diff --git a/leptonica/src/numafunc1.c b/leptonica/src/numafunc1.c new file mode 100644 index 00000000..acff5559 --- /dev/null +++ b/leptonica/src/numafunc1.c @@ -0,0 +1,3697 @@ +/*====================================================================* + - Copyright (C) 2001 Leptonica. All rights reserved. + - + - Redistribution and use in source and binary forms, with or without + - modification, are permitted provided that the following conditions + - are met: + - 1. Redistributions of source code must retain the above copyright + - notice, this list of conditions and the following disclaimer. + - 2. Redistributions in binary form must reproduce the above + - copyright notice, this list of conditions and the following + - disclaimer in the documentation and/or other materials + - provided with the distribution. + - + - THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + - ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + - LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + - A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL ANY + - CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + - EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + - PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + - PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY + - OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + - NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + - SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + *====================================================================*/ + +/*! + * \file numafunc1.c + * <pre> + * + * -------------------------------------- + * This file has these Numa utilities: + * - arithmetic operations + * - simple data analysis + * - generation of special sequences + * - permutations + * - interpolation + * - sorting + * - data analysis requiring sorting + * - joins and rearrangements + * -------------------------------------- + * + * Arithmetic and logic + * NUMA *numaArithOp() + * NUMA *numaLogicalOp() + * NUMA *numaInvert() + * l_int32 numaSimilar() + * l_int32 numaAddToNumber() + * + * Simple extractions + * l_int32 numaGetMin() + * l_int32 numaGetMax() + * l_int32 numaGetSum() + * NUMA *numaGetPartialSums() + * l_int32 numaGetSumOnInterval() + * l_int32 numaHasOnlyIntegers() + * NUMA *numaSubsample() + * NUMA *numaMakeDelta() + * NUMA *numaMakeSequence() + * NUMA *numaMakeConstant() + * NUMA *numaMakeAbsValue() + * NUMA *numaAddBorder() + * NUMA *numaAddSpecifiedBorder() + * NUMA *numaRemoveBorder() + * l_int32 numaCountNonzeroRuns() + * l_int32 numaGetNonzeroRange() + * l_int32 numaGetCountRelativeToZero() + * NUMA *numaClipToInterval() + * NUMA *numaMakeThresholdIndicator() + * NUMA *numaUniformSampling() + * NUMA *numaReverse() + * + * Signal feature extraction + * NUMA *numaLowPassIntervals() + * NUMA *numaThresholdEdges() + * NUMA *numaGetSpanValues() + * NUMA *numaGetEdgeValues() + * + * Interpolation + * l_int32 numaInterpolateEqxVal() + * l_int32 numaInterpolateEqxInterval() + * l_int32 numaInterpolateArbxVal() + * l_int32 numaInterpolateArbxInterval() + * + * Functions requiring interpolation + * l_int32 numaFitMax() + * l_int32 numaDifferentiateInterval() + * l_int32 numaIntegrateInterval() + * + * Sorting + * NUMA *numaSortGeneral() + * NUMA *numaSortAutoSelect() + * NUMA *numaSortIndexAutoSelect() + * l_int32 numaChooseSortType() + * NUMA *numaSort() + * NUMA *numaBinSort() + * NUMA *numaGetSortIndex() + * NUMA *numaGetBinSortIndex() + * NUMA *numaSortByIndex() + * l_int32 numaIsSorted() + * l_int32 numaSortPair() + * NUMA *numaInvertMap() + * l_int32 numaAddSorted() + * l_int32 numaFindSortedLoc() + * + * Random permutation + * NUMA *numaPseudorandomSequence() + * NUMA *numaRandomPermutation() + * + * Functions requiring sorting + * l_int32 numaGetRankValue() + * l_int32 numaGetMedian() + * l_int32 numaGetBinnedMedian() + * l_int32 numaGetMeanDevFromMedian() + * l_int32 numaGetMedianDevFromMedian() + * l_int32 numaGetMode() + * + * Rearrangements + * l_int32 numaJoin() + * l_int32 numaaJoin() + * NUMA *numaaFlattenToNuma() + * + * Things to remember when using the Numa: + * + * (1) The numa is a struct, not an array. Always use accessors + * (see numabasic.c), never the fields directly. + * + * (2) The number array holds l_float32 values. It can also + * be used to store l_int32 values. See numabasic.c for + * details on using the accessors. + * + * (3) If you use numaCreate(), no numbers are stored and the size is 0. + * You have to add numbers to increase the size. + * If you want to start with a numa of a fixed size, with each + * entry initialized to the same value, use numaMakeConstant(). + * + * (4) Occasionally, in the comments we denote the i-th element of a + * numa by na[i]. This is conceptual only -- the numa is not an array! + * </pre> + */ + +#ifdef HAVE_CONFIG_H +#include <config_auto.h> +#endif /* HAVE_CONFIG_H */ + +#include <math.h> +#include "allheaders.h" + +/*----------------------------------------------------------------------* + * Arithmetic and logical ops on Numas * + *----------------------------------------------------------------------*/ +/*! + * \brief numaArithOp() + * + * \param[in] nad [optional] can be null or equal to na1 (in-place + * \param[in] na1 + * \param[in] na2 + * \param[in] op L_ARITH_ADD, L_ARITH_SUBTRACT, + * L_ARITH_MULTIPLY, L_ARITH_DIVIDE + * \return nad always: operation applied to na1 and na2 + * + * <pre> + * Notes: + * (1) The sizes of na1 and na2 must be equal. + * (2) nad can only null or equal to na1. + * (3) To add a constant to a numa, or to multipy a numa by + * a constant, use numaTransform(). + * </pre> + */ +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; +} + + +/*! + * \brief numaLogicalOp() + * + * \param[in] nad [optional] can be null or equal to na1 (in-place + * \param[in] na1 + * \param[in] na2 + * \param[in] op L_UNION, L_INTERSECTION, L_SUBTRACTION, L_EXCLUSIVE_OR + * \return nad always: operation applied to na1 and na2 + * + * <pre> + * Notes: + * (1) The sizes of na1 and na2 must be equal. + * (2) nad can only be null or equal to na1. + * (3) This is intended for use with indicator arrays (0s and 1s). + * Input data is extracted as integers (0 == false, anything + * else == true); output results are 0 and 1. + * (4) L_SUBTRACTION is subtraction of val2 from val1. For bit logical + * arithmetic this is (val1 & ~val2), but because these values + * are integers, we use (val1 && !val2). + * </pre> + */ +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; +} + + +/*! + * \brief numaInvert() + * + * \param[in] nad [optional] can be null or equal to nas (in-place + * \param[in] nas + * \return nad always: 'inverts' nas + * + * <pre> + * Notes: + * (1) This is intended for use with indicator arrays (0s and 1s). + * It gives a boolean-type output, taking the input as + * an integer and inverting it: + * 0 --> 1 + * anything else --> 0 + * </pre> + */ +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; +} + + +/*! + * \brief numaSimilar() + * + * \param[in] na1 + * \param[in] na2 + * \param[in] maxdiff use 0.0 for exact equality + * \param[out] psimilar 1 if similar; 0 if different + * \return 0 if OK, 1 on error + * + * <pre> + * Notes: + * (1) Float values can differ slightly due to roundoff and + * accumulated errors. Using %maxdiff > 0.0 allows similar + * arrays to be identified. + * </pre> +*/ +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; +} + + +/*! + * \brief numaAddToNumber() + * + * \param[in] na source numa + * \param[in] index element to be changed + * \param[in] val new value to be added + * \return 0 if OK, 1 on error + * + * <pre> + * Notes: + * (1) This is useful for accumulating sums, regardless of the index + * order in which the values are made available. + * (2) Before use, the numa has to be filled up to %index. This would + * typically be used by creating the numa with the full sized + * array, initialized to 0.0, using numaMakeConstant(). + * </pre> + */ +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 * + *----------------------------------------------------------------------*/ +/*! + * \brief numaGetMin() + * + * \param[in] na source numa + * \param[out] pminval [optional] min value + * \param[out] piminloc [optional] index of min location + * \return 0 if OK; 1 on error + */ +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; +} + + +/*! + * \brief numaGetMax() + * + * \param[in] na source numa + * \param[out] pmaxval [optional] max value + * \param[out] pimaxloc [optional] index of max location + * \return 0 if OK; 1 on error + */ +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; +} + + +/*! + * \brief numaGetSum() + * + * \param[in] na source numa + * \param[out] psum sum of values + * \return 0 if OK, 1 on error + */ +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; +} + + +/*! + * \brief numaGetPartialSums() + * + * \param[in] na source numa + * \return nasum, or NULL on error + * + * <pre> + * Notes: + * (1) nasum[i] is the sum for all j <= i of na[j]. + * So nasum[0] = na[0]. + * (2) If you want to generate a rank function, where rank[0] - 0.0, + * insert a 0.0 at the beginning of the nasum array. + * </pre> + */ +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; +} + + +/*! + * \brief numaGetSumOnInterval() + * + * \param[in] na source numa + * \param[in] first beginning index + * \param[in] last final index; use -1 to go to the end + * \param[out] psum sum of values in the index interval range + * \return 0 if OK, 1 on error + */ +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; +} + + +/*! + * \brief numaHasOnlyIntegers() + * + * \param[in] na source numa + * \param[out] pallints 1 if all sampled values are ints; else 0 + * \return 0 if OK, 1 on error + */ +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; +} + + +/*! + * \brief numaSubsample() + * + * \param[in] nas + * \param[in] subfactor subsample factor, >= 1 + * \return nad evenly sampled values from nas, or NULL on error + */ +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; +} + + +/*! + * \brief numaMakeDelta() + * + * \param[in] nas input numa + * \return numa of difference values val[i+1] - val[i], + * or NULL on error + */ +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; +} + + +/*! + * \brief numaMakeSequence() + * + * \param[in] startval + * \param[in] increment + * \param[in] size of sequence + * \return numa of sequence of evenly spaced values, or NULL on error + */ +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; +} + + +/*! + * \brief numaMakeConstant() + * + * \param[in] val + * \param[in] size of numa + * \return numa of given size with all entries equal to 'val', + * or NULL on error + */ +NUMA * +numaMakeConstant(l_float32 val, + l_int32 size) +{ + return numaMakeSequence(val, 0.0, size); +} + + +/*! + * \brief numaMakeAbsValue() + * + * \param[in] nad can be null for new array, or the same as nas for inplace + * \param[in] nas input numa + * \return nad with all numbers being the absval of the input, + * or NULL on error + */ +NUMA * +numaMakeAbsValue(NUMA *nad, + NUMA *nas) +{ +l_int32 i, n; +l_float32 val; + + PROCNAME("numaMakeAbsValue"); + + 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; +} + + +/*! + * \brief numaAddBorder() + * + * \param[in] nas + * \param[in] left number of elements to add before the start + * \param[in] right number of elements to add after the end + * \param[in] val initialize border elements + * \return nad with added elements at left and right, or NULL on error + */ +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; +} + + +/*! + * \brief numaAddSpecifiedBorder() + * + * \param[in] nas + * \param[in] left number of elements to add before the start + * \param[in] right number of elements to add after the end + * \param[in] type L_CONTINUED_BORDER, L_MIRRORED_BORDER + * \return nad with added elements at left and right, or NULL on error + */ +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; +} + + +/*! + * \brief numaRemoveBorder() + * + * \param[in] nas + * \param[in] left number of elements to remove from the start + * \param[in] right number of elements to remove up to the end + * \return nad with removed elements at left and right, or NULL on error + */ +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; +} + + +/*! + * \brief numaCountNonzeroRuns() + * + * \param[in] na e.g., of pixel counts in rows or columns + * \param[out] pcount number of nonzero runs + * \return 0 if OK, 1 on error + */ +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; +} + + +/*! + * \brief numaGetNonzeroRange() + * + * \param[in] na source numa + * \param[in] eps largest value considered to be zero + * \param[out] pfirst, plast interval of array indices + * where values are nonzero + * \return 0 if OK, 1 on error or if no nonzero range is found. + */ +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; +} + + +/*! + * \brief numaGetCountRelativeToZero() + * + * \param[in] na source numa + * \param[in] type L_LESS_THAN_ZERO, L_EQUAL_TO_ZERO, L_GREATER_THAN_ZERO + * \param[out] pcount count of values of given type + * \return 0 if OK, 1 on error + */ +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; +} + + +/*! + * \brief numaClipToInterval() + * + * \param[in] nas + * \param[in] first >= 0; <= last + * \param[in] last + * \return numa with the same values as the input, but clipped + * to the specified interval + * + * <pre> + * Notes: + * If you want the indices of the array values to be unchanged, + * use first = 0. + * Usage: + * This is useful to clip a histogram that has a few nonzero + * values to its nonzero range. + * </pre> + */ +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; +} + + +/*! + * \brief numaMakeThresholdIndicator() + * + * \param[in] nas input numa + * \param[in] thresh threshold value + * \param[in] type L_SELECT_IF_LT, L_SELECT_IF_GT, + * L_SELECT_IF_LTE, L_SELECT_IF_GTE + * \return nad : indicator array: values are 0 and 1 + * + * <pre> + * Notes: + * (1) For each element in nas, if the constraint given by 'type' + * correctly specifies its relation to thresh, a value of 1 + * is recorded in nad. + * </pre> + */ +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; +} + + +/*! + * \brief numaUniformSampling() + * + * \param[in] nas input numa + * \param[in] nsamp number of samples + * \return nad : resampled array, or NULL on error + * + * <pre> + * Notes: + * (1) This resamples the values in the array, using %nsamp + * equal divisions. + * </pre> + */ +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; +} + + +/*! + * \brief numaReverse() + * + * \param[in] nad [optional] can be null or equal to nas + * \param[in] nas input numa + * \return nad : reversed, or NULL on error + * + * <pre> + * Notes: + * (1) Usage: + * numaReverse(nas, nas); // in-place + * nad = numaReverse(NULL, nas); // makes a new one + * </pre> + */ +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 * + *----------------------------------------------------------------------*/ +/*! + * \brief numaLowPassIntervals() + * + * \param[in] nas input numa + * \param[in] thresh threshold fraction of max; in [0.0 ... 1.0] + * \param[in] maxn for normalizing; set maxn = 0.0 to use the max in nas + * \return nad : interval abscissa pairs, or NULL on error + * + * <pre> + * Notes: + * (1) For each interval where the value is less than a specified + * fraction of the maximum, this records the left and right "x" + * value. + * </pre> + */ +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; +} + + +/*! + * \brief numaThresholdEdges() + * + * \param[in] nas input numa + * \param[in] thresh1 low threshold as fraction of max; in [0.0 ... 1.0] + * \param[in] thresh2 high threshold as fraction of max; in [0.0 ... 1.0] + * \param[in] maxn for normalizing; set maxn = 0.0 to use the max in nas + * \return nad edge interval triplets, or NULL on error + * + * <pre> + * Notes: + * (1) For each edge interval, where where the value is less + * than %thresh1 on one side, greater than %thresh2 on + * the other, and between these thresholds throughout the + * interval, this records a triplet of values: the + * 'left' and 'right' edges, and either +1 or -1, depending + * on whether the edge is rising or falling. + * (2) No assumption is made about the value outside the array, + * so if the value at the array edge is between the threshold + * values, it is not considered part of an edge. We start + * looking for edge intervals only after leaving the thresholded + * band. + * </pre> + */ +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; +} + + +/*! + * \brief numaGetSpanValues() + * + * \param[in] na numa that is output of numaLowPassIntervals() + * \param[in] span span number, zero-based + * \param[out] pstart [optional] location of start of transition + * \param[out] pend [optional] location of end of transition + * \return 0 if OK, 1 on error + */ +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; +} + + +/*! + * \brief numaGetEdgeValues() + * + * \param[in] na numa that is output of numaThresholdEdges() + * \param[in] edge edge number, zero-based + * \param[out] pstart [optional] location of start of transition + * \param[out] pend [optional] location of end of transition + * \param[out] psign [optional] transition sign: +1 is rising, + * -1 is falling + * \return 0 if OK, 1 on error + */ +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 * + *----------------------------------------------------------------------*/ +/*! + * \brief numaInterpolateEqxVal() + * + * \param[in] startx xval corresponding to first element in array + * \param[in] deltax x increment between array elements + * \param[in] nay numa of ordinate values, assumed equally spaced + * \param[in] type L_LINEAR_INTERP, L_QUADRATIC_INTERP + * \param[in] xval + * \param[out] pyval interpolated value + * \return 0 if OK, 1 on error e.g., if xval is outside range + * + * <pre> + * Notes: + * (1) Considering nay as a function of x, the x values + * are equally spaced + * (2) Caller should check for valid return. + * + * For linear Lagrangian interpolation (through 2 data pts): + * y(x) = y1(x-x2)/(x1-x2) + y2(x-x1)/(x2-x1) + * + * For quadratic Lagrangian interpolation (through 3 data pts): + * y(x) = y1(x-x2)(x-x3)/((x1-x2)(x1-x3)) + + * y2(x-x1)(x-x3)/((x2-x1)(x2-x3)) + + * y3(x-x1)(x-x2)/((x3-x1)(x3-x2)) + * + * </pre> + */ +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; +} + + +/*! + * \brief numaInterpolateArbxVal() + * + * \param[in] nax numa of abscissa values + * \param[in] nay numa of ordinate values, corresponding to nax + * \param[in] type L_LINEAR_INTERP, L_QUADRATIC_INTERP + * \param[in] xval + * \param[out] pyval interpolated value + * \return 0 if OK, 1 on error e.g., if xval is outside range + * + * <pre> + * Notes: + * (1) The values in nax must be sorted in increasing order. + * If, additionally, they are equally spaced, you can use + * numaInterpolateEqxVal(). + * (2) Caller should check for valid return. + * (3) Uses lagrangian interpolation. See numaInterpolateEqxVal() + * for formulas. + * </pre> + */ +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; +} + + +/*! + * \brief numaInterpolateEqxInterval() + * + * \param[in] startx xval corresponding to first element in nas + * \param[in] deltax x increment between array elements in nas + * \param[in] nasy numa of ordinate values, assumed equally spaced + * \param[in] type L_LINEAR_INTERP, L_QUADRATIC_INTERP + * \param[in] x0 start value of interval + * \param[in] x1 end value of interval + * \param[in] npts number of points to evaluate function in interval + * \param[out] pnax [optional] array of x values in interval + * \param[out] pnay array of y values in interval + * \return 0 if OK, 1 on error + * + * <pre> + * Notes: + * (1) Considering nasy as a function of x, the x values + * are equally spaced. + * (2) This creates nay (and optionally nax) of interpolated + * values over the specified interval (x0, x1). + * (3) If the interval (x0, x1) lies partially outside the array + * nasy (as interpreted by startx and deltax), it is an + * error and returns 1. + * (4) Note that deltax is the intrinsic x-increment for the input + * array nasy, whereas delx is the intrinsic x-increment for the + * output interpolated array nay. + * </pre> + */ +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; +} + + +/*! + * \brief numaInterpolateArbxInterval() + * + * \param[in] nax numa of abscissa values + * \param[in] nay numa of ordinate values, corresponding to nax + * \param[in] type L_LINEAR_INTERP, L_QUADRATIC_INTERP + * \param[in] x0 start value of interval + * \param[in] x1 end value of interval + * \param[in] npts number of points to evaluate function in interval + * \param[out] pnadx [optional] array of x values in interval + * \param[out] pnady array of y values in interval + * \return 0 if OK, 1 on error e.g., if x0 or x1 is outside range + * + * <pre> + * Notes: + * (1) The values in nax must be sorted in increasing order. + * If they are not sorted, we do it here, and complain. + * (2) If the values in nax are equally spaced, you can use + * numaInterpolateEqxInterval(). + * (3) Caller should check for valid return. + * (4) We don't call numaInterpolateArbxVal() for each output + * point, because that requires an O(n) search for + * each point. Instead, we do a single O(n) pass through + * nax, saving the indices to be used for each output yval. + * (5) Uses lagrangian interpolation. See numaInterpolateEqxVal() + * for formulas. + * </pre> + */ +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 * + *----------------------------------------------------------------------*/ +/*! + * \brief numaFitMax() + * + * \param[in] na numa of ordinate values, to fit a max to + * \param[out] pmaxval max value + * \param[in] naloc [optional] associated numa of abscissa values + * \param[out] pmaxloc abscissa value that gives max value in na; + * if naloc == null, this is given as an interpolated + * index value + * \return 0 if OK; 1 on error + * + * <pre> + * Notes: + * If %naloc is given, there is no requirement that the + * data points are evenly spaced. Lagrangian interpolation + * handles that. The only requirement is that the + * data points are ordered so that the values in naloc + * are either increasing or decreasing. We test to make + * sure that the sizes of na and naloc are equal, and it + * is assumed that the correspondences %na[i] as a function + * of %naloc[i] are properly arranged for all i. + * + * The formula for Lagrangian interpolation through 3 data pts is: + * y(x) = y1(x-x2)(x-x3)/((x1-x2)(x1-x3)) + + * y2(x-x1)(x-x3)/((x2-x1)(x2-x3)) + + * y3(x-x1)(x-x2)/((x3-x1)(x3-x2)) + * + * Then the derivative, using the constants (c1,c2,c3) defined below, + * is set to 0: + * y'(x) = 2x(c1+c2+c3) - c1(x2+x3) - c2(x1+x3) - c3(x1+x2) = 0 + * </pre> + */ +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; +} + + +/*! + * \brief numaDifferentiateInterval() + * + * \param[in] nax numa of abscissa values + * \param[in] nay numa of ordinate values, corresponding to nax + * \param[in] x0 start value of interval + * \param[in] x1 end value of interval + * \param[in] npts number of points to evaluate function in interval + * \param[out] pnadx [optional] array of x values in interval + * \param[out] pnady array of derivatives in interval + * \return 0 if OK, 1 on error e.g., if x0 or x1 is outside range + * + * <pre> + * Notes: + * (1) The values in nax must be sorted in increasing order. + * If they are not sorted, it is done in the interpolation + * step, and a warning is issued. + * (2) Caller should check for valid return. + * </pre> + */ +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; +} + + +/*! + * \brief numaIntegrateInterval() + * + * \param[in] nax numa of abscissa values + * \param[in] nay numa of ordinate values, corresponding to nax + * \param[in] x0 start value of interval + * \param[in] x1 end value of interval + * \param[in] npts number of points to evaluate function in interval + * \param[out] psum integral of function over interval + * \return 0 if OK, 1 on error e.g., if x0 or x1 is outside range + * + * <pre> + * Notes: + * (1) The values in nax must be sorted in increasing order. + * If they are not sorted, it is done in the interpolation + * step, and a warning is issued. + * (2) Caller should check for valid return. + * </pre> + */ +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 * + *----------------------------------------------------------------------*/ +/*! + * \brief numaSortGeneral() + * + * \param[in] na source numa + * \param[out] pnasort [optional] sorted numa + * \param[out] pnaindex [optional] index of elements in na associated + * with each element of nasort + * \param[out] pnainvert [optional] index of elements in nasort associated + * with each element of na + * \param[in] sortorder L_SORT_INCREASING or L_SORT_DECREASING + * \param[in] sorttype L_SHELL_SORT or L_BIN_SORT + * \return 0 if OK, 1 on error + * + * <pre> + * Notes: + * (1) Sorting can be confusing. Here's an array of five values with + * the results shown for the 3 output arrays. + * + * na nasort naindex nainvert + * ----------------------------------- + * 3 9 2 3 + * 4 6 3 2 + * 9 4 1 0 + * 6 3 0 1 + * 1 1 4 4 + * + * Note that naindex is a LUT into na for the sorted array values, + * and nainvert directly gives the sorted index values for the + * input array. It is useful to view naindex is as a map: + * 0 --> 2 + * 1 --> 3 + * 2 --> 1 + * 3 --> 0 + * 4 --> 4 + * and nainvert, the inverse of this map: + * 0 --> 3 + * 1 --> 2 + * 2 --> 0 + * 3 --> 1 + * 4 --> 4 + * + * We can write these relations symbolically as: + * nasort[i] = na[naindex[i]] + * na[i] = nasort[nainvert[i]] + * </pre> + */ +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; +} + + +/*! + * \brief numaSortAutoSelect() + * + * \param[in] nas + * \param[in] sortorder L_SORT_INCREASING or L_SORT_DECREASING + * \return naout output sorted numa, or NULL on error + * + * <pre> + * Notes: + * (1) This does either a shell sort or a bin sort, depending on + * the number of elements in nas and the dynamic range. + * </pre> + */ +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); +} + + +/*! + * \brief numaSortIndexAutoSelect() + * + * \param[in] nas + * \param[in] sortorder L_SORT_INCREASING or L_SORT_DECREASING + * \return nad indices of nas, sorted by value in nas, or NULL on error + * + * <pre> + * Notes: + * (1) This does either a shell sort or a bin sort, depending on + * the number of elements in nas and the dynamic range. + * </pre> + */ +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); +} + + +/*! + * \brief numaChooseSortType() + * + * \param[in] nas to be sorted + * \return sorttype L_SHELL_SORT or L_BIN_SORT, or UNDEF on error. + * + * <pre> + * Notes: + * (1) This selects either a shell sort or a bin sort, depending on + * the number of elements in nas and the dynamic range. + * (2) If there are negative values in nas, it selects shell sort. + * </pre> + */ +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; +} + + +/*! + * \brief numaSort() + * + * \param[in] naout output numa; can be NULL or equal to nain + * \param[in] nain input numa + * \param[in] sortorder L_SORT_INCREASING or L_SORT_DECREASING + * \return naout output sorted numa, or NULL on error + * + * <pre> + * Notes: + * (1) Set naout = nain for in-place; otherwise, set naout = NULL. + * (2) Source: Shell sort, modified from K&R, 2nd edition, p.62. + * Slow but simple O(n logn) sort. + * </pre> + */ +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; +} + + +/*! + * \brief numaBinSort() + * + * \param[in] nas of non-negative integers with a max that can + * not exceed (MaxInitPtraSize - 1) + * \param[in] sortorder L_SORT_INCREASING or L_SORT_DECREASING + * \return na sorted, or NULL on error + * + * <pre> + * Notes: + * (1) Because this uses a bin sort with buckets of size 1, it + * is not appropriate for sorting either small arrays or + * arrays containing very large integer values. For such + * arrays, use a standard general sort function like + * numaSort(). + * (2) You can use numaSortAutoSelect() to decide which sorting + * method to use. + * </pre> + */ +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; +} + + +/*! + * \brief numaGetSortIndex() + * + * \param[in] na source numa + * \param[in] sortorder L_SORT_INCREASING or L_SORT_DECREASING + * \return na giving an array of indices that would sort + * the input array, or NULL on error + */ +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; +} + + +/*! + * \brief numaGetBinSortIndex() + * + * \param[in] nas of non-negative integers with a max that can + * not exceed (MaxInitPtraSize - 1) + * \param[in] sortorder L_SORT_INCREASING or L_SORT_DECREASING + * \return na sorted, or NULL on error + * + * <pre> + * Notes: + * (1) This creates an array (or lookup table) that contains + * the sorted position of the elements in the input Numa. + * (2) Because it uses a bin sort with buckets of size 1, it + * is not appropriate for sorting either small arrays or + * arrays containing very large integer values. For such + * arrays, use a standard general sort function like + * numaGetSortIndex(). + * (3) You can use numaSortIndexAutoSelect() to decide which + * sorting method to use. + * </pre> + */ +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; +} + + +/*! + * \brief numaSortByIndex() + * + * \param[in] nas + * \param[in] naindex na that maps from the new numa to the input numa + * \return nad sorted, or NULL on error + */ +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; +} + + +/*! + * \brief numaIsSorted() + * + * \param[in] nas + * \param[in] sortorder L_SORT_INCREASING or L_SORT_DECREASING + * \param[out] psorted 1 if sorted; 0 if not + * \return 1 if OK; 0 on error + * + * <pre> + * Notes: + * (1) This is a quick O(n) test if nas is sorted. It is useful + * in situations where the array is likely to be already + * sorted, and a sort operation can be avoided. + * </pre> + */ +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; +} + + +/*! + * \brief numaSortPair() + * + * \param[in] nax, nay input arrays + * \param[in] sortorder L_SORT_INCREASING or L_SORT_DECREASING + * \param[out] pnasx sorted + * \param[out] pnasy sorted exactly in order of nasx + * \return 0 if OK, 1 on error + * + * <pre> + * Notes: + * (1) This function sorts the two input arrays, nax and nay, + * together, using nax as the key for sorting. + * </pre> + */ +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; +} + + +/*! + * \brief numaInvertMap() + * + * \param[in] nas + * \return nad the inverted map, or NULL on error or if not invertible + * + * <pre> + * Notes: + * (1) This requires that nas contain each integer from 0 to n-1. + * The array is typically an index array into a sort or permutation + * of another array. + * </pre> + */ +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; +} + +/*! + * \brief numaAddSorted() + * + * \param[in] na sorted input + * \param[in] val value to be inserted in sorted order + * \return 0 if OK, 1 on error + * + * <pre> + * Notes: + * (1) The input %na is sorted. This function determines the + * sort order of %na and inserts %val into the array. + * </pre> + */ +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; +} + + +/*! + * \brief numaFindSortedLoc() + * + * \param[in] na sorted input + * \param[in] val value to be inserted in sorted order + * \param[out] *ploc index location to insert @val + * \return 0 if OK, 1 on error + * + * <pre> + * Notes: + * (1) The input %na is sorted. This determines the sort order of @na, + * either increasing or decreasing, and does a binary search for the + * location to insert %val into the array. The search is O(log n). + * (2) The index returned is the location to insert into the array. + * The value at the index, and all values to the right, are + * moved to the right (increasing their index location by 1). + * (3) If n is the size of %na, *ploc can be anything in [0 ... n]. + * if *ploc == 0, the value is inserted at the beginning of the + * array; if *ploc == n, it is inserted at the end. + * (4) If the size of %na is 1, insert with an increasing sort. + * </pre> + */ +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 * + *----------------------------------------------------------------------*/ +/*! + * \brief numaPseudorandomSequence() + * + * \param[in] size of sequence + * \param[in] seed for random number generation + * \return na pseudorandom on [0,...,size - 1], or NULL on error + * + * <pre> + * Notes: + * (1) This uses the Durstenfeld shuffle. + * See: http://en.wikipedia.org/wiki/Fisher–Yates_shuffle. + * Result is a pseudorandom permutation of the sequence of integers + * from 0 to size - 1. + * </pre> + */ +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; +} + + +/*! + * \brief numaRandomPermutation() + * + * \param[in] nas input array + * \param[in] seed for random number generation + * \return nas randomly shuffled array, or NULL on error + */ +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 * + *----------------------------------------------------------------------*/ +/*! + * \brief numaGetRankValue() + * + * \param[in] na source numa + * \param[in] fract use 0.0 for smallest, 1.0 for largest + * \param[in] nasort [optional] increasing sorted version of na + * \param[in] usebins 0 for general sort; 1 for bin sort + * \param[out] pval rank val + * \return 0 if OK; 1 on error + * + * <pre> + * Notes: + * (1) Computes the rank value of a number in the %na, which is + * the number that is a fraction %fract from the small + * end of the sorted version of %na. + * (2) If you do this multiple times for different rank values, + * sort the array in advance and use that for %nasort; + * if you're only calling this once, input %nasort == NULL. + * (3) If %usebins == 1, this uses a bin sorting method. + * Use this only where: + * * the numbers are non-negative integers + * * there are over 100 numbers + * * the maximum value is less than about 50,000 + * (4) The advantage of using a bin sort is that it is O(n), + * instead of O(nlogn) for general sort routines. + * </pre> + */ +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; +} + + +/*! + * \brief numaGetMedian() + * + * \param[in] na source numa + * \param[out] pval median value + * \return 0 if OK; 1 on error + * + * <pre> + * Notes: + * (1) Computes the median value of the numbers in the numa, by + * sorting and finding the middle value in the sorted array. + * </pre> + */ +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); +} + + +/*! + * \brief numaGetBinnedMedian() + * + * \param[in] na source numa + * \param[out] pval integer median value + * \return 0 if OK; 1 on error + * + * <pre> + * Notes: + * (1) Computes the median value of the numbers in the numa, + * using bin sort and finding the middle value in the sorted array. + * (2) See numaGetRankValue() for conditions on na for which + * this should be used. Otherwise, use numaGetMedian(). + * </pre> + */ +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; +} + + +/*! + * \brief numaGetMeanDevFromMedian() + * + * \param[in] na source numa + * \param[in] med median value + * \param[out] pdev average absolute value deviation from median value + * \return 0 if OK; 1 on error + */ +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; +} + + +/*! + * \brief numaGetMedianDevFromMedian() + * + * \param[in] na source numa + * \param[out] pmed [optional] median value + * \param[out] pdev median deviation from median val + * \return 0 if OK; 1 on error + * + * <pre> + * Notes: + * (1) Finds the median of the absolute value of the deviation from + * the median value in the array. Why take the absolute value? + * Consider the case where you have values equally distributed + * about both sides of a median value. Without taking the absolute + * value of the differences, you will get 0 for the deviation, + * and this is not useful. + * </pre> + */ +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; +} + + +/*! + * \brief numaGetMode() + * + * \param[in] na source numa + * \param[out] pval mode val + * \param[out] pcount [optional] mode count + * \return 0 if OK; 1 on error + * + * <pre> + * Notes: + * (1) Computes the mode value of the numbers in the numa, by + * sorting and finding the value of the number with the + * largest count. + * (2) Optionally, also returns that count. + * </pre> + */ +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 * + *----------------------------------------------------------------------*/ +/*! + * \brief numaJoin() + * + * \param[in] nad dest numa; add to this one + * \param[in] nas [optional] source numa; add from this one + * \param[in] istart starting index in nas + * \param[in] iend ending index in nas; use -1 to cat all + * \return 0 if OK, 1 on error + * + * <pre> + * Notes: + * (1) istart < 0 is taken to mean 'read from the start' (istart = 0) + * (2) iend < 0 means 'read to the end' + * (3) if nas == NULL, this is a no-op + * </pre> + */ +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; +} + + +/*! + * \brief numaaJoin() + * + * \param[in] naad dest naa; add to this one + * \param[in] naas [optional] source naa; add from this one + * \param[in] istart starting index in nas + * \param[in] iend ending index in naas; use -1 to cat all + * \return 0 if OK, 1 on error + * + * <pre> + * Notes: + * (1) istart < 0 is taken to mean 'read from the start' (istart = 0) + * (2) iend < 0 means 'read to the end' + * (3) if naas == NULL, this is a no-op + * </pre> + */ +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; +} + + +/*! + * \brief numaaFlattenToNuma() + * + * \param[in] naa + * \return numa, or NULL on error + * + * <pre> + * Notes: + * (1) This 'flattens' the Numaa to a Numa, by joining successively + * each Numa in the Numaa. + * (2) It doesn't make any assumptions about the location of the + * Numas in the Numaa array, unlike most Numaa functions. + * (3) It leaves the input Numaa unchanged. + * </pre> + */ +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; +} + |