IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 4, 2021, 6:02:11 PM (5 years ago)
Author:
eugene
Message:

merge changes from eam_branches/ipp-dev-20210817 (add chebyshev, add fit2d_full, add inner_fraction to medimage, change dvoImageOverlap test to CERSTD)

Location:
trunk/Ohana
Files:
8 edited
2 copied

Legend:

Unmodified
Added
Removed

  • trunk/Ohana

  • trunk/Ohana/src/opihi/cmd.data/Makefile

    r41432 r41891  
    2828$(SRC)/bisection.$(ARCH).o              \
    2929$(SRC)/cast.$(ARCH).o           \
     30$(SRC)/chebyshev.$(ARCH).o              \
     31$(SRC)/chebyshev_commands.$(ARCH).o     \
    3032$(SRC)/center.$(ARCH).o \
    3133$(SRC)/clear.$(ARCH).o          \

  • trunk/Ohana/src/opihi/cmd.data/applyfit2d.c

    r29001 r41891  
    116116
    117117
     118int applyfit2d_full (int argc, char **argv) {
     119 
     120  int i, j, n, order;
     121  char *c, name[64];
     122  double **C, X, Y;
     123  opihi_flt *z;
     124  Vector *xvec, *yvec, *zvec;
     125
     126  if (argc != 4) {
     127    gprint (GP_ERR, "USAGE: applyfit2d_full x y z\n");
     128    return (FALSE);
     129  }
     130
     131  c = get_variable ("Cnn");
     132  if (c == NULL) {
     133    gprint (GP_ERR, "no fit available\n");
     134    return (FALSE);
     135  }
     136  order = atof (c);
     137  free (c);
     138
     139  if ((xvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);   
     140  if ((yvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);   
     141  if ((zvec = SelectVector (argv[3], ANYVECTOR, TRUE)) == NULL) return (FALSE);   
     142
     143  ALLOCATE (C, double *, order+1);
     144  for (i = 0; i <= order; i++) {
     145    ALLOCATE (C[i], double, order+1);
     146    for (j = 0; j <= order; j++) {
     147      sprintf (name, "CX%dY%d", i, j);
     148      c = get_variable (name);
     149      if (c == NULL) {
     150        gprint (GP_ERR, "missing fit term %d,%d\n", i, j);
     151        for (j = 0; j < i; j++) free (C[j]);
     152        free (C);
     153        return (FALSE);
     154      }
     155      C[i][j] = atof (c);
     156      free (c);
     157    }
     158  }
     159
     160  ResetVector (zvec, OPIHI_FLT, xvec[0].Nelements);
     161  bzero (zvec[0].elements.Flt, sizeof(opihi_flt)*zvec[0].Nelements);
     162  z = zvec[0].elements.Flt;
     163
     164  // we have four cases (x == flt or int) and (y == flt or int)
     165  if ((xvec[0].type == OPIHI_FLT) && (yvec[0].type == OPIHI_FLT)) {
     166    opihi_flt *x = xvec[0].elements.Flt;
     167    opihi_flt *y = yvec[0].elements.Flt;
     168    for (n = 0; n < xvec[0].Nelements; n++, x++, y++, z++) {
     169      Y = X = 1;
     170      for (j = 0; j <= order; j++) {
     171        X = Y;
     172        for (i = 0; i <= order; i++) {
     173          *z += C[i][j]*X;
     174          X = X * (*x);
     175        }
     176        Y = Y * (*y);
     177      }
     178    }
     179  }
     180  if ((xvec[0].type == OPIHI_FLT) && (yvec[0].type == OPIHI_INT)) {
     181    opihi_flt *x = xvec[0].elements.Flt;
     182    opihi_int *y = yvec[0].elements.Int;
     183    for (n = 0; n < xvec[0].Nelements; n++, x++, y++, z++) {
     184      Y = X = 1;
     185      for (j = 0; j <= order; j++) {
     186        X = Y;
     187        for (i = 0; i <= order; i++) {
     188          *z += C[i][j]*X;
     189          X = X * (*x);
     190        }
     191        Y = Y * (*y);
     192      }
     193    }
     194  }
     195  if ((xvec[0].type == OPIHI_INT) && (yvec[0].type == OPIHI_FLT)) {
     196    opihi_int *x = xvec[0].elements.Int;
     197    opihi_flt *y = yvec[0].elements.Flt;
     198    for (n = 0; n < xvec[0].Nelements; n++, x++, y++, z++) {
     199      Y = X = 1;
     200      for (j = 0; j <= order; j++) {
     201        X = Y;
     202        for (i = 0; i <= order; i++) {
     203          *z += C[i][j]*X;
     204          X = X * (*x);
     205        }
     206        Y = Y * (*y);
     207      }
     208    }
     209  }
     210  if ((xvec[0].type == OPIHI_INT) && (yvec[0].type == OPIHI_INT)) {
     211    opihi_int *x = xvec[0].elements.Int;
     212    opihi_int *y = yvec[0].elements.Int;
     213    for (n = 0; n < xvec[0].Nelements; n++, x++, y++, z++) {
     214      Y = X = 1;
     215      for (j = 0; j <= order; j++) {
     216        X = Y;
     217        for (i = 0; i <= order; i++) {
     218          *z += C[i][j]*X;
     219          X = X * (*x);
     220        }
     221        Y = Y * (*y);
     222      }
     223    }
     224  }
     225
     226  for (i = 0; i < order + 1; i++) free (C[i]);
     227  free (C);
     228
     229  return (TRUE);
     230}
     231
     232

  • trunk/Ohana/src/opihi/cmd.data/fit2d.c

    r39020 r41891  
    294294  return (TRUE);
    295295}
     296
     297int fit2d_full (int argc, char **argv) {
     298 
     299  int N;
     300  char name[64];
     301
     302  int VERBOSE = FALSE;
     303  if ((N = get_argument (argc, argv, "-v"))) {
     304    VERBOSE = TRUE;
     305    remove_argument (N, &argc, argv);
     306  }
     307
     308  int Quiet = FALSE;
     309  if ((N = get_argument (argc, argv, "-q"))) {
     310    Quiet = TRUE;
     311    remove_argument (N, &argc, argv);
     312  }
     313  if ((N = get_argument (argc, argv, "-quiet"))) {
     314    Quiet = TRUE;
     315    remove_argument (N, &argc, argv);
     316  }
     317
     318  float ClipNSigma = 0;
     319  int ClipNiter  = 1;
     320  if ((N = get_argument (argc, argv, "-clip"))) {
     321    remove_argument (N, &argc, argv);
     322    ClipNSigma = atof(argv[N]);
     323    remove_argument (N, &argc, argv);
     324    ClipNiter  = atof(argv[N]);
     325    remove_argument (N, &argc, argv);
     326  }
     327
     328  opihi_flt *dz = NULL;
     329  Vector *dzvec = NULL;
     330  int Weight = FALSE;
     331  if ((N = get_argument (argc, argv, "-dz"))) {
     332    remove_argument (N, &argc, argv);
     333    if ((dzvec = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE);   
     334    remove_argument (N, &argc, argv);
     335    Weight = TRUE;
     336  }
     337
     338  if (argc != 5) {
     339    gprint (GP_ERR, "USAGE: fit2d_full x y z order [-dz wt] [-quiet/-q] [-clip Nsigma Niter]\n");
     340    return (FALSE);
     341  }
     342
     343  Vector *xvec, *yvec, *zvec;
     344  if ((xvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);   
     345  if ((yvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);   
     346  if ((zvec = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE);   
     347
     348  if (xvec[0].Nelements != yvec[0].Nelements) {
     349    gprint (GP_ERR, "vectors must have same length\n");
     350    return (FALSE);
     351  }
     352  if (xvec[0].Nelements != zvec[0].Nelements) {
     353    gprint (GP_ERR, "vectors must have same length\n");
     354    return (FALSE);
     355  }
     356  CastVector (xvec, OPIHI_FLT);
     357  CastVector (yvec, OPIHI_FLT);
     358  CastVector (zvec, OPIHI_FLT);
     359
     360  if (Weight) {
     361    CastVector (dzvec, OPIHI_FLT);
     362    if (xvec[0].Nelements != dzvec[0].Nelements) {
     363      gprint (GP_ERR, "vectors must have same length\n");
     364      return (FALSE);
     365    }
     366  }
     367 
     368  int order = atof (argv[4]);
     369  int nterm = SQ(order + 1);
     370
     371  ALLOCATE_PTR (zfit, opihi_flt, xvec[0].Nelements);
     372  ALLOCATE_PTR (mask, char, xvec[0].Nelements);
     373  memset (mask, 0, xvec[0].Nelements);
     374
     375  /* allocate the summation matrices */
     376  ALLOCATE_PTR (b, double *, nterm);
     377  ALLOCATE_PTR (c, double *, nterm);
     378  for (int i = 0; i < nterm; i++) {
     379    ALLOCATE (c[i], double, nterm);
     380    ALLOCATE (b[i], double, 1);
     381  }
     382
     383  for (N = 0; N < ClipNiter; N++) {
     384
     385    /* init registers for current pass */
     386    for (int i = 0; i < nterm; i++) {
     387      memset (c[i], 0, nterm*sizeof(double));
     388      memset (b[i], 0, sizeof(double));
     389    }
     390
     391    opihi_flt *x = xvec[0].elements.Flt;
     392    opihi_flt *y = yvec[0].elements.Flt;
     393    opihi_flt *z = zvec[0].elements.Flt;
     394    if (Weight) dz = dzvec[0].elements.Flt;
     395
     396    /* add up the x,y values */
     397    for (int i = 0; i < xvec[0].Nelements; i++, x++, y++, z++) {
     398      if (mask[i]) continue;
     399      if (!finite(*x) || !finite(*y)) continue;
     400      double dZ = 1.0;
     401      if (Weight) {
     402        dZ = 1.0 / SQ(*dz);
     403        dz ++;
     404      }
     405
     406      int n = 0;
     407      double iY = 1*dZ;
     408      double iX = iY;
     409
     410      for (int iy = 0; iy <= order; iy++) {
     411        iX = iY;
     412        for (int ix = 0; ix <= order; ix++) {
     413          b[n][0] += iX * (*z);
     414
     415          int m = 0;
     416          double jY = 1;
     417          double jX = jY;
     418
     419          for (int jy = 0; jy <= order; jy++) {
     420            jX = jY;
     421            for (int jx = 0; jx <= order; jx++) {
     422              c[n][m] += iX*jX;
     423              jX = jX * (*x);
     424              m ++;
     425            }
     426            jY = jY * (*y);
     427          }
     428          iX = iX * (*x);
     429          n ++;
     430        }
     431        iY = iY * (*y);
     432      }
     433    }
     434
     435    /** test print **/
     436    if (VERBOSE) {
     437      for (int i = 0; i < nterm; i++) {
     438        for (int j = 0; j < nterm; j++) {
     439          gprint (GP_ERR, "%g  ", c[i][j]);
     440        }
     441        gprint (GP_ERR, "\n");
     442      }
     443      gprint (GP_ERR, "-----\n");
     444    }
     445
     446    dgaussjordan (c, b, nterm, 1);
     447
     448    /** test print **/
     449    if (VERBOSE) {
     450      int n = 0;
     451      for (int ny = 0; ny <= order; ny++) {
     452        for (int nx = 0; nx <= order; nx++) {
     453          gprint (GP_ERR, "x^%d y^%d: %g\n", nx, ny, b[n][0]);
     454          n ++;
     455        }
     456      }
     457    }
     458
     459    /** test print **/
     460    if (VERBOSE) {
     461      for (int i = 0; i < nterm; i++) {
     462        for (int j = 0; j < nterm; j++) {
     463          gprint (GP_ERR, "%g  ", c[i][j]);
     464        }
     465        gprint (GP_ERR, "\n");
     466      }
     467      gprint (GP_ERR, "-----\n");
     468    }
     469
     470    /* the b[][0] terms are in the following order:
     471       y^0 x^0, y^0 x^1, ... y^0 x^N
     472       y^1 x^0, y^1 x^1, ... y^1 x^N
     473       ...
     474       y^N x^0, y^N x^1, ... y^N x^N
     475    */
     476
     477    /* generate fitted values */
     478    x = xvec[0].elements.Flt;
     479    y = yvec[0].elements.Flt;
     480    opihi_flt *zf = zfit;
     481    for (int n = 0; n < xvec[0].Nelements; n++, x++, y++, zf++) {
     482      if (!finite(*x) || !finite(*y)) continue;
     483      *zf = 0;
     484      double Y = 1;
     485      double X = Y;
     486
     487      int n = 0;
     488      for (int ny = 0; ny <= order; ny++) {
     489        X = Y;
     490        for (int nx = 0; nx <= order; nx++) {
     491          *zf += b[n][0]*X;
     492          X = X * (*x);
     493          n ++;
     494        }
     495        Y = Y * (*y);
     496      }
     497    }
     498
     499    /* measure fit residual scatter */
     500    x  = xvec[0].elements.Flt;
     501    y  = yvec[0].elements.Flt;
     502    z  = zvec[0].elements.Flt;
     503    zf = zfit;
     504    double dZ = 0.0;
     505    double dZ2 = 0.0;
     506    int Npt = 0;
     507    for (int i = 0, Npt = 0; i < xvec[0].Nelements; i++, x++, y++, z++, zf++) {
     508      if (mask[i]) continue;
     509      if (!finite(*x) || !finite(*y)) continue;
     510      dZ  += (*z - *zf);
     511      dZ2 += SQ(*z - *zf);
     512      Npt ++;
     513    }
     514    double mean  = dZ / Npt;
     515    double sigma = sqrt (fabs(dZ2/Npt - SQ(mean)));
     516    double maxsigma = ClipNSigma * sigma;
     517
     518    if (VERBOSE) gprint (GP_ERR, "mean: %g, sigma: %g, maxsigma: %g\n", mean, sigma, maxsigma);
     519
     520    /* mask outlier points */
     521    x  = xvec[0].elements.Flt;
     522    y  = yvec[0].elements.Flt;
     523    z  = zvec[0].elements.Flt;
     524    zf = zfit;
     525    int Nmask = 0;
     526    for (int i = 0; ClipNSigma && (i < xvec[0].Nelements); i++, x++, y++, z++, zf++) {
     527      dZ = (*z - *zf);
     528      if (fabs(dZ) > maxsigma) {
     529        mask[i] = TRUE;
     530        Nmask ++;
     531      } else {
     532        mask[i] = FALSE;
     533      }
     534    }
     535    if (VERBOSE) gprint (GP_ERR, "pass: %d, Nmask: %d\n", N, Nmask);
     536  }
     537
     538  if (!Quiet) gprint (GP_ERR, "z = ");
     539  int n = 0;
     540  for (int ny = 0; ny <= order; ny++) {
     541    for (int nx = 0; nx <= order; nx++) {
     542      sprintf (name, "CX%dY%d", nx, ny);
     543      set_variable (name, b[n][0]);
     544      if (!Quiet) gprint (GP_ERR, "%f x^%d y^%d  ", b[n][0], nx, ny);
     545      n++;
     546    }
     547  }
     548  sprintf (name, "Cnn");
     549  set_variable (name, (double) order);
     550 
     551  if (!Quiet) gprint (GP_ERR, "\n");
     552
     553  /* free internal data */
     554  free (zfit);
     555  free (mask);
     556
     557  for (int i = 0; i < nterm; i++) {
     558    free (c[i]);
     559    free (b[i]);
     560  }
     561  free (b);
     562  free (c);
     563
     564  return (TRUE);
     565}

  • trunk/Ohana/src/opihi/cmd.data/init.c

    r41432 r41891  
    77int applyfit1d       PROTO((int, char **));
    88int applyfit2d       PROTO((int, char **));
     9int applyfit2d_full  PROTO((int, char **));
    910int applyfit3d       PROTO((int, char **));
    1011int box              PROTO((int, char **));
    … …  
    1314int center           PROTO((int, char **));
    1415int cast             PROTO((int, char **));
     16int chebyshev_command PROTO((int, char **));
    1517int circstats        PROTO((int, char **));
    1618int clear            PROTO((int, char **));
    … …  
    4547int fit1d_irls       PROTO((int, char **));
    4648int fit2d            PROTO((int, char **));
     49int fit2d_full       PROTO((int, char **));
    4750int fit3d            PROTO((int, char **));
    4851int gaussjordan      PROTO((int, char **));
    … …  
    204207  {1, "applyfit1d",   applyfit1d,       "apply 1-d fit to new vector"},
    205208  {1, "applyfit2d",   applyfit2d,       "apply 2-d fit to new vector"},
     209  {1, "applyfit2d_full", applyfit2d_full, "apply 2-d fit to new vector"},
    206210  {1, "applyfit3d",   applyfit3d,       "apply 3-d fit to new vector"},
    207211  {1, "bisection",    bisection,        "use bisection to find threshold in vector"},
    … …  
    211215  {1, "center",       center,           "center image on coords"},
    212216  {1, "cast",         cast,             "cast input vector to specified type"},
     217  {1, "chebyshev",    chebyshev_command, "Chebyshev polynomial commands"},
    213218  {1, "circstats",    circstats,        "circular statistics"},
    214219  {1, "clear",        clear,            "erase plot"},
    … …  
    244249  {1, "fit1d_irls",   fit1d_irls,       "fit polynomial to vector pair"},
    245250  {1, "fit2d",        fit2d,            "fit 2-d polynomial to vector triplet"},
     251  {1, "fit2d_full",   fit2d_full,       "fit 2-d polynomial to vector triplet (all X^n Y^m coeffs)"},
    246252  {1, "fit3d",        fit3d,            "fit 3-d polynomial to vector quad"},
    247253  {1, "gaussdev",     gaussdeviate,     "generate a gaussian deviate vector"},

  • trunk/Ohana/src/opihi/cmd.data/medimage_calc.c

    r41341 r41891  
    1919static float *irls_testval = NULL;
    2020
    21 enum {CALC_MEDIAN, CALC_MEAN, CALC_IRLS, CALC_WTMEAN};
     21enum {CALC_MEDIAN, CALC_MEAN, CALC_IRLS, CALC_WTMEAN, CALC_INNER_FRACTION};
    2222
    2323int medimage_calc (int argc, char **argv) {
    … …  
    3939  }
    4040
     41  float minRange = NAN;
     42  float maxRange = NAN;
     43
    4144  int mode = CALC_MEDIAN;
    4245  if ((N = get_argument (argc, argv, "-mean"))) {
    … …  
    4548  }
    4649  if ((N = get_argument (argc, argv, "-irls"))) {
    47     if (mode != CALC_MEDIAN) { gprint (GP_ERR, "supply only one of -mean, -median, -irls, -wtmean\n"); return FALSE; }
     50    if (mode != CALC_MEDIAN) goto choose_one;
    4851    mode = CALC_IRLS;
    4952    remove_argument (N, &argc, argv);
    5053  }
    5154  if ((N = get_argument (argc, argv, "-wtmean"))) {
    52     if (mode != CALC_MEDIAN) { gprint (GP_ERR, "supply only one of -mean, -median, -irls, -wtmean\n"); return FALSE; }
     55    if (mode != CALC_MEDIAN) goto choose_one;
    5356    mode = CALC_WTMEAN;
    5457    remove_argument (N, &argc, argv);
    5558  }
     59  if ((N = get_argument (argc, argv, "-inner-fraction"))) {
     60    if (mode != CALC_MEDIAN) goto choose_one;
     61    if (argc < N + 3) goto fraction_options;
     62    mode = CALC_INNER_FRACTION;
     63    remove_argument (N, &argc, argv);
     64    minRange = atof(argv[N]);
     65    if ((minRange < 0.0) || (minRange > 1.0)) goto fraction_options;
     66    remove_argument (N, &argc, argv);
     67    maxRange = atof(argv[N]);
     68    if ((maxRange < 0.0) || (maxRange > 1.0)) goto fraction_options;
     69    if (maxRange <= minRange) goto fraction_options;
     70    remove_argument (N, &argc, argv);
     71  }
    5672  if ((N = get_argument (argc, argv, "-median"))) {
    57     if (mode != CALC_MEDIAN) { gprint (GP_ERR, "supply only one of -mean, -median, -irls, -wtmean\n"); return FALSE; }
     73    if (mode != CALC_MEDIAN) goto choose_one;
    5874    mode = CALC_MEDIAN;
    5975    remove_argument (N, &argc, argv);
    … …  
    157173          }
    158174          break;
     175        case CALC_INNER_FRACTION:
     176          fsort (val, N);
     177          int Ns = MIN(MAX(0, minRange * N),N); // e.g., 0.1 * 50 = 5
     178          int Ne = MIN(MAX(0, maxRange * N),N); // e.g., 0.9 * 50 = 45
     179          int Npt = Ne - Ns;
     180
     181          float sum = 0.0;
     182          for (n = Ns; n < Ne; n++) {
     183            sum += val[n];
     184          }
     185          outvalue[Npix] = sum / (float) Npt;
     186
     187          if (varvalue) {
     188            float sum = 0.0;
     189            for (n = Ns; n < Ne; n++) {
     190              sum += SQ(val[n] - outvalue[Npix]);
     191            }
     192            // variance on the mean (stdev / sqrt(N))^2
     193            varvalue[Npix] = sum / (Npt - 1) / Npt;
     194          }
     195          break;
    159196        case CALC_MEAN: {
    160197          float sum = 0.0;
    … …  
    195232  irls_free();
    196233  return TRUE;
     234
     235 choose_one:
     236  gprint (GP_ERR, "supply only one of -mean, -median, -irls, -wtmean, -inner-fraction\n");
     237  return FALSE;
     238
     239 fraction_options:
     240  gprint (GP_ERR, "USE: -inner-fraction (min) (max) : min & max in range 0.0 to 1.0\n");
     241  return FALSE;
    197242}
    198243

  • trunk/Ohana/src/opihi/cmd.data/tdhistogram.c

    r39233 r41891  
    4545  int valid = (!reuse && (argc == 11)) || (reuse && (argc == 5));
    4646  if (!valid) {
    47     gprint (GP_ERR, "USAGE: tdhistogram buffer x y z (Xmin) (Xmax) (Ymin) (Ymax) (Zmin) (Zmax) [-dx dx] [-dy dy] [-dz dz]\n");
     47    gprint (GP_ERR, "USAGE: tdhistogram buffer x y z (Xmin) (Xmax) (Ymin) (Ymax) (Zmin) (Zmax) [-dx dx] [-dy dy] [-dz dz] [-range range]\n");
    4848    gprint (GP_ERR, "   OR: tdhistogram buffer x y z -reuse\n");
     49    gprint (GP_ERR, "   -dx, -dy, -dy specify the bin size in these directions\n");
     50    gprint (GP_ERR, "   -range : generate an output vector corresponding to the elements in the z-direction\n");
    4951    gprint (GP_ERR, " output buffer is 3D\n");
    5052    return (FALSE);

  • trunk/Ohana/src/opihi/cmd.data/vstats.c

    r38441 r41891  
    4848  }
    4949
    50   if (argc != 2) {
    51     gprint (GP_ERR, "USAGE: vstat (vector) [-ignore value] [-q] [-quiet] [-iter Niter] [-sigma-clip Nsigma]\n");
     50  int valid = (argc == 2);
     51  valid |= (argc > 3) && !strcmp (argv[2], "where");
     52  if (!valid) {
     53    gprint (GP_ERR, "USAGE: vstat (vector) [-ignore value] [-q] [-quiet] [-iter Niter] [-sigma-clip Nsigma] [where (logical expression)]\n");
    5254    gprint (GP_ERR, "  default is 1 iteration without sigma clipping or 3 with sigma clipping\n");
    5355    return (FALSE);
    … …  
    5860  /* we need two passes, one for max, min, mean, sum, one for median, stdev, etc */
    5961
    60   // set a good / bad mask
     62  // tvec is used for logical test (truth vector)
     63  Vector *tvec = NULL;
     64  if (argc > 3) {
     65    int size;
     66    char *out = dvomath (argc - 3, &argv[3], &size, 1);
     67    if (out == NULL) {
     68      print_error ();
     69      return FALSE;
     70    }
     71    if ((tvec = SelectVector (out, OLDVECTOR, TRUE)) == NULL) {
     72      gprint (GP_ERR, " invalid logic result\n");
     73      DeleteNamedVector (out);
     74      free (out);
     75      return (FALSE);
     76    }
     77    if (tvec->Nelements != vec->Nelements) {
     78      gprint (GP_ERR, "logical vector does not match in length\n");
     79      return FALSE;
     80    }
     81  }
     82
     83  // set a good / bad mask (mask == 1 means 'bad')
    6184  ALLOCATE (mask, char, vec[0].Nelements);
    6285  if (vec[0].type == OPIHI_FLT) {
    6386    opihi_flt *X = vec[0].elements.Flt;
    6487    for (i = 0; i < vec[0].Nelements; i++, X++) {
    65       mask[i] = 1;
    66       if (!finite (*X)) continue;
    67       if (Ignore && (*X == IgnoreValue)) continue;
    68       mask[i] = 0;
     88      mask[i] = 0; // do not mask unless we have a reason below
     89      if (tvec) {
     90        // note the logical vector is 0 == bad (ignore) while mask[i] has the opposite sense (0 is good)
     91        mask[i] = (tvec->type == OPIHI_FLT) ? (tvec->elements.Flt[i] == 0.0) : (tvec->elements.Int[i] == 0);
     92      }
     93      if (!finite (*X)) { mask[i] = 1; }
     94      if (Ignore && (*X == IgnoreValue)) { mask[i] = 1; }
    6995    }     
    7096  } else {
    7197    opihi_int *X = vec[0].elements.Int;
    7298    for (i = 0; i < vec[0].Nelements; i++, X++) {
    73       mask[i] = 1;
    74       if (!finite (*X)) continue;
    75       if (Ignore && (*X == IgnoreValue)) continue;
    76       mask[i] = 0;
     99      mask[i] = 0; // do not mask unless we have a reason below
     100      if (tvec) {
     101        // note the logical vector is 0 == bad (ignore) while mask[i] has the opposite sense (0 is good)
     102        mask[i] = (tvec->type == OPIHI_FLT) ? (tvec->elements.Flt[i] == 0.0) : (tvec->elements.Int[i] == 0);
     103      }
     104      if (!finite (*X)) { mask[i] = 1; }
     105      if (Ignore && (*X == IgnoreValue)) { mask[i] = 1; }
    77106    }     
    78107  }
Note: See TracChangeset for help on using the changeset viewer.