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