IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 30610


Ignore:
Timestamp:
Feb 13, 2011, 11:20:31 AM (15 years ago)
Author:
eugene
Message:

add +=, -= operators to opihi; limits -by-image; -closest option for match2d; fix reindex function; output precision for write_vectors.c; add sigma-clip and iteration to vstats; fix possible crash in subset (allow 0-length result)

Location:
trunk/Ohana/src/opihi/cmd.data
Files:
6 edited

Legend:

Unmodified
Added
Removed

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

    r29540 r30610  
    2626
    2727  // XXX need an option to set the limits based on the current image bounds
     28  if ((N = get_argument (argc, argv, "-image"))) {
     29    remove_argument (N, &argc, argv);
     30    KapaGetImageRange (kapa, &graphmode.xmin, &graphmode.xmax, &graphmode.ymax, &graphmode.ymin);
     31
     32    set_variable ("XMIN", graphmode.xmin);
     33    set_variable ("XMAX", graphmode.xmax);
     34    set_variable ("YMIN", graphmode.ymin);
     35    set_variable ("YMAX", graphmode.ymax);
     36
     37    // if (!NoClear) KapaClearSections (kapa);
     38    SetGraph (&graphmode);
     39    KapaSetLimits (kapa, &graphmode);
     40    return (TRUE);
     41    // Set Region based on image
     42  }
    2843
    2944  if (argc == 1) {

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

    r27817 r30610  
    11# include "data.h"
     2
    23int find_matches2d (Vector *X1, Vector *Y1, Vector *X2, Vector *Y2, double Radius, Vector *index1, Vector *index2);
     4int find_matches2d_closest (Vector *X1, Vector *Y1, Vector *X2, Vector *Y2, double Radius, Vector *index);
    35
    46// match2d (X1) (Y1) (X2) (Y2) (Radius) [-index1 (index1)] [-index2 (index2)] [-nomatch1 nomatch1] [-nomatch2 nomatch2]
    … …  
    68int match2d (int argc, char **argv) {
    79 
    8   int N;
     10  int N, CLOSEST;
    911  double Radius;
    1012  char *endptr;
    … …  
    1214  Vector *index1, *index2;
    1315
     16  CLOSEST = FALSE;
     17  if ((N = get_argument (argc, argv, "-closest"))) {
     18    remove_argument (N, &argc, argv);
     19    CLOSEST = TRUE;
     20  }
     21
    1422  if ((N = get_argument (argc, argv, "-index1"))) {
    1523    remove_argument (N, &argc, argv);
    … …  
    3038  if (argc != 6) {
    3139    gprint (GP_ERR, "USAGE: match2d X1 Y1 X2 Y2 Radius [-index1 (index1)] [-index2 (index2)] [-closest]\n");
    32     // gprint (GP_ERR, "USAGE: match2d X1 Y1 X2 Y2 Radius [-index1 (index1)] [-index2 (index2)] [-closest]\n");
    33     // gprint (GP_ERR, "  if -closest is provided, index1 & index2 will have the same length as X1 and X2 (respectively)\n");
    34     // gprint (GP_ERR, "    with either the index of the match or a value of -1 for non-matches\n");
    35     return (FALSE);
    36   }
     40    gprint (GP_ERR, "  if -closest is provided, index1 & index2 will have the same length as X1 and X2 (respectively)\n");
     41    gprint (GP_ERR, "    with either the index of the match or a value of -1 for non-matches\n");
     42    return (FALSE);
     43  }
     44
     45  /*
     46    we have two modes of operation: 
     47
     48    without -closest, we are finding all matched pairs within the match radius.  in this
     49    case, the two index vectors have the same length, one entry per matched pair.
     50    x1[index1],y1[index1] matches to x2[index2],y2[index2].
     51
     52    with -closest selected, we are finding the closest element of set 1 to each of set 2
     53    and vice versa.  in this case, index1 is always the same length as x1,y1, while index2
     54    is the same lengths as x2,y2.  x2[index1],y2[index1] matches x1,y1 while
     55    x1[index2],y1[index2] matches x2,y2
     56
     57   */
    3758
    3859  if ((X1vec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);   
    … …  
    6182  }
    6283
    63   find_matches2d (X1vec, Y1vec, X2vec, Y2vec, Radius, index1, index2);
     84  if (CLOSEST) {
     85      find_matches2d_closest (X1vec, Y1vec, X2vec, Y2vec, Radius, index1);
     86      find_matches2d_closest (X2vec, Y2vec, X1vec, Y1vec, Radius, index2);
     87  } else {
     88      find_matches2d (X1vec, Y1vec, X2vec, Y2vec, Radius, index1, index2);
     89  }
    6490
    6591  return (TRUE);
    … …  
    132158  return (TRUE);
    133159}
     160
     161// find the elements of X2,Y2 which are closest to each element of X1,Y1 (-1 if no match)
     162int find_matches2d_closest (Vector *X1, Vector *Y1, Vector *X2, Vector *Y2, double Radius, Vector *index) {
     163 
     164  off_t i, j, Jmin, Ji, I, J, *N1, *N2, NMATCH;
     165  double dX, dY, dR, Radius2, Rmin;
     166
     167  NMATCH = X1->Nelements;
     168  ResetVector (index, OPIHI_INT, NMATCH);
     169
     170  ALLOCATE (N1, off_t, X1->Nelements);
     171  ALLOCATE (N2, off_t, X2->Nelements);
     172
     173  for (i = 0; i < X1->Nelements; i++) { N1[i] = i; }
     174  for (i = 0; i < X2->Nelements; i++) { N2[i] = i; }
     175
     176  sort_coords_indexonly (X1->elements.Flt, Y1->elements.Flt, N1, X1->Nelements);
     177  sort_coords_indexonly (X2->elements.Flt, Y2->elements.Flt, N2, X2->Nelements);
     178
     179  Radius2 = Radius*Radius;
     180
     181  // find the closest entry in list 2 to the current entry in list 1:
     182  for (i = j = 0; (i < X1->Nelements) && (j < X2->Nelements);) {
     183    I = N1[i];
     184    J = N2[j];
     185
     186    dX = X1->elements.Flt[I] - X2->elements.Flt[J];
     187
     188    if (dX <= -1.02*Radius) {
     189      // no match in list 2 to this entry
     190      index->elements.Int[I] = -1;
     191      i++;
     192      continue;
     193    }
     194    if (dX >= +1.02*Radius) { j++; continue; }
     195
     196    // look for closest matches of list2() to list1(i)
     197    Jmin = -1;
     198    Rmin = Radius2;
     199    for (Ji = j; (dX > -1.02*Radius) && (Ji < X2->Nelements); Ji++) {
     200      J = N2[Ji];
     201      dX = X1->elements.Flt[I] - X2->elements.Flt[J];
     202      dY = Y1->elements.Flt[I] - Y2->elements.Flt[J];
     203      dR = dX*dX + dY*dY;
     204      if (dR > Radius2) continue;
     205      if (dR < Rmin) {
     206        Rmin = dR;
     207        Jmin  = J;
     208      }
     209    }
     210
     211    // no match in list 2 to this entry
     212    if (Jmin == -1) {
     213      index->elements.Int[I] = -1;
     214      i++;
     215      continue;
     216    }
     217    index->elements.Int[I] = Jmin;
     218    i++;
     219  }
     220  index->Nelements = NMATCH;
     221
     222  free (N1);
     223  free (N2);
     224
     225  return (TRUE);
     226}

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

    r29540 r30610  
    55int reindex (int argc, char **argv) {
    66 
    7   int  i, Npts, Nmax, valid;
     7  int  i, Npts, Nmax;
    88  Vector *ivec, *ovec, *xvec;
    99
    … …  
    1111  Npts = 0;
    1212
    13   valid = TRUE;
    14   valid &= (argc >= 6);
    15   valid &= !strcmp(argv[2], "=");
    16   valid &= !strcmp(argv[4], "using");
    17   if (!valid) {
    18     gprint (GP_ERR, "USAGE: reindex (out) = (in) using (index)\n");
    19     gprint (GP_ERR, "  creates a new vectors (out) from (in) based on sequence in (index)\n");
    20     return (FALSE);
    21   }
     13  if (argc != 6) goto usage;
     14  if (strcmp(argv[2], "=")) goto usage;
     15  if (strcmp(argv[4], "using")) goto usage;
    2216
    2317  if ((ovec = SelectVector (argv[1], ANYVECTOR, TRUE)) == NULL) goto error;
    … …  
    6357  DeleteVector (ovec);
    6458  return (FALSE);
     59
     60usage:
     61    gprint (GP_ERR, "USAGE: reindex (out) = (in) using (index)\n");
     62    gprint (GP_ERR, "  creates a new vectors (out) from (in) based on sequence in (index)\n");
     63    return (FALSE);
    6564}
    6665

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

    r27491 r30610  
    3434
    3535  // ovec matches ivec in type
    36   ResetVector (ovec, ivec->type, MAX (tvec[0].Nelements, 1));
     36  ResetVector (ovec, ivec->type, tvec[0].Nelements);
    3737
    3838  // we have four cases: (ivec == flt or int) and (tvec == flt or int)
    … …  
    7575
    7676  // free up unused memory
    77   ResetVector (ovec, ivec->type, MAX (Npts, 1));
     77  ResetVector (ovec, ivec->type, Npts);
    7878
    7979  DeleteVector (tvec);

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

    r21508 r30610  
    33int vstats (int argc, char **argv) {
    44 
    5   int i, N;
    6   double max, min, sum, var, dvar, mean, stdev;
    7   float IgnoreValue;
     5  int i, iter, N, Nbin, Niter;
     6  double max, min, pmin, pmax, sum, var, dvar, mean, stdev;
     7  float IgnoreValue, Nsigma;
    88  int Ignore, Quiet;
    99
    10   int *Nval, bin, Nmode, Nmed;
    11   double dx, mode, median;
     10  int *Nval, bin, Nmode, Nmed, Nused;
     11  double dx, mode, median, threshold;
    1212  Vector *vec;
     13  char *mask = NULL;
    1314
    1415  IgnoreValue = 0;
    … …  
    2122  }
    2223
     24  Niter = 1;
     25  Nsigma = 3.0;
     26
     27  if ((N = get_argument (argc, argv, "-sigma-clip"))) {
     28    remove_argument (N, &argc, argv);
     29    Nsigma = atof(argv[N]);
     30    Niter = 3;
     31    remove_argument (N, &argc, argv);
     32  }
     33
     34  if ((N = get_argument (argc, argv, "-iter"))) {
     35    remove_argument (N, &argc, argv);
     36    Niter = atoi(argv[N]);
     37    remove_argument (N, &argc, argv);
     38  }
     39
    2340  Quiet = FALSE;
    2441  if ((N = get_argument (argc, argv, "-q"))) {
    … …  
    3249
    3350  if (argc != 2) {
    34     gprint (GP_ERR, "USAGE: vstat (vector)\n");
     51    gprint (GP_ERR, "USAGE: vstat (vector) [-ignore value] [-q] [-quiet] [-iter Niter] [-sigma-clip Nsigma]\n");
     52    gprint (GP_ERR, "  default is 1 iteration without sigma clipping or 3 with sigma clipping\n");
    3553    return (FALSE);
    3654  }
    … …  
    3957
    4058  /* we need two passes, one for max, min, mean, sum, one for median, stdev, etc */
     59
     60  // set a good / bad mask
     61  ALLOCATE (mask, char, vec[0].Nelements);
     62  if (vec[0].type == OPIHI_FLT) {
     63    opihi_flt *X = vec[0].elements.Flt;
     64    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;
     69    }     
     70  } else {
     71    opihi_int *X = vec[0].elements.Int;
     72    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;
     77    }     
     78  }
    4179
    4280  /* calculate max, min, mean, sum, npix */
    … …  
    4785    opihi_flt *X = vec[0].elements.Flt;
    4886    for (i = 0; i < vec[0].Nelements; i++, X++) {
    49       if (!finite (*X)) continue;
    50       if (Ignore && (*X == IgnoreValue)) continue;
     87      if (mask[i]) continue;
    5188      max = MAX (*X, max);
    5289      min = MIN (*X, min);
    … …  
    5794    opihi_int *X = vec[0].elements.Int;
    5895    for (i = 0; i < vec[0].Nelements; i++, X++) {
    59       if (!finite (*X)) continue;
    60       if (Ignore && (*X == IgnoreValue)) continue;
     96      if (mask[i]) continue;
    6197      max = MAX (*X, max);
    6298      min = MIN (*X, min);
    … …  
    67103  mean = sum / N;
    68104
    69   /* calculate median and mode with resolution of (max - min) / 1000 */
    70   dx = (max - min) / 1000;
    71   if (dx == 0) {
    72     median = mode = min;
    73     stdev = 0.0;
    74     goto skip;
    75   }
    76 
    77   ALLOCATE (Nval, int, 1002);
    78   bzero (Nval, 1000*sizeof(int));
    79   var = 0;
    80   if (vec[0].type == OPIHI_FLT) {
    81     opihi_flt *X = vec[0].elements.Flt;
    82     for (i = 0; i < vec[0].Nelements; i++, X++) {
    83       if (!finite (*X)) continue;
    84       if (Ignore && (*X == IgnoreValue)) continue;
    85       bin = MAX (0, MIN (1000, (*X - min) / dx));
    86       Nval[bin] ++;
    87       dvar = (*X - mean);
    88       var += dvar*dvar;
    89     }     
    90   } else {
    91     opihi_int *X = vec[0].elements.Int;
    92     for (i = 0; i < vec[0].Nelements; i++, X++) {
    93       if (!finite (*X)) continue;
    94       if (Ignore && (*X == IgnoreValue)) continue;
    95       bin = MAX (0, MIN (1000, (*X - min) / dx));
    96       Nval[bin] ++;
    97       dvar = (*X - mean);
    98       var += dvar*dvar;
    99     }     
    100   }
    101   stdev = sqrt (var / N);
    102 
    103   Nmode = 0;
    104   mode = Nval[Nmode];
    105   median = 0;
    106   Nmed = -1;
    107   for (i = 0; i < 1001; i++) {
    108     if (Nmed == -1) {
    109       median += Nval[i];
    110       if (median >= N / 2.0) {
    111         Nmed = i;
    112         median = i * dx + min;
     105  // we do Niter passes; on each pass, we exclude entries > Nsigma from the median
     106  pmin = min;
     107  pmax = max;
     108
     109  stdev = mode = median = Nused = 0.0;
     110
     111  for (iter = 0; iter < Niter; iter ++) {
     112
     113    // reduce Nbin after the first iteration?
     114    Nbin = 1000;
     115
     116    /* calculate median and mode with resolution of (max - min) / 1000 */
     117    dx = (pmax - pmin) / Nbin;
     118    if (dx == 0) {
     119      median = mode = min;
     120      stdev = 0.0;
     121      goto skip;
     122    }
     123
     124    ALLOCATE (Nval, int, Nbin + 2);
     125    bzero (Nval, Nbin*sizeof(int));
     126    var = 0;
     127    if (vec[0].type == OPIHI_FLT) {
     128      opihi_flt *X = vec[0].elements.Flt;
     129      for (i = 0; i < vec[0].Nelements; i++, X++) {
     130        if (mask[i]) continue;
     131        bin = MAX (0, MIN (Nbin, (*X - pmin) / dx));
     132        Nval[bin] ++;
     133        dvar = (*X - mean);
     134        var += dvar*dvar;
     135      }     
     136    } else {
     137      opihi_int *X = vec[0].elements.Int;
     138      for (i = 0; i < vec[0].Nelements; i++, X++) {
     139        if (mask[i]) continue;
     140        bin = MAX (0, MIN (Nbin, (*X - pmin) / dx));
     141        Nval[bin] ++;
     142        dvar = (*X - mean);
     143        var += dvar*dvar;
     144      }     
     145    }
     146    stdev = sqrt (var / N);
     147
     148    Nmode = 0;
     149    mode = Nval[Nmode];
     150    median = 0;
     151    Nmed = -1;
     152    for (i = 0; i < Nbin + 1; i++) {
     153      if (Nmed == -1) {
     154        median += Nval[i];
     155        if (median >= N / 2.0) {
     156          Nmed = i;
     157          median = i * dx + pmin;
     158        }
    113159      }
    114     }
    115     if (mode < Nval[i]) {
    116       Nmode = i;
    117       mode = Nval[Nmode];
    118     }
    119   }
    120   mode = Nmode * dx + min;
    121   free (Nval);
     160      if (mode < Nval[i]) {
     161        Nmode = i;
     162        mode = Nval[Nmode];
     163      }
     164    }
     165    mode = Nmode * dx + pmin;
     166    free (Nval);
     167
     168    threshold = Nsigma * stdev;
     169
     170    // we are going to do another pass: mark the entries to skip
     171    pmin = max;
     172    pmax = min;
     173    if (iter < Niter - 1) {
     174      if (vec[0].type == OPIHI_FLT) {
     175        opihi_flt *X = vec[0].elements.Flt;
     176        for (i = 0; i < vec[0].Nelements; i++, X++) {
     177          if (mask[i]) continue;
     178          if (fabs(*X - median) > threshold) {
     179            mask[i] = 1;
     180          } else {
     181            pmin = MIN (*X, pmin);
     182            pmax = MAX (*X, pmax);
     183          }     
     184        }
     185      } else {
     186        opihi_int *X = vec[0].elements.Int;
     187        for (i = 0; i < vec[0].Nelements; i++, X++) {
     188          if (mask[i]) continue;
     189          if (fabs(*X - median) > threshold) {
     190            mask[i] = 1;
     191          } else {
     192            pmin = MIN (*X, pmin);
     193            pmax = MAX (*X, pmax);
     194          }
     195        }     
     196      }
     197    }
     198    // gprint (GP_ERR, "iter %d, mean: %g, stdev: %g, min: %g, max: %g, median: %g, mode: %g, Npts: %d\n",
     199    // iter, mean, stdev, min, max, median, mode, N);
     200  }
     201
     202  Nused = 0;
     203  for (i = 0; i < vec[0].Nelements; i++) {
     204    if (mask[i]) continue;
     205    Nused ++;
     206  }
    122207
    123208skip:
     209  FREE(mask);
     210
    124211  if (!Quiet) {
    125212    gprint (GP_ERR, "mean: %g, stdev: %g, min: %g, max: %g, median: %g, mode: %g, Npts: %d\n",
    126              mean, stdev, min, max, median, mode, N);
    127   }
    128 
     213            mean, stdev, min, max, median, mode, N);
     214  }
     215
     216  set_variable ("PMIN",      pmin);
     217  set_variable ("PMAX",      pmax);
    129218  set_variable ("MIN",      min);
    130219  set_variable ("MAX",      max);
    … …  
    134223  set_variable ("TOTAL",    sum);
    135224  set_int_variable ("NPIX", N);
     225  set_int_variable ("NPTS", N);
     226  set_int_variable ("NUSED", Nused);
    136227  set_variable ("SIGMA",    stdev);
    137228  return (TRUE);

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

    r20936 r30610  
    140140      if (fmttype[j] == 'd') {
    141141        if (vec[j][0].type == OPIHI_FLT) {
    142           fprintf (f, fmtlist[j], (int)(vec[j][0].elements.Flt[i]));
     142          fprintf (f, fmtlist[j], (opihi_int)(vec[j][0].elements.Flt[i]));
    143143        } else {
    144           fprintf (f, fmtlist[j], (int)(vec[j][0].elements.Int[i]));
     144          fprintf (f, fmtlist[j], (opihi_int)(vec[j][0].elements.Int[i]));
    145145        }
    146146      }
    147147      if (fmttype[j] == 'f') {
    148148        if (vec[j][0].type == OPIHI_FLT) {
    149           fprintf (f, fmtlist[j], (float)(vec[j][0].elements.Flt[i]));
     149          fprintf (f, fmtlist[j], (opihi_flt)(vec[j][0].elements.Flt[i]));
    150150        } else {
    151           fprintf (f, fmtlist[j], (float)(vec[j][0].elements.Int[i]));
     151          fprintf (f, fmtlist[j], (opihi_flt)(vec[j][0].elements.Int[i]));
    152152        }
    153153      }
Note: See TracChangeset for help on using the changeset viewer.