Changeset 29714 for branches/eam_branches/ipp-20101103/Ohana/src/dvomerge

Timestamp:

Nov 8, 2010, 6:20:37 PM (16 years ago)

Author:

eugene

Message:

add the ability to read a broken fits table (pad the missing areas); adjust dvorepair to count missing detections per image

Location:

branches/eam_branches/ipp-20101103/Ohana/src/dvomerge

Files:

• doc (added)
• doc/notes.txt (added)
• src/dvorepair.c (modified) (1 diff)
• src/dvorepairCPT.c (added)

branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/src/dvorepair.c

@@ -4 +4 @@
 # define myAssert(LOGIC,MSG) { if (!LOGIC) { fprintf (stderr, "%s\n", MSG); abort(); } }
 
-// broken cpt file, valid cpm file: we can recover everything in the cpt file from the cpm file:
-// * load the full cpm file
-// * loop over detections
-// * if ave_ref is new: create a new object
-//   * determine RA & DEC from ext_id
-//   * obj_id, cat_id are defined in detection
-//   * 
+// find images which are missing detections:
+// * load the Images.dat file
+// * create an index for imageID (seq = imageIDindex[i])
+// * scan over all catalog files
+// * load the cpm file (pad if short, identify the padded section)
+// * loop over detections: increment detection count for each imageID
 
 int main (int argc, char **argv) {
 
-  off_t Nmeasure, Nimage;
-  int i, Nbytes, NaveMax, Naverage, NAVERAGE, Nave, Nold;
-  int *found;
+  off_t i, j, Nmeasure, Nvalid, Nimage, Nindex, *imageIdx, index, imageSeq;
+  int Nbytes, NbytesPerRow, Nbad, Ncheck, Ntol;
+  int *detCounts;
 
-  Image *image, *thisImage;
-  Average *average;
+  Image *image;
   Measure *measure;
 
-  Matrix matrix;
+  SkyTable *insky;
+  SkyList *inlist;
 
-  char *cpmFilename, *cptFilenameSrc, *cptFilenameTgt, *imageFilename;
+  char *catdir = NULL;
+  char *cpmFilename = NULL;
+  char *imageFilename = NULL;
 
-  Header cptHeaderPHU, cptHeaderTBL;
-  Header cpmHeaderPHU, cpmHeaderTBL;
-  FTable cpmFtable, cptFtable;
+  Header cpmHeaderPHU;
+  Header cpmHeaderTBL;
+  FTable cpmFtable;
 
-  FILE *cptFileSrc = NULL;
-  FILE *cptFileTgt = NULL;
   FILE *cpmFile = NULL;
 
   char catformat;
 
-  if (argc != 5) {
-    fprintf (stderr, "USAGE: dvorepair (images) (cpm) (cptInput) (cptOutput)\n");
+  if (argc != 3) {
+    fprintf (stderr, "USAGE: dvorepair (catdir) (Ntol)\n");
+    fprintf (stderr, "  catdir : database of interest\n");
+    fprintf (stderr, "  Ntol : allow Ntol missing detections\n");
     exit (2);
   }
 
-  imageFilename  = argv[1];
-  cpmFilename    = argv[2];
-  cptFilenameSrc = argv[3];
-  cptFilenameTgt = argv[4];
-
-  cpmFtable.header = &cpmHeaderTBL;
-  cptFtable.header = &cptHeaderTBL;
-
-  // XXX don't bother locking for now: this is totally manual..
+  catdir = argv[1];
+  Ntol = atoi(argv[2]);
 
   // load the image data
+  ALLOCATE(imageFilename, char, strlen(catdir) + 12);
+  sprintf (imageFilename, "%s/Images.dat", catdir);
   if ((image = LoadImages (imageFilename, &Nimage)) == NULL) return (FALSE);
   BuildChipMatch (image, Nimage);
 
-  // open cpm file
-  cpmFile = fopen(cpmFilename, "r");
-  myAssert(cpmFile, "failed to open cpm file");
+  // generate an index for imageIDs:
+  Nindex = 0;
+  for (i = 0; i < Nimage; i++) {
+    Nindex = MAX(image[i].imageID, Nindex);
+  }
+  ALLOCATE (imageIdx, off_t, (Nindex + 1));
+  memset (imageIdx, 0, (Nindex + 1)*sizeof(off_t));
 
-  // load the cpm header
-  if (!gfits_fread_header (cpmFile, &cpmHeaderPHU)) {
-    myAbort("failure to cpm header");
+  for (i = 0; i < Nimage; i++) {
+    index = image[i].imageID;
+    if (index == 0) {
+      fprintf (stderr, "?");
+      continue;
+    }
+    if (imageIdx[index]) {
+      fprintf (stderr, "!");
+      continue;
+    }
+    imageIdx[index] = i;
   }
 
-  // move to TBL header
-  Nbytes = cpmHeaderPHU.datasize + gfits_data_size (&cpmHeaderPHU);
-  fseeko (cpmFile, Nbytes, SEEK_SET);
+  // generate a list of the detection counts:
+  ALLOCATE (detCounts, int, Nimage);
+  memset (detCounts, 0, Nimage*sizeof(int));
 
-  // read cpm TBL header
-  if (!gfits_fread_header (cpmFile, &cpmHeaderTBL)) { 
-    myAbort("can't read header for cpm table");
+  // load the sky table for the existing database
+  insky = SkyTableLoadOptimal (catdir, NULL, NULL, FALSE, SKY_DEPTH_HST, VERBOSE);
+  myAssert(insky, "can't read SkyTable");
+  SkyTableSetFilenames (insky, catdir, "cpt");
+
+  // XXX apply this...generate the subset matching the user-selected region
+  SkyRegion UserPatch;
+  UserPatch.Rmin = 0.0;
+  UserPatch.Rmax = 360.0;
+  UserPatch.Dmin = -90.0;
+  UserPatch.Dmax = +90.0;
+  inlist = SkyListByPatch (insky, -1, &UserPatch);
+
+  // SkyListPopulatedRange (&Ns, &Ne, inlist, 0);
+  // depth = inlist[0].regions[Ns][0].depth;
+  
+  ALLOCATE(cpmFilename, char, strlen(catdir) + 64);
+
+  // loop over the populated input regions
+  Ncheck = 0;
+  for (i = 0; i < inlist[0].Nregions; i++) {
+    if (!inlist[0].regions[i][0].table) continue;
+
+    if (1) {
+      sprintf (cpmFilename, "%s/%s.cpm", catdir, inlist[0].regions[i][0].name);
+      cpmFtable.header = &cpmHeaderTBL;
+
+      // open cpm file
+      cpmFile = fopen(cpmFilename, "r");
+      if (!cpmFile) continue;
+      // myAssert(cpmFile, "failed to open cpm file");
+    
+      // fprintf (stderr, "input: %s\n", inlist[0].regions[i][0].name);
+
+      // load the cpm header
+      if (!gfits_fread_header (cpmFile, &cpmHeaderPHU)) {
+        myAbort("failure to cpm header");
+      }
+
+      // move to TBL header
+      Nbytes = cpmHeaderPHU.datasize + gfits_data_size (&cpmHeaderPHU);
+      fseeko (cpmFile, Nbytes, SEEK_SET);
+
+      // read cpm TBL header
+      if (!gfits_fread_header (cpmFile, &cpmHeaderTBL)) { 
+        myAbort("can't read header for cpm table");
+      }
+
+      // read Measure table data : format is irrelevant here */
+      if (!gfits_fread_ftable_data (cpmFile, &cpmFtable, TRUE)) { 
+        myAbort("can't read data for cpm table");
+      }
+
+      gfits_scan(&cpmHeaderTBL, "NAXIS1", "%d", 1, &NbytesPerRow);
+
+      measure = FtableToMeasure (&cpmFtable, &Nmeasure, &catformat);
+      myAssert(measure, "failed to convert ftable to measure data");
+    
+      Nvalid = (int)(cpmFtable.validsize / NbytesPerRow);
+      Nvalid = MIN(Nmeasure, Nvalid);
+
+      // examine all measurements and new objects as needed
+      for (j = 0; j < Nvalid; j++) {
+        index = measure[j].imageID;
+        if (!index) {
+          fprintf (stderr, "?");
+          continue;
+        }
+
+        imageSeq = imageIdx[index];
+        // XXX check the range?
+
+        detCounts[imageSeq] ++;
+      }
+      fclose(cpmFile);
+      gfits_free_header (&cpmHeaderPHU);
+      gfits_free_header (&cpmHeaderTBL);
+      free (measure);
+
+      Ncheck ++;
+      if (Ncheck % 1000 == 0) {
+        fprintf (stderr, "%s...", inlist[0].regions[i][0].name);
+      }
+    }
   }
-  // read Measure table data : format is irrelevant here */
-  if (!gfits_fread_ftable_data (cpmFile, &cpmFtable)) { 
-    myAbort("can't read data for cpm table");
+  fprintf (stderr, "\n");
+
+  Nbad = 0;
+  for (i = 0; i < Nimage; i++) {
+    if (detCounts[i] < image[i].nstar - Ntol) {
+      fprintf (stdout, "image %s : %d vs %d\n", image[i].name, detCounts[i], image[i].nstar);
+      Nbad ++;
+    }
   }
-
-  measure = FtableToMeasure (&cpmFtable, &Nmeasure, &catformat);
-  myAssert(measure, "failed to convert ftable to measure data");
-
-  NaveMax = 0;
-  NAVERAGE = 1000;
-  ALLOCATE (average, Average, NAVERAGE);
-  memset (average, 0, NAVERAGE*sizeof(Average));
-
-  ALLOCATE (found, int, NAVERAGE);
-  memset (found, 0, NAVERAGE*sizeof(int));
-
-  // examine all measurements and new objects as needed
-  for (i = 0; i < Nmeasure; i++) {
-    Nave = measure[i].averef;
-    if (found[Nave]) {
-      average[Nave].Nmeasure ++;
-      myAssert(average[Nave].objID == measure[i].objID, "objIDs do not match!");
-      myAssert(average[Nave].catID == measure[i].catID, "catIDs do not match!");
-      continue;
-    }
-
-    if (Nave >= NAVERAGE) {
-      Nold = NAVERAGE;
-      NAVERAGE = MAX(Nave + 1000, NAVERAGE + 1000);
-      REALLOCATE (average, Average, NAVERAGE);
-      memset (&average[Nold], 0, (NAVERAGE - Nold)*sizeof(Average));
-
-      REALLOCATE (found, int, NAVERAGE);
-      memset (&found[Nold], 0, (NAVERAGE - Nold)*sizeof(int));
-    }
-
-    NaveMax = MAX(Nave, NaveMax);
-
-    found[Nave] = TRUE;
-
-    // we are going to leave most of the elements of average unset: they are the result of 
-    // the relastro analysis for this object and can be recreated with a call to relastro
-
-    // fields we have to set:
-
-    // need to find image so we can use ccd coordinates to determine RA & DEC
-    thisImage = MatchImage (image, Nimage, measure[i].t, measure[i].photcode, measure[i].imageID);
-    XY_to_RD (&average[Nave].R, &average[Nave].D, measure[i].Xccd, measure[i].Yccd, &thisImage[0].coords);
-
-    average[Nave].Nmeasure = 1;
-    average[Nave].Nmissing = 0;
-
-    // assume the resulting table set is unsorted
-    average[Nave].measureOffset = -1;
-    average[Nave].missingOffset = -1;
-    average[Nave].extendOffset = -1;
-
-    average[Nave].objID = measure[i].objID;
-    average[Nave].catID = measure[i].catID;
-    average[Nave].extID = CreatePSPSObjectID(average[Nave].R, average[Nave].D);
+  if (!Nbad) {
+    fprintf (stderr, "no bad images found\n");
   }
-  Naverage = NaveMax + 1;
-
-  // have we created all objects in the range 0 - Naverage?
-  for (i = 0; i < Naverage; i++) {
-    myAssert(found[i], "failed to find one");
-  }
-
-  // open source cpt file
-  cptFileSrc = fopen(cptFilenameSrc, "r");
-  myAssert(cptFileSrc, "failed to open cpt file");
-
-  // load the cpt header (use for CATID, RA,DEC range, filenames)
-  if (!gfits_fread_header (cptFileSrc, &cptHeaderPHU)) {
-    myAbort("failure to cpt header");
-  }
-
-  // update the output header
-  gfits_modify (&cptHeaderPHU, "NSTARS",     "%d",      1,  Naverage);
-  gfits_modify (&cptHeaderPHU, "NMEAS",      OFF_T_FMT, 1,  Nmeasure);
-  gfits_modify (&cptHeaderPHU, "NMISS",      "%d",      1,  0);
-  gfits_modify_alt (&cptHeaderPHU, "SORTED", "%t",      1,  FALSE);
-
-  /* convert internal to external format */
-  if (!AverageToFtable (&cptFtable, average, Naverage, catformat, NULL)) {
-    myAbort("trouble converting format");
-  }
-
-  // create and write the output file
-  cptFileTgt = fopen(cptFilenameTgt, "w");
-  myAssert(cptFileTgt, "failed to open cpt file");
-    
-  // write PHU header
-  if (!gfits_fwrite_header (cptFileTgt, &cptHeaderPHU)) {
-    myAbort("can't write primary header");
-  }
-
-  // write the PHU matrix; this is probably a NOP, do I have to keep it in?
-  gfits_create_matrix (&cptHeaderPHU, &matrix);
-  if (!gfits_fwrite_matrix  (cptFileTgt, &matrix)) {
-    myAbort("can't write primary matrix");
-  }
-  gfits_free_matrix (&matrix);
-
-  // write the table data
-  if (!gfits_fwrite_ftable_range (cptFileTgt, &cptFtable, 0, Naverage, 0, Naverage)) {
-    myAbort("can't write table data");
-  }
-  
-  fclose(cptFileTgt);
-  fclose(cptFileSrc);
-  fclose(cpmFile);
 
   exit (0);