1 /*====================================================================*
2 - Copyright (C) 2001 Leptonica. All rights reserved.
3 - This software is distributed in the hope that it will be
4 - useful, but with NO WARRANTY OF ANY KIND.
5 - No author or distributor accepts responsibility to anyone for the
6 - consequences of using this software, or for whether it serves any
7 - particular purpose or works at all, unless he or she says so in
8 - writing. Everyone is granted permission to copy, modify and
9 - redistribute this source code, for commercial or non-commercial
10 - purposes, with the following restrictions: (1) the origin of this
11 - source code must not be misrepresented; (2) modified versions must
12 - be plainly marked as such; and (3) this notice may not be removed
13 - or altered from any source or modified source distribution.
14 *====================================================================*/
15
16 /*
17 * numafunc1.c
18 *
19 * Arithmetic and logic
20 * NUMA *numaArithOp()
21 * NUMA *numaLogicalOp()
22 * NUMA *numaInvert()
23 *
24 * Simple extractions
25 * l_int32 numaGetMin()
26 * l_int32 numaGetMax()
27 * l_int32 numaGetSum()
28 * NUMA *numaGetPartialSums()
29 * l_int32 numaGetSumOnInterval()
30 * l_int32 numaHasOnlyIntegers()
31 * NUMA *numaSubsample()
32 * NUMA *numaMakeSequence()
33 * NUMA *numaMakeConstant()
34 * l_int32 numaGetNonzeroRange()
35 * l_int32 numaGetCountRelativeToZero()
36 * NUMA *numaClipToInterval()
37 * NUMA *numaMakeThresholdIndicator()
38 * NUMA *numaUniformSampling()
39 *
40 * Interpolation
41 * l_int32 numaInterpolateEqxVal()
42 * l_int32 numaInterpolateEqxInterval()
43 * l_int32 numaInterpolateArbxVal()
44 * l_int32 numaInterpolateArbxInterval()
45 *
46 * Functions requiring interpolation
47 * l_int32 numaFitMax()
48 * l_int32 numaDifferentiateInterval()
49 * l_int32 numaIntegrateInterval()
50 *
51 * Sorting
52 * NUMA *numaSort()
53 * NUMA *numaGetSortIndex()
54 * NUMA *numaSortByIndex()
55 * l_int32 numaIsSorted()
56 * l_int32 numaSortPair()
57 * NUMA *numaPseudorandomSequence()
58 *
59 * Functions requiring sorting
60 * l_int32 numaGetRankValue()
61 * l_int32 numaGetMedian()
62 * l_int32 numaGetMode()
63 *
64 * Numa combination
65 * l_int32 numaJoin()
66 * NUMA *numaaFlattenToNuma()
67 *
68 *
69 * Things to remember when using the Numa:
70 *
71 * (1) The numa is a struct, not an array. Always use accessors
72 * (see numabasic.c), never the fields directly.
73 *
74 * (2) The number array holds l_float32 values. It can also
75 * be used to store l_int32 values. See numabasic.c for
76 * details on using the accessors.
77 *
78 * (3) If you use numaCreate(), no numbers are stored and the size is 0.
79 * You have to add numbers to increase the size.
80 * If you want to start with a numa of a fixed size, with each
81 * entry initialized to the same value, use numaMakeConstant().
82 *
83 * (4) Occasionally, in the comments we denote the i-th element of a
84 * numa by na[i]. This is conceptual only -- the numa is not an array!
85 */
86
87 #include <stdio.h>
88 #include <string.h>
89 #include <stdlib.h>
90 #include <math.h>
91 #include "allheaders.h"
92
93
94 /*----------------------------------------------------------------------*
95 * Arithmetic and logical ops on Numas *
96 *----------------------------------------------------------------------*/
97 /*!
98 * numaArithOp()
99 *
100 * Input: nad (<optional> can be null or equal to na1 (in-place)
101 * na1
102 * na2
103 * op (L_ARITH_ADD, L_ARITH_SUBTRACT,
104 * L_ARITH_MULTIPLY, L_ARITH_DIVIDE)
105 * Return: nad (always: operation applied to na1 and na2)
106 *
107 * Notes:
108 * (1) The sizes of na1 and na2 must be equal.
109 * (2) nad can only null or equal to na1.
110 * (3) To add a constant to a numa, or to multipy a numa by
111 * a constant, use numaTransform().
112 */
113 NUMA *
numaArithOp(NUMA * nad,NUMA * na1,NUMA * na2,l_int32 op)114 numaArithOp(NUMA *nad,
115 NUMA *na1,
116 NUMA *na2,
117 l_int32 op)
118 {
119 l_int32 i, n;
120 l_float32 val1, val2;
121
122 PROCNAME("numaArithOp");
123
124 if (!na1 || !na2)
125 return (NUMA *)ERROR_PTR("na1, na2 not both defined", procName, nad);
126 n = numaGetCount(na1);
127 if (n != numaGetCount(na2))
128 return (NUMA *)ERROR_PTR("na1, na2 sizes differ", procName, nad);
129 if (nad && nad != na1)
130 return (NUMA *)ERROR_PTR("nad defined but not in-place", procName, nad);
131 if (op != L_ARITH_ADD && op != L_ARITH_SUBTRACT &&
132 op != L_ARITH_MULTIPLY && op != L_ARITH_DIVIDE)
133 return (NUMA *)ERROR_PTR("invalid op", procName, nad);
134 if (op == L_ARITH_DIVIDE) {
135 for (i = 0; i < n; i++) {
136 numaGetFValue(na2, i, &val2);
137 if (val2 == 0.0)
138 return (NUMA *)ERROR_PTR("na2 has 0 element", procName, nad);
139 }
140 }
141
142 /* If nad is not identical to na1, make it an identical copy */
143 if (!nad)
144 nad = numaCopy(na1);
145
146 for (i = 0; i < n; i++) {
147 numaGetFValue(nad, i, &val1);
148 numaGetFValue(na2, i, &val2);
149 switch (op) {
150 case L_ARITH_ADD:
151 numaSetValue(nad, i, val1 + val2);
152 break;
153 case L_ARITH_SUBTRACT:
154 numaSetValue(nad, i, val1 - val2);
155 break;
156 case L_ARITH_MULTIPLY:
157 numaSetValue(nad, i, val1 * val2);
158 break;
159 case L_ARITH_DIVIDE:
160 numaSetValue(nad, i, val1 / val2);
161 break;
162 default:
163 fprintf(stderr, " Unknown arith op: %d\n", op);
164 return nad;
165 }
166 }
167
168 return nad;
169 }
170
171
172 /*!
173 * numaLogicalOp()
174 *
175 * Input: nad (<optional> can be null or equal to na1 (in-place)
176 * na1
177 * na2
178 * op (L_UNION, L_INTERSECTION, L_SUBTRACTION, L_EXCLUSIVE_OR)
179 * Return: nad (always: operation applied to na1 and na2)
180 *
181 * Notes:
182 * (1) The sizes of na1 and na2 must be equal.
183 * (2) nad can only null or equal to na1.
184 * (3) This is intended for use with indicator arrays (0s and 1s).
185 * Input data is extracted as integers (0 == false, anything
186 * else == true); output results are 0 and 1.
187 * (4) L_SUBTRACTION is subtraction of val2 from val1. For bit logical
188 * arithmetic this is (val1 & ~val2), but because these values
189 * are integers, we use (val1 && !val2).
190 */
191 NUMA *
numaLogicalOp(NUMA * nad,NUMA * na1,NUMA * na2,l_int32 op)192 numaLogicalOp(NUMA *nad,
193 NUMA *na1,
194 NUMA *na2,
195 l_int32 op)
196 {
197 l_int32 i, n, val1, val2, val;
198
199 PROCNAME("numaLogicalOp");
200
201 if (!na1 || !na2)
202 return (NUMA *)ERROR_PTR("na1, na2 not both defined", procName, nad);
203 n = numaGetCount(na1);
204 if (n != numaGetCount(na2))
205 return (NUMA *)ERROR_PTR("na1, na2 sizes differ", procName, nad);
206 if (nad && nad != na1)
207 return (NUMA *)ERROR_PTR("nad defined; not in-place", procName, nad);
208 if (op != L_UNION && op != L_INTERSECTION &&
209 op != L_SUBTRACTION && op != L_EXCLUSIVE_OR)
210 return (NUMA *)ERROR_PTR("invalid op", procName, nad);
211
212 /* If nad is not identical to na1, make it an identical copy */
213 if (!nad)
214 nad = numaCopy(na1);
215
216 for (i = 0; i < n; i++) {
217 numaGetIValue(nad, i, &val1);
218 numaGetIValue(na2, i, &val2);
219 switch (op) {
220 case L_UNION:
221 val = (val1 || val2) ? 1 : 0;
222 numaSetValue(nad, i, val);
223 break;
224 case L_INTERSECTION:
225 val = (val1 && val2) ? 1 : 0;
226 numaSetValue(nad, i, val);
227 break;
228 case L_SUBTRACTION:
229 val = (val1 && !val2) ? 1 : 0;
230 numaSetValue(nad, i, val);
231 break;
232 case L_EXCLUSIVE_OR:
233 val = ((val1 && !val2) || (!val1 && val2)) ? 1 : 0;
234 numaSetValue(nad, i, val);
235 break;
236 default:
237 fprintf(stderr, " Unknown logical op: %d\n", op);
238 return nad;
239 }
240 }
241
242 return nad;
243 }
244
245
246 /*!
247 * numaInvert()
248 *
249 * Input: nad (<optional> can be null or equal to nas (in-place)
250 * nas
251 * Return: nad (always: 'inverts' nas)
252 *
253 * Notes:
254 * (1) This is intended for use with indicator arrays (0s and 1s).
255 * It gives a boolean-type output, taking the input as
256 * an integer and inverting it:
257 * 0 --> 1
258 * anything else --> 0
259 */
260 NUMA *
numaInvert(NUMA * nad,NUMA * nas)261 numaInvert(NUMA *nad,
262 NUMA *nas)
263 {
264 l_int32 i, n, val;
265
266 PROCNAME("numaInvert");
267
268 if (!nas)
269 return (NUMA *)ERROR_PTR("nas not defined", procName, nad);
270 if (nad && nad != nas)
271 return (NUMA *)ERROR_PTR("nad defined; not in-place", procName, nad);
272
273 if (!nad)
274 nad = numaCopy(nas);
275 n = numaGetCount(nad);
276 for (i = 0; i < n; i++) {
277 numaGetIValue(nad, i, &val);
278 if (!val)
279 val = 1;
280 else
281 val = 0;
282 numaSetValue(nad, i, val);
283 }
284
285 return nad;
286 }
287
288
289 /*----------------------------------------------------------------------*
290 * Simple extractions *
291 *----------------------------------------------------------------------*/
292 /*!
293 * numaGetMin()
294 *
295 * Input: na
296 * &minval (<optional return> min value)
297 * &iminloc (<optional return> index of min location)
298 * Return: 0 if OK; 1 on error
299 */
300 l_int32
numaGetMin(NUMA * na,l_float32 * pminval,l_int32 * piminloc)301 numaGetMin(NUMA *na,
302 l_float32 *pminval,
303 l_int32 *piminloc)
304 {
305 l_int32 i, n, iminloc;
306 l_float32 val, minval;
307
308 PROCNAME("numaGetMin");
309
310 if (!pminval && !piminloc)
311 return ERROR_INT("nothing to do", procName, 1);
312 if (pminval) *pminval = 0.0;
313 if (piminloc) *piminloc = 0;
314 if (!na)
315 return ERROR_INT("na not defined", procName, 1);
316
317 minval = +1000000000.;
318 iminloc = 0;
319 n = numaGetCount(na);
320 for (i = 0; i < n; i++) {
321 numaGetFValue(na, i, &val);
322 if (val < minval) {
323 minval = val;
324 iminloc = i;
325 }
326 }
327
328 if (pminval) *pminval = minval;
329 if (piminloc) *piminloc = iminloc;
330 return 0;
331 }
332
333
334 /*!
335 * numaGetMax()
336 *
337 * Input: na
338 * &maxval (<optional return> max value)
339 * &imaxloc (<optional return> index of max location)
340 * Return: 0 if OK; 1 on error
341 */
342 l_int32
numaGetMax(NUMA * na,l_float32 * pmaxval,l_int32 * pimaxloc)343 numaGetMax(NUMA *na,
344 l_float32 *pmaxval,
345 l_int32 *pimaxloc)
346 {
347 l_int32 i, n, imaxloc;
348 l_float32 val, maxval;
349
350 PROCNAME("numaGetMax");
351
352 if (!pmaxval && !pimaxloc)
353 return ERROR_INT("nothing to do", procName, 1);
354 if (pmaxval) *pmaxval = 0.0;
355 if (pimaxloc) *pimaxloc = 0;
356 if (!na)
357 return ERROR_INT("na not defined", procName, 1);
358
359 maxval = -1000000000.;
360 imaxloc = 0;
361 n = numaGetCount(na);
362 for (i = 0; i < n; i++) {
363 numaGetFValue(na, i, &val);
364 if (val > maxval) {
365 maxval = val;
366 imaxloc = i;
367 }
368 }
369
370 if (pmaxval) *pmaxval = maxval;
371 if (pimaxloc) *pimaxloc = imaxloc;
372 return 0;
373 }
374
375
376 /*!
377 * numaGetSum()
378 *
379 * Input: na
380 * &sum (<return> sum of values)
381 * Return: 0 if OK, 1 on error
382 */
383 l_int32
numaGetSum(NUMA * na,l_float32 * psum)384 numaGetSum(NUMA *na,
385 l_float32 *psum)
386 {
387 l_int32 i, n;
388 l_float32 val, sum;
389
390 PROCNAME("numaGetSum");
391
392 if (!na)
393 return ERROR_INT("na not defined", procName, 1);
394 if (!psum)
395 return ERROR_INT("&sum not defined", procName, 1);
396
397 sum = 0.0;
398 n = numaGetCount(na);
399 for (i = 0; i < n; i++) {
400 numaGetFValue(na, i, &val);
401 sum += val;
402 }
403 *psum = sum;
404 return 0;
405 }
406
407
408 /*!
409 * numaGetPartialSums()
410 *
411 * Input: na
412 * Return: nasum, or null on error
413 *
414 * Notes:
415 * (1) nasum[i] is the sum for all j <= i of na[j].
416 * So nasum[0] = na[0].
417 * (2) If you want to generate a rank function, where rank[0] - 0.0,
418 * insert a 0.0 at the beginning of the nasum array.
419 */
420 NUMA *
numaGetPartialSums(NUMA * na)421 numaGetPartialSums(NUMA *na)
422 {
423 l_int32 i, n;
424 l_float32 val, sum;
425 NUMA *nasum;
426
427 PROCNAME("numaGetPartialSums");
428
429 if (!na)
430 return (NUMA *)ERROR_PTR("na not defined", procName, NULL);
431
432 n = numaGetCount(na);
433 nasum = numaCreate(n);
434 sum = 0.0;
435 for (i = 0; i < n; i++) {
436 numaGetFValue(na, i, &val);
437 sum += val;
438 numaAddNumber(nasum, sum);
439 }
440 return nasum;
441 }
442
443
444 /*!
445 * numaGetSumOnInterval()
446 *
447 * Input: na
448 * first (beginning index)
449 * last (final index)
450 * &sum (<return> sum of values in the index interval range)
451 * Return: 0 if OK, 1 on error
452 */
453 l_int32
numaGetSumOnInterval(NUMA * na,l_int32 first,l_int32 last,l_float32 * psum)454 numaGetSumOnInterval(NUMA *na,
455 l_int32 first,
456 l_int32 last,
457 l_float32 *psum)
458 {
459 l_int32 i, n, truelast;
460 l_float32 val, sum;
461
462 PROCNAME("numaGetSumOnInterval");
463
464 if (!na)
465 return ERROR_INT("na not defined", procName, 1);
466 if (!psum)
467 return ERROR_INT("&sum not defined", procName, 1);
468 *psum = 0.0;
469
470 sum = 0.0;
471 n = numaGetCount(na);
472 if (first >= n) /* not an error */
473 return 0;
474 truelast = L_MIN(last, n - 1);
475
476 for (i = first; i <= truelast; i++) {
477 numaGetFValue(na, i, &val);
478 sum += val;
479 }
480 *psum = sum;
481 return 0;
482 }
483
484
485 /*!
486 * numaHasOnlyIntegers()
487 *
488 * Input: na
489 * maxsamples (maximum number of samples to check)
490 * &allints (<return> 1 if all sampled values are ints; else 0)
491 * Return: 0 if OK, 1 on error
492 *
493 * Notes:
494 * (1) Set @maxsamples == 0 to check every integer in na. Otherwise,
495 * this samples no more than @maxsamples.
496 */
497 l_int32
numaHasOnlyIntegers(NUMA * na,l_int32 maxsamples,l_int32 * pallints)498 numaHasOnlyIntegers(NUMA *na,
499 l_int32 maxsamples,
500 l_int32 *pallints)
501 {
502 l_int32 i, n, incr;
503 l_float32 val;
504
505 PROCNAME("numaHasOnlyIntegers");
506
507 if (!pallints)
508 return ERROR_INT("&allints not defined", procName, 1);
509 *pallints = TRUE;
510 if (!na)
511 return ERROR_INT("na not defined", procName, 1);
512
513 if ((n = numaGetCount(na)) == 0)
514 return ERROR_INT("na empty", procName, 1);
515 if (maxsamples <= 0)
516 incr = 1;
517 else
518 incr = (l_int32)((n + maxsamples - 1) / maxsamples);
519 for (i = 0; i < n; i += incr) {
520 numaGetFValue(na, i, &val);
521 if (val != (l_int32)val) {
522 *pallints = FALSE;
523 return 0;
524 }
525 }
526
527 return 0;
528 }
529
530
531 /*!
532 * numaSubsample()
533 *
534 * Input: nas
535 * subfactor (subsample factor, >= 1)
536 * Return: nad (evenly sampled values from nas), or null on error
537 */
538 NUMA *
numaSubsample(NUMA * nas,l_int32 subfactor)539 numaSubsample(NUMA *nas,
540 l_int32 subfactor)
541 {
542 l_int32 i, n;
543 l_float32 val;
544 NUMA *nad;
545
546 PROCNAME("numaSubsample");
547
548 if (!nas)
549 return (NUMA *)ERROR_PTR("nas not defined", procName, NULL);
550 if (subfactor < 1)
551 return (NUMA *)ERROR_PTR("subfactor < 1", procName, NULL);
552
553 nad = numaCreate(0);
554 n = numaGetCount(nas);
555 for (i = 0; i < n; i++) {
556 if (i % subfactor != 0) continue;
557 numaGetFValue(nas, i, &val);
558 numaAddNumber(nad, val);
559 }
560
561 return nad;
562 }
563
564
565 /*!
566 * numaMakeSequence()
567 *
568 * Input: startval
569 * increment
570 * size (of sequence)
571 * Return: numa of sequence of evenly spaced values, or null on error
572 */
573 NUMA *
numaMakeSequence(l_float32 startval,l_float32 increment,l_int32 size)574 numaMakeSequence(l_float32 startval,
575 l_float32 increment,
576 l_int32 size)
577 {
578 l_int32 i;
579 l_float32 val;
580 NUMA *na;
581
582 PROCNAME("numaMakeSequence");
583
584 if ((na = numaCreate(size)) == NULL)
585 return (NUMA *)ERROR_PTR("na not made", procName, NULL);
586
587 for (i = 0; i < size; i++) {
588 val = startval + i * increment;
589 numaAddNumber(na, val);
590 }
591
592 return na;
593 }
594
595
596 /*!
597 * numaMakeConstant()
598 *
599 * Input: val
600 * size (of numa)
601 * Return: numa of given size with all entries equal to 'val',
602 * or null on error
603 */
604 NUMA *
numaMakeConstant(l_float32 val,l_int32 size)605 numaMakeConstant(l_float32 val,
606 l_int32 size)
607 {
608 return numaMakeSequence(val, 0.0, size);
609 }
610
611
612 /*!
613 * numaGetNonzeroRange()
614 *
615 * Input: numa
616 * eps (largest value considered to be zero)
617 * &first, &last (<return> interval of array indices
618 * where values are nonzero)
619 * Return: 0 if OK, 1 on error or if no nonzero range is found.
620 */
621 l_int32
numaGetNonzeroRange(NUMA * na,l_float32 eps,l_int32 * pfirst,l_int32 * plast)622 numaGetNonzeroRange(NUMA *na,
623 l_float32 eps,
624 l_int32 *pfirst,
625 l_int32 *plast)
626 {
627 l_int32 n, i, found;
628 l_float32 val;
629
630 PROCNAME("numaGetNonzeroRange");
631
632 if (!na)
633 return ERROR_INT("na not defined", procName, 1);
634 if (!pfirst || !plast)
635 return ERROR_INT("pfirst and plast not both defined", procName, 1);
636 n = numaGetCount(na);
637 found = FALSE;
638 for (i = 0; i < n; i++) {
639 numaGetFValue(na, i, &val);
640 if (val > eps) {
641 found = TRUE;
642 break;
643 }
644 }
645 if (!found) {
646 *pfirst = n - 1;
647 *plast = 0;
648 return 1;
649 }
650
651 *pfirst = i;
652 for (i = n - 1; i >= 0; i--) {
653 numaGetFValue(na, i, &val);
654 if (val > eps)
655 break;
656 }
657 *plast = i;
658 return 0;
659 }
660
661
662 /*!
663 * numaGetCountRelativeToZero()
664 *
665 * Input: numa
666 * type (L_LESS_THAN_ZERO, L_EQUAL_TO_ZERO, L_GREATER_THAN_ZERO)
667 * &count (<return> count of values of given type)
668 * Return: 0 if OK, 1 on error
669 */
670 l_int32
numaGetCountRelativeToZero(NUMA * na,l_int32 type,l_int32 * pcount)671 numaGetCountRelativeToZero(NUMA *na,
672 l_int32 type,
673 l_int32 *pcount)
674 {
675 l_int32 n, i, count;
676 l_float32 val;
677
678 PROCNAME("numaGetCountRelativeToZero");
679
680 if (!pcount)
681 return ERROR_INT("&count not defined", procName, 1);
682 *pcount = 0;
683 if (!na)
684 return ERROR_INT("na not defined", procName, 1);
685 n = numaGetCount(na);
686 for (i = 0, count = 0; i < n; i++) {
687 numaGetFValue(na, i, &val);
688 if (type == L_LESS_THAN_ZERO && val < 0.0)
689 count++;
690 else if (type == L_EQUAL_TO_ZERO && val == 0.0)
691 count++;
692 else if (type == L_GREATER_THAN_ZERO && val > 0.0)
693 count++;
694 }
695
696 *pcount = count;
697 return 0;
698 }
699
700
701 /*!
702 * numaClipToInterval()
703 *
704 * Input: numa
705 * first, last (clipping interval)
706 * Return: numa with the same values as the input, but clipped
707 * to the specified interval
708 *
709 * Note: If you want the indices of the array values to be unchanged,
710 * use first = 0.
711 * Usage: This is useful to clip a histogram that has a few nonzero
712 * values to its nonzero range.
713 */
714 NUMA *
numaClipToInterval(NUMA * nas,l_int32 first,l_int32 last)715 numaClipToInterval(NUMA *nas,
716 l_int32 first,
717 l_int32 last)
718 {
719 l_int32 n, i, truelast;
720 l_float32 val;
721 NUMA *nad;
722
723 PROCNAME("numaClipToInterval");
724
725 if (!nas)
726 return (NUMA *)ERROR_PTR("nas not defined", procName, NULL);
727 if (first > last)
728 return (NUMA *)ERROR_PTR("range not valid", procName, NULL);
729
730 n = numaGetCount(nas);
731 if (first >= n)
732 return (NUMA *)ERROR_PTR("no elements in range", procName, NULL);
733 truelast = L_MIN(last, n - 1);
734 if ((nad = numaCreate(truelast - first + 1)) == NULL)
735 return (NUMA *)ERROR_PTR("nad not made", procName, NULL);
736 for (i = first; i <= truelast; i++) {
737 numaGetFValue(nas, i, &val);
738 numaAddNumber(nad, val);
739 }
740
741 return nad;
742 }
743
744
745 /*!
746 * numaMakeThresholdIndicator()
747 *
748 * Input: nas (input numa)
749 * thresh (threshold value)
750 * type (L_SELECT_IF_LT, L_SELECT_IF_GT,
751 * L_SELECT_IF_LTE, L_SELECT_IF_GTE)
752 * Output: nad (indicator array: values are 0 and 1)
753 *
754 * Notes:
755 * (1) For each element in nas, if the constraint given by 'type'
756 * correctly specifies its relation to thresh, a value of 1
757 * is recorded in nad.
758 */
759 NUMA *
numaMakeThresholdIndicator(NUMA * nas,l_float32 thresh,l_int32 type)760 numaMakeThresholdIndicator(NUMA *nas,
761 l_float32 thresh,
762 l_int32 type)
763 {
764 l_int32 n, i, ival;
765 l_float32 fval;
766 NUMA *nai;
767
768 PROCNAME("numaMakeThresholdIndicator");
769
770 if (!nas)
771 return (NUMA *)ERROR_PTR("nas not defined", procName, NULL);
772 n = numaGetCount(nas);
773 nai = numaCreate(n);
774 for (i = 0; i < n; i++) {
775 numaGetFValue(nas, i, &fval);
776 ival = 0;
777 switch (type)
778 {
779 case L_SELECT_IF_LT:
780 if (fval < thresh) ival = 1;
781 break;
782 case L_SELECT_IF_GT:
783 if (fval > thresh) ival = 1;
784 break;
785 case L_SELECT_IF_LTE:
786 if (fval <= thresh) ival = 1;
787 break;
788 case L_SELECT_IF_GTE:
789 if (fval >= thresh) ival = 1;
790 break;
791 default:
792 numaDestroy(&nai);
793 return (NUMA *)ERROR_PTR("invalid type", procName, NULL);
794 }
795 numaAddNumber(nai, ival);
796 }
797
798 return nai;
799 }
800
801
802 /*!
803 * numaUniformSampling()
804 *
805 * Input: nas (input numa)
806 * nsamp (number of samples)
807 * Output: nad (resampled array), or null on error
808 *
809 * Notes:
810 * (1) This resamples the values in the array, using @nsamp
811 * equal divisions.
812 */
813 NUMA *
numaUniformSampling(NUMA * nas,l_int32 nsamp)814 numaUniformSampling(NUMA *nas,
815 l_int32 nsamp)
816 {
817 l_int32 n, i, j, ileft, iright;
818 l_float32 left, right, binsize, lfract, rfract, sum, startx, delx;
819 l_float32 *array;
820 NUMA *nad;
821
822 PROCNAME("numaUniformSampling");
823
824 if (!nas)
825 return (NUMA *)ERROR_PTR("nas not defined", procName, NULL);
826 if (nsamp <= 0)
827 return (NUMA *)ERROR_PTR("nsamp must be > 0", procName, NULL);
828
829 n = numaGetCount(nas);
830 nad = numaCreate(nsamp);
831 array = numaGetFArray(nas, L_NOCOPY);
832 binsize = (l_float32)n / (l_float32)nsamp;
833 numaGetXParameters(nas, &startx, &delx);
834 numaSetXParameters(nad, startx, binsize * delx);
835 left = 0.0;
836 for (i = 0; i < nsamp; i++) {
837 sum = 0.0;
838 right = left + binsize;
839 ileft = (l_int32)left;
840 lfract = 1.0 - left + ileft;
841 if (lfract >= 1.0) /* on left bin boundary */
842 lfract = 0.0;
843 iright = (l_int32)right;
844 rfract = right - iright;
845 iright = L_MIN(iright, n - 1);
846 if (ileft == iright) { /* both are within the same original sample */
847 sum += (lfract + rfract - 1.0) * array[ileft];
848 }
849 else {
850 if (lfract > 0.0001) /* left fraction */
851 sum += lfract * array[ileft];
852 if (rfract > 0.0001) /* right fraction */
853 sum += rfract * array[iright];
854 for (j = ileft + 1; j < iright; j++) /* entire pixels */
855 sum += array[j];
856 }
857
858 numaAddNumber(nad, sum);
859 left = right;
860 }
861 return nad;
862 }
863
864
865
866 /*----------------------------------------------------------------------*
867 * Interpolation *
868 *----------------------------------------------------------------------*/
869 /*!
870 * numaInterpolateEqxVal()
871 *
872 * Input: startx (xval corresponding to first element in array)
873 * deltax (x increment between array elements)
874 * nay (numa of ordinate values, assumed equally spaced)
875 * type (L_LINEAR_INTERP, L_QUADRATIC_INTERP)
876 * xval
877 * &yval (<return> interpolated value)
878 * Return: 0 if OK, 1 on error (e.g., if xval is outside range)
879 *
880 * Notes:
881 * (1) Considering nay as a function of x, the x values
882 * are equally spaced
883 * (2) Caller should check for valid return.
884 *
885 * For linear Lagrangian interpolation (through 2 data pts):
886 * y(x) = y1(x-x2)/(x1-x2) + y2(x-x1)/(x2-x1)
887 *
888 * For quadratic Lagrangian interpolation (through 3 data pts):
889 * y(x) = y1(x-x2)(x-x3)/((x1-x2)(x1-x3)) +
890 * y2(x-x1)(x-x3)/((x2-x1)(x2-x3)) +
891 * y3(x-x1)(x-x2)/((x3-x1)(x3-x2))
892 *
893 */
894 l_int32
numaInterpolateEqxVal(l_float32 startx,l_float32 deltax,NUMA * nay,l_int32 type,l_float32 xval,l_float32 * pyval)895 numaInterpolateEqxVal(l_float32 startx,
896 l_float32 deltax,
897 NUMA *nay,
898 l_int32 type,
899 l_float32 xval,
900 l_float32 *pyval)
901 {
902 l_int32 i, n, i1, i2, i3;
903 l_float32 x1, x2, x3, fy1, fy2, fy3, d1, d2, d3, del, fi, maxx;
904 l_float32 *fa;
905
906 PROCNAME("numaInterpolateEqxVal");
907
908 if (!pyval)
909 return ERROR_INT("&yval not defined", procName, 1);
910 *pyval = 0.0;
911 if (!nay)
912 return ERROR_INT("nay not defined", procName, 1);
913 if (deltax <= 0.0)
914 return ERROR_INT("deltax not > 0", procName, 1);
915 if (type != L_LINEAR_INTERP && type != L_QUADRATIC_INTERP)
916 return ERROR_INT("invalid interp type", procName, 1);
917 n = numaGetCount(nay);
918 if (n < 2)
919 return ERROR_INT("not enough points", procName, 1);
920 if (type == L_QUADRATIC_INTERP && n == 2) {
921 type = L_LINEAR_INTERP;
922 L_WARNING("only 2 points; using linear interp", procName);
923 }
924 maxx = startx + deltax * (n - 1);
925 if (xval < startx || xval > maxx)
926 return ERROR_INT("xval is out of bounds", procName, 1);
927
928 fa = numaGetFArray(nay, L_NOCOPY);
929 fi = (xval - startx) / deltax;
930 i = (l_int32)fi;
931 del = fi - i;
932 if (del == 0.0) { /* no interpolation required */
933 *pyval = fa[i];
934 return 0;
935 }
936
937 if (type == L_LINEAR_INTERP) {
938 *pyval = fa[i] + del * (fa[i + 1] - fa[i]);
939 return 0;
940 }
941
942 /* Quadratic interpolation */
943 d1 = d3 = 0.5 / (deltax * deltax);
944 d2 = -2. * d1;
945 if (i == 0) {
946 i1 = i;
947 i2 = i + 1;
948 i3 = i + 2;
949 }
950 else {
951 i1 = i - 1;
952 i2 = i;
953 i3 = i + 1;
954 }
955 x1 = startx + i1 * deltax;
956 x2 = startx + i2 * deltax;
957 x3 = startx + i3 * deltax;
958 fy1 = d1 * fa[i1];
959 fy2 = d2 * fa[i2];
960 fy3 = d3 * fa[i3];
961 *pyval = fy1 * (xval - x2) * (xval - x3) +
962 fy2 * (xval - x1) * (xval - x3) +
963 fy3 * (xval - x1) * (xval - x2);
964 return 0;
965 }
966
967
968 /*!
969 * numaInterpolateArbxVal()
970 *
971 * Input: nax (numa of abscissa values)
972 * nay (numa of ordinate values, corresponding to nax)
973 * type (L_LINEAR_INTERP, L_QUADRATIC_INTERP)
974 * xval
975 * &yval (<return> interpolated value)
976 * Return: 0 if OK, 1 on error (e.g., if xval is outside range)
977 *
978 * Notes:
979 * (1) The values in nax must be sorted in increasing order.
980 * If, additionally, they are equally spaced, you can use
981 * numaInterpolateEqxVal().
982 * (2) Caller should check for valid return.
983 * (3) Uses lagrangian interpolation. See numaInterpolateEqxVal()
984 * for formulas.
985 */
986 l_int32
numaInterpolateArbxVal(NUMA * nax,NUMA * nay,l_int32 type,l_float32 xval,l_float32 * pyval)987 numaInterpolateArbxVal(NUMA *nax,
988 NUMA *nay,
989 l_int32 type,
990 l_float32 xval,
991 l_float32 *pyval)
992 {
993 l_int32 i, im, nx, ny, i1, i2, i3;
994 l_float32 delu, dell, fract, d1, d2, d3;
995 l_float32 minx, maxx;
996 l_float32 *fax, *fay;
997
998 PROCNAME("numaInterpolateArbxVal");
999
1000 if (!pyval)
1001 return ERROR_INT("&yval not defined", procName, 1);
1002 *pyval = 0.0;
1003 if (!nax)
1004 return ERROR_INT("nax not defined", procName, 1);
1005 if (!nay)
1006 return ERROR_INT("nay not defined", procName, 1);
1007 if (type != L_LINEAR_INTERP && type != L_QUADRATIC_INTERP)
1008 return ERROR_INT("invalid interp type", procName, 1);
1009 ny = numaGetCount(nay);
1010 nx = numaGetCount(nax);
1011 if (nx != ny)
1012 return ERROR_INT("nax and nay not same size arrays", procName, 1);
1013 if (ny < 2)
1014 return ERROR_INT("not enough points", procName, 1);
1015 if (type == L_QUADRATIC_INTERP && ny == 2) {
1016 type = L_LINEAR_INTERP;
1017 L_WARNING("only 2 points; using linear interp", procName);
1018 }
1019 numaGetFValue(nax, 0, &minx);
1020 numaGetFValue(nax, nx - 1, &maxx);
1021 if (xval < minx || xval > maxx)
1022 return ERROR_INT("xval is out of bounds", procName, 1);
1023
1024 fax = numaGetFArray(nax, L_NOCOPY);
1025 fay = numaGetFArray(nay, L_NOCOPY);
1026
1027 /* Linear search for interval. We are guaranteed
1028 * to either return or break out of the loop.
1029 * In addition, we are assured that fax[i] - fax[im] > 0.0 */
1030 if (xval == fax[0]) {
1031 *pyval = fay[0];
1032 return 0;
1033 }
1034 for (i = 1; i < nx; i++) {
1035 delu = fax[i] - xval;
1036 if (delu >= 0.0) { /* we've passed it */
1037 if (delu == 0.0) {
1038 *pyval = fay[i];
1039 return 0;
1040 }
1041 im = i - 1;
1042 dell = xval - fax[im]; /* >= 0 */
1043 break;
1044 }
1045 }
1046 fract = dell / (fax[i] - fax[im]);
1047
1048 if (type == L_LINEAR_INTERP) {
1049 *pyval = fay[i] + fract * (fay[i + 1] - fay[i]);
1050 return 0;
1051 }
1052
1053 /* Quadratic interpolation */
1054 if (im == 0) {
1055 i1 = im;
1056 i2 = im + 1;
1057 i3 = im + 2;
1058 }
1059 else {
1060 i1 = im - 1;
1061 i2 = im;
1062 i3 = im + 1;
1063 }
1064 d1 = (fax[i1] - fax[i2]) * (fax[i1] - fax[i3]);
1065 d2 = (fax[i2] - fax[i1]) * (fax[i2] - fax[i3]);
1066 d3 = (fax[i3] - fax[i1]) * (fax[i3] - fax[i2]);
1067 *pyval = fay[i1] * (xval - fax[i2]) * (xval - fax[i3]) / d1 +
1068 fay[i2] * (xval - fax[i1]) * (xval - fax[i3]) / d2 +
1069 fay[i3] * (xval - fax[i1]) * (xval - fax[i2]) / d3;
1070 return 0;
1071 }
1072
1073
1074 /*!
1075 * numaInterpolateEqxInterval()
1076 *
1077 * Input: startx (xval corresponding to first element in nas)
1078 * deltax (x increment between array elements in nas)
1079 * nasy (numa of ordinate values, assumed equally spaced)
1080 * type (L_LINEAR_INTERP, L_QUADRATIC_INTERP)
1081 * x0 (start value of interval)
1082 * x1 (end value of interval)
1083 * npts (number of points to evaluate function in interval)
1084 * &nax (<optional return> array of x values in interval)
1085 * &nay (<return> array of y values in interval)
1086 * Return: 0 if OK, 1 on error
1087 *
1088 * Notes:
1089 * (1) Considering nasy as a function of x, the x values
1090 * are equally spaced.
1091 * (2) This creates nay (and optionally nax) of interpolated
1092 * values over the specified interval (x0, x1).
1093 * (3) If the interval (x0, x1) lies partially outside the array
1094 * nasy (as interpreted by startx and deltax), it is an
1095 * error and returns 1.
1096 * (4) Note that deltax is the intrinsic x-increment for the input
1097 * array nasy, whereas delx is the intrinsic x-increment for the
1098 * output interpolated array nay.
1099 */
1100 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)1101 numaInterpolateEqxInterval(l_float32 startx,
1102 l_float32 deltax,
1103 NUMA *nasy,
1104 l_int32 type,
1105 l_float32 x0,
1106 l_float32 x1,
1107 l_int32 npts,
1108 NUMA **pnax,
1109 NUMA **pnay)
1110 {
1111 l_int32 i, n;
1112 l_float32 x, yval, maxx, delx;
1113 NUMA *nax, *nay;
1114
1115 PROCNAME("numaInterpolateEqxInterval");
1116
1117 if (pnax) *pnax = NULL;
1118 if (!pnay)
1119 return ERROR_INT("&nay not defined", procName, 1);
1120 *pnay = NULL;
1121 if (!nasy)
1122 return ERROR_INT("nasy not defined", procName, 1);
1123 if (deltax <= 0.0)
1124 return ERROR_INT("deltax not > 0", procName, 1);
1125 if (type != L_LINEAR_INTERP && type != L_QUADRATIC_INTERP)
1126 return ERROR_INT("invalid interp type", procName, 1);
1127 n = numaGetCount(nasy);
1128 if (type == L_QUADRATIC_INTERP && n == 2) {
1129 type = L_LINEAR_INTERP;
1130 L_WARNING("only 2 points; using linear interp", procName);
1131 }
1132 maxx = startx + deltax * (n - 1);
1133 if (x0 < startx || x1 > maxx || x1 <= x0)
1134 return ERROR_INT("[x0 ... x1] is not valid", procName, 1);
1135 if (npts < 3)
1136 return ERROR_INT("npts < 3", procName, 1);
1137 delx = (x1 - x0) / (l_float32)(npts - 1); /* delx is for output nay */
1138
1139 if ((nay = numaCreate(npts)) == NULL)
1140 return ERROR_INT("nay not made", procName, 1);
1141 numaSetXParameters(nay, x0, delx);
1142 *pnay = nay;
1143 if (pnax) {
1144 nax = numaCreate(npts);
1145 *pnax = nax;
1146 }
1147
1148 for (i = 0; i < npts; i++) {
1149 x = x0 + i * delx;
1150 if (pnax)
1151 numaAddNumber(nax, x);
1152 numaInterpolateEqxVal(startx, deltax, nasy, type, x, &yval);
1153 numaAddNumber(nay, yval);
1154 }
1155
1156 return 0;
1157 }
1158
1159
1160 /*!
1161 * numaInterpolateArbxInterval()
1162 *
1163 * Input: nax (numa of abscissa values)
1164 * nay (numa of ordinate values, corresponding to nax)
1165 * type (L_LINEAR_INTERP, L_QUADRATIC_INTERP)
1166 * x0 (start value of interval)
1167 * x1 (end value of interval)
1168 * npts (number of points to evaluate function in interval)
1169 * &nadx (<optional return> array of x values in interval)
1170 * &nady (<return> array of y values in interval)
1171 * Return: 0 if OK, 1 on error (e.g., if x0 or x1 is outside range)
1172 *
1173 * Notes:
1174 * (1) The values in nax must be sorted in increasing order.
1175 * If they are not sorted, we do it here, and complain.
1176 * (2) If the values in nax are equally spaced, you can use
1177 * numaInterpolateEqxInterval().
1178 * (3) Caller should check for valid return.
1179 * (4) We don't call numaInterpolateArbxVal() for each output
1180 * point, because that requires an O(n) search for
1181 * each point. Instead, we do a single O(n) pass through
1182 * nax, saving the indices to be used for each output yval.
1183 * (5) Uses lagrangian interpolation. See numaInterpolateEqxVal()
1184 * for formulas.
1185 */
1186 l_int32
numaInterpolateArbxInterval(NUMA * nax,NUMA * nay,l_int32 type,l_float32 x0,l_float32 x1,l_int32 npts,NUMA ** pnadx,NUMA ** pnady)1187 numaInterpolateArbxInterval(NUMA *nax,
1188 NUMA *nay,
1189 l_int32 type,
1190 l_float32 x0,
1191 l_float32 x1,
1192 l_int32 npts,
1193 NUMA **pnadx,
1194 NUMA **pnady)
1195 {
1196 l_int32 i, im, j, nx, ny, i1, i2, i3, sorted;
1197 l_int32 *index;
1198 l_float32 del, xval, yval, excess, fract, minx, maxx, d1, d2, d3;
1199 l_float32 *fax, *fay;
1200 NUMA *nasx, *nasy, *nadx, *nady;
1201
1202 PROCNAME("numaInterpolateArbxInterval");
1203
1204 if (pnadx) *pnadx = NULL;
1205 if (!pnady)
1206 return ERROR_INT("&nady not defined", procName, 1);
1207 *pnady = NULL;
1208 if (!nay)
1209 return ERROR_INT("nay not defined", procName, 1);
1210 if (!nax)
1211 return ERROR_INT("nax not defined", procName, 1);
1212 if (type != L_LINEAR_INTERP && type != L_QUADRATIC_INTERP)
1213 return ERROR_INT("invalid interp type", procName, 1);
1214 if (x0 > x1)
1215 return ERROR_INT("x0 > x1", procName, 1);
1216 ny = numaGetCount(nay);
1217 nx = numaGetCount(nax);
1218 if (nx != ny)
1219 return ERROR_INT("nax and nay not same size arrays", procName, 1);
1220 if (ny < 2)
1221 return ERROR_INT("not enough points", procName, 1);
1222 if (type == L_QUADRATIC_INTERP && ny == 2) {
1223 type = L_LINEAR_INTERP;
1224 L_WARNING("only 2 points; using linear interp", procName);
1225 }
1226 numaGetMin(nax, &minx, NULL);
1227 numaGetMax(nax, &maxx, NULL);
1228 if (x0 < minx || x1 > maxx)
1229 return ERROR_INT("xval is out of bounds", procName, 1);
1230
1231 /* Make sure that nax is sorted in increasing order */
1232 numaIsSorted(nax, L_SORT_INCREASING, &sorted);
1233 if (!sorted) {
1234 L_WARNING("we are sorting nax in increasing order", procName);
1235 numaSortPair(nax, nay, L_SORT_INCREASING, &nasx, &nasy);
1236 }
1237 else {
1238 nasx = numaClone(nax);
1239 nasy = numaClone(nay);
1240 }
1241
1242 fax = numaGetFArray(nasx, L_NOCOPY);
1243 fay = numaGetFArray(nasy, L_NOCOPY);
1244
1245 /* Get array of indices into fax for interpolated locations */
1246 if ((index = (l_int32 *)CALLOC(npts, sizeof(l_int32))) == NULL)
1247 return ERROR_INT("ind not made", procName, 1);
1248 del = (x1 - x0) / (npts - 1.0);
1249 for (i = 0, j = 0; j < nx && i < npts; i++) {
1250 xval = x0 + i * del;
1251 while (j < nx - 1 && xval > fax[j])
1252 j++;
1253 if (xval == fax[j])
1254 index[i] = L_MIN(j, nx - 1);
1255 else /* the index of fax[] is just below xval */
1256 index[i] = L_MAX(j - 1, 0);
1257 }
1258
1259 /* For each point to be interpolated, get the y value */
1260 nady = numaCreate(npts);
1261 *pnady = nady;
1262 if (pnadx) {
1263 nadx = numaCreate(npts);
1264 *pnadx = nadx;
1265 }
1266 for (i = 0; i < npts; i++) {
1267 xval = x0 + i * del;
1268 if (pnadx)
1269 numaAddNumber(nadx, xval);
1270 im = index[i];
1271 excess = xval - fax[im];
1272 if (excess == 0.0) {
1273 numaAddNumber(nady, fay[im]);
1274 continue;
1275 }
1276 fract = excess / (fax[im + 1] - fax[im]);
1277
1278 if (type == L_LINEAR_INTERP) {
1279 yval = fay[im] + fract * (fay[im + 1] - fay[im]);
1280 numaAddNumber(nady, yval);
1281 continue;
1282 }
1283
1284 /* Quadratic interpolation */
1285 if (im == 0) {
1286 i1 = im;
1287 i2 = im + 1;
1288 i3 = im + 2;
1289 }
1290 else {
1291 i1 = im - 1;
1292 i2 = im;
1293 i3 = im + 1;
1294 }
1295 d1 = (fax[i1] - fax[i2]) * (fax[i1] - fax[i3]);
1296 d2 = (fax[i2] - fax[i1]) * (fax[i2] - fax[i3]);
1297 d3 = (fax[i3] - fax[i1]) * (fax[i3] - fax[i2]);
1298 yval = fay[i1] * (xval - fax[i2]) * (xval - fax[i3]) / d1 +
1299 fay[i2] * (xval - fax[i1]) * (xval - fax[i3]) / d2 +
1300 fay[i3] * (xval - fax[i1]) * (xval - fax[i2]) / d3;
1301 numaAddNumber(nady, yval);
1302 }
1303
1304 FREE(index);
1305 numaDestroy(&nasx);
1306 numaDestroy(&nasy);
1307 return 0;
1308 }
1309
1310
1311 /*----------------------------------------------------------------------*
1312 * Functions requiring interpolation *
1313 *----------------------------------------------------------------------*/
1314 /*!
1315 * numaFitMax()
1316 *
1317 * Input: na (numa of ordinate values, to fit a max to)
1318 * &maxval (<return> max value)
1319 * naloc (<optional> associated numa of abscissa values)
1320 * &maxloc (<return> abscissa value that gives max value in na;
1321 * if naloc == null, this is given as an interpolated
1322 * index value)
1323 * Return: 0 if OK; 1 on error
1324 *
1325 * Note: if naloc is given, there is no requirement that the
1326 * data points are evenly spaced. Lagrangian interpolation
1327 * handles that. The only requirement is that the
1328 * data points are ordered so that the values in naloc
1329 * are either increasing or decreasing. We test to make
1330 * sure that the sizes of na and naloc are equal, and it
1331 * is assumed that the correspondences na[i] as a function
1332 * of naloc[i] are properly arranged for all i.
1333 *
1334 * The formula for Lagrangian interpolation through 3 data pts is:
1335 * y(x) = y1(x-x2)(x-x3)/((x1-x2)(x1-x3)) +
1336 * y2(x-x1)(x-x3)/((x2-x1)(x2-x3)) +
1337 * y3(x-x1)(x-x2)/((x3-x1)(x3-x2))
1338 *
1339 * Then the derivative, using the constants (c1,c2,c3) defined below,
1340 * is set to 0:
1341 * y'(x) = 2x(c1+c2+c3) - c1(x2+x3) - c2(x1+x3) - c3(x1+x2) = 0
1342 */
1343 l_int32
numaFitMax(NUMA * na,l_float32 * pmaxval,NUMA * naloc,l_float32 * pmaxloc)1344 numaFitMax(NUMA *na,
1345 l_float32 *pmaxval,
1346 NUMA *naloc,
1347 l_float32 *pmaxloc)
1348 {
1349 l_float32 val;
1350 l_float32 smaxval; /* start value of maximum sample, before interpolating */
1351 l_int32 n, imaxloc;
1352 l_float32 x1, x2, x3, y1, y2, y3, c1, c2, c3, a, b, xmax, ymax;
1353
1354 PROCNAME("numaFitMax");
1355
1356 *pmaxval = *pmaxloc = 0.0; /* init */
1357
1358 if (!na)
1359 return ERROR_INT("na not defined", procName, 1);
1360 if (!pmaxval)
1361 return ERROR_INT("&maxval not defined", procName, 1);
1362 if (!pmaxloc)
1363 return ERROR_INT("&maxloc not defined", procName, 1);
1364
1365 n = numaGetCount(na);
1366 if (naloc) {
1367 if (n != numaGetCount(naloc))
1368 return ERROR_INT("na and naloc of unequal size", procName, 1);
1369 }
1370
1371 numaGetMax(na, &smaxval, &imaxloc);
1372
1373 /* Simple case: max is at end point */
1374 if (imaxloc == 0 || imaxloc == n - 1) {
1375 *pmaxval = smaxval;
1376 if (naloc) {
1377 numaGetFValue(naloc, imaxloc, &val);
1378 *pmaxloc = val;
1379 }
1380 else
1381 *pmaxloc = imaxloc;
1382 return 0;
1383 }
1384
1385 /* Interior point; use quadratic interpolation */
1386 y2 = smaxval;
1387 numaGetFValue(na, imaxloc - 1, &val);
1388 y1 = val;
1389 numaGetFValue(na, imaxloc + 1, &val);
1390 y3 = val;
1391 if (naloc) {
1392 numaGetFValue(naloc, imaxloc - 1, &val);
1393 x1 = val;
1394 numaGetFValue(naloc, imaxloc, &val);
1395 x2 = val;
1396 numaGetFValue(naloc, imaxloc + 1, &val);
1397 x3 = val;
1398 }
1399 else {
1400 x1 = imaxloc - 1;
1401 x2 = imaxloc;
1402 x3 = imaxloc + 1;
1403 }
1404
1405 /* Can't interpolate; just use the max val in na
1406 * and the corresponding one in naloc */
1407 if (x1 == x2 || x1 == x3 || x2 == x3) {
1408 *pmaxval = y2;
1409 *pmaxloc = x2;
1410 return 0;
1411 }
1412
1413 /* Use lagrangian interpolation; set dy/dx = 0 */
1414 c1 = y1 / ((x1 - x2) * (x1 - x3));
1415 c2 = y2 / ((x2 - x1) * (x2 - x3));
1416 c3 = y3 / ((x3 - x1) * (x3 - x2));
1417 a = c1 + c2 + c3;
1418 b = c1 * (x2 + x3) + c2 * (x1 + x3) + c3 * (x1 + x2);
1419 xmax = b / (2 * a);
1420 ymax = c1 * (xmax - x2) * (xmax - x3) +
1421 c2 * (xmax - x1) * (xmax - x3) +
1422 c3 * (xmax - x1) * (xmax - x2);
1423 *pmaxval = ymax;
1424 *pmaxloc = xmax;
1425
1426 return 0;
1427 }
1428
1429
1430 /*!
1431 * numaDifferentiateInterval()
1432 *
1433 * Input: nax (numa of abscissa values)
1434 * nay (numa of ordinate values, corresponding to nax)
1435 * x0 (start value of interval)
1436 * x1 (end value of interval)
1437 * npts (number of points to evaluate function in interval)
1438 * &nadx (<optional return> array of x values in interval)
1439 * &nady (<return> array of derivatives in interval)
1440 * Return: 0 if OK, 1 on error (e.g., if x0 or x1 is outside range)
1441 *
1442 * Notes:
1443 * (1) The values in nax must be sorted in increasing order.
1444 * If they are not sorted, it is done in the interpolation
1445 * step, and a warning is issued.
1446 * (2) Caller should check for valid return.
1447 */
1448 l_int32
numaDifferentiateInterval(NUMA * nax,NUMA * nay,l_float32 x0,l_float32 x1,l_int32 npts,NUMA ** pnadx,NUMA ** pnady)1449 numaDifferentiateInterval(NUMA *nax,
1450 NUMA *nay,
1451 l_float32 x0,
1452 l_float32 x1,
1453 l_int32 npts,
1454 NUMA **pnadx,
1455 NUMA **pnady)
1456 {
1457 l_int32 i, nx, ny;
1458 l_float32 minx, maxx, der, invdel;
1459 l_float32 *fay;
1460 NUMA *nady, *naiy;
1461
1462 PROCNAME("numaDifferentiateInterval");
1463
1464 if (pnadx) *pnadx = NULL;
1465 if (!pnady)
1466 return ERROR_INT("&nady not defined", procName, 1);
1467 *pnady = NULL;
1468 if (!nay)
1469 return ERROR_INT("nay not defined", procName, 1);
1470 if (!nax)
1471 return ERROR_INT("nax not defined", procName, 1);
1472 if (x0 > x1)
1473 return ERROR_INT("x0 > x1", procName, 1);
1474 ny = numaGetCount(nay);
1475 nx = numaGetCount(nax);
1476 if (nx != ny)
1477 return ERROR_INT("nax and nay not same size arrays", procName, 1);
1478 if (ny < 2)
1479 return ERROR_INT("not enough points", procName, 1);
1480 numaGetMin(nax, &minx, NULL);
1481 numaGetMax(nax, &maxx, NULL);
1482 if (x0 < minx || x1 > maxx)
1483 return ERROR_INT("xval is out of bounds", procName, 1);
1484 if (npts < 2)
1485 return ERROR_INT("npts < 2", procName, 1);
1486
1487 /* Generate interpolated array over specified interval */
1488 if (numaInterpolateArbxInterval(nax, nay, L_LINEAR_INTERP, x0, x1,
1489 npts, pnadx, &naiy))
1490 return ERROR_INT("interpolation failed", procName, 1);
1491
1492 nady = numaCreate(npts);
1493 *pnady = nady;
1494 invdel = 0.5 * ((l_float32)npts - 1.0) / (x1 - x0);
1495 fay = numaGetFArray(naiy, L_NOCOPY);
1496
1497 /* Compute and save derivatives */
1498 der = 0.5 * invdel * (fay[1] - fay[0]);
1499 numaAddNumber(nady, der);
1500 for (i = 1; i < npts - 1; i++) {
1501 der = invdel * (fay[i + 1] - fay[i - 1]);
1502 numaAddNumber(nady, der);
1503 }
1504 der = 0.5 * invdel * (fay[npts - 1] - fay[npts - 2]);
1505 numaAddNumber(nady, der);
1506
1507 numaDestroy(&naiy);
1508 return 0;
1509 }
1510
1511
1512 /*!
1513 * numaIntegrateInterval()
1514 *
1515 * Input: nax (numa of abscissa values)
1516 * nay (numa of ordinate values, corresponding to nax)
1517 * x0 (start value of interval)
1518 * x1 (end value of interval)
1519 * npts (number of points to evaluate function in interval)
1520 * &sum (<return> integral of function over interval)
1521 * Return: 0 if OK, 1 on error (e.g., if x0 or x1 is outside range)
1522 *
1523 * Notes:
1524 * (1) The values in nax must be sorted in increasing order.
1525 * If they are not sorted, it is done in the interpolation
1526 * step, and a warning is issued.
1527 * (2) Caller should check for valid return.
1528 */
1529 l_int32
numaIntegrateInterval(NUMA * nax,NUMA * nay,l_float32 x0,l_float32 x1,l_int32 npts,l_float32 * psum)1530 numaIntegrateInterval(NUMA *nax,
1531 NUMA *nay,
1532 l_float32 x0,
1533 l_float32 x1,
1534 l_int32 npts,
1535 l_float32 *psum)
1536 {
1537 l_int32 i, nx, ny;
1538 l_float32 minx, maxx, sum, del;
1539 l_float32 *fay;
1540 NUMA *naiy;
1541
1542 PROCNAME("numaIntegrateInterval");
1543
1544 if (!psum)
1545 return ERROR_INT("&sum not defined", procName, 1);
1546 *psum = 0.0;
1547 if (!nay)
1548 return ERROR_INT("nay not defined", procName, 1);
1549 if (!nax)
1550 return ERROR_INT("nax not defined", procName, 1);
1551 if (x0 > x1)
1552 return ERROR_INT("x0 > x1", procName, 1);
1553 if (npts < 2)
1554 return ERROR_INT("npts < 2", procName, 1);
1555 ny = numaGetCount(nay);
1556 nx = numaGetCount(nax);
1557 if (nx != ny)
1558 return ERROR_INT("nax and nay not same size arrays", procName, 1);
1559 if (ny < 2)
1560 return ERROR_INT("not enough points", procName, 1);
1561 numaGetMin(nax, &minx, NULL);
1562 numaGetMax(nax, &maxx, NULL);
1563 if (x0 < minx || x1 > maxx)
1564 return ERROR_INT("xval is out of bounds", procName, 1);
1565
1566 /* Generate interpolated array over specified interval */
1567 if (numaInterpolateArbxInterval(nax, nay, L_LINEAR_INTERP, x0, x1,
1568 npts, NULL, &naiy))
1569 return ERROR_INT("interpolation failed", procName, 1);
1570
1571 del = (x1 - x0) / ((l_float32)npts - 1.0);
1572 fay = numaGetFArray(naiy, L_NOCOPY);
1573
1574 /* Compute integral (simple trapezoid) */
1575 sum = 0.5 * (fay[0] + fay[npts - 1]);
1576 for (i = 1; i < npts - 1; i++)
1577 sum += fay[i];
1578 *psum = del * sum;
1579
1580 numaDestroy(&naiy);
1581 return 0;
1582 }
1583
1584
1585 /*----------------------------------------------------------------------*
1586 * Sorting *
1587 *----------------------------------------------------------------------*/
1588 /*!
1589 * numaSort()
1590 *
1591 * Input: naout (output numa; can be NULL or equal to nain)
1592 * nain (input numa)
1593 * sortorder (L_SORT_INCREASING or L_SORT_DECREASING)
1594 * Return: naout (output sorted numa), or null on error
1595 *
1596 * Notes:
1597 * (1) Set naout = nain for in-place; otherwise, set naout = NULL.
1598 * (2) Source: Shell sort, modified from K&R, 2nd edition, p.62.
1599 * Slow but simple O(n logn) sort.
1600 */
1601 NUMA *
numaSort(NUMA * naout,NUMA * nain,l_int32 sortorder)1602 numaSort(NUMA *naout,
1603 NUMA *nain,
1604 l_int32 sortorder)
1605 {
1606 l_int32 i, n, gap, j;
1607 l_float32 tmp;
1608 l_float32 *array;
1609
1610 PROCNAME("numaSort");
1611
1612 if (!nain)
1613 return (NUMA *)ERROR_PTR("nain not defined", procName, NULL);
1614
1615 /* Make naout if necessary; otherwise do in-place */
1616 if (!naout)
1617 naout = numaCopy(nain);
1618 else if (nain != naout)
1619 return (NUMA *)ERROR_PTR("invalid: not in-place", procName, NULL);
1620 array = naout->array; /* operate directly on the array */
1621 n = numaGetCount(naout);
1622
1623 /* Shell sort */
1624 for (gap = n/2; gap > 0; gap = gap / 2) {
1625 for (i = gap; i < n; i++) {
1626 for (j = i - gap; j >= 0; j -= gap) {
1627 if ((sortorder == L_SORT_INCREASING &&
1628 array[j] > array[j + gap]) ||
1629 (sortorder == L_SORT_DECREASING &&
1630 array[j] < array[j + gap]))
1631 {
1632 tmp = array[j];
1633 array[j] = array[j + gap];
1634 array[j + gap] = tmp;
1635 }
1636 }
1637 }
1638 }
1639
1640 return naout;
1641 }
1642
1643
1644 /*!
1645 * numaGetSortIndex()
1646 *
1647 * Input: na
1648 * sortorder (L_SORT_INCREASING or L_SORT_DECREASING)
1649 * Return: na giving an array of indices that would sort
1650 * the input array, or null on error
1651 */
1652 NUMA *
numaGetSortIndex(NUMA * na,l_int32 sortorder)1653 numaGetSortIndex(NUMA *na,
1654 l_int32 sortorder)
1655 {
1656 l_int32 i, n, gap, j;
1657 l_float32 tmp;
1658 l_float32 *array; /* copy of input array */
1659 l_float32 *iarray; /* array of indices */
1660 NUMA *naisort;
1661
1662 PROCNAME("numaGetSortIndex");
1663
1664 if (!na)
1665 return (NUMA *)ERROR_PTR("na not defined", procName, NULL);
1666 if (sortorder != L_SORT_INCREASING && sortorder != L_SORT_DECREASING)
1667 return (NUMA *)ERROR_PTR("invalid sortorder", procName, NULL);
1668
1669 n = numaGetCount(na);
1670 if ((array = numaGetFArray(na, L_COPY)) == NULL)
1671 return (NUMA *)ERROR_PTR("array not made", procName, NULL);
1672 if ((iarray = (l_float32 *)CALLOC(n, sizeof(l_float32))) == NULL)
1673 return (NUMA *)ERROR_PTR("iarray not made", procName, NULL);
1674 for (i = 0; i < n; i++)
1675 iarray[i] = i;
1676
1677 /* Shell sort */
1678 for (gap = n/2; gap > 0; gap = gap / 2) {
1679 for (i = gap; i < n; i++) {
1680 for (j = i - gap; j >= 0; j -= gap) {
1681 if ((sortorder == L_SORT_INCREASING &&
1682 array[j] > array[j + gap]) ||
1683 (sortorder == L_SORT_DECREASING &&
1684 array[j] < array[j + gap]))
1685 {
1686 tmp = array[j];
1687 array[j] = array[j + gap];
1688 array[j + gap] = tmp;
1689 tmp = iarray[j];
1690 iarray[j] = iarray[j + gap];
1691 iarray[j + gap] = tmp;
1692 }
1693 }
1694 }
1695 }
1696
1697 naisort = numaCreate(n);
1698 for (i = 0; i < n; i++)
1699 numaAddNumber(naisort, iarray[i]);
1700
1701 FREE(array);
1702 FREE(iarray);
1703 return naisort;
1704 }
1705
1706
1707 /*!
1708 * numaSortByIndex()
1709 *
1710 * Input: nas
1711 * naindex (na that maps from the new numa to the input numa)
1712 * Return: nad (sorted), or null on error
1713 */
1714 NUMA *
numaSortByIndex(NUMA * nas,NUMA * naindex)1715 numaSortByIndex(NUMA *nas,
1716 NUMA *naindex)
1717 {
1718 l_int32 i, n, index;
1719 l_float32 val;
1720 NUMA *nad;
1721
1722 PROCNAME("numaSortByIndex");
1723
1724 if (!nas)
1725 return (NUMA *)ERROR_PTR("nas not defined", procName, NULL);
1726 if (!naindex)
1727 return (NUMA *)ERROR_PTR("naindex not defined", procName, NULL);
1728
1729 n = numaGetCount(nas);
1730 nad = numaCreate(n);
1731 for (i = 0; i < n; i++) {
1732 numaGetIValue(naindex, i, &index);
1733 numaGetFValue(nas, index, &val);
1734 numaAddNumber(nad, val);
1735 }
1736
1737 return nad;
1738 }
1739
1740
1741 /*!
1742 * numaIsSorted()
1743 *
1744 * Input: nas
1745 * sortorder (L_SORT_INCREASING or L_SORT_DECREASING)
1746 * &sorted (<return> 1 if sorted; 0 if not)
1747 * Return: 1 if OK; 0 on error
1748 *
1749 * Notes:
1750 * (1) This is a quick O(n) test if nas is sorted. It is useful
1751 * in situations where the array is likely to be already
1752 * sorted, and a sort operation can be avoided.
1753 */
1754 l_int32
numaIsSorted(NUMA * nas,l_int32 sortorder,l_int32 * psorted)1755 numaIsSorted(NUMA *nas,
1756 l_int32 sortorder,
1757 l_int32 *psorted)
1758 {
1759 l_int32 i, n;
1760 l_float32 preval, val;
1761
1762 PROCNAME("numaIsSorted");
1763
1764 if (!psorted)
1765 return ERROR_INT("&sorted not defined", procName, 1);
1766 *psorted = FALSE;
1767 if (!nas)
1768 return ERROR_INT("nas not defined", procName, 1);
1769 if (sortorder != L_SORT_INCREASING && sortorder != L_SORT_DECREASING)
1770 return ERROR_INT("invalid sortorder", procName, 1);
1771
1772 n = numaGetCount(nas);
1773 numaGetFValue(nas, 0, &preval);
1774 for (i = 1; i < n; i++) {
1775 numaGetFValue(nas, i, &val);
1776 if ((sortorder == L_SORT_INCREASING && val < preval) ||
1777 (sortorder == L_SORT_DECREASING && val > preval))
1778 return 0;
1779 }
1780
1781 *psorted = TRUE;
1782 return 0;
1783 }
1784
1785
1786 /*!
1787 * numaSortPair()
1788 *
1789 * Input: nax, nay (input arrays)
1790 * sortorder (L_SORT_INCREASING or L_SORT_DECREASING)
1791 * &nasx (<return> sorted)
1792 * &naxy (<return> sorted exactly in order of nasx)
1793 * Return: 0 if OK, 1 on error
1794 *
1795 * Notes:
1796 * (1) This function sorts the two input arrays, nax and nay,
1797 * together, using nax as the key for sorting.
1798 */
1799 l_int32
numaSortPair(NUMA * nax,NUMA * nay,l_int32 sortorder,NUMA ** pnasx,NUMA ** pnasy)1800 numaSortPair(NUMA *nax,
1801 NUMA *nay,
1802 l_int32 sortorder,
1803 NUMA **pnasx,
1804 NUMA **pnasy)
1805 {
1806 l_int32 sorted;
1807 NUMA *naindex;
1808
1809 PROCNAME("numaSortPair");
1810
1811 if (!pnasx)
1812 return ERROR_INT("&nasx not defined", procName, 1);
1813 if (!pnasy)
1814 return ERROR_INT("&nasy not defined", procName, 1);
1815 *pnasx = *pnasy = NULL;
1816 if (!nax)
1817 return ERROR_INT("nax not defined", procName, 1);
1818 if (!nay)
1819 return ERROR_INT("nay not defined", procName, 1);
1820 if (sortorder != L_SORT_INCREASING && sortorder != L_SORT_DECREASING)
1821 return ERROR_INT("invalid sortorder", procName, 1);
1822
1823 numaIsSorted(nax, sortorder, &sorted);
1824 if (sorted == TRUE) {
1825 *pnasx = numaCopy(nax);
1826 *pnasy = numaCopy(nay);
1827 }
1828 else {
1829 naindex = numaGetSortIndex(nax, sortorder);
1830 *pnasx = numaSortByIndex(nax, naindex);
1831 *pnasy = numaSortByIndex(nay, naindex);
1832 numaDestroy(&naindex);
1833 }
1834
1835 return 0;
1836 }
1837
1838
1839 /*!
1840 * numaPseudorandomSequence()
1841 *
1842 * Input: size (of sequence)
1843 * seed (prime number; use 0 for default)
1844 * Return: na (pseudorandom on {0,...,size - 1}), or null on error
1845 *
1846 * Notes:
1847 * (1) Result is a permutation of the sequence of integers
1848 * from 0 to size - 1, where (seed % size) is repeatedly
1849 * added to the previous result, and the result is taken mod size.
1850 * This is not particularly random!
1851 */
1852 NUMA *
numaPseudorandomSequence(l_int32 size,l_int32 seed)1853 numaPseudorandomSequence(l_int32 size,
1854 l_int32 seed)
1855 {
1856 l_int32 i, val;
1857 NUMA *na;
1858
1859 PROCNAME("numaPseudorandomSequence");
1860
1861 if (size <= 0)
1862 return (NUMA *)ERROR_PTR("size <= 0", procName, NULL);
1863 if (seed == 0)
1864 seed = 165653;
1865
1866 if ((na = numaCreate(size)) == NULL)
1867 return (NUMA *)ERROR_PTR("na not made", procName, NULL);
1868 val = seed / 7;
1869 for (i = 0; i < size; i++) {
1870 val = (val + seed) % size;
1871 numaAddNumber(na, val);
1872 }
1873
1874 return na;
1875 }
1876
1877
1878 /*----------------------------------------------------------------------*
1879 * Functions requiring sorting *
1880 *----------------------------------------------------------------------*/
1881 /*!
1882 * numaGetRankValue()
1883 *
1884 * Input: na
1885 * fract (use 0.0 for smallest, 1.0 for largest)
1886 * &val (<return> rank val)
1887 * Return: 0 if OK; 1 on error
1888 *
1889 * Notes:
1890 * (1) Computes the rank value of the numbers in the numa, by
1891 * sorting and finding the value a fraction @fract from
1892 * the small end.
1893 */
1894 l_int32
numaGetRankValue(NUMA * na,l_float32 fract,l_float32 * pval)1895 numaGetRankValue(NUMA *na,
1896 l_float32 fract,
1897 l_float32 *pval)
1898 {
1899 l_int32 n, index;
1900 NUMA *nasort;
1901
1902 PROCNAME("numaGetRankValue");
1903
1904 if (!pval)
1905 return ERROR_INT("&val not defined", procName, 1);
1906 *pval = 0.0; /* init */
1907 if (!na)
1908 return ERROR_INT("na not defined", procName, 1);
1909 if (fract < 0.0 || fract > 1.0)
1910 return ERROR_INT("fract not in [0.0 ... 1.0]", procName, 1);
1911 n = numaGetCount(na);
1912 if (n == 0)
1913 return ERROR_INT("na empty", procName, 1);
1914
1915 if ((nasort = numaSort(NULL, na, L_SORT_INCREASING)) == NULL)
1916 return ERROR_INT("nasort not made", procName, 1);
1917 index = (l_int32)(fract * (l_float32)(n - 1) + 0.5);
1918 numaGetFValue(nasort, index, pval);
1919
1920 numaDestroy(&nasort);
1921 return 0;
1922 }
1923
1924
1925 /*!
1926 * numaGetMedian()
1927 *
1928 * Input: na
1929 * &val (<return> median val)
1930 * Return: 0 if OK; 1 on error
1931 *
1932 * Notes:
1933 * (1) Computes the median value of the numbers in the numa, by
1934 * sorting and finding the middle value in the sorted array.
1935 */
1936 l_int32
numaGetMedian(NUMA * na,l_float32 * pval)1937 numaGetMedian(NUMA *na,
1938 l_float32 *pval)
1939 {
1940 PROCNAME("numaGetMedian");
1941
1942 if (!pval)
1943 return ERROR_INT("&val not defined", procName, 1);
1944 *pval = 0.0; /* init */
1945 if (!na)
1946 return ERROR_INT("na not defined", procName, 1);
1947
1948 return numaGetRankValue(na, 0.5, pval);
1949 }
1950
1951
1952 /*!
1953 * numaGetMode()
1954 *
1955 * Input: na
1956 * &val (<return> mode val)
1957 * &count (<optional return> mode count)
1958 * Return: 0 if OK; 1 on error
1959 *
1960 * Notes:
1961 * (1) Computes the mode value of the numbers in the numa, by
1962 * sorting and finding the value of the number with the
1963 * largest count.
1964 * (2) Optionally, also returns that count.
1965 */
1966 l_int32
numaGetMode(NUMA * na,l_float32 * pval,l_int32 * pcount)1967 numaGetMode(NUMA *na,
1968 l_float32 *pval,
1969 l_int32 *pcount)
1970 {
1971 l_int32 i, n, maxcount, prevcount;
1972 l_float32 val, maxval, prevval;
1973 l_float32 *array;
1974 NUMA *nasort;
1975
1976 PROCNAME("numaGetMode");
1977
1978 if (!na)
1979 return ERROR_INT("na not defined", procName, 1);
1980 if (!pval)
1981 return ERROR_INT("&val not defined", procName, 1);
1982
1983 *pval = 0.0;
1984 if (pcount) *pcount = 0;
1985 if ((n = numaGetCount(na)) == 0)
1986 return 1;
1987
1988 if ((nasort = numaSort(NULL, na, L_SORT_DECREASING)) == NULL)
1989 return ERROR_INT("nas not made", procName, 1);
1990 array = numaGetFArray(nasort, L_NOCOPY);
1991
1992 /* Initialize with array[0] */
1993 prevval = array[0];
1994 prevcount = 1;
1995 maxval = prevval;
1996 maxcount = prevcount;
1997
1998 /* Scan the sorted array, aggregating duplicates */
1999 for (i = 1; i < n; i++) {
2000 val = array[i];
2001 if (val == prevval)
2002 prevcount++;
2003 else { /* new value */
2004 if (prevcount > maxcount) { /* new max */
2005 maxcount = prevcount;
2006 maxval = prevval;
2007 }
2008 prevval = val;
2009 prevcount = 1;
2010 }
2011 }
2012
2013 /* Was the mode the last run of elements? */
2014 if (prevcount > maxcount) {
2015 maxcount = prevcount;
2016 maxval = prevval;
2017 }
2018
2019 *pval = maxval;
2020 if (pcount)
2021 *pcount = maxcount;
2022
2023 numaDestroy(&nasort);
2024 return 0;
2025 }
2026
2027
2028 /*----------------------------------------------------------------------*
2029 * Numa combination *
2030 *----------------------------------------------------------------------*/
2031 /*!
2032 * numaJoin()
2033 *
2034 * Input: nad (dest numa; add to this one)
2035 * nas (<optional> source numa; add from this one)
2036 * istart (starting index in nas)
2037 * iend (ending index in nas; use 0 to cat all)
2038 * Return: 0 if OK, 1 on error
2039 *
2040 * Notes:
2041 * (1) istart < 0 is taken to mean 'read from the start' (istart = 0)
2042 * (2) iend <= 0 means 'read to the end'
2043 * (3) if nas == NULL, this is a no-op
2044 */
2045 l_int32
numaJoin(NUMA * nad,NUMA * nas,l_int32 istart,l_int32 iend)2046 numaJoin(NUMA *nad,
2047 NUMA *nas,
2048 l_int32 istart,
2049 l_int32 iend)
2050 {
2051 l_int32 ns, i;
2052 l_float32 val;
2053
2054 PROCNAME("numaJoin");
2055
2056 if (!nad)
2057 return ERROR_INT("nad not defined", procName, 1);
2058 if (!nas)
2059 return 0;
2060 ns = numaGetCount(nas);
2061 if (istart < 0)
2062 istart = 0;
2063 if (istart >= ns)
2064 return ERROR_INT("istart out of bounds", procName, 1);
2065 if (iend <= 0)
2066 iend = ns - 1;
2067 if (iend >= ns)
2068 return ERROR_INT("iend out of bounds", procName, 1);
2069 if (istart > iend)
2070 return ERROR_INT("istart > iend; nothing to add", procName, 1);
2071
2072 for (i = istart; i <= iend; i++) {
2073 numaGetFValue(nas, i, &val);
2074 numaAddNumber(nad, val);
2075 }
2076
2077 return 0;
2078 }
2079
2080
2081 /*!
2082 * numaaFlattenToNuma()
2083 *
2084 * Input: numaa
2085 * Return: numa, or null on error
2086 *
2087 * Notes:
2088 * (1) This 'flattens' the Numaa to a Numa, by joining successively
2089 * each Numa in the Numaa.
2090 * (2) It doesn't make any assumptions about the location of the
2091 * Numas in the Numaa array, unlike most Numaa functions.
2092 * (3) It leaves the input Numaa unchanged.
2093 */
2094 NUMA *
numaaFlattenToNuma(NUMAA * naa)2095 numaaFlattenToNuma(NUMAA *naa)
2096 {
2097 l_int32 i, nalloc;
2098 NUMA *na, *nad;
2099 NUMA **array;
2100
2101 PROCNAME("numaaFlattenToNuma");
2102
2103 if (!naa)
2104 return (NUMA *)ERROR_PTR("naa not defined", procName, NULL);
2105
2106 nalloc = naa->nalloc;
2107 array = numaaGetPtrArray(naa);
2108 nad = numaCreate(0);
2109 for (i = 0; i < nalloc; i++) {
2110 na = array[i];
2111 if (!na) continue;
2112 numaJoin(nad, na, 0, 0);
2113 }
2114
2115 return nad;
2116 }
2117
2118
2119