Changeset 30201

Timestamp:

Jan 5, 2011, 7:30:11 AM (15 years ago)

Author:

eugene

Message:

little improvements

Location:

branches/eam_branches/ipp-20101205/Ohana/src

Files:

• libkapa/src/DrawRotString.c (modified) (3 diffs)
• libkapa/src/PSRotFont.c (modified) (3 diffs)
• libkapa/src/bDrawRotFont.c (modified) (3 diffs)
• opihi/cmd.data/match2d.c (modified) (6 diffs)
• opihi/cmd.data/vstats.c (modified) (8 diffs)

branches/eam_branches/ipp-20101205/Ohana/src/libkapa/src/DrawRotString.c

@@ -27 +27 @@
   unsigned char *bitmap;
   char *currentname, basename[64]; 
-  int i, dy, dx, N, X, Y, code;
+  int i, dy, dx, N, X, Y, code, protect;
   int dX, Xoff, dY, Yoff, YoffBase;
   int currentsize, basesize;
@@ -63 +63 @@
 
   code = FALSE;
+  protect = FALSE;
 
   YoffBase = Yoff;
@@ -70 +71 @@
     if ((N < 0) || (N >= NROTCHARS)) continue;
 
+    if (N == 39) { // single-quote
+      protect = protect ? FALSE : TRUE;
+    } 
+
     /* check for special characters */
-    if (!code) {
+    if (!code && !protect) {
       if (N == 94) {
         SetRotFont (currentname, (int)(0.8*currentsize));

branches/eam_branches/ipp-20101205/Ohana/src/libkapa/src/PSRotFont.c

@@ -7 +7 @@
 
   char *segment, basename[64], *currentname;
-  int i, N, code;
+  int i, N, code, protect;
   int dX, dY, Xoff, Yoff, X, Y, Nseg, NSEG, YoffBase;
   double cs, sn, fscale, currentscale;
@@ -60 +60 @@
 
   code = FALSE;
+  protect = FALSE;
+
   YoffBase = 0;
   /* accumulate string segments with common state */
@@ -66 +68 @@
     if ((N < 0) || (N >= NROTCHARS)) continue;
 
+    if (N == 39) { // single-quote
+      protect = protect ? FALSE : TRUE;
+    } 
+
     /* check for special characters */
-    if (!code) {
+    if (!code && !protect) {
       /* superscript character (^) */
       if (N == 94) {

branches/eam_branches/ipp-20101205/Ohana/src/libkapa/src/bDrawRotFont.c

@@ -14 +14 @@
   unsigned char *bitmap;
   char *currentname, basename[64]; 
-  int i, dy, dx, N, X, Y, code;
+  int i, dy, dx, N, X, Y, code, protect;
   int dX, Xoff, dY, Yoff, YoffBase;
   int currentsize, basesize;
@@ -53 +53 @@
 
   code = FALSE;
+  protect = FALSE;
 
   YoffBase = Yoff;
@@ -60 +61 @@
     if ((N < 0) || (N >= NROTCHARS)) continue;
 
+    if (N == 39) { // single-quote
+      protect = protect ? FALSE : TRUE;
+    } 
+
     /* check for special characters */
-    if (!code) {
+    if (!code && !protect) {
       if (N == 94) {
         SetRotFont (currentname, (int)(0.8*currentsize));

branches/eam_branches/ipp-20101205/Ohana/src/opihi/cmd.data/match2d.c

@@ -1 +1 @@
 # include "data.h"
+
 int find_matches2d (Vector *X1, Vector *Y1, Vector *X2, Vector *Y2, double Radius, Vector *index1, Vector *index2);
+int find_matches2d_closest (Vector *X1, Vector *Y1, Vector *X2, Vector *Y2, double Radius, Vector *index);
 
 // match2d (X1) (Y1) (X2) (Y2) (Radius) [-index1 (index1)] [-index2 (index2)] [-nomatch1 nomatch1] [-nomatch2 nomatch2]
@@ -6 +8 @@
 int match2d (int argc, char **argv) {
   
-  int N;
+  int N, CLOSEST;
   double Radius;
   char *endptr;
@@ -12 +14 @@
   Vector *index1, *index2;
 
+  CLOSEST = FALSE;
+  if ((N = get_argument (argc, argv, "-closest"))) {
+    remove_argument (N, &argc, argv);
+    CLOSEST = TRUE;
+  }
+
   if ((N = get_argument (argc, argv, "-index1"))) {
     remove_argument (N, &argc, argv);
@@ -30 +38 @@
   if (argc != 6) {
     gprint (GP_ERR, "USAGE: match2d X1 Y1 X2 Y2 Radius [-index1 (index1)] [-index2 (index2)] [-closest]\n");
-    // gprint (GP_ERR, "USAGE: match2d X1 Y1 X2 Y2 Radius [-index1 (index1)] [-index2 (index2)] [-closest]\n");
-    // gprint (GP_ERR, "  if -closest is provided, index1 & index2 will have the same length as X1 and X2 (respectively)\n");
-    // gprint (GP_ERR, "    with either the index of the match or a value of -1 for non-matches\n");
-    return (FALSE);
-  }
+    gprint (GP_ERR, "  if -closest is provided, index1 & index2 will have the same length as X1 and X2 (respectively)\n");
+    gprint (GP_ERR, "    with either the index of the match or a value of -1 for non-matches\n");
+    return (FALSE);
+  }
+
+  /*
+    we have two modes of operation:  
+
+    without -closest, we are finding all matched pairs within the match radius.  in this
+    case, the two index vectors have the same length, one entry per matched pair.
+    x1[index1],y1[index1] matches to x2[index2],y2[index2].
+
+    with -closest selected, we are finding the closest element of set 1 to each of set 2
+    and vice versa.  in this case, index1 is always the same length as x1,y1, while index2
+    is the same lengths as x2,y2.  x2[index1],y2[index1] matches x1,y1 while
+    x1[index2],y1[index2] matches x2,y2
+
+   */
 
   if ((X1vec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);    
@@ -61 +82 @@
   }
 
-  find_matches2d (X1vec, Y1vec, X2vec, Y2vec, Radius, index1, index2);
+  if (CLOSEST) {
+      find_matches2d_closest (X1vec, Y1vec, X2vec, Y2vec, Radius, index1);
+      find_matches2d_closest (X2vec, Y2vec, X1vec, Y1vec, Radius, index2);
+  } else {
+      find_matches2d (X1vec, Y1vec, X2vec, Y2vec, Radius, index1, index2);
+  }
 
   return (TRUE);
@@ -132 +158 @@
   return (TRUE);
 }
+
+// find the elements of X2,Y2 which are closest to each element of X1,Y1 (-1 if no match)
+int find_matches2d_closest (Vector *X1, Vector *Y1, Vector *X2, Vector *Y2, double Radius, Vector *index) {
+  
+  off_t i, j, Jmin, Ji, I, J, *N1, *N2, NMATCH;
+  double dX, dY, dR, Radius2, Rmin;
+
+  NMATCH = X1->Nelements;
+  ResetVector (index, OPIHI_INT, NMATCH);
+
+  ALLOCATE (N1, off_t, X1->Nelements);
+  ALLOCATE (N2, off_t, X2->Nelements);
+
+  for (i = 0; i < X1->Nelements; i++) { N1[i] = i; }
+  for (i = 0; i < X2->Nelements; i++) { N2[i] = i; }
+
+  sort_coords_indexonly (X1->elements.Flt, Y1->elements.Flt, N1, X1->Nelements);
+  sort_coords_indexonly (X2->elements.Flt, Y2->elements.Flt, N2, X2->Nelements);
+
+  Radius2 = Radius*Radius;
+
+  // find the closest entry in list 2 to the current entry in list 1:
+  for (i = j = 0; (i < X1->Nelements) && (j < X2->Nelements);) {
+    I = N1[i];
+    J = N2[j];
+
+    dX = X1->elements.Flt[I] - X2->elements.Flt[J];
+
+    if (dX <= -1.02*Radius) { 
+      // no match in list 2 to this entry
+      index->elements.Int[I] = -1;
+      i++; 
+      continue; 
+    }
+    if (dX >= +1.02*Radius) { j++; continue; }
+
+    // look for closest matches of list2() to list1(i)
+    Jmin = -1;
+    Rmin = Radius2;
+    for (Ji = j; (dX > -1.02*Radius) && (Ji < X2->Nelements); Ji++) {
+      J = N2[Ji];
+      dX = X1->elements.Flt[I] - X2->elements.Flt[J];
+      dY = Y1->elements.Flt[I] - Y2->elements.Flt[J];
+      dR = dX*dX + dY*dY;
+      if (dR > Radius2) continue;
+      if (dR < Rmin) {
+        Rmin = dR;
+        Jmin  = J;
+      }
+    }
+
+    // no match in list 2 to this entry
+    if (Jmin == -1) {
+      index->elements.Int[I] = -1;
+      i++;
+      continue;
+    }
+    index->elements.Int[I] = Jmin;
+    i++;
+  }
+  index->Nelements = NMATCH;
+
+  free (N1);
+  free (N2);
+
+  return (TRUE);
+}

branches/eam_branches/ipp-20101205/Ohana/src/opihi/cmd.data/vstats.c

@@ -3 +3 @@
 int vstats (int argc, char **argv) {
   
-  int i, N;
-  double max, min, sum, var, dvar, mean, stdev;
-  float IgnoreValue;
+  int i, iter, N, Nbin, Niter;
+  double max, min, pmin, pmax, sum, var, dvar, mean, stdev;
+  float IgnoreValue, Nsigma;
   int Ignore, Quiet;
 
-  int *Nval, bin, Nmode, Nmed;
-  double dx, mode, median;
+  int *Nval, bin, Nmode, Nmed, Nused;
+  double dx, mode, median, threshold;
   Vector *vec;
+  char *mask = NULL;
 
   IgnoreValue = 0;
@@ -21 +22 @@
   }
 
+  Niter = 1;
+  Nsigma = 3.0;
+
+  if ((N = get_argument (argc, argv, "-sigma-clip"))) {
+    remove_argument (N, &argc, argv);
+    Nsigma = atof(argv[N]);
+    Niter = 3;
+    remove_argument (N, &argc, argv);
+  }
+
+  if ((N = get_argument (argc, argv, "-iter"))) {
+    remove_argument (N, &argc, argv);
+    Niter = atoi(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
   Quiet = FALSE;
   if ((N = get_argument (argc, argv, "-q"))) {
@@ -32 +49 @@
 
   if (argc != 2) {
-    gprint (GP_ERR, "USAGE: vstat (vector)\n");
+    gprint (GP_ERR, "USAGE: vstat (vector) [-ignore value] [-q] [-quiet] [-iter Niter] [-sigma-clip Nsigma]\n");
+    gprint (GP_ERR, "  default is 1 iteration without sigma clipping or 3 with sigma clipping\n");
     return (FALSE);
   }
@@ -39 +57 @@
 
   /* we need two passes, one for max, min, mean, sum, one for median, stdev, etc */
+
+  // set a good / bad mask
+  ALLOCATE (mask, char, vec[0].Nelements);
+  if (vec[0].type == OPIHI_FLT) {
+    opihi_flt *X = vec[0].elements.Flt;
+    for (i = 0; i < vec[0].Nelements; i++, X++) {
+      mask[i] = 1;
+      if (!finite (*X)) continue;
+      if (Ignore && (*X == IgnoreValue)) continue;
+      mask[i] = 0;
+    }      
+  } else {
+    opihi_int *X = vec[0].elements.Int;
+    for (i = 0; i < vec[0].Nelements; i++, X++) {
+      mask[i] = 1;
+      if (!finite (*X)) continue;
+      if (Ignore && (*X == IgnoreValue)) continue;
+      mask[i] = 0;
+    }      
+  }
 
   /* calculate max, min, mean, sum, npix */
@@ -47 +85 @@
     opihi_flt *X = vec[0].elements.Flt;
     for (i = 0; i < vec[0].Nelements; i++, X++) {
-      if (!finite (*X)) continue;
-      if (Ignore && (*X == IgnoreValue)) continue;
+      if (mask[i]) continue;
       max = MAX (*X, max);
       min = MIN (*X, min);
@@ -57 +94 @@
     opihi_int *X = vec[0].elements.Int;
     for (i = 0; i < vec[0].Nelements; i++, X++) {
-      if (!finite (*X)) continue;
-      if (Ignore && (*X == IgnoreValue)) continue;
+      if (mask[i]) continue;
       max = MAX (*X, max);
       min = MIN (*X, min);
@@ -67 +103 @@
   mean = sum / N;
 
-  /* calculate median and mode with resolution of (max - min) / 1000 */ 
-  dx = (max - min) / 1000;
-  if (dx == 0) {
-    median = mode = min;
-    stdev = 0.0;
-    goto skip;
-  }
-
-  ALLOCATE (Nval, int, 1002);
-  bzero (Nval, 1000*sizeof(int));
-  var = 0;
-  if (vec[0].type == OPIHI_FLT) {
-    opihi_flt *X = vec[0].elements.Flt;
-    for (i = 0; i < vec[0].Nelements; i++, X++) {
-      if (!finite (*X)) continue;
-      if (Ignore && (*X == IgnoreValue)) continue;
-      bin = MAX (0, MIN (1000, (*X - min) / dx));
-      Nval[bin] ++;
-      dvar = (*X - mean);
-      var += dvar*dvar;
-    }      
-  } else {
-    opihi_int *X = vec[0].elements.Int;
-    for (i = 0; i < vec[0].Nelements; i++, X++) {
-      if (!finite (*X)) continue;
-      if (Ignore && (*X == IgnoreValue)) continue;
-      bin = MAX (0, MIN (1000, (*X - min) / dx));
-      Nval[bin] ++;
-      dvar = (*X - mean);
-      var += dvar*dvar;
-    }      
-  }
-  stdev = sqrt (var / N);
-
-  Nmode = 0;
-  mode = Nval[Nmode];
-  median = 0;
-  Nmed = -1;
-  for (i = 0; i < 1001; i++) {
-    if (Nmed == -1) {
-      median += Nval[i];
-      if (median >= N / 2.0) {
-        Nmed = i;
-        median = i * dx + min;
+  // we do Niter passes; on each pass, we exclude entries > Nsigma from the median
+  pmin = min;
+  pmax = max;
+
+  for (iter = 0; iter < Niter; iter ++) {
+
+    // reduce Nbin after the first iteration?
+    Nbin = 1000;
+
+    /* calculate median and mode with resolution of (max - min) / 1000 */ 
+    dx = (pmax - pmin) / Nbin;
+    if (dx == 0) {
+      median = mode = min;
+      stdev = 0.0;
+      goto skip;
+    }
+
+    ALLOCATE (Nval, int, Nbin + 2);
+    bzero (Nval, Nbin*sizeof(int));
+    var = 0;
+    if (vec[0].type == OPIHI_FLT) {
+      opihi_flt *X = vec[0].elements.Flt;
+      for (i = 0; i < vec[0].Nelements; i++, X++) {
+        if (mask[i]) continue;
+        bin = MAX (0, MIN (Nbin, (*X - pmin) / dx));
+        Nval[bin] ++;
+        dvar = (*X - mean);
+        var += dvar*dvar;
+      }      
+    } else {
+      opihi_int *X = vec[0].elements.Int;
+      for (i = 0; i < vec[0].Nelements; i++, X++) {
+        if (mask[i]) continue;
+        bin = MAX (0, MIN (Nbin, (*X - pmin) / dx));
+        Nval[bin] ++;
+        dvar = (*X - mean);
+        var += dvar*dvar;
+      }      
+    }
+    stdev = sqrt (var / N);
+
+    Nmode = 0;
+    mode = Nval[Nmode];
+    median = 0;
+    Nmed = -1;
+    for (i = 0; i < Nbin + 1; i++) {
+      if (Nmed == -1) {
+        median += Nval[i];
+        if (median >= N / 2.0) {
+          Nmed = i;
+          median = i * dx + pmin;
+        }
       }
-    }
-    if (mode < Nval[i]) {
-      Nmode = i;
-      mode = Nval[Nmode];
-    }
-  }
-  mode = Nmode * dx + min;
-  free (Nval);
+      if (mode < Nval[i]) {
+        Nmode = i;
+        mode = Nval[Nmode];
+      }
+    }
+    mode = Nmode * dx + pmin;
+    free (Nval);
+
+    threshold = Nsigma * stdev;
+
+    // we are going to do another pass: mark the entries to skip
+    pmin = max;
+    pmax = min;
+    if (iter < Niter - 1) {
+      if (vec[0].type == OPIHI_FLT) {
+        opihi_flt *X = vec[0].elements.Flt;
+        for (i = 0; i < vec[0].Nelements; i++, X++) {
+          if (mask[i]) continue;
+          if (fabs(*X - median) > threshold) {
+            mask[i] = 1;
+          } else {
+            pmin = MIN (*X, pmin);
+            pmax = MAX (*X, pmax);
+          }      
+        }
+      } else {
+        opihi_int *X = vec[0].elements.Int;
+        for (i = 0; i < vec[0].Nelements; i++, X++) {
+          if (mask[i]) continue;
+          if (fabs(*X - median) > threshold) {
+            mask[i] = 1;
+          } else {
+            pmin = MIN (*X, pmin);
+            pmax = MAX (*X, pmax);
+          }
+        }      
+      }
+    }
+    // gprint (GP_ERR, "iter %d, mean: %g, stdev: %g, min: %g, max: %g, median: %g, mode: %g, Npts: %d\n", 
+    // iter, mean, stdev, min, max, median, mode, N);
+  }
+
+  Nused = 0;
+  for (i = 0; i < vec[0].Nelements; i++) {
+    if (mask[i]) continue;
+    Nused ++;
+  }
 
 skip:
+  FREE(mask);
+
   if (!Quiet) {
     gprint (GP_ERR, "mean: %g, stdev: %g, min: %g, max: %g, median: %g, mode: %g, Npts: %d\n", 
-             mean, stdev, min, max, median, mode, N);
-  }
-
+            mean, stdev, min, max, median, mode, N);
+  }
+
+  set_variable ("PMIN",      pmin);
+  set_variable ("PMAX",      pmax);
   set_variable ("MIN",      min);
   set_variable ("MAX",      max);
@@ -134 +221 @@
   set_variable ("TOTAL",    sum);
   set_int_variable ("NPIX", N);
+  set_int_variable ("NPTS", N);
+  set_int_variable ("NUSED", Nused);
   set_variable ("SIGMA",    stdev);
   return (TRUE);