IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 8, 2019, 8:21:47 PM (7 years ago)
Author:
eugene
Message:

add reverse option for threshold; median filter for vectors

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Ohana/src/opihi/cmd.data/threshold.c

    r33963 r40623  
    33int threshold (int argc, char **argv) {
    44 
    5   int N, QUIET;
    6   double value;
     5  int N;
    76  Vector *vecx, *vecy;
    87
    9   QUIET = FALSE;
     8  int QUIET = FALSE;
    109  if ((N = get_argument (argc, argv, "-q"))) {
    1110    QUIET = TRUE;
     
    1312  }
    1413
     14  int REVERSE = FALSE;
     15  if ((N = get_argument (argc, argv, "-r"))) {
     16    REVERSE = TRUE;
     17    remove_argument (N, &argc, argv);
     18  }
     19
     20  int BinMin = -1;
     21  int BinMax = -1;
     22  if ((N = get_argument (argc, argv, "-range"))) {
     23    remove_argument (N, &argc, argv);
     24    BinMin = atoi (argv[N]);
     25    remove_argument (N, &argc, argv);
     26    BinMax = atoi (argv[N]);
     27    remove_argument (N, &argc, argv);
     28  }
     29
    1530  if (argc != 4) {
    16     gprint (GP_ERR, "USAGE: threshold <x> <y> (value)\n");
     31    gprint (GP_ERR, "USAGE: threshold <x> <y> (value) [-q] [-r] [-range BinMin BinMax] \n");
    1732    gprint (GP_ERR, "  find the x coordinate at which we pass the specified value\n");
    1833    gprint (GP_ERR, "  by default, y must be monotonically increasing\n");
     34    gprint (GP_ERR, "  -q : quiet mode\n");
     35    gprint (GP_ERR, "  -r : reverse (find downward transition)\n");
    1936    return (FALSE);
    2037  }
     
    2239  if ((vecx = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
    2340  if ((vecy = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
     41  if (vecx[0].Nelements != vecy[0].Nelements) return FALSE;
     42  int Nelements = vecx[0].Nelements;
    2443 
    25   value = atof (argv[3]);
     44  double value = atof (argv[3]);
    2645 
    27   // if (argc == 6) {
    28   //   start = atof (argv[4]);
    29   //   end   = atof (argv[5]);
    30   // } else {
    31   //   start = (vecx[0].type == OPIHI_FLT) ? vecx[0].elements.Flt[0] : vecx[0].elements.Int[0];
    32   //   end   = (vecx[0].type == OPIHI_FLT) ? vecx[0].elements.Flt[vecx[0].Nelements - 1] : vecx[0].elements.Int[vecx[0].Nelements - 1];
    33   // }
    34 
    3546  int isFltX = (vecx[0].type == OPIHI_FLT);
    3647  int isFltY = (vecy[0].type == OPIHI_FLT);
    3748
     49  // BinMin = -1 -> range not provide
     50  if (BinMin == -1) {
     51    BinMin = 0;
     52    BinMax = Nelements - 1;
     53  } else {
     54    if ((BinMin >= vecx[0].Nelements) || (BinMin < -1*vecx[0].Nelements)) {
     55      gprint (GP_ERR, "vector subscript out of range\n");
     56      return FALSE;
     57    }
     58    if (BinMin < 0) BinMin += vecx[0].Nelements;
     59    if ((BinMax >= vecx[0].Nelements) || (BinMax < -1*vecx[0].Nelements)) {
     60      gprint (GP_ERR, "vector subscript out of range\n");
     61      return FALSE;
     62    }
     63    if (BinMax < 0) BinMax += vecx[0].Nelements;
     64  }
     65
    3866  // use bisection to find the value
    39   int Nlo, Nhi;
    40   int Nelements = vecx[0].Nelements;
    4167
    4268  // find the last entry before start
    43   Nlo = 0;
    44   Nhi = Nelements - 1;
    45   while (Nhi - Nlo > 10) {
    46     N = 0.5*(Nlo + Nhi);
    47     double testval = isFltY ? vecy[0].elements.Flt[N] : vecy[0].elements.Int[N];
    48     if (testval < value) {
    49       Nlo = MAX(N, 0);
    50     } else {
    51       Nhi = MIN(N, Nelements - 1);
     69  int Nlo = BinMin;
     70  int Nhi = BinMax;
     71
     72  if (!REVERSE) {
     73    while (Nhi - Nlo > 10) {
     74      N = 0.5*(Nlo + Nhi);
     75      double testval = isFltY ? vecy[0].elements.Flt[N] : vecy[0].elements.Int[N];
     76      if (testval < value) {
     77        Nlo = MAX(N, 0);
     78      } else {
     79        Nhi = MIN(N, Nelements - 1);
     80      }
    5281    }
    53   }
    54   // v[Nlo] < value <= v[Nhi]
    55   for (N = Nlo; N <= Nhi; N++) {
    56     double testval = isFltY ? vecy[0].elements.Flt[N] : vecy[0].elements.Int[N];
    57     if (testval > value) {
    58       Nhi = N;
    59       break;
     82    // v[Nlo] < value <= v[Nhi]
     83    for (N = Nlo; N <= Nhi; N++) {
     84      double testval = isFltY ? vecy[0].elements.Flt[N] : vecy[0].elements.Int[N];
     85      if (testval > value) {
     86        Nhi = N;
     87        break;
     88      }
     89    }
     90  } else {
     91    while (Nhi - Nlo > 10) {
     92      N = 0.5*(Nlo + Nhi);
     93      double testval = isFltY ? vecy[0].elements.Flt[N] : vecy[0].elements.Int[N];
     94      if (testval > value) {
     95        Nlo = MAX(N, 0);
     96      } else {
     97        Nhi = MIN(N, Nelements - 1);
     98      }
     99    }
     100    // v[Nlo] > value >= v[Nhi]
     101    for (N = Nlo; N <= Nhi; N++) {
     102      double testval = isFltY ? vecy[0].elements.Flt[N] : vecy[0].elements.Int[N];
     103      if (testval < value) {
     104        Nhi = N;
     105        break;
     106      }
    60107    }
    61108  }
Note: See TracChangeset for help on using the changeset viewer.