Index: /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/include/relastro.h
===================================================================
--- /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/include/relastro.h	(revision 23648)
+++ /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/include/relastro.h	(revision 23649)
@@ -120,4 +120,6 @@
 int TimeSelect;
 time_t TSTART, TSTOP;
+
+int FlagOutlier;
 
 FitMode FIT_MODE;
@@ -266,4 +268,6 @@
 int UpdateMeasures (Catalog *catalog, int Ncatalog);
 void fixImageRaw (Catalog *catalog, int Ncatalog, int im);
+void FlagOutliers(Catalog *catalog, int Ncatalog);
+int MeasFilterTest(Measure *measure);
 
 int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon);
Index: /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/ImageOps.c
===================================================================
--- /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/ImageOps.c	(revision 23648)
+++ /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/ImageOps.c	(revision 23649)
@@ -386,2 +386,165 @@
   return (ref);
 }
+
+/** lifted from relphot/StarOps.clean_measures */
+void FlagOutliers (Catalog *catalog, int Ncatalog) {
+
+  int i, j, k, m, N, Ndel, Nave, Nmax, TOOFEW, Nsecfilt;
+  double Ns;
+  double *R, *D, *dR, *dD;
+  StatType statsR, statsD;
+
+  Nsecfilt = GetPhotcodeNsecfilt();
+  assert(catalog[0].Nsecfilt == Nsecfilt);
+
+  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);
+    }
+  }
+  ALLOCATE (R, double, Nmax);
+  ALLOCATE (D, double, Nmax);
+  ALLOCATE (dR, double, Nmax);
+  ALLOCATE (dD, double, Nmax);
+
+  /* it makes no sense to mark 3-sigma outliers with <5 measurements */
+  TOOFEW = MAX (5, SRC_MEAS_TOOFEW);
+
+  Ns = 3;
+  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++;
+      }
+      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++;
+	}
+      }
+      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);
+    }
+  }
+}
+
+ 
+/** Determine whether a measurement should be included in the analysis, based on supplied filter criteria */ 
+int MeasFilterTest(Measure *measure) {
+  int found, k;
+  long mask;
+  PhotCode *code;
+  float mag;
+
+  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) {
+    found = FALSE;
+    for (k = 0; (k < NphotcodesKeep) && !found; k++) {
+      if (photcodesKeep[k][0].code == measure[0].photcode) found = TRUE;
+      if (photcodesKeep[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE;
+    }
+    if (!found) return FALSE;
+  }
+  
+  if (NphotcodesSkip > 0) {
+    found = FALSE;
+    for (k = 0; (k < NphotcodesSkip) && !found; k++) {
+      if (photcodesSkip[k][0].code == measure[0].photcode) found = TRUE;
+      if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE;
+    }
+    if (found) return FALSE;
+  }  
+  
+  /* select measurements by time */
+  if (TimeSelect) {
+    if (measure[0].t < TSTART) return FALSE;
+    if (measure[0].t > TSTOP) return FALSE;
+  }
+  
+  /* select measurements by quality */
+  if (PhotFlagSelect) {
+    if (PhotFlagBad) {
+      mask = PhotFlagBad;
+    } else {
+      code = GetPhotcodebyCode (measure[0].photcode);
+      mask = code[0].astromBadMask;
+    }
+    if (mask & measure[0].photFlags) return FALSE;
+  }
+
+  /* select measurements by measurement error */
+  if ((SIGMA_LIM > 0) && (measure[0].dM > SIGMA_LIM)) return FALSE;
+  
+  /* select measurements by mag limit */
+  if (ImagSelect) {
+    mag = PhotInst (measure);
+    if (mag < ImagMin || mag > ImagMax) return FALSE;
+  }
+  
+  return TRUE;
+}
Index: /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/UpdateObjects.c
===================================================================
--- /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/UpdateObjects.c	(revision 23648)
+++ /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/UpdateObjects.c	(revision 23649)
@@ -97,6 +97,11 @@
       for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
 
-	//exclude measurements which have non-finite astrometry
-	if (~finite(catalog[i].measure[m].dR) || ~finite(catalog[i].measure[m].dD)) continue;
+	//does the measurement pass the supplied filtering constraints?
+	if (!MeasFilterTest(&catalog[i].measure[m])) continue;
+
+	//outlier rejection
+	if (FlagOutliers && (catalog[i].measure[m].dbFlags & ID_MEAS_POOR_ASTROM)) {
+	  continue;
+	}
 
 	// exclude measurements by previous outlier detection
@@ -109,49 +114,4 @@
 	# endif
 
-	/* exclude measurements by quality */
-	if (PhotFlagSelect) {
-	  if (PhotFlagBad) {
-	    mask = PhotFlagBad;
-	  } else {
-	    code = GetPhotcodebyCode (catalog[i].measure[m].photcode);
-	    mask = code[0].astromBadMask;
-	  }
-	  if (mask & catalog[i].measure[m].photFlags) {
-	    catalog[i].measure[m].dbFlags |= ID_MEAS_SKIP_ASTROM;
-	    continue;
-	  }
-	}
-	
-	/* exclude measurements by mag limit */
-	if (ImagSelect) {
-	  mag = PhotInst (&catalog[i].measure[m]);
-	  if (mag < ImagMin) {
-	    catalog[i].measure[m].dbFlags |= ID_MEAS_SKIP_ASTROM;
-	    continue;
-	  }
-	  if (mag > ImagMax) {
-	    catalog[i].measure[m].dbFlags |= ID_MEAS_SKIP_ASTROM;
-	    continue;
-	  }
-	}
-
-	/* select or exclude measurements by photcode, or equiv photcode, if specified */
-	if (NphotcodesKeep > 0) {
-	  found = FALSE;
-	  for (kp = 0; (kp < NphotcodesKeep) && !found; kp++) {
-	    if (photcodesKeep[kp][0].code == catalog[i].measure[m].photcode) found = TRUE;
-	    if (photcodesKeep[kp][0].code == GetPhotcodeEquivCodebyCode(catalog[i].measure[m].photcode)) found = TRUE;
-	  }
-	  if (!found) continue;
-	}
-	if (NphotcodesSkip > 0) {
-	  found = FALSE;
-	  for (kp = 0; (kp < NphotcodesSkip) && !found; kp++) {
-	    if (photcodesSkip[kp][0].code == catalog[i].measure[m].photcode) found = TRUE;
-	    if (photcodesSkip[kp][0].code == GetPhotcodeEquivCodebyCode(catalog[i].measure[m].photcode)) found = TRUE;
-	  }
-	  if (found) 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]);
@@ -179,5 +139,5 @@
 
       // too few measurements for average position (require 2 values)
-      if (N < 2) {
+      if (N < SRC_MEAS_TOOFEW) {
 	// XXX need to define PHOTOM and ASTROM object flags
 	// catalog[i].average[j].code |= ID_STAR_FEW;
Index: /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/args.c
===================================================================
--- /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/args.c	(revision 23648)
+++ /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/args.c	(revision 23649)
@@ -39,4 +39,8 @@
     remove_argument (N, &argc, argv);
     FIT_TARGET = TARGET_MOSAICS;
+  }
+  if ((N = get_argument (argc, argv, "-clip"))) {
+    remove_argument (N, &argc, argv);
+    FlagOutlier = TRUE;
   }
   if (FIT_TARGET == TARGET_NONE) usage();
Index: /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/bcatalog.c
===================================================================
--- /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/bcatalog.c	(revision 23648)
+++ /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/bcatalog.c	(revision 23649)
@@ -8,5 +8,4 @@
   int mask;
   PhotCode *code;
-  int skipFew = skipPhotCodeKeep = skipPhotCodeSkip = skipTime = skipPhotFlag = skipSigmaLim = skipImagSelect = 0;
 
   // XXX in the future, use catalog[0].Nsecfilt only?  allow catalogs to have variable Nsecfilt?
@@ -25,5 +24,4 @@
   for (i = 0; i < catalog[0].Naverage; i++) {
     if (catalog[0].average[i].Nmeasure <= SRC_MEAS_TOOFEW) {
-      skipFew++;
       continue;
     }
@@ -48,73 +46,11 @@
 
       offset = catalog[0].average[i].measureOffset + j;
+      
+      //filter objects based on user supplied criteria
+      if (!MeasFilterTest(&catalog[0].measure[offset])) continue;
 
-      /* select measurements by photcode, or equiv photcode, if specified */
-      if (NphotcodesKeep > 0) {
-        found = FALSE;
-        for (k = 0; (k < NphotcodesKeep) && !found; k++) {
-          if (photcodesKeep[k][0].code == catalog[0].measure[offset].photcode) found = TRUE;
-          if (photcodesKeep[k][0].code == GetPhotcodeEquivCodebyCode(catalog[0].measure[offset].photcode)) found = TRUE;
-        }
-        if (!found) {
-	  skipPhotCodeKeep++;
-	  continue;
-	}
-      }
-      if (NphotcodesSkip > 0) {
-        found = FALSE;
-        for (k = 0; (k < NphotcodesSkip) && !found; k++) {
-          if (photcodesSkip[k][0].code == catalog[0].measure[offset].photcode) found = TRUE;
-          if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(catalog[0].measure[offset].photcode)) found = TRUE;
-        }
-        if (found) {
-	  skipPhotCodeSkip++;
-	  continue;
-	}
-      }
-
-      /* select measurements by time */
-      if (TimeSelect) {
-        if (catalog[0].measure[offset].t < TSTART) {
-	  skipTime++;
-	  continue;
-	}
-        if (catalog[0].measure[offset].t > TSTOP) {
-	  skipTime++;
-	  continue;
-	}
-      }
+      //filter out outliers
+      if (FlagOutlier && (catalog[0].measure[offset].dbFlags & ID_MEAS_POOR_ASTROM)) continue;
       
-      /* select measurements by quality */
-      // XXX FIX THIS!!
-      // if (DophotSelect && (catalog[0].measure[offset].dophot != DophotValue)) continue;
-      
-      /* select measurements by quality */
-      if (PhotFlagSelect) {
-        if (PhotFlagBad) {
-          mask = PhotFlagBad;
-        } else {
-          code = GetPhotcodebyCode (catalog[0].measure[offset].photcode);
-          mask = code[0].astromBadMask;
-        }
-        if (mask & catalog[0].measure[offset].photFlags) {
-	  skipPhotFlag++;
-	  continue;
-	}
-      }
-      
-      /* select measurements by measurement error */
-      if ((SIGMA_LIM > 0) && (catalog[0].measure[offset].dM > SIGMA_LIM)) {
-	skipSigmaLim++;
-	continue;
-      }
-      
-      /* select measurements by mag limit */
-      if (ImagSelect) {
-        mag = PhotInst (&catalog[0].measure[offset]);
-        if (mag < ImagMin || mag > ImagMax) {
-	  skipImagSelect++;
-	  continue;
-	}
-      }
       // re-assess on each run of relastro if a measurement should be used
 
@@ -162,12 +98,4 @@
 
   if (VERBOSE) {
-    fprintf(stderr, "Reasons for exclusion in bcatalog:\n");
-    fprintf(stderr, "\ntoo few: %d \nphotCodeSkip: %d \nphotCodeKeep: %d \ntime: %d\n",
-	    skipfew, skipPhotCodeKeep, skipPhotCodeSkip, skipTime);
-    fprintf(stderr, "photFlag: %d \nmagnitude error: %d \nImag: %d\n",
-	    skipPhotFlag, skipSigmaLim, skipImagSelect);
-  }
-
-  if (VERBOSE) {
     fprintf (stderr, "%d: using %d stars (%d measures) for catalog\n", i,
              subcatalog[0].Naverage, subcatalog[0].Nmeasure);
Index: /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/load_catalogs.c
===================================================================
--- /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/load_catalogs.c	(revision 23648)
+++ /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/load_catalogs.c	(revision 23649)
@@ -29,4 +29,9 @@
     }
     if (VERBOSE && !pcatalog[0].Naves_disk) fprintf (stderr, "no data in %s, skipping\n", pcatalog[0].filename);
+
+    //outlier rejection
+    if (FlagOutlier) {
+      FlagOutliers(&tcatalog, 1);
+    }
 
     // select only the brighter stars
Index: /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/relastro_objects.c
===================================================================
--- /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/relastro_objects.c	(revision 23648)
+++ /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/relastro_objects.c	(revision 23649)
@@ -37,4 +37,8 @@
     }
 
+    if (FlagOutliers) {
+      FlagOutliers(&catalog, 1);
+    }
+
     // XXX consider what gets reset (only ASTROM flags)
     if (RESET) {