Changeset 33568

Timestamp:

Mar 19, 2012, 6:09:12 PM (14 years ago)

Author:

eugene

Message:

count images with too few, nearly too few objects; require images to be completely within the selected region (or we have gradients across the images); use psastro style order selection based on number of objects; free some possible leaks; clean up bcatalog messages (and list counts for selected objects at each stage); pass SIGMA_LIM to relastro_client -load-objects

Location:

branches/eam_branches/ipp-20111122/Ohana/src/relastro/src

Files:

• FitChip.c (modified) (7 diffs)
• ImageOps.c (modified) (1 diff)
• UpdateChips.c (modified) (2 diffs)
• bcatalog.c (modified) (10 diffs)
• fitpoly.c (modified) (1 diff)
• load_catalogs.c (modified) (3 diffs)
• select_images.c (modified) (1 diff)

branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/FitChip.c

@@ -4 +4 @@
 int FitChip (StarData *raw, StarData *ref, int Nmatch, Image *image) {
 
-  int i, NstatFull, Nstat, Niter, skip;
+  int i, NstatFull, Nstat, Niter;
   float dLsig, dMsig, dRsig;
   float dLsigFull, dMsigFull, dRsigFull;
@@ -52 +52 @@
     // fit the requested order polynomial
     if (CHIPORDER > 0) {
-      image[0].coords.Npolyterms = CHIPORDER;
-    }
+      int Ndof_min = 3;
+      int order_max = 0.5*(sqrt(4*Nmatch - 4*Ndof_min + 1) - 3);
+      int order_use = MIN (CHIPORDER, order_max);
+      if (order_use < 1) {
+        if (VERBOSE2) fprintf (stderr, "insufficient measurements (%d) for linear fit\n", Nmatch);
+        image[0].flags |= ID_IMAGE_ASTROM_FEW;
+        return FALSE;
+      }
+      image[0].coords.Npolyterms = order_use;
+    }
+
     if (fit) fit_free (fit);
     fit = fit_init (image[0].coords.Npolyterms);
@@ -63 +72 @@
     }
 
+# if (0)
     // check if the fit has enough data points for the polynomial order
     skip = FALSE;
@@ -86 +96 @@
       return FALSE;
     }
+# endif
 
     // measure the fit, update the coords & object coordinates
@@ -91 +102 @@
       fprintf (stderr, "failed to fit new model\n");
       image[0].flags |= ID_IMAGE_ASTROM_FAIL;
+      if (fit) fit_free (fit);
       return FALSE;
     }
@@ -97 +109 @@
       fprintf (stderr, "failed to fit new model\n");
       image[0].flags |= ID_IMAGE_ASTROM_FAIL;
+      if (fit) fit_free (fit);
       return FALSE;
     }
@@ -133 +146 @@
   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);
+  if (VERBOSE2) 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;

branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/ImageOps.c

@@ -181 +181 @@
   }
 
-  for (i = 0; VERBOSE && (i < Nimage); i++) {
-    name = GetPhotcodeNamebyCode (image[i].photcode);
-    fprintf (stderr, "image "OFF_T_FMT" (%d, %s) has "OFF_T_FMT" of %d measures (%s, %s)\n",  i,  image[i].imageID, image[i].name, N_onImage[i], image[i].nstar,
-             ohana_sec_to_date(image[i].tzero), name);
-  }
+  // watch for chips with few stars
+  int Nfew = 0;
+  int Nbad = 0;
+  for (i = 0; i < Nimage; i++) {
+    // ignore the PHU entries
+    if (!strcmp(&image[i].coords.ctype[4], "-DIS")) continue;
+    if (VERBOSE2 || (MATCHCAT && (N_onImage[i] < 20))) {
+      name = GetPhotcodeNamebyCode (image[i].photcode);
+      fprintf (stderr, "image "OFF_T_FMT" (%d, %s) has "OFF_T_FMT" of %d measures (%s, %s)",  i,  image[i].imageID, image[i].name, N_onImage[i], image[i].nstar,
+               ohana_sec_to_date(image[i].tzero), name);
+      if (N_onImage[i] < 20) {
+        fprintf (stderr, "*");
+        Nfew ++;
+      } 
+      if (N_onImage[i] < 15) {
+        fprintf (stderr, "*** warning ***");
+        Nbad ++;
+      } 
+      fprintf (stderr, "\n");
+    }
+  }
+  fprintf (stderr, OFF_T_FMT" total images, %d with < 20 measurements, %d with < 15\n", Nimage, Nfew, Nbad);
 }
 

branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/UpdateChips.c

@@ -25 +25 @@
   ALLOCATE (mode, char, Nimage);
 
+  // XXX for faster processing in the future, this can be easily run in parallel
+  // each chip is fitted independently, so we could do N at once in parallel
   for (i = 0; i < Nimage; i++) {
 
@@ -47 +49 @@
       Nskip ++;
       mode[i] = 0;
+      free (raw);
       continue;
     }

branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/bcatalog.c

@@ -1 +1 @@
 # include "relastro.h"
+
+static int Nskip1 = 0;
+static int Nskip2 = 0;
+static unsigned int Tref = 1323353985;
+static short Cref = 10355;
 
 int bcatalog (Catalog *subcatalog, Catalog *catalog) {
@@ -34 +39 @@
     if (catalog[0].average[i].Nmeasure <= SRC_MEAS_TOOFEW) {
       continue;
+      // XXX can't (easily) count this for a given chip
     }
     
@@ -55 +61 @@
       if (!MeasFilterTest(&catalog[0].measure[offset], TRUE)) {
         catalog[0].measure[offset].dbFlags &= ~ID_MEAS_USED_CHIP;
+        if (FALSE && (abs(catalog[0].measure[offset].t - Tref) < 10) && (catalog[0].measure[offset].photcode == Cref)) {
+          Nskip1 ++;
+        }
         continue;
       }
@@ -61 +70 @@
       if (FlagOutlier && (catalog[0].measure[offset].dbFlags & ID_MEAS_POOR_ASTROM)) {
         catalog[0].measure[offset].dbFlags &= ~ID_MEAS_USED_CHIP;
+        if (FALSE && (abs(catalog[0].measure[offset].t - Tref) < 10) && (catalog[0].measure[offset].photcode == Cref)) {
+          Nskip2 ++;
+        }
         continue;
       }
@@ -115 +127 @@
   if (MaxDensityUse) {
     LimitDensityCatalog_ByNmeasure (subcatalog, catalog);
-  }
-
-  if (VERBOSE) {
-    fprintf (stderr, OFF_T_FMT": using "OFF_T_FMT" stars ("OFF_T_FMT" measures) for catalog\n",  i,  subcatalog[0].Naverage,  subcatalog[0].Nmeasure);
-   }
+  } else {
+    if (VERBOSE) {
+      char *basename = filebasename (catalog[0].filename);
+      fprintf (stderr, "subset of "OFF_T_FMT" ("OFF_T_FMT" total) stars, "OFF_T_FMT" ("OFF_T_FMT" total) measures for catalog %s\n", 
+               subcatalog[0].Naverage, catalog[0].Naverage, subcatalog[0].Nmeasure, catalog[0].Nmeasure, basename);
+      free (basename);
+    }
+  }
   return (TRUE);
+}
+
+void bcatalog_show_skips () {
+  fprintf (stderr, "Nskip: %d, %d\n", Nskip1, Nskip2);
 }
 
@@ -140 +159 @@
 
 /* this version does NOT use AverageTiny, MeasureTiny */ 
-int LimitDensityCatalog_ByNmeasure (Catalog *subcatalog, Catalog *catalog) {
+int LimitDensityCatalog_ByNmeasure (Catalog *subcatalog, Catalog *oldcatalog) {
 
   Catalog tmpcatalog;
@@ -148 +167 @@
   int Nsecfilt = GetPhotcodeNsecfilt ();
 
-  gfits_scan (&catalog[0].header, "RA0",  "%lf", 1, &Rmin);
-  gfits_scan (&catalog[0].header, "DEC0", "%lf", 1, &Dmin);
-  gfits_scan (&catalog[0].header, "RA1",  "%lf", 1, &Rmax);
-  gfits_scan (&catalog[0].header, "DEC1", "%lf", 1, &Dmax);
+  gfits_scan (&oldcatalog[0].header, "RA0",  "%lf", 1, &Rmin);
+  gfits_scan (&oldcatalog[0].header, "DEC0", "%lf", 1, &Dmin);
+  gfits_scan (&oldcatalog[0].header, "RA1",  "%lf", 1, &Rmax);
+  gfits_scan (&oldcatalog[0].header, "DEC1", "%lf", 1, &Dmax);
 
   if (VERBOSE2) fprintf (stderr, "extracting from catalog covering region %f,%f to %f,%f\n", Rmin, Dmin, Rmax, Dmax);
@@ -177 +196 @@
   for (i = 0; i < Naverage; i++) {
     index[i] = i;
-    value[i] = catalog[0].average[i].Nmeasure;
+    value[i] = subcatalog[0].average[i].Nmeasure;
   }
   sort_by_Nmeasure (value, index, Naverage);
@@ -212 +231 @@
 
   if (VERBOSE) {
-    fprintf (stderr, "limited to "OFF_T_FMT" of "OFF_T_FMT" stars ("OFF_T_FMT" of "OFF_T_FMT" measures) for catalog %s\n", 
-             Nmax, subcatalog[0].Naverage, Nmeasure, subcatalog[0].Nmeasure,  catalog[0].filename);
+    char *basename = filebasename (oldcatalog[0].filename);
+    fprintf (stderr, "limited to "OFF_T_FMT" ("OFF_T_FMT" subset, "OFF_T_FMT" total) stars, "OFF_T_FMT" ("OFF_T_FMT" subset, "OFF_T_FMT" total) measures for catalog %s\n", 
+             Nmax, subcatalog[0].Naverage, oldcatalog[0].Naverage, Nmeasure, subcatalog[0].Nmeasure,  oldcatalog[0].Nmeasure, basename);
+    free (basename);
   }
 
@@ -227 +248 @@
   subcatalog[0].Naverage = Nmax;
   subcatalog[0].Nmeasure = Nmeasure;
-  subcatalog[0].Nsecfilt = catalog[0].Nsecfilt;
-  subcatalog[0].Nsecf_mem = Naverage * catalog[0].Nsecfilt;
+  subcatalog[0].Nsecfilt = oldcatalog[0].Nsecfilt;
+  subcatalog[0].Nsecf_mem = Naverage * oldcatalog[0].Nsecfilt;
 
   return (TRUE);

branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/fitpoly.c

@@ -149 +149 @@
   
   if (!dgaussjordan (matrix, vector, fit[0].Nelems, 2)) {
+    array_free (matrix, fit[0].Nelems);
+    array_free (vector, fit[0].Nelems);
     return (FALSE);
   }

branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/load_catalogs.c

@@ -2 +2 @@
 
 Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog);
+void bcatalog_show_skips ();
 
 Catalog *load_catalogs (SkyList *skylist, int *Ncatalog, int subselect, int hostID, char *hostpath) {
@@ -83 +84 @@
   }
 
+  // XXX TEST : bcatalog_show_skips();
+
   Nstar = 0;
   for (i = 0; i < Ncat; i++) {
@@ -139 +142 @@
 
     char command[1024];
-    snprintf (command, 1024, "relastro_client -load-objects %s -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -statmode %s -minerror %f", 
-              table->hosts[i].results, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, STATMODE, MIN_ERROR);
+    snprintf (command, 1024, "relastro_client -load-objects %s -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -statmode %s -minerror %f -D RELASTRO_SIGMA_LIM %f", 
+              table->hosts[i].results, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, STATMODE, MIN_ERROR, SIGMA_LIM);
 
     char tmpline[1024];

branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/select_images.c

@@ -196 +196 @@
     if (RmaxImage < RminSkyRegion) continue;
 
+    if (!strncmp(timage[i].name, "o5903g0638o", 10) && (timage[i].photcode == 10355)) {
+      fprintf (stderr, "test image\n");
+    }
+
+    // require that the full image be inside the region of interest (only for non PHU
+    // images) XXX : if we calibrate the mosaic, require that the full mosaic be inside
+    // the region and all of its chips..
+# define FULL_OVERLAP 1
+    if (FULL_OVERLAP && strcmp(&timage[i].coords.ctype[4], "-DIS")) {
+      if (RmaxImage > UserPatch.Rmax) continue;
+      if (RminImage < UserPatch.Rmin) continue;
+      if (DmaxImage > UserPatch.Dmax) continue;
+      if (DminImage < UserPatch.Dmin) continue;
+    }
+
     // image overlaps region, keep it
     if (USE_BASIC_CHECK) goto found_it;