Changeset 39600

Timestamp:

Jun 16, 2016, 4:55:58 PM (10 years ago)

Author:

eugene

Message:

for repair mode: do not save unless changes are made; make a backup for changed catalogs; add header metadata to record modifications

Location:

trunk/Ohana/src/relastro

Files:

• include/relastro.h (modified) (3 diffs)
• src/RepairWarpMeasures.c (modified) (11 diffs)
• src/RepairWarps.c (modified) (2 diffs)

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

@@ -244 +244 @@
 } FrameCorrectionSet;
 
-
 typedef struct {
   double *Rave;
@@ -252 +251 @@
   int Nicrfobj;
 } ICRFobj;
+
+typedef struct {
+  int NfixChipID;
+  int NfixStackID;
+  int NfixWarpID;
+  int NfixWarpImageID;
+  int NfixWarpCoord;
+  int NmissWarp;
+  int NmissStack;
+  int NbadWarp;
+  int NbadWarpTime;
+  int NwarpNoImage;
+} WarpRepairResult;
 
 # define ID_MEAS_OBJECT_HAS_2MASS ID_MEAS_POOR_PHOTOM
@@ -802 +814 @@
 
 int RepairWarps (SkyList *skylist, int hostID, char *hostpath);
-int RepairWarpMeasures (Catalog *catalog);
+WarpRepairResult RepairWarpMeasures (Catalog *catalog);
 
 int FindWarpGroups (void);

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

@@ -2 +2 @@
 
 // XXX need to load_images first and generate some lookup tables
-int RepairWarpMeasures (Catalog *catalog) {
+WarpRepairResult RepairWarpMeasures (Catalog *catalog) {
 
   off_t Nimage;
@@ -9 +9 @@
   myAssert (!catalog->Nmeasure || catalog->measure, "programming error");
 
-  int NfixChipID = 0, NfixStackID = 0, NfixWarpID = 0, NfixWarpImageID = 0, NmissWarp = 0, NmissStack = 0, NbadWarp = 0, NwarpNoImage = 0;
+  WarpRepairResult result;
+  result.NfixChipID = 0;
+  result.NfixStackID = 0;
+  result.NfixWarpID = 0;
+  result.NfixWarpImageID = 0;
+  result.NfixWarpCoord = 0;
+  result.NmissWarp = 0;
+  result.NmissStack = 0;
+  result.NbadWarp = 0;
+  result.NbadWarpTime = 0;
+  result.NwarpNoImage = 0;
 
   int onePercent = catalog->Nmeasure / 100;
@@ -25 +35 @@
       if (extID != measure->extID) {
         measure->extID = extID;
-        NfixChipID ++;
+        result.NfixChipID ++;
       }
       continue;
@@ -34 +44 @@
       int im = getImageByID (measure->imageID);
       if (im < 0) {
-        if (NmissStack < 10) fprintf (stderr, "missing stack exposure: %f %f : %d %d\n", measure->R, measure->D, measure->imageID, measure->photcode);
-        NmissStack ++;
+        if (result.NmissStack < 10) fprintf (stderr, "missing stack exposure: %f %f : %d %d\n", measure->R, measure->D, measure->imageID, measure->photcode);
+        result.NmissStack ++;
         continue;
       }
@@ -41 +51 @@
       if (extID != measure->extID) {
         measure->extID = extID;
-        NfixStackID ++;
+        result.NfixStackID ++;
       }
       continue;
@@ -64 +74 @@
     double Dwrp = measure->D;
 
+    /* 
+       we have a warp detection.  X,Y are correct (never modified).  things which might be wrong:
+       * imageID (thus the Mcal, R,D, extID may all be broken)
+       */
+
+    int coordsDiffer = FALSE;
     if (USE_IMAGE_COORDS_FOR_REPAIR) {
       off_t idx = getImageByID (measure->imageID);
       if (idx == -1) {
         if (VERBOSE2) fprintf (stderr, "can't match detection to image?\n");
-        NwarpNoImage ++;
+        result.NwarpNoImage ++;
         continue;
       }
 
       Coords *imcoords = &image[idx].coords;
+
+      // save the original warp coordinates for comparison
+      double Rold = ohana_normalize_angle_to_midpoint(Rwrp, 180.0);
+      double Dold = ohana_normalize_angle_to_midpoint(Dwrp,   0.0);
 
       XY_to_RD (&Rwrp, &Dwrp, X, Y, imcoords);
@@ -78 +98 @@
       Rwrp = ohana_normalize_angle_to_midpoint(Rwrp, 180.0);
       Dwrp = ohana_normalize_angle_to_midpoint(Dwrp,   0.0);
+
+      // check if the new Rwrp, Dwrp is different from the current R,D:
+      float csdec = cos(Dold * RAD_DEG);
+
+      // find the ra,dec displacement in arcsec:
+      double dR = 3600.0*(Rwrp - Rold)*csdec;
+      double dD = 3600.0*(Dwrp - Dold);
+      
+      // skip detections which are within a small distance of the expected location
+      // NOTE: this actually works surprisingly well near the pole
+      double dPos = hypot(dR,dD);
+      if (dPos > 0.5) { 
+        // in this case, we should double-check measure.R,D values below
+        coordsDiffer = TRUE;
+      }
     }
 
@@ -93 +128 @@
     double dPos = hypot(dR,dD);
     if (dPos < 2.0) { 
+      // image coord agrees with average coord
+      if (coordsDiffer) {
+        // measure coord disagrees with image coord: replace with image coord
+        measure->R = Rwrp;
+        measure->D = Dwrp;
+        result.NfixWarpCoord ++;
+      }
+      // check that the current extID value is good for warp detection:
       int im = getImageByID (measure->imageID);
       uint64_t extID = CreatePSPSStackDetectionID (34, image[im].externID, measure->detID);
       if (extID != measure->extID) {
         measure->extID = extID;
-        NfixWarpID ++;
-      }
-      continue;
+        result.NfixWarpID ++;
+      }
+      // check match of obstime:
+      if ((measure->t != image[im].tzero) && (result.NbadWarpTime < 10)) {
+        fprintf (stderr, "inconsistent warp measurement times: %d vs %d for %f, %f (%d, %ld)", measure->t, image[im].tzero, Rwrp, Dwrp, measure->photcode, extID); 
+        result.NbadWarpTime ++;
+        continue;
+      }
     }
 
@@ -105 +153 @@
     int warpSeq = GetWarpSeq (image, measure->t, measure->photcode, Rave, Dave, X, Y);
     if (warpSeq < 0) {
-      if (NmissWarp < 10) fprintf (stderr, "missing warp exposure: %f %f : %d %d\n", Rave, Dave, measure->t, measure->photcode);
-      NmissWarp ++;
+      if (result.NmissWarp < 10) fprintf (stderr, "missing warp exposure: %f %f : %d %d\n", Rave, Dave, measure->t, measure->photcode);
+      result.NmissWarp ++;
       continue;
     }
@@ -127 +175 @@
     if (dPos > 5.0) {
       fprintf (stderr, "measurement still far from average location: %f %f vs %f %f : %d %d\n", Rave, Dave, Rnew, Dnew, measure->t, measure->photcode);
-      NbadWarp ++;
+      result.NbadWarp ++;
     }
 
@@ -137 +185 @@
     if (extID != measure->extID) {
       measure->extID = extID;
-      NfixWarpID ++;
-    }
-    NfixWarpImageID ++;
+      result.NfixWarpID ++;
+    }
+    result.NfixWarpImageID ++;
   }
 
   fprintf (stderr, "\n");
-  fprintf (stderr, "repaired %s : warp image ID %d : ext ID chip %d stack %d warp %d : missed warp %d stack %d : bad warp: %d : warp no image: %d\n", catalog->filename, NfixWarpImageID, NfixChipID, NfixStackID, NfixWarpID, NmissWarp, NmissStack, NbadWarp, NwarpNoImage);
-
-  return (TRUE);
+  fprintf (stderr, "repaired %s : warp image ID %d : ext ID chip %d stack %d warp %d : missed warp %d stack %d : bad warp: %d : warp no image: %d : badWarpTime: %d, fixWarpCoord: %d\n", 
+           catalog->filename, 
+           result.NfixWarpImageID, result.NfixChipID, result.NfixStackID, 
+           result.NfixWarpID, result.NmissWarp, result.NmissStack, 
+           result.NbadWarp, result.NwarpNoImage, result.NbadWarpTime, 
+           result.NfixWarpCoord);
+
+  return (result);
 }

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

@@ -45 +45 @@
 
     // update the detection coordinates using the new image parameters
-    RepairWarpMeasures (&catalog);
+    WarpRepairResult result = RepairWarpMeasures (&catalog);
 
+    // track number of corrections and only update if corrections are made
+    int Nmods = 0;
+    Nmods += result.NfixChipID;
+    Nmods += result.NfixStackID;
+    Nmods += result.NfixWarpID;
+    Nmods += result.NfixWarpImageID;
+    Nmods += result.NfixWarpCoord;
+    Nmods += result.NmissWarp;
+    Nmods += result.NmissStack;
+    Nmods += result.NbadWarp;
+    Nmods += result.NbadWarpTime;
+    Nmods += result.NwarpNoImage;
+    if (!Nmods) {
+      fprintf (stderr, "no changes to %s, no output\n", catalog.filename);
+      dvo_catalog_unlock (&catalog);
+      dvo_catalog_free (&catalog);
+      continue;
+    }
+
+    if (!UPDATE) {
+      dvo_catalog_unlock (&catalog);
+      dvo_catalog_free (&catalog);
+      continue;
+    }
+    
     char history[128];
     struct timeval now;
@@ -55 +80 @@
     free (moddate);
 
-    if (!UPDATE) {
-      dvo_catalog_unlock (&catalog);
-      dvo_catalog_free (&catalog);
-      continue;
+    // add metadata to define number of corrections
+    char line[128];
+    snprintf (line, 60, "ChipID %d, StackID %d, warpID %d", result.NfixChipID, result.NfixStackID, result.NfixWarpID);
+    gfits_modify (&catalog.header, "REPAIR_1", "%s", 1, line);
+    snprintf (line, 60, "WarpImageID %d, WarpCoord %d, WarpNoImage %d", result.NfixWarpImageID, result.NfixWarpCoord, result.NwarpNoImage);
+    gfits_modify (&catalog.header, "REPAIR_2", "%s", 1, line);
+    snprintf (line, 60, "missWarp %d, badWarp %d, badWarpTime %d, missStack %d", result.NmissWarp, result.NbadWarp, result.NbadWarpTime, result.NmissStack);
+    gfits_modify (&catalog.header, "REPAIR_3", "%s", 1, line);
+
+    // save backup of original cpm file
+    if (!dvo_catalog_backup (&catalog, "rp0", TRUE)) {
+      fprintf (stderr, "ERROR: failed to make backup for catalog %s\n", catalog.filename);
+      exit (1);
     }
-    
+
     // write the updated detections to disk
     save_catalogs (&catalog, 1);