Changeset 24243 for branches/cnb_branches/cnb_branch_20090301

Timestamp:

May 26, 2009, 1:47:40 PM (17 years ago)

Author:

beaumont

Message:

enhancements to relastro, relastrovisual

Location:

branches/cnb_branches/cnb_branch_20090301

Files:

• . (modified) (1 prop)
• Ohana (modified) (1 prop)
• Ohana/src/libdvo/include/dvo.h (modified) (1 diff)
• Ohana/src/opihi/dimm/camera_cmds.c (modified) (1 prop)
• Ohana/src/opihi/dimm/telescope_cmds.c (modified) (1 prop)
• Ohana/src/relastro/include/relastro.h (modified) (2 diffs)
• Ohana/src/relastro/src/FitChip.c (modified) (1 diff)
• Ohana/src/relastro/src/GetAstromError.c (modified) (1 diff)
• Ohana/src/relastro/src/ImageOps.c (modified) (8 diffs)
• Ohana/src/relastro/src/UpdateObjects.c (modified) (4 diffs)
• Ohana/src/relastro/src/args.c (modified) (1 diff)
• Ohana/src/relastro/src/bcatalog.c (modified) (2 diffs)
• Ohana/src/relastro/src/load_catalogs.c (modified) (1 diff)
• Ohana/src/relastro/src/relastroVisual.c (modified) (3 diffs)
• Ohana/src/relastro/src/relastroVisual.h (modified) (1 diff)
• Ohana/src/relastro/src/relastro_objects.c (modified) (4 diffs)
• ppImage (modified) (1 prop)
• ppStack (modified) (1 prop)
• ppSub (modified) (1 prop)
• psModules (modified) (1 prop)
• psModules/src/astrom/pmAstrometryModel.c (modified) (1 diff)
• psastro (modified) (1 prop)
• psastro/src/psastroModel.c (modified) (2 diffs)
• psphot (modified) (1 prop)
• pswarp/src/pswarpVersion.c (modified) (2 diffs)

branches/cnb_branches/cnb_branch_20090301/Ohana/src/libdvo/include/dvo.h

@@ -108 +108 @@
 # define ID_MEAS_POOR_ASTROM  0x0010  /* detection is astrometry outlier                                             */   
 # define ID_MEAS_SKIP_ASTROM  0x0020  /* detection was ignored for astrometry measurement                            */   
-# define ID_MEAS_UNDEF_1      0x0040  /* unused                                                                      */
-# define ID_MEAS_UNDEF_2      0x0080  /* unused                                                                      */
+# define ID_MEAS_USED_OBJ     0x0040  /* detection was used during opdate objects                                    */
+# define ID_MEAS_USED_CHIP    0x0080  /* detection was used during update chips                                      */
 # define ID_MEAS_BLEND_MEAS   0x0100  /* detection is within radius of multiple objects                              */
 # define ID_MEAS_BLEND_OBJ    0x0200  /* multiple detections within radius of object                                 */

branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/include/relastro.h

@@ -122 +122 @@
 
 int FlagOutlier;
+int    CLIP_THRESH;
 
 FitMode FIT_MODE;
@@ -268 +269 @@
 int UpdateMeasures (Catalog *catalog, int Ncatalog);
 void fixImageRaw (Catalog *catalog, int Ncatalog, int im);
-void FlagOutliers(Catalog *catalog, int Ncatalog);
+void FlagOutliers(Catalog *catalog);
 int MeasFilterTest(Measure *measure);
 

branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/FitChip.c

@@ -53 +53 @@
     // generate the fit matches
     for (i = 0; i < Nmatch; i++) {
-      if (raw[i].mask) continue;
-      if (isnan(raw[i].dMag) || raw[i].dMag > FIT_CHIP_MAX_ERROR) continue;
+      if (raw[i].mask) {
+        continue;
+      }
+      if (isnan(raw[i].dMag) || raw[i].dMag > FIT_CHIP_MAX_ERROR) {
+        continue;
+      }
 
       // only keep objects within dRmax

branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/GetAstromError.c

@@ -30 +30 @@
   
   dPtotal = sqrt(SQ(dPsys) + AS*SQ(dPobs) + MS*SQ(dM));
+
+  //XXX dXccd, dYccd are now working correctly
+  //dPtotal = AS * dPobs;
+
   dPtotal = MAX (dPtotal, MIN_ERROR);
 

branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/ImageOps.c

@@ -1 +1 @@
 # include "relastro.h"
+#define TESTING
 
 static unsigned int *start;
@@ -163 +164 @@
   plot_list (&graphdata, Mlist, dlist, Nimage, "Mcal vs dMcal", NULL);
 
-# define NBIN 200
+# define NBIN 20
   REALLOCATE (xlist, double, NBIN);
   REALLOCATE (Mlist, double, NBIN);
@@ -292 +293 @@
 
     // an object with only one detection provides no information about the image calibration
+    //XXX this is already taken care of in bcatalog
     raw[i].mask = FALSE;
     int mask = FALSE;
@@ -297 +299 @@
       mask = TRUE;
     }
-    if (~finite(catalog[c].measure[m].dR) || ~finite(catalog[c].measure[m].dD)) mask = TRUE;
+    if (!finite(catalog[c].measure[m].dR) || !finite(catalog[c].measure[m].dD)) {
+      mask = TRUE;
+    }
 
     raw[i].mask = mask;
@@ -388 +392 @@
 
 /** lifted from relphot/StarOps.clean_measures */
-void FlagOutliers (Catalog *catalog, int Ncatalog) {
+void FlagOutliers2D(Catalog *catalog);
+
+void FlagOutliers (Catalog *catalog) {
+  // XXX FlagOutliers is just using FlagOutliers2D
+  FlagOutliers2D(catalog);
+  return;
 
   int i, j, k, m, N, Ndel, Nave, Nmax, TOOFEW, Nsecfilt;
-  double Ns;
+  double Ns, theta, x, y;
   double *R, *D, *dR, *dD;
   StatType statsR, statsD;
@@ -399 +408 @@
 
   if (VERBOSE) fprintf (stderr, "marking poor measures\n");
-  /* Nmeasure is now different, need to reallocate */
   Nmax = 0;
-  for (i = 0; i < Ncatalog; i++) {
-    for (j = 0; j < catalog[i].Naverage; j++) {
-      Nmax = MAX (Nmax, catalog[i].average[j].Nmeasure);
-    }
-  }
+  for (i = 0; i < catalog[0].Naverage; i++) {
+    Nmax = MAX (Nmax, catalog[0].average[i].Nmeasure);
+  }
+
   ALLOCATE (R, double, Nmax);
   ALLOCATE (D, double, Nmax);
@@ -414 +421 @@
   TOOFEW = MAX (5, SRC_MEAS_TOOFEW);
 
-  Ns = 3;
+  Ns = CLIP_THRESH;
   Ndel = Nave = 0;
-  for (i = 0; i < Ncatalog; i++) {
-    for (j = 0; j < catalog[i].Naverage; j++) {
       
-      /* accumulate list of valid measurements */
-      m = catalog[i].average[j].measureOffset;
-      N = 0;
-      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
-
-        // skip measurements based on user selected criteria
-        if (!MeasFilterTest(&catalog[i].measure[m])) continue;
-
-        R[N] = getMeanR(&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
-        D[N] = getMeanD(&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
-        dR[N] = GetAstromError( &catalog[i].measure[m], ERROR_MODE_RA);
-        dD[N] = GetAstromError( &catalog[i].measure[m], ERROR_MODE_DEC);
-        if (isnan(R[N]) || isnan(D[N])) continue;
-        N++;
+  /* loop over each object in the catalog */
+  for (j = 0; j < catalog[0].Naverage; j++) {
+    
+    /* accumulate list of valid measurements */
+    m = catalog[0].average[j].measureOffset;
+    N = 0;
+    for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
+      
+      // skip measurements based on user selected criteria
+      if (!MeasFilterTest(&catalog[0].measure[m])) continue;
+      R[N] = catalog[0].measure[m].dR;
+      D[N] = catalog[0].measure[m].dD;
+      dR[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_RA);
+      dD[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_DEC);
+      if (isnan(R[N]) || isnan(D[N])) continue;
+      N++;
+    }
+    if (N <= TOOFEW) continue;
+    
+    /* 3-sigma clip based on stats of inner 50% */
+    initstats ("MEAN");
+    liststats (R, dR, N, &statsR);
+    liststats (D, dD, N, &statsD);
+    
+    statsR.sigma = MAX (MIN_ERROR, statsR.sigma);
+    statsD.sigma = MAX (MIN_ERROR, statsD.sigma);
+    
+    /* compare per-object distance to this standard deviation, and flag outliers*/
+    m = catalog[0].average[j].measureOffset;
+    N = 0;
+    for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
+      //reset flag on each invocation
+      catalog[0].measure[m].dbFlags &= ~ID_MEAS_POOR_ASTROM;
+
+      // skip measurements based on user selected criteria
+      if (!MeasFilterTest(&catalog[0].measure[m])) continue;
+      
+      x = catalog[0].measure[m].dR - statsR.median;
+      y = catalog[0].measure[m].dD - statsD.median;
+      theta = atan2(y,x);
+      if ((x*x + y*y) > (SQR(statsR.sigma * Ns * cos(theta)) + 
+                         SQR(statsD.sigma * Ns * sin(theta)))) {   
+        catalog[0].measure[m].dbFlags |= ID_MEAS_POOR_ASTROM;
+        Ndel++;
       }
-      if (N <= TOOFEW) continue;
-
-      /* 3-sigma clip based on stats of inner 50% */
-      initstats ("INNER_MEAN");
-      liststats (R, dR, N, &statsR);
-      liststats (R, dR, N, &statsD);
-
-      statsR.sigma = MAX (MIN_ERROR, statsR.sigma);
-      statsD.sigma = MAX (MIN_ERROR, statsD.sigma);
-
-      for (k = m = 0; k < N; k++) {
-        if ((fabs (R[k] - statsR.median) < Ns*statsR.sigma) && 
-            (fabs (D[k] - statsD.median) < Ns*statsD.sigma)) {
-          R[m] = R[k];
-          dR[m] = dR[k];
-          D[m] = D[k];
-          dD[m] = dD[k];
-          m++;
-        }
+      N++;
+      Nave ++;
+    }
+
+    //examine results
+    relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset, 
+                               catalog[0].average[j].Nmeasure, 
+                               statsR, statsD, Ns);
+  }
+  
+  if (VERBOSE) fprintf (stderr, "%d measures marked poor, %d total\n", Ndel, Nave);
+  free (R);
+  free(dR);
+  free(D);
+  free(dD); 
+}
+
+
+/** an alternative outlier rejection scheme */
+void FlagOutliers2D (Catalog *catalog) {
+
+  int i, j, k, m, N, Ndel, Nave, Nmax, TOOFEW, Nsecfilt;
+  double *index;
+  double Ns, theta, x, y;
+  double *R, *D, *dR, *dD, *d2;
+  StatType statsR, statsD;
+
+  Nsecfilt = GetPhotcodeNsecfilt();
+  assert(catalog[0].Nsecfilt == Nsecfilt);
+
+  if (VERBOSE) fprintf (stderr, "marking poor measures\n");
+  Nmax = 0;
+  for (i = 0; i < catalog[0].Naverage; i++) {
+    Nmax = MAX (Nmax, catalog[0].average[i].Nmeasure);
+  }
+
+  ALLOCATE (R, double, Nmax);
+  ALLOCATE (D, double, Nmax);
+  ALLOCATE (dR, double, Nmax);
+  ALLOCATE (dD, double, Nmax);
+  ALLOCATE (d2, double, Nmax);
+  ALLOCATE (index, double, Nmax);
+
+  /* it makes no sense to mark 3-sigma outliers with <5 measurements */
+  TOOFEW = MAX (5, SRC_MEAS_TOOFEW);
+
+  Ns = CLIP_THRESH;
+  Ndel = Nave = 0;
+      
+  /* loop over each object in the catalog */
+  for (j = 0; j < catalog[0].Naverage; j++) {
+    
+    /* accumulate list of valid measurements */
+    m = catalog[0].average[j].measureOffset;
+    N = 0;
+    
+    for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
+      
+      //reset flag on each invocation
+      catalog[0].measure[m].dbFlags &= ~ID_MEAS_POOR_ASTROM;
+      
+      // skip measurements based on user selected criteria
+      if (!MeasFilterTest(&catalog[0].measure[m])) continue;
+      R[N] = catalog[0].measure[m].dR;
+      D[N] = catalog[0].measure[m].dD;
+      dR[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_RA);
+      dD[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_DEC);
+      if (isnan(R[N]) || isnan(D[N])) continue;
+      N++;
+    }
+    if (N <= TOOFEW) continue;
+    
+    /* calculate mean of all points*/
+    initstats ("MEAN");
+    liststats (R, dR, N, &statsR);
+    liststats (D, dD, N, &statsD);
+    statsR.sigma = MAX (MIN_ERROR, statsR.sigma);
+    statsD.sigma = MAX (MIN_ERROR, statsD.sigma);
+    
+    /* calculate deviations of all points*/
+    m = catalog[0].average[j].measureOffset;
+    N = 0;
+    for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
+      //skip bad measurements
+      if (!MeasFilterTest(&catalog[0].measure[m])) continue;  
+      x = catalog[0].measure[m].dR - statsR.median;
+      y = catalog[0].measure[m].dD - statsD.median;
+      theta = atan2(y,x);
+      d2[N] = (x*x + y*y) / (SQR(statsR.sigma * Ns * cos(theta)) + 
+                             SQR(statsD.sigma * Ns * sin(theta)));      
+      index[N] = m;
+      N++;
+    }
+    
+    // sort d2
+    dsortpair(d2, index, N);
+    N = (N/2 > (N-1)) ? N/2 : N-1;
+
+    // recalculate image center, sigma based on closest 50% of points
+    for(k = 0;  k < N; k++) {
+      int ind = (int) index[k];
+      R[k] = catalog[0].measure[ind].dR;
+      D[k] = catalog[0].measure[ind].dD;
+      dR[k] = GetAstromError( &catalog[0].measure[ind], ERROR_MODE_RA);
+      dD[k] = GetAstromError( &catalog[0].measure[ind], ERROR_MODE_DEC);
+    }
+    liststats (R, dR, N, &statsR);
+    liststats (D, dD, N, &statsD);
+    statsR.sigma = MAX (MIN_ERROR, statsR.sigma);
+    statsD.sigma = MAX (MIN_ERROR, statsD.sigma);
+    
+    // use these new statistics to flag outliers 
+    m = catalog[0].average[j].measureOffset;
+    N = 0;
+    for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
+      //skip bad measurements
+      if (!MeasFilterTest(&catalog[0].measure[m])) continue;  
+      x = catalog[0].measure[m].dR - statsR.median;
+      y = catalog[0].measure[m].dD - statsD.median;
+      theta = atan2(y,x);
+      d2[N] = (x*x + y*y) / (SQR(statsR.sigma * Ns * cos(theta)) + 
+                             SQR(statsD.sigma * Ns * sin(theta)));      
+      if ((d2[N]) > 1) {
+        catalog[0].measure[m].dbFlags |= ID_MEAS_POOR_ASTROM;
+        Ndel ++;
       }
-      initstats ("MEAN");
-      liststats (R, dR, m, &statsR);
-      liststats (D, dD, m, &statsD);
-      statsR.sigma = MAX (MIN_ERROR, statsR.sigma);
-      statsD.sigma = MAX (MIN_ERROR, statsD.sigma);
-
-      /* apply to list of all relevant measurements*/
-      m = catalog[i].average[j].measureOffset;
-      N = 0;
-      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
-        // skip measurements based on user selected criteria
-        if (!MeasFilterTest(&catalog[i].measure[m])) continue;
-
-        R[N] = getMeanR(&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
-        D[N] = getMeanD(&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
-        dR[N] = GetAstromError( &catalog[i].measure[m], ERROR_MODE_RA);
-        dD[N] = GetAstromError( &catalog[i].measure[m], ERROR_MODE_DEC);
-        
-        if((fabs(R[N] - statsR.median) > Ns * statsR.sigma) ||
-           (fabs(D[N] - statsR.median) > Ns * statsD.sigma)) {
-          catalog[i].measure[m].dbFlags |= ID_MEAS_POOR_ASTROM;
-          Ndel++;
-        }
-        N++;
-        Nave ++;
-      }
-      initstats (STATMODE);
-   
-      if (TRUE) fprintf (stderr, "%d measures marked poor, %d total\n", Ndel, Nave);
-      free (R);
-      free(dR);
-      free(D);
-      free(dD);
-    }
-  }
-}
-
- 
+      N++;
+      Nave++;
+    }  //done rejecting outliers
+
+    //examine results
+    relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset, 
+                               catalog[0].average[j].Nmeasure, 
+                               statsR, statsD, Ns);
+    
+  } //done looping over objects
+  
+  if (VERBOSE) fprintf (stderr, "%d measures marked poor, %d total\n", Ndel, Nave);
+  free (R);
+  free(dR);
+  free(D);
+  free(dD); 
+  free(d2);
+  free(index);
+}
+
+
+
 /** Determine whether a measurement should be included in the analysis, based on supplied filter criteria */ 
 int MeasFilterTest(Measure *measure) {
@@ -501 +626 @@
   if (!finite(measure[0].dR) || !finite(measure[0].dD)) return FALSE;
   if (!finite(measure[0].M)) return FALSE; //XXX is this necessary for all relastro tasks?
-
+  
   /* select measurements by photcode, or equiv photcode, if specified */
   if (NphotcodesKeep > 0) {

branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/UpdateObjects.c

@@ -98 +98 @@
 
         //does the measurement pass the supplied filtering constraints?
-        if (!MeasFilterTest(&catalog[i].measure[m])) continue;
+        if (!MeasFilterTest(&catalog[i].measure[m])) {
+          catalog[i].measure[m].dbFlags &= ~ID_MEAS_USED_OBJ;
+          continue;
+        }
 
         //outlier rejection
-        if (FlagOutliers && (catalog[i].measure[m].dbFlags & ID_MEAS_POOR_ASTROM)) {
+        if (FlagOutlier && (catalog[i].measure[m].dbFlags & ID_MEAS_POOR_ASTROM)) {
+          catalog[i].measure[m].dbFlags &= ~ID_MEAS_USED_OBJ;
           continue;
         }
@@ -114 +118 @@
         # endif
 
+        catalog[i].measure[m].dbFlags |= ID_MEAS_USED_OBJ;
+
         R[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
         D[N] = getMeanD (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
@@ -141 +147 @@
       if (N < SRC_MEAS_TOOFEW) {
         // XXX need to define PHOTOM and ASTROM object flags
-        // catalog[i].average[j].code |= ID_STAR_FEW;
-        continue;
+        catalog[i].average[j].flags |= ID_STAR_FEW;
+        if (N < 2) continue;
       }
 
@@ -236 +242 @@
       }
 
+      //make sure that the fit succeeded
+      assert(finite(fit.Ro) && finite(fit.Do) && 
+             finite(fit.dRo) && finite(fit.dDo) &&
+             finite(fit.uR) && finite(fit.uD) &&
+             finite(fit.duR) && finite(fit.duD) &&
+             finite(fit.p) && finite(fit.dp));
+
       // the measure fields must be updated before the average fields
       m = catalog[i].average[j].measureOffset;

branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/args.c

@@ -40 +40 @@
     FIT_TARGET = TARGET_MOSAICS;
   }
+
+  FlagOutlier = FALSE;
   if ((N = get_argument (argc, argv, "-clip"))) {
     remove_argument (N, &argc, argv);
+    CLIP_THRESH = atof (argv[N]);
+    remove_argument (N, &argc, argv);
     FlagOutlier = TRUE;
   }
+
   if (FIT_TARGET == TARGET_NONE) usage();
 

branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/bcatalog.c

@@ -48 +48 @@
       
       //filter objects based on user supplied criteria
-      if (!MeasFilterTest(&catalog[0].measure[offset])) continue;
+      if (!MeasFilterTest(&catalog[0].measure[offset])) {
+        catalog[0].measure[offset].dbFlags &= ~ID_MEAS_USED_CHIP;
+        continue;
+      }
 
       //filter out outliers
-      if (FlagOutlier && (catalog[0].measure[offset].dbFlags & ID_MEAS_POOR_ASTROM)) continue;
-      
+      if (FlagOutlier && (catalog[0].measure[offset].dbFlags & ID_MEAS_POOR_ASTROM)) {
+        catalog[0].measure[offset].dbFlags &= ~ID_MEAS_USED_CHIP;
+        continue;
+      }
+      catalog[0].measure[offset].dbFlags |= ID_MEAS_USED_CHIP;
+
       // re-assess on each run of relastro if a measurement should be used
 
@@ -73 +80 @@
         subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_AREA;
       }
+     
       Nmeasure ++;
       Nm ++;

branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/load_catalogs.c

@@ -32 +32 @@
     //outlier rejection
     if (FlagOutlier) {
-      FlagOutliers(&tcatalog, 1);
+      FlagOutliers(&tcatalog);
     }
 

branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/relastroVisual.c

@@ -16 +16 @@
 
 static int isVisual = FALSE;
-static int plotRawRef = TRUE;
-static int plotScatter = TRUE;
-static int plotResid = TRUE;
-static int plotVector = TRUE;
+static int plotRawRef = FALSE;
+static int plotScatter = FALSE;
+static int plotResid = FALSE;
+static int plotVector = FALSE;
+static int plotOutliers = TRUE;
 
 /** Spawn a kapa window */
@@ -75 +76 @@
                      float xVec[], float yVec[],
                      int npts, int *kapaID) {
-    if (!plotResid) return 1;
+    if (!isVisual || !plotResid) return TRUE;
     Graphdata graphdata;
     KapaSection section;
 
-    if (!initWindow(kapaID)) return 0;
+    if (!initWindow(kapaID)) return FALSE;
 
     KapaInitGraph(&graphdata);
@@ -389 +390 @@
   return 1;
 }
+
+int relastroVisualPlotOutliers(Catalog *catalog, int offset, int Nmeasure,
+                               StatType statsR, StatType statsD, double thresh) {
+    
+  float *Din, *Rin, *Dout, *Rout;
+  double xmin, xmax, ymin, ymax, range;
+  float xCirc[100], yCirc[100];
+  int m, i;
+  int Nin, Nout;
+  Measure meas;
+  Graphdata graphdata;
+  KapaSection section;
+  
+  if (!isVisual || !plotOutliers) return 1;
+  if (!initWindow(&kapa)) return 0;
+  
+  // populate vectors
+  ALLOCATE(Din, float, Nmeasure);
+  ALLOCATE(Rin, float, Nmeasure);
+  ALLOCATE(Dout, float, Nmeasure);
+  ALLOCATE(Rout, float, Nmeasure);
+  
+  //create the threshhold ellipse
+  for(i = 0; i < 100; i++) {
+    xCirc[i] = statsR.median + thresh * statsR.sigma * cos(2 * 3.14 / 99. * i);
+    yCirc[i] = statsD.median + thresh * statsD.sigma * sin(2 * 3.14 / 99. * i);  
+  }
+
+  m = offset;
+  Nin = Nout = 0;
+  xmin = +FLT_MAX;
+  xmax = -FLT_MAX;
+  ymax = -FLT_MAX;
+  ymin = +FLT_MAX;
+  for(i = 0; i < Nmeasure; i++, m++) {
+    meas = catalog[0].measure[m];
+    if (!MeasFilterTest(&meas)) continue;
+    xmin = MIN(xmin, meas.dR);
+    xmax = MAX(xmax, meas.dR);
+    ymin = MIN(ymin, meas.dD);
+    ymax = MAX(ymax, meas.dD);
+    
+    if (meas.dbFlags & ID_MEAS_POOR_ASTROM) {
+      Rout[Nout] = (float)(meas.dR);
+      Dout[Nout] = (float)(meas.dD);
+      fprintf(stderr, "r: %f\td: %f\t outlier: 1\n", Rout[Nout], Dout[Nout]);
+      Nout++;
+    } else {
+      Rin[Nin] = (float)(meas.dR);
+      Din[Nin] = (float)(meas.dD);
+      fprintf(stderr, "r: %f\td: %f\t outlier: 0\n", Rin[Nin], Din[Nin]);
+      Nin++;
+    }
+  }
+ 
+  range = (xmax - xmin);
+  xmin -= .1 * range;
+  xmax += .1 * range;
+  range = (ymax - ymin);
+  ymax += .1 * range;
+  ymin -= .1 * range;
+
+  //temporary fix
+  xmin = -1; xmax = 1; ymin = -1; ymax = 1;
+
+
+  //initialize graph info
+  section.x = 0; section.y = 0; section.dx = 1; section.dy = 1;
+  section.name = "junk";
+  
+  KapaInitGraph(&graphdata);
+  KapaClearPlots(kapa);
+  KapaSetFont(kapa, "helvetica", 14);
+ 
+  graphdata.ptype = 7;
+  graphdata.style = 2;
+  graphdata.size = 3;
+  graphdata.xmin = xmin;
+  graphdata.xmax = xmax;
+  graphdata.ymin = ymin;
+  graphdata.ymax = ymax;
+
+  KapaSetSection(kapa, &section);
+  KapaSetLimits(kapa, &graphdata);
+  KapaBox(kapa, &graphdata);
+
+  KapaSendLabel( kapa, "RA (arcsec)", KAPA_LABEL_XM);
+  KapaSendLabel( kapa, "Dec (arcsec)", KAPA_LABEL_YM);
+  KapaSendLabel( kapa, "Points flagged as outliers (red)",
+                 KAPA_LABEL_XP);
+
+  graphdata.color = KapaColorByName("green");
+  KapaPrepPlot(kapa, Nin, &graphdata);
+  KapaPlotVector(kapa, Nin, Rin, "x");
+  KapaPlotVector(kapa, Nin, Din, "y");
+
+  graphdata.color = KapaColorByName("red");
+  KapaPrepPlot(kapa, Nout, &graphdata);
+  KapaPlotVector(kapa, Nout, Rout, "x");
+  KapaPlotVector(kapa, Nout, Dout, "y");
+
+  graphdata.color = KapaColorByName("black");
+  graphdata.ptype = 0;
+  graphdata.style = 0;
+  KapaPrepPlot(kapa, 100, &graphdata);
+  KapaPlotVector(kapa, 100, xCirc, "x");
+  KapaPlotVector(kapa, 100, yCirc, "y");
+
+
+  
+  askUser(&plotOutliers);
+
+  FREE(Rout);
+  FREE(Dout);
+  FREE(Rin);
+  FREE(Din);
+
+  return TRUE;
+}
+  
+  

branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/relastroVisual.h

@@ -4 +4 @@
 int relastroVisualPlotRawRef(StarData *raw, StarData *ref, double dRmax, int numObj);
 int relastroVisualPlotScatter(double values[], double thresh, int npts);
+int relastroVisualPlotOutliers(Catalog *catalog, int offset, int Nmeasure, 
+                               StatType statsR, StatType statsD, double thresh);
 
 #endif

branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/relastro_objects.c

@@ -9 +9 @@
   Catalog catalog;
 
-  //CNB
-  fprintf(stderr, "beginning relastro objects\n");
 
   // load the current sky table (layout of all SkyRegions) 
@@ -19 +17 @@
   skylist = SkyListByPatch (sky, -1, &UserPatch);
 
-  //CNB
-  fprintf(stder, "looping over skylist regions\n");
   // load data from each region file, only use bright stars
   for (i = 0; i < skylist[0].Nregions; i++) {
@@ -42 +38 @@
     }
 
-    if (FlagOutliers) {
-      FlagOutliers(&catalog, 1);
+    if (FlagOutlier) {
+      FlagOutliers(&catalog);
     }
 
@@ -57 +53 @@
     }
 
-    //CNB
-    fprintf(stderr, "updating objects\n");
     UpdateObjects (&catalog, 1);
 

branches/cnb_branches/cnb_branch_20090301/psModules/src/astrom/pmAstrometryModel.c

@@ -471 +471 @@
         char *thisone = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME");
         if (!thisone) continue;
-        if (!strcmp (name, thisone)) return (chip);
+        if (!strcmp (name, thisone)) return (chip);
     }
     return NULL;

branches/cnb_branches/cnb_branch_20090301/psastro/src/psastroModel.c

@@ -46 +46 @@
         exit (1);
     }
-    
+
     // run the full astrometry analysis (chip and/or mosaic)
     if (!psastroModelAdjust (config)) {
@@ -52 +52 @@
         exit (1);
     }
-    
+
     // save the model
     if (!psastroModelDataSave (config)) {

branches/cnb_branches/cnb_branch_20090301/pswarp/src/pswarpVersion.c

@@ -10 +10 @@
  *  Copyright 2009 Institute for Astronomy, University of Hawaii
  */
-
 #ifdef HAVE_CONFIG_H
 #include <config.h>
@@ -20 +19 @@
 #include <psphot.h>
 #include <ppStats.h>
-#include "pswarp.h"
-
 
 psString pswarpVersion(void)