Index: trunk/Ohana/src/relphot/src/args.c
===================================================================
--- trunk/Ohana/src/relphot/src/args.c	(revision 30616)
+++ trunk/Ohana/src/relphot/src/args.c	(revision 31160)
@@ -78,4 +78,5 @@
     remove_argument (N, &argc, argv);
     PLOTDELAY = 1e6*atof(argv[N]);
+    PLOTSTUFF = TRUE; // always turn on plotting if i request a plot delay
     remove_argument (N, &argc, argv);
   }
Index: trunk/Ohana/src/relphot/src/bcatalog.c
===================================================================
--- trunk/Ohana/src/relphot/src/bcatalog.c	(revision 30616)
+++ trunk/Ohana/src/relphot/src/bcatalog.c	(revision 31160)
@@ -7,5 +7,5 @@
   off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm;
   float mag;
-  int Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew;
+  int Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Ngalaxy, Npsfqf;
 
   // XXX PhotNsec as a global is a bad idea; either get it from catalog
@@ -23,5 +23,5 @@
   Nmeasure = Naverage = 0;
 
-  Ncode = Ntime = Ndophot = Nmag = Nsigma = Nimag = Nfew = 0;
+  Ncode = Ntime = Ndophot = Nmag = Nsigma = Nimag = Nfew = Npsfqf = Ngalaxy = 0;
 
   /* exclude stars not in range or with too few measurements */
@@ -44,4 +44,6 @@
 
     Nm = 0;
+    int nEXT = 0;
+    int nPSF = 0;
     for (j = 0; j < catalog[0].average[i].Nmeasure; j++) {
 
@@ -62,4 +64,16 @@
       // XXX chnage this to select by bitflags
       if (DophotSelect && ((catalog[0].measure[offset].photFlags >> 16) != DophotValue)) { Ndophot ++; continue; }
+
+      // skip garbage measurements
+      if (catalog[0].measure[offset].psfQual < 0.85) { Npsfqf ++; continue; }
+
+      // check for galaxies
+      if (!isnan(catalog[0].measure[offset].Map)) {
+	  if (catalog[0].measure[offset].M - catalog[0].measure[offset].Map > 0.15) {
+	      nEXT ++;
+	  } else {
+	      nPSF ++;
+	  }
+      }
 
       /* select measurements by mag limit */
@@ -95,4 +109,11 @@
     }
 
+    // skip object if it is likely to be a galaxy
+    if (nEXT >= nPSF) {
+      Nmeasure -= Nm;
+      Ngalaxy ++;
+      continue; 
+    }
+
     // XXXX test : what checks do I need to make elsewhere to avoid problems here?
     if (Nm <= STAR_TOOFEW) { /* enough measurements in band? */
@@ -121,6 +142,6 @@
     fprintf (stderr, "using "OFF_T_FMT" stars ("OFF_T_FMT" measures) of "OFF_T_FMT" for catalog %s\n", 
 	     subcatalog[0].Naverage,  subcatalog[0].Nmeasure,  i, catalog[0].filename);
-    fprintf (stderr, "rejections: %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d few\n", 
-	     Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew);
+    fprintf (stderr, "rejections: %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d few, %d psfqf, %d galaxies\n", 
+	     Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Npsfqf, Ngalaxy);
   }
   return (TRUE);
Index: trunk/Ohana/src/relphot/src/load_images.c
===================================================================
--- trunk/Ohana/src/relphot/src/load_images.c	(revision 30616)
+++ trunk/Ohana/src/relphot/src/load_images.c	(revision 31160)
@@ -7,4 +7,6 @@
   fprintf (stderr, MSG, __VA_ARGS__); }
 
+// This function generates a subset of the images based on selections.  Input db has already
+// been loaded with the raw fits table data
 SkyList *load_images (FITS_DB *db, char *regionName, SkyRegion *region, int RegionSelect) {
 
@@ -32,5 +34,5 @@
   }
 
-  // convert database table to internal structure
+  // convert database table to internal structure (binary to Image)
   image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped);
   if (!image) {
@@ -44,10 +46,15 @@
   MARKTIME("selected images: %f sec\n", dtime);
 
+  // generate db->vtable from db->ftable based on the selection
+  // XXX does this simply duplicate the memory needlessly?  we recreate these lines
+  // in reload_images.  If we had saved the line numbers, we could avoid this
   gfits_vtable_from_ftable (&db[0].ftable, &db[0].vtable, LineNumber, Nsubset);
   MARKTIME("converted ftable to vtable: %f sec\n", dtime);
 
+  // save the subset of images in the static reference in ImageOps, set up indexes
   initImages (subset, Nsubset);
   MARKTIME("init images: %f sec\n", dtime);
 
+  // match chips to mosaics (if applicable)
   initMosaics (subset, Nsubset);
   MARKTIME("init mosaics: %f sec\n", dtime);
Index: trunk/Ohana/src/relphot/src/select_images.c
===================================================================
--- trunk/Ohana/src/relphot/src/select_images.c	(revision 30616)
+++ trunk/Ohana/src/relphot/src/select_images.c	(revision 31160)
@@ -33,5 +33,5 @@
   struct timeval start, stop;
   
-  double RmaxSkyRegion, RminSkyRegion, DminSkyRegion, DmaxSkyRegion, RmidSkyRegion;
+  double RmaxSkyRegion, RminSkyRegion, RmidSkyRegion, DminSkyRegion, DmaxSkyRegion;
 
   double *RmaxSky;
@@ -131,13 +131,21 @@
     }
 
-    /* define image corners */
-    Xi[0] = 0;            Yi[0] = 0;
-    Xi[1] = timage[i].NX; Yi[1] = 0;
-    Xi[2] = timage[i].NX; Yi[2] = timage[i].NY;
-    Xi[3] = 0;            Yi[3] = timage[i].NY;
-    Xi[4] = 0;            Yi[4] = 0;
+    /* define image corners - note the DIS images (mosaic phu) are special */
+    if (!strcmp(&timage[i].coords.ctype[4], "-DIS")) {
+      Xi[0] = -0.5*timage[i].NX; Yi[0] = -0.5*timage[i].NY;
+      Xi[1] = +0.5*timage[i].NX; Yi[1] = -0.5*timage[i].NY;
+      Xi[2] = +0.5*timage[i].NX; Yi[2] = +0.5*timage[i].NY;
+      Xi[3] = -0.5*timage[i].NX; Yi[3] = +0.5*timage[i].NY;
+      Xi[4] = -0.5*timage[i].NX; Yi[4] = -0.5*timage[i].NY;
+    } else {
+      Xi[0] = 0;            Yi[0] = 0;
+      Xi[1] = timage[i].NX; Yi[1] = 0;
+      Xi[2] = timage[i].NX; Yi[2] = timage[i].NY;
+      Xi[3] = 0;            Yi[3] = timage[i].NY;
+      Xi[4] = 0;            Yi[4] = 0;
+    }
     found = FALSE;
 
-    /* transform corners to ra,dec -- costs ~3sec for 3M images */
+    /* transform corners to ra,dec -- costs ~3sec for 3M images (pikake) */
     double RminImage = 360.0;
     double RmaxImage =   0.0;
@@ -147,5 +155,5 @@
       XY_to_RD (&Ri[j], &Di[j], Xi[j], Yi[j], &timage[i].coords);
       Ri[j] = ohana_normalize_angle_to_midpoint (Ri[j], RmidSkyRegion);
-      
+
       RminImage = MIN(RminImage, Ri[j]);
       RmaxImage = MAX(RmaxImage, Ri[j]);
@@ -158,5 +166,5 @@
 	RminImage = tmp - 360.0;
     }
-    
+
     // check that this image is even in range of the searched region
     if (DminImage > DmaxSkyRegion) continue;
@@ -182,5 +190,5 @@
       tcoords.crval2 = skycoords[m].Dc;
 
-      /* transform to ra,dec */
+      /* transform corner coords to X,Y in this catalog system */
       InRange = TRUE;
       for (j = 0; (j < 5) && InRange; j++) {
@@ -352,21 +360,2 @@
   return (Nlo);
 }
-
-off_t getRegionStopByRA (double R, double *Rref, off_t Nregions) {
-
-  // use bisection to find the overlapping mosaic
-
-  off_t Nlo, Nhi, N;
-
-  // find the last mosaic before start
-  Nlo = 0; Nhi = Nregions;
-  while (Nhi - Nlo > 10) {
-    N = 0.5*(Nlo + Nhi);
-    if (Rref[N] < R) {
-      Nlo = MAX(N, 0);
-    } else {
-      Nhi = MIN(N, Nregions);
-    }
-  }
-  return (Nlo);
-}