Leptonica 1.68
C Image Processing Library

numafunc1.c

Go to the documentation of this file.
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 
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Defines