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