Changeset 37446

Timestamp:

Sep 29, 2014, 8:42:00 AM (12 years ago)

Author:

eugene

Message:

initial testing of the astrom map IO / application is now working in relastro

Location:

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

Files:

• include/relastro.h (modified) (2 diffs)
• src/FitChip.c (modified) (5 diffs)
• src/args.c (modified) (1 diff)
• src/high_speed_utils.c (modified) (2 diffs)
• src/load_images.c (modified) (3 diffs)
• src/relastro_images.c (modified) (2 diffs)

branches/eam_branches/ipp-20140904/Ohana/src/relastro/include/relastro.h

@@ -204 +204 @@
 int    PLOTDELAY;
 int    CHIPORDER;
+int    CHIPMAP;
 
 int MaxDensityUse;
@@ -547 +548 @@
 int isGPC1stack (int photcode);
 int isGPC1warp (int photcode);
+
+int save_astrom_table ();
+AstromOffsetTable *get_astrom_table ();
+
+int fit_map (AstromOffsetMap *map, StarData *raw, StarData *ref, int Npts);
+

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

@@ -57 +57 @@
     relastroVisualPlotChipFit(raw, ref, dRmax, Nmatch);
 
+    // for polynomial fits, I used (5,30,60) = (5,6,7) stars per term
+
     // set the maximum order for the polynomial (based on number of stars kept above)
     int order_use = 0;
-    if (Nkeep >  5) order_use = 1; // 4 stars per polynomial term (per dimension)
-    if (Nkeep > 30) order_use = 2; // 5 stars per polynomial term (per dimension)
-    if (Nkeep > 60) order_use = 3; // 6 stars per polynomial term (per dimension)
+    if (Nkeep >   5) order_use = 1; //  5 stars per cell
+    if (Nkeep >  24) order_use = 2; //  6 stars per cell
+    if (Nkeep >  63) order_use = 3; //  7 stars per cell
+    if (Nkeep > 128) order_use = 4; //  8 stars per cell
+    if (Nkeep > 225) order_use = 5; //  9 stars per cell
+    if (Nkeep > 360) order_use = 6; // 10 stars per cell
     if (order_use < 1) {
       if (VERBOSE2) fprintf (stderr, "insufficient measurements (%d) for linear fit\n", Nkeep);
@@ -67 +72 @@
       return FALSE;
     }
+    // fprintf (stderr, "using %d for %s\n", order_use, image[0].name);
 
     // when fitting the map, first fit a linear model (below? change Npolyterms to -1)
-    if (MAP) {
+    if (CHIPMAP) {
       image[0].coords.Npolyterms = 1;
     } else {
@@ -99 +105 @@
     }
 
+    // apply fit to get the fitted X,Y coordinates.  we need these to fit the residual map below
     for (i = 0; i < Nmatch; i++) {
-      // we have not yet fitted the map, so image[0].coords.Npolyterms needs to 
-      // be 1 here...
+      // we have not yet fitted the map, so Npolyterms needs to be 1 here:
+      LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[0].coords);
+    }
+
+    if (CHIPMAP) {
+      if (image[0].coords.offsetMap == NULL) {
+        // allocate a new table and assign to this image
+        AstromOffsetTable *table = get_astrom_table ();
+        AstromOffsetTableNewMap(table, order_use, image);
+      }
+      fit_map (image[0].coords.offsetMap, raw, ref, Nmatch);
+      image[0].coords.Npolyterms = -1;
+    }
+
+    for (i = 0; i < Nmatch; i++) {
+      // if we have fitted the map above, Npolyterms needs to be -1 here:
       XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[0].coords);
-    }
-
-    // where do we apply the fit all the way b
-
-    if (MAP) {
-      if (image[0].coords.offsetMap == NULL) {
-        AstromOffsetMap *map = AstromOffsetTableNewMap(table, order);
-        image[0].coords.offsetMap = map;
-      }
-      fit_map (image[0].coords.offsetMap, raw, ref);
-      image[0].coords.Npolyterm = -1;
     }
   }
@@ -274 +284 @@
 */
 
-int fit_map (AstromOffsetMap *map) {
+int fit_map (AstromOffsetMap *map, StarData *raw, StarData *ref, int Npts) {
 
   // we are actually fitting the residual after the linear fit has been taken off
@@ -280 +290 @@
   // fit the linear terms as above
   // calculate dX (raw.X - ref.X) and dY 
-  
-  AstromOffsetMapFit (map, x, y, dX, Npts, TRUE);
-  AstromOffsetMapFit (map, x, y, dY, Npts, FALSE);
+
+  int i, N;
+
+  float *x, *y, *dX, *dY;
+  ALLOCATE (x,  float, Npts);
+  ALLOCATE (y,  float, Npts);
+  ALLOCATE (dX, float, Npts);
+  ALLOCATE (dY, float, Npts);
+
+  N = 0;
+  for (i = 0; i < Npts; i++) {
+    if (raw[i].mask) continue;
+    x[N] = raw[i].X;
+    y[N] = raw[i].Y;
+    dX[N] = ref[i].X - raw[i].X;
+    dY[N] = ref[i].Y - raw[i].Y;
+    N++;
+  }
+
+  // in coordsops.c:XY_to_LM, the map is defined to carry dX,dY so that:
+  // (L,M) = f(X',Y') : (X',Y') = (X,Y) + (dX,dY)
+
+  AstromOffsetMapFit (map, x, y, dX, N, TRUE);
+  AstromOffsetMapFit (map, x, y, dY, N, FALSE);
+
+  free (x);
+  free (y);
+  free (dX);
+  free (dY);
 
   return TRUE;

branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/args.c

@@ -370 +370 @@
     CHIPORDER = atoi(argv[N]);
     remove_argument (N, &argc, argv);
+  }
+
+  CHIPMAP = FALSE;
+  if ((N = get_argument (argc, argv, "-chipmap"))) {
+    remove_argument (N, &argc, argv);
+    CHIPMAP = TRUE;
   }
 

branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/high_speed_utils.c

@@ -141 +141 @@
     off_t n;
     for (n = 0; n < NfieldsA; n++) {
-        valuesA[n] = dbExtractAverages (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], &catalog[0].measure[m], &fieldsA[n]);
+      valuesA[n] = dbExtractAverages (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], &catalog[0].measure[m], NULL, NULL, &fieldsA[n]);
     }
     return dbBooleanCond(stackA, NstackA, valuesA);
@@ -155 +155 @@
     off_t n;
     for (n = 0; n < NfieldsB; n++) {
-        valuesB[n] = dbExtractAverages (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], &catalog[0].measure[m], &fieldsB[n]);
+        valuesB[n] = dbExtractAverages (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], &catalog[0].measure[m], NULL, NULL, &fieldsB[n]);
     }
     return dbBooleanCond(stackB, NstackB, valuesB);

branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/load_images.c

@@ -1 +1 @@
 # include "relastro.h"
+
+static AstromOffsetTable *table = NULL;
+
+/* AstromOffsetMap is the map of Astrometry Offsets for non-polynomial astrometric
+ * corrections for each chip.  
+ */
 
 int load_images (FITS_DB *db, SkyList *skylist, int UseFullOverlap) {
@@ -17 +23 @@
   MARKTIME("  convert image table: %f sec\n", dtime);
 
+  // assign image->parent and image->coords.mosaic 
   BuildChipMatch (image, Nimage);
   MARKTIME("build chip match: %f sec\n", dtime);
+
+  char mapfile[DVO_MAX_PATH];
+  snprintf (mapfile, DVO_MAX_PATH, "%s/AstroMap.fits", CATDIR);
+  table = AstromOffsetMapLoad (mapfile, VERBOSE);
+
+  // assign images.coords.offsetMap -> table->map[i]
+  if (table) {
+    AstromOffsetTableMatchChips (image, Nimage, table);
+  } else {
+    table = AstromOffsetTableInit ();
+  }
 
   // select the images which overlap the selected sky regions
@@ -54 +72 @@
   return TRUE;
 }
+
+int save_astrom_table () {
+
+  char mapfile[DVO_MAX_PATH];
+  snprintf (mapfile, DVO_MAX_PATH, "%s/AstroMap.fits", CATDIR);
+  AstromOffsetMapSave (table, mapfile);
+
+  return TRUE;
+}
+
+AstromOffsetTable *get_astrom_table () {
+  return table;
+}
+

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

@@ -102 +102 @@
   if (PARALLEL) {
     // save the updated image parameters
+    // need to also save the image map table...
+    save_astrom_table ();
     dvo_image_update (&db, VERBOSE);
     dvo_image_unlock (&db); 
@@ -114 +116 @@
   if (!PARALLEL) {
     // save the updated image parameters
+    save_astrom_table ();
     dvo_image_update (&db, VERBOSE);
     dvo_image_unlock (&db);