Coverage Report

Created: 2024-07-27 06:27

/src/leptonica/src/numafunc1.c
Line
Count
Source (jump to first uncovered line)
1
/*====================================================================*
2
 -  Copyright (C) 2001 Leptonica.  All rights reserved.
3
 -
4
 -  Redistribution and use in source and binary forms, with or without
5
 -  modification, are permitted provided that the following conditions
6
 -  are met:
7
 -  1. Redistributions of source code must retain the above copyright
8
 -     notice, this list of conditions and the following disclaimer.
9
 -  2. Redistributions in binary form must reproduce the above
10
 -     copyright notice, this list of conditions and the following
11
 -     disclaimer in the documentation and/or other materials
12
 -     provided with the distribution.
13
 -
14
 -  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
15
 -  ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
16
 -  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
17
 -  A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL ANY
18
 -  CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
19
 -  EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
20
 -  PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
21
 -  PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
22
 -  OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
23
 -  NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24
 -  SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25
 *====================================================================*/
26
27
/*!
28
 * \file  numafunc1.c
29
 * <pre>
30
 *
31
 *      --------------------------------------
32
 *      This file has these Numa utilities:
33
 *         - arithmetic operations
34
 *         - simple data analysis
35
 *         - generation of special sequences
36
 *         - permutations
37
 *         - interpolation
38
 *         - sorting
39
 *         - data analysis requiring sorting
40
 *         - joins and rearrangements
41
 *      --------------------------------------
42
 *
43
 *      Arithmetic and logic
44
 *          NUMA        *numaArithOp()
45
 *          NUMA        *numaLogicalOp()
46
 *          NUMA        *numaInvert()
47
 *          l_int32      numaSimilar()
48
 *          l_int32      numaAddToNumber()
49
 *
50
 *      Simple extractions
51
 *          l_int32      numaGetMin()
52
 *          l_int32      numaGetMax()
53
 *          l_int32      numaGetSum()
54
 *          NUMA        *numaGetPartialSums()
55
 *          l_int32      numaGetSumOnInterval()
56
 *          l_int32      numaHasOnlyIntegers()
57
 *          l_int32      numaGetMean()
58
 *          l_int32      numaGetMeanAbsval()
59
 *          NUMA        *numaSubsample()
60
 *          NUMA        *numaMakeDelta()
61
 *          NUMA        *numaMakeSequence()
62
 *          NUMA        *numaMakeConstant()
63
 *          NUMA        *numaMakeAbsval()
64
 *          NUMA        *numaAddBorder()
65
 *          NUMA        *numaAddSpecifiedBorder()
66
 *          NUMA        *numaRemoveBorder()
67
 *          l_int32      numaCountNonzeroRuns()
68
 *          l_int32      numaGetNonzeroRange()
69
 *          l_int32      numaGetCountRelativeToZero()
70
 *          NUMA        *numaClipToInterval()
71
 *          NUMA        *numaMakeThresholdIndicator()
72
 *          NUMA        *numaUniformSampling()
73
 *          NUMA        *numaReverse()
74
 *
75
 *      Signal feature extraction
76
 *          NUMA        *numaLowPassIntervals()
77
 *          NUMA        *numaThresholdEdges()
78
 *          NUMA        *numaGetSpanValues()
79
 *          NUMA        *numaGetEdgeValues()
80
 *
81
 *      Interpolation
82
 *          l_int32      numaInterpolateEqxVal()
83
 *          l_int32      numaInterpolateEqxInterval()
84
 *          l_int32      numaInterpolateArbxVal()
85
 *          l_int32      numaInterpolateArbxInterval()
86
 *
87
 *      Functions requiring interpolation
88
 *          l_int32      numaFitMax()
89
 *          l_int32      numaDifferentiateInterval()
90
 *          l_int32      numaIntegrateInterval()
91
 *
92
 *      Sorting
93
 *          NUMA        *numaSortGeneral()
94
 *          NUMA        *numaSortAutoSelect()
95
 *          NUMA        *numaSortIndexAutoSelect()
96
 *          l_int32      numaChooseSortType()
97
 *          NUMA        *numaSort()
98
 *          NUMA        *numaBinSort()
99
 *          NUMA        *numaGetSortIndex()
100
 *          NUMA        *numaGetBinSortIndex()
101
 *          NUMA        *numaSortByIndex()
102
 *          l_int32      numaIsSorted()
103
 *          l_int32      numaSortPair()
104
 *          NUMA        *numaInvertMap()
105
 *          l_int32      numaAddSorted()
106
 *          l_int32      numaFindSortedLoc()
107
 *
108
 *      Random permutation
109
 *          NUMA        *numaPseudorandomSequence()
110
 *          NUMA        *numaRandomPermutation()
111
 *
112
 *      Functions requiring sorting
113
 *          l_int32      numaGetRankValue()
114
 *          l_int32      numaGetMedian()
115
 *          l_int32      numaGetBinnedMedian()
116
 *          l_int32      numaGetMeanDevFromMedian()
117
 *          l_int32      numaGetMedianDevFromMedian()
118
 *          l_int32      numaGetMode()
119
 *
120
 *      Rearrangements
121
 *          l_int32      numaJoin()
122
 *          l_int32      numaaJoin()
123
 *          NUMA        *numaaFlattenToNuma()
124
 *
125
 *    Things to remember when using the Numa:
126
 *
127
 *    (1) The numa is a struct, not an array.  Always use accessors
128
 *        (see numabasic.c), never the fields directly.
129
 *
130
 *    (2) The number array holds l_float32 values.  It can also
131
 *        be used to store l_int32 values.  See numabasic.c for
132
 *        details on using the accessors.
133
 *
134
 *    (3) If you use numaCreate(), no numbers are stored and the size is 0.
135
 *        You have to add numbers to increase the size.
136
 *        If you want to start with a numa of a fixed size, with each
137
 *        entry initialized to the same value, use numaMakeConstant().
138
 *
139
 *    (4) Occasionally, in the comments we denote the i-th element of a
140
 *        numa by na[i].  This is conceptual only -- the numa is not an array!
141
 * </pre>
142
 */
143
144
#ifdef HAVE_CONFIG_H
145
#include <config_auto.h>
146
#endif  /* HAVE_CONFIG_H */
147
148
#include <math.h>
149
#include "allheaders.h"
150
#include "array_internal.h"
151
152
/*----------------------------------------------------------------------*
153
 *                Arithmetic and logical ops on Numas                   *
154
 *----------------------------------------------------------------------*/
155
/*!
156
 * \brief   numaArithOp()
157
 *
158
 * \param[in]    nad     [optional] can be null or equal to na1 (in-place
159
 * \param[in]    na1
160
 * \param[in]    na2
161
 * \param[in]    op      L_ARITH_ADD, L_ARITH_SUBTRACT,
162
 *                       L_ARITH_MULTIPLY, L_ARITH_DIVIDE
163
 * \return  nad always: operation applied to na1 and na2
164
 *
165
 * <pre>
166
 * Notes:
167
 *      (1) The sizes of na1 and na2 must be equal.
168
 *      (2) nad can only null or equal to na1.
169
 *      (3) To add a constant to a numa, or to multiply a numa by
170
 *          a constant, use numaTransform().
171
 * </pre>
172
 */
173
NUMA *
174
numaArithOp(NUMA    *nad,
175
            NUMA    *na1,
176
            NUMA    *na2,
177
            l_int32  op)
178
0
{
179
0
l_int32    i, n;
180
0
l_float32  val1, val2;
181
182
0
    if (!na1 || !na2)
183
0
        return (NUMA *)ERROR_PTR("na1, na2 not both defined", __func__, nad);
184
0
    n = numaGetCount(na1);
185
0
    if (n != numaGetCount(na2))
186
0
        return (NUMA *)ERROR_PTR("na1, na2 sizes differ", __func__, nad);
187
0
    if (nad && nad != na1)
188
0
        return (NUMA *)ERROR_PTR("nad defined but not in-place", __func__, nad);
189
0
    if (op != L_ARITH_ADD && op != L_ARITH_SUBTRACT &&
190
0
        op != L_ARITH_MULTIPLY && op != L_ARITH_DIVIDE)
191
0
        return (NUMA *)ERROR_PTR("invalid op", __func__, nad);
192
0
    if (op == L_ARITH_DIVIDE) {
193
0
        for (i = 0; i < n; i++) {
194
0
            numaGetFValue(na2, i, &val2);
195
0
            if (val2 == 0.0)
196
0
                return (NUMA *)ERROR_PTR("na2 has 0 element", __func__, nad);
197
0
        }
198
0
    }
199
200
        /* If nad is not identical to na1, make it an identical copy */
201
0
    if (!nad)
202
0
        nad = numaCopy(na1);
203
204
0
    for (i = 0; i < n; i++) {
205
0
        numaGetFValue(nad, i, &val1);
206
0
        numaGetFValue(na2, i, &val2);
207
0
        switch (op) {
208
0
        case L_ARITH_ADD:
209
0
            numaSetValue(nad, i, val1 + val2);
210
0
            break;
211
0
        case L_ARITH_SUBTRACT:
212
0
            numaSetValue(nad, i, val1 - val2);
213
0
            break;
214
0
        case L_ARITH_MULTIPLY:
215
0
            numaSetValue(nad, i, val1 * val2);
216
0
            break;
217
0
        case L_ARITH_DIVIDE:
218
0
            numaSetValue(nad, i, val1 / val2);
219
0
            break;
220
0
        default:
221
0
            lept_stderr(" Unknown arith op: %d\n", op);
222
0
            return nad;
223
0
        }
224
0
    }
225
226
0
    return nad;
227
0
}
228
229
230
/*!
231
 * \brief   numaLogicalOp()
232
 *
233
 * \param[in]    nad     [optional] can be null or equal to na1 (in-place
234
 * \param[in]    na1
235
 * \param[in]    na2
236
 * \param[in]    op      L_UNION, L_INTERSECTION, L_SUBTRACTION, L_EXCLUSIVE_OR
237
 * \return  nad always: operation applied to na1 and na2
238
 *
239
 * <pre>
240
 * Notes:
241
 *      (1) The sizes of na1 and na2 must be equal.
242
 *      (2) nad can only be null or equal to na1.
243
 *      (3) This is intended for use with indicator arrays (0s and 1s).
244
 *          Input data is extracted as integers (0 == false, anything
245
 *          else == true); output results are 0 and 1.
246
 *      (4) L_SUBTRACTION is subtraction of val2 from val1.  For bit logical
247
 *          arithmetic this is (val1 & ~val2), but because these values
248
 *          are integers, we use (val1 && !val2).
249
 * </pre>
250
 */
251
NUMA *
252
numaLogicalOp(NUMA    *nad,
253
              NUMA    *na1,
254
              NUMA    *na2,
255
              l_int32  op)
256
0
{
257
0
l_int32  i, n, val1, val2, val;
258
259
0
    if (!na1 || !na2)
260
0
        return (NUMA *)ERROR_PTR("na1, na2 not both defined", __func__, nad);
261
0
    n = numaGetCount(na1);
262
0
    if (n != numaGetCount(na2))
263
0
        return (NUMA *)ERROR_PTR("na1, na2 sizes differ", __func__, nad);
264
0
    if (nad && nad != na1)
265
0
        return (NUMA *)ERROR_PTR("nad defined; not in-place", __func__, nad);
266
0
    if (op != L_UNION && op != L_INTERSECTION &&
267
0
        op != L_SUBTRACTION && op != L_EXCLUSIVE_OR)
268
0
        return (NUMA *)ERROR_PTR("invalid op", __func__, nad);
269
270
        /* If nad is not identical to na1, make it an identical copy */
271
0
    if (!nad)
272
0
        nad = numaCopy(na1);
273
274
0
    for (i = 0; i < n; i++) {
275
0
        numaGetIValue(nad, i, &val1);
276
0
        numaGetIValue(na2, i, &val2);
277
0
        val1 = (val1 == 0) ? 0 : 1;
278
0
        val2 = (val2 == 0) ? 0 : 1;
279
0
        switch (op) {
280
0
        case L_UNION:
281
0
            val = (val1 || val2) ? 1 : 0;
282
0
            numaSetValue(nad, i, val);
283
0
            break;
284
0
        case L_INTERSECTION:
285
0
            val = (val1 && val2) ? 1 : 0;
286
0
            numaSetValue(nad, i, val);
287
0
            break;
288
0
        case L_SUBTRACTION:
289
0
            val = (val1 && !val2) ? 1 : 0;
290
0
            numaSetValue(nad, i, val);
291
0
            break;
292
0
        case L_EXCLUSIVE_OR:
293
0
            val = (val1 != val2) ? 1 : 0;
294
0
            numaSetValue(nad, i, val);
295
0
            break;
296
0
        default:
297
0
            lept_stderr(" Unknown logical op: %d\n", op);
298
0
            return nad;
299
0
        }
300
0
    }
301
302
0
    return nad;
303
0
}
304
305
306
/*!
307
 * \brief   numaInvert()
308
 *
309
 * \param[in]    nad    [optional] can be null or equal to nas (in-place
310
 * \param[in]    nas
311
 * \return  nad always: 'inverts' nas
312
 *
313
 * <pre>
314
 * Notes:
315
 *      (1) This is intended for use with indicator arrays (0s and 1s).
316
 *          It gives a boolean-type output, taking the input as
317
 *          an integer and inverting it:
318
 *              0              -->  1
319
 *              anything else  -->   0
320
 * </pre>
321
 */
322
NUMA *
323
numaInvert(NUMA  *nad,
324
           NUMA  *nas)
325
0
{
326
0
l_int32  i, n, val;
327
328
0
    if (!nas)
329
0
        return (NUMA *)ERROR_PTR("nas not defined", __func__, nad);
330
0
    if (nad && nad != nas)
331
0
        return (NUMA *)ERROR_PTR("nad defined; not in-place", __func__, nad);
332
333
0
    if (!nad)
334
0
        nad = numaCopy(nas);
335
0
    n = numaGetCount(nad);
336
0
    for (i = 0; i < n; i++) {
337
0
        numaGetIValue(nad, i, &val);
338
0
        if (!val)
339
0
            val = 1;
340
0
        else
341
0
            val = 0;
342
0
        numaSetValue(nad, i, val);
343
0
    }
344
0
    return nad;
345
0
}
346
347
348
/*!
349
 * \brief   numaSimilar()
350
 *
351
 * \param[in]    na1
352
 * \param[in]    na2
353
 * \param[in]    maxdiff    use 0.0 for exact equality
354
 * \param[out]   psimilar   1 if similar; 0 if different
355
 * \return  0 if OK, 1 on error
356
 *
357
 * <pre>
358
 * Notes:
359
 *      (1) Float values can differ slightly due to roundoff and
360
 *          accumulated errors.  Using %maxdiff > 0.0 allows similar
361
 *          arrays to be identified.
362
 * </pre>
363
*/
364
l_int32
365
numaSimilar(NUMA      *na1,
366
            NUMA      *na2,
367
            l_float32  maxdiff,
368
            l_int32   *psimilar)
369
0
{
370
0
l_int32    i, n;
371
0
l_float32  val1, val2;
372
373
0
    if (!psimilar)
374
0
        return ERROR_INT("&similar not defined", __func__, 1);
375
0
    *psimilar = 0;
376
0
    if (!na1 || !na2)
377
0
        return ERROR_INT("na1 and na2 not both defined", __func__, 1);
378
0
    maxdiff = L_ABS(maxdiff);
379
380
0
    n = numaGetCount(na1);
381
0
    if (n != numaGetCount(na2)) return 0;
382
383
0
    for (i = 0; i < n; i++) {
384
0
        numaGetFValue(na1, i, &val1);
385
0
        numaGetFValue(na2, i, &val2);
386
0
        if (L_ABS(val1 - val2) > maxdiff) return 0;
387
0
    }
388
389
0
    *psimilar = 1;
390
0
    return 0;
391
0
}
392
393
394
/*!
395
 * \brief   numaAddToNumber()
396
 *
397
 * \param[in]    na       source numa
398
 * \param[in]    index    element to be changed
399
 * \param[in]    val      new value to be added
400
 * \return  0 if OK, 1 on error
401
 *
402
 * <pre>
403
 * Notes:
404
 *      (1) This is useful for accumulating sums, regardless of the index
405
 *          order in which the values are made available.
406
 *      (2) Before use, the numa has to be filled up to %index.  This would
407
 *          typically be used by creating the numa with the full sized
408
 *          array, initialized to 0.0, using numaMakeConstant().
409
 * </pre>
410
 */
411
l_ok
412
numaAddToNumber(NUMA      *na,
413
                l_int32    index,
414
                l_float32  val)
415
0
{
416
0
l_int32  n;
417
418
0
    if (!na)
419
0
        return ERROR_INT("na not defined", __func__, 1);
420
0
    if ((n = numaGetCount(na)) == 0)
421
0
        return ERROR_INT("na is empty", __func__, 1);
422
0
    if (index < 0 || index >= n) {
423
0
        L_ERROR("index %d not in [0,...,%d]\n", __func__, index, n - 1);
424
0
        return 1;
425
0
    }
426
427
0
    na->array[index] += val;
428
0
    return 0;
429
0
}
430
431
432
/*----------------------------------------------------------------------*
433
 *                         Simple extractions                           *
434
 *----------------------------------------------------------------------*/
435
/*!
436
 * \brief   numaGetMin()
437
 *
438
 * \param[in]    na        source numa
439
 * \param[out]   pminval   [optional] min value
440
 * \param[out]   piminloc  [optional] index of min location
441
 * \return  0 if OK; 1 on error
442
 */
443
l_ok
444
numaGetMin(NUMA       *na,
445
           l_float32  *pminval,
446
           l_int32    *piminloc)
447
0
{
448
0
l_int32    i, n, iminloc;
449
0
l_float32  val, minval;
450
451
0
    if (!pminval && !piminloc)
452
0
        return ERROR_INT("nothing to do", __func__, 1);
453
0
    if (pminval) *pminval = 0.0;
454
0
    if (piminloc) *piminloc = 0;
455
0
    if (!na)
456
0
        return ERROR_INT("na not defined", __func__, 1);
457
0
    if ((n = numaGetCount(na)) == 0)
458
0
        return ERROR_INT("na is empty", __func__, 1);
459
460
0
    minval = +1000000000.;
461
0
    iminloc = 0;
462
0
    for (i = 0; i < n; i++) {
463
0
        numaGetFValue(na, i, &val);
464
0
        if (val < minval) {
465
0
            minval = val;
466
0
            iminloc = i;
467
0
        }
468
0
    }
469
470
0
    if (pminval) *pminval = minval;
471
0
    if (piminloc) *piminloc = iminloc;
472
0
    return 0;
473
0
}
474
475
476
/*!
477
 * \brief   numaGetMax()
478
 *
479
 * \param[in]    na        source numa
480
 * \param[out]   pmaxval   [optional] max value
481
 * \param[out]   pimaxloc  [optional] index of max location
482
 * \return  0 if OK; 1 on error
483
 */
484
l_ok
485
numaGetMax(NUMA       *na,
486
           l_float32  *pmaxval,
487
           l_int32    *pimaxloc)
488
0
{
489
0
l_int32    i, n, imaxloc;
490
0
l_float32  val, maxval;
491
492
0
    if (!pmaxval && !pimaxloc)
493
0
        return ERROR_INT("nothing to do", __func__, 1);
494
0
    if (pmaxval) *pmaxval = 0.0;
495
0
    if (pimaxloc) *pimaxloc = 0;
496
0
    if (!na)
497
0
        return ERROR_INT("na not defined", __func__, 1);
498
0
    if ((n = numaGetCount(na)) == 0)
499
0
        return ERROR_INT("na is empty", __func__, 1);
500
501
0
    maxval = -1000000000.;
502
0
    imaxloc = 0;
503
0
    for (i = 0; i < n; i++) {
504
0
        numaGetFValue(na, i, &val);
505
0
        if (val > maxval) {
506
0
            maxval = val;
507
0
            imaxloc = i;
508
0
        }
509
0
    }
510
511
0
    if (pmaxval) *pmaxval = maxval;
512
0
    if (pimaxloc) *pimaxloc = imaxloc;
513
0
    return 0;
514
0
}
515
516
517
/*!
518
 * \brief   numaGetSum()
519
 *
520
 * \param[in]    na     source numa
521
 * \param[out]   psum   sum of values
522
 * \return  0 if OK, 1 on error
523
 */
524
l_ok
525
numaGetSum(NUMA       *na,
526
           l_float32  *psum)
527
0
{
528
0
l_int32    i, n;
529
0
l_float32  val, sum;
530
531
0
    if (!psum)
532
0
        return ERROR_INT("&sum not defined", __func__, 1);
533
0
    *psum = 0;
534
0
    if (!na)
535
0
        return ERROR_INT("na not defined", __func__, 1);
536
537
0
    if ((n = numaGetCount(na)) == 0)
538
0
        return ERROR_INT("na is empty", __func__, 1);
539
0
    sum = 0.0;
540
0
    for (i = 0; i < n; i++) {
541
0
        numaGetFValue(na, i, &val);
542
0
        sum += val;
543
0
    }
544
0
    *psum = sum;
545
0
    return 0;
546
0
}
547
548
549
/*!
550
 * \brief   numaGetPartialSums()
551
 *
552
 * \param[in]    na    source numa
553
 * \return  nasum, or NULL on error
554
 *
555
 * <pre>
556
 * Notes:
557
 *      (1) nasum[i] is the sum for all j <= i of na[j].
558
 *          So nasum[0] = na[0].
559
 *      (2) If you want to generate a rank function, where rank[0] - 0.0,
560
 *          insert a 0.0 at the beginning of the nasum array.
561
 * </pre>
562
 */
563
NUMA *
564
numaGetPartialSums(NUMA  *na)
565
0
{
566
0
l_int32    i, n;
567
0
l_float32  val, sum;
568
0
NUMA      *nasum;
569
570
0
    if (!na)
571
0
        return (NUMA *)ERROR_PTR("na not defined", __func__, NULL);
572
573
0
    if ((n = numaGetCount(na)) == 0)
574
0
        L_WARNING("na is empty\n", __func__);
575
0
    nasum = numaCreate(n);
576
0
    sum = 0.0;
577
0
    for (i = 0; i < n; i++) {
578
0
        numaGetFValue(na, i, &val);
579
0
        sum += val;
580
0
        numaAddNumber(nasum, sum);
581
0
    }
582
0
    return nasum;
583
0
}
584
585
586
/*!
587
 * \brief   numaGetSumOnInterval()
588
 *
589
 * \param[in]    na      source numa
590
 * \param[in]    first   beginning index
591
 * \param[in]    last    final index; use -1 to go to the end
592
 * \param[out]   psum    sum of values in the index interval range
593
 * \return  0 if OK, 1 on error
594
 */
595
l_ok
596
numaGetSumOnInterval(NUMA       *na,
597
                     l_int32     first,
598
                     l_int32     last,
599
                     l_float32  *psum)
600
0
{
601
0
l_int32    i, n;
602
0
l_float32  val, sum;
603
604
0
    if (!psum)
605
0
        return ERROR_INT("&sum not defined", __func__, 1);
606
0
    *psum = 0.0;
607
0
    if (!na)
608
0
        return ERROR_INT("na not defined", __func__, 1);
609
610
0
    sum = 0.0;
611
0
    if ((n = numaGetCount(na)) == 0)
612
0
        return ERROR_INT("na is empty", __func__, 1);
613
0
    if (first < 0) first = 0;
614
0
    if (first >= n || last < -1)  /* not an error */
615
0
        return 0;
616
0
    if (last == -1)
617
0
        last = n - 1;
618
0
    last = L_MIN(last, n - 1);
619
620
0
    for (i = first; i <= last; i++) {
621
0
        numaGetFValue(na, i, &val);
622
0
        sum += val;
623
0
    }
624
0
    *psum = sum;
625
0
    return 0;
626
0
}
627
628
629
/*!
630
 * \brief   numaHasOnlyIntegers()
631
 *
632
 * \param[in]    na           source numa
633
 * \param[out]   pallints     1 if all sampled values are ints; else 0
634
 * \return  0 if OK, 1 on error
635
 */
636
l_ok
637
numaHasOnlyIntegers(NUMA     *na,
638
                    l_int32  *pallints)
639
0
{
640
0
l_int32    i, n;
641
0
l_float32  val;
642
643
0
    if (!pallints)
644
0
        return ERROR_INT("&allints not defined", __func__, 1);
645
0
    *pallints = TRUE;
646
0
    if (!na)
647
0
        return ERROR_INT("na not defined", __func__, 1);
648
649
0
    if ((n = numaGetCount(na)) == 0)
650
0
        return ERROR_INT("na is empty", __func__, 1);
651
0
    for (i = 0; i < n; i ++) {
652
0
        numaGetFValue(na, i, &val);
653
0
        if (val != (l_int32)val) {
654
0
            *pallints = FALSE;
655
0
            return 0;
656
0
        }
657
0
    }
658
0
    return 0;
659
0
}
660
661
662
/*!
663
 * \brief   numaGetMean()
664
 *
665
 * \param[in]    na     source numa
666
 * \param[out]   pave   average of values
667
 * \return  0 if OK, 1 on error
668
 */
669
l_ok
670
numaGetMean(NUMA       *na,
671
            l_float32  *pave)
672
0
{
673
0
l_int32    n;
674
0
l_float32  sum;
675
676
0
    if (!pave)
677
0
        return ERROR_INT("&ave not defined", __func__, 1);
678
0
    *pave = 0;
679
0
    if (!na)
680
0
        return ERROR_INT("na not defined", __func__, 1);
681
0
    if ((n = numaGetCount(na)) == 0)
682
0
        return ERROR_INT("na is empty", __func__, 1);
683
684
0
    numaGetSum(na, &sum);
685
0
    *pave = sum / n;
686
0
    return 0;
687
0
}
688
689
/*!
690
 * \brief   numaGetMeanAbsval()
691
 *
692
 * \param[in]    na         source numa
693
 * \param[out]   paveabs    average of absolute values
694
 * \return  0 if OK, 1 on error
695
 */
696
l_ok
697
numaGetMeanAbsval(NUMA       *na,
698
                  l_float32  *paveabs)
699
0
{
700
0
l_int32  n;
701
0
NUMA    *na1;
702
703
0
    if (!paveabs)
704
0
        return ERROR_INT("&aveabs not defined", __func__, 1);
705
0
    *paveabs = 0;
706
0
    if (!na)
707
0
        return ERROR_INT("na not defined", __func__, 1);
708
0
    if ((n = numaGetCount(na)) == 0)
709
0
        return ERROR_INT("na is empty", __func__, 1);
710
711
0
    na1 = numaMakeAbsval(NULL, na);
712
0
    numaGetMean(na1, paveabs);
713
0
    numaDestroy(&na1);
714
0
    return 0;
715
0
}
716
717
718
/*!
719
 * \brief   numaSubsample()
720
 *
721
 * \param[in]    nas
722
 * \param[in]    subfactor    subsample factor, >= 1
723
 * \return  nad evenly sampled values from nas, or NULL on error
724
 */
725
NUMA *
726
numaSubsample(NUMA    *nas,
727
              l_int32  subfactor)
728
0
{
729
0
l_int32    i, n;
730
0
l_float32  val;
731
0
NUMA      *nad;
732
733
0
    if (!nas)
734
0
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
735
0
    if (subfactor < 1)
736
0
        return (NUMA *)ERROR_PTR("subfactor < 1", __func__, NULL);
737
738
0
    nad = numaCreate(0);
739
0
    if ((n = numaGetCount(nas)) == 0)
740
0
        L_WARNING("nas is empty\n", __func__);
741
0
    for (i = 0; i < n; i++) {
742
0
        if (i % subfactor != 0) continue;
743
0
        numaGetFValue(nas, i, &val);
744
0
        numaAddNumber(nad, val);
745
0
    }
746
747
0
    return nad;
748
0
}
749
750
751
/*!
752
 * \brief   numaMakeDelta()
753
 *
754
 * \param[in]    nas    input numa
755
 * \return  numa of difference values val[i+1] - val[i],
756
 *                    or NULL on error
757
 */
758
NUMA *
759
numaMakeDelta(NUMA  *nas)
760
0
{
761
0
l_int32    i, n;
762
0
l_float32  prev, cur;
763
0
NUMA      *nad;
764
765
0
    if (!nas)
766
0
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
767
0
    if ((n = numaGetCount(nas)) < 2) {
768
0
        L_WARNING("n < 2; returning empty numa\n", __func__);
769
0
        return numaCreate(1);
770
0
    }
771
772
0
    nad = numaCreate(n - 1);
773
0
    numaGetFValue(nas, 0, &prev);
774
0
    for (i = 1; i < n; i++) {
775
0
        numaGetFValue(nas, i, &cur);
776
0
        numaAddNumber(nad, cur - prev);
777
0
        prev = cur;
778
0
    }
779
0
    return nad;
780
0
}
781
782
783
/*!
784
 * \brief   numaMakeSequence()
785
 *
786
 * \param[in]    startval
787
 * \param[in]    increment
788
 * \param[in]    size        of sequence
789
 * \return  numa of sequence of evenly spaced values, or NULL on error
790
 */
791
NUMA *
792
numaMakeSequence(l_float32  startval,
793
                 l_float32  increment,
794
                 l_int32    size)
795
0
{
796
0
l_int32    i;
797
0
l_float32  val;
798
0
NUMA      *na;
799
800
0
    if ((na = numaCreate(size)) == NULL)
801
0
        return (NUMA *)ERROR_PTR("na not made", __func__, NULL);
802
803
0
    for (i = 0; i < size; i++) {
804
0
        val = startval + i * increment;
805
0
        numaAddNumber(na, val);
806
0
    }
807
0
    return na;
808
0
}
809
810
811
/*!
812
 * \brief   numaMakeConstant()
813
 *
814
 * \param[in]    val
815
 * \param[in]    size     of numa
816
 * \return  numa of given size with all entries equal to 'val',
817
 *              or NULL on error
818
 */
819
NUMA *
820
numaMakeConstant(l_float32  val,
821
                 l_int32    size)
822
0
{
823
0
    return numaMakeSequence(val, 0.0, size);
824
0
}
825
826
827
/*!
828
 * \brief   numaMakeAbsval()
829
 *
830
 * \param[in]    nad   can be null for new array, or the same as nas for inplace
831
 * \param[in]    nas   input numa
832
 * \return  nad with all numbers being the absval of the input,
833
 *              or NULL on error
834
 */
835
NUMA *
836
numaMakeAbsval(NUMA  *nad,
837
               NUMA  *nas)
838
0
{
839
0
l_int32    i, n;
840
0
l_float32  val;
841
842
0
    if (!nas)
843
0
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
844
0
    if (nad && nad != nas)
845
0
        return (NUMA *)ERROR_PTR("nad and not in-place", __func__, NULL);
846
847
0
    if (!nad)
848
0
        nad = numaCopy(nas);
849
0
    n = numaGetCount(nad);
850
0
    for (i = 0; i < n; i++) {
851
0
        val = nad->array[i];
852
0
        nad->array[i] = L_ABS(val);
853
0
    }
854
855
0
    return nad;
856
0
}
857
858
859
/*!
860
 * \brief   numaAddBorder()
861
 *
862
 * \param[in]    nas
863
 * \param[in]    left    number of elements to add before the start
864
 * \param[in]    right   number of elements to add after the end
865
 * \param[in]    val     initialize border elements
866
 * \return  nad with added elements at left and right, or NULL on error
867
 */
868
NUMA *
869
numaAddBorder(NUMA      *nas,
870
              l_int32    left,
871
              l_int32    right,
872
              l_float32  val)
873
0
{
874
0
l_int32     i, n, len;
875
0
l_float32   startx, delx;
876
0
l_float32  *fas, *fad;
877
0
NUMA       *nad;
878
879
0
    if (!nas)
880
0
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
881
0
    if (left < 0) left = 0;
882
0
    if (right < 0) right = 0;
883
0
    if (left == 0 && right == 0)
884
0
        return numaCopy(nas);
885
886
0
    n = numaGetCount(nas);
887
0
    len = n + left + right;
888
0
    nad = numaMakeConstant(val, len);
889
0
    numaGetParameters(nas, &startx, &delx);
890
0
    numaSetParameters(nad, startx - delx * left, delx);
891
0
    fas = numaGetFArray(nas, L_NOCOPY);
892
0
    fad = numaGetFArray(nad, L_NOCOPY);
893
0
    for (i = 0; i < n; i++)
894
0
        fad[left + i] = fas[i];
895
896
0
    return nad;
897
0
}
898
899
900
/*!
901
 * \brief   numaAddSpecifiedBorder()
902
 *
903
 * \param[in]    nas
904
 * \param[in]    left    number of elements to add before the start
905
 * \param[in]    right   number of elements to add after the end
906
 * \param[in]    type    L_CONTINUED_BORDER, L_MIRRORED_BORDER
907
 * \return  nad with added elements at left and right, or NULL on error
908
 */
909
NUMA *
910
numaAddSpecifiedBorder(NUMA    *nas,
911
                       l_int32  left,
912
                       l_int32  right,
913
                       l_int32  type)
914
0
{
915
0
l_int32     i, n;
916
0
l_float32  *fa;
917
0
NUMA       *nad;
918
919
0
    if (!nas)
920
0
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
921
0
    if (left < 0) left = 0;
922
0
    if (right < 0) right = 0;
923
0
    if (left == 0 && right == 0)
924
0
        return numaCopy(nas);
925
0
    if (type != L_CONTINUED_BORDER && type != L_MIRRORED_BORDER)
926
0
        return (NUMA *)ERROR_PTR("invalid type", __func__, NULL);
927
0
    n = numaGetCount(nas);
928
0
    if (type == L_MIRRORED_BORDER && (left > n || right > n))
929
0
        return (NUMA *)ERROR_PTR("border too large", __func__, NULL);
930
931
0
    nad = numaAddBorder(nas, left, right, 0);
932
0
    n = numaGetCount(nad);
933
0
    fa = numaGetFArray(nad, L_NOCOPY);
934
0
    if (type == L_CONTINUED_BORDER) {
935
0
        for (i = 0; i < left; i++)
936
0
            fa[i] = fa[left];
937
0
        for (i = n - right; i < n; i++)
938
0
            fa[i] = fa[n - right - 1];
939
0
    } else {  /* type == L_MIRRORED_BORDER */
940
0
        for (i = 0; i < left; i++)
941
0
            fa[i] = fa[2 * left - 1 - i];
942
0
        for (i = 0; i < right; i++)
943
0
            fa[n - right + i] = fa[n - right - i - 1];
944
0
    }
945
946
0
    return nad;
947
0
}
948
949
950
/*!
951
 * \brief   numaRemoveBorder()
952
 *
953
 * \param[in]    nas
954
 * \param[in]    left    number of elements to remove from the start
955
 * \param[in]    right   number of elements to remove up to the end
956
 * \return  nad with removed elements at left and right, or NULL on error
957
 */
958
NUMA *
959
numaRemoveBorder(NUMA      *nas,
960
                 l_int32    left,
961
                 l_int32    right)
962
0
{
963
0
l_int32     i, n, len;
964
0
l_float32   startx, delx;
965
0
l_float32  *fas, *fad;
966
0
NUMA       *nad;
967
968
0
    if (!nas)
969
0
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
970
0
    if (left < 0) left = 0;
971
0
    if (right < 0) right = 0;
972
0
    if (left == 0 && right == 0)
973
0
        return numaCopy(nas);
974
975
0
    n = numaGetCount(nas);
976
0
    if ((len = n - left - right) < 0)
977
0
        return (NUMA *)ERROR_PTR("len < 0 after removal", __func__, NULL);
978
0
    nad = numaMakeConstant(0, len);
979
0
    numaGetParameters(nas, &startx, &delx);
980
0
    numaSetParameters(nad, startx + delx * left, delx);
981
0
    fas = numaGetFArray(nas, L_NOCOPY);
982
0
    fad = numaGetFArray(nad, L_NOCOPY);
983
0
    for (i = 0; i < len; i++)
984
0
        fad[i] = fas[left + i];
985
986
0
    return nad;
987
0
}
988
989
990
/*!
991
 * \brief   numaCountNonzeroRuns()
992
 *
993
 * \param[in]    na      e.g., of pixel counts in rows or columns
994
 * \param[out]   pcount  number of nonzero runs
995
 * \return  0 if OK, 1 on error
996
 */
997
l_ok
998
numaCountNonzeroRuns(NUMA     *na,
999
                     l_int32  *pcount)
1000
0
{
1001
0
l_int32  n, i, val, count, inrun;
1002
1003
0
    if (!pcount)
1004
0
        return ERROR_INT("&count not defined", __func__, 1);
1005
0
    *pcount = 0;
1006
0
    if (!na)
1007
0
        return ERROR_INT("na not defined", __func__, 1);
1008
0
    if ((n = numaGetCount(na)) == 0)
1009
0
        return ERROR_INT("na is empty", __func__, 1);
1010
1011
0
    count = 0;
1012
0
    inrun = FALSE;
1013
0
    for (i = 0; i < n; i++) {
1014
0
        numaGetIValue(na, i, &val);
1015
0
        if (!inrun && val > 0) {
1016
0
            count++;
1017
0
            inrun = TRUE;
1018
0
        } else if (inrun && val == 0) {
1019
0
            inrun = FALSE;
1020
0
        }
1021
0
    }
1022
0
    *pcount = count;
1023
0
    return 0;
1024
0
}
1025
1026
1027
/*!
1028
 * \brief   numaGetNonzeroRange()
1029
 *
1030
 * \param[in]    na              source numa
1031
 * \param[in]    eps             largest value considered to be zero
1032
 * \param[out]   pfirst, plast   interval of array indices
1033
 *                               where values are nonzero
1034
 * \return  0 if OK, 1 on error or if no nonzero range is found.
1035
 */
1036
l_ok
1037
numaGetNonzeroRange(NUMA      *na,
1038
                    l_float32  eps,
1039
                    l_int32   *pfirst,
1040
                    l_int32   *plast)
1041
0
{
1042
0
l_int32    n, i, found;
1043
0
l_float32  val;
1044
1045
0
    if (pfirst) *pfirst = 0;
1046
0
    if (plast) *plast = 0;
1047
0
    if (!pfirst || !plast)
1048
0
        return ERROR_INT("pfirst and plast not both defined", __func__, 1);
1049
0
    if (!na)
1050
0
        return ERROR_INT("na not defined", __func__, 1);
1051
0
    if ((n = numaGetCount(na)) == 0)
1052
0
        return ERROR_INT("na is empty", __func__, 1);
1053
1054
0
    found = FALSE;
1055
0
    for (i = 0; i < n; i++) {
1056
0
        numaGetFValue(na, i, &val);
1057
0
        if (val > eps) {
1058
0
            found = TRUE;
1059
0
            break;
1060
0
        }
1061
0
    }
1062
0
    if (!found) {
1063
0
        *pfirst = n - 1;
1064
0
        *plast = 0;
1065
0
        return 1;
1066
0
    }
1067
1068
0
    *pfirst = i;
1069
0
    for (i = n - 1; i >= 0; i--) {
1070
0
        numaGetFValue(na, i, &val);
1071
0
        if (val > eps)
1072
0
            break;
1073
0
    }
1074
0
    *plast = i;
1075
0
    return 0;
1076
0
}
1077
1078
1079
/*!
1080
 * \brief   numaGetCountRelativeToZero()
1081
 *
1082
 * \param[in]    na      source numa
1083
 * \param[in]    type    L_LESS_THAN_ZERO, L_EQUAL_TO_ZERO, L_GREATER_THAN_ZERO
1084
 * \param[out]   pcount  count of values of given type
1085
 * \return  0 if OK, 1 on error
1086
 */
1087
l_ok
1088
numaGetCountRelativeToZero(NUMA     *na,
1089
                           l_int32   type,
1090
                           l_int32  *pcount)
1091
0
{
1092
0
l_int32    n, i, count;
1093
0
l_float32  val;
1094
1095
0
    if (!pcount)
1096
0
        return ERROR_INT("&count not defined", __func__, 1);
1097
0
    *pcount = 0;
1098
0
    if (!na)
1099
0
        return ERROR_INT("na not defined", __func__, 1);
1100
0
    if ((n = numaGetCount(na)) == 0)
1101
0
        return ERROR_INT("na is empty", __func__, 1);
1102
1103
0
    for (i = 0, count = 0; i < n; i++) {
1104
0
        numaGetFValue(na, i, &val);
1105
0
        if (type == L_LESS_THAN_ZERO && val < 0.0)
1106
0
            count++;
1107
0
        else if (type == L_EQUAL_TO_ZERO && val == 0.0)
1108
0
            count++;
1109
0
        else if (type == L_GREATER_THAN_ZERO && val > 0.0)
1110
0
            count++;
1111
0
    }
1112
1113
0
    *pcount = count;
1114
0
    return 0;
1115
0
}
1116
1117
1118
/*!
1119
 * \brief   numaClipToInterval()
1120
 *
1121
 * \param[in]    nas
1122
 * \param[in]    first    >= 0; <= last
1123
 * \param[in]    last
1124
 * \return  numa with the same values as the input, but clipped
1125
 *              to the specified interval
1126
 *
1127
 * <pre>
1128
 * Notes:
1129
 *        If you want the indices of the array values to be unchanged,
1130
 *        use first = 0.
1131
 *  Usage:
1132
 *        This is useful to clip a histogram that has a few nonzero
1133
 *        values to its nonzero range.
1134
 * </pre>
1135
 */
1136
NUMA *
1137
numaClipToInterval(NUMA    *nas,
1138
                   l_int32  first,
1139
                   l_int32  last)
1140
0
{
1141
0
l_int32    n, i;
1142
0
l_float32  val, startx, delx;
1143
0
NUMA      *nad;
1144
1145
0
    if (!nas)
1146
0
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
1147
0
    if ((n = numaGetCount(nas)) == 0)
1148
0
        return (NUMA *)ERROR_PTR("nas is empty", __func__, NULL);
1149
0
    if (first < 0 || first > last)
1150
0
        return (NUMA *)ERROR_PTR("range not valid", __func__, NULL);
1151
0
    if (first >= n)
1152
0
        return (NUMA *)ERROR_PTR("no elements in range", __func__, NULL);
1153
1154
0
    last = L_MIN(last, n - 1);
1155
0
    if ((nad = numaCreate(last - first + 1)) == NULL)
1156
0
        return (NUMA *)ERROR_PTR("nad not made", __func__, NULL);
1157
0
    for (i = first; i <= last; i++) {
1158
0
        numaGetFValue(nas, i, &val);
1159
0
        numaAddNumber(nad, val);
1160
0
    }
1161
0
    numaGetParameters(nas, &startx, &delx);
1162
0
    numaSetParameters(nad, startx + first * delx, delx);
1163
0
    return nad;
1164
0
}
1165
1166
1167
/*!
1168
 * \brief   numaMakeThresholdIndicator()
1169
 *
1170
 * \param[in]    nas      input numa
1171
 * \param[in]    thresh   threshold value
1172
 * \param[in]    type     L_SELECT_IF_LT, L_SELECT_IF_GT,
1173
 *                        L_SELECT_IF_LTE, L_SELECT_IF_GTE
1174
 * \return   nad : indicator array: values are 0 and 1
1175
 *
1176
 * <pre>
1177
 * Notes:
1178
 *      (1) For each element in nas, if the constraint given by 'type'
1179
 *          correctly specifies its relation to thresh, a value of 1
1180
 *          is recorded in nad.
1181
 * </pre>
1182
 */
1183
NUMA *
1184
numaMakeThresholdIndicator(NUMA      *nas,
1185
                           l_float32  thresh,
1186
                           l_int32    type)
1187
0
{
1188
0
l_int32    n, i, ival;
1189
0
l_float32  fval;
1190
0
NUMA      *nai;
1191
1192
0
    if (!nas)
1193
0
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
1194
0
    if ((n = numaGetCount(nas)) == 0)
1195
0
        return (NUMA *)ERROR_PTR("nas is empty", __func__, NULL);
1196
1197
0
    nai = numaCreate(n);
1198
0
    for (i = 0; i < n; i++) {
1199
0
        numaGetFValue(nas, i, &fval);
1200
0
        ival = 0;
1201
0
        switch (type)
1202
0
        {
1203
0
        case L_SELECT_IF_LT:
1204
0
            if (fval < thresh) ival = 1;
1205
0
            break;
1206
0
        case L_SELECT_IF_GT:
1207
0
            if (fval > thresh) ival = 1;
1208
0
            break;
1209
0
        case L_SELECT_IF_LTE:
1210
0
            if (fval <= thresh) ival = 1;
1211
0
            break;
1212
0
        case L_SELECT_IF_GTE:
1213
0
            if (fval >= thresh) ival = 1;
1214
0
            break;
1215
0
        default:
1216
0
            numaDestroy(&nai);
1217
0
            return (NUMA *)ERROR_PTR("invalid type", __func__, NULL);
1218
0
        }
1219
0
        numaAddNumber(nai, ival);
1220
0
    }
1221
1222
0
    return nai;
1223
0
}
1224
1225
1226
/*!
1227
 * \brief   numaUniformSampling()
1228
 *
1229
 * \param[in]    nas     input numa
1230
 * \param[in]    nsamp   number of samples
1231
 * \return  nad : resampled array, or NULL on error
1232
 *
1233
 * <pre>
1234
 * Notes:
1235
 *      (1) This resamples the values in the array, using %nsamp
1236
 *          equal divisions.
1237
 * </pre>
1238
 */
1239
NUMA *
1240
numaUniformSampling(NUMA    *nas,
1241
                    l_int32  nsamp)
1242
0
{
1243
0
l_int32     n, i, j, ileft, iright;
1244
0
l_float32   left, right, binsize, lfract, rfract, sum, startx, delx;
1245
0
l_float32  *array;
1246
0
NUMA       *nad;
1247
1248
0
    if (!nas)
1249
0
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
1250
0
    if ((n = numaGetCount(nas)) == 0)
1251
0
        return (NUMA *)ERROR_PTR("nas is empty", __func__, NULL);
1252
0
    if (nsamp <= 0)
1253
0
        return (NUMA *)ERROR_PTR("nsamp must be > 0", __func__, NULL);
1254
1255
0
    nad = numaCreate(nsamp);
1256
0
    array = numaGetFArray(nas, L_NOCOPY);
1257
0
    binsize = (l_float32)n / (l_float32)nsamp;
1258
0
    numaGetParameters(nas, &startx, &delx);
1259
0
    numaSetParameters(nad, startx, binsize * delx);
1260
0
    left = 0.0;
1261
0
    for (i = 0; i < nsamp; i++) {
1262
0
        sum = 0.0;
1263
0
        right = left + binsize;
1264
0
        ileft = (l_int32)left;
1265
0
        lfract = 1.0f - left + ileft;
1266
0
        if (lfract >= 1.0)  /* on left bin boundary */
1267
0
            lfract = 0.0;
1268
0
        iright = (l_int32)right;
1269
0
        rfract = right - iright;
1270
0
        iright = L_MIN(iright, n - 1);
1271
0
        if (ileft == iright) {  /* both are within the same original sample */
1272
0
            sum += (lfract + rfract - 1.0f) * array[ileft];
1273
0
        } else {
1274
0
            if (lfract > 0.0001)  /* left fraction */
1275
0
                sum += lfract * array[ileft];
1276
0
            if (rfract > 0.0001)  /* right fraction */
1277
0
                sum += rfract * array[iright];
1278
0
            for (j = ileft + 1; j < iright; j++)  /* entire pixels */
1279
0
                sum += array[j];
1280
0
        }
1281
1282
0
        numaAddNumber(nad, sum);
1283
0
        left = right;
1284
0
    }
1285
0
    return nad;
1286
0
}
1287
1288
1289
/*!
1290
 * \brief   numaReverse()
1291
 *
1292
 * \param[in]    nad    [optional] can be null or equal to nas
1293
 * \param[in]    nas    input numa
1294
 * \return  nad : reversed, or NULL on error
1295
 *
1296
 * <pre>
1297
 * Notes:
1298
 *      (1) Usage:
1299
 *            numaReverse(nas, nas);   // in-place
1300
 *            nad = numaReverse(NULL, nas);  // makes a new one
1301
 * </pre>
1302
 */
1303
NUMA *
1304
numaReverse(NUMA  *nad,
1305
            NUMA  *nas)
1306
0
{
1307
0
l_int32    n, i;
1308
0
l_float32  val1, val2;
1309
1310
0
    if (!nas)
1311
0
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
1312
0
    if (nad && nas != nad)
1313
0
        return (NUMA *)ERROR_PTR("nad defined but != nas", __func__, NULL);
1314
1315
0
    n = numaGetCount(nas);
1316
0
    if (nad) {  /* in-place */
1317
0
        for (i = 0; i < n / 2; i++) {
1318
0
            numaGetFValue(nad, i, &val1);
1319
0
            numaGetFValue(nad, n - i - 1, &val2);
1320
0
            numaSetValue(nad, i, val2);
1321
0
            numaSetValue(nad, n - i - 1, val1);
1322
0
        }
1323
0
    } else {
1324
0
        nad = numaCreate(n);
1325
0
        for (i = n - 1; i >= 0; i--) {
1326
0
            numaGetFValue(nas, i, &val1);
1327
0
            numaAddNumber(nad, val1);
1328
0
        }
1329
0
    }
1330
1331
        /* Reverse the startx and delx fields */
1332
0
    nad->startx = nas->startx + (n - 1) * nas->delx;
1333
0
    nad->delx = -nas->delx;
1334
0
    return nad;
1335
0
}
1336
1337
1338
/*----------------------------------------------------------------------*
1339
 *                       Signal feature extraction                      *
1340
 *----------------------------------------------------------------------*/
1341
/*!
1342
 * \brief   numaLowPassIntervals()
1343
 *
1344
 * \param[in]    nas      input numa
1345
 * \param[in]    thresh   threshold fraction of max; in [0.0 ... 1.0]
1346
 * \param[in]    maxn     for normalizing; set maxn = 0.0 to use the max in nas
1347
 * \return  nad : interval abscissa pairs, or NULL on error
1348
 *
1349
 * <pre>
1350
 * Notes:
1351
 *      (1) For each interval where the value is less than a specified
1352
 *          fraction of the maximum, this records the left and right "x"
1353
 *          value.
1354
 * </pre>
1355
 */
1356
NUMA *
1357
numaLowPassIntervals(NUMA      *nas,
1358
                     l_float32  thresh,
1359
                     l_float32  maxn)
1360
0
{
1361
0
l_int32    n, i, inrun;
1362
0
l_float32  maxval, threshval, fval, startx, delx, x0, x1;
1363
0
NUMA      *nad;
1364
1365
0
    if (!nas)
1366
0
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
1367
0
    if ((n = numaGetCount(nas)) == 0)
1368
0
        return (NUMA *)ERROR_PTR("nas is empty", __func__, NULL);
1369
0
    if (thresh < 0.0 || thresh > 1.0)
1370
0
        return (NUMA *)ERROR_PTR("invalid thresh", __func__, NULL);
1371
1372
        /* The input threshold is a fraction of the max.
1373
         * The first entry in nad is the value of the max. */
1374
0
    if (maxn == 0.0)
1375
0
        numaGetMax(nas, &maxval, NULL);
1376
0
    else
1377
0
        maxval = maxn;
1378
0
    numaGetParameters(nas, &startx, &delx);
1379
0
    threshval = thresh * maxval;
1380
0
    nad = numaCreate(0);
1381
0
    numaAddNumber(nad, maxval);
1382
1383
        /* Write pairs of pts (x0, x1) for the intervals */
1384
0
    inrun = FALSE;
1385
0
    for (i = 0; i < n; i++) {
1386
0
        numaGetFValue(nas, i, &fval);
1387
0
        if (fval < threshval && inrun == FALSE) {  /* start a new run */
1388
0
            inrun = TRUE;
1389
0
            x0 = startx + i * delx;
1390
0
        } else if (fval > threshval && inrun == TRUE) {  /* end the run */
1391
0
            inrun = FALSE;
1392
0
            x1 = startx + i * delx;
1393
0
            numaAddNumber(nad, x0);
1394
0
            numaAddNumber(nad, x1);
1395
0
        }
1396
0
    }
1397
0
    if (inrun == TRUE) {  /* must end the last run */
1398
0
        x1 = startx + (n - 1) * delx;
1399
0
        numaAddNumber(nad, x0);
1400
0
        numaAddNumber(nad, x1);
1401
0
    }
1402
1403
0
    return nad;
1404
0
}
1405
1406
1407
/*!
1408
 * \brief   numaThresholdEdges()
1409
 *
1410
 * \param[in]    nas      input numa
1411
 * \param[in]    thresh1  low threshold as fraction of max; in [0.0 ... 1.0]
1412
 * \param[in]    thresh2  high threshold as fraction of max; in [0.0 ... 1.0]
1413
 * \param[in]    maxn     for normalizing; set maxn = 0.0 to use the max in nas
1414
 * \return  nad   edge interval triplets, or NULL on error
1415
 *
1416
 * <pre>
1417
 * Notes:
1418
 *      (1) For each edge interval, where the value is less
1419
 *          than %thresh1 on one side, greater than %thresh2 on
1420
 *          the other, and between these thresholds throughout the
1421
 *          interval, this records a triplet of values: the
1422
 *          'left' and 'right' edges, and either +1 or -1, depending
1423
 *          on whether the edge is rising or falling.
1424
 *      (2) No assumption is made about the value outside the array,
1425
 *          so if the value at the array edge is between the threshold
1426
 *          values, it is not considered part of an edge.  We start
1427
 *          looking for edge intervals only after leaving the thresholded
1428
 *          band.
1429
 * </pre>
1430
 */
1431
NUMA *
1432
numaThresholdEdges(NUMA      *nas,
1433
                   l_float32  thresh1,
1434
                   l_float32  thresh2,
1435
                   l_float32  maxn)
1436
0
{
1437
0
l_int32    n, i, istart, inband, output, sign;
1438
0
l_int32    startbelow, below, above, belowlast, abovelast;
1439
0
l_float32  maxval, threshval1, threshval2, fval, startx, delx, x0, x1;
1440
0
NUMA      *nad;
1441
1442
0
    if (!nas)
1443
0
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
1444
0
    if ((n = numaGetCount(nas)) == 0)
1445
0
        return (NUMA *)ERROR_PTR("nas is empty", __func__, NULL);
1446
0
    if (thresh1 < 0.0 || thresh1 > 1.0 || thresh2 < 0.0 || thresh2 > 1.0)
1447
0
        return (NUMA *)ERROR_PTR("invalid thresholds", __func__, NULL);
1448
0
    if (thresh2 < thresh1)
1449
0
        return (NUMA *)ERROR_PTR("thresh2 < thresh1", __func__, NULL);
1450
1451
        /* The input thresholds are fractions of the max.
1452
         * The first entry in nad is the value of the max used
1453
         * here for normalization. */
1454
0
    if (maxn == 0.0)
1455
0
        numaGetMax(nas, &maxval, NULL);
1456
0
    else
1457
0
        maxval = maxn;
1458
0
    numaGetMax(nas, &maxval, NULL);
1459
0
    numaGetParameters(nas, &startx, &delx);
1460
0
    threshval1 = thresh1 * maxval;
1461
0
    threshval2 = thresh2 * maxval;
1462
0
    nad = numaCreate(0);
1463
0
    numaAddNumber(nad, maxval);
1464
1465
        /* Write triplets of pts (x0, x1, sign) for the edges.
1466
         * First make sure we start search from outside the band.
1467
         * Only one of {belowlast, abovelast} is true. */
1468
0
    for (i = 0; i < n; i++) {
1469
0
        istart = i;
1470
0
        numaGetFValue(nas, i, &fval);
1471
0
        belowlast = (fval < threshval1) ? TRUE : FALSE;
1472
0
        abovelast = (fval > threshval2) ? TRUE : FALSE;
1473
0
        if (belowlast == TRUE || abovelast == TRUE)
1474
0
            break;
1475
0
    }
1476
0
    if (istart == n)  /* no intervals found */
1477
0
        return nad;
1478
1479
        /* x0 and x1 can only be set from outside the edge.
1480
         * They are the values just before entering the band,
1481
         * and just after entering the band.  We can jump through
1482
         * the band, in which case they differ by one index in nas. */
1483
0
    inband = FALSE;
1484
0
    startbelow = belowlast; /* one of these is true */
1485
0
    output = FALSE;
1486
0
    x0 = startx + istart * delx;
1487
0
    for (i = istart + 1; i < n; i++) {
1488
0
        numaGetFValue(nas, i, &fval);
1489
0
        below = (fval < threshval1) ? TRUE : FALSE;
1490
0
        above = (fval > threshval2) ? TRUE : FALSE;
1491
0
        if (!inband && belowlast && above) {  /* full jump up */
1492
0
            x1 = startx + i * delx;
1493
0
            sign = 1;
1494
0
            startbelow = FALSE;  /* for the next transition */
1495
0
            output = TRUE;
1496
0
        } else if (!inband && abovelast && below) {  /* full jump down */
1497
0
            x1 = startx + i * delx;
1498
0
            sign = -1;
1499
0
            startbelow = TRUE;  /* for the next transition */
1500
0
            output = TRUE;
1501
0
        } else if (inband && startbelow && above) {  /* exit rising; success */
1502
0
            x1 = startx + i * delx;
1503
0
            sign = 1;
1504
0
            inband = FALSE;
1505
0
            startbelow = FALSE;  /* for the next transition */
1506
0
            output = TRUE;
1507
0
        } else if (inband && !startbelow && below) {
1508
                /* exit falling; success */
1509
0
            x1 = startx + i * delx;
1510
0
            sign = -1;
1511
0
            inband = FALSE;
1512
0
            startbelow = TRUE;  /* for the next transition */
1513
0
            output = TRUE;
1514
0
        } else if (inband && !startbelow && above) {  /* exit rising; failure */
1515
0
            x0 = startx + i * delx;
1516
0
            inband = FALSE;
1517
0
        } else if (inband && startbelow && below) {  /* exit falling; failure */
1518
0
            x0 = startx + i * delx;
1519
0
            inband = FALSE;
1520
0
        } else if (!inband && !above && !below) {  /* enter */
1521
0
            inband = TRUE;
1522
0
            startbelow = belowlast;
1523
0
        } else if (!inband && (above || below)) {  /* outside and remaining */
1524
0
            x0 = startx + i * delx;  /* update position */
1525
0
        }
1526
0
        belowlast = below;
1527
0
        abovelast = above;
1528
0
        if (output) {  /* we have exited; save new x0 */
1529
0
            numaAddNumber(nad, x0);
1530
0
            numaAddNumber(nad, x1);
1531
0
            numaAddNumber(nad, sign);
1532
0
            output = FALSE;
1533
0
            x0 = startx + i * delx;
1534
0
        }
1535
0
    }
1536
1537
0
    return nad;
1538
0
}
1539
1540
1541
/*!
1542
 * \brief   numaGetSpanValues()
1543
 *
1544
 * \param[in]    na       numa that is output of numaLowPassIntervals()
1545
 * \param[in]    span     span number, zero-based
1546
 * \param[out]   pstart   [optional] location of start of transition
1547
 * \param[out]   pend     [optional] location of end of transition
1548
 * \return  0 if OK, 1 on error
1549
 */
1550
l_int32
1551
numaGetSpanValues(NUMA    *na,
1552
                  l_int32  span,
1553
                  l_int32 *pstart,
1554
                  l_int32 *pend)
1555
0
{
1556
0
l_int32  n, nspans;
1557
1558
0
    if (!na)
1559
0
        return ERROR_INT("na not defined", __func__, 1);
1560
0
    if ((n = numaGetCount(na)) == 0)
1561
0
        return ERROR_INT("na is empty", __func__, 1);
1562
0
    if (n % 2 != 1)
1563
0
        return ERROR_INT("n is not odd", __func__, 1);
1564
0
    nspans = n / 2;
1565
0
    if (nspans < 0 || span >= nspans)
1566
0
        return ERROR_INT("invalid span", __func__, 1);
1567
1568
0
    if (pstart) numaGetIValue(na, 2 * span + 1, pstart);
1569
0
    if (pend) numaGetIValue(na, 2 * span + 2, pend);
1570
0
    return 0;
1571
0
}
1572
1573
1574
/*!
1575
 * \brief   numaGetEdgeValues()
1576
 *
1577
 * \param[in]    na       numa that is output of numaThresholdEdges()
1578
 * \param[in]    edge     edge number, zero-based
1579
 * \param[out]   pstart   [optional] location of start of transition
1580
 * \param[out]   pend     [optional] location of end of transition
1581
 * \param[out]   psign    [optional] transition sign: +1 is rising,
1582
 *                        -1 is falling
1583
 * \return  0 if OK, 1 on error
1584
 */
1585
l_int32
1586
numaGetEdgeValues(NUMA    *na,
1587
                  l_int32  edge,
1588
                  l_int32 *pstart,
1589
                  l_int32 *pend,
1590
                  l_int32 *psign)
1591
0
{
1592
0
l_int32  n, nedges;
1593
1594
0
    if (!na)
1595
0
        return ERROR_INT("na not defined", __func__, 1);
1596
0
    if ((n = numaGetCount(na)) == 0)
1597
0
        return ERROR_INT("na is empty", __func__, 1);
1598
0
    if (n % 3 != 1)
1599
0
        return ERROR_INT("n % 3 is not 1", __func__, 1);
1600
0
    nedges = (n - 1) / 3;
1601
0
    if (edge < 0 || edge >= nedges)
1602
0
        return ERROR_INT("invalid edge", __func__, 1);
1603
1604
0
    if (pstart) numaGetIValue(na, 3 * edge + 1, pstart);
1605
0
    if (pend) numaGetIValue(na, 3 * edge + 2, pend);
1606
0
    if (psign) numaGetIValue(na, 3 * edge + 3, psign);
1607
0
    return 0;
1608
0
}
1609
1610
1611
/*----------------------------------------------------------------------*
1612
 *                             Interpolation                            *
1613
 *----------------------------------------------------------------------*/
1614
/*!
1615
 * \brief   numaInterpolateEqxVal()
1616
 *
1617
 * \param[in]    startx   xval corresponding to first element in array
1618
 * \param[in]    deltax   x increment between array elements
1619
 * \param[in]    nay      numa of ordinate values, assumed equally spaced
1620
 * \param[in]    type     L_LINEAR_INTERP, L_QUADRATIC_INTERP
1621
 * \param[in]    xval
1622
 * \param[out]   pyval    interpolated value
1623
 * \return  0 if OK, 1 on error e.g., if xval is outside range
1624
 *
1625
 * <pre>
1626
 * Notes:
1627
 *      (1) Considering nay as a function of x, the x values
1628
 *          are equally spaced
1629
 *      (2) Caller should check for valid return.
1630
 *
1631
 *  For linear Lagrangian interpolation (through 2 data pts):
1632
 *         y(x) = y1(x-x2)/(x1-x2) + y2(x-x1)/(x2-x1)
1633
 *
1634
 *  For quadratic Lagrangian interpolation (through 3 data pts):
1635
 *         y(x) = y1(x-x2)(x-x3)/((x1-x2)(x1-x3)) +
1636
 *                y2(x-x1)(x-x3)/((x2-x1)(x2-x3)) +
1637
 *                y3(x-x1)(x-x2)/((x3-x1)(x3-x2))
1638
 *
1639
 * </pre>
1640
 */
1641
l_ok
1642
numaInterpolateEqxVal(l_float32   startx,
1643
                      l_float32   deltax,
1644
                      NUMA       *nay,
1645
                      l_int32     type,
1646
                      l_float32   xval,
1647
                      l_float32  *pyval)
1648
0
{
1649
0
l_int32     i, n, i1, i2, i3;
1650
0
l_float32   x1, x2, x3, fy1, fy2, fy3, d1, d2, d3, del, fi, maxx;
1651
0
l_float32  *fa;
1652
1653
0
    if (!pyval)
1654
0
        return ERROR_INT("&yval not defined", __func__, 1);
1655
0
    *pyval = 0.0;
1656
0
    if (!nay)
1657
0
        return ERROR_INT("nay not defined", __func__, 1);
1658
0
    if (deltax <= 0.0)
1659
0
        return ERROR_INT("deltax not > 0", __func__, 1);
1660
0
    if (type != L_LINEAR_INTERP && type != L_QUADRATIC_INTERP)
1661
0
        return ERROR_INT("invalid interp type", __func__, 1);
1662
0
    if ((n = numaGetCount(nay)) < 2)
1663
0
        return ERROR_INT("not enough points", __func__, 1);
1664
0
    if (type == L_QUADRATIC_INTERP && n == 2) {
1665
0
        type = L_LINEAR_INTERP;
1666
0
        L_WARNING("only 2 points; using linear interp\n", __func__);
1667
0
    }
1668
0
    maxx = startx + deltax * (n - 1);
1669
0
    if (xval < startx || xval > maxx)
1670
0
        return ERROR_INT("xval is out of bounds", __func__, 1);
1671
1672
0
    fa = numaGetFArray(nay, L_NOCOPY);
1673
0
    fi = (xval - startx) / deltax;
1674
0
    i = (l_int32)fi;
1675
0
    del = fi - i;
1676
0
    if (del == 0.0) {  /* no interpolation required */
1677
0
        *pyval = fa[i];
1678
0
        return 0;
1679
0
    }
1680
1681
0
    if (type == L_LINEAR_INTERP) {
1682
0
        *pyval = fa[i] + del * (fa[i + 1] - fa[i]);
1683
0
        return 0;
1684
0
    }
1685
1686
        /* Quadratic interpolation */
1687
0
    d1 = d3 = 0.5f / (deltax * deltax);
1688
0
    d2 = -2.f * d1;
1689
0
    if (i == 0) {
1690
0
        i1 = i;
1691
0
        i2 = i + 1;
1692
0
        i3 = i + 2;
1693
0
    } else {
1694
0
        i1 = i - 1;
1695
0
        i2 = i;
1696
0
        i3 = i + 1;
1697
0
    }
1698
0
    x1 = startx + i1 * deltax;
1699
0
    x2 = startx + i2 * deltax;
1700
0
    x3 = startx + i3 * deltax;
1701
0
    fy1 = d1 * fa[i1];
1702
0
    fy2 = d2 * fa[i2];
1703
0
    fy3 = d3 * fa[i3];
1704
0
    *pyval = fy1 * (xval - x2) * (xval - x3) +
1705
0
             fy2 * (xval - x1) * (xval - x3) +
1706
0
             fy3 * (xval - x1) * (xval - x2);
1707
0
    return 0;
1708
0
}
1709
1710
1711
/*!
1712
 * \brief   numaInterpolateArbxVal()
1713
 *
1714
 * \param[in]    nax    numa of abscissa values
1715
 * \param[in]    nay    numa of ordinate values, corresponding to nax
1716
 * \param[in]    type   L_LINEAR_INTERP, L_QUADRATIC_INTERP
1717
 * \param[in]    xval
1718
 * \param[out]   pyval  interpolated value
1719
 * \return  0 if OK, 1 on error e.g., if xval is outside range
1720
 *
1721
 * <pre>
1722
 * Notes:
1723
 *      (1) The values in nax must be sorted in increasing order.
1724
 *          If, additionally, they are equally spaced, you can use
1725
 *          numaInterpolateEqxVal().
1726
 *      (2) Caller should check for valid return.
1727
 *      (3) Uses lagrangian interpolation.  See numaInterpolateEqxVal()
1728
 *          for formulas.
1729
 * </pre>
1730
 */
1731
l_ok
1732
numaInterpolateArbxVal(NUMA       *nax,
1733
                       NUMA       *nay,
1734
                       l_int32     type,
1735
                       l_float32   xval,
1736
                       l_float32  *pyval)
1737
0
{
1738
0
l_int32     i, im, nx, ny, i1, i2, i3;
1739
0
l_float32   delu, dell, fract, d1, d2, d3;
1740
0
l_float32   minx, maxx;
1741
0
l_float32  *fax, *fay;
1742
1743
0
    if (!pyval)
1744
0
        return ERROR_INT("&yval not defined", __func__, 1);
1745
0
    *pyval = 0.0;
1746
0
    if (!nax)
1747
0
        return ERROR_INT("nax not defined", __func__, 1);
1748
0
    if (!nay)
1749
0
        return ERROR_INT("nay not defined", __func__, 1);
1750
0
    if (type != L_LINEAR_INTERP && type != L_QUADRATIC_INTERP)
1751
0
        return ERROR_INT("invalid interp type", __func__, 1);
1752
0
    ny = numaGetCount(nay);
1753
0
    nx = numaGetCount(nax);
1754
0
    if (nx != ny)
1755
0
        return ERROR_INT("nax and nay not same size arrays", __func__, 1);
1756
0
    if (ny < 2)
1757
0
        return ERROR_INT("not enough points", __func__, 1);
1758
0
    if (type == L_QUADRATIC_INTERP && ny == 2) {
1759
0
        type = L_LINEAR_INTERP;
1760
0
        L_WARNING("only 2 points; using linear interp\n", __func__);
1761
0
    }
1762
0
    numaGetFValue(nax, 0, &minx);
1763
0
    numaGetFValue(nax, nx - 1, &maxx);
1764
0
    if (xval < minx || xval > maxx)
1765
0
        return ERROR_INT("xval is out of bounds", __func__, 1);
1766
1767
0
    fax = numaGetFArray(nax, L_NOCOPY);
1768
0
    fay = numaGetFArray(nay, L_NOCOPY);
1769
1770
        /* Linear search for interval.  We are guaranteed
1771
         * to either return or break out of the loop.
1772
         * In addition, we are assured that fax[i] - fax[im] > 0.0 */
1773
0
    if (xval == fax[0]) {
1774
0
        *pyval = fay[0];
1775
0
        return 0;
1776
0
    }
1777
0
    im = 0;
1778
0
    dell = 0.0;
1779
0
    for (i = 1; i < nx; i++) {
1780
0
        delu = fax[i] - xval;
1781
0
        if (delu >= 0.0) {  /* we've passed it */
1782
0
            if (delu == 0.0) {
1783
0
                *pyval = fay[i];
1784
0
                return 0;
1785
0
            }
1786
0
            im = i - 1;
1787
0
            dell = xval - fax[im];  /* >= 0 */
1788
0
            break;
1789
0
        }
1790
0
    }
1791
0
    fract = dell / (fax[i] - fax[im]);
1792
1793
0
    if (type == L_LINEAR_INTERP) {
1794
0
        *pyval = fay[i] + fract * (fay[i + 1] - fay[i]);
1795
0
        return 0;
1796
0
    }
1797
1798
        /* Quadratic interpolation */
1799
0
    if (im == 0) {
1800
0
        i1 = im;
1801
0
        i2 = im + 1;
1802
0
        i3 = im + 2;
1803
0
    } else {
1804
0
        i1 = im - 1;
1805
0
        i2 = im;
1806
0
        i3 = im + 1;
1807
0
    }
1808
0
    d1 = (fax[i1] - fax[i2]) * (fax[i1] - fax[i3]);
1809
0
    d2 = (fax[i2] - fax[i1]) * (fax[i2] - fax[i3]);
1810
0
    d3 = (fax[i3] - fax[i1]) * (fax[i3] - fax[i2]);
1811
0
    *pyval = fay[i1] * (xval - fax[i2]) * (xval - fax[i3]) / d1 +
1812
0
             fay[i2] * (xval - fax[i1]) * (xval - fax[i3]) / d2 +
1813
0
             fay[i3] * (xval - fax[i1]) * (xval - fax[i2]) / d3;
1814
0
    return 0;
1815
0
}
1816
1817
1818
/*!
1819
 * \brief   numaInterpolateEqxInterval()
1820
 *
1821
 * \param[in]    startx    xval corresponding to first element in nas
1822
 * \param[in]    deltax    x increment between array elements in nas
1823
 * \param[in]    nasy      numa of ordinate values, assumed equally spaced
1824
 * \param[in]    type      L_LINEAR_INTERP, L_QUADRATIC_INTERP
1825
 * \param[in]    x0        start value of interval
1826
 * \param[in]    x1        end value of interval
1827
 * \param[in]    npts      number of points to evaluate function in interval
1828
 * \param[out]   pnax      [optional] array of x values in interval
1829
 * \param[out]   pnay      array of y values in interval
1830
 * \return  0 if OK, 1 on error
1831
 *
1832
 * <pre>
1833
 * Notes:
1834
 *      (1) Considering nasy as a function of x, the x values
1835
 *          are equally spaced.
1836
 *      (2) This creates nay (and optionally nax) of interpolated
1837
 *          values over the specified interval (x0, x1).
1838
 *      (3) If the interval (x0, x1) lies partially outside the array
1839
 *          nasy (as interpreted by startx and deltax), it is an
1840
 *          error and returns 1.
1841
 *      (4) Note that deltax is the intrinsic x-increment for the input
1842
 *          array nasy, whereas delx is the intrinsic x-increment for the
1843
 *          output interpolated array nay.
1844
 * </pre>
1845
 */
1846
l_ok
1847
numaInterpolateEqxInterval(l_float32  startx,
1848
                           l_float32  deltax,
1849
                           NUMA      *nasy,
1850
                           l_int32    type,
1851
                           l_float32  x0,
1852
                           l_float32  x1,
1853
                           l_int32    npts,
1854
                           NUMA     **pnax,
1855
                           NUMA     **pnay)
1856
0
{
1857
0
l_int32     i, n;
1858
0
l_float32   x, yval, maxx, delx;
1859
0
NUMA       *nax = NULL, *nay;
1860
1861
0
    if (pnax) *pnax = NULL;
1862
0
    if (!pnay)
1863
0
        return ERROR_INT("&nay not defined", __func__, 1);
1864
0
    *pnay = NULL;
1865
0
    if (!nasy)
1866
0
        return ERROR_INT("nasy not defined", __func__, 1);
1867
0
    if ((n = numaGetCount(nasy)) < 2)
1868
0
        return ERROR_INT("n < 2", __func__, 1);
1869
0
    if (deltax <= 0.0)
1870
0
        return ERROR_INT("deltax not > 0", __func__, 1);
1871
0
    if (type != L_LINEAR_INTERP && type != L_QUADRATIC_INTERP)
1872
0
        return ERROR_INT("invalid interp type", __func__, 1);
1873
0
    if (type == L_QUADRATIC_INTERP && n == 2) {
1874
0
        type = L_LINEAR_INTERP;
1875
0
        L_WARNING("only 2 points; using linear interp\n", __func__);
1876
0
    }
1877
0
    maxx = startx + deltax * (n - 1);
1878
0
    if (x0 < startx || x1 > maxx || x1 <= x0)
1879
0
        return ERROR_INT("[x0 ... x1] is not valid", __func__, 1);
1880
0
    if (npts < 3)
1881
0
        return ERROR_INT("npts < 3", __func__, 1);
1882
0
    delx = (x1 - x0) / (l_float32)(npts - 1);  /* delx is for output nay */
1883
1884
0
    if ((nay = numaCreate(npts)) == NULL)
1885
0
        return ERROR_INT("nay not made", __func__, 1);
1886
0
    numaSetParameters(nay, x0, delx);
1887
0
    *pnay = nay;
1888
0
    if (pnax) {
1889
0
        nax = numaCreate(npts);
1890
0
        *pnax = nax;
1891
0
    }
1892
1893
0
    for (i = 0; i < npts; i++) {
1894
0
        x = x0 + i * delx;
1895
0
        if (pnax)
1896
0
            numaAddNumber(nax, x);
1897
0
        numaInterpolateEqxVal(startx, deltax, nasy, type, x, &yval);
1898
0
        numaAddNumber(nay, yval);
1899
0
    }
1900
1901
0
    return 0;
1902
0
}
1903
1904
1905
/*!
1906
 * \brief   numaInterpolateArbxInterval()
1907
 *
1908
 * \param[in]    nax     numa of abscissa values
1909
 * \param[in]    nay     numa of ordinate values, corresponding to nax
1910
 * \param[in]    type    L_LINEAR_INTERP, L_QUADRATIC_INTERP
1911
 * \param[in]    x0      start value of interval
1912
 * \param[in]    x1      end value of interval
1913
 * \param[in]    npts    number of points to evaluate function in interval
1914
 * \param[out]   pnadx   [optional] array of x values in interval
1915
 * \param[out]   pnady   array of y values in interval
1916
 * \return  0 if OK, 1 on error e.g., if x0 or x1 is outside range
1917
 *
1918
 * <pre>
1919
 * Notes:
1920
 *      (1) The values in nax must be sorted in increasing order.
1921
 *          If they are not sorted, we do it here, and complain.
1922
 *      (2) If the values in nax are equally spaced, you can use
1923
 *          numaInterpolateEqxInterval().
1924
 *      (3) Caller should check for valid return.
1925
 *      (4) We don't call numaInterpolateArbxVal() for each output
1926
 *          point, because that requires an O(n) search for
1927
 *          each point.  Instead, we do a single O(n) pass through
1928
 *          nax, saving the indices to be used for each output yval.
1929
 *      (5) Uses lagrangian interpolation.  See numaInterpolateEqxVal()
1930
 *          for formulas.
1931
 * </pre>
1932
 */
1933
l_ok
1934
numaInterpolateArbxInterval(NUMA       *nax,
1935
                            NUMA       *nay,
1936
                            l_int32     type,
1937
                            l_float32   x0,
1938
                            l_float32   x1,
1939
                            l_int32     npts,
1940
                            NUMA      **pnadx,
1941
                            NUMA      **pnady)
1942
0
{
1943
0
l_int32     i, im, j, nx, ny, i1, i2, i3, sorted;
1944
0
l_int32    *index;
1945
0
l_float32   del, xval, yval, excess, fract, minx, maxx, d1, d2, d3;
1946
0
l_float32  *fax, *fay;
1947
0
NUMA       *nasx, *nasy, *nadx = NULL, *nady;
1948
1949
0
    if (pnadx) *pnadx = NULL;
1950
0
    if (!pnady)
1951
0
        return ERROR_INT("&nady not defined", __func__, 1);
1952
0
    *pnady = NULL;
1953
0
    if (!nay)
1954
0
        return ERROR_INT("nay not defined", __func__, 1);
1955
0
    if (!nax)
1956
0
        return ERROR_INT("nax not defined", __func__, 1);
1957
0
    if (type != L_LINEAR_INTERP && type != L_QUADRATIC_INTERP)
1958
0
        return ERROR_INT("invalid interp type", __func__, 1);
1959
0
    if (x0 > x1)
1960
0
        return ERROR_INT("x0 > x1", __func__, 1);
1961
0
    ny = numaGetCount(nay);
1962
0
    nx = numaGetCount(nax);
1963
0
    if (nx != ny)
1964
0
        return ERROR_INT("nax and nay not same size arrays", __func__, 1);
1965
0
    if (ny < 2)
1966
0
        return ERROR_INT("not enough points", __func__, 1);
1967
0
    if (type == L_QUADRATIC_INTERP && ny == 2) {
1968
0
        type = L_LINEAR_INTERP;
1969
0
        L_WARNING("only 2 points; using linear interp\n", __func__);
1970
0
    }
1971
0
    numaGetMin(nax, &minx, NULL);
1972
0
    numaGetMax(nax, &maxx, NULL);
1973
0
    if (x0 < minx || x1 > maxx)
1974
0
        return ERROR_INT("xval is out of bounds", __func__, 1);
1975
1976
        /* Make sure that nax is sorted in increasing order */
1977
0
    numaIsSorted(nax, L_SORT_INCREASING, &sorted);
1978
0
    if (!sorted) {
1979
0
        L_WARNING("we are sorting nax in increasing order\n", __func__);
1980
0
        numaSortPair(nax, nay, L_SORT_INCREASING, &nasx, &nasy);
1981
0
    } else {
1982
0
        nasx = numaClone(nax);
1983
0
        nasy = numaClone(nay);
1984
0
    }
1985
1986
0
    fax = numaGetFArray(nasx, L_NOCOPY);
1987
0
    fay = numaGetFArray(nasy, L_NOCOPY);
1988
1989
        /* Get array of indices into fax for interpolated locations */
1990
0
    if ((index = (l_int32 *)LEPT_CALLOC(npts, sizeof(l_int32))) == NULL) {
1991
0
        numaDestroy(&nasx);
1992
0
        numaDestroy(&nasy);
1993
0
        return ERROR_INT("ind not made", __func__, 1);
1994
0
    }
1995
0
    del = (x1 - x0) / (npts - 1.0f);
1996
0
    for (i = 0, j = 0; j < nx && i < npts; i++) {
1997
0
        xval = x0 + i * del;
1998
0
        while (j < nx - 1 && xval > fax[j])
1999
0
            j++;
2000
0
        if (xval == fax[j])
2001
0
            index[i] = L_MIN(j, nx - 1);
2002
0
        else    /* the index of fax[] is just below xval */
2003
0
            index[i] = L_MAX(j - 1, 0);
2004
0
    }
2005
2006
        /* For each point to be interpolated, get the y value */
2007
0
    nady = numaCreate(npts);
2008
0
    *pnady = nady;
2009
0
    if (pnadx) {
2010
0
        nadx = numaCreate(npts);
2011
0
        *pnadx = nadx;
2012
0
    }
2013
0
    for (i = 0; i < npts; i++) {
2014
0
        xval = x0 + i * del;
2015
0
        if (pnadx)
2016
0
            numaAddNumber(nadx, xval);
2017
0
        im = index[i];
2018
0
        excess = xval - fax[im];
2019
0
        if (excess == 0.0) {
2020
0
            numaAddNumber(nady, fay[im]);
2021
0
            continue;
2022
0
        }
2023
0
        fract = excess / (fax[im + 1] - fax[im]);
2024
2025
0
        if (type == L_LINEAR_INTERP) {
2026
0
            yval = fay[im] + fract * (fay[im + 1] - fay[im]);
2027
0
            numaAddNumber(nady, yval);
2028
0
            continue;
2029
0
        }
2030
2031
            /* Quadratic interpolation */
2032
0
        if (im == 0) {
2033
0
            i1 = im;
2034
0
            i2 = im + 1;
2035
0
            i3 = im + 2;
2036
0
        } else {
2037
0
            i1 = im - 1;
2038
0
            i2 = im;
2039
0
            i3 = im + 1;
2040
0
        }
2041
0
        d1 = (fax[i1] - fax[i2]) * (fax[i1] - fax[i3]);
2042
0
        d2 = (fax[i2] - fax[i1]) * (fax[i2] - fax[i3]);
2043
0
        d3 = (fax[i3] - fax[i1]) * (fax[i3] - fax[i2]);
2044
0
        yval = fay[i1] * (xval - fax[i2]) * (xval - fax[i3]) / d1 +
2045
0
               fay[i2] * (xval - fax[i1]) * (xval - fax[i3]) / d2 +
2046
0
               fay[i3] * (xval - fax[i1]) * (xval - fax[i2]) / d3;
2047
0
        numaAddNumber(nady, yval);
2048
0
    }
2049
2050
0
    LEPT_FREE(index);
2051
0
    numaDestroy(&nasx);
2052
0
    numaDestroy(&nasy);
2053
0
    return 0;
2054
0
}
2055
2056
2057
/*----------------------------------------------------------------------*
2058
 *                     Functions requiring interpolation                *
2059
 *----------------------------------------------------------------------*/
2060
/*!
2061
 * \brief   numaFitMax()
2062
 *
2063
 * \param[in]    na       numa of ordinate values, to fit a max to
2064
 * \param[out]   pmaxval  max value
2065
 * \param[in]    naloc    [optional] associated numa of abscissa values
2066
 * \param[out]   pmaxloc  abscissa value that gives max value in na;
2067
 *                        if naloc == null, this is given as an interpolated
2068
 *                        index value
2069
 * \return  0 if OK; 1 on error
2070
 *
2071
 * <pre>
2072
 * Notes:
2073
 *        If %naloc is given, there is no requirement that the
2074
 *        data points are evenly spaced.  Lagrangian interpolation
2075
 *        handles that.  The only requirement is that the
2076
 *        data points are ordered so that the values in naloc
2077
 *        are either increasing or decreasing.  We test to make
2078
 *        sure that the sizes of na and naloc are equal, and it
2079
 *        is assumed that the correspondences %na[i] as a function
2080
 *        of %naloc[i] are properly arranged for all i.
2081
 *
2082
 *  The formula for Lagrangian interpolation through 3 data pts is:
2083
 *       y(x) = y1(x-x2)(x-x3)/((x1-x2)(x1-x3)) +
2084
 *              y2(x-x1)(x-x3)/((x2-x1)(x2-x3)) +
2085
 *              y3(x-x1)(x-x2)/((x3-x1)(x3-x2))
2086
 *
2087
 *  Then the derivative, using the constants (c1,c2,c3) defined below,
2088
 *  is set to 0:
2089
 *       y'(x) = 2x(c1+c2+c3) - c1(x2+x3) - c2(x1+x3) - c3(x1+x2) = 0
2090
 * </pre>
2091
 */
2092
l_ok
2093
numaFitMax(NUMA       *na,
2094
           l_float32  *pmaxval,
2095
           NUMA       *naloc,
2096
           l_float32  *pmaxloc)
2097
0
{
2098
0
l_float32  val;
2099
0
l_float32  smaxval;  /* start value of maximum sample, before interpolating */
2100
0
l_int32    n, imaxloc;
2101
0
l_float32  x1, x2, x3, y1, y2, y3, c1, c2, c3, a, b, xmax, ymax;
2102
2103
0
    if (pmaxval) *pmaxval = 0.0;
2104
0
    if (pmaxloc) *pmaxloc = 0.0;
2105
0
    if (!na)
2106
0
        return ERROR_INT("na not defined", __func__, 1);
2107
0
    if ((n = numaGetCount(na)) == 0)
2108
0
        return ERROR_INT("na is empty", __func__, 1);
2109
0
    if (!pmaxval)
2110
0
        return ERROR_INT("&maxval not defined", __func__, 1);
2111
0
    if (!pmaxloc)
2112
0
        return ERROR_INT("&maxloc not defined", __func__, 1);
2113
0
    if (naloc) {
2114
0
        if (n != numaGetCount(naloc))
2115
0
            return ERROR_INT("na and naloc of unequal size", __func__, 1);
2116
0
    }
2117
2118
0
    numaGetMax(na, &smaxval, &imaxloc);
2119
2120
        /* Simple case: max is at end point */
2121
0
    if (imaxloc == 0 || imaxloc == n - 1) {
2122
0
        *pmaxval = smaxval;
2123
0
        if (naloc) {
2124
0
            numaGetFValue(naloc, imaxloc, &val);
2125
0
            *pmaxloc = val;
2126
0
        } else {
2127
0
            *pmaxloc = imaxloc;
2128
0
        }
2129
0
        return 0;
2130
0
    }
2131
2132
        /* Interior point; use quadratic interpolation */
2133
0
    y2 = smaxval;
2134
0
    numaGetFValue(na, imaxloc - 1, &val);
2135
0
    y1 = val;
2136
0
    numaGetFValue(na, imaxloc + 1, &val);
2137
0
    y3 = val;
2138
0
    if (naloc) {
2139
0
        numaGetFValue(naloc, imaxloc - 1, &val);
2140
0
        x1 = val;
2141
0
        numaGetFValue(naloc, imaxloc, &val);
2142
0
        x2 = val;
2143
0
        numaGetFValue(naloc, imaxloc + 1, &val);
2144
0
        x3 = val;
2145
0
    } else {
2146
0
        x1 = imaxloc - 1;
2147
0
        x2 = imaxloc;
2148
0
        x3 = imaxloc + 1;
2149
0
    }
2150
2151
        /* Can't interpolate; just use the max val in na
2152
         * and the corresponding one in naloc */
2153
0
    if (x1 == x2 || x1 == x3 || x2 == x3) {
2154
0
        *pmaxval = y2;
2155
0
        *pmaxloc = x2;
2156
0
        return 0;
2157
0
    }
2158
2159
        /* Use lagrangian interpolation; set dy/dx = 0 */
2160
0
    c1 = y1 / ((x1 - x2) * (x1 - x3));
2161
0
    c2 = y2 / ((x2 - x1) * (x2 - x3));
2162
0
    c3 = y3 / ((x3 - x1) * (x3 - x2));
2163
0
    a = c1 + c2 + c3;
2164
0
    b = c1 * (x2 + x3) + c2 * (x1 + x3) + c3 * (x1 + x2);
2165
0
    xmax = b / (2 * a);
2166
0
    ymax = c1 * (xmax - x2) * (xmax - x3) +
2167
0
           c2 * (xmax - x1) * (xmax - x3) +
2168
0
           c3 * (xmax - x1) * (xmax - x2);
2169
0
    *pmaxval = ymax;
2170
0
    *pmaxloc = xmax;
2171
2172
0
    return 0;
2173
0
}
2174
2175
2176
/*!
2177
 * \brief   numaDifferentiateInterval()
2178
 *
2179
 * \param[in]    nax     numa of abscissa values
2180
 * \param[in]    nay     numa of ordinate values, corresponding to nax
2181
 * \param[in]    x0      start value of interval
2182
 * \param[in]    x1      end value of interval
2183
 * \param[in]    npts    number of points to evaluate function in interval
2184
 * \param[out]   pnadx   [optional] array of x values in interval
2185
 * \param[out]   pnady   array of derivatives in interval
2186
 * \return  0 if OK, 1 on error e.g., if x0 or x1 is outside range
2187
 *
2188
 * <pre>
2189
 * Notes:
2190
 *      (1) The values in nax must be sorted in increasing order.
2191
 *          If they are not sorted, it is done in the interpolation
2192
 *          step, and a warning is issued.
2193
 *      (2) Caller should check for valid return.
2194
 * </pre>
2195
 */
2196
l_ok
2197
numaDifferentiateInterval(NUMA       *nax,
2198
                          NUMA       *nay,
2199
                          l_float32   x0,
2200
                          l_float32   x1,
2201
                          l_int32     npts,
2202
                          NUMA      **pnadx,
2203
                          NUMA      **pnady)
2204
0
{
2205
0
l_int32     i, nx, ny;
2206
0
l_float32   minx, maxx, der, invdel;
2207
0
l_float32  *fay;
2208
0
NUMA       *nady, *naiy;
2209
2210
0
    if (pnadx) *pnadx = NULL;
2211
0
    if (!pnady)
2212
0
        return ERROR_INT("&nady not defined", __func__, 1);
2213
0
    *pnady = NULL;
2214
0
    if (!nay)
2215
0
        return ERROR_INT("nay not defined", __func__, 1);
2216
0
    if (!nax)
2217
0
        return ERROR_INT("nax not defined", __func__, 1);
2218
0
    if (x0 > x1)
2219
0
        return ERROR_INT("x0 > x1", __func__, 1);
2220
0
    ny = numaGetCount(nay);
2221
0
    nx = numaGetCount(nax);
2222
0
    if (nx != ny)
2223
0
        return ERROR_INT("nax and nay not same size arrays", __func__, 1);
2224
0
    if (ny < 2)
2225
0
        return ERROR_INT("not enough points", __func__, 1);
2226
0
    numaGetMin(nax, &minx, NULL);
2227
0
    numaGetMax(nax, &maxx, NULL);
2228
0
    if (x0 < minx || x1 > maxx)
2229
0
        return ERROR_INT("xval is out of bounds", __func__, 1);
2230
0
    if (npts < 2)
2231
0
        return ERROR_INT("npts < 2", __func__, 1);
2232
2233
        /* Generate interpolated array over specified interval */
2234
0
    if (numaInterpolateArbxInterval(nax, nay, L_LINEAR_INTERP, x0, x1,
2235
0
                                    npts, pnadx, &naiy))
2236
0
        return ERROR_INT("interpolation failed", __func__, 1);
2237
2238
0
    nady = numaCreate(npts);
2239
0
    *pnady = nady;
2240
0
    invdel = 0.5f * ((l_float32)npts - 1.0f) / (x1 - x0);
2241
0
    fay = numaGetFArray(naiy, L_NOCOPY);
2242
2243
        /* Compute and save derivatives */
2244
0
    der = 0.5f * invdel * (fay[1] - fay[0]);
2245
0
    numaAddNumber(nady, der);
2246
0
    for (i = 1; i < npts - 1; i++)  {
2247
0
        der = invdel * (fay[i + 1] - fay[i - 1]);
2248
0
        numaAddNumber(nady, der);
2249
0
    }
2250
0
    der = 0.5f * invdel * (fay[npts - 1] - fay[npts - 2]);
2251
0
    numaAddNumber(nady, der);
2252
2253
0
    numaDestroy(&naiy);
2254
0
    return 0;
2255
0
}
2256
2257
2258
/*!
2259
 * \brief   numaIntegrateInterval()
2260
 *
2261
 * \param[in]    nax     numa of abscissa values
2262
 * \param[in]    nay     numa of ordinate values, corresponding to nax
2263
 * \param[in]    x0      start value of interval
2264
 * \param[in]    x1      end value of interval
2265
 * \param[in]    npts    number of points to evaluate function in interval
2266
 * \param[out]   psum    integral of function over interval
2267
 * \return  0 if OK, 1 on error e.g., if x0 or x1 is outside range
2268
 *
2269
 * <pre>
2270
 * Notes:
2271
 *      (1) The values in nax must be sorted in increasing order.
2272
 *          If they are not sorted, it is done in the interpolation
2273
 *          step, and a warning is issued.
2274
 *      (2) Caller should check for valid return.
2275
 * </pre>
2276
 */
2277
l_ok
2278
numaIntegrateInterval(NUMA       *nax,
2279
                      NUMA       *nay,
2280
                      l_float32   x0,
2281
                      l_float32   x1,
2282
                      l_int32     npts,
2283
                      l_float32  *psum)
2284
0
{
2285
0
l_int32     i, nx, ny;
2286
0
l_float32   minx, maxx, sum, del;
2287
0
l_float32  *fay;
2288
0
NUMA       *naiy;
2289
2290
0
    if (!psum)
2291
0
        return ERROR_INT("&sum not defined", __func__, 1);
2292
0
    *psum = 0.0;
2293
0
    if (!nay)
2294
0
        return ERROR_INT("nay not defined", __func__, 1);
2295
0
    if (!nax)
2296
0
        return ERROR_INT("nax not defined", __func__, 1);
2297
0
    if (x0 > x1)
2298
0
        return ERROR_INT("x0 > x1", __func__, 1);
2299
0
    if (npts < 2)
2300
0
        return ERROR_INT("npts < 2", __func__, 1);
2301
0
    ny = numaGetCount(nay);
2302
0
    nx = numaGetCount(nax);
2303
0
    if (nx != ny)
2304
0
        return ERROR_INT("nax and nay not same size arrays", __func__, 1);
2305
0
    if (ny < 2)
2306
0
        return ERROR_INT("not enough points", __func__, 1);
2307
0
    numaGetMin(nax, &minx, NULL);
2308
0
    numaGetMax(nax, &maxx, NULL);
2309
0
    if (x0 < minx || x1 > maxx)
2310
0
        return ERROR_INT("xval is out of bounds", __func__, 1);
2311
2312
        /* Generate interpolated array over specified interval */
2313
0
    if (numaInterpolateArbxInterval(nax, nay, L_LINEAR_INTERP, x0, x1,
2314
0
                                    npts, NULL, &naiy))
2315
0
        return ERROR_INT("interpolation failed", __func__, 1);
2316
2317
0
    del = (x1 - x0) / ((l_float32)npts - 1.0f);
2318
0
    fay = numaGetFArray(naiy, L_NOCOPY);
2319
2320
        /* Compute integral (simple trapezoid) */
2321
0
    sum = 0.5f * (fay[0] + fay[npts - 1]);
2322
0
    for (i = 1; i < npts - 1; i++)
2323
0
        sum += fay[i];
2324
0
    *psum = del * sum;
2325
2326
0
    numaDestroy(&naiy);
2327
0
    return 0;
2328
0
}
2329
2330
2331
/*----------------------------------------------------------------------*
2332
 *                                Sorting                               *
2333
 *----------------------------------------------------------------------*/
2334
/*!
2335
 * \brief   numaSortGeneral()
2336
 *
2337
 * \param[in]    na          source numa
2338
 * \param[out]   pnasort     [optional] sorted numa
2339
 * \param[out]   pnaindex    [optional] index of elements in na associated
2340
 *                           with each element of nasort
2341
 * \param[out]   pnainvert   [optional] index of elements in nasort associated
2342
 *                           with each element of na
2343
 * \param[in]    sortorder   L_SORT_INCREASING or L_SORT_DECREASING
2344
 * \param[in]    sorttype    L_SHELL_SORT or L_BIN_SORT
2345
 * \return  0 if OK, 1 on error
2346
 *
2347
 * <pre>
2348
 * Notes:
2349
 *      (1) Sorting can be confusing.  Here's an array of five values with
2350
 *          the results shown for the 3 output arrays.
2351
 *
2352
 *          na      nasort   naindex   nainvert
2353
 *          -----------------------------------
2354
 *          3         9         2         3
2355
 *          4         6         3         2
2356
 *          9         4         1         0
2357
 *          6         3         0         1
2358
 *          1         1         4         4
2359
 *
2360
 *          Note that naindex is a LUT into na for the sorted array values,
2361
 *          and nainvert directly gives the sorted index values for the
2362
 *          input array.  It is useful to view naindex is as a map:
2363
 *                 0  -->  2
2364
 *                 1  -->  3
2365
 *                 2  -->  1
2366
 *                 3  -->  0
2367
 *                 4  -->  4
2368
 *          and nainvert, the inverse of this map:
2369
 *                 0  -->  3
2370
 *                 1  -->  2
2371
 *                 2  -->  0
2372
 *                 3  -->  1
2373
 *                 4  -->  4
2374
 *
2375
 *          We can write these relations symbolically as:
2376
 *              nasort[i] = na[naindex[i]]
2377
 *              na[i] = nasort[nainvert[i]]
2378
 * </pre>
2379
 */
2380
l_ok
2381
numaSortGeneral(NUMA    *na,
2382
                NUMA   **pnasort,
2383
                NUMA   **pnaindex,
2384
                NUMA   **pnainvert,
2385
                l_int32  sortorder,
2386
                l_int32  sorttype)
2387
0
{
2388
0
l_int32    isize;
2389
0
l_float32  size;
2390
0
NUMA      *naindex = NULL;
2391
2392
0
    if (pnasort) *pnasort = NULL;
2393
0
    if (pnaindex) *pnaindex = NULL;
2394
0
    if (pnainvert) *pnainvert = NULL;
2395
0
    if (!na)
2396
0
        return ERROR_INT("na not defined", __func__, 1);
2397
0
    if (sortorder != L_SORT_INCREASING && sortorder != L_SORT_DECREASING)
2398
0
        return ERROR_INT("invalid sort order", __func__, 1);
2399
0
    if (sorttype != L_SHELL_SORT && sorttype != L_BIN_SORT)
2400
0
        return ERROR_INT("invalid sort type", __func__, 1);
2401
0
    if (!pnasort && !pnaindex && !pnainvert)
2402
0
        return ERROR_INT("nothing to do", __func__, 1);
2403
2404
0
    if (sorttype == L_BIN_SORT) {
2405
0
        numaGetMax(na, &size, NULL);
2406
0
        isize = (l_int32)size;
2407
0
        if (isize > MaxInitPtraSize - 1) {
2408
0
            L_WARNING("array too large; using shell sort\n", __func__);
2409
0
            sorttype = L_SHELL_SORT;
2410
0
        } else {
2411
0
            naindex = numaGetBinSortIndex(na, sortorder);
2412
0
        }
2413
0
    }
2414
2415
0
    if (sorttype == L_SHELL_SORT)
2416
0
        naindex = numaGetSortIndex(na, sortorder);
2417
2418
0
    if (pnasort)
2419
0
        *pnasort = numaSortByIndex(na, naindex);
2420
0
    if (pnainvert)
2421
0
        *pnainvert = numaInvertMap(naindex);
2422
0
    if (pnaindex)
2423
0
        *pnaindex = naindex;
2424
0
    else
2425
0
        numaDestroy(&naindex);
2426
0
    return 0;
2427
0
}
2428
2429
2430
/*!
2431
 * \brief   numaSortAutoSelect()
2432
 *
2433
 * \param[in]    nas
2434
 * \param[in]    sortorder   L_SORT_INCREASING or L_SORT_DECREASING
2435
 * \return  naout output sorted numa, or NULL on error
2436
 *
2437
 * <pre>
2438
 * Notes:
2439
 *      (1) This does either a shell sort or a bin sort, depending on
2440
 *          the number of elements in nas and the dynamic range.
2441
 * </pre>
2442
 */
2443
NUMA *
2444
numaSortAutoSelect(NUMA    *nas,
2445
                   l_int32  sortorder)
2446
0
{
2447
0
l_int32  type;
2448
2449
0
    if (!nas)
2450
0
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
2451
0
    if (numaGetCount(nas) == 0) {
2452
0
        L_WARNING("nas is empty; returning copy\n", __func__);
2453
0
        return numaCopy(nas);
2454
0
    }
2455
0
    if (sortorder != L_SORT_INCREASING && sortorder != L_SORT_DECREASING)
2456
0
        return (NUMA *)ERROR_PTR("invalid sort order", __func__, NULL);
2457
2458
0
    type = numaChooseSortType(nas);
2459
0
    if (type != L_SHELL_SORT && type != L_BIN_SORT)
2460
0
        return (NUMA *)ERROR_PTR("invalid sort type", __func__, NULL);
2461
2462
0
    if (type == L_BIN_SORT)
2463
0
        return numaBinSort(nas, sortorder);
2464
0
    else  /* shell sort */
2465
0
        return numaSort(NULL, nas, sortorder);
2466
0
}
2467
2468
2469
/*!
2470
 * \brief   numaSortIndexAutoSelect()
2471
 *
2472
 * \param[in]    nas
2473
 * \param[in]    sortorder     L_SORT_INCREASING or L_SORT_DECREASING
2474
 * \return  nad indices of nas, sorted by value in nas, or NULL on error
2475
 *
2476
 * <pre>
2477
 * Notes:
2478
 *      (1) This does either a shell sort or a bin sort, depending on
2479
 *          the number of elements in nas and the dynamic range.
2480
 * </pre>
2481
 */
2482
NUMA *
2483
numaSortIndexAutoSelect(NUMA    *nas,
2484
                        l_int32  sortorder)
2485
0
{
2486
0
l_int32  type;
2487
2488
0
    if (!nas)
2489
0
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
2490
0
    if (numaGetCount(nas) == 0) {
2491
0
        L_WARNING("nas is empty; returning copy\n", __func__);
2492
0
        return numaCopy(nas);
2493
0
    }
2494
0
    if (sortorder != L_SORT_INCREASING && sortorder != L_SORT_DECREASING)
2495
0
        return (NUMA *)ERROR_PTR("invalid sort order", __func__, NULL);
2496
0
    type = numaChooseSortType(nas);
2497
0
    if (type != L_SHELL_SORT && type != L_BIN_SORT)
2498
0
        return (NUMA *)ERROR_PTR("invalid sort type", __func__, NULL);
2499
2500
0
    if (type == L_BIN_SORT)
2501
0
        return numaGetBinSortIndex(nas, sortorder);
2502
0
    else  /* shell sort */
2503
0
        return numaGetSortIndex(nas, sortorder);
2504
0
}
2505
2506
2507
/*!
2508
 * \brief   numaChooseSortType()
2509
 *
2510
 * \param[in]    nas     to be sorted
2511
 * \return  sorttype  L_SHELL_SORT or L_BIN_SORT, or UNDEF on error.
2512
 *
2513
 * <pre>
2514
 * Notes:
2515
 *      (1) This selects either a shell sort or a bin sort, depending on
2516
 *          the number of elements in nas and the dynamic range.
2517
 *      (2) If there are negative values in nas, it selects shell sort.
2518
 * </pre>
2519
 */
2520
l_int32
2521
numaChooseSortType(NUMA  *nas)
2522
0
{
2523
0
l_int32    n;
2524
0
l_float32  minval, maxval;
2525
2526
0
    if (!nas)
2527
0
        return ERROR_INT("nas not defined", __func__, UNDEF);
2528
2529
        /* If small histogram or negative values; use shell sort */
2530
0
    numaGetMin(nas, &minval, NULL);
2531
0
    n = numaGetCount(nas);
2532
0
    if (minval < 0.0 || n < 200)
2533
0
        return L_SHELL_SORT;
2534
2535
        /* If large maxval, use shell sort */
2536
0
    numaGetMax(nas, &maxval, NULL);
2537
0
    if (maxval > MaxInitPtraSize - 1)
2538
0
        return L_SHELL_SORT;
2539
2540
        /* Otherwise, need to compare nlog(n) with maxval.
2541
         * The factor of 0.003 was determined by comparing times for
2542
         * different histogram sizes and maxval.  It is very small
2543
         * because binsort is fast and shell sort gets slow for large n. */
2544
0
    if (n * log((l_float32)n) < 0.003 * maxval)
2545
0
        return L_SHELL_SORT;
2546
0
    else
2547
0
        return L_BIN_SORT;
2548
0
}
2549
2550
2551
/*!
2552
 * \brief   numaSort()
2553
 *
2554
 * \param[in]    naout       output numa; can be NULL or equal to nain
2555
 * \param[in]    nain        input numa
2556
 * \param[in]    sortorder   L_SORT_INCREASING or L_SORT_DECREASING
2557
 * \return  naout   output sorted numa, or NULL on error
2558
 *
2559
 * <pre>
2560
 * Notes:
2561
 *      (1) Set naout = nain for in-place; otherwise, set naout = NULL.
2562
 *      (2) Source: Shell sort, modified from K&R, 2nd edition, p.62.
2563
 *          Slow but simple O(n logn) sort.
2564
 * </pre>
2565
 */
2566
NUMA *
2567
numaSort(NUMA    *naout,
2568
         NUMA    *nain,
2569
         l_int32  sortorder)
2570
0
{
2571
0
l_int32     i, n, gap, j;
2572
0
l_float32   tmp;
2573
0
l_float32  *array;
2574
2575
0
    if (!nain)
2576
0
        return (NUMA *)ERROR_PTR("nain not defined", __func__, NULL);
2577
0
    if (sortorder != L_SORT_INCREASING && sortorder != L_SORT_DECREASING)
2578
0
        return (NUMA *)ERROR_PTR("invalid sort order", __func__, NULL);
2579
2580
        /* Make naout if necessary; otherwise do in-place */
2581
0
    if (!naout)
2582
0
        naout = numaCopy(nain);
2583
0
    else if (nain != naout)
2584
0
        return (NUMA *)ERROR_PTR("invalid: not in-place", __func__, NULL);
2585
0
    if ((n = numaGetCount(naout)) == 0) {
2586
0
        L_WARNING("naout is empty\n", __func__);
2587
0
        return naout;
2588
0
    }
2589
0
    array = naout->array;  /* operate directly on the array */
2590
0
    n = numaGetCount(naout);
2591
2592
        /* Shell sort */
2593
0
    for (gap = n/2; gap > 0; gap = gap / 2) {
2594
0
        for (i = gap; i < n; i++) {
2595
0
            for (j = i - gap; j >= 0; j -= gap) {
2596
0
                if ((sortorder == L_SORT_INCREASING &&
2597
0
                     array[j] > array[j + gap]) ||
2598
0
                    (sortorder == L_SORT_DECREASING &&
2599
0
                     array[j] < array[j + gap]))
2600
0
                {
2601
0
                    tmp = array[j];
2602
0
                    array[j] = array[j + gap];
2603
0
                    array[j + gap] = tmp;
2604
0
                }
2605
0
            }
2606
0
        }
2607
0
    }
2608
2609
0
    return naout;
2610
0
}
2611
2612
2613
/*!
2614
 * \brief   numaBinSort()
2615
 *
2616
 * \param[in]    nas         of non-negative integers with a max that can
2617
 *                           not exceed (MaxInitPtraSize - 1)
2618
 * \param[in]    sortorder   L_SORT_INCREASING or L_SORT_DECREASING
2619
 * \return  na   sorted, or NULL on error
2620
 *
2621
 * <pre>
2622
 * Notes:
2623
 *      (1) Because this uses a bin sort with buckets of size 1, it
2624
 *          is not appropriate for sorting either small arrays or
2625
 *          arrays containing very large integer values.  For such
2626
 *          arrays, use a standard general sort function like
2627
 *          numaSort().
2628
 *      (2) You can use numaSortAutoSelect() to decide which sorting
2629
 *          method to use.
2630
 * </pre>
2631
 */
2632
NUMA *
2633
numaBinSort(NUMA    *nas,
2634
            l_int32  sortorder)
2635
0
{
2636
0
NUMA  *nat, *nad;
2637
2638
0
    if (!nas)
2639
0
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
2640
0
    if (numaGetCount(nas) == 0) {
2641
0
        L_WARNING("nas is empty; returning copy\n", __func__);
2642
0
        return numaCopy(nas);
2643
0
    }
2644
0
    if (sortorder != L_SORT_INCREASING && sortorder != L_SORT_DECREASING)
2645
0
        return (NUMA *)ERROR_PTR("invalid sort order", __func__, NULL);
2646
2647
0
    if ((nat = numaGetBinSortIndex(nas, sortorder)) == NULL)
2648
0
        return (NUMA *)ERROR_PTR("bin sort failed", __func__, NULL);
2649
0
    nad = numaSortByIndex(nas, nat);
2650
0
    numaDestroy(&nat);
2651
0
    return nad;
2652
0
}
2653
2654
2655
/*!
2656
 * \brief   numaGetSortIndex()
2657
 *
2658
 * \param[in]    na          source numa
2659
 * \param[in]    sortorder   L_SORT_INCREASING or L_SORT_DECREASING
2660
 * \return  na  giving an array of indices that would sort
2661
 *              the input array, or NULL on error
2662
 */
2663
NUMA *
2664
numaGetSortIndex(NUMA    *na,
2665
                 l_int32  sortorder)
2666
0
{
2667
0
l_int32     i, n, gap, j;
2668
0
l_float32   tmp;
2669
0
l_float32  *array;   /* copy of input array */
2670
0
l_float32  *iarray;  /* array of indices */
2671
0
NUMA       *naisort;
2672
2673
0
    if (!na)
2674
0
        return (NUMA *)ERROR_PTR("na not defined", __func__, NULL);
2675
0
    if (numaGetCount(na) == 0) {
2676
0
        L_WARNING("na is empty\n", __func__);
2677
0
        return numaCreate(1);
2678
0
    }
2679
0
    if (sortorder != L_SORT_INCREASING && sortorder != L_SORT_DECREASING)
2680
0
        return (NUMA *)ERROR_PTR("invalid sortorder", __func__, NULL);
2681
2682
0
    n = numaGetCount(na);
2683
0
    if ((array = numaGetFArray(na, L_COPY)) == NULL)
2684
0
        return (NUMA *)ERROR_PTR("array not made", __func__, NULL);
2685
0
    if ((iarray = (l_float32 *)LEPT_CALLOC(n, sizeof(l_float32))) == NULL) {
2686
0
        LEPT_FREE(array);
2687
0
        return (NUMA *)ERROR_PTR("iarray not made", __func__, NULL);
2688
0
    }
2689
0
    for (i = 0; i < n; i++)
2690
0
        iarray[i] = i;
2691
2692
        /* Shell sort */
2693
0
    for (gap = n/2; gap > 0; gap = gap / 2) {
2694
0
        for (i = gap; i < n; i++) {
2695
0
            for (j = i - gap; j >= 0; j -= gap) {
2696
0
                if ((sortorder == L_SORT_INCREASING &&
2697
0
                     array[j] > array[j + gap]) ||
2698
0
                    (sortorder == L_SORT_DECREASING &&
2699
0
                     array[j] < array[j + gap]))
2700
0
                {
2701
0
                    tmp = array[j];
2702
0
                    array[j] = array[j + gap];
2703
0
                    array[j + gap] = tmp;
2704
0
                    tmp = iarray[j];
2705
0
                    iarray[j] = iarray[j + gap];
2706
0
                    iarray[j + gap] = tmp;
2707
0
                }
2708
0
            }
2709
0
        }
2710
0
    }
2711
2712
0
    naisort = numaCreate(n);
2713
0
    for (i = 0; i < n; i++)
2714
0
        numaAddNumber(naisort, iarray[i]);
2715
2716
0
    LEPT_FREE(array);
2717
0
    LEPT_FREE(iarray);
2718
0
    return naisort;
2719
0
}
2720
2721
2722
/*!
2723
 * \brief   numaGetBinSortIndex()
2724
 *
2725
 * \param[in]    nas         of non-negative integers with a max that can
2726
 *                           not exceed (MaxInitPtraSize - 1)
2727
 * \param[in]    sortorder   L_SORT_INCREASING or L_SORT_DECREASING
2728
 * \return  na  sorted, or NULL on error
2729
 *
2730
 * <pre>
2731
 * Notes:
2732
 *      (1) This creates an array (or lookup table) that contains
2733
 *          the sorted position of the elements in the input Numa.
2734
 *      (2) Because it uses a bin sort with buckets of size 1, it
2735
 *          is not appropriate for sorting either small arrays or
2736
 *          arrays containing very large integer values.  For such
2737
 *          arrays, use a standard general sort function like
2738
 *          numaGetSortIndex().
2739
 *      (3) You can use numaSortIndexAutoSelect() to decide which
2740
 *          sorting method to use.
2741
 * </pre>
2742
 */
2743
NUMA *
2744
numaGetBinSortIndex(NUMA    *nas,
2745
                    l_int32  sortorder)
2746
0
{
2747
0
l_int32    i, n, isize, ival, imax;
2748
0
l_float32  minsize, size;
2749
0
NUMA      *na, *nai, *nad;
2750
0
L_PTRA    *paindex;
2751
2752
0
    if (!nas)
2753
0
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
2754
0
    if (numaGetCount(nas) == 0) {
2755
0
        L_WARNING("nas is empty\n", __func__);
2756
0
        return numaCreate(1);
2757
0
    }
2758
0
    if (sortorder != L_SORT_INCREASING && sortorder != L_SORT_DECREASING)
2759
0
        return (NUMA *)ERROR_PTR("invalid sort order", __func__, NULL);
2760
0
    numaGetMin(nas, &minsize, NULL);
2761
0
    if (minsize < 0)
2762
0
        return (NUMA *)ERROR_PTR("nas has negative numbers", __func__, NULL);
2763
0
    numaGetMax(nas, &size, NULL);
2764
0
    isize = (l_int32)size;
2765
0
    if (isize > MaxInitPtraSize - 1) {
2766
0
        L_ERROR("array too large: %d elements > max size = %d\n",
2767
0
                __func__, isize, MaxInitPtraSize - 1);
2768
0
        return NULL;
2769
0
    }
2770
2771
        /* Set up a ptra holding numa at indices for which there
2772
         * are values in nas.  Suppose nas has the value 230 at index
2773
         * 7355.  A numa holding the index 7355 is created and stored
2774
         * at the ptra index 230.  If there is another value of 230
2775
         * in nas, its index is added to the same numa (at index 230
2776
         * in the ptra).  When finished, the ptra can be scanned for numa,
2777
         * and the original indices in the nas can be read out.  In this
2778
         * way, the ptra effectively sorts the input numbers in the nas. */
2779
0
    paindex = ptraCreate(isize + 1);
2780
0
    n = numaGetCount(nas);
2781
0
    for (i = 0; i < n; i++) {
2782
0
        numaGetIValue(nas, i, &ival);
2783
0
        nai = (NUMA *)ptraGetPtrToItem(paindex, ival);
2784
0
        if (!nai) {  /* make it; no shifting will occur */
2785
0
            nai = numaCreate(1);
2786
0
            ptraInsert(paindex, ival, nai, L_MIN_DOWNSHIFT);
2787
0
        }
2788
0
        numaAddNumber(nai, i);
2789
0
    }
2790
2791
        /* Sort by scanning the ptra, extracting numas and pulling
2792
         * the (index into nas) numbers out of each numa, taken
2793
         * successively in requested order. */
2794
0
    ptraGetMaxIndex(paindex, &imax);
2795
0
    nad = numaCreate(0);
2796
0
    if (sortorder == L_SORT_INCREASING) {
2797
0
        for (i = 0; i <= imax; i++) {
2798
0
            na = (NUMA *)ptraRemove(paindex, i, L_NO_COMPACTION);
2799
0
            if (!na) continue;
2800
0
            numaJoin(nad, na, 0, -1);
2801
0
            numaDestroy(&na);
2802
0
        }
2803
0
    } else {  /* L_SORT_DECREASING */
2804
0
        for (i = imax; i >= 0; i--) {
2805
0
            na = (NUMA *)ptraRemoveLast(paindex);
2806
0
            if (!na) break;  /* they've all been removed */
2807
0
            numaJoin(nad, na, 0, -1);
2808
0
            numaDestroy(&na);
2809
0
        }
2810
0
    }
2811
2812
0
    ptraDestroy(&paindex, FALSE, FALSE);
2813
0
    return nad;
2814
0
}
2815
2816
2817
/*!
2818
 * \brief   numaSortByIndex()
2819
 *
2820
 * \param[in]    nas
2821
 * \param[in]    naindex     na that maps from the new numa to the input numa
2822
 * \return  nad  sorted, or NULL on error
2823
 */
2824
NUMA *
2825
numaSortByIndex(NUMA  *nas,
2826
                NUMA  *naindex)
2827
0
{
2828
0
l_int32    i, n, ni, index;
2829
0
l_float32  val;
2830
0
NUMA      *nad;
2831
2832
0
    if (!nas)
2833
0
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
2834
0
    if (!naindex)
2835
0
        return (NUMA *)ERROR_PTR("naindex not defined", __func__, NULL);
2836
0
    n = numaGetCount(nas);
2837
0
    ni = numaGetCount(naindex);
2838
0
    if (n != ni)
2839
0
        return (NUMA *)ERROR_PTR("numa sizes differ", __func__, NULL);
2840
0
    if (n == 0) {
2841
0
        L_WARNING("nas is empty\n", __func__);
2842
0
        return numaCopy(nas);
2843
0
    }
2844
2845
0
    nad = numaCreate(n);
2846
0
    for (i = 0; i < n; i++) {
2847
0
        numaGetIValue(naindex, i, &index);
2848
0
        numaGetFValue(nas, index, &val);
2849
0
        numaAddNumber(nad, val);
2850
0
    }
2851
2852
0
    return nad;
2853
0
}
2854
2855
2856
/*!
2857
 * \brief   numaIsSorted()
2858
 *
2859
 * \param[in]    nas
2860
 * \param[in]    sortorder   L_SORT_INCREASING or L_SORT_DECREASING
2861
 * \param[out]   psorted     1 if sorted; 0 if not
2862
 * \return  1 if OK; 0 on error
2863
 *
2864
 * <pre>
2865
 * Notes:
2866
 *      (1) This is a quick O(n) test if nas is sorted.  It is useful
2867
 *          in situations where the array is likely to be already
2868
 *          sorted, and a sort operation can be avoided.
2869
 * </pre>
2870
 */
2871
l_int32
2872
numaIsSorted(NUMA     *nas,
2873
             l_int32   sortorder,
2874
             l_int32  *psorted)
2875
0
{
2876
0
l_int32    i, n;
2877
0
l_float32  prevval, val;
2878
2879
0
    if (!psorted)
2880
0
        return ERROR_INT("&sorted not defined", __func__, 1);
2881
0
    *psorted = FALSE;
2882
0
    if (!nas)
2883
0
        return ERROR_INT("nas not defined", __func__, 1);
2884
0
    if ((n = numaGetCount(nas))== 0) {
2885
0
        L_WARNING("nas is empty\n", __func__);
2886
0
        *psorted = TRUE;
2887
0
        return 0;
2888
0
    }
2889
0
    if (sortorder != L_SORT_INCREASING && sortorder != L_SORT_DECREASING)
2890
0
        return ERROR_INT("invalid sortorder", __func__, 1);
2891
2892
0
    n = numaGetCount(nas);
2893
0
    numaGetFValue(nas, 0, &prevval);
2894
0
    for (i = 1; i < n; i++) {
2895
0
        numaGetFValue(nas, i, &val);
2896
0
        if ((sortorder == L_SORT_INCREASING && val < prevval) ||
2897
0
            (sortorder == L_SORT_DECREASING && val > prevval))
2898
0
            return 0;
2899
0
    }
2900
2901
0
    *psorted = TRUE;
2902
0
    return 0;
2903
0
}
2904
2905
2906
/*!
2907
 * \brief   numaSortPair()
2908
 *
2909
 * \param[in]    nax, nay     input arrays
2910
 * \param[in]    sortorder    L_SORT_INCREASING or L_SORT_DECREASING
2911
 * \param[out]   pnasx        sorted
2912
 * \param[out]   pnasy        sorted exactly in order of nasx
2913
 * \return  0 if OK, 1 on error
2914
 *
2915
 * <pre>
2916
 * Notes:
2917
 *      (1) This function sorts the two input arrays, nax and nay,
2918
 *          together, using nax as the key for sorting.
2919
 * </pre>
2920
 */
2921
l_ok
2922
numaSortPair(NUMA    *nax,
2923
             NUMA    *nay,
2924
             l_int32  sortorder,
2925
             NUMA   **pnasx,
2926
             NUMA   **pnasy)
2927
0
{
2928
0
l_int32  sorted;
2929
0
NUMA    *naindex;
2930
2931
0
    if (pnasx) *pnasx = NULL;
2932
0
    if (pnasy) *pnasy = NULL;
2933
0
    if (!pnasx || !pnasy)
2934
0
        return ERROR_INT("&nasx and/or &nasy not defined", __func__, 1);
2935
0
    if (!nax)
2936
0
        return ERROR_INT("nax not defined", __func__, 1);
2937
0
    if (!nay)
2938
0
        return ERROR_INT("nay not defined", __func__, 1);
2939
0
    if (sortorder != L_SORT_INCREASING && sortorder != L_SORT_DECREASING)
2940
0
        return ERROR_INT("invalid sortorder", __func__, 1);
2941
2942
0
    numaIsSorted(nax, sortorder, &sorted);
2943
0
    if (sorted == TRUE) {
2944
0
        *pnasx = numaCopy(nax);
2945
0
        *pnasy = numaCopy(nay);
2946
0
    } else {
2947
0
        naindex = numaGetSortIndex(nax, sortorder);
2948
0
        *pnasx = numaSortByIndex(nax, naindex);
2949
0
        *pnasy = numaSortByIndex(nay, naindex);
2950
0
        numaDestroy(&naindex);
2951
0
    }
2952
2953
0
    return 0;
2954
0
}
2955
2956
2957
/*!
2958
 * \brief   numaInvertMap()
2959
 *
2960
 * \param[in]    nas
2961
 * \return  nad  the inverted map, or NULL on error or if not invertible
2962
 *
2963
 * <pre>
2964
 * Notes:
2965
 *      (1) This requires that nas contain each integer from 0 to n-1.
2966
 *          The array is typically an index array into a sort or permutation
2967
 *          of another array.
2968
 * </pre>
2969
 */
2970
NUMA *
2971
numaInvertMap(NUMA  *nas)
2972
0
{
2973
0
l_int32   i, n, val, error;
2974
0
l_int32  *test;
2975
0
NUMA     *nad;
2976
2977
0
    if (!nas)
2978
0
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
2979
0
    if ((n = numaGetCount(nas)) == 0) {
2980
0
        L_WARNING("nas is empty\n", __func__);
2981
0
        return numaCopy(nas);
2982
0
    }
2983
2984
0
    nad = numaMakeConstant(0.0, n);
2985
0
    test = (l_int32 *)LEPT_CALLOC(n, sizeof(l_int32));
2986
0
    error = 0;
2987
0
    for (i = 0; i < n; i++) {
2988
0
        numaGetIValue(nas, i, &val);
2989
0
        if (val >= n) {
2990
0
            error = 1;
2991
0
            break;
2992
0
        }
2993
0
        numaReplaceNumber(nad, val, i);
2994
0
        if (test[val] == 0) {
2995
0
            test[val] = 1;
2996
0
        } else {
2997
0
            error = 1;
2998
0
            break;
2999
0
        }
3000
0
    }
3001
3002
0
    LEPT_FREE(test);
3003
0
    if (error) {
3004
0
        numaDestroy(&nad);
3005
0
        return (NUMA *)ERROR_PTR("nas not invertible", __func__, NULL);
3006
0
    }
3007
3008
0
    return nad;
3009
0
}
3010
3011
/*!
3012
 * \brief   numaAddSorted()
3013
 *
3014
 * \param[in]    na     sorted input
3015
 * \param[in]    val    value to be inserted in sorted order
3016
 * \return  0 if OK, 1 on error
3017
 *
3018
 * <pre>
3019
 * Notes:
3020
 *      (1) The input %na is sorted.  This function determines the
3021
 *          sort order of %na and inserts %val into the array.
3022
 * </pre>
3023
 */
3024
l_ok
3025
numaAddSorted(NUMA      *na,
3026
              l_float32  val)
3027
0
{
3028
0
l_int32  index;
3029
3030
0
    if (!na)
3031
0
        return ERROR_INT("na not defined", __func__, 1);
3032
3033
0
    if (numaFindSortedLoc(na, val, &index) == 1)
3034
0
        return ERROR_INT("insert failure", __func__, 1);
3035
0
    numaInsertNumber(na, index, val);
3036
0
    return 0;
3037
0
}
3038
3039
3040
/*!
3041
 * \brief   numaFindSortedLoc()
3042
 *
3043
 * \param[in]    na     sorted input
3044
 * \param[in]    val    value to be inserted in sorted order
3045
 * \param[out]  *ploc   index location to insert @val
3046
 * \return  0 if OK, 1 on error
3047
 *
3048
 * <pre>
3049
 * Notes:
3050
 *      (1) The input %na is sorted.  This determines the sort order of @na,
3051
 *          either increasing or decreasing, and does a binary search for the
3052
 *          location to insert %val into the array.  The search is O(log n).
3053
 *      (2) The index returned is the location to insert into the array.
3054
 *          The value at the index, and all values to the right, are
3055
 *          moved to the right (increasing their index location by 1).
3056
 *      (3) If n is the size of %na, *ploc can be anything in [0 ... n].
3057
 *          if *ploc == 0, the value is inserted at the beginning of the
3058
 *          array; if *ploc == n, it is inserted at the end.
3059
 *      (4) If the size of %na is 1, insert with an increasing sort.
3060
 * </pre>
3061
 */
3062
l_ok
3063
numaFindSortedLoc(NUMA      *na,
3064
                  l_float32  val,
3065
                  l_int32   *pindex)
3066
0
{
3067
0
l_int32    n, increasing, lindex, rindex, midindex;
3068
0
l_float32  val0, valn, valmid;
3069
3070
0
    if (!pindex)
3071
0
        return ERROR_INT("&index not defined", __func__, 1);
3072
0
    *pindex = 0;
3073
0
    if (!na)
3074
0
        return ERROR_INT("na not defined", __func__, 1);
3075
3076
0
    n = numaGetCount(na);
3077
0
    if (n == 0) return 0;
3078
0
    numaGetFValue(na, 0, &val0);
3079
0
    if (n == 1) {  /* use increasing sort order */
3080
0
        if (val >= val0)
3081
0
            *pindex = 1;
3082
0
        return 0;
3083
0
    }
3084
3085
        /* -----------------  n >= 2 ----------------- */
3086
0
    numaGetFValue(na, n - 1, &valn);
3087
0
    increasing = (valn >= val0) ? 1 : 0;  /* sort order */
3088
3089
        /* Check if outside bounds of existing array */
3090
0
    if (increasing) {
3091
0
        if (val < val0) {
3092
0
            *pindex = 0;
3093
0
            return 0;
3094
0
        } else if (val > valn) {
3095
0
            *pindex = n;
3096
0
            return 0;
3097
0
        }
3098
0
    } else {  /* decreasing */
3099
0
        if (val > val0) {
3100
0
            *pindex = 0;
3101
0
            return 0;
3102
0
        } else if (val < valn) {
3103
0
            *pindex = n;
3104
0
            return 0;
3105
0
        }
3106
0
    }
3107
3108
        /* Within bounds of existing array; search */
3109
0
    lindex = 0;
3110
0
    rindex = n - 1;
3111
0
    while (1) {
3112
0
        midindex = (lindex + rindex) / 2;
3113
0
        if (midindex == lindex || midindex == rindex) break;
3114
0
        numaGetFValue(na, midindex, &valmid);
3115
0
        if (increasing) {
3116
0
            if (val > valmid)
3117
0
                lindex = midindex;
3118
0
            else
3119
0
                rindex = midindex;
3120
0
        } else {  /* decreasing */
3121
0
            if (val > valmid)
3122
0
                rindex = midindex;
3123
0
            else
3124
0
                lindex = midindex;
3125
0
        }
3126
0
    }
3127
0
    *pindex = rindex;
3128
0
    return 0;
3129
0
}
3130
3131
3132
/*----------------------------------------------------------------------*
3133
 *                          Random permutation                          *
3134
 *----------------------------------------------------------------------*/
3135
/*!
3136
 * \brief   numaPseudorandomSequence()
3137
 *
3138
 * \param[in]    size     of sequence
3139
 * \param[in]    seed     for random number generation
3140
 * \return  na  pseudorandom on [0,...,size - 1], or NULL on error
3141
 *
3142
 * <pre>
3143
 * Notes:
3144
 *      (1) This uses the Durstenfeld shuffle.
3145
 *          See: http://en.wikipedia.org/wiki/Fisher–Yates_shuffle.
3146
 *          Result is a pseudorandom permutation of the sequence of integers
3147
 *          from 0 to size - 1.
3148
 * </pre>
3149
 */
3150
NUMA *
3151
numaPseudorandomSequence(l_int32  size,
3152
                         l_int32  seed)
3153
0
{
3154
0
l_int32   i, index, temp;
3155
0
l_int32  *array;
3156
0
NUMA     *na;
3157
3158
0
    if (size <= 0)
3159
0
        return (NUMA *)ERROR_PTR("size <= 0", __func__, NULL);
3160
3161
0
    if ((array = (l_int32 *)LEPT_CALLOC(size, sizeof(l_int32))) == NULL)
3162
0
        return (NUMA *)ERROR_PTR("array not made", __func__, NULL);
3163
0
    for (i = 0; i < size; i++)
3164
0
        array[i] = i;
3165
0
    srand(seed);
3166
0
    for (i = size - 1; i > 0; i--) {
3167
0
        index = (l_int32)((i + 1) * ((l_float64)rand() / (l_float64)RAND_MAX));
3168
0
        index = L_MIN(index, i);
3169
0
        temp = array[i];
3170
0
        array[i] = array[index];
3171
0
        array[index] = temp;
3172
0
    }
3173
3174
0
    na = numaCreateFromIArray(array, size);
3175
0
    LEPT_FREE(array);
3176
0
    return na;
3177
0
}
3178
3179
3180
/*!
3181
 * \brief   numaRandomPermutation()
3182
 *
3183
 * \param[in]    nas    input array
3184
 * \param[in]    seed   for random number generation
3185
 * \return  nas  randomly shuffled array, or NULL on error
3186
 */
3187
NUMA *
3188
numaRandomPermutation(NUMA    *nas,
3189
                      l_int32  seed)
3190
0
{
3191
0
l_int32    i, index, size;
3192
0
l_float32  val;
3193
0
NUMA      *naindex, *nad;
3194
3195
0
    if (!nas)
3196
0
        return (NUMA *)ERROR_PTR("nas not defined", __func__, NULL);
3197
0
    if ((size = numaGetCount(nas)) == 0) {
3198
0
        L_WARNING("nas is empty\n", __func__);
3199
0
        return numaCopy(nas);
3200
0
    }
3201
3202
0
    naindex = numaPseudorandomSequence(size, seed);
3203
0
    nad = numaCreate(size);
3204
0
    for (i = 0; i < size; i++) {
3205
0
        numaGetIValue(naindex, i, &index);
3206
0
        numaGetFValue(nas, index, &val);
3207
0
        numaAddNumber(nad, val);
3208
0
    }
3209
0
    numaDestroy(&naindex);
3210
0
    return nad;
3211
0
}
3212
3213
3214
/*----------------------------------------------------------------------*
3215
 *                     Functions requiring sorting                      *
3216
 *----------------------------------------------------------------------*/
3217
/*!
3218
 * \brief   numaGetRankValue()
3219
 *
3220
 * \param[in]    na        source numa
3221
 * \param[in]    fract     use 0.0 for smallest, 1.0 for largest
3222
 * \param[in]    nasort    [optional] increasing sorted version of na
3223
 * \param[in]    usebins   0 for general sort; 1 for bin sort
3224
 * \param[out]   pval      rank val
3225
 * \return  0 if OK; 1 on error
3226
 *
3227
 * <pre>
3228
 * Notes:
3229
 *      (1) Computes the rank value of a number in the %na, which is
3230
 *          the number that is a fraction %fract from the small
3231
 *          end of the sorted version of %na.
3232
 *      (2) If you do this multiple times for different rank values,
3233
 *          sort the array in advance and use that for %nasort;
3234
 *          if you're only calling this once, input %nasort == NULL.
3235
 *      (3) If %usebins == 1, this uses a bin sorting method.
3236
 *          Use this only where:
3237
 *           * the numbers are non-negative integers
3238
 *           * there are over 100 numbers
3239
 *           * the maximum value is less than about 50,000
3240
 *      (4) The advantage of using a bin sort is that it is O(n),
3241
 *          instead of O(nlogn) for general sort routines.
3242
 * </pre>
3243
 */
3244
l_ok
3245
numaGetRankValue(NUMA       *na,
3246
                 l_float32   fract,
3247
                 NUMA       *nasort,
3248
                 l_int32     usebins,
3249
                 l_float32  *pval)
3250
0
{
3251
0
l_int32  n, index;
3252
0
NUMA    *nas;
3253
3254
0
    if (!pval)
3255
0
        return ERROR_INT("&val not defined", __func__, 1);
3256
0
    *pval = 0.0;  /* init */
3257
0
    if (!na)
3258
0
        return ERROR_INT("na not defined", __func__, 1);
3259
0
    if ((n = numaGetCount(na)) == 0)
3260
0
        return ERROR_INT("na empty", __func__, 1);
3261
0
    if (fract < 0.0 || fract > 1.0)
3262
0
        return ERROR_INT("fract not in [0.0 ... 1.0]", __func__, 1);
3263
3264
0
    if (nasort) {
3265
0
        nas = nasort;
3266
0
    } else {
3267
0
        if (usebins == 0)
3268
0
            nas = numaSort(NULL, na, L_SORT_INCREASING);
3269
0
        else
3270
0
            nas = numaBinSort(na, L_SORT_INCREASING);
3271
0
        if (!nas)
3272
0
            return ERROR_INT("nas not made", __func__, 1);
3273
0
    }
3274
0
    index = (l_int32)(fract * (l_float32)(n - 1) + 0.5);
3275
0
    numaGetFValue(nas, index, pval);
3276
3277
0
    if (!nasort) numaDestroy(&nas);
3278
0
    return 0;
3279
0
}
3280
3281
3282
/*!
3283
 * \brief   numaGetMedian()
3284
 *
3285
 * \param[in]    na     source numa
3286
 * \param[out]   pval   median value
3287
 * \return  0 if OK; 1 on error
3288
 *
3289
 * <pre>
3290
 * Notes:
3291
 *      (1) Computes the median value of the numbers in the numa, by
3292
 *          sorting and finding the middle value in the sorted array.
3293
 * </pre>
3294
 */
3295
l_ok
3296
numaGetMedian(NUMA       *na,
3297
              l_float32  *pval)
3298
0
{
3299
0
    if (!pval)
3300
0
        return ERROR_INT("&val not defined", __func__, 1);
3301
0
    *pval = 0.0;  /* init */
3302
0
    if (!na || numaGetCount(na) == 0)
3303
0
        return ERROR_INT("na not defined or empty", __func__, 1);
3304
3305
0
    return numaGetRankValue(na, 0.5, NULL, 0, pval);
3306
0
}
3307
3308
3309
/*!
3310
 * \brief   numaGetBinnedMedian()
3311
 *
3312
 * \param[in]    na      source numa
3313
 * \param[out]   pval    integer median value
3314
 * \return  0 if OK; 1 on error
3315
 *
3316
 * <pre>
3317
 * Notes:
3318
 *      (1) Computes the median value of the numbers in the numa,
3319
 *          using bin sort and finding the middle value in the sorted array.
3320
 *      (2) See numaGetRankValue() for conditions on na for which
3321
 *          this should be used.  Otherwise, use numaGetMedian().
3322
 * </pre>
3323
 */
3324
l_ok
3325
numaGetBinnedMedian(NUMA     *na,
3326
                    l_int32  *pval)
3327
0
{
3328
0
l_int32    ret;
3329
0
l_float32  fval;
3330
3331
0
    if (!pval)
3332
0
        return ERROR_INT("&val not defined", __func__, 1);
3333
0
    *pval = 0;  /* init */
3334
0
    if (!na || numaGetCount(na) == 0)
3335
0
        return ERROR_INT("na not defined or empty", __func__, 1);
3336
3337
0
    ret = numaGetRankValue(na, 0.5, NULL, 1, &fval);
3338
0
    *pval = lept_roundftoi(fval);
3339
0
    return ret;
3340
0
}
3341
3342
3343
/*!
3344
 * \brief   numaGetMeanDevFromMedian()
3345
 *
3346
 * \param[in]      na     source numa
3347
 * \param[in]      med    median value
3348
 * \param[out]     pdev   average absolute value deviation from median value
3349
 * \return  0 if OK; 1 on error
3350
 */
3351
l_ok
3352
numaGetMeanDevFromMedian(NUMA       *na,
3353
                         l_float32   med,
3354
                         l_float32  *pdev)
3355
0
{
3356
0
l_int32    i, n;
3357
0
l_float32  val, dev;
3358
3359
0
    if (!pdev)
3360
0
        return ERROR_INT("&dev not defined", __func__, 1);
3361
0
    *pdev = 0.0;  /* init */
3362
0
    if (!na)
3363
0
        return ERROR_INT("na not defined", __func__, 1);
3364
0
    if ((n = numaGetCount(na)) == 0)
3365
0
        return ERROR_INT("na is empty", __func__, 1);
3366
3367
0
    dev = 0.0;
3368
0
    for (i = 0; i < n; i++) {
3369
0
        numaGetFValue(na, i, &val);
3370
0
        dev += L_ABS(val - med);
3371
0
    }
3372
0
    *pdev = dev / (l_float32)n;
3373
0
    return 0;
3374
0
}
3375
3376
3377
/*!
3378
 * \brief   numaGetMedianDevFromMedian()
3379
 *
3380
 * \param[in]    na        source numa
3381
 * \param[out]   pmed      [optional] median value
3382
 * \param[out]   pdev      median deviation from median val
3383
 * \return  0 if OK; 1 on error
3384
 *
3385
 * <pre>
3386
 * Notes:
3387
 *      (1) Finds the median of the absolute value of the deviation from
3388
 *          the median value in the array.  Why take the absolute value?
3389
 *          Consider the case where you have values equally distributed
3390
 *          about both sides of a median value.  Without taking the absolute
3391
 *          value of the differences, you will get 0 for the deviation,
3392
 *          and this is not useful.
3393
 * </pre>
3394
 */
3395
l_ok
3396
numaGetMedianDevFromMedian(NUMA       *na,
3397
                           l_float32  *pmed,
3398
                           l_float32  *pdev)
3399
0
{
3400
0
l_int32    n, i;
3401
0
l_float32  val, med;
3402
0
NUMA      *nadev;
3403
3404
0
    if (pmed) *pmed = 0.0;
3405
0
    if (!pdev)
3406
0
        return ERROR_INT("&dev not defined", __func__, 1);
3407
0
    *pdev = 0.0;
3408
0
    if (!na || numaGetCount(na) == 0)
3409
0
        return ERROR_INT("na not defined or empty", __func__, 1);
3410
3411
0
    numaGetMedian(na, &med);
3412
0
    if (pmed) *pmed = med;
3413
0
    n = numaGetCount(na);
3414
0
    nadev = numaCreate(n);
3415
0
    for (i = 0; i < n; i++) {
3416
0
        numaGetFValue(na, i, &val);
3417
0
        numaAddNumber(nadev, L_ABS(val - med));
3418
0
    }
3419
0
    numaGetMedian(nadev, pdev);
3420
3421
0
    numaDestroy(&nadev);
3422
0
    return 0;
3423
0
}
3424
3425
3426
/*!
3427
 * \brief   numaGetMode()
3428
 *
3429
 * \param[in]    na      source numa
3430
 * \param[out]   pval    mode val
3431
 * \param[out]   pcount  [optional] mode count
3432
 * \return  0 if OK; 1 on error
3433
 *
3434
 * <pre>
3435
 * Notes:
3436
 *      (1) Computes the mode value of the numbers in the numa, by
3437
 *          sorting and finding the value of the number with the
3438
 *          largest count.
3439
 *      (2) Optionally, also returns that count.
3440
 * </pre>
3441
 */
3442
l_ok
3443
numaGetMode(NUMA       *na,
3444
            l_float32  *pval,
3445
            l_int32    *pcount)
3446
0
{
3447
0
l_int32     i, n, maxcount, prevcount;
3448
0
l_float32   val, maxval, prevval;
3449
0
l_float32  *array;
3450
0
NUMA       *nasort;
3451
3452
0
    if (pcount) *pcount = 0;
3453
0
    if (!pval)
3454
0
        return ERROR_INT("&val not defined", __func__, 1);
3455
0
    *pval = 0.0;
3456
0
    if (!na)
3457
0
        return ERROR_INT("na not defined", __func__, 1);
3458
0
    if ((n = numaGetCount(na)) == 0)
3459
0
        return ERROR_INT("na is empty", __func__, 1);
3460
3461
0
    if ((nasort = numaSort(NULL, na, L_SORT_DECREASING)) == NULL)
3462
0
        return ERROR_INT("nas not made", __func__, 1);
3463
0
    array = numaGetFArray(nasort, L_NOCOPY);
3464
3465
        /* Initialize with array[0] */
3466
0
    prevval = array[0];
3467
0
    prevcount = 1;
3468
0
    maxval = prevval;
3469
0
    maxcount = prevcount;
3470
3471
        /* Scan the sorted array, aggregating duplicates */
3472
0
    for (i = 1; i < n; i++) {
3473
0
        val = array[i];
3474
0
        if (val == prevval) {
3475
0
            prevcount++;
3476
0
        } else {  /* new value */
3477
0
            if (prevcount > maxcount) {  /* new max */
3478
0
                maxcount = prevcount;
3479
0
                maxval = prevval;
3480
0
            }
3481
0
            prevval = val;
3482
0
            prevcount = 1;
3483
0
        }
3484
0
    }
3485
3486
        /* Was the mode the last run of elements? */
3487
0
    if (prevcount > maxcount) {
3488
0
        maxcount = prevcount;
3489
0
        maxval = prevval;
3490
0
    }
3491
3492
0
    *pval = maxval;
3493
0
    if (pcount)
3494
0
        *pcount = maxcount;
3495
3496
0
    numaDestroy(&nasort);
3497
0
    return 0;
3498
0
}
3499
3500
3501
/*----------------------------------------------------------------------*
3502
 *                            Rearrangements                            *
3503
 *----------------------------------------------------------------------*/
3504
/*!
3505
 * \brief   numaJoin()
3506
 *
3507
 * \param[in]    nad      dest numa; add to this one
3508
 * \param[in]    nas      [optional] source numa; add from this one
3509
 * \param[in]    istart   starting index in nas
3510
 * \param[in]    iend     ending index in nas; use -1 to cat all
3511
 * \return  0 if OK, 1 on error
3512
 *
3513
 * <pre>
3514
 * Notes:
3515
 *      (1) istart < 0 is taken to mean 'read from the start' (istart = 0)
3516
 *      (2) iend < 0 means 'read to the end'
3517
 *      (3) if nas == NULL, this is a no-op
3518
 * </pre>
3519
 */
3520
l_ok
3521
numaJoin(NUMA    *nad,
3522
         NUMA    *nas,
3523
         l_int32  istart,
3524
         l_int32  iend)
3525
0
{
3526
0
l_int32    n, i;
3527
0
l_float32  val;
3528
3529
0
    if (!nad)
3530
0
        return ERROR_INT("nad not defined", __func__, 1);
3531
0
    if (!nas)
3532
0
        return 0;
3533
3534
0
    if (istart < 0)
3535
0
        istart = 0;
3536
0
    n = numaGetCount(nas);
3537
0
    if (iend < 0 || iend >= n)
3538
0
        iend = n - 1;
3539
0
    if (istart > iend)
3540
0
        return ERROR_INT("istart > iend; nothing to add", __func__, 1);
3541
3542
0
    for (i = istart; i <= iend; i++) {
3543
0
        numaGetFValue(nas, i, &val);
3544
0
        numaAddNumber(nad, val);
3545
0
    }
3546
3547
0
    return 0;
3548
0
}
3549
3550
3551
/*!
3552
 * \brief   numaaJoin()
3553
 *
3554
 * \param[in]    naad     dest naa; add to this one
3555
 * \param[in]    naas     [optional] source naa; add from this one
3556
 * \param[in]    istart   starting index in nas
3557
 * \param[in]    iend     ending index in naas; use -1 to cat all
3558
 * \return  0 if OK, 1 on error
3559
 *
3560
 * <pre>
3561
 * Notes:
3562
 *      (1) istart < 0 is taken to mean 'read from the start' (istart = 0)
3563
 *      (2) iend < 0 means 'read to the end'
3564
 *      (3) if naas == NULL, this is a no-op
3565
 * </pre>
3566
 */
3567
l_ok
3568
numaaJoin(NUMAA   *naad,
3569
          NUMAA   *naas,
3570
          l_int32  istart,
3571
          l_int32  iend)
3572
0
{
3573
0
l_int32  n, i;
3574
0
NUMA    *na;
3575
3576
0
    if (!naad)
3577
0
        return ERROR_INT("naad not defined", __func__, 1);
3578
0
    if (!naas)
3579
0
        return 0;
3580
3581
0
    if (istart < 0)
3582
0
        istart = 0;
3583
0
    n = numaaGetCount(naas);
3584
0
    if (iend < 0 || iend >= n)
3585
0
        iend = n - 1;
3586
0
    if (istart > iend)
3587
0
        return ERROR_INT("istart > iend; nothing to add", __func__, 1);
3588
3589
0
    for (i = istart; i <= iend; i++) {
3590
0
        na = numaaGetNuma(naas, i, L_CLONE);
3591
0
        numaaAddNuma(naad, na, L_INSERT);
3592
0
    }
3593
3594
0
    return 0;
3595
0
}
3596
3597
3598
/*!
3599
 * \brief   numaaFlattenToNuma()
3600
 *
3601
 * \param[in]    naa
3602
 * \return  numa, or NULL on error
3603
 *
3604
 * <pre>
3605
 * Notes:
3606
 *      (1) This 'flattens' the Numaa to a Numa, by joining successively
3607
 *          each Numa in the Numaa.
3608
 *      (2) It doesn't make any assumptions about the location of the
3609
 *          Numas in the Numaa array, unlike most Numaa functions.
3610
 *      (3) It leaves the input Numaa unchanged.
3611
 * </pre>
3612
 */
3613
NUMA *
3614
numaaFlattenToNuma(NUMAA  *naa)
3615
0
{
3616
0
l_int32  i, nalloc;
3617
0
NUMA    *na, *nad;
3618
0
NUMA   **array;
3619
3620
0
    if (!naa)
3621
0
        return (NUMA *)ERROR_PTR("naa not defined", __func__, NULL);
3622
3623
0
    nalloc = naa->nalloc;
3624
0
    array = numaaGetPtrArray(naa);
3625
0
    nad = numaCreate(0);
3626
0
    for (i = 0; i < nalloc; i++) {
3627
0
        na = array[i];
3628
0
        if (!na) continue;
3629
0
        numaJoin(nad, na, 0, -1);
3630
0
    }
3631
3632
0
    return nad;
3633
0
}
3634