Changeset 36769

Timestamp:

May 23, 2014, 9:42:18 AM (12 years ago)

Author:

eugene

Message:

add option to reset image astrometry to a specific photcode (eg, 2MASS_J); add DensityLimit by small grid cell (require N per 200x200 arcsec^{2 box); require more stars for higher order fits}

Location:

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

Files:

• FitChip.c (modified) (2 diffs)
• bcatalog.c (modified) (8 diffs)
• launch_region_hosts.c (modified) (2 diffs)
• relastro_images.c (modified) (1 diff)
• relastro_parallel_images.c (modified) (1 diff)

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

@@ -50 +50 @@
     relastroVisualPlotChipFit(raw, ref, dRmax, Nmatch);
 
-    // fit the requested order polynomial
+    // set the maximum order for the polynomial (based on number of stars)
+    int order_use = 0;
+    if (Nmatch >  5) order_use = 1; // 4 stars per polynomial term (per dimension)
+    if (Nmatch > 30) order_use = 2; // 5 stars per polynomial term (per dimension)
+    if (Nmatch > 60) order_use = 3; // 6 stars per polynomial term (per dimension)
+    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 (0)
     if (CHIPORDER > 0) {
       int Ndof_min = 3;
@@ -62 +74 @@
       image[0].coords.Npolyterms = order_use;
     }
+# endif
 
     if (fit) fit_free (fit);

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

@@ -1 +1 @@
 # include "relastro.h"
 
+static int Nkeep1 = 0;
+static int Nkeep2 = 0;
 static int Nskip1 = 0;
 static int Nskip2 = 0;
-static unsigned int Tref = 1323353985;
-static short Cref = 10355;
+
+// test image: 2013/06/15,13:25:51, GPC1.r.XY50
+static int CHECK_TEST_IMAGE = FALSE;
+static unsigned int Tref = 1378812312; 
+static short Cref = 10001;
+
+int LimitDensityCatalog_ByNmeasureGrid (Catalog *subcatalog, Catalog *oldcatalog);
 
 int bcatalog (Catalog *subcatalog, Catalog *catalog) {
@@ -60 +67 @@
       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)) {
+        if (CHECK_TEST_IMAGE && (abs(catalog[0].measure[offset].t - Tref) < 10) && (catalog[0].measure[offset].photcode == Cref)) {
           Nskip1 ++;
         }
         continue;
+      }
+      if (CHECK_TEST_IMAGE && (abs(catalog[0].measure[offset].t - Tref) < 10) && (catalog[0].measure[offset].photcode == Cref)) {
+        Nkeep1 ++;
       }
 
@@ -87 +97 @@
       // but before the final average properties are calculated, these measurements may be
       // allowed.
+
+      if (CHECK_TEST_IMAGE && (abs(catalog[0].measure[offset].t - Tref) < 10) && (catalog[0].measure[offset].photcode == Cref)) {
+        Nkeep2 ++;
+      }
 
       CopyMeasureToTiny (&subcatalog[0].measureT[Nmeasure], &catalog[0].measure[offset]);
@@ -125 +139 @@
   // limit the total number of stars in the catalog
   if (MaxDensityUse) {
-    LimitDensityCatalog_ByNmeasure (subcatalog, catalog);
+    LimitDensityCatalog_ByNmeasureGrid (subcatalog, catalog);
   } else {
     if (VERBOSE2) {
@@ -134 +148 @@
     }
   }
+  if (CHECK_TEST_IMAGE && (Nkeep1 + Nkeep2 + Nskip1 + Nskip2 > 0)) {
+    fprintf (stderr, "kept %d %d, skipped %d %d\n", Nkeep1, Nkeep2, Nskip1, Nskip2);
+  }
+
   return (TRUE);
 }
@@ -139 +157 @@
 void bcatalog_show_skips () {
   fprintf (stderr, "Nskip: %d, %d\n", Nskip1, Nskip2);
+  fprintf (stderr, "Nkeep: %d, %d\n", Nkeep1, Nkeep2);
 }
 
@@ -155 +174 @@
 # undef COMPARE
 
+}
+
+/* generate a grid in a locally-projected coordinate system, try to select average entries 
+   from each grid cell in decending Nmeasure order.
+ */ 
+int LimitDensityCatalog_ByNmeasureGrid (Catalog *subcatalog, Catalog *oldcatalog) {
+
+  off_t i, j;
+  int ix, iy;
+
+  Catalog tmpcatalog;
+
+  double Rmin, Rmax, Dmin, Dmax;
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  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);
+
+  float AREA = fabs(Dmax - Dmin) * fabs(Rmax - Rmin) * cos (0.5*RAD_DEG*(Dmax + Dmin));
+  assert (AREA > 0);
+
+  off_t Nmax = MaxDensityValue * AREA;
+  if (subcatalog[0].Naverage <= Nmax) {
+    if (VERBOSE) {
+      fprintf (stderr, "subcatalog has less than the max density\n");
+    }
+    return (TRUE);
+  }
+
+  off_t Naverage = subcatalog[0].Naverage;
+
+  // generate a grid in locally projected space
+  double Rc = 0.5*(Rmin + Rmax);
+  double Dc = 0.5*(Dmin + Dmax);
+
+  /* project coordinates to a plane centered on the object with units of arcsec */
+  Coords coords;
+  coords.crval1 = Rc;
+  coords.crval2 = Dc;
+  coords.crpix1 = 0;
+  coords.crpix2 = 0;
+  coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0;
+  coords.pc1_1  = coords.pc2_2 = 1.0;
+  coords.pc1_2  = coords.pc2_1 = 0.0;
+  coords.Npolyterms = 1;
+  strcpy (coords.ctype, "DEC--SIN");
+
+  // convert all average R,D values to X,Y:
+  double *X, *Y;
+  ALLOCATE (X, double, Naverage);
+  ALLOCATE (Y, double, Naverage);
+  float Xmin = +10000.0, Ymin = +10000.0;
+  float Xmax = -10000.0, Ymax = -10000.0;
+  for (i = 0; i < Naverage; i++) {
+    X[i] = NAN;
+    Y[i] = NAN;
+    // skip any stars which are outside of nominal catalog range
+    if (subcatalog[0].average[i].R < Rmin) continue;
+    if (subcatalog[0].average[i].R > Rmax) continue;
+    if (subcatalog[0].average[i].D < Dmin) continue;
+    if (subcatalog[0].average[i].D > Dmax) continue;
+    RD_to_XY (&X[i], &Y[i], subcatalog[0].average[i].R, subcatalog[0].average[i].D, &coords);
+    Xmin = MIN (Xmin, X[i]);
+    Xmax = MAX (Xmax, X[i]);
+    Ymin = MIN (Ymin, Y[i]);
+    Ymax = MAX (Ymax, Y[i]);
+  }
+
+  // how many grid cells? what is the grid spacing? 
+  float dX = Xmax - Xmin;
+  float dY = Ymax - Ymin;
+
+  // *** XXX for the moment, I'm using a hard-wired cell size (200 arcsec ~ 3.3 arcmin)
+  int NX = (int)(dX / 200) + 1;
+  int NY = (int)(dY / 200) + 1;
+  // fprintf (stderr, "Density Grid: %d x %d\n", NX, NY);
+  // XXX check that NX,NY are sensible (5 degrees / 200 arcsec seems like the absolute max)
+  if (NX > 1000) { 
+    fprintf (stderr, "serious problem with %s: NX = %d\n", subcatalog[0].filename, NX); 
+    exit (3); 
+  }
+  if (NY > 1000) { 
+    fprintf (stderr, "serious problem with %s: NY = %d\n", subcatalog[0].filename, NY); 
+    exit (3); 
+  }
+  
+  // kind of ugly : generate a grid of index, Nmeasure arrays
+  // to be filled below (I also need NN and Nn to track the number of 
+  // entries in each).
+  int    **NN_grid;
+  int    **Nn_grid;
+  int   ***Nm_grid;
+  off_t ***idxgrid;
+  ALLOCATE (NN_grid, int *, NX);
+  ALLOCATE (Nn_grid, int *, NX);
+  ALLOCATE (Nm_grid, int **, NX);
+  ALLOCATE (idxgrid, off_t **, NX);
+
+  for (ix = 0; ix < NX; ix++) {
+    ALLOCATE (NN_grid[ix], int, NY);
+    ALLOCATE (Nn_grid[ix], int, NY);
+    ALLOCATE (Nm_grid[ix], int *, NY);
+    ALLOCATE (idxgrid[ix], off_t *, NY);
+    for (iy = 0; iy < NY; iy++) {
+      Nn_grid[ix][iy] = 0;
+      NN_grid[ix][iy] = 100;
+      ALLOCATE (Nm_grid[ix][iy], int,   NN_grid[ix][iy]);
+      ALLOCATE (idxgrid[ix][iy], off_t, NN_grid[ix][iy]);
+    }
+  }
+
+  // assign all of the average entries to a grid cell
+  for (i = 0; i < Naverage; i++) {
+    if (isnan(X[i])) continue;
+    if (isnan(Y[i])) continue;
+    ix = (X[i] - Xmin) / 200.0;
+    iy = (Y[i] - Ymin) / 200.0;
+    int Nn = Nn_grid[ix][iy];
+    Nm_grid[ix][iy][Nn] = subcatalog[0].average[i].Nmeasure;
+    
+    // if we are resetting to a given photcode, we need to have that photcode...
+    if (NphotcodesReset) {
+      int k;
+      int foundReset = FALSE;
+      int m = subcatalog[0].average[i].measureOffset;
+      MeasureTiny *measure = &subcatalog[0].measureT[m];
+      for (j = 0; (j < subcatalog[0].average[i].Nmeasure) && !foundReset; j++) {
+        if (CHECK_TEST_IMAGE && (abs(measure[j].t - Tref) < 10) && (measure[j].photcode == Cref)) {
+          fprintf (stderr, ".");
+        }
+        for (k = 0; (k < NphotcodesReset) && !foundReset; k++) {
+          if (photcodesReset[k][0].code == measure[j].photcode) foundReset = TRUE;
+        }
+      }
+      if (!foundReset) {
+        Nm_grid[ix][iy][Nn] = 0;
+      }
+    }
+
+    idxgrid[ix][iy][Nn] = i;
+    Nn_grid[ix][iy] ++;
+    if (Nn_grid[ix][iy] >= NN_grid[ix][iy]) {
+      NN_grid[ix][iy] += 100;
+      REALLOCATE (Nm_grid[ix][iy], int,   NN_grid[ix][iy]);
+      REALLOCATE (idxgrid[ix][iy], off_t, NN_grid[ix][iy]);
+    }
+  }
+    
+  // sort all of the grid cells
+  for (ix = 0; ix < NX; ix++) {
+    for (iy = 0; iy < NY; iy++) {
+      sort_by_Nmeasure (Nm_grid[ix][iy], idxgrid[ix][iy], Nn_grid[ix][iy]);
+      NN_grid[ix][iy] = 0; // I'm going to use this array to track which element I've already selected
+    }
+  }
+
+  // cycle over the grid until we ready Nmax
+  off_t *keepidx = NULL;
+  ALLOCATE (keepidx, off_t, Naverage);
+  memset (keepidx, 0, Naverage*sizeof(off_t));
+  int Nkeep = 0;
+
+  for (i = 0; (i < 20) && (Nkeep < Nmax); i++) {
+    for (ix = 0; (ix < NX) && (Nkeep < Nmax); ix++) {
+      for (iy = 0; (iy < NY) && (Nkeep < Nmax); iy++) {
+        if (NN_grid[ix][iy] >= Nn_grid[ix][iy]) continue; // all used up!
+        int Nn = NN_grid[ix][iy];
+        keepidx[Nkeep] = idxgrid[ix][iy][Nn];
+        Nkeep ++;
+        NN_grid[ix][iy] ++;
+      }
+    }
+  }
+
+  // count the number of measurements this selection will yield
+  off_t ave, NMEASURE = 0;
+  for (i = 0; i < Nkeep; i++) {
+    ave = keepidx[i];
+    NMEASURE += subcatalog[0].average[ave].Nmeasure;
+  }
+
+  // test catID : 37262 37261 37257 37258
+  int dumpit = FALSE;
+  dumpit |= (oldcatalog[0].catID == 37007);
+  // dumpit |= (oldcatalog[0].catID == 37261);
+  // dumpit |= (oldcatalog[0].catID == 37257);
+  // dumpit |= (oldcatalog[0].catID == 37258);
+  if (dumpit) {
+    char name[64];
+    snprintf (name, 64, "cat.%05d.dump.dat", oldcatalog[0].catID);
+    FILE *fdump = fopen (name, "w");
+    for (i = 0; i < Nkeep; i++) {
+      ave = keepidx[i];
+      fprintf (fdump, "%10.6f %10.6f %d\n", subcatalog[0].average[ave].R, subcatalog[0].average[ave].D, subcatalog[0].average[ave].Nmeasure);
+    }
+    fclose (fdump);
+  }
+
+  // allocate the output data 
+  ALLOCATE (tmpcatalog.average,  Average,     Nkeep);
+  ALLOCATE (tmpcatalog.measureT, MeasureTiny, NMEASURE);
+  ALLOCATE (tmpcatalog.secfilt,  SecFilt,     Nkeep * Nsecfilt);
+
+  off_t Nmeasure = 0;
+
+  // copy the Nkeep selected entries from subcatalog to tmpcatalog (adjusting links)
+  for (i = 0; i < Nkeep; i++) {
+    ave = keepidx[i];
+    tmpcatalog.average[i] = subcatalog[0].average[ave];
+    tmpcatalog.average[i].measureOffset = Nmeasure;
+    for (j = 0; j < tmpcatalog.average[i].Nmeasure; j++) {
+      off_t offset = subcatalog[0].average[ave].measureOffset + j;
+      tmpcatalog.measureT[Nmeasure] = subcatalog[0].measureT[offset];
+      tmpcatalog.measureT[Nmeasure].averef = i;
+      Nmeasure ++;
+    }
+    for (j = 0; j < Nsecfilt; j++) {
+      tmpcatalog.secfilt[i*Nsecfilt + j] = subcatalog[0].secfilt[ave*Nsecfilt + j];
+    }
+  }
+
+  if (VERBOSE2) {
+    char *basename = filebasename (oldcatalog[0].filename);
+    fprintf (stderr, "limited to %d ("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", 
+             Nkeep, subcatalog[0].Naverage, oldcatalog[0].Naverage, Nmeasure, subcatalog[0].Nmeasure,  oldcatalog[0].Nmeasure, basename);
+    free (basename);
+  }
+
+  free (X);
+  free (Y);
+
+  for (ix = 0; ix < NX; ix++) {
+    for (iy = 0; iy < NY; iy++) {
+      free (Nm_grid[ix][iy]);
+      free (idxgrid[ix][iy]);
+    }
+    free (NN_grid[ix]);
+    free (Nn_grid[ix]);
+    free (Nm_grid[ix]);
+    free (idxgrid[ix]);
+  }
+  free (NN_grid);
+  free (Nn_grid);
+  free (Nm_grid);
+  free (idxgrid);
+
+  free (keepidx);
+
+  free (subcatalog[0].average);
+  free (subcatalog[0].measureT);
+  free (subcatalog[0].secfilt);
+
+  subcatalog[0].average = tmpcatalog.average;
+  subcatalog[0].measureT = tmpcatalog.measureT;
+  subcatalog[0].secfilt = tmpcatalog.secfilt;
+  subcatalog[0].Naverage = Nkeep;
+  subcatalog[0].Nmeasure = Nmeasure;
+  subcatalog[0].Nsecfilt = oldcatalog[0].Nsecfilt;
+  subcatalog[0].Nsecf_mem = Naverage * oldcatalog[0].Nsecfilt;
+
+  return (TRUE);
 }
 
@@ -205 +490 @@
     NMEASURE += subcatalog[0].average[ave].Nmeasure;
   }
+
+# if (0)
+  if (oldcatalog[0].catID == 59962) {
+    FILE *fdump = fopen ("cat.dump.dat", "w");
+    for (i = 0; i < Nmax; i++) {
+      ave = index[i];
+      fprintf (fdump, "%10.6f %10.6f %d\n", subcatalog[0].average[ave].R, subcatalog[0].average[ave].D, subcatalog[0].average[ave].Nmeasure);
+    }
+    fclose (fdump);
+  }
+# endif
 
   // allocate the output data 

branches/eam_branches/ipp-20140423/Ohana/src/relastro/src/launch_region_hosts.c

@@ -88 +88 @@
     if (PHOTCODE_KEEP_LIST) strextend (command, "+photcode %s", PHOTCODE_KEEP_LIST); 
     if (PHOTCODE_SKIP_LIST) strextend (command, "-photcode %s", PHOTCODE_SKIP_LIST);
+    if (PHOTCODE_RESET_LIST) strextend (command, "-reset-to-photcode %s", PHOTCODE_RESET_LIST);
 
     if (MaxDensityUse)      strextend (command, "-max-density %f", MaxDensityValue);
@@ -97 +98 @@
     if (MinBadQF > 0.0)        strextend (command, "-min-bad-psfqf %f", MinBadQF);
     if (MaxMeanOffset != 10.0) strextend (command, "-max-mean-offset  %f", MaxMeanOffset);
+
+    strextend (command, "-D RELASTRO_SRC_MEAS_TOOFEW %d", SRC_MEAS_TOOFEW);
+    strextend (command, "-D RELASTRO_SIGMA_LIM %f", SIGMA_LIM);
+    strextend (command, "-D RELASTRO_DPOS_MAX %f", DPOS_MAX);
+    strextend (command, "-D ADDSTAR_RADIUS %d", ADDSTAR_RADIUS);
 
     if (TimeSelect) { 

branches/eam_branches/ipp-20140423/Ohana/src/relastro/src/relastro_images.c

@@ -34 +34 @@
   MARKTIME("load catalog data: %f sec\n", dtime);
 
-   photcodesKeep =  photcodesReset;
-  NphotcodesKeep = NphotcodesReset;
+  if (photcodesReset) {
+    photcodesKeep  = photcodesReset;
+    NphotcodesKeep = NphotcodesReset;
+  }
 
   if (Ncatalog == 0) {

branches/eam_branches/ipp-20140423/Ohana/src/relastro/src/relastro_parallel_images.c

@@ -48 +48 @@
   MARKTIME("-- load catalog data: %f sec\n", dtime);
   
+  if (photcodesReset) {
+    photcodesKeep  = photcodesReset;
+    NphotcodesKeep = NphotcodesReset;
+  }
+
   // NOTE: if I let all hosts load blindly, I saturate the data clients with too many
   // relastro_client requests.  I need to have the master mediate this.  the master