Changeset 30616 for trunk/Ohana/src/relastro

Timestamp:

Feb 13, 2011, 11:30:02 AM (15 years ago)

Author:

eugene

Message:

various fixes to photdbc, relastro and relphot based on analysis of the first reference database

Location:

trunk/Ohana/src/relastro

Files:

• include/relastro.h (modified) (6 diffs)
• src/ConfigInit.c (modified) (1 diff)
• src/FitChip.c (modified) (6 diffs)
• src/GetAstromError.c (modified) (3 diffs)
• src/ImageOps.c (modified) (17 diffs)
• src/StarMaps.c (modified) (3 diffs)
• src/UpdateChips.c (modified) (3 diffs)
• src/UpdateObjectOffsets.c (modified) (1 diff)
• src/UpdateObjects.c (modified) (2 diffs)
• src/args.c (modified) (3 diffs)
• src/bcatalog.c (modified) (2 diffs)
• src/load_catalogs.c (modified) (1 diff)
• src/plotstuff.c (modified) (1 diff)
• src/relastro.c (modified) (1 diff)
• src/relastroVisual.c (modified) (21 diffs)
• src/select_images.c (modified) (8 diffs)

trunk/Ohana/src/relastro/include/relastro.h

@@ -26 +26 @@
   double dPos;
   int mask;
+  int Nmeas;
 } StarData;
 
@@ -94 +95 @@
 
 double SIGMA_LIM;
-int SRC_MEAS_TOOFEW; //catalog objects wich fewer detections then this are ignored
+int    SRC_MEAS_TOOFEW; //catalog objects wich fewer detections then this are ignored
 double MIN_ERROR;
+
+int    IMFIT_CLIP_NITER; // number of clipping iterations to perform in FitChip
+double IMFIT_CLIP_NSIGMA; // number of sigma to clip in FitChip
+double IMFIT_SYS_SIGMA_LIM; // max dMag for objects used to measure systematic scatter
 
 double RADIUS; // match radius for high-speed objects
@@ -141 +146 @@
 int FlagOutlier;
 int    CLIP_THRESH;
+int USE_BASIC_CHECK;
 
 FitMode FIT_MODE;
@@ -296 +302 @@
 void fixImageRaw (Catalog *catalog, int Ncatalog, off_t im);
 void FlagOutliers(Catalog *catalog);
-int MeasFilterTest(Measure *measure);
+int MeasFilterTest(Measure *measure, int applySigmaLim);
 
 int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon);
@@ -325 +331 @@
 int UpdateObjectOffsets (SkyList *skylist);
 
-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);
-
-
+int relastroVisualPlotChipFit(StarData *raw, StarData *ref, double dRmax, int numObj);
+// 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);
+void relastroSetVisual(int state);
+int relastroGetVisual(void);
 
 int FixProblemImages (SkyList *skylist);
@@ -352 +358 @@
 int createStarMapPoints();
 int checkStarMap(int N);
+
+int GetScatterRawRef(float *dLsig, float *dMsig, float *dRsig, int *nKeep, StarData *raw, StarData *ref, int Nstars, float SigmaLimit);

trunk/Ohana/src/relastro/src/ConfigInit.c

@@ -19 +19 @@
   if (VERBOSE) fprintf (stderr, "loaded config file: %s\n", file);
 
-  GetConfig (config, "RELASTRO_SIGMA_LIM",     "%lf", 0, &SIGMA_LIM);
+  GetConfig (config, "RELASTRO_SIGMA_LIM",         "%lf", 0, &SIGMA_LIM); // exclude measurements on this basis
   GetConfig (config, "RELASTRO_SRC_MEAS_TOOFEW",   "%d",  0, &SRC_MEAS_TOOFEW);
+
+  if (!ScanConfig (config, "RELASTRO_IMFIT_CLIP_NITER",    "%d",  0, &IMFIT_CLIP_NITER))    IMFIT_CLIP_NITER  = 3;
+  if (!ScanConfig (config, "RELASTRO_IMFIT_CLIP_NSIGMA",   "%lf", 0, &IMFIT_CLIP_NSIGMA))   IMFIT_CLIP_NSIGMA = 3.0;
+  if (!ScanConfig (config, "RELASTRO_IMFIT_SYS_SIGMA_LIM", "%lf", 0, &IMFIT_SYS_SIGMA_LIM)) IMFIT_SYS_SIGMA_LIM = 0.01;
 
   // XXX these are used in relphot to identify poor stars / images -- define

trunk/Ohana/src/relastro/src/FitChip.c

@@ -2 +2 @@
 # include "relastroVisual.h"
 
-// XXX make these user parameters
-# define FIT_CHIP_NITER     3
-# define FIT_CHIP_NSIGMA    3.0
-
-// XXX save the fit[0].Npts value in the image table?
-
-// XXX save measurements of the fit quality (scatter, chisq) in the image table
-
 int FitChip (StarData *raw, StarData *ref, int Nmatch, Image *image) {
 
-  int i, Nscatter, Niter, skip;
-  CoordFit *fit;
-  double dL, dM, dR, dRmax, *values;
-
-  ALLOCATE (values, double, Nmatch);
-  for (Niter = 0; Niter < FIT_CHIP_NITER; Niter ++) {
+  int i, NstatFull, Nstat, Niter, skip;
+  float dLsig, dMsig, dRsig;
+  float dLsigFull, dMsigFull, dRsigFull;
+  float dL, dM, dR, dRmax;
+
+  CoordFit *fit = NULL;
+
+  dRmax = 0.0;
+
+  for (Niter = 0; Niter < IMFIT_CLIP_NITER; Niter ++) {
+    
+    // measure the scatter for the unmarked (good) measurements with dM < limit
+    // SIGMA_LIM here is redundant with ImageOps.c:915
+    GetScatterRawRef(&dLsig, &dMsig, &dRsig, &Nstat, raw, ref, Nmatch, SIGMA_LIM);
+    dRmax = IMFIT_CLIP_NSIGMA*dRsig;
+
+    // (a) skip unfittable points
+    // (b) mark good and bad points for fit (in and outliers)
 
     // measure the scatter distribution (use only the bright end detections)
-    for (i = Nscatter = 0; i < Nmatch; i++) {
+    for (i = 0; i < Nmatch; i++) {
       if (raw[i].mask) continue;
-      if (isnan(raw[i].dMag) || raw[i].dMag > SIGMA_LIM) continue;
+      if (isnan(raw[i].dMag)) { 
+        raw[i].mask |= 0x0004; 
+        continue; 
+      }
+      if (raw[i].dMag > SIGMA_LIM) { 
+        raw[i].mask |= 0x0008; 
+        continue; 
+      } // is this redundant with ImageOps.c:915?
 
       dL = raw[i].L - ref[i].L;
       dM = raw[i].M - ref[i].M;
       dR = hypot (dL, dM);
-
-      values[Nscatter] = dR;
-      Nscatter++;
-    }
-
-    if (Nscatter > 5) {
-      // for a 2D Gaussian, 40% of the points are within 1 sigma; dRmax is ~ 3 sigma
-      // XXX this test is not sensible for Nscatter < XXX (5?)
-      dsort (values, Nscatter);
-      dRmax = FIT_CHIP_NSIGMA*values[(int)(0.40*Nscatter)];
-      relastroVisualPlotScatter(values, dRmax, Nscatter);
-      relastroVisualPlotRawRef(raw, ref, dRmax, Nmatch);
-    } else {
-      dRmax = 0;
-    }
+      if (dR > dRmax) {
+        raw[i].mask |= 0x0010;
+        continue;
+      }
+    }
+
+    // figures to assess the fitting process:
+    // x vs dx, x vs dy, y vs dx, y vs dy : 
+    // residual vector field 
+    // dx vs mag, dy vs mag
+    relastroVisualPlotChipFit(raw, ref, dRmax, Nmatch);
 
     // fit the requested order polynomial
@@ -47 +54 @@
       image[0].coords.Npolyterms = CHIPORDER;
     }
+    if (fit) fit_free (fit);
     fit = fit_init (image[0].coords.Npolyterms);
 
     // generate the fit matches
     for (i = 0; i < Nmatch; i++) {
-      if (raw[i].mask) {
-        continue;
-      }
-      if (isnan(raw[i].dMag) || raw[i].dMag > SIGMA_LIM) {
-        continue;
-      }
-
-      // only keep objects within dRmax
-      dL = raw[i].L - ref[i].L;
-      dM = raw[i].M - ref[i].M;
-      dR = hypot (dL, dM);
-
-      // fprintf (stderr, "fit %f %f -> %f %f : %f %f (%f vs %f) wt: %f\n", raw[i].X, raw[i].Y, ref[i].L, ref[i].M, raw[i].L, raw[i].M, dR, dRmax, raw[i].dPos);
-
-      if ((dRmax > 0.0) && (dR > dRmax)) continue;
-
+      if (raw[i].mask) continue;
       fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, raw[i].dPos);
     }
@@ -88 +81 @@
     }
     if (skip) {
-      if (VERBOSE) fprintf (stderr, "insufficient measurements (%d) for requested order (%d)\n", fit[0].Npts, image[0].coords.Npolyterms);
+      if (VERBOSE2) fprintf (stderr, "insufficient measurements (%d) for requested order (%d)\n", fit[0].Npts, image[0].coords.Npolyterms);
       fit_free (fit);
-      free (values);
       image[0].flags |= ID_IMAGE_ASTROM_FEW;
       return FALSE;
     }
-
-    // fprintf (stderr, "scatter limit: %f based on %d detections; using %d of %d for fit\n", dRmax, Nscatter, fit[0].Npts, Nmatch);
 
     // measure the fit, update the coords & object coordinates
@@ -110 +100 @@
     }
 
-    fit_free (fit);
-
     for (i = 0; i < Nmatch; i++) {
       XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[0].coords);
@@ -118 +106 @@
   }
 
-  free (values);
+  // count mask classes
+  int nMask1 = 0;
+  int nMask2 = 0;
+  int nMask3 = 0;
+  int nMask4 = 0;
+  int nMask5 = 0;
+  for (i = 0; i < Nmatch; i++) {
+    if (!raw[i].mask) continue;
+    if (raw[i].mask & 0x01) nMask1 ++;
+    if (raw[i].mask & 0x02) nMask2 ++;
+    if (raw[i].mask & 0x04) nMask3 ++;
+    if (raw[i].mask & 0x08) nMask4 ++;
+    if (raw[i].mask & 0x10) nMask5 ++;
+  }
+
+  GetScatterRawRef(&dLsigFull, &dMsigFull, &dRsigFull, &NstatFull, raw, ref, Nmatch, SIGMA_LIM);
+  GetScatterRawRef(&dLsig, &dMsig, &dRsig, &Nstat, raw, ref, Nmatch, IMFIT_SYS_SIGMA_LIM);
+
+  int Nm = 0;
+  int Ns = 0;
+  for (i = 0; i < Nmatch; i++) {
+    if (raw[i].mask) continue;
+    Nm += raw[i].Nmeas;
+    Ns++;
+  }
+  image[0].nLinkAstrom = (Nm / Ns);
+
+  if (VERBOSE) fprintf (stderr, "fit sigma: %f (%f, %f) : full: %f (%f, %f), scatter limit: %f (%d full, %d bright, %d fit, %d all) (%d %d %d %d %d)\n", dRsig, dLsig, dMsig, dRsigFull, dLsigFull, dMsigFull, dRmax, NstatFull, Nstat, fit[0].Npts, Nmatch, nMask1, nMask2, nMask3, nMask4, nMask5);
+
+  image[0].dXpixSys = dLsig;
+  image[0].dYpixSys = dMsig;
+  image[0].nFitAstrom = fit[0].Npts;
+
+  if (fit) fit_free (fit);
+
   return TRUE;
+}
+
+// measure the scatter distribution (limit selected objects by dMag)
+int GetScatterRawRef(float *dLsig, float *dMsig, float *dRsig, int *nKeep, StarData *raw, StarData *ref, int Nstars, float SigmaLimit) {
+
+  int i, Ns;
+  float dL, dM, dLsum, dLsum2, dMsum, dMsum2, *dR;
+
+  ALLOCATE (dR, float, Nstars);
+
+  Ns = 0;
+  dLsum = dLsum2 = dMsum = dMsum2 = 0.0;
+  for (i = 0; i < Nstars; i++) {
+    if (raw[i].mask) continue;
+    if (isnan(raw[i].dMag)) continue;
+    if ((SigmaLimit > 0.0) && (raw[i].dMag > SigmaLimit)) continue;
+    
+    dL = raw[i].L - ref[i].L;
+    dM = raw[i].M - ref[i].M;
+
+    dLsum  += dL;
+    dLsum2 += dL*dL;
+    dMsum  += dM;
+    dMsum2 += dM*dM;
+
+    dR[Ns] = hypot (dL, dM);
+    Ns++;
+  }
+
+  *dLsig = sqrt((dLsum2 / Ns) - SQ(dLsum/Ns));
+  *dMsig = sqrt((dMsum2 / Ns) - SQ(dMsum/Ns));
+  *nKeep = Ns;
+
+  if (Ns < 5) {
+    *dRsig = NAN;
+    free  (dR);
+    return (FALSE);
+  }
+
+  // for a 2D Gaussian, 40% of the points are within R = 1 sigma
+  fsort (dR, Ns);
+  *dRsig = dR[(int)(0.40*Ns)];
+  
+  free  (dR);
+  return (TRUE);
 }
 
@@ -162 +229 @@
   fclose (f);
 */
+

trunk/Ohana/src/relastro/src/GetAstromError.c

@@ -1 +1 @@
 # include "relastro.h"
+# define WEIGHTED_ERRORS 1
 
 float GetAstromError (Measure *measure, int mode) {
-  //BIG HACKXXXXXXXX
-  return 0.1;
+
   PhotCode *code;
-  float dPobs, dPsys, dPtotal, dM, AS, MS;
+  float dPobs, dPsys, dPtotal, dM, AS, MS, dX, dY;
+
+  if (!WEIGHTED_ERRORS) {
+    // if we don't understand the errors at all, this at least lets us get things roughly
+    // right:
+    return 0.1;
+  }
+
   switch (mode) {
     case ERROR_MODE_RA:
-      dPobs = measure[0].dXccd / 100.0;  // need to redefine this as RAerr
+      dPobs = FromShortPixels(measure[0].dXccd);  // dXccd is a value in pixels
       break;
     case ERROR_MODE_DEC:
-      dPobs = measure[0].dYccd / 100.0;  // need to redefine this as RAerr
+      dPobs = FromShortPixels(measure[0].dYccd);  // dYccd is a value in pixels
       break;
     case ERROR_MODE_POS:
-      dPobs = hypot (measure[0].dXccd, measure[0].dYccd) / 100.0;  // need to redefine this as RAerr
+      dX = FromShortPixels(measure[0].dXccd);  // dXccd is a value in pixels
+      dY = FromShortPixels(measure[0].dYccd);  // dYccd is a value in pixels
+      dPobs = hypot (dX, dY);
       break;
     default:
       abort();
   }
-  /* the astrometric errors are not being carried yet (but should be!) */
-  /* we use the photometric mag error as a weighting term */
 
   code  = GetPhotcodebyCode (measure[0].photcode);
@@ -27 +34 @@
   dPsys = code[0].astromErrSys;
   dM    = measure[0].dM;
-  dPtotal = sqrt(SQ(dPsys) + AS*SQ(dPobs) + MS*SQ(dM));
-
-  //XXX dXccd, dYccd are now working correctly
-  //dPtotal = AS * dPobs;
+  dPtotal = sqrt(SQ(dPsys) + SQ(AS*dPobs) + SQ(MS*dM));
 
   dPtotal = MAX (dPtotal, MIN_ERROR);
@@ -36 +40 @@
 }
 
+/* for a long time, psphot was either not reported position errors, or was reporting
+ * completely wrong astrometry errors.  This function lets us handle, in the
+ * configuration, different strategies to generating a position error 
+ *
+ * astrometry systematic error : this is the minimum expected per-position error.  You
+ * should probably measure this from you data.  watch out for the chicken and egg problem!
+ *
+ * astrometry error scale : this field lets you accept position errors in (say) pixels and
+ * convert them with this term to arcsec.  AS : pixel scale in arcsec
+ *
+ * astrometry mag scale : this field lets you define position errors based on the
+ * photometry error.  the scale factor should be something like a typical seeing number
+ * (in arcsec) for the given instrument
+ *
+ */

trunk/Ohana/src/relastro/src/ImageOps.c

@@ -146 +146 @@
   free (clist);
   free (mlist);
+  free (Nlist);
+  free (NLIST);
 }
 
@@ -325 +327 @@
   for (j = 0; j < average[0].Nmeasure; j++) {
     off = average[0].measureOffset + j;
-    fprintf (stderr, "%f, %f\n", measure[off].dR, measure[off].dD);
+    fprintf (stderr, "dR, dD, mag, dMag: %f, %f, %f, %f\n", measure[off].dR, measure[off].dD, measure[off].M, measure[off].dM);
   }
   return;
@@ -336 +338 @@
   off_t i, m, c, n, nPos;
   double X, Y, L, M, P, Q, R, D, dR, dD;
-  double dPos, DPOS_MAX_ASEC;
+  double dPos, dPosSys, DPOS_MAX_ASEC;
 
   Mosaic *mosaic;
@@ -359 +361 @@
   }
 
-  // accumulate the rms position offsets.  if this value, or any specific entry, is too
-  // large, we will reset the image to the original coords at the end of the analysis
-
+  // these are used to accumulate the rms position offsets.  if this value, or any
+  // specific entry, is too large, we will reset the image to the original coords at the
+  // end of the analysis
   dPos = 0.0;
   nPos = 0;
+
+  // convert the image systematic error in pixels to a value in arcsec
+  { 
+    double dLsig, dMsig;
+    double Ro, Do, Rx, Dx, dP0, dP1;
+    Coords *coords;
+
+    // these values are in pixels, but we to convert to arcsec
+    dLsig = image[0].dXpixSys;
+    dMsig = image[0].dYpixSys;
+
+    if (moscoords == NULL) {
+      coords = imcoords;
+    } else {
+      coords = moscoords;
+    }
+    XY_to_LM (&Ro, &Do, 0.0, 0.0, coords);
+    XY_to_LM (&Rx, &Dx, dLsig, 0.0, coords);
+    dP0 = 3600.0 * hypot(Rx - Ro, Dx - Do); // convert to arcsec
+    XY_to_LM (&Rx, &Dx, 0.0, dLsig, coords);
+    dP1 = 3600.0 * hypot(Rx - Ro, Dx - Do); // convert to arcsec
+    dPosSys = 0.5 * (dP0 + dP1);
+  }      
 
   for (i = 0; i < Nlist[im]; i++) {
@@ -390 +415 @@
     // complain if the new location is far from the average location
     // NOTE: This should never happen, or our StarMap tests are not working
-    if (fabs(dR) > 1.5*ADDSTAR_RADIUS) {
+    if (fabs(dR) > 3.0*ADDSTAR_RADIUS) {
       fprintf (stderr, "measurement is far from average location (R): %f %f (%f %f)\n", catalog[c].average[n].R, catalog[c].average[n].D, dR, dD);
       dump_measures (&catalog[c].average[n], catalog[c].measure);
       // abort ();
     }
-    if (fabs(dD) > 1.5*ADDSTAR_RADIUS) {
+    if (fabs(dD) > 3.0*ADDSTAR_RADIUS) {
       fprintf (stderr, "measurement is far from average location (D): %f %f (%f %f)\n", catalog[c].average[n].R, catalog[c].average[n].D, dR, dD);
       dump_measures (&catalog[c].average[n], catalog[c].measure);
@@ -418 +443 @@
     catalog[c].measure[m].dR = dR;
     catalog[c].measure[m].dD = dD;
-
+    
     if (catalog[c].measure[m].dR > +180.0*3600.0) {
       // average on high end of boundary, move star up
@@ -429 +454 @@
       catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R);
     }
+
+    // set the systematic error for this image:
+    catalog[c].measure[m].dRsys = ToShortPixels(dPosSys);
   }
 
@@ -544 +572 @@
 
     // 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;
+    // XXX this is already taken care of in bcatalog
+    raw[i].mask = 0x0000;
     if (catalog[c].average[n].Nmeasure <= SRC_MEAS_TOOFEW) {
-      mask = TRUE;
+      raw[i].mask |= 0x0001;
     }
     if (!finite(catalog[c].measure[m].dR) || !finite(catalog[c].measure[m].dD)) {
-      mask = TRUE;
-    }
-
-    raw[i].mask = mask;
-
+      raw[i].mask |= 0x0002;
+    }
+    raw[i].Nmeas = catalog[c].average[n].Nmeasure; // record so we can check how well connected an image is
 
     switch (mode) {
@@ -685 +710 @@
       
       // skip measurements based on user selected criteria
-      if (!MeasFilterTest(&catalog[0].measure[m])) continue;
+      if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue;
       R[N] = catalog[0].measure[m].dR;
       D[N] = catalog[0].measure[m].dD;
@@ -711 +736 @@
 
       // skip measurements based on user selected criteria
-      if (!MeasFilterTest(&catalog[0].measure[m])) continue;
+      if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue;
       
       x = catalog[0].measure[m].dR - statsR.median;
@@ -725 +750 @@
     }
 
-    //examine results
-    relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset, 
-                               catalog[0].average[j].Nmeasure, 
-                               statsR, statsD, Ns);
+    // examine results
+    // relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset, catalog[0].average[j].Nmeasure, statsR, statsD, Ns);
   }
   
@@ -784 +807 @@
       
       // skip measurements based on user selected criteria
-      if (!MeasFilterTest(&catalog[0].measure[m])) continue;
+      if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue;
       R[N] = catalog[0].measure[m].dR;
       D[N] = catalog[0].measure[m].dD;
@@ -806 +829 @@
     for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
       //skip bad measurements
-      if (!MeasFilterTest(&catalog[0].measure[m])) continue;  
+      if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue;  
       x = catalog[0].measure[m].dR - statsR.median;
       y = catalog[0].measure[m].dD - statsD.median;
@@ -838 +861 @@
     for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
       //skip bad measurements
-      if (!MeasFilterTest(&catalog[0].measure[m])) continue;  
+      if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue;  
       x = catalog[0].measure[m].dR - statsR.median;
       y = catalog[0].measure[m].dD - statsD.median;
@@ -852 +875 @@
     }  //done rejecting outliers
 
-    //examine results
-    relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset, 
-                               catalog[0].average[j].Nmeasure, 
-                               statsR, statsD, Ns);
+    // examine results
+    // relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset, catalog[0].average[j].Nmeasure, statsR, statsD, Ns);
     
   } //done looping over objects
@@ -871 +892 @@
 
 /** Determine whether a measurement should be included in the analysis, based on supplied filter criteria */ 
-int MeasFilterTest(Measure *measure) {
+// we only optionally apply the sigma limit: for object averages, this should not be used (should it?)
+int MeasFilterTest(Measure *measure, int applySigmaLim) {
   int found, k;
   long mask;
@@ -917 +939 @@
 
   /* select measurements by measurement error */
-  if ((SIGMA_LIM > 0) && (measure[0].dM > SIGMA_LIM)) return FALSE;
+  if (applySigmaLim && (SIGMA_LIM > 0) && (measure[0].dM > SIGMA_LIM)) {
+    return FALSE;
+  }
   
   /* select measurements by mag limit */

trunk/Ohana/src/relastro/src/StarMaps.c

@@ -79 +79 @@
     starmap[N].stars[ybin*NX_MAP + xbin] ++;
   }
-  MARKTIME("assign stars to starmap bins: %f sec\n", dtime);
+  if (VERBOSE2) { MARKTIME("assign stars to starmap bins: %f sec\n", dtime); }
 
   return (TRUE);
@@ -117 +117 @@
     }
 
-    if (VERBOSE) fprintf (stderr, "starmap: %d points for image %s\n", starmap[i].Npoints, images[i].name);
+    if (VERBOSE2) fprintf (stderr, "starmap: %d points for image %s\n", starmap[i].Npoints, images[i].name);
   }
 
@@ -145 +145 @@
   }
 
-  if (VERBOSE) fprintf (stderr, "max deviations for %s using %d pts : %f, %f\n", images[N].name, starmap[N].Npoints, dLmax, dMmax);
+  if (VERBOSE2) fprintf (stderr, "max deviations for %s using %d pts : %f, %f\n", images[N].name, starmap[N].Npoints, dLmax, dMmax);
 
-  if (dLmax > DPOS_MAX) return (FALSE);
-  if (dMmax > DPOS_MAX) return (FALSE);
+  if (dLmax > DPOS_MAX) {
+      if (VERBOSE) fprintf (stderr, "max deviations for %s using %d pts : %f, %f\n", images[N].name, starmap[N].Npoints, dLmax, dMmax);
+      return (FALSE);
+  }
+  if (dMmax > DPOS_MAX) {
+      if (VERBOSE) fprintf (stderr, "max deviations for %s using %d pts : %f, %f\n", images[N].name, starmap[N].Npoints, dLmax, dMmax);
+      return (FALSE);
+  }
   return (TRUE);
 }

trunk/Ohana/src/relastro/src/UpdateChips.c

@@ -1 +1 @@
 # include "relastro.h"
+int plotChipFits (double *Ro, double *Do, char *mode, int Nimage);
+int saveCenter (Image *image, double *Ro, double *Do, int im);
 
 int UpdateChips (Catalog *catalog, int Ncatalog) {
 
+  int Nskip, Nmosaic, NnewFit, NoldFit;
+
   /* we can measure new image parameters for each non-mosaic chip independently */
-  off_t i, Nimage, Nraw, Nref;
+  off_t i, Nimage, Nraw, Nref, nFitAstr;
   Image *image;
   StarData *raw, *ref;
   Coords *oldCoords;
+  float dXpixSys, dYpixSys;
+  double *Ro, *Do;
+  char *mode;
+
+  Nskip = Nmosaic = NnewFit = NoldFit = 0;
 
   image = getimages (&Nimage);
 
+  // save fit results for summary plot
+  ALLOCATE (Ro, double, Nimage);
+  ALLOCATE (Do, double, Nimage);
+  ALLOCATE (mode, char, Nimage);
+
   for (i = 0; i < Nimage; i++) {
 
     /* skip all except WRP images */
-    if (strcmp(&image[i].coords.ctype[4], "-WRP")) continue;
+    if (strcmp(&image[i].coords.ctype[4], "-WRP")) {
+      Nmosaic ++;
+      mode[i] = 0;
+      continue;
+    }
 
     /* convert measure coordinates to raw entries */
     raw = getImageRaw (catalog, Ncatalog, i, &Nraw, MODE_MOSAIC);
-    if (!raw) continue;
+    if (!raw) {
+      Nskip ++;
+      mode[i] = 0;
+      continue;
+    }
 
     /* convert average coordinates to ref entries */
     ref = getImageRef (catalog, Ncatalog, i, &Nref, MODE_MOSAIC);
-    if (!ref) continue;
+    if (!ref) {
+      Nskip ++;
+      mode[i] = 0;
+      continue;
+    }
 
     // note that Nraw & Nref must be equal: if not, we made a programming error in one of these two functions.
     assert (Nraw == Nref);
 
+    // save these in case of failure
     saveCoords (&image[i].coords, i);
+    dXpixSys = image[i].dXpixSys;
+    dYpixSys = image[i].dYpixSys;
+    nFitAstr = image[i].nFitAstrom;
 
     // FitChip does iterative, clipped fitting
@@ -33 +63 @@
     if (!FitChip (raw, ref, Nraw, &image[i])) {
       if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT"\n", image[i].name,  i,  Nraw);
+
+      // restore status quo ante
       oldCoords = getCoords (i);
       memcpy (&image[i].coords, oldCoords, sizeof(Coords));
+      image[i].dXpixSys = dXpixSys;
+      image[i].dYpixSys = dYpixSys; 
+      image[i].nFitAstrom = nFitAstr;
+
+      saveCenter (image, &Ro[i], &Do[i], i);
+      mode[i] = 1;
+      NoldFit ++;
       free (raw);
       free (ref);
@@ -42 +81 @@
     if (!checkStarMap (i)) {
       if (VERBOSE) fprintf (stderr, "fit diverges too much for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT"\n", image[i].name,  i,  Nraw);
+      // restore status quo ante
       oldCoords = getCoords (i);
       memcpy (&image[i].coords, oldCoords, sizeof(Coords));
+      image[i].dXpixSys = dXpixSys;
+      image[i].dYpixSys = dYpixSys; 
+      image[i].nFitAstrom = nFitAstr;
+
+      saveCenter (image, &Ro[i], &Do[i], i);
+      mode[i] = 2;
       image[i].flags |= ID_IMAGE_ASTROM_POOR;
-    }
-
+      NoldFit ++;
+      free (raw);
+      free (ref);
+      continue;
+    } 
+
+    saveCenter (image, &Ro[i], &Do[i], i);
+    mode[i] = 3;
+    NnewFit ++;
     free (raw);
     free (ref);
   }
 
+  plotChipFits (Ro, Do, mode, Nimage);
+
+  fprintf (stderr, "UpdateChips: %d fitted, %d keep old, %d skipped, %d mosaic (skipped)\n", NnewFit, NoldFit, Nskip, Nmosaic);
   return (TRUE);
 }
 
+int saveCenter (Image *image, double *Ro, double *Do, int im) {
+
+  Mosaic *mosaic;
+  Coords *moscoords, *imcoords;
+  double X, Y, L, M, P, Q, R, D;
+
+  moscoords = NULL;
+  if (!strcmp(&image[im].coords.ctype[4], "-WRP")) {
+    mosaic = getMosaicForImage (im);
+    if (mosaic == NULL) return FALSE;  // if we cannot find the associated image, skip it
+    moscoords = &mosaic[0].coords;
+  }
+  imcoords = &image[im].coords;
+  
+  if (!strcmp(&image[im].coords.ctype[4], "-WRP")) {
+    X = 0.5*image[im].NX; 
+    Y = 0.5*image[im].NY;
+  } else {
+    X = 0.0; 
+    Y = 0.0;
+  }
+
+  if (moscoords == NULL) {
+    // this is a Simple image (not a mosaic)
+    // note that for a Simple image, L,M = P,Q
+    XY_to_LM (&L, &M, X, Y, imcoords);
+    LM_to_RD (&R, &D, L, M, imcoords);
+  } else {
+    XY_to_LM (&L, &M, X, Y, imcoords);
+    XY_to_LM (&P, &Q, L, M, moscoords);
+    LM_to_RD (&R, &D, P, Q, moscoords);
+  }
+
+  double Rmid;
+  if (UserCatalog) {
+      Rmid = UserCatalogRA;
+  } else {
+      Rmid = 0.5*(UserPatch.Rmin + UserPatch.Rmax);
+  }
+
+  R = ohana_normalize_angle_to_midpoint (R, Rmid);
+
+  *Ro = R;
+  *Do = D;
+
+  return (TRUE);
+}
+ 
+int plotChipFits (double *Ro, double *Do, char *mode, int Nimage) {
+
+  static int kapa = -1; 
+
+  int i, N;
+  double Rmin, Rmax, Dmin, Dmax;
+  float *xvec, *yvec;
+  Graphdata graphdata;
+
+  if (!relastroGetVisual()) return (TRUE);
+
+  if (kapa == -1) {
+    kapa = KapaOpenNamedSocket("kapa", "relastro");
+    if (kapa == -1) {
+      fprintf (stderr, "can't open kapa window\n");
+      return FALSE;
+    }
+  }
+
+  Rmin = +720;
+  Rmax = -720;
+  Dmin = +90;
+  Dmax = -90;
+
+  // find the R, D range
+  for (i = 0; i < Nimage; i++) {
+    if (!mode[i]) continue;
+
+    Rmin = MIN(Rmin, Ro[i]);
+    Rmax = MAX(Rmax, Ro[i]);
+    Dmin = MIN(Dmin, Do[i]);
+    Dmax = MAX(Dmax, Do[i]);
+  }
+
+  ALLOCATE (xvec, float, Nimage);
+  ALLOCATE (yvec, float, Nimage);
+
+  bzero (&graphdata, sizeof(Graphdata));
+  plot_defaults (&graphdata);
+  graphdata.xmin = Rmin;
+  graphdata.xmax = Rmax;
+  graphdata.ymin = Dmin;
+  graphdata.ymax = Dmax;
+  graphdata.style = 2;
+  graphdata.size = 1;
+
+  KapaSetFont (kapa, "helvetica", 14);
+  KapaSetLimits (kapa, &graphdata);
+  KapaBox (kapa, &graphdata);
+
+  // *** good images ***
+  N = 0;
+  for (i = 0; i < Nimage; i++) {
+    if (mode[i] != 3) continue;
+    xvec[N] = Ro[i];
+    yvec[N] = Do[i];
+    N++;
+  }
+  graphdata.ptype = 7;
+  graphdata.color = KapaColorByName("black");
+  KapaPrepPlot (kapa, N, &graphdata);
+  KapaPlotVector (kapa, N, xvec, "x");
+  KapaPlotVector (kapa, N, yvec, "y");
+  
+  // *** reject fit ***
+  N = 0;
+  for (i = 0; i < Nimage; i++) {
+    if (mode[i] != 1) continue;
+    xvec[N] = Ro[i];
+    yvec[N] = Do[i];
+    N++;
+  }
+  graphdata.ptype = 3;
+  graphdata.color = KapaColorByName("red");
+  KapaPrepPlot (kapa, N, &graphdata);
+  KapaPlotVector (kapa, N, xvec, "x");
+  KapaPlotVector (kapa, N, yvec, "y");
+  
+  // *** divergent fit ***
+  N = 0;
+  for (i = 0; i < Nimage; i++) {
+    if (mode[i] != 2) continue;
+    xvec[N] = Ro[i];
+    yvec[N] = Do[i];
+    N++;
+  }
+  graphdata.ptype = 2;
+  graphdata.color = KapaColorByName("blue");
+  KapaPrepPlot (kapa, N, &graphdata);
+  KapaPlotVector (kapa, N, xvec, "x");
+  KapaPlotVector (kapa, N, yvec, "y");
+  
+  free (xvec);
+  free (yvec);
+
+  return (TRUE);
+}
+
+// XXX if (!FindMosaicForImage (image, Nimage, i)) { }

trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c

@@ -40 +40 @@
     UpdateMeasures (&catalog, 1);
 
+    UpdateObjects (&catalog, 1);
+
     freeImageBins (1);
 

trunk/Ohana/src/relastro/src/UpdateObjects.c

@@ -103 +103 @@
 
         //does the measurement pass the supplied filtering constraints?
-        if (!MeasFilterTest(&catalog[i].measure[m])) {
+        if (!MeasFilterTest(&catalog[i].measure[m], FALSE)) {
           catalog[i].measure[m].dbFlags &= ~ID_MEAS_USED_OBJ;
           continue;
@@ -134 +134 @@
 
         // dX, dY : error in arcsec -- 
-        // dX[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_RA);
-        // dY[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_DEC);
-
-        dX[N] = 0.1;
-        dY[N] = 0.1;
+        dX[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_RA);
+        dY[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_DEC);
+
+        // add systematic error in quadrature, if desired
+        // only do this after the fit has converged (or you will never improve the poor images)
+        // if (INCLUDE_SYS_ERR) {
+        // float dRsys = FromShortPixels(catalog[i].measure[m].dRsys);
+        // dX[N] = hypot(dX[N], dRsys);
+        // dY[N] = hypot(dY[N], dRsys);
+        // }
+
+        // dX[N] = 0.1;
+        // dY[N] = 0.1;
+
         dT[N] = catalog[i].measure[m].dt;
 

trunk/Ohana/src/relastro/src/args.c

@@ -13 +13 @@
     remove_argument (N, &argc, argv);
     FIT_TARGET = TARGET_OBJECTS;
-
-    // check for object fitting modes (not valid for images)
-    if ((N = get_argument (argc, argv, "-pm"))) {
-        remove_argument (N, &argc, argv);
-        FIT_MODE = FIT_PM_ONLY;
-    }
-    if ((N = get_argument (argc, argv, "-par"))) {
-        remove_argument (N, &argc, argv);
-        FIT_MODE = FIT_PAR_ONLY;
-    }
-    if ((N = get_argument (argc, argv, "-pmpar"))) {
-        remove_argument (N, &argc, argv);
-        FIT_MODE = FIT_PM_AND_PAR;
-    }
-  }
+  }
+
+  // check for object fitting modes
+  if ((N = get_argument (argc, argv, "-pm"))) {
+    remove_argument (N, &argc, argv);
+    FIT_MODE = FIT_PM_ONLY;
+  }
+  if ((N = get_argument (argc, argv, "-par"))) {
+    remove_argument (N, &argc, argv);
+    FIT_MODE = FIT_PAR_ONLY;
+  }
+  if ((N = get_argument (argc, argv, "-pmpar"))) {
+    remove_argument (N, &argc, argv);
+    FIT_MODE = FIT_PM_AND_PAR;
+  }
+
   if ((N = get_argument (argc, argv, "-high-speed"))) {
     // XXX include a parallax / no-parallax option
@@ -93 +94 @@
       usage ();
     }
+  }
+
+  USE_BASIC_CHECK = FALSE;
+  if ((N = get_argument (argc, argv, "-basic-image-search"))) {
+    remove_argument (N, &argc, argv);
+    USE_BASIC_CHECK = TRUE;
   }
 
@@ -144 +151 @@
     VERBOSE = VERBOSE2 = TRUE;
     remove_argument (N, &argc, argv);
+  }
+
+  if ((N = get_argument (argc, argv, "-visual"))) {
+    remove_argument (N, &argc, argv);
+    relastroSetVisual(TRUE);
   }
 

trunk/Ohana/src/relastro/src/bcatalog.c

@@ -45 +45 @@
       offset = catalog[0].average[i].measureOffset + j;
       
-      // filter objects based on user supplied criteria
-      if (!MeasFilterTest(&catalog[0].measure[offset])) {
+      // filter objects based on user supplied criteria, including SIGMA_LIM
+      if (!MeasFilterTest(&catalog[0].measure[offset], TRUE)) {
         catalog[0].measure[offset].dbFlags &= ~ID_MEAS_USED_CHIP;
         continue;
@@ -78 +78 @@
         subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_AREA;
       }
-     
+
       Nmeasure ++;
       Nm ++;

trunk/Ohana/src/relastro/src/load_catalogs.c

@@ -26 +26 @@
     pcatalog[0].Nsecfilt  = GetPhotcodeNsecfilt ();
 
-    if (!dvo_catalog_open (pcatalog, skylist[0].regions[i], VERBOSE, "w")) {
+    if (!dvo_catalog_open (pcatalog, skylist[0].regions[i], VERBOSE2, "w")) {
       fprintf (stderr, "ERROR: failure reading catalog %s\n", pcatalog[0].filename);
       exit (1);
     }
-    if (VERBOSE && !pcatalog[0].Naves_disk) fprintf (stderr, "no data in %s, skipping\n", pcatalog[0].filename);
+    if (VERBOSE2 && !pcatalog[0].Naves_disk) fprintf (stderr, "no data in %s, skipping\n", pcatalog[0].filename);
 
     //outlier rejection

trunk/Ohana/src/relastro/src/plotstuff.c

@@ -146 +146 @@
   graphdata[0].ymax = dUNDEF;
    
+  graphdata[0].ticktextPad = NAN;
+  graphdata[0].labelPadXm = NAN;
+  graphdata[0].labelPadXp = NAN;
+  graphdata[0].labelPadYm = NAN;
+  graphdata[0].labelPadYp = NAN;
+  graphdata[0].padXm = NAN;
+  graphdata[0].padXp = NAN;
+  graphdata[0].padYm = NAN;
+  graphdata[0].padYp = NAN;
 }

trunk/Ohana/src/relastro/src/relastro.c

@@ -78 +78 @@
         MARKTIME("update chips: %f sec\n", dtime);
       }
+      // create summary plots of the process
+      // relastroVisualSummaryChips();
       break;
 

trunk/Ohana/src/relastro/src/relastroVisual.c

@@ -11 +11 @@
 #define KAPAY 700
 
-static int kapa = -1;
+static int kapa1 = -1;
 static int kapa2 = -1;
 static int kapa3 = -1;
+static int kapa4 = -1;
 
 static int isVisual = FALSE;
-static int plotRawRef = FALSE;
-static int plotScatter = FALSE;
-static int plotResid = FALSE;
-static int plotVector = FALSE;
+static int plotRawRef = TRUE;
+// static int plotScatter = TRUE;
+// static int plotResid = TRUE;
+// static int plotVector = TRUE;
 static int plotOutliers = TRUE;
 
@@ -29 +30 @@
       fprintf(stderr, "Failure to open kapa.\n");
       isVisual = 0;
-      return 0;
+      return FALSE;
     }
-    //    KapaResize (*kapid, KAPAX, KAPAY);
-  }
-  return 1;
+  }
+  return TRUE;
 }
 
@@ -49 +49 @@
     isVisual = 0;
   }
-  return 1;
+  return TRUE;
+}
+
+void relastroSetVisual(int state) {
+  isVisual = state;
+}
+
+int relastroGetVisual(void) {
+  return isVisual;
 }
 
@@ -70 +78 @@
     graphdata->xmax = xhi;
     graphdata->ymax = yhi;
-    return 1;
-}
-
-static int residPlot(float x[], float y[],
-                     float xVec[], float yVec[],
-                     int npts, int *kapaID) {
-    if (!isVisual || !plotResid) return TRUE;
+    return TRUE;
+}
+
+/** 4-panel plot of x vs dx, y vs dx, x vs dy, y vs dy **/
+int relastroVisualRawRef(int *kapaID, float *x, float *y, float *dx, float *dy, float *dPos, int npts) {
+
     Graphdata graphdata;
     KapaSection section;
@@ -83 +90 @@
 
     KapaInitGraph(&graphdata);
-    KapaClearPlots(*kapaID);
+    KapaClearSections(*kapaID);
     KapaSetFont(*kapaID, "helvetica", 14);
 
@@ -89 +96 @@
     section.x = 0.0; section.y = 0.0;
     section.dx = .45, section.dy = .45;
+    section.bg = KapaColorByName("white");
     graphdata.ptype = 7;
     graphdata.style = 2;
+    graphdata.etype |= 0x01;
 
     KapaSetSection(*kapaID, &section);
-    if(!scaleGraphdata(x, xVec, &graphdata, npts)) return 0;
+    if(!scaleGraphdata(x, dx, &graphdata, npts)) return 0;
     KapaSetLimits(*kapaID, &graphdata);
     KapaBox(*kapaID, &graphdata);
@@ -100 +109 @@
     KapaPrepPlot(*kapaID, npts, &graphdata);
     KapaPlotVector(*kapaID, npts, x, "x");
-    KapaPlotVector(*kapaID, npts, xVec, "y");
+    KapaPlotVector(*kapaID, npts, dx, "y");
+    KapaPlotVector(*kapaID, npts, dPos, "dym");
+    KapaPlotVector(*kapaID, npts, dPos, "dyp");
 
     section.x = .5; section.y = 0; section.name="1";
     KapaSetSection(*kapaID, &section);
-    if(!scaleGraphdata(x, yVec, &graphdata, npts)) return 0;
+    if(!scaleGraphdata(x, dy, &graphdata, npts)) return 0;
     KapaSetLimits(*kapaID, &graphdata);
     KapaBox(*kapaID, &graphdata);
@@ -111 +122 @@
     KapaPrepPlot(*kapaID, npts, &graphdata);
     KapaPlotVector(*kapaID, npts, x, "x");
-    KapaPlotVector(*kapaID, npts, yVec, "y");
+    KapaPlotVector(*kapaID, npts, dy, "y");
+    KapaPlotVector(*kapaID, npts, dPos, "dym");
+    KapaPlotVector(*kapaID, npts, dPos, "dyp");
 
     section.x = .0; section.y = .5; section.name="2";
     KapaSetSection(*kapaID, &section);
-    if(!scaleGraphdata(y, xVec, &graphdata, npts)) return 0;;
+    if(!scaleGraphdata(y, dx, &graphdata, npts)) return 0;;
     KapaSetLimits(*kapaID, &graphdata);
     KapaBox(*kapaID, &graphdata);
@@ -122 +135 @@
     KapaPrepPlot(*kapaID, npts, &graphdata);
     KapaPlotVector(*kapaID, npts, y, "x");
-    KapaPlotVector(*kapaID, npts, xVec, "y");
+    KapaPlotVector(*kapaID, npts, dx, "y");
+    KapaPlotVector(*kapaID, npts, dPos, "dym");
+    KapaPlotVector(*kapaID, npts, dPos, "dyp");
 
     section.x = .5; section.y = .5; section.name="3";
     KapaSetSection(*kapaID, &section);
-    if(!scaleGraphdata(y, yVec, &graphdata, npts)) return 0;
+    if(!scaleGraphdata(y, dy, &graphdata, npts)) return 0;
     KapaSetLimits(*kapaID, &graphdata);
     KapaBox(*kapaID, &graphdata);
@@ -133 +148 @@
     KapaPrepPlot(*kapaID, npts, &graphdata);
     KapaPlotVector(*kapaID, npts, y, "x");
-    KapaPlotVector(*kapaID, npts, yVec, "y");
-
-    return 1;
-}
-
-
-/**Plot a vector field*/
-static int plotVectorField(float x[], float y[],
-                           float xVec[], float yVec[],
-                           int npts, int *kapaID, double maxVecLength) {
-
-    if(!plotVector) return 1;
+    KapaPlotVector(*kapaID, npts, dy, "y");
+    KapaPlotVector(*kapaID, npts, dPos, "dym");
+    KapaPlotVector(*kapaID, npts, dPos, "dyp");
+
+    return TRUE;
+}
+
+
+/** Plot a vector field (circles at vector origin, scaled lines giving vector directions */
+int relastroVisualVectorField(int *kapaID, float *x, float *y, float *dx, float *dy, int npts, double maxVecLength) {
 
     Graphdata graphdata;
-    float singleX[2], singleY[2];
+    float *xVec, *yVec;
     float vecScaleFactor;
     float graphSize;
     int i;
     char plotTitle[50];
-    sprintf(plotTitle, "Maximum Vector Size = %5.1e", maxVecLength);
-
+
+    if (!npts) return FALSE;
     if (!initWindow(kapaID)) return 0;
 
+    snprintf(plotTitle, 50, "Maximum Vector Size = %5.1e", maxVecLength);
+
     KapaInitGraph(&graphdata);
-    KapaClearPlots(*kapaID);
-    if(!scaleGraphdata(x, y, &graphdata, npts)) return 0;
+    KapaClearSections(*kapaID);
+    KapaSetFont(*kapaID, "helvetica", 14);
+
+    if (!scaleGraphdata(x, y, &graphdata, npts)) return 0;
 
     graphSize = graphdata.xmax - graphdata.xmin;
@@ -164 +181 @@
         graphSize = graphdata.ymax - graphdata.ymin;
     }
-
     vecScaleFactor = graphSize * 0.02 / (float)maxVecLength;
-#ifdef TESTING
-    fprintf(stderr, "GraphSize: %e\n maxVecLength: %e\n vecScaleFactor: %e\n",
-            graphSize, maxVecLength, vecScaleFactor);
-#endif
-
+
+    // fprintf(stderr, "GraphSize: %e\n maxVecLength: %e\n vecScaleFactor: %e\n", graphSize, maxVecLength, vecScaleFactor);
 
     KapaSetFont (*kapaID, "helvetica", 14);
@@ -179 +192 @@
     graphdata.ptype = 7;
     graphdata.style = 2;
+    graphdata.color = KapaColorByName("black");
     KapaPrepPlot(*kapaID, npts, &graphdata);
     KapaPlotVector(*kapaID, npts, x, "x");
     KapaPlotVector(*kapaID, npts, y, "y");
 
-    //plot each vector individually
-    graphdata.ptype = 0;
-    graphdata.style = 0;
+    ALLOCATE (xVec, float, 2*npts);
+    ALLOCATE (yVec, float, 2*npts);
+    for (i = 0; i < npts; i++) {
+      xVec[2*i + 0] = x[i];
+      yVec[2*i + 0] = y[i];
+      xVec[2*i + 1] = x[i] + dx[i] * vecScaleFactor;
+      yVec[2*i + 1] = y[i] + dy[i] * vecScaleFactor;
+    }
+
+    graphdata.ptype = 100; // line segements by point pair
+    graphdata.style = 2;
     graphdata.color = KapaColorByName("blue");
-    for(i = 0; i < npts; i++) {
-        singleX[0] = x[i];
-        singleY[0] = y[i];
-        singleX[1] = x[i] + xVec[i] * vecScaleFactor;
-        singleY[1] = y[i] + yVec[i] * vecScaleFactor;
-        KapaPrepPlot(*kapaID, 2, &graphdata);
-        KapaPlotVector(*kapaID, 2, singleX, "x");
-        KapaPlotVector(*kapaID, 2, singleY, "y");
-    }
-    return 1;
-}
-
-int relastroVisualPlotScatter(double values[], double thresh, int npts) {
-    float *x, *data;
-    int i;
-    float xline[2], yline[2];
-    if (!isVisual || !plotScatter) return 1;
-    if (!initWindow(&kapa2)) return 0;
-
-    ALLOCATE(x, float, npts);
-    ALLOCATE(data, float, npts);
-
-    for(i = 0; i < npts; i++) {
-        x[i] = i;
-        data[i] = (float) values[i];
-    }
+    KapaPrepPlot  (*kapaID, 2*npts, &graphdata);
+    KapaPlotVector(*kapaID, 2*npts, xVec, "x");
+    KapaPlotVector(*kapaID, 2*npts, yVec, "y");
+
+    free (xVec);
+    free (yVec);
+    return TRUE;
+}
+
+int relastroVisualPlotScatter(int *kapaID, float *dXfit, float *dYfit, float *mag, int npts) {
 
     Graphdata graphdata;
     KapaSection section;
-    section.x = 0; section.y = 0; section.dx = 1; section.dy = 1;
-    section.name = "junk";
+
+    if (!initWindow(kapaID)) return 0;
 
     KapaInitGraph(&graphdata);
-    KapaClearPlots(kapa2);
-    KapaSetSection(kapa2, &section);
-
-    graphdata.ptype = 0;
-    graphdata.style = 0;
-    graphdata.xmin = 0;
-    graphdata.xmax = npts;
-    graphdata.ymin = 0;
-    graphdata.ymax = data[npts-1];
-
-    KapaSetFont(kapa2, "helvetica", 14);
-    KapaSetLimits(kapa2, &graphdata);
-    KapaBox(kapa2, &graphdata);
-    KapaSendLabel( kapa2, "Object", KAPA_LABEL_XM);
-    KapaSendLabel( kapa2, "Offset(pixels)", KAPA_LABEL_YM);
-    KapaSendLabel( kapa2, "Astrometric Offset with fit cutoff",
-                   KAPA_LABEL_XP);
-    KapaPrepPlot(kapa2, npts, &graphdata);
-    KapaPlotVector(kapa2, npts, x, "x");
-    KapaPlotVector(kapa2, npts, data, "y");
-
-    graphdata.color = KapaColorByName("red");
-    KapaPrepPlot(kapa2, 2, &graphdata);
-    yline[0] = (float) thresh;
-    yline[1] = (float) thresh;
-    xline[0] = graphdata.xmin;
-    xline[1] = graphdata.xmax;
-    KapaPlotVector(kapa2, 2, xline, "x");
-    KapaPlotVector(kapa2, 2, yline, "y");
-
-    askUser(&plotScatter);
-    return 1;
-}
-
-
+    KapaClearSections(*kapaID);
+    KapaSetFont(*kapaID, "helvetica", 14);
+
+    section.name = "a";
+    section.x = 0.0; section.y = 0.0; section.dx = 1.0; section.dy = 0.5;
+    section.bg = KapaColorByName("white");
+    KapaSetSection(*kapaID, &section);
+
+    graphdata.ptype = 2;
+    graphdata.style = 2;
+
+    if(!scaleGraphdata(mag, dXfit, &graphdata, npts)) return 0;
+    KapaSetLimits(*kapaID, &graphdata);
+    KapaBox(*kapaID, &graphdata);
+    KapaSendLabel(*kapaID, "mag", KAPA_LABEL_XM);
+    KapaSendLabel(*kapaID, "dXfit", KAPA_LABEL_YM);
+    KapaPrepPlot(*kapaID, npts, &graphdata);
+    KapaPlotVector(*kapaID, npts, mag, "x");
+    KapaPlotVector(*kapaID, npts, dXfit, "y");
+
+    section.name = "b";
+    section.x = 0.0; section.y = 0.5; section.dx = 1.0; section.dy = 0.5;
+    section.bg = KapaColorByName("white");
+    KapaSetSection(*kapaID, &section);
+
+    graphdata.ptype = 2;
+    graphdata.style = 2;
+
+    if(!scaleGraphdata(mag, dYfit, &graphdata, npts)) return 0;
+    KapaSetLimits(*kapaID, &graphdata);
+    KapaBox(*kapaID, &graphdata);
+    KapaSendLabel(*kapaID, "mag", KAPA_LABEL_XM);
+    KapaSendLabel(*kapaID, "dYfit", KAPA_LABEL_YM);
+    KapaPrepPlot(*kapaID, npts, &graphdata);
+    KapaPlotVector(*kapaID, npts, mag, "x");
+    KapaPlotVector(*kapaID, npts, dYfit, "y");
+
+    return TRUE;
+}
 
 /** plot raw vs ref (L, M). Only those whose distance is < drMax are used in fit*/
-int relastroVisualPlotRawRef(StarData *raw, StarData *ref, double dRmax, int numObj) {
-
-  if( !isVisual || !plotRawRef) return 1;
-  if( !initWindow(&kapa)) return 0;
+int relastroVisualPlotFittedStars(int *kapaID, float *rawX, float *rawY, float *refX, float *refY, int numNoFit, float *rawXfit, float *rawYfit, float *refXfit, float *refYfit, int numFit) {
+
+  Graphdata graphdata;
+
+  if (!initWindow(kapaID)) return 0;
+
+  KapaInitGraph(&graphdata);
+  KapaClearPlots(*kapaID);
+
+  graphdata.ptype = 7;
+  graphdata.style = 2;
+
+  if(!scaleGraphdata(rawXfit, rawYfit, &graphdata, numFit)) {
+      fprintf(stderr, "Not enough finite points for plotting");
+      return 0;
+  }
+
+  KapaSetFont(*kapaID, "helvetica", 14);
+  KapaSetLimits(*kapaID, &graphdata);
+  KapaBox(*kapaID, &graphdata);
+  KapaSendLabel( *kapaID, "X", KAPA_LABEL_XM);
+  KapaSendLabel( *kapaID, "Y", KAPA_LABEL_YM);
+  KapaSendLabel( *kapaID, "orange, red, green, blue: (raw, ref), (nofit, fit)",
+                 KAPA_LABEL_XP);
+
+  graphdata.color = KapaColorByName("orange");
+  graphdata.size = 1;
+  KapaPrepPlot(*kapaID, numNoFit, &graphdata);
+  KapaPlotVector(*kapaID, numNoFit, rawX, "x");
+  KapaPlotVector(*kapaID, numNoFit, rawY, "y");
+
+  graphdata.color = KapaColorByName("red");
+  graphdata.size = 2;
+  KapaPrepPlot(*kapaID, numNoFit, &graphdata);
+  KapaPlotVector(*kapaID, numNoFit, refX, "x");
+  KapaPlotVector(*kapaID, numNoFit, refY, "y");
+
+  graphdata.color = KapaColorByName("green");
+  graphdata.size = 1;
+  KapaPrepPlot(*kapaID, numFit,  &graphdata);
+  KapaPlotVector(*kapaID, numFit, rawXfit, "x");
+  KapaPlotVector(*kapaID, numFit, rawYfit, "y");
+
+  graphdata.color = KapaColorByName("blue");
+  graphdata.size = 2;
+  KapaPrepPlot(*kapaID, numFit, &graphdata);
+  KapaPlotVector(*kapaID, numFit, refXfit, "x");
+  KapaPlotVector(*kapaID, numFit, refYfit, "y");
+
+  return TRUE;
+}
+
+/** create various plots using the data in raw and ref **/
+int relastroVisualPlotChipFit(StarData *raw, StarData *ref, double dRmax, int numObj) {
 
   float *rawX, *rawY,  *refX, *refY;
   float *rawXfit, *rawYfit, *refXfit, *refYfit;
   float *magRaw, *magRef, *magRawfit, *magReffit;
-  float *xVec, *yVec;
+  float *dXfit, *dYfit, *dPos;
   int numFit = 0, numNoFit = 0;
   double dL, dM, dR;
+
+  if (!isVisual) return TRUE;
 
   ALLOCATE(rawX,      float, numObj);
@@ -273 +335 @@
   ALLOCATE(refX,      float, numObj);
   ALLOCATE(refY,      float, numObj);
+  ALLOCATE(magRaw,    float, numObj);
+  ALLOCATE(magRef,    float, numObj);
+
   ALLOCATE(rawXfit,   float, numObj);
   ALLOCATE(rawYfit,   float, numObj);
   ALLOCATE(refXfit,   float, numObj);
   ALLOCATE(refYfit,   float, numObj);
-  ALLOCATE(magRaw,    float, numObj);
-  ALLOCATE(magRef,    float, numObj);
   ALLOCATE(magRawfit, float, numObj);
   ALLOCATE(magReffit, float, numObj);
+  ALLOCATE(dXfit,     float, numObj);
+  ALLOCATE(dYfit,     float, numObj);
+  ALLOCATE(dPos,      float, numObj);
 
   int i;
   for(i = 0; i < numObj; i++) {
-    if  (raw[i].mask) continue;
+    if (raw[i].mask) continue; // XXX 
 
     dL = raw[i].L - ref[i].L;
     dM = raw[i].M - ref[i].M;
     dR = hypot (dL, dM);
+
+    // XXX change the selection to a mask-based thing
     if (dR > dRmax) {
+      // UNFITTED values
       rawX[numNoFit] = raw[i].X;
       rawY[numNoFit] = raw[i].Y;
@@ -298 +367 @@
       numNoFit++;
     } else {
-      rawXfit[numFit] = raw[i].X;
-      rawYfit[numFit] = raw[i].Y;
-      refXfit[numFit] = ref[i].X;
-      refYfit[numFit] = ref[i].Y;
-      magRaw[numFit] = raw[i].Mag;
-      magRef[numFit] = ref[i].Mag;
+      // FITTED values
+      rawXfit[numFit]   = raw[i].X;
+      rawYfit[numFit]   = raw[i].Y;
+      refXfit[numFit]   = ref[i].X;
+      refYfit[numFit]   = ref[i].Y;
+      magRawfit[numFit] = raw[i].Mag;
+      magReffit[numFit] = ref[i].Mag;
+      dPos[numFit]      = ref[i].dPos;
+      dXfit[numFit]     = raw[i].X - ref[i].X;
+      dYfit[numFit]     = raw[i].Y - ref[i].Y;
       numFit++;
     }
@@ -310 +383 @@
   if (numFit == 0) return 0;
 
-  Graphdata graphdata;
-
-  KapaInitGraph(&graphdata);
-  KapaClearPlots(kapa);
-
-  graphdata.ptype = 7;
-  graphdata.style = 2;
-
-  if(!scaleGraphdata(rawXfit, rawYfit, &graphdata, numFit)) {
-      fprintf(stderr, "Not enough finite points for plotting");
-      return 0;
-  }
-
-  KapaSetFont(kapa, "helvetica", 14);
-  KapaSetLimits(kapa, &graphdata);
-  KapaBox(kapa, &graphdata);
-  KapaSendLabel( kapa, "X", KAPA_LABEL_XM);
-  KapaSendLabel( kapa, "Y", KAPA_LABEL_YM);
-  KapaSendLabel( kapa, "orange, red, green, blue: (raw, ref), (nofit, fit)",
-                 KAPA_LABEL_XP);
-
-  graphdata.color = KapaColorByName("orange");
-  graphdata.size = 1;
-  KapaPrepPlot(kapa, numNoFit, &graphdata);
-  KapaPlotVector(kapa, numNoFit, rawX, "x");
-  KapaPlotVector(kapa, numNoFit, rawY, "y");
-
-  graphdata.color = KapaColorByName("red");
-  graphdata.size = 2;
-  KapaPrepPlot(kapa, numNoFit, &graphdata);
-  KapaPlotVector(kapa, numNoFit, refX, "x");
-  KapaPlotVector(kapa, numNoFit, refY, "y");
-
-  graphdata.color = KapaColorByName("green");
-  graphdata.size = 1;
-  KapaPrepPlot(kapa, numFit,  &graphdata);
-  KapaPlotVector(kapa, numFit, rawXfit, "x");
-  KapaPlotVector(kapa, numFit, rawYfit, "y");
-
-  graphdata.color = KapaColorByName("blue");
-  graphdata.size = 2;
-  KapaPrepPlot(kapa, numFit, &graphdata);
-  KapaPlotVector(kapa, numFit, refXfit, "x");
-  KapaPlotVector(kapa, numFit, refYfit, "y");
-
-  ALLOCATE(xVec, float, numFit);
-  ALLOCATE(yVec, float, numFit);
-
-  //plot the fitted objects as vectors
-  for(i = 0; i < numFit; i++) {
-      xVec[i] = rawXfit[i] - refXfit[i];
-      yVec[i] = rawYfit[i] - refYfit[i];
-  }
-
-  plotVectorField(rawXfit, rawYfit, xVec, yVec, numFit, &kapa3, dRmax);
-  if(!residPlot(rawXfit, rawYfit, xVec, yVec, numFit, &kapa2)) {
-      fprintf(stderr, "Unable to plot residuals");
-      return 0;
-  }
-
-  FREE(xVec);
-  FREE(yVec);
+  // 4-panel plot of x vs dx, y vs dx, x vs dy, y vs dy 
+  relastroVisualRawRef(&kapa1, rawXfit, rawYfit, dXfit, dYfit, dPos, numFit);
+
+  // vector line plot for the fit
+  relastroVisualVectorField(&kapa2, rawXfit, rawYfit, dXfit, dYfit, numFit, dRmax);
+
+  // dXfit, dYfit vs mag
+  relastroVisualPlotScatter(&kapa3, dXfit, dYfit, magRawfit, numFit);
+
+  // plot the positions of the fitted stars on the chip
+  relastroVisualPlotFittedStars(&kapa4, rawX, rawY, refX, refY, numNoFit, rawXfit, rawYfit, refXfit, refYfit, numFit);
 
   askUser(&plotRawRef);
 
+  FREE(dXfit);
+  FREE(dYfit);
+  FREE(dPos);
   FREE(rawX);
   FREE(rawY);
@@ -388 +413 @@
   FREE(magReffit);
 
-  return 1;
+  return TRUE;
 }
 
@@ -403 +428 @@
   KapaSection section;
   
-  if (!isVisual || !plotOutliers) return 1;
-  if (!initWindow(&kapa)) return 0;
+  if (!isVisual || !plotOutliers) return TRUE;
+  if (!initWindow(&kapa1)) return 0;
   
   // populate vectors
@@ -426 +451 @@
   for(i = 0; i < Nmeasure; i++, m++) {
     meas = catalog[0].measure[m];
-    if (!MeasFilterTest(&meas)) continue;
+    if (!MeasFilterTest(&meas, FALSE)) continue;
     xmin = MIN(xmin, meas.dR);
     xmax = MAX(xmax, meas.dR);
@@ -459 +484 @@
   section.x = 0; section.y = 0; section.dx = 1; section.dy = 1;
   section.name = "junk";
+  section.bg = KapaColorByName("white");
   
   KapaInitGraph(&graphdata);
-  KapaClearPlots(kapa);
-  KapaSetFont(kapa, "helvetica", 14);
+  KapaClearPlots(kapa1);
+  KapaSetFont(kapa1, "helvetica", 14);
  
   graphdata.ptype = 7;
@@ -472 +498 @@
   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)",
+  KapaSetSection(kapa1, &section);
+  KapaSetLimits(kapa1, &graphdata);
+  KapaBox(kapa1, &graphdata);
+
+  KapaSendLabel( kapa1, "RA (arcsec)", KAPA_LABEL_XM);
+  KapaSendLabel( kapa1, "Dec (arcsec)", KAPA_LABEL_YM);
+  KapaSendLabel( kapa1, "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");
+  KapaPrepPlot(kapa1, Nin, &graphdata);
+  KapaPlotVector(kapa1, Nin, Rin, "x");
+  KapaPlotVector(kapa1, Nin, Din, "y");
 
   graphdata.color = KapaColorByName("red");
-  KapaPrepPlot(kapa, Nout, &graphdata);
-  KapaPlotVector(kapa, Nout, Rout, "x");
-  KapaPlotVector(kapa, Nout, Dout, "y");
+  KapaPrepPlot(kapa1, Nout, &graphdata);
+  KapaPlotVector(kapa1, Nout, Rout, "x");
+  KapaPlotVector(kapa1, 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");
-
-
-  
+  KapaPrepPlot(kapa1, 100, &graphdata);
+  KapaPlotVector(kapa1, 100, xCirc, "x");
+  KapaPlotVector(kapa1, 100, yCirc, "y");
+
   askUser(&plotOutliers);
 
@@ -510 +534 @@
 }
   
-  
+# if (0)
+int relastroVisualSummaryChips() {
+
+  // plot the dXsys, dYsys histograms
+
+  // plot x vs dx, y vs dy, etc for all mosaics
+
+  // plot a map of median star scatter
+
+}
+# endif
+

trunk/Ohana/src/relastro/src/select_images.c

@@ -15 +15 @@
 void dsortindex (double *X, off_t *Y, int N);
 off_t getRegionStartByRA (double R, double *Rref, off_t Nregions);
+
+# define MARKTIME(MSG,...) { \
+  float dtime; \
+  gettimeofday (&stop, (void *) NULL); \
+  dtime = DTIME (stop, start); \
+  fprintf (stderr, MSG, __VA_ARGS__); }
 
 Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage) {
@@ -25 +31 @@
   Coords tcoords;
   SkyRegionCoords *skycoords;
-  
+  struct timeval start, stop;
+  
+  double RmaxSkyRegion, RminSkyRegion, RmidSkyRegion, DminSkyRegion, DmaxSkyRegion;
+
   double *RmaxSky;
   off_t *index;
@@ -40 +49 @@
     return NULL;
   }
+
+  gettimeofday (&start, (void *) NULL);
 
   // the comparison is made in the catalog local projection. below we set crval1,2
@@ -52 +63 @@
   ALLOCATE (RmaxSky, double, skylist[0].Nregions);
   ALLOCATE (index, off_t, skylist[0].Nregions);
+
+  RminSkyRegion = +360.0;
+  RmaxSkyRegion = -360.0;
+  DminSkyRegion = +90.0;
+  DmaxSkyRegion = -90.0;
 
   /* compare with each region file */
@@ -80 +96 @@
     skycoords[i].Xc[3] -= dx; skycoords[i].Yc[3] += dy;
     skycoords[i].Xc[4] -= dx; skycoords[i].Yc[4] -= dy;
-  }
+
+    RminSkyRegion = MIN(RminSkyRegion, skylist[0].regions[i][0].Rmin);
+    RmaxSkyRegion = MAX(RmaxSkyRegion, skylist[0].regions[i][0].Rmax);
+    DminSkyRegion = MIN(DminSkyRegion, skylist[0].regions[i][0].Dmin);
+    DmaxSkyRegion = MAX(DmaxSkyRegion, skylist[0].regions[i][0].Dmax);
+  }
+  RmidSkyRegion = 0.5*(RminSkyRegion + RmaxSkyRegion);
+  MARKTIME("create sky region coords: %f sec\n", dtime);
 
   dsortindex (RmaxSky, index, skylist[0].Nregions);
+  MARKTIME("sort sky coords: %f sec\n", dtime);
 
   if (VERBOSE) fprintf (stderr, "finding images\n");
   BuildChipMatch (timage, Ntimage);
+  MARKTIME("build chip match: %f sec\n", dtime);
 
   nimage = 0;
@@ -124 +149 @@
     }
 
-    /* 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 */
     double RminImage = 360.0;
+    double RmaxImage =   0.0;
+    double DminImage = +90.0;
+    double DmaxImage = -90.0;
+    // int leftside = FALSE;
     for (j = 0; j < 5; j++) {
       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]);
+      DminImage = MIN(DminImage, Di[j]);
+      DmaxImage = MAX(DmaxImage, Di[j]);
+    }
+    if (RmaxImage - RminImage > 180.0) {
+        double tmp = RminImage;
+        RmaxImage = RminImage;
+        RminImage = tmp - 360.0;
+    }
+
+    // check that this image is even in range of the searched region
+    if (DminImage > DmaxSkyRegion) continue;
+    if (DmaxImage < DminSkyRegion) continue;
+    
+    // the sky region RA is defined to be 0 - 360.0
+    if (RminImage > RmaxSkyRegion) continue;
+    if (RmaxImage < RminSkyRegion) continue;
+
+    // image overlaps region, keep it
+    if (USE_BASIC_CHECK) goto found_it;
 
     // RA(nStart) is guaranteed to be < RminImage:
@@ -161 +219 @@
       for (j = 0; (j < 4) && !found; j++) {
         found = corner_check (&Xi[j], &Yi[j], &skycoords[m].Xc[0], &skycoords[m].Yc[0]);
+        if (found) goto found_it;
       }
       /* check if catalog corner inside image */
       for (j = 0; (j < 4) && !found; j++) {
         found = corner_check (&skycoords[m].Xc[j], &skycoords[m].Yc[j], &Xi[0], &Yi[0]);
+        if (found) goto found_it;
       }
       /* check if edges cross */
@@ -170 +230 @@
         for (k = 0; (k < 4) && !found; k++) {
           found = edge_check (&Xi[j], &Yi[j], &skycoords[m].Xc[k], &skycoords[m].Yc[k]);
+          if (found) goto found_it;
         }
       }
-      if (!found) continue;
-
-      image[nimage] = timage[i]; 
-      /* always allow 'few' images to succeed, if possible */
-      if (image[nimage].flags & ID_IMAGE_ASTROM_FEW) { 
-        image[nimage].flags &= ~ID_IMAGE_ASTROM_FEW;
-      }
-      if (RESET) {
-        // XXX do we need / want to do this in relastro?
-        assignMcal (&image[nimage], (double *) NULL, -1);
-        image[nimage].dMcal = NAN;
-        image[nimage].flags &= ~badImage;
-      }
-      line_number[nimage] = i;
-      nimage ++;
-      if (nimage == NIMAGE) {
-        NIMAGE += 100;
-        REALLOCATE (image, Image, NIMAGE);
-        REALLOCATE (line_number, off_t, NIMAGE);
-      }
-    }
-  }
-      
+    }
+    if (!found) continue;
+
+  found_it:
+    image[nimage] = timage[i]; 
+    /* always allow 'few' images to succeed, if possible */
+    if (image[nimage].flags & ID_IMAGE_ASTROM_FEW) { 
+      image[nimage].flags &= ~ID_IMAGE_ASTROM_FEW;
+    }
+    if (RESET) {
+      // XXX do we need / want to do this in relastro?
+      // assignMcal (&image[nimage], (double *) NULL, -1);
+      // image[nimage].Mcal = NAN;
+      // image[nimage].dMcal = NAN;
+      image[nimage].flags &= ~badImage;
+    }
+    line_number[nimage] = i;
+    nimage ++;
+    if (nimage == NIMAGE) {
+      NIMAGE += 100;
+      REALLOCATE (image, Image, NIMAGE);
+      REALLOCATE (line_number, off_t, NIMAGE);
+    }
+  }
+  MARKTIME("finish image selection: %f sec\n", dtime);
+    
   if (VERBOSE) fprintf (stderr, "found "OFF_T_FMT" images\n",  nimage);