Leptonica 1.68
C Image Processing Library
|
00001 /*====================================================================* 00002 - Copyright (C) 2001 Leptonica. All rights reserved. 00003 - This software is distributed in the hope that it will be 00004 - useful, but with NO WARRANTY OF ANY KIND. 00005 - No author or distributor accepts responsibility to anyone for the 00006 - consequences of using this software, or for whether it serves any 00007 - particular purpose or works at all, unless he or she says so in 00008 - writing. Everyone is granted permission to copy, modify and 00009 - redistribute this source code, for commercial or non-commercial 00010 - purposes, with the following restrictions: (1) the origin of this 00011 - source code must not be misrepresented; (2) modified versions must 00012 - be plainly marked as such; and (3) this notice may not be removed 00013 - or altered from any source or modified source distribution. 00014 *====================================================================*/ 00015 00016 /* 00017 * numafunc1.c 00018 * 00019 * Arithmetic and logic 00020 * NUMA *numaArithOp() 00021 * NUMA *numaLogicalOp() 00022 * NUMA *numaInvert() 00023 * 00024 * Simple extractions 00025 * l_int32 numaGetMin() 00026 * l_int32 numaGetMax() 00027 * l_int32 numaGetSum() 00028 * NUMA *numaGetPartialSums() 00029 * l_int32 numaGetSumOnInterval() 00030 * l_int32 numaHasOnlyIntegers() 00031 * NUMA *numaSubsample() 00032 * NUMA *numaMakeDelta() 00033 * NUMA *numaMakeSequence() 00034 * NUMA *numaMakeConstant() 00035 * NUMA *numaAddBorder() 00036 * NUMA *numaAddSpecifiedBorder() 00037 * NUMA *numaRemoveBorder() 00038 * l_int32 numaGetNonzeroRange() 00039 * l_int32 numaGetCountRelativeToZero() 00040 * NUMA *numaClipToInterval() 00041 * NUMA *numaMakeThresholdIndicator() 00042 * NUMA *numaUniformSampling() 00043 * 00044 * Signal feature extraction 00045 * NUMA *numaLowPassIntervals() 00046 * NUMA *numaThresholdEdges() 00047 * NUMA *numaGetSpanValues() 00048 * NUMA *numaGetEdgeValues() 00049 * 00050 * Interpolation 00051 * l_int32 numaInterpolateEqxVal() 00052 * l_int32 numaInterpolateEqxInterval() 00053 * l_int32 numaInterpolateArbxVal() 00054 * l_int32 numaInterpolateArbxInterval() 00055 * 00056 * Functions requiring interpolation 00057 * l_int32 numaFitMax() 00058 * l_int32 numaDifferentiateInterval() 00059 * l_int32 numaIntegrateInterval() 00060 * 00061 * Sorting 00062 * NUMA *numaSort() 00063 * NUMA *numaGetSortIndex() 00064 * NUMA *numaSortByIndex() 00065 * l_int32 numaIsSorted() 00066 * l_int32 numaSortPair() 00067 * 00068 * Random permutation 00069 * NUMA *numaPseudorandomSequence() 00070 * NUMA *numaRandomPermutation() 00071 * 00072 * Functions requiring sorting 00073 * l_int32 numaGetRankValue() 00074 * l_int32 numaGetMedian() 00075 * l_int32 numaGetMode() 00076 * 00077 * Numa combination 00078 * l_int32 numaJoin() 00079 * NUMA *numaaFlattenToNuma() 00080 * 00081 * 00082 * Things to remember when using the Numa: 00083 * 00084 * (1) The numa is a struct, not an array. Always use accessors 00085 * (see numabasic.c), never the fields directly. 00086 * 00087 * (2) The number array holds l_float32 values. It can also 00088 * be used to store l_int32 values. See numabasic.c for 00089 * details on using the accessors. 00090 * 00091 * (3) If you use numaCreate(), no numbers are stored and the size is 0. 00092 * You have to add numbers to increase the size. 00093 * If you want to start with a numa of a fixed size, with each 00094 * entry initialized to the same value, use numaMakeConstant(). 00095 * 00096 * (4) Occasionally, in the comments we denote the i-th element of a 00097 * numa by na[i]. This is conceptual only -- the numa is not an array! 00098 */ 00099 00100 #include <stdlib.h> 00101 #include <math.h> 00102 #include "allheaders.h" 00103 00104 00105 /*----------------------------------------------------------------------* 00106 * Arithmetic and logical ops on Numas * 00107 *----------------------------------------------------------------------*/ 00108 /*! 00109 * numaArithOp() 00110 * 00111 * Input: nad (<optional> can be null or equal to na1 (in-place) 00112 * na1 00113 * na2 00114 * op (L_ARITH_ADD, L_ARITH_SUBTRACT, 00115 * L_ARITH_MULTIPLY, L_ARITH_DIVIDE) 00116 * Return: nad (always: operation applied to na1 and na2) 00117 * 00118 * Notes: 00119 * (1) The sizes of na1 and na2 must be equal. 00120 * (2) nad can only null or equal to na1. 00121 * (3) To add a constant to a numa, or to multipy a numa by 00122 * a constant, use numaTransform(). 00123 */ 00124 NUMA * 00125 numaArithOp(NUMA *nad, 00126 NUMA *na1, 00127 NUMA *na2, 00128 l_int32 op) 00129 { 00130 l_int32 i, n; 00131 l_float32 val1, val2; 00132 00133 PROCNAME("numaArithOp"); 00134 00135 if (!na1 || !na2) 00136 return (NUMA *)ERROR_PTR("na1, na2 not both defined", procName, nad); 00137 n = numaGetCount(na1); 00138 if (n != numaGetCount(na2)) 00139 return (NUMA *)ERROR_PTR("na1, na2 sizes differ", procName, nad); 00140 if (nad && nad != na1) 00141 return (NUMA *)ERROR_PTR("nad defined but not in-place", procName, nad); 00142 if (op != L_ARITH_ADD && op != L_ARITH_SUBTRACT && 00143 op != L_ARITH_MULTIPLY && op != L_ARITH_DIVIDE) 00144 return (NUMA *)ERROR_PTR("invalid op", procName, nad); 00145 if (op == L_ARITH_DIVIDE) { 00146 for (i = 0; i < n; i++) { 00147 numaGetFValue(na2, i, &val2); 00148 if (val2 == 0.0) 00149 return (NUMA *)ERROR_PTR("na2 has 0 element", procName, nad); 00150 } 00151 } 00152 00153 /* If nad is not identical to na1, make it an identical copy */ 00154 if (!nad) 00155 nad = numaCopy(na1); 00156 00157 for (i = 0; i < n; i++) { 00158 numaGetFValue(nad, i, &val1); 00159 numaGetFValue(na2, i, &val2); 00160 switch (op) { 00161 case L_ARITH_ADD: 00162 numaSetValue(nad, i, val1 + val2); 00163 break; 00164 case L_ARITH_SUBTRACT: 00165 numaSetValue(nad, i, val1 - val2); 00166 break; 00167 case L_ARITH_MULTIPLY: 00168 numaSetValue(nad, i, val1 * val2); 00169 break; 00170 case L_ARITH_DIVIDE: 00171 numaSetValue(nad, i, val1 / val2); 00172 break; 00173 default: 00174 fprintf(stderr, " Unknown arith op: %d\n", op); 00175 return nad; 00176 } 00177 } 00178 00179 return nad; 00180 } 00181 00182 00183 /*! 00184 * numaLogicalOp() 00185 * 00186 * Input: nad (<optional> can be null or equal to na1 (in-place) 00187 * na1 00188 * na2 00189 * op (L_UNION, L_INTERSECTION, L_SUBTRACTION, L_EXCLUSIVE_OR) 00190 * Return: nad (always: operation applied to na1 and na2) 00191 * 00192 * Notes: 00193 * (1) The sizes of na1 and na2 must be equal. 00194 * (2) nad can only null or equal to na1. 00195 * (3) This is intended for use with indicator arrays (0s and 1s). 00196 * Input data is extracted as integers (0 == false, anything 00197 * else == true); output results are 0 and 1. 00198 * (4) L_SUBTRACTION is subtraction of val2 from val1. For bit logical 00199 * arithmetic this is (val1 & ~val2), but because these values 00200 * are integers, we use (val1 && !val2). 00201 */ 00202 NUMA * 00203 numaLogicalOp(NUMA *nad, 00204 NUMA *na1, 00205 NUMA *na2, 00206 l_int32 op) 00207 { 00208 l_int32 i, n, val1, val2, val; 00209 00210 PROCNAME("numaLogicalOp"); 00211 00212 if (!na1 || !na2) 00213 return (NUMA *)ERROR_PTR("na1, na2 not both defined", procName, nad); 00214 n = numaGetCount(na1); 00215 if (n != numaGetCount(na2)) 00216 return (NUMA *)ERROR_PTR("na1, na2 sizes differ", procName, nad); 00217 if (nad && nad != na1) 00218 return (NUMA *)ERROR_PTR("nad defined; not in-place", procName, nad); 00219 if (op != L_UNION && op != L_INTERSECTION && 00220 op != L_SUBTRACTION && op != L_EXCLUSIVE_OR) 00221 return (NUMA *)ERROR_PTR("invalid op", procName, nad); 00222 00223 /* If nad is not identical to na1, make it an identical copy */ 00224 if (!nad) 00225 nad = numaCopy(na1); 00226 00227 for (i = 0; i < n; i++) { 00228 numaGetIValue(nad, i, &val1); 00229 numaGetIValue(na2, i, &val2); 00230 switch (op) { 00231 case L_UNION: 00232 val = (val1 || val2) ? 1 : 0; 00233 numaSetValue(nad, i, val); 00234 break; 00235 case L_INTERSECTION: 00236 val = (val1 && val2) ? 1 : 0; 00237 numaSetValue(nad, i, val); 00238 break; 00239 case L_SUBTRACTION: 00240 val = (val1 && !val2) ? 1 : 0; 00241 numaSetValue(nad, i, val); 00242 break; 00243 case L_EXCLUSIVE_OR: 00244 val = ((val1 && !val2) || (!val1 && val2)) ? 1 : 0; 00245 numaSetValue(nad, i, val); 00246 break; 00247 default: 00248 fprintf(stderr, " Unknown logical op: %d\n", op); 00249 return nad; 00250 } 00251 } 00252 00253 return nad; 00254 } 00255 00256 00257 /*! 00258 * numaInvert() 00259 * 00260 * Input: nad (<optional> can be null or equal to nas (in-place) 00261 * nas 00262 * Return: nad (always: 'inverts' nas) 00263 * 00264 * Notes: 00265 * (1) This is intended for use with indicator arrays (0s and 1s). 00266 * It gives a boolean-type output, taking the input as 00267 * an integer and inverting it: 00268 * 0 --> 1 00269 * anything else --> 0 00270 */ 00271 NUMA * 00272 numaInvert(NUMA *nad, 00273 NUMA *nas) 00274 { 00275 l_int32 i, n, val; 00276 00277 PROCNAME("numaInvert"); 00278 00279 if (!nas) 00280 return (NUMA *)ERROR_PTR("nas not defined", procName, nad); 00281 if (nad && nad != nas) 00282 return (NUMA *)ERROR_PTR("nad defined; not in-place", procName, nad); 00283 00284 if (!nad) 00285 nad = numaCopy(nas); 00286 n = numaGetCount(nad); 00287 for (i = 0; i < n; i++) { 00288 numaGetIValue(nad, i, &val); 00289 if (!val) 00290 val = 1; 00291 else 00292 val = 0; 00293 numaSetValue(nad, i, val); 00294 } 00295 00296 return nad; 00297 } 00298 00299 00300 /*----------------------------------------------------------------------* 00301 * Simple extractions * 00302 *----------------------------------------------------------------------*/ 00303 /*! 00304 * numaGetMin() 00305 * 00306 * Input: na 00307 * &minval (<optional return> min value) 00308 * &iminloc (<optional return> index of min location) 00309 * Return: 0 if OK; 1 on error 00310 */ 00311 l_int32 00312 numaGetMin(NUMA *na, 00313 l_float32 *pminval, 00314 l_int32 *piminloc) 00315 { 00316 l_int32 i, n, iminloc; 00317 l_float32 val, minval; 00318 00319 PROCNAME("numaGetMin"); 00320 00321 if (!pminval && !piminloc) 00322 return ERROR_INT("nothing to do", procName, 1); 00323 if (pminval) *pminval = 0.0; 00324 if (piminloc) *piminloc = 0; 00325 if (!na) 00326 return ERROR_INT("na not defined", procName, 1); 00327 00328 minval = +1000000000.; 00329 iminloc = 0; 00330 n = numaGetCount(na); 00331 for (i = 0; i < n; i++) { 00332 numaGetFValue(na, i, &val); 00333 if (val < minval) { 00334 minval = val; 00335 iminloc = i; 00336 } 00337 } 00338 00339 if (pminval) *pminval = minval; 00340 if (piminloc) *piminloc = iminloc; 00341 return 0; 00342 } 00343 00344 00345 /*! 00346 * numaGetMax() 00347 * 00348 * Input: na 00349 * &maxval (<optional return> max value) 00350 * &imaxloc (<optional return> index of max location) 00351 * Return: 0 if OK; 1 on error 00352 */ 00353 l_int32 00354 numaGetMax(NUMA *na, 00355 l_float32 *pmaxval, 00356 l_int32 *pimaxloc) 00357 { 00358 l_int32 i, n, imaxloc; 00359 l_float32 val, maxval; 00360 00361 PROCNAME("numaGetMax"); 00362 00363 if (!pmaxval && !pimaxloc) 00364 return ERROR_INT("nothing to do", procName, 1); 00365 if (pmaxval) *pmaxval = 0.0; 00366 if (pimaxloc) *pimaxloc = 0; 00367 if (!na) 00368 return ERROR_INT("na not defined", procName, 1); 00369 00370 maxval = -1000000000.; 00371 imaxloc = 0; 00372 n = numaGetCount(na); 00373 for (i = 0; i < n; i++) { 00374 numaGetFValue(na, i, &val); 00375 if (val > maxval) { 00376 maxval = val; 00377 imaxloc = i; 00378 } 00379 } 00380 00381 if (pmaxval) *pmaxval = maxval; 00382 if (pimaxloc) *pimaxloc = imaxloc; 00383 return 0; 00384 } 00385 00386 00387 /*! 00388 * numaGetSum() 00389 * 00390 * Input: na 00391 * &sum (<return> sum of values) 00392 * Return: 0 if OK, 1 on error 00393 */ 00394 l_int32 00395 numaGetSum(NUMA *na, 00396 l_float32 *psum) 00397 { 00398 l_int32 i, n; 00399 l_float32 val, sum; 00400 00401 PROCNAME("numaGetSum"); 00402 00403 if (!na) 00404 return ERROR_INT("na not defined", procName, 1); 00405 if (!psum) 00406 return ERROR_INT("&sum not defined", procName, 1); 00407 00408 sum = 0.0; 00409 n = numaGetCount(na); 00410 for (i = 0; i < n; i++) { 00411 numaGetFValue(na, i, &val); 00412 sum += val; 00413 } 00414 *psum = sum; 00415 return 0; 00416 } 00417 00418 00419 /*! 00420 * numaGetPartialSums() 00421 * 00422 * Input: na 00423 * Return: nasum, or null on error 00424 * 00425 * Notes: 00426 * (1) nasum[i] is the sum for all j <= i of na[j]. 00427 * So nasum[0] = na[0]. 00428 * (2) If you want to generate a rank function, where rank[0] - 0.0, 00429 * insert a 0.0 at the beginning of the nasum array. 00430 */ 00431 NUMA * 00432 numaGetPartialSums(NUMA *na) 00433 { 00434 l_int32 i, n; 00435 l_float32 val, sum; 00436 NUMA *nasum; 00437 00438 PROCNAME("numaGetPartialSums"); 00439 00440 if (!na) 00441 return (NUMA *)ERROR_PTR("na not defined", procName, NULL); 00442 00443 n = numaGetCount(na); 00444 nasum = numaCreate(n); 00445 sum = 0.0; 00446 for (i = 0; i < n; i++) { 00447 numaGetFValue(na, i, &val); 00448 sum += val; 00449 numaAddNumber(nasum, sum); 00450 } 00451 return nasum; 00452 } 00453 00454 00455 /*! 00456 * numaGetSumOnInterval() 00457 * 00458 * Input: na 00459 * first (beginning index) 00460 * last (final index) 00461 * &sum (<return> sum of values in the index interval range) 00462 * Return: 0 if OK, 1 on error 00463 */ 00464 l_int32 00465 numaGetSumOnInterval(NUMA *na, 00466 l_int32 first, 00467 l_int32 last, 00468 l_float32 *psum) 00469 { 00470 l_int32 i, n, truelast; 00471 l_float32 val, sum; 00472 00473 PROCNAME("numaGetSumOnInterval"); 00474 00475 if (!na) 00476 return ERROR_INT("na not defined", procName, 1); 00477 if (!psum) 00478 return ERROR_INT("&sum not defined", procName, 1); 00479 *psum = 0.0; 00480 00481 sum = 0.0; 00482 n = numaGetCount(na); 00483 if (first >= n) /* not an error */ 00484 return 0; 00485 truelast = L_MIN(last, n - 1); 00486 00487 for (i = first; i <= truelast; i++) { 00488 numaGetFValue(na, i, &val); 00489 sum += val; 00490 } 00491 *psum = sum; 00492 return 0; 00493 } 00494 00495 00496 /*! 00497 * numaHasOnlyIntegers() 00498 * 00499 * Input: na 00500 * maxsamples (maximum number of samples to check) 00501 * &allints (<return> 1 if all sampled values are ints; else 0) 00502 * Return: 0 if OK, 1 on error 00503 * 00504 * Notes: 00505 * (1) Set @maxsamples == 0 to check every integer in na. Otherwise, 00506 * this samples no more than @maxsamples. 00507 */ 00508 l_int32 00509 numaHasOnlyIntegers(NUMA *na, 00510 l_int32 maxsamples, 00511 l_int32 *pallints) 00512 { 00513 l_int32 i, n, incr; 00514 l_float32 val; 00515 00516 PROCNAME("numaHasOnlyIntegers"); 00517 00518 if (!pallints) 00519 return ERROR_INT("&allints not defined", procName, 1); 00520 *pallints = TRUE; 00521 if (!na) 00522 return ERROR_INT("na not defined", procName, 1); 00523 00524 if ((n = numaGetCount(na)) == 0) 00525 return ERROR_INT("na empty", procName, 1); 00526 if (maxsamples <= 0) 00527 incr = 1; 00528 else 00529 incr = (l_int32)((n + maxsamples - 1) / maxsamples); 00530 for (i = 0; i < n; i += incr) { 00531 numaGetFValue(na, i, &val); 00532 if (val != (l_int32)val) { 00533 *pallints = FALSE; 00534 return 0; 00535 } 00536 } 00537 00538 return 0; 00539 } 00540 00541 00542 /*! 00543 * numaSubsample() 00544 * 00545 * Input: nas 00546 * subfactor (subsample factor, >= 1) 00547 * Return: nad (evenly sampled values from nas), or null on error 00548 */ 00549 NUMA * 00550 numaSubsample(NUMA *nas, 00551 l_int32 subfactor) 00552 { 00553 l_int32 i, n; 00554 l_float32 val; 00555 NUMA *nad; 00556 00557 PROCNAME("numaSubsample"); 00558 00559 if (!nas) 00560 return (NUMA *)ERROR_PTR("nas not defined", procName, NULL); 00561 if (subfactor < 1) 00562 return (NUMA *)ERROR_PTR("subfactor < 1", procName, NULL); 00563 00564 nad = numaCreate(0); 00565 n = numaGetCount(nas); 00566 for (i = 0; i < n; i++) { 00567 if (i % subfactor != 0) continue; 00568 numaGetFValue(nas, i, &val); 00569 numaAddNumber(nad, val); 00570 } 00571 00572 return nad; 00573 } 00574 00575 00576 /*! 00577 * numaMakeDelta() 00578 * 00579 * Input: nas (input numa) 00580 * Return: numa (of difference values val[i+1] - val[i]), 00581 * or null on error 00582 */ 00583 NUMA * 00584 numaMakeDelta(NUMA *nas) 00585 { 00586 l_int32 i, n, prev, cur; 00587 NUMA *nad; 00588 00589 PROCNAME("numaMakeDelta"); 00590 00591 if (!nas) 00592 return (NUMA *)ERROR_PTR("nas not defined", procName, NULL); 00593 n = numaGetCount(nas); 00594 nad = numaCreate(n - 1); 00595 prev = 0; 00596 for (i = 1; i < n; i++) { 00597 numaGetIValue(nas, i, &cur); 00598 numaAddNumber(nad, cur - prev); 00599 prev = cur; 00600 } 00601 return nad; 00602 } 00603 00604 00605 /*! 00606 * numaMakeSequence() 00607 * 00608 * Input: startval 00609 * increment 00610 * size (of sequence) 00611 * Return: numa of sequence of evenly spaced values, or null on error 00612 */ 00613 NUMA * 00614 numaMakeSequence(l_float32 startval, 00615 l_float32 increment, 00616 l_int32 size) 00617 { 00618 l_int32 i; 00619 l_float32 val; 00620 NUMA *na; 00621 00622 PROCNAME("numaMakeSequence"); 00623 00624 if ((na = numaCreate(size)) == NULL) 00625 return (NUMA *)ERROR_PTR("na not made", procName, NULL); 00626 00627 for (i = 0; i < size; i++) { 00628 val = startval + i * increment; 00629 numaAddNumber(na, val); 00630 } 00631 00632 return na; 00633 } 00634 00635 00636 /*! 00637 * numaMakeConstant() 00638 * 00639 * Input: val 00640 * size (of numa) 00641 * Return: numa of given size with all entries equal to 'val', 00642 * or null on error 00643 */ 00644 NUMA * 00645 numaMakeConstant(l_float32 val, 00646 l_int32 size) 00647 { 00648 return numaMakeSequence(val, 0.0, size); 00649 } 00650 00651 00652 /*! 00653 * numaAddBorder() 00654 * 00655 * Input: nas 00656 * left, right (number of elements to add on each side) 00657 * val (initialize border elements) 00658 * Return: nad (with added elements at left and right), or null on error 00659 */ 00660 NUMA * 00661 numaAddBorder(NUMA *nas, 00662 l_int32 left, 00663 l_int32 right, 00664 l_float32 val) 00665 { 00666 l_int32 i, n, len; 00667 l_float32 startx, delx; 00668 l_float32 *fas, *fad; 00669 NUMA *nad; 00670 00671 PROCNAME("numaAddBorder"); 00672 00673 if (!nas) 00674 return (NUMA *)ERROR_PTR("nas not defined", procName, NULL); 00675 if (left < 0) left = 0; 00676 if (right < 0) right = 0; 00677 if (left == 0 && right == 0) 00678 return numaCopy(nas); 00679 00680 n = numaGetCount(nas); 00681 len = n + left + right; 00682 nad = numaMakeConstant(val, len); 00683 numaGetXParameters(nas, &startx, &delx); 00684 numaSetXParameters(nad, startx - delx * left, delx); 00685 fas = numaGetFArray(nas, L_NOCOPY); 00686 fad = numaGetFArray(nad, L_NOCOPY); 00687 for (i = 0; i < n; i++) 00688 fad[left + i] = fas[i]; 00689 00690 return nad; 00691 } 00692 00693 00694 /*! 00695 * numaAddSpecifiedBorder() 00696 * 00697 * Input: nas 00698 * left, right (number of elements to add on each side) 00699 * type (L_EXTENDED_BORDER, L_MIRRORED_BORDER) 00700 * Return: nad (with added elements at left and right), or null on error 00701 */ 00702 NUMA * 00703 numaAddSpecifiedBorder(NUMA *nas, 00704 l_int32 left, 00705 l_int32 right, 00706 l_int32 type) 00707 { 00708 l_int32 i, n; 00709 l_float32 *fa; 00710 NUMA *nad; 00711 00712 PROCNAME("numaAddSpecifiedBorder"); 00713 00714 if (!nas) 00715 return (NUMA *)ERROR_PTR("nas not defined", procName, NULL); 00716 if (left < 0) left = 0; 00717 if (right < 0) right = 0; 00718 if (left == 0 && right == 0) 00719 return numaCopy(nas); 00720 if (type != L_EXTENDED_BORDER && type != L_MIRRORED_BORDER) 00721 return (NUMA *)ERROR_PTR("invalid type", procName, NULL); 00722 n = numaGetCount(nas); 00723 if (type == L_MIRRORED_BORDER && (left > n || right > n)) 00724 return (NUMA *)ERROR_PTR("border too large", procName, NULL); 00725 00726 nad = numaAddBorder(nas, left, right, 0); 00727 n = numaGetCount(nad); 00728 fa = numaGetFArray(nad, L_NOCOPY); 00729 if (type == L_EXTENDED_BORDER) { 00730 for (i = 0; i < left; i++) 00731 fa[i] = fa[left]; 00732 for (i = n - right; i < n; i++) 00733 fa[i] = fa[n - right - 1]; 00734 } 00735 else { /* type == L_MIRRORED_BORDER */ 00736 for (i = 0; i < left; i++) 00737 fa[i] = fa[2 * left - 1 - i]; 00738 for (i = 0; i < right; i++) 00739 fa[n - right + i] = fa[n - right - i - 1]; 00740 } 00741 00742 return nad; 00743 } 00744 00745 00746 /*! 00747 * numaRemoveBorder() 00748 * 00749 * Input: nas 00750 * left, right (number of elements to remove from each side) 00751 * Return: nad (with removed elements at left and right), or null on error 00752 */ 00753 NUMA * 00754 numaRemoveBorder(NUMA *nas, 00755 l_int32 left, 00756 l_int32 right) 00757 { 00758 l_int32 i, n, len; 00759 l_float32 startx, delx; 00760 l_float32 *fas, *fad; 00761 NUMA *nad; 00762 00763 PROCNAME("numaRemoveBorder"); 00764 00765 if (!nas) 00766 return (NUMA *)ERROR_PTR("nas not defined", procName, NULL); 00767 if (left < 0) left = 0; 00768 if (right < 0) right = 0; 00769 if (left == 0 && right == 0) 00770 return numaCopy(nas); 00771 00772 n = numaGetCount(nas); 00773 if ((len = n - left - right) < 0) 00774 return (NUMA *)ERROR_PTR("len < 0 after removal", procName, NULL); 00775 nad = numaMakeConstant(0, len); 00776 numaGetXParameters(nas, &startx, &delx); 00777 numaSetXParameters(nad, startx + delx * left, delx); 00778 fas = numaGetFArray(nas, L_NOCOPY); 00779 fad = numaGetFArray(nad, L_NOCOPY); 00780 for (i = 0; i < len; i++) 00781 fad[i] = fas[left + i]; 00782 00783 return nad; 00784 } 00785 00786 00787 /*! 00788 * numaGetNonzeroRange() 00789 * 00790 * Input: numa 00791 * eps (largest value considered to be zero) 00792 * &first, &last (<return> interval of array indices 00793 * where values are nonzero) 00794 * Return: 0 if OK, 1 on error or if no nonzero range is found. 00795 */ 00796 l_int32 00797 numaGetNonzeroRange(NUMA *na, 00798 l_float32 eps, 00799 l_int32 *pfirst, 00800 l_int32 *plast) 00801 { 00802 l_int32 n, i, found; 00803 l_float32 val; 00804 00805 PROCNAME("numaGetNonzeroRange"); 00806 00807 if (!na) 00808 return ERROR_INT("na not defined", procName, 1); 00809 if (!pfirst || !plast) 00810 return ERROR_INT("pfirst and plast not both defined", procName, 1); 00811 n = numaGetCount(na); 00812 found = FALSE; 00813 for (i = 0; i < n; i++) { 00814 numaGetFValue(na, i, &val); 00815 if (val > eps) { 00816 found = TRUE; 00817 break; 00818 } 00819 } 00820 if (!found) { 00821 *pfirst = n - 1; 00822 *plast = 0; 00823 return 1; 00824 } 00825 00826 *pfirst = i; 00827 for (i = n - 1; i >= 0; i--) { 00828 numaGetFValue(na, i, &val); 00829 if (val > eps) 00830 break; 00831 } 00832 *plast = i; 00833 return 0; 00834 } 00835 00836 00837 /*! 00838 * numaGetCountRelativeToZero() 00839 * 00840 * Input: numa 00841 * type (L_LESS_THAN_ZERO, L_EQUAL_TO_ZERO, L_GREATER_THAN_ZERO) 00842 * &count (<return> count of values of given type) 00843 * Return: 0 if OK, 1 on error 00844 */ 00845 l_int32 00846 numaGetCountRelativeToZero(NUMA *na, 00847 l_int32 type, 00848 l_int32 *pcount) 00849 { 00850 l_int32 n, i, count; 00851 l_float32 val; 00852 00853 PROCNAME("numaGetCountRelativeToZero"); 00854 00855 if (!pcount) 00856 return ERROR_INT("&count not defined", procName, 1); 00857 *pcount = 0; 00858 if (!na) 00859 return ERROR_INT("na not defined", procName, 1); 00860 n = numaGetCount(na); 00861 for (i = 0, count = 0; i < n; i++) { 00862 numaGetFValue(na, i, &val); 00863 if (type == L_LESS_THAN_ZERO && val < 0.0) 00864 count++; 00865 else if (type == L_EQUAL_TO_ZERO && val == 0.0) 00866 count++; 00867 else if (type == L_GREATER_THAN_ZERO && val > 0.0) 00868 count++; 00869 } 00870 00871 *pcount = count; 00872 return 0; 00873 } 00874 00875 00876 /*! 00877 * numaClipToInterval() 00878 * 00879 * Input: numa 00880 * first, last (clipping interval) 00881 * Return: numa with the same values as the input, but clipped 00882 * to the specified interval 00883 * 00884 * Note: If you want the indices of the array values to be unchanged, 00885 * use first = 0. 00886 * Usage: This is useful to clip a histogram that has a few nonzero 00887 * values to its nonzero range. 00888 */ 00889 NUMA * 00890 numaClipToInterval(NUMA *nas, 00891 l_int32 first, 00892 l_int32 last) 00893 { 00894 l_int32 n, i, truelast; 00895 l_float32 val; 00896 NUMA *nad; 00897 00898 PROCNAME("numaClipToInterval"); 00899 00900 if (!nas) 00901 return (NUMA *)ERROR_PTR("nas not defined", procName, NULL); 00902 if (first > last) 00903 return (NUMA *)ERROR_PTR("range not valid", procName, NULL); 00904 00905 n = numaGetCount(nas); 00906 if (first >= n) 00907 return (NUMA *)ERROR_PTR("no elements in range", procName, NULL); 00908 truelast = L_MIN(last, n - 1); 00909 if ((nad = numaCreate(truelast - first + 1)) == NULL) 00910 return (NUMA *)ERROR_PTR("nad not made", procName, NULL); 00911 for (i = first; i <= truelast; i++) { 00912 numaGetFValue(nas, i, &val); 00913 numaAddNumber(nad, val); 00914 } 00915 00916 return nad; 00917 } 00918 00919 00920 /*! 00921 * numaMakeThresholdIndicator() 00922 * 00923 * Input: nas (input numa) 00924 * thresh (threshold value) 00925 * type (L_SELECT_IF_LT, L_SELECT_IF_GT, 00926 * L_SELECT_IF_LTE, L_SELECT_IF_GTE) 00927 * Output: nad (indicator array: values are 0 and 1) 00928 * 00929 * Notes: 00930 * (1) For each element in nas, if the constraint given by 'type' 00931 * correctly specifies its relation to thresh, a value of 1 00932 * is recorded in nad. 00933 */ 00934 NUMA * 00935 numaMakeThresholdIndicator(NUMA *nas, 00936 l_float32 thresh, 00937 l_int32 type) 00938 { 00939 l_int32 n, i, ival; 00940 l_float32 fval; 00941 NUMA *nai; 00942 00943 PROCNAME("numaMakeThresholdIndicator"); 00944 00945 if (!nas) 00946 return (NUMA *)ERROR_PTR("nas not defined", procName, NULL); 00947 n = numaGetCount(nas); 00948 nai = numaCreate(n); 00949 for (i = 0; i < n; i++) { 00950 numaGetFValue(nas, i, &fval); 00951 ival = 0; 00952 switch (type) 00953 { 00954 case L_SELECT_IF_LT: 00955 if (fval < thresh) ival = 1; 00956 break; 00957 case L_SELECT_IF_GT: 00958 if (fval > thresh) ival = 1; 00959 break; 00960 case L_SELECT_IF_LTE: 00961 if (fval <= thresh) ival = 1; 00962 break; 00963 case L_SELECT_IF_GTE: 00964 if (fval >= thresh) ival = 1; 00965 break; 00966 default: 00967 numaDestroy(&nai); 00968 return (NUMA *)ERROR_PTR("invalid type", procName, NULL); 00969 } 00970 numaAddNumber(nai, ival); 00971 } 00972 00973 return nai; 00974 } 00975 00976 00977 /*! 00978 * numaUniformSampling() 00979 * 00980 * Input: nas (input numa) 00981 * nsamp (number of samples) 00982 * Output: nad (resampled array), or null on error 00983 * 00984 * Notes: 00985 * (1) This resamples the values in the array, using @nsamp 00986 * equal divisions. 00987 */ 00988 NUMA * 00989 numaUniformSampling(NUMA *nas, 00990 l_int32 nsamp) 00991 { 00992 l_int32 n, i, j, ileft, iright; 00993 l_float32 left, right, binsize, lfract, rfract, sum, startx, delx; 00994 l_float32 *array; 00995 NUMA *nad; 00996 00997 PROCNAME("numaUniformSampling"); 00998 00999 if (!nas) 01000 return (NUMA *)ERROR_PTR("nas not defined", procName, NULL); 01001 if (nsamp <= 0) 01002 return (NUMA *)ERROR_PTR("nsamp must be > 0", procName, NULL); 01003 01004 n = numaGetCount(nas); 01005 nad = numaCreate(nsamp); 01006 array = numaGetFArray(nas, L_NOCOPY); 01007 binsize = (l_float32)n / (l_float32)nsamp; 01008 numaGetXParameters(nas, &startx, &delx); 01009 numaSetXParameters(nad, startx, binsize * delx); 01010 left = 0.0; 01011 for (i = 0; i < nsamp; i++) { 01012 sum = 0.0; 01013 right = left + binsize; 01014 ileft = (l_int32)left; 01015 lfract = 1.0 - left + ileft; 01016 if (lfract >= 1.0) /* on left bin boundary */ 01017 lfract = 0.0; 01018 iright = (l_int32)right; 01019 rfract = right - iright; 01020 iright = L_MIN(iright, n - 1); 01021 if (ileft == iright) { /* both are within the same original sample */ 01022 sum += (lfract + rfract - 1.0) * array[ileft]; 01023 } 01024 else { 01025 if (lfract > 0.0001) /* left fraction */ 01026 sum += lfract * array[ileft]; 01027 if (rfract > 0.0001) /* right fraction */ 01028 sum += rfract * array[iright]; 01029 for (j = ileft + 1; j < iright; j++) /* entire pixels */ 01030 sum += array[j]; 01031 } 01032 01033 numaAddNumber(nad, sum); 01034 left = right; 01035 } 01036 return nad; 01037 } 01038 01039 01040 /*----------------------------------------------------------------------* 01041 * Signal feature extraction * 01042 *----------------------------------------------------------------------*/ 01043 /*! 01044 * numaLowPassIntervals() 01045 * 01046 * Input: nas (input numa) 01047 * thresh (threshold fraction of max; in [0.0 ... 1.0]) 01048 * maxn (for normalizing; set maxn = 0.0 to use the max in nas) 01049 * Output: nad (interval abscissa pairs), or null on error 01050 * 01051 * Notes: 01052 * (1) For each interval where the value is less than a specified 01053 * fraction of the maximum, this records the left and right "x" 01054 * value. 01055 */ 01056 NUMA * 01057 numaLowPassIntervals(NUMA *nas, 01058 l_float32 thresh, 01059 l_float32 maxn) 01060 { 01061 l_int32 n, i, inrun; 01062 l_float32 maxval, threshval, fval, startx, delx, x0, x1; 01063 NUMA *nad; 01064 01065 PROCNAME("numaLowPassIntervals"); 01066 01067 if (!nas) 01068 return (NUMA *)ERROR_PTR("nas not defined", procName, NULL); 01069 if (thresh < 0.0 || thresh > 1.0) 01070 return (NUMA *)ERROR_PTR("invalid thresh", procName, NULL); 01071 01072 /* The input threshold is a fraction of the max. 01073 * The first entry in nad is the value of the max. */ 01074 n = numaGetCount(nas); 01075 if (maxn == 0.0) 01076 numaGetMax(nas, &maxval, NULL); 01077 else 01078 maxval = maxn; 01079 numaGetXParameters(nas, &startx, &delx); 01080 threshval = thresh * maxval; 01081 nad = numaCreate(0); 01082 numaAddNumber(nad, maxval); 01083 01084 /* Write pairs of pts (x0, x1) for the intervals */ 01085 inrun = FALSE; 01086 for (i = 0; i < n; i++) { 01087 numaGetFValue(nas, i, &fval); 01088 if (fval < threshval && inrun == FALSE) { /* start a new run */ 01089 inrun = TRUE; 01090 x0 = startx + i * delx; 01091 } 01092 else if (fval > threshval && inrun == TRUE) { /* end the run */ 01093 inrun = FALSE; 01094 x1 = startx + i * delx; 01095 numaAddNumber(nad, x0); 01096 numaAddNumber(nad, x1); 01097 } 01098 } 01099 if (inrun == TRUE) { /* must end the last run */ 01100 x1 = startx + (n - 1) * delx; 01101 numaAddNumber(nad, x0); 01102 numaAddNumber(nad, x1); 01103 } 01104 01105 return nad; 01106 } 01107 01108 01109 /*! 01110 * numaThresholdEdges() 01111 * 01112 * Input: nas (input numa) 01113 * thresh1 (low threshold as fraction of max; in [0.0 ... 1.0]) 01114 * thresh2 (high threshold as fraction of max; in [0.0 ... 1.0]) 01115 * maxn (for normalizing; set maxn = 0.0 to use the max in nas) 01116 * Output: nad (edge interval triplets), or null on error 01117 * 01118 * Notes: 01119 * (1) For each edge interval, where where the value is less 01120 * than @thresh1 on one side, greater than @thresh2 on 01121 * the other, and between these thresholds throughout the 01122 * interval, this records a triplet of values: the 01123 * 'left' and 'right' edges, and either +1 or -1, depending 01124 * on whether the edge is rising or falling. 01125 * (2) No assumption is made about the value outside the array, 01126 * so if the value at the array edge is between the threshold 01127 * values, it is not considered part of an edge. We start 01128 * looking for edge intervals only after leaving the thresholded 01129 * band. 01130 */ 01131 NUMA * 01132 numaThresholdEdges(NUMA *nas, 01133 l_float32 thresh1, 01134 l_float32 thresh2, 01135 l_float32 maxn) 01136 { 01137 l_int32 n, i, istart, inband, output, sign; 01138 l_int32 startbelow, below, above, belowlast, abovelast; 01139 l_float32 maxval, threshval1, threshval2, fval, startx, delx, x0, x1; 01140 NUMA *nad; 01141 01142 PROCNAME("numaThresholdEdges"); 01143 01144 if (!nas) 01145 return (NUMA *)ERROR_PTR("nas not defined", procName, NULL); 01146 if (thresh1 < 0.0 || thresh1 > 1.0 || thresh2 < 0.0 || thresh2 > 1.0) 01147 return (NUMA *)ERROR_PTR("invalid thresholds", procName, NULL); 01148 if (thresh2 < thresh1) 01149 return (NUMA *)ERROR_PTR("thresh2 < thresh1", procName, NULL); 01150 01151 /* The input thresholds are fractions of the max. 01152 * The first entry in nad is the value of the max used 01153 * here for normalization. */ 01154 n = numaGetCount(nas); 01155 if (maxn == 0.0) 01156 numaGetMax(nas, &maxval, NULL); 01157 else 01158 maxval = maxn; 01159 numaGetMax(nas, &maxval, NULL); 01160 numaGetXParameters(nas, &startx, &delx); 01161 threshval1 = thresh1 * maxval; 01162 threshval2 = thresh2 * maxval; 01163 nad = numaCreate(0); 01164 numaAddNumber(nad, maxval); 01165 01166 /* Write triplets of pts (x0, x1, sign) for the edges. 01167 * First make sure we start search from outside the band. 01168 * Only one of {belowlast, abovelast} is true. */ 01169 for (i = 0; i < n; i++) { 01170 istart = i; 01171 numaGetFValue(nas, i, &fval); 01172 belowlast = (fval < threshval1) ? TRUE : FALSE; 01173 abovelast = (fval > threshval2) ? TRUE : FALSE; 01174 if (belowlast == TRUE || abovelast == TRUE) 01175 break; 01176 } 01177 if (istart == n) /* no intervals found */ 01178 return nad; 01179 01180 /* x0 and x1 can only be set from outside the edge. 01181 * They are the values just before entering the band, 01182 * and just after entering the band. We can jump through 01183 * the band, in which case they differ by one index in nas. */ 01184 inband = FALSE; 01185 startbelow = belowlast; /* one of these is true */ 01186 output = FALSE; 01187 x0 = startx + istart * delx; 01188 for (i = istart + 1; i < n; i++) { 01189 numaGetFValue(nas, i, &fval); 01190 below = (fval < threshval1) ? TRUE : FALSE; 01191 above = (fval > threshval2) ? TRUE : FALSE; 01192 if (!inband && belowlast && above) { /* full jump up */ 01193 x1 = startx + i * delx; 01194 sign = 1; 01195 startbelow = FALSE; /* for the next transition */ 01196 output = TRUE; 01197 } 01198 else if (!inband && abovelast && below) { /* full jump down */ 01199 x1 = startx + i * delx; 01200 sign = -1; 01201 startbelow = TRUE; /* for the next transition */ 01202 output = TRUE; 01203 } 01204 else if (inband && startbelow && above) { /* exit rising; success */ 01205 x1 = startx + i * delx; 01206 sign = 1; 01207 inband = FALSE; 01208 startbelow = FALSE; /* for the next transition */ 01209 output = TRUE; 01210 } 01211 else if (inband && !startbelow && below) { /* exit falling; success */ 01212 x1 = startx + i * delx; 01213 sign = -1; 01214 inband = FALSE; 01215 startbelow = TRUE; /* for the next transition */ 01216 output = TRUE; 01217 } 01218 else if (inband && !startbelow && above) { /* exit rising; failure */ 01219 x0 = startx + i * delx; 01220 inband = FALSE; 01221 } 01222 else if (inband && startbelow && below) { /* exit falling; failure */ 01223 x0 = startx + i * delx; 01224 inband = FALSE; 01225 } 01226 else if (!inband && !above && !below) { /* enter */ 01227 inband = TRUE; 01228 startbelow = belowlast; 01229 } 01230 else if (!inband && (above || below)) { /* outside and remaining */ 01231 x0 = startx + i * delx; /* update position */ 01232 } 01233 belowlast = below; 01234 abovelast = above; 01235 if (output) { /* we have exited; save new x0 */ 01236 numaAddNumber(nad, x0); 01237 numaAddNumber(nad, x1); 01238 numaAddNumber(nad, sign); 01239 output = FALSE; 01240 x0 = startx + i * delx; 01241 } 01242 } 01243 01244 return nad; 01245 } 01246 01247 01248 /*! 01249 * numaGetSpanValues() 01250 * 01251 * Input: na (numa that is output of numaLowPassIntervals()) 01252 * span (span number, zero-based) 01253 * &start (<optional return> location of start of transition) 01254 * &end (<optional return> location of end of transition) 01255 * Output: 0 if OK, 1 on error 01256 */ 01257 l_int32 01258 numaGetSpanValues(NUMA *na, 01259 l_int32 span, 01260 l_int32 *pstart, 01261 l_int32 *pend) 01262 { 01263 l_int32 n, nspans; 01264 01265 PROCNAME("numaGetSpanValues"); 01266 01267 if (!na) 01268 return ERROR_INT("na not defined", procName, 1); 01269 n = numaGetCount(na); 01270 if (n % 2 != 1) 01271 return ERROR_INT("n is not odd", procName, 1); 01272 nspans = n / 2; 01273 if (nspans < 0 || span >= nspans) 01274 return ERROR_INT("invalid span", procName, 1); 01275 01276 if (pstart) numaGetIValue(na, 2 * span + 1, pstart); 01277 if (pend) numaGetIValue(na, 2 * span + 2, pend); 01278 return 0; 01279 } 01280 01281 01282 /*! 01283 * numaGetEdgeValues() 01284 * 01285 * Input: na (numa that is output of numaThresholdEdges()) 01286 * edge (edge number, zero-based) 01287 * &start (<optional return> location of start of transition) 01288 * &end (<optional return> location of end of transition) 01289 * &sign (<optional return> transition sign: +1 is rising, 01290 * -1 is falling) 01291 * Output: 0 if OK, 1 on error 01292 */ 01293 l_int32 01294 numaGetEdgeValues(NUMA *na, 01295 l_int32 edge, 01296 l_int32 *pstart, 01297 l_int32 *pend, 01298 l_int32 *psign) 01299 { 01300 l_int32 n, nedges; 01301 01302 PROCNAME("numaGetEdgeValues"); 01303 01304 if (!na) 01305 return ERROR_INT("na not defined", procName, 1); 01306 n = numaGetCount(na); 01307 if (n % 3 != 1) 01308 return ERROR_INT("n % 3 is not 1", procName, 1); 01309 nedges = (n - 1) / 3; 01310 if (edge < 0 || edge >= nedges) 01311 return ERROR_INT("invalid edge", procName, 1); 01312 01313 if (pstart) numaGetIValue(na, 3 * edge + 1, pstart); 01314 if (pend) numaGetIValue(na, 3 * edge + 2, pend); 01315 if (psign) numaGetIValue(na, 3 * edge + 3, psign); 01316 return 0; 01317 } 01318 01319 01320 /*----------------------------------------------------------------------* 01321 * Interpolation * 01322 *----------------------------------------------------------------------*/ 01323 /*! 01324 * numaInterpolateEqxVal() 01325 * 01326 * Input: startx (xval corresponding to first element in array) 01327 * deltax (x increment between array elements) 01328 * nay (numa of ordinate values, assumed equally spaced) 01329 * type (L_LINEAR_INTERP, L_QUADRATIC_INTERP) 01330 * xval 01331 * &yval (<return> interpolated value) 01332 * Return: 0 if OK, 1 on error (e.g., if xval is outside range) 01333 * 01334 * Notes: 01335 * (1) Considering nay as a function of x, the x values 01336 * are equally spaced 01337 * (2) Caller should check for valid return. 01338 * 01339 * For linear Lagrangian interpolation (through 2 data pts): 01340 * y(x) = y1(x-x2)/(x1-x2) + y2(x-x1)/(x2-x1) 01341 * 01342 * For quadratic Lagrangian interpolation (through 3 data pts): 01343 * y(x) = y1(x-x2)(x-x3)/((x1-x2)(x1-x3)) + 01344 * y2(x-x1)(x-x3)/((x2-x1)(x2-x3)) + 01345 * y3(x-x1)(x-x2)/((x3-x1)(x3-x2)) 01346 * 01347 */ 01348 l_int32 01349 numaInterpolateEqxVal(l_float32 startx, 01350 l_float32 deltax, 01351 NUMA *nay, 01352 l_int32 type, 01353 l_float32 xval, 01354 l_float32 *pyval) 01355 { 01356 l_int32 i, n, i1, i2, i3; 01357 l_float32 x1, x2, x3, fy1, fy2, fy3, d1, d2, d3, del, fi, maxx; 01358 l_float32 *fa; 01359 01360 PROCNAME("numaInterpolateEqxVal"); 01361 01362 if (!pyval) 01363 return ERROR_INT("&yval not defined", procName, 1); 01364 *pyval = 0.0; 01365 if (!nay) 01366 return ERROR_INT("nay not defined", procName, 1); 01367 if (deltax <= 0.0) 01368 return ERROR_INT("deltax not > 0", procName, 1); 01369 if (type != L_LINEAR_INTERP && type != L_QUADRATIC_INTERP) 01370 return ERROR_INT("invalid interp type", procName, 1); 01371 n = numaGetCount(nay); 01372 if (n < 2) 01373 return ERROR_INT("not enough points", procName, 1); 01374 if (type == L_QUADRATIC_INTERP && n == 2) { 01375 type = L_LINEAR_INTERP; 01376 L_WARNING("only 2 points; using linear interp", procName); 01377 } 01378 maxx = startx + deltax * (n - 1); 01379 if (xval < startx || xval > maxx) 01380 return ERROR_INT("xval is out of bounds", procName, 1); 01381 01382 fa = numaGetFArray(nay, L_NOCOPY); 01383 fi = (xval - startx) / deltax; 01384 i = (l_int32)fi; 01385 del = fi - i; 01386 if (del == 0.0) { /* no interpolation required */ 01387 *pyval = fa[i]; 01388 return 0; 01389 } 01390 01391 if (type == L_LINEAR_INTERP) { 01392 *pyval = fa[i] + del * (fa[i + 1] - fa[i]); 01393 return 0; 01394 } 01395 01396 /* Quadratic interpolation */ 01397 d1 = d3 = 0.5 / (deltax * deltax); 01398 d2 = -2. * d1; 01399 if (i == 0) { 01400 i1 = i; 01401 i2 = i + 1; 01402 i3 = i + 2; 01403 } 01404 else { 01405 i1 = i - 1; 01406 i2 = i; 01407 i3 = i + 1; 01408 } 01409 x1 = startx + i1 * deltax; 01410 x2 = startx + i2 * deltax; 01411 x3 = startx + i3 * deltax; 01412 fy1 = d1 * fa[i1]; 01413 fy2 = d2 * fa[i2]; 01414 fy3 = d3 * fa[i3]; 01415 *pyval = fy1 * (xval - x2) * (xval - x3) + 01416 fy2 * (xval - x1) * (xval - x3) + 01417 fy3 * (xval - x1) * (xval - x2); 01418 return 0; 01419 } 01420 01421 01422 /*! 01423 * numaInterpolateArbxVal() 01424 * 01425 * Input: nax (numa of abscissa values) 01426 * nay (numa of ordinate values, corresponding to nax) 01427 * type (L_LINEAR_INTERP, L_QUADRATIC_INTERP) 01428 * xval 01429 * &yval (<return> interpolated value) 01430 * Return: 0 if OK, 1 on error (e.g., if xval is outside range) 01431 * 01432 * Notes: 01433 * (1) The values in nax must be sorted in increasing order. 01434 * If, additionally, they are equally spaced, you can use 01435 * numaInterpolateEqxVal(). 01436 * (2) Caller should check for valid return. 01437 * (3) Uses lagrangian interpolation. See numaInterpolateEqxVal() 01438 * for formulas. 01439 */ 01440 l_int32 01441 numaInterpolateArbxVal(NUMA *nax, 01442 NUMA *nay, 01443 l_int32 type, 01444 l_float32 xval, 01445 l_float32 *pyval) 01446 { 01447 l_int32 i, im, nx, ny, i1, i2, i3; 01448 l_float32 delu, dell, fract, d1, d2, d3; 01449 l_float32 minx, maxx; 01450 l_float32 *fax, *fay; 01451 01452 PROCNAME("numaInterpolateArbxVal"); 01453 01454 if (!pyval) 01455 return ERROR_INT("&yval not defined", procName, 1); 01456 *pyval = 0.0; 01457 if (!nax) 01458 return ERROR_INT("nax not defined", procName, 1); 01459 if (!nay) 01460 return ERROR_INT("nay not defined", procName, 1); 01461 if (type != L_LINEAR_INTERP && type != L_QUADRATIC_INTERP) 01462 return ERROR_INT("invalid interp type", procName, 1); 01463 ny = numaGetCount(nay); 01464 nx = numaGetCount(nax); 01465 if (nx != ny) 01466 return ERROR_INT("nax and nay not same size arrays", procName, 1); 01467 if (ny < 2) 01468 return ERROR_INT("not enough points", procName, 1); 01469 if (type == L_QUADRATIC_INTERP && ny == 2) { 01470 type = L_LINEAR_INTERP; 01471 L_WARNING("only 2 points; using linear interp", procName); 01472 } 01473 numaGetFValue(nax, 0, &minx); 01474 numaGetFValue(nax, nx - 1, &maxx); 01475 if (xval < minx || xval > maxx) 01476 return ERROR_INT("xval is out of bounds", procName, 1); 01477 01478 fax = numaGetFArray(nax, L_NOCOPY); 01479 fay = numaGetFArray(nay, L_NOCOPY); 01480 01481 /* Linear search for interval. We are guaranteed 01482 * to either return or break out of the loop. 01483 * In addition, we are assured that fax[i] - fax[im] > 0.0 */ 01484 if (xval == fax[0]) { 01485 *pyval = fay[0]; 01486 return 0; 01487 } 01488 for (i = 1; i < nx; i++) { 01489 delu = fax[i] - xval; 01490 if (delu >= 0.0) { /* we've passed it */ 01491 if (delu == 0.0) { 01492 *pyval = fay[i]; 01493 return 0; 01494 } 01495 im = i - 1; 01496 dell = xval - fax[im]; /* >= 0 */ 01497 break; 01498 } 01499 } 01500 fract = dell / (fax[i] - fax[im]); 01501 01502 if (type == L_LINEAR_INTERP) { 01503 *pyval = fay[i] + fract * (fay[i + 1] - fay[i]); 01504 return 0; 01505 } 01506 01507 /* Quadratic interpolation */ 01508 if (im == 0) { 01509 i1 = im; 01510 i2 = im + 1; 01511 i3 = im + 2; 01512 } 01513 else { 01514 i1 = im - 1; 01515 i2 = im; 01516 i3 = im + 1; 01517 } 01518 d1 = (fax[i1] - fax[i2]) * (fax[i1] - fax[i3]); 01519 d2 = (fax[i2] - fax[i1]) * (fax[i2] - fax[i3]); 01520 d3 = (fax[i3] - fax[i1]) * (fax[i3] - fax[i2]); 01521 *pyval = fay[i1] * (xval - fax[i2]) * (xval - fax[i3]) / d1 + 01522 fay[i2] * (xval - fax[i1]) * (xval - fax[i3]) / d2 + 01523 fay[i3] * (xval - fax[i1]) * (xval - fax[i2]) / d3; 01524 return 0; 01525 } 01526 01527 01528 /*! 01529 * numaInterpolateEqxInterval() 01530 * 01531 * Input: startx (xval corresponding to first element in nas) 01532 * deltax (x increment between array elements in nas) 01533 * nasy (numa of ordinate values, assumed equally spaced) 01534 * type (L_LINEAR_INTERP, L_QUADRATIC_INTERP) 01535 * x0 (start value of interval) 01536 * x1 (end value of interval) 01537 * npts (number of points to evaluate function in interval) 01538 * &nax (<optional return> array of x values in interval) 01539 * &nay (<return> array of y values in interval) 01540 * Return: 0 if OK, 1 on error 01541 * 01542 * Notes: 01543 * (1) Considering nasy as a function of x, the x values 01544 * are equally spaced. 01545 * (2) This creates nay (and optionally nax) of interpolated 01546 * values over the specified interval (x0, x1). 01547 * (3) If the interval (x0, x1) lies partially outside the array 01548 * nasy (as interpreted by startx and deltax), it is an 01549 * error and returns 1. 01550 * (4) Note that deltax is the intrinsic x-increment for the input 01551 * array nasy, whereas delx is the intrinsic x-increment for the 01552 * output interpolated array nay. 01553 */ 01554 l_int32 01555 numaInterpolateEqxInterval(l_float32 startx, 01556 l_float32 deltax, 01557 NUMA *nasy, 01558 l_int32 type, 01559 l_float32 x0, 01560 l_float32 x1, 01561 l_int32 npts, 01562 NUMA **pnax, 01563 NUMA **pnay) 01564 { 01565 l_int32 i, n; 01566 l_float32 x, yval, maxx, delx; 01567 NUMA *nax, *nay; 01568 01569 PROCNAME("numaInterpolateEqxInterval"); 01570 01571 if (pnax) *pnax = NULL; 01572 if (!pnay) 01573 return ERROR_INT("&nay not defined", procName, 1); 01574 *pnay = NULL; 01575 if (!nasy) 01576 return ERROR_INT("nasy not defined", procName, 1); 01577 if (deltax <= 0.0) 01578 return ERROR_INT("deltax not > 0", procName, 1); 01579 if (type != L_LINEAR_INTERP && type != L_QUADRATIC_INTERP) 01580 return ERROR_INT("invalid interp type", procName, 1); 01581 n = numaGetCount(nasy); 01582 if (type == L_QUADRATIC_INTERP && n == 2) { 01583 type = L_LINEAR_INTERP; 01584 L_WARNING("only 2 points; using linear interp", procName); 01585 } 01586 maxx = startx + deltax * (n - 1); 01587 if (x0 < startx || x1 > maxx || x1 <= x0) 01588 return ERROR_INT("[x0 ... x1] is not valid", procName, 1); 01589 if (npts < 3) 01590 return ERROR_INT("npts < 3", procName, 1); 01591 delx = (x1 - x0) / (l_float32)(npts - 1); /* delx is for output nay */ 01592 01593 if ((nay = numaCreate(npts)) == NULL) 01594 return ERROR_INT("nay not made", procName, 1); 01595 numaSetXParameters(nay, x0, delx); 01596 *pnay = nay; 01597 if (pnax) { 01598 nax = numaCreate(npts); 01599 *pnax = nax; 01600 } 01601 01602 for (i = 0; i < npts; i++) { 01603 x = x0 + i * delx; 01604 if (pnax) 01605 numaAddNumber(nax, x); 01606 numaInterpolateEqxVal(startx, deltax, nasy, type, x, &yval); 01607 numaAddNumber(nay, yval); 01608 } 01609 01610 return 0; 01611 } 01612 01613 01614 /*! 01615 * numaInterpolateArbxInterval() 01616 * 01617 * Input: nax (numa of abscissa values) 01618 * nay (numa of ordinate values, corresponding to nax) 01619 * type (L_LINEAR_INTERP, L_QUADRATIC_INTERP) 01620 * x0 (start value of interval) 01621 * x1 (end value of interval) 01622 * npts (number of points to evaluate function in interval) 01623 * &nadx (<optional return> array of x values in interval) 01624 * &nady (<return> array of y values in interval) 01625 * Return: 0 if OK, 1 on error (e.g., if x0 or x1 is outside range) 01626 * 01627 * Notes: 01628 * (1) The values in nax must be sorted in increasing order. 01629 * If they are not sorted, we do it here, and complain. 01630 * (2) If the values in nax are equally spaced, you can use 01631 * numaInterpolateEqxInterval(). 01632 * (3) Caller should check for valid return. 01633 * (4) We don't call numaInterpolateArbxVal() for each output 01634 * point, because that requires an O(n) search for 01635 * each point. Instead, we do a single O(n) pass through 01636 * nax, saving the indices to be used for each output yval. 01637 * (5) Uses lagrangian interpolation. See numaInterpolateEqxVal() 01638 * for formulas. 01639 */ 01640 l_int32 01641 numaInterpolateArbxInterval(NUMA *nax, 01642 NUMA *nay, 01643 l_int32 type, 01644 l_float32 x0, 01645 l_float32 x1, 01646 l_int32 npts, 01647 NUMA **pnadx, 01648 NUMA **pnady) 01649 { 01650 l_int32 i, im, j, nx, ny, i1, i2, i3, sorted; 01651 l_int32 *index; 01652 l_float32 del, xval, yval, excess, fract, minx, maxx, d1, d2, d3; 01653 l_float32 *fax, *fay; 01654 NUMA *nasx, *nasy, *nadx, *nady; 01655 01656 PROCNAME("numaInterpolateArbxInterval"); 01657 01658 if (pnadx) *pnadx = NULL; 01659 if (!pnady) 01660 return ERROR_INT("&nady not defined", procName, 1); 01661 *pnady = NULL; 01662 if (!nay) 01663 return ERROR_INT("nay not defined", procName, 1); 01664 if (!nax) 01665 return ERROR_INT("nax not defined", procName, 1); 01666 if (type != L_LINEAR_INTERP && type != L_QUADRATIC_INTERP) 01667 return ERROR_INT("invalid interp type", procName, 1); 01668 if (x0 > x1) 01669 return ERROR_INT("x0 > x1", procName, 1); 01670 ny = numaGetCount(nay); 01671 nx = numaGetCount(nax); 01672 if (nx != ny) 01673 return ERROR_INT("nax and nay not same size arrays", procName, 1); 01674 if (ny < 2) 01675 return ERROR_INT("not enough points", procName, 1); 01676 if (type == L_QUADRATIC_INTERP && ny == 2) { 01677 type = L_LINEAR_INTERP; 01678 L_WARNING("only 2 points; using linear interp", procName); 01679 } 01680 numaGetMin(nax, &minx, NULL); 01681 numaGetMax(nax, &maxx, NULL); 01682 if (x0 < minx || x1 > maxx) 01683 return ERROR_INT("xval is out of bounds", procName, 1); 01684 01685 /* Make sure that nax is sorted in increasing order */ 01686 numaIsSorted(nax, L_SORT_INCREASING, &sorted); 01687 if (!sorted) { 01688 L_WARNING("we are sorting nax in increasing order", procName); 01689 numaSortPair(nax, nay, L_SORT_INCREASING, &nasx, &nasy); 01690 } 01691 else { 01692 nasx = numaClone(nax); 01693 nasy = numaClone(nay); 01694 } 01695 01696 fax = numaGetFArray(nasx, L_NOCOPY); 01697 fay = numaGetFArray(nasy, L_NOCOPY); 01698 01699 /* Get array of indices into fax for interpolated locations */ 01700 if ((index = (l_int32 *)CALLOC(npts, sizeof(l_int32))) == NULL) 01701 return ERROR_INT("ind not made", procName, 1); 01702 del = (x1 - x0) / (npts - 1.0); 01703 for (i = 0, j = 0; j < nx && i < npts; i++) { 01704 xval = x0 + i * del; 01705 while (j < nx - 1 && xval > fax[j]) 01706 j++; 01707 if (xval == fax[j]) 01708 index[i] = L_MIN(j, nx - 1); 01709 else /* the index of fax[] is just below xval */ 01710 index[i] = L_MAX(j - 1, 0); 01711 } 01712 01713 /* For each point to be interpolated, get the y value */ 01714 nady = numaCreate(npts); 01715 *pnady = nady; 01716 if (pnadx) { 01717 nadx = numaCreate(npts); 01718 *pnadx = nadx; 01719 } 01720 for (i = 0; i < npts; i++) { 01721 xval = x0 + i * del; 01722 if (pnadx) 01723 numaAddNumber(nadx, xval); 01724 im = index[i]; 01725 excess = xval - fax[im]; 01726 if (excess == 0.0) { 01727 numaAddNumber(nady, fay[im]); 01728 continue; 01729 } 01730 fract = excess / (fax[im + 1] - fax[im]); 01731 01732 if (type == L_LINEAR_INTERP) { 01733 yval = fay[im] + fract * (fay[im + 1] - fay[im]); 01734 numaAddNumber(nady, yval); 01735 continue; 01736 } 01737 01738 /* Quadratic interpolation */ 01739 if (im == 0) { 01740 i1 = im; 01741 i2 = im + 1; 01742 i3 = im + 2; 01743 } 01744 else { 01745 i1 = im - 1; 01746 i2 = im; 01747 i3 = im + 1; 01748 } 01749 d1 = (fax[i1] - fax[i2]) * (fax[i1] - fax[i3]); 01750 d2 = (fax[i2] - fax[i1]) * (fax[i2] - fax[i3]); 01751 d3 = (fax[i3] - fax[i1]) * (fax[i3] - fax[i2]); 01752 yval = fay[i1] * (xval - fax[i2]) * (xval - fax[i3]) / d1 + 01753 fay[i2] * (xval - fax[i1]) * (xval - fax[i3]) / d2 + 01754 fay[i3] * (xval - fax[i1]) * (xval - fax[i2]) / d3; 01755 numaAddNumber(nady, yval); 01756 } 01757 01758 FREE(index); 01759 numaDestroy(&nasx); 01760 numaDestroy(&nasy); 01761 return 0; 01762 } 01763 01764 01765 /*----------------------------------------------------------------------* 01766 * Functions requiring interpolation * 01767 *----------------------------------------------------------------------*/ 01768 /*! 01769 * numaFitMax() 01770 * 01771 * Input: na (numa of ordinate values, to fit a max to) 01772 * &maxval (<return> max value) 01773 * naloc (<optional> associated numa of abscissa values) 01774 * &maxloc (<return> abscissa value that gives max value in na; 01775 * if naloc == null, this is given as an interpolated 01776 * index value) 01777 * Return: 0 if OK; 1 on error 01778 * 01779 * Note: if naloc is given, there is no requirement that the 01780 * data points are evenly spaced. Lagrangian interpolation 01781 * handles that. The only requirement is that the 01782 * data points are ordered so that the values in naloc 01783 * are either increasing or decreasing. We test to make 01784 * sure that the sizes of na and naloc are equal, and it 01785 * is assumed that the correspondences na[i] as a function 01786 * of naloc[i] are properly arranged for all i. 01787 * 01788 * The formula for Lagrangian interpolation through 3 data pts is: 01789 * y(x) = y1(x-x2)(x-x3)/((x1-x2)(x1-x3)) + 01790 * y2(x-x1)(x-x3)/((x2-x1)(x2-x3)) + 01791 * y3(x-x1)(x-x2)/((x3-x1)(x3-x2)) 01792 * 01793 * Then the derivative, using the constants (c1,c2,c3) defined below, 01794 * is set to 0: 01795 * y'(x) = 2x(c1+c2+c3) - c1(x2+x3) - c2(x1+x3) - c3(x1+x2) = 0 01796 */ 01797 l_int32 01798 numaFitMax(NUMA *na, 01799 l_float32 *pmaxval, 01800 NUMA *naloc, 01801 l_float32 *pmaxloc) 01802 { 01803 l_float32 val; 01804 l_float32 smaxval; /* start value of maximum sample, before interpolating */ 01805 l_int32 n, imaxloc; 01806 l_float32 x1, x2, x3, y1, y2, y3, c1, c2, c3, a, b, xmax, ymax; 01807 01808 PROCNAME("numaFitMax"); 01809 01810 *pmaxval = *pmaxloc = 0.0; /* init */ 01811 01812 if (!na) 01813 return ERROR_INT("na not defined", procName, 1); 01814 if (!pmaxval) 01815 return ERROR_INT("&maxval not defined", procName, 1); 01816 if (!pmaxloc) 01817 return ERROR_INT("&maxloc not defined", procName, 1); 01818 01819 n = numaGetCount(na); 01820 if (naloc) { 01821 if (n != numaGetCount(naloc)) 01822 return ERROR_INT("na and naloc of unequal size", procName, 1); 01823 } 01824 01825 numaGetMax(na, &smaxval, &imaxloc); 01826 01827 /* Simple case: max is at end point */ 01828 if (imaxloc == 0 || imaxloc == n - 1) { 01829 *pmaxval = smaxval; 01830 if (naloc) { 01831 numaGetFValue(naloc, imaxloc, &val); 01832 *pmaxloc = val; 01833 } 01834 else 01835 *pmaxloc = imaxloc; 01836 return 0; 01837 } 01838 01839 /* Interior point; use quadratic interpolation */ 01840 y2 = smaxval; 01841 numaGetFValue(na, imaxloc - 1, &val); 01842 y1 = val; 01843 numaGetFValue(na, imaxloc + 1, &val); 01844 y3 = val; 01845 if (naloc) { 01846 numaGetFValue(naloc, imaxloc - 1, &val); 01847 x1 = val; 01848 numaGetFValue(naloc, imaxloc, &val); 01849 x2 = val; 01850 numaGetFValue(naloc, imaxloc + 1, &val); 01851 x3 = val; 01852 } 01853 else { 01854 x1 = imaxloc - 1; 01855 x2 = imaxloc; 01856 x3 = imaxloc + 1; 01857 } 01858 01859 /* Can't interpolate; just use the max val in na 01860 * and the corresponding one in naloc */ 01861 if (x1 == x2 || x1 == x3 || x2 == x3) { 01862 *pmaxval = y2; 01863 *pmaxloc = x2; 01864 return 0; 01865 } 01866 01867 /* Use lagrangian interpolation; set dy/dx = 0 */ 01868 c1 = y1 / ((x1 - x2) * (x1 - x3)); 01869 c2 = y2 / ((x2 - x1) * (x2 - x3)); 01870 c3 = y3 / ((x3 - x1) * (x3 - x2)); 01871 a = c1 + c2 + c3; 01872 b = c1 * (x2 + x3) + c2 * (x1 + x3) + c3 * (x1 + x2); 01873 xmax = b / (2 * a); 01874 ymax = c1 * (xmax - x2) * (xmax - x3) + 01875 c2 * (xmax - x1) * (xmax - x3) + 01876 c3 * (xmax - x1) * (xmax - x2); 01877 *pmaxval = ymax; 01878 *pmaxloc = xmax; 01879 01880 return 0; 01881 } 01882 01883 01884 /*! 01885 * numaDifferentiateInterval() 01886 * 01887 * Input: nax (numa of abscissa values) 01888 * nay (numa of ordinate values, corresponding to nax) 01889 * x0 (start value of interval) 01890 * x1 (end value of interval) 01891 * npts (number of points to evaluate function in interval) 01892 * &nadx (<optional return> array of x values in interval) 01893 * &nady (<return> array of derivatives in interval) 01894 * Return: 0 if OK, 1 on error (e.g., if x0 or x1 is outside range) 01895 * 01896 * Notes: 01897 * (1) The values in nax must be sorted in increasing order. 01898 * If they are not sorted, it is done in the interpolation 01899 * step, and a warning is issued. 01900 * (2) Caller should check for valid return. 01901 */ 01902 l_int32 01903 numaDifferentiateInterval(NUMA *nax, 01904 NUMA *nay, 01905 l_float32 x0, 01906 l_float32 x1, 01907 l_int32 npts, 01908 NUMA **pnadx, 01909 NUMA **pnady) 01910 { 01911 l_int32 i, nx, ny; 01912 l_float32 minx, maxx, der, invdel; 01913 l_float32 *fay; 01914 NUMA *nady, *naiy; 01915 01916 PROCNAME("numaDifferentiateInterval"); 01917 01918 if (pnadx) *pnadx = NULL; 01919 if (!pnady) 01920 return ERROR_INT("&nady not defined", procName, 1); 01921 *pnady = NULL; 01922 if (!nay) 01923 return ERROR_INT("nay not defined", procName, 1); 01924 if (!nax) 01925 return ERROR_INT("nax not defined", procName, 1); 01926 if (x0 > x1) 01927 return ERROR_INT("x0 > x1", procName, 1); 01928 ny = numaGetCount(nay); 01929 nx = numaGetCount(nax); 01930 if (nx != ny) 01931 return ERROR_INT("nax and nay not same size arrays", procName, 1); 01932 if (ny < 2) 01933 return ERROR_INT("not enough points", procName, 1); 01934 numaGetMin(nax, &minx, NULL); 01935 numaGetMax(nax, &maxx, NULL); 01936 if (x0 < minx || x1 > maxx) 01937 return ERROR_INT("xval is out of bounds", procName, 1); 01938 if (npts < 2) 01939 return ERROR_INT("npts < 2", procName, 1); 01940 01941 /* Generate interpolated array over specified interval */ 01942 if (numaInterpolateArbxInterval(nax, nay, L_LINEAR_INTERP, x0, x1, 01943 npts, pnadx, &naiy)) 01944 return ERROR_INT("interpolation failed", procName, 1); 01945 01946 nady = numaCreate(npts); 01947 *pnady = nady; 01948 invdel = 0.5 * ((l_float32)npts - 1.0) / (x1 - x0); 01949 fay = numaGetFArray(naiy, L_NOCOPY); 01950 01951 /* Compute and save derivatives */ 01952 der = 0.5 * invdel * (fay[1] - fay[0]); 01953 numaAddNumber(nady, der); 01954 for (i = 1; i < npts - 1; i++) { 01955 der = invdel * (fay[i + 1] - fay[i - 1]); 01956 numaAddNumber(nady, der); 01957 } 01958 der = 0.5 * invdel * (fay[npts - 1] - fay[npts - 2]); 01959 numaAddNumber(nady, der); 01960 01961 numaDestroy(&naiy); 01962 return 0; 01963 } 01964 01965 01966 /*! 01967 * numaIntegrateInterval() 01968 * 01969 * Input: nax (numa of abscissa values) 01970 * nay (numa of ordinate values, corresponding to nax) 01971 * x0 (start value of interval) 01972 * x1 (end value of interval) 01973 * npts (number of points to evaluate function in interval) 01974 * &sum (<return> integral of function over interval) 01975 * Return: 0 if OK, 1 on error (e.g., if x0 or x1 is outside range) 01976 * 01977 * Notes: 01978 * (1) The values in nax must be sorted in increasing order. 01979 * If they are not sorted, it is done in the interpolation 01980 * step, and a warning is issued. 01981 * (2) Caller should check for valid return. 01982 */ 01983 l_int32 01984 numaIntegrateInterval(NUMA *nax, 01985 NUMA *nay, 01986 l_float32 x0, 01987 l_float32 x1, 01988 l_int32 npts, 01989 l_float32 *psum) 01990 { 01991 l_int32 i, nx, ny; 01992 l_float32 minx, maxx, sum, del; 01993 l_float32 *fay; 01994 NUMA *naiy; 01995 01996 PROCNAME("numaIntegrateInterval"); 01997 01998 if (!psum) 01999 return ERROR_INT("&sum not defined", procName, 1); 02000 *psum = 0.0; 02001 if (!nay) 02002 return ERROR_INT("nay not defined", procName, 1); 02003 if (!nax) 02004 return ERROR_INT("nax not defined", procName, 1); 02005 if (x0 > x1) 02006 return ERROR_INT("x0 > x1", procName, 1); 02007 if (npts < 2) 02008 return ERROR_INT("npts < 2", procName, 1); 02009 ny = numaGetCount(nay); 02010 nx = numaGetCount(nax); 02011 if (nx != ny) 02012 return ERROR_INT("nax and nay not same size arrays", procName, 1); 02013 if (ny < 2) 02014 return ERROR_INT("not enough points", procName, 1); 02015 numaGetMin(nax, &minx, NULL); 02016 numaGetMax(nax, &maxx, NULL); 02017 if (x0 < minx || x1 > maxx) 02018 return ERROR_INT("xval is out of bounds", procName, 1); 02019 02020 /* Generate interpolated array over specified interval */ 02021 if (numaInterpolateArbxInterval(nax, nay, L_LINEAR_INTERP, x0, x1, 02022 npts, NULL, &naiy)) 02023 return ERROR_INT("interpolation failed", procName, 1); 02024 02025 del = (x1 - x0) / ((l_float32)npts - 1.0); 02026 fay = numaGetFArray(naiy, L_NOCOPY); 02027 02028 /* Compute integral (simple trapezoid) */ 02029 sum = 0.5 * (fay[0] + fay[npts - 1]); 02030 for (i = 1; i < npts - 1; i++) 02031 sum += fay[i]; 02032 *psum = del * sum; 02033 02034 numaDestroy(&naiy); 02035 return 0; 02036 } 02037 02038 02039 /*----------------------------------------------------------------------* 02040 * Sorting * 02041 *----------------------------------------------------------------------*/ 02042 /*! 02043 * numaSort() 02044 * 02045 * Input: naout (output numa; can be NULL or equal to nain) 02046 * nain (input numa) 02047 * sortorder (L_SORT_INCREASING or L_SORT_DECREASING) 02048 * Return: naout (output sorted numa), or null on error 02049 * 02050 * Notes: 02051 * (1) Set naout = nain for in-place; otherwise, set naout = NULL. 02052 * (2) Source: Shell sort, modified from K&R, 2nd edition, p.62. 02053 * Slow but simple O(n logn) sort. 02054 */ 02055 NUMA * 02056 numaSort(NUMA *naout, 02057 NUMA *nain, 02058 l_int32 sortorder) 02059 { 02060 l_int32 i, n, gap, j; 02061 l_float32 tmp; 02062 l_float32 *array; 02063 02064 PROCNAME("numaSort"); 02065 02066 if (!nain) 02067 return (NUMA *)ERROR_PTR("nain not defined", procName, NULL); 02068 02069 /* Make naout if necessary; otherwise do in-place */ 02070 if (!naout) 02071 naout = numaCopy(nain); 02072 else if (nain != naout) 02073 return (NUMA *)ERROR_PTR("invalid: not in-place", procName, NULL); 02074 array = naout->array; /* operate directly on the array */ 02075 n = numaGetCount(naout); 02076 02077 /* Shell sort */ 02078 for (gap = n/2; gap > 0; gap = gap / 2) { 02079 for (i = gap; i < n; i++) { 02080 for (j = i - gap; j >= 0; j -= gap) { 02081 if ((sortorder == L_SORT_INCREASING && 02082 array[j] > array[j + gap]) || 02083 (sortorder == L_SORT_DECREASING && 02084 array[j] < array[j + gap])) 02085 { 02086 tmp = array[j]; 02087 array[j] = array[j + gap]; 02088 array[j + gap] = tmp; 02089 } 02090 } 02091 } 02092 } 02093 02094 return naout; 02095 } 02096 02097 02098 /*! 02099 * numaGetSortIndex() 02100 * 02101 * Input: na 02102 * sortorder (L_SORT_INCREASING or L_SORT_DECREASING) 02103 * Return: na giving an array of indices that would sort 02104 * the input array, or null on error 02105 */ 02106 NUMA * 02107 numaGetSortIndex(NUMA *na, 02108 l_int32 sortorder) 02109 { 02110 l_int32 i, n, gap, j; 02111 l_float32 tmp; 02112 l_float32 *array; /* copy of input array */ 02113 l_float32 *iarray; /* array of indices */ 02114 NUMA *naisort; 02115 02116 PROCNAME("numaGetSortIndex"); 02117 02118 if (!na) 02119 return (NUMA *)ERROR_PTR("na not defined", procName, NULL); 02120 if (sortorder != L_SORT_INCREASING && sortorder != L_SORT_DECREASING) 02121 return (NUMA *)ERROR_PTR("invalid sortorder", procName, NULL); 02122 02123 n = numaGetCount(na); 02124 if ((array = numaGetFArray(na, L_COPY)) == NULL) 02125 return (NUMA *)ERROR_PTR("array not made", procName, NULL); 02126 if ((iarray = (l_float32 *)CALLOC(n, sizeof(l_float32))) == NULL) 02127 return (NUMA *)ERROR_PTR("iarray not made", procName, NULL); 02128 for (i = 0; i < n; i++) 02129 iarray[i] = i; 02130 02131 /* Shell sort */ 02132 for (gap = n/2; gap > 0; gap = gap / 2) { 02133 for (i = gap; i < n; i++) { 02134 for (j = i - gap; j >= 0; j -= gap) { 02135 if ((sortorder == L_SORT_INCREASING && 02136 array[j] > array[j + gap]) || 02137 (sortorder == L_SORT_DECREASING && 02138 array[j] < array[j + gap])) 02139 { 02140 tmp = array[j]; 02141 array[j] = array[j + gap]; 02142 array[j + gap] = tmp; 02143 tmp = iarray[j]; 02144 iarray[j] = iarray[j + gap]; 02145 iarray[j + gap] = tmp; 02146 } 02147 } 02148 } 02149 } 02150 02151 naisort = numaCreate(n); 02152 for (i = 0; i < n; i++) 02153 numaAddNumber(naisort, iarray[i]); 02154 02155 FREE(array); 02156 FREE(iarray); 02157 return naisort; 02158 } 02159 02160 02161 /*! 02162 * numaSortByIndex() 02163 * 02164 * Input: nas 02165 * naindex (na that maps from the new numa to the input numa) 02166 * Return: nad (sorted), or null on error 02167 */ 02168 NUMA * 02169 numaSortByIndex(NUMA *nas, 02170 NUMA *naindex) 02171 { 02172 l_int32 i, n, index; 02173 l_float32 val; 02174 NUMA *nad; 02175 02176 PROCNAME("numaSortByIndex"); 02177 02178 if (!nas) 02179 return (NUMA *)ERROR_PTR("nas not defined", procName, NULL); 02180 if (!naindex) 02181 return (NUMA *)ERROR_PTR("naindex not defined", procName, NULL); 02182 02183 n = numaGetCount(nas); 02184 nad = numaCreate(n); 02185 for (i = 0; i < n; i++) { 02186 numaGetIValue(naindex, i, &index); 02187 numaGetFValue(nas, index, &val); 02188 numaAddNumber(nad, val); 02189 } 02190 02191 return nad; 02192 } 02193 02194 02195 /*! 02196 * numaIsSorted() 02197 * 02198 * Input: nas 02199 * sortorder (L_SORT_INCREASING or L_SORT_DECREASING) 02200 * &sorted (<return> 1 if sorted; 0 if not) 02201 * Return: 1 if OK; 0 on error 02202 * 02203 * Notes: 02204 * (1) This is a quick O(n) test if nas is sorted. It is useful 02205 * in situations where the array is likely to be already 02206 * sorted, and a sort operation can be avoided. 02207 */ 02208 l_int32 02209 numaIsSorted(NUMA *nas, 02210 l_int32 sortorder, 02211 l_int32 *psorted) 02212 { 02213 l_int32 i, n; 02214 l_float32 preval, val; 02215 02216 PROCNAME("numaIsSorted"); 02217 02218 if (!psorted) 02219 return ERROR_INT("&sorted not defined", procName, 1); 02220 *psorted = FALSE; 02221 if (!nas) 02222 return ERROR_INT("nas not defined", procName, 1); 02223 if (sortorder != L_SORT_INCREASING && sortorder != L_SORT_DECREASING) 02224 return ERROR_INT("invalid sortorder", procName, 1); 02225 02226 n = numaGetCount(nas); 02227 numaGetFValue(nas, 0, &preval); 02228 for (i = 1; i < n; i++) { 02229 numaGetFValue(nas, i, &val); 02230 if ((sortorder == L_SORT_INCREASING && val < preval) || 02231 (sortorder == L_SORT_DECREASING && val > preval)) 02232 return 0; 02233 } 02234 02235 *psorted = TRUE; 02236 return 0; 02237 } 02238 02239 02240 /*! 02241 * numaSortPair() 02242 * 02243 * Input: nax, nay (input arrays) 02244 * sortorder (L_SORT_INCREASING or L_SORT_DECREASING) 02245 * &nasx (<return> sorted) 02246 * &naxy (<return> sorted exactly in order of nasx) 02247 * Return: 0 if OK, 1 on error 02248 * 02249 * Notes: 02250 * (1) This function sorts the two input arrays, nax and nay, 02251 * together, using nax as the key for sorting. 02252 */ 02253 l_int32 02254 numaSortPair(NUMA *nax, 02255 NUMA *nay, 02256 l_int32 sortorder, 02257 NUMA **pnasx, 02258 NUMA **pnasy) 02259 { 02260 l_int32 sorted; 02261 NUMA *naindex; 02262 02263 PROCNAME("numaSortPair"); 02264 02265 if (!pnasx) 02266 return ERROR_INT("&nasx not defined", procName, 1); 02267 if (!pnasy) 02268 return ERROR_INT("&nasy not defined", procName, 1); 02269 *pnasx = *pnasy = NULL; 02270 if (!nax) 02271 return ERROR_INT("nax not defined", procName, 1); 02272 if (!nay) 02273 return ERROR_INT("nay not defined", procName, 1); 02274 if (sortorder != L_SORT_INCREASING && sortorder != L_SORT_DECREASING) 02275 return ERROR_INT("invalid sortorder", procName, 1); 02276 02277 numaIsSorted(nax, sortorder, &sorted); 02278 if (sorted == TRUE) { 02279 *pnasx = numaCopy(nax); 02280 *pnasy = numaCopy(nay); 02281 } 02282 else { 02283 naindex = numaGetSortIndex(nax, sortorder); 02284 *pnasx = numaSortByIndex(nax, naindex); 02285 *pnasy = numaSortByIndex(nay, naindex); 02286 numaDestroy(&naindex); 02287 } 02288 02289 return 0; 02290 } 02291 02292 02293 /*----------------------------------------------------------------------* 02294 * Random permutation * 02295 *----------------------------------------------------------------------*/ 02296 /*! 02297 * numaPseudorandomSequence() 02298 * 02299 * Input: size (of sequence) 02300 * seed (for random number generation) 02301 * Return: na (pseudorandom on {0,...,size - 1}), or null on error 02302 * 02303 * Notes: 02304 * (1) This uses the Durstenfeld shuffle. 02305 * See: http://en.wikipedia.org/wiki/Fisher–Yates_shuffle. 02306 * Result is a pseudorandom permutation of the sequence of integers 02307 * from 0 to size - 1. 02308 */ 02309 NUMA * 02310 numaPseudorandomSequence(l_int32 size, 02311 l_int32 seed) 02312 { 02313 l_int32 i, index, temp; 02314 l_int32 *array; 02315 NUMA *na; 02316 02317 PROCNAME("numaPseudorandomSequence"); 02318 02319 if (size <= 0) 02320 return (NUMA *)ERROR_PTR("size <= 0", procName, NULL); 02321 02322 if ((array = (l_int32 *)CALLOC(size, sizeof(l_int32))) == NULL) 02323 return (NUMA *)ERROR_PTR("array not made", procName, NULL); 02324 for (i = 0; i < size; i++) 02325 array[i] = i; 02326 srand(seed); 02327 for (i = size - 1; i > 0; i--) { 02328 index = (l_int32)((i + 1) * ((l_float64)rand() / (l_float64)RAND_MAX)); 02329 index = L_MIN(index, i); 02330 temp = array[i]; 02331 array[i] = array[index]; 02332 array[index] = temp; 02333 } 02334 02335 na = numaCreateFromIArray(array, size); 02336 FREE(array); 02337 return na; 02338 } 02339 02340 02341 /*! 02342 * numaRandomPermutation() 02343 * 02344 * Input: nas (input array) 02345 * seed (for random number generation) 02346 * Return: nas (randomly shuggled array), or null on error 02347 */ 02348 NUMA * 02349 numaRandomPermutation(NUMA *nas, 02350 l_int32 seed) 02351 { 02352 l_int32 i, index, size; 02353 l_float32 val; 02354 NUMA *naindex, *nad; 02355 02356 PROCNAME("numaRandomPermutation"); 02357 02358 if (!nas) 02359 return (NUMA *)ERROR_PTR("nas not defined", procName, NULL); 02360 02361 size = numaGetCount(nas); 02362 naindex = numaPseudorandomSequence(size, seed); 02363 nad = numaCreate(size); 02364 for (i = 0; i < size; i++) { 02365 numaGetIValue(naindex, i, &index); 02366 numaGetFValue(nas, index, &val); 02367 numaAddNumber(nad, val); 02368 } 02369 02370 numaDestroy(&naindex); 02371 return nad; 02372 } 02373 02374 02375 /*----------------------------------------------------------------------* 02376 * Functions requiring sorting * 02377 *----------------------------------------------------------------------*/ 02378 /*! 02379 * numaGetRankValue() 02380 * 02381 * Input: na 02382 * fract (use 0.0 for smallest, 1.0 for largest) 02383 * &val (<return> rank val) 02384 * Return: 0 if OK; 1 on error 02385 * 02386 * Notes: 02387 * (1) Computes the rank value of the numbers in the numa, by 02388 * sorting and finding the value a fraction @fract from 02389 * the small end. 02390 */ 02391 l_int32 02392 numaGetRankValue(NUMA *na, 02393 l_float32 fract, 02394 l_float32 *pval) 02395 { 02396 l_int32 n, index; 02397 NUMA *nasort; 02398 02399 PROCNAME("numaGetRankValue"); 02400 02401 if (!pval) 02402 return ERROR_INT("&val not defined", procName, 1); 02403 *pval = 0.0; /* init */ 02404 if (!na) 02405 return ERROR_INT("na not defined", procName, 1); 02406 if (fract < 0.0 || fract > 1.0) 02407 return ERROR_INT("fract not in [0.0 ... 1.0]", procName, 1); 02408 n = numaGetCount(na); 02409 if (n == 0) 02410 return ERROR_INT("na empty", procName, 1); 02411 02412 if ((nasort = numaSort(NULL, na, L_SORT_INCREASING)) == NULL) 02413 return ERROR_INT("nasort not made", procName, 1); 02414 index = (l_int32)(fract * (l_float32)(n - 1) + 0.5); 02415 numaGetFValue(nasort, index, pval); 02416 02417 numaDestroy(&nasort); 02418 return 0; 02419 } 02420 02421 02422 /*! 02423 * numaGetMedian() 02424 * 02425 * Input: na 02426 * &val (<return> median val) 02427 * Return: 0 if OK; 1 on error 02428 * 02429 * Notes: 02430 * (1) Computes the median value of the numbers in the numa, by 02431 * sorting and finding the middle value in the sorted array. 02432 */ 02433 l_int32 02434 numaGetMedian(NUMA *na, 02435 l_float32 *pval) 02436 { 02437 PROCNAME("numaGetMedian"); 02438 02439 if (!pval) 02440 return ERROR_INT("&val not defined", procName, 1); 02441 *pval = 0.0; /* init */ 02442 if (!na) 02443 return ERROR_INT("na not defined", procName, 1); 02444 02445 return numaGetRankValue(na, 0.5, pval); 02446 } 02447 02448 02449 /*! 02450 * numaGetMode() 02451 * 02452 * Input: na 02453 * &val (<return> mode val) 02454 * &count (<optional return> mode count) 02455 * Return: 0 if OK; 1 on error 02456 * 02457 * Notes: 02458 * (1) Computes the mode value of the numbers in the numa, by 02459 * sorting and finding the value of the number with the 02460 * largest count. 02461 * (2) Optionally, also returns that count. 02462 */ 02463 l_int32 02464 numaGetMode(NUMA *na, 02465 l_float32 *pval, 02466 l_int32 *pcount) 02467 { 02468 l_int32 i, n, maxcount, prevcount; 02469 l_float32 val, maxval, prevval; 02470 l_float32 *array; 02471 NUMA *nasort; 02472 02473 PROCNAME("numaGetMode"); 02474 02475 if (!na) 02476 return ERROR_INT("na not defined", procName, 1); 02477 if (!pval) 02478 return ERROR_INT("&val not defined", procName, 1); 02479 02480 *pval = 0.0; 02481 if (pcount) *pcount = 0; 02482 if ((n = numaGetCount(na)) == 0) 02483 return 1; 02484 02485 if ((nasort = numaSort(NULL, na, L_SORT_DECREASING)) == NULL) 02486 return ERROR_INT("nas not made", procName, 1); 02487 array = numaGetFArray(nasort, L_NOCOPY); 02488 02489 /* Initialize with array[0] */ 02490 prevval = array[0]; 02491 prevcount = 1; 02492 maxval = prevval; 02493 maxcount = prevcount; 02494 02495 /* Scan the sorted array, aggregating duplicates */ 02496 for (i = 1; i < n; i++) { 02497 val = array[i]; 02498 if (val == prevval) 02499 prevcount++; 02500 else { /* new value */ 02501 if (prevcount > maxcount) { /* new max */ 02502 maxcount = prevcount; 02503 maxval = prevval; 02504 } 02505 prevval = val; 02506 prevcount = 1; 02507 } 02508 } 02509 02510 /* Was the mode the last run of elements? */ 02511 if (prevcount > maxcount) { 02512 maxcount = prevcount; 02513 maxval = prevval; 02514 } 02515 02516 *pval = maxval; 02517 if (pcount) 02518 *pcount = maxcount; 02519 02520 numaDestroy(&nasort); 02521 return 0; 02522 } 02523 02524 02525 /*----------------------------------------------------------------------* 02526 * Numa combination * 02527 *----------------------------------------------------------------------*/ 02528 /*! 02529 * numaJoin() 02530 * 02531 * Input: nad (dest numa; add to this one) 02532 * nas (<optional> source numa; add from this one) 02533 * istart (starting index in nas) 02534 * iend (ending index in nas; use 0 to cat all) 02535 * Return: 0 if OK, 1 on error 02536 * 02537 * Notes: 02538 * (1) istart < 0 is taken to mean 'read from the start' (istart = 0) 02539 * (2) iend <= 0 means 'read to the end' 02540 * (3) if nas == NULL, this is a no-op 02541 */ 02542 l_int32 02543 numaJoin(NUMA *nad, 02544 NUMA *nas, 02545 l_int32 istart, 02546 l_int32 iend) 02547 { 02548 l_int32 ns, i; 02549 l_float32 val; 02550 02551 PROCNAME("numaJoin"); 02552 02553 if (!nad) 02554 return ERROR_INT("nad not defined", procName, 1); 02555 if (!nas) 02556 return 0; 02557 ns = numaGetCount(nas); 02558 if (istart < 0) 02559 istart = 0; 02560 if (istart >= ns) 02561 return ERROR_INT("istart out of bounds", procName, 1); 02562 if (iend <= 0) 02563 iend = ns - 1; 02564 if (iend >= ns) 02565 return ERROR_INT("iend out of bounds", procName, 1); 02566 if (istart > iend) 02567 return ERROR_INT("istart > iend; nothing to add", procName, 1); 02568 02569 for (i = istart; i <= iend; i++) { 02570 numaGetFValue(nas, i, &val); 02571 numaAddNumber(nad, val); 02572 } 02573 02574 return 0; 02575 } 02576 02577 02578 /*! 02579 * numaaFlattenToNuma() 02580 * 02581 * Input: numaa 02582 * Return: numa, or null on error 02583 * 02584 * Notes: 02585 * (1) This 'flattens' the Numaa to a Numa, by joining successively 02586 * each Numa in the Numaa. 02587 * (2) It doesn't make any assumptions about the location of the 02588 * Numas in the Numaa array, unlike most Numaa functions. 02589 * (3) It leaves the input Numaa unchanged. 02590 */ 02591 NUMA * 02592 numaaFlattenToNuma(NUMAA *naa) 02593 { 02594 l_int32 i, nalloc; 02595 NUMA *na, *nad; 02596 NUMA **array; 02597 02598 PROCNAME("numaaFlattenToNuma"); 02599 02600 if (!naa) 02601 return (NUMA *)ERROR_PTR("naa not defined", procName, NULL); 02602 02603 nalloc = naa->nalloc; 02604 array = numaaGetPtrArray(naa); 02605 nad = numaCreate(0); 02606 for (i = 0; i < nalloc; i++) { 02607 na = array[i]; 02608 if (!na) continue; 02609 numaJoin(nad, na, 0, 0); 02610 } 02611 02612 return nad; 02613 } 02614