Changeset 29734

Timestamp:

Nov 10, 2010, 9:35:40 AM (15 years ago)

Author:

eugene

Message:

dvorepair: now able to remove identified detections by image ID

Location:

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

Files:

• dvorepair.c (modified) (5 diffs)
• dvorepairCPT.c (modified) (1 diff)
• dvorepairImageVsMeasure.c (added)

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

@@ -2 +2 @@
 
 # define myAbort(MSG) { fprintf (stderr, "%s\n", MSG); abort(); }
-# define myAssert(LOGIC,MSG) { if (!LOGIC) { fprintf (stderr, "%s\n", MSG); abort(); } }
-
-// find images which are missing detections:
+# define myAssert(LOGIC,MSG) { if (!(LOGIC)) { fprintf (stderr, "%s\n", MSG); abort(); } }
+
+// delete images based on a text table of the target IDs
+// * load the delete table file (validate format)
 // * load the Images.dat file
 // * create an index for imageID (seq = imageIDindex[i])
-// * scan over all catalog files
+// * create an index of the imageIDs to be deleted (delete (T/F) = imageDELindex[i])
+// * determine the full RA & DEC range of the images to be deleted
+// * scan over all catalog files in the specified RA & DEC range
 // * load the cpm file (pad if short, identify the padded section)
-// * loop over detections: increment detection count for each imageID
+// * create a new cpm file
+// * loop over detections : only keep those not from images to be deleted
+// * load the cpt file
+// * rebuild the cpt file or delete the matching entries?
+
+int *ReadDeleteList(char *filename, int *nindex);
+int RepairTableCPT(char *cptFilenameSrc, char *cptFilenameTgt, char *cpsFilenameSrc, char *cpsFilenameTgt, Measure *measure, off_t Nmeasure, Image *image, off_t Nimage, char catformat);
 
 int main (int argc, char **argv) {
 
-  off_t i, j, Nmeasure, Nvalid, Nimage, Nindex, *imageIdx, index, imageSeq;
-  int Nbytes, NbytesPerRow, Nbad, Ncheck, Ntol;
-  int *detCounts;
+  off_t i, j, Nmeasure, NmeasureNew, Ndelete, Nvalid, Nimage, Nindex, *imageIdx, index;
+  int seq, Nbytes, NbytesPerRow, Ncheck, nPass, raPass;
+  double Rthis, Dthis, Rmin, Rmax, Dmin, Dmax, Qthis, Qmin, Qmax, dR, dQ;
 
   Image *image;
   Measure *measure;
+  Measure *measureNew;
+
+  Matrix matrix;
 
   SkyTable *insky;
   SkyList *inlist;
 
-  char *catdir = NULL;
-  char *cpmFilename = NULL;
+  char *cpmFilenameSrc = NULL;
+  char *cptFilenameSrc = NULL;
+  char *cpsFilenameSrc = NULL;
+  char *cpmFilenameTgt = NULL;
+  char *cptFilenameTgt = NULL;
+  char *cpsFilenameTgt = NULL;
+
   char *imageFilename = NULL;
 
@@ -36 +53 @@
 
   if (argc != 3) {
-    fprintf (stderr, "USAGE: dvorepair (catdir) (Ntol)\n");
+    fprintf (stderr, "USAGE: dvorepair (catdir) (deleteList)\n");
     fprintf (stderr, "  catdir : database of interest\n");
-    fprintf (stderr, "  Ntol : allow Ntol missing detections\n");
+    fprintf (stderr, "  deleteList : table from dvorepair giving images to be deleted\n");
     exit (2);
   }
 
-  catdir = argv[1];
-  Ntol = atoi(argv[2]);
+  char *catdir = argv[1];
+  char *delList = argv[2];
 
   // load the image data
@@ -51 +68 @@
   BuildChipMatch (image, Nimage);
 
-  // generate an index for imageIDs:
+  // generate an index for imageIDs
+
+  // find the full range of the imageID values
   Nindex = 0;
   for (i = 0; i < Nimage; i++) {
     Nindex = MAX(image[i].imageID, Nindex);
   }
+
+  // create the index vector
   ALLOCATE (imageIdx, off_t, (Nindex + 1));
   memset (imageIdx, 0, (Nindex + 1)*sizeof(off_t));
 
+  // assign the imageIDs to the index vector
   for (i = 0; i < Nimage; i++) {
     index = image[i].imageID;
-    if (index == 0) {
-      fprintf (stderr, "?");
+    myAssert(index, "image index cannot be 0");
+    myAssert(!imageIdx[index], "stomping on an earlier index");
+    imageIdx[index] = i;
+  }
+
+  // generate an index of the deletion image
+  int *deleteIndex;
+  ALLOCATE (deleteIndex, int, Nindex + 1);
+  memset (deleteIndex, 0, (Nindex + 1)*sizeof(int));
+
+  int NdeleteIDs;
+  int *deleteIDs = ReadDeleteList(delList, &NdeleteIDs);
+
+  for (i = 0; i < NdeleteIDs; i++) {
+    index = deleteIDs[i];
+    myAssert(index > 0, "image index cannot be <= 0");
+    myAssert(index <= Nindex, "delete index out of range of real data");
+    myAssert(!deleteIndex[index], "stomping on an earlier index");
+    deleteIndex[index] = TRUE;
+  }
+
+  // find the RA & DEC range of the images we want to delete
+  Rmin = 360.0; 
+  Rmax = 0.0;
+  Dmin = +90.0;
+  Dmax = -90.0;
+  for (i = 0; i < NdeleteIDs; i++) {
+    index = deleteIDs[i];
+    seq = imageIdx[index];
+    if (!FindMosaicForImage (image, Nimage, seq)) {
+      fprintf (stderr, "cannot find mosaic for %s, skipping\n", image[seq].name);
       continue;
     }
-    if (imageIdx[index]) {
-      fprintf (stderr, "!");
-      continue;
-    }
-    imageIdx[index] = i;
-  }
-
-  // generate a list of the detection counts:
-  ALLOCATE (detCounts, int, Nimage);
-  memset (detCounts, 0, Nimage*sizeof(int));
+
+    XY_to_RD(&Rthis, &Dthis, 0, 0, &image[seq].coords);
+    Rmin = MIN(Rthis, Rmin);
+    Rmax = MAX(Rthis, Rmax);
+    Dmin = MIN(Dthis, Dmin);
+    Dmax = MAX(Dthis, Dmax);
+    Qthis = (Rthis < 180.0) ? Rthis : Rthis - 360.0;
+    Qmin = MIN(Qthis, Qmin);
+    Qmax = MAX(Qthis, Qmax);
+
+    fprintf (stderr, "image @ %f,%f\n", Rthis, Dthis);
+
+    XY_to_RD(&Rthis, &Dthis, image[seq].NX, 0, &image[seq].coords);
+    Rmin = MIN(Rthis, Rmin);
+    Rmax = MAX(Rthis, Rmax);
+    Dmin = MIN(Dthis, Dmin);
+    Dmax = MAX(Dthis, Dmax);
+    Qthis = (Rthis < 180.0) ? Rthis : Rthis - 360.0;
+    Qmin = MIN(Qthis, Qmin);
+    Qmax = MAX(Qthis, Qmax);
+
+    XY_to_RD(&Rthis, &Dthis, 0, image[seq].NY, &image[seq].coords);
+    Rmin = MIN(Rthis, Rmin);
+    Rmax = MAX(Rthis, Rmax);
+    Dmin = MIN(Dthis, Dmin);
+    Dmax = MAX(Dthis, Dmax);
+    Qthis = (Rthis < 180.0) ? Rthis : Rthis - 360.0;
+    Qmin = MIN(Qthis, Qmin);
+    Qmax = MAX(Qthis, Qmax);
+
+    XY_to_RD(&Rthis, &Dthis, image[seq].NX, image[seq].NY, &image[seq].coords);
+    Rmin = MIN(Rthis, Rmin);
+    Rmax = MAX(Rthis, Rmax);
+    Dmin = MIN(Dthis, Dmin);
+    Dmax = MAX(Dthis, Dmax);
+    Qthis = (Rthis < 180.0) ? Rthis : Rthis - 360.0;
+    Qmin = MIN(Qthis, Qmin);
+    Qmax = MAX(Qthis, Qmax);
+  }
+
+  dQ = Qmax - Qmin;
+  dR = Rmax - Rmin;
 
   // load the sky table for the existing database
@@ -80 +163 @@
   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;
+
+  if (dR < dQ) {
+    fprintf (stderr, "R,D range: %f - %f, %f - %f\n", Rmin, Rmax, Dmin, Dmax);
+    nPass = 1;
+  } else {
+    fprintf (stderr, "R,D range: %f - %f, %f - %f\n", Qmin, Qmax, Dmin, Dmax);
+    nPass = 2;
+  }
   
-  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);
+  ALLOCATE(cpmFilenameSrc, char, strlen(catdir) + 64);
+  ALLOCATE(cptFilenameSrc, char, strlen(catdir) + 64);
+  ALLOCATE(cpsFilenameSrc, char, strlen(catdir) + 64);
+  ALLOCATE(cpmFilenameTgt, char, strlen(catdir) + 64);
+  ALLOCATE(cptFilenameTgt, char, strlen(catdir) + 64);
+  ALLOCATE(cpsFilenameTgt, char, strlen(catdir) + 64);
+
+  for (raPass = 0; raPass < nPass; raPass++) {
+    
+    UserPatch.Dmin = Dmin + 0.1;
+    UserPatch.Dmax = Dmax + 0.1;
+
+    if (dR < dQ) {
+      UserPatch.Rmin = Rmin;
+      UserPatch.Rmax = Rmax;
+    } else {
+      if (raPass == 0) {
+        UserPatch.Rmin = 0.0;
+        UserPatch.Rmax = Qmax;
+      } else {
+        UserPatch.Rmin = Qmin + 360.0;
+        UserPatch.Rmax = 360.0;
+      }
+    }
+
+    inlist = SkyListByPatch (insky, -1, &UserPatch);
+  
+    // SkyListPopulatedRange (&Ns, &Ne, inlist, 0);
+    // depth = inlist[0].regions[Ns][0].depth;
+  
+    // loop over the populated input regions
+    Ncheck = 0;
+    for (i = 0; i < inlist[0].Nregions; i++) {
+      if (!inlist[0].regions[i][0].table) continue;
+
+      sprintf (cpmFilenameSrc, "%s/%s.cpm", catdir, inlist[0].regions[i][0].name);
+      sprintf (cptFilenameSrc, "%s/%s.cpt", catdir, inlist[0].regions[i][0].name);
+      sprintf (cpsFilenameSrc, "%s/%s.cps", catdir, inlist[0].regions[i][0].name);
+      sprintf (cpmFilenameTgt, "%s/%s.cpm.fixed", catdir, inlist[0].regions[i][0].name);
+      sprintf (cptFilenameTgt, "%s/%s.cpt.fixed", catdir, inlist[0].regions[i][0].name);
+      sprintf (cpsFilenameTgt, "%s/%s.cps.fixed", 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");
+      /*** read and examine the CPM file ***/
+      {
+        fprintf (stderr, "check %s\n", cpmFilenameSrc);
+
+        // open cpm file
+        cpmFile = fopen(cpmFilenameSrc, "r");
+        if (!cpmFile) continue;
+        // myAssert(cpmFile, "failed to open cpm file");
+    
+        // 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);
+
+        // close the input cpm file
+        fclose(cpmFile);
+
+        // allocate an output array of measures (to replace, if needed)
+        ALLOCATE (measureNew, Measure, Nvalid);
+
+        NmeasureNew = 0;
+        Ndelete = 0;
+
+        // examine all measurements: find ones that need to be deleted
+        for (j = 0; j < Nvalid; j++) {
+          index = measure[j].imageID;
+          myAssert(index, "measure is missing an image ID");
+
+          if (deleteIndex[index]) {
+            Ndelete ++;
+            continue;
+          }
+          measureNew[NmeasureNew] = measure[j];
+          NmeasureNew ++;
+        }
       }
 
-      // 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");
+      fprintf (stderr, "deleting %d of %d (keep %d) detections from %s -> %s\n", (int) Ndelete, (int) Nvalid, (int) NmeasureNew, cpmFilenameSrc, cpmFilenameTgt);
+
+      // if we actually want to delete any measurements, write out a new cpm file
+      if (Ndelete > 0) {
+
+        // the CPT and CPS tables need to be regenerated.  This must happen first because, in the process, we also update measure->averef
+        RepairTableCPT(cptFilenameSrc, cptFilenameTgt, cpsFilenameSrc, cpsFilenameTgt, measureNew, NmeasureNew, image, Nimage, catformat);
+
+        // convert internal to external format 
+        if (!MeasureToFtable (&cpmFtable, measureNew, NmeasureNew, catformat)) {
+          myAbort("trouble converting format");
+        }
+
+        // create and write the output file
+        cpmFile = fopen(cpmFilenameTgt, "w");
+        myAssert(cpmFile, "failed to open cpt file");
+        
+        // write PHU header
+        if (!gfits_fwrite_header (cpmFile, &cpmHeaderPHU)) {
+          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 (&cpmHeaderPHU, &matrix);
+        if (!gfits_fwrite_matrix  (cpmFile, &matrix)) {
+          myAbort("can't write primary matrix");
+        }
+        gfits_free_matrix (&matrix);
+        
+        // write the table data
+        if (!gfits_fwrite_ftable_range (cpmFile, &cpmFtable, 0, NmeasureNew, 0, NmeasureNew)) {
+          myAbort("can't write table data");
+        }
+        fclose (cpmFile);
       }
-
-      // 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);
+      free (measureNew);
 
       Ncheck ++;
@@ -160 +322 @@
       }
     }
-  }
-  fprintf (stderr, "\n");
-
-  Nbad = 0;
-  for (i = 0; i < Nimage; i++) {
-    // careful: off_t math does not do well with subtractions...
-    if (detCounts[i] + Ntol < image[i].nstar) {
-      fprintf (stdout, "image %s (%d) : %d vs %d\n", image[i].name, image[i].imageID, detCounts[i], image[i].nstar);
-      Nbad ++;
-    }
-  }
-  if (!Nbad) {
-    fprintf (stderr, "no bad images found\n");
-  }
+    fprintf (stderr, "\n");
+    SkyListFree(inlist);
+  }
+
+  free (imageFilename);
+  free (imageIdx);
+  free (deleteIndex);
+  free (deleteIDs);
+
+  free(cpmFilenameSrc);
+  free(cptFilenameSrc);
+  free(cpsFilenameSrc);
+  free(cpmFilenameTgt);
+  free(cptFilenameTgt);
+  free(cpsFilenameTgt);
+
+  SkyTableFree(insky);
 
   exit (0);
 }
+
+int *ReadDeleteList(char *filename, int *nindex) {
+
+  int j, index, Nindex, NINDEX, *indexList, Nline;
+  int Nstart, Nbytes, Nread, status;
+  char *c0, *c1, *space, *indexPoint, *buffer;
+
+  FILE *f = fopen (filename, "r");
+  myAssert(f, "failed to open delete list");
+  
+  Nindex = 0;
+  NINDEX = 1000;
+  ALLOCATE(indexList, int, NINDEX);
+
+  ALLOCATE (buffer, char, 0x10001);
+  bzero (buffer, 0x10001);
+
+  Nstart = 0;
+  Nline = 0;
+  while (1) {
+    Nbytes = 0x10000 - Nstart;
+    bzero (&buffer[Nstart], Nbytes);
+    Nread = fread (&buffer[Nstart], 1, Nbytes, f);
+    if (ferror (f)) {
+      perror ("error reading data file");
+      exit (1);
+    }
+    if (Nread == 0) break; // end of data
+    // nbytes = Nread + Nstart;
+    
+    status = TRUE;
+    c0 = buffer; 
+    while (status) {
+      // identify the range of the line
+      c1 = strchr (c0, '\n');
+      if (c1 == (char *) NULL) {
+        Nstart = strlen (c0);
+        memmove (buffer, c0, Nstart);
+        break;
+      } else {
+        *c1 = 0;
+      }      
+      if (*c0 == '#') goto next_line;
+
+      // confirm 'image' as first token:
+      if (strncmp(c0, "image", 5)) {
+        fprintf (stderr, "error on line %d: missing 'image' token\n", Nline);
+        goto next_line;
+      }
+
+      // confirm we have 7 fields broken by 6 spaces:
+      space = c0;
+      for (j = 0; j < 6; j++) {
+        space = strchr(space, ' '); 
+        if (!space) {
+          fprintf (stderr, "line only has %d word(s)\n", j + 1);
+          goto next_line;
+        }
+        space ++;
+        if (j == 1) {
+          indexPoint = space + 1;
+        }
+      }
+
+      // save the value we have found:
+      index = atoi(indexPoint);
+      indexList[Nindex] = index;
+
+      fprintf (stderr, "index: %d, line: %s\n", index, c0);
+
+      Nindex ++;
+      if (Nindex == NINDEX) {
+        NINDEX += 1000;
+        REALLOCATE (indexList, int, NINDEX);
+      }
+      
+    next_line:
+      Nline ++;
+      c0 = c1 + 1;
+    }
+  }
+
+  *nindex = Nindex;
+  return (indexList);
+}
+
+/* Delete List Format:
+   image o5252g0035o.140669.cm.51017.smf[XY11.hdr] (219974) : 368 vs 375
+   image o5252g0051o.140685.cm.51123.smf[XY11.hdr] (226417) : 350 vs 356
+   image o5399g0207o.194950.cm.98351.smf[XY46.hdr] (1683277) : 10341 vs 10347
+   image o5464g0284o.230438.cm.123774.smf[XY21.hdr] (2634913) : 214 vs 242
+   image o5464g0305o.230459.cm.123795.smf[XY21.hdr] (2635758) : 232 vs 257
+   image o5465g0306o.231171.cm.124201.smf[XY01.hdr] (2654941) : 0 vs 494
+   image o5465g0306o.231171.cm.124201.smf[XY02.hdr] (2654942) : 10 vs 813
+   image o5465g0306o.231171.cm.124201.smf[XY10.hdr] (2654947) : 0 vs 534
+   image o5465g0306o.231171.cm.124201.smf[XY11.hdr] (2654948) : 0 vs 447
+*/
+
+int RepairTableCPT(char *cptFilenameSrc, char *cptFilenameTgt, char *cpsFilenameSrc, char *cpsFilenameTgt, Measure *measure, off_t Nmeasure, Image *image, off_t Nimage, char catformat) {
+
+  off_t *averefMatch;
+  off_t i, NaveMax, Naverage, NAVERAGE, NaverageOut, Nave, Nout, Nold;
+  int *found, Nsecfilt;
+
+  Image *thisImage;
+  Average *average, *averageOut;
+
+  Matrix matrix;
+
+  Header cptHeaderPHU;
+  Header cptHeaderTBL;
+  FTable cptFtable;
+
+  FILE *cptFileSrc = NULL;
+  FILE *cptFileTgt = NULL;
+
+  cptFtable.header = &cptHeaderTBL;
+
+  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 (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));
+    }
+
+    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;
+    }
+
+    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
+
+    // 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);
+  }
+  Naverage = NaveMax + 1;
+
+  // we now have an average table, but there will be holes due to deleted measurements 
+  // create a new average table with only existing entries
+
+  ALLOCATE (averageOut, Average, Naverage);
+  memset (averageOut, 0, Naverage*sizeof(Average));
+
+  ALLOCATE (averefMatch, off_t, Naverage);
+  memset (averefMatch, 0, Naverage*sizeof(int));
+
+  Nave = 0;
+  for (i = 0; i < Naverage; i++) {
+    if (!found[i]) continue;
+    averageOut[Nave] = average[i]; // use a memcpy?
+    averefMatch[i] = Nave;
+    Nave ++;
+  }
+  NaverageOut = Nave;
+
+  for (i = 0; i < Nmeasure; i++) {
+    Nave = measure[i].averef;
+    Nout = averefMatch[Nave];
+    myAssert(Nout < NaverageOut, "output averef is wrong");
+    
+    myAssert(average[Nave].objID == measure[i].objID, "objIDs do not match");
+    myAssert(average[Nave].catID == measure[i].catID, "objIDs do not match");
+    myAssert(averageOut[Nout].objID == measure[i].objID, "objIDs do not match");
+    myAssert(averageOut[Nout].catID == measure[i].catID, "objIDs do not match");
+
+    measure[i].averef = Nout;
+  }
+
+  fprintf (stderr, "cpt file : %d obj -> %d obj (%s -> %s)\n", (int) Naverage, (int) NaverageOut, cptFilenameSrc, cptFilenameTgt);
+
+  // 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",     OFF_T_FMT, 1,  NaverageOut);
+  gfits_modify (&cptHeaderPHU, "NMEAS",      OFF_T_FMT, 1,  Nmeasure);
+  gfits_modify (&cptHeaderPHU, "NMISS",      "%d",      1,  0);
+  gfits_modify_alt (&cptHeaderPHU, "SORTED", "%t",      1,  FALSE);
+
+  gfits_scan (&cptHeaderPHU, "NSECFILT",     "%d",      1,  &Nsecfilt);
+
+  /* convert internal to external format */
+  if (!AverageToFtable (&cptFtable, averageOut, NaverageOut, 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, NaverageOut, 0, NaverageOut)) {
+    myAbort("can't write table data");
+  }
+  
+  fclose(cptFileTgt);
+  fclose(cptFileSrc);
+
+  gfits_free_header (&cptHeaderPHU);
+  gfits_free_header (&cptHeaderTBL);
+  free (average);
+  free (averageOut);
+
+  free (found);
+  free (averefMatch);
+
+  { 
+    Header cpsHeaderPHU;
+    Header cpsHeaderTBL;
+    FTable cpsFtable;
+
+    FILE *cpsFileSrc = NULL;
+    FILE *cpsFileTgt = NULL;
+
+    SecFilt *secfilt = NULL;
+
+    cpsFtable.header = &cpsHeaderTBL;
+
+    // open source cpt file
+    cpsFileSrc = fopen(cpsFilenameSrc, "r");
+    myAssert(cpsFileSrc, "failed to open cps file");
+    
+    // load the cps header (use for CATID, RA,DEC range, filenames)
+    if (!gfits_fread_header (cpsFileSrc, &cpsHeaderPHU)) {
+      myAbort("failure to cps header");
+    }
+
+    int Nrows = Nsecfilt*NaverageOut;
+    ALLOCATE (secfilt, SecFilt, Nrows);
+
+    /* convert internal to external format */
+    if (!SecFiltToFtable (&cpsFtable, secfilt, Nrows, catformat)) {
+      myAbort("trouble converting format");
+    }
+
+    // create and write the output file
+    cpsFileTgt = fopen(cpsFilenameTgt, "w");
+    myAssert(cpsFileTgt, "failed to open cps file");
+    
+    // write PHU header
+    if (!gfits_fwrite_header (cpsFileTgt, &cpsHeaderPHU)) {
+      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 (&cpsHeaderPHU, &matrix);
+    if (!gfits_fwrite_matrix  (cpsFileTgt, &matrix)) {
+      myAbort("can't write primary matrix");
+    }
+    gfits_free_matrix (&matrix);
+
+    // write the table data
+    if (!gfits_fwrite_ftable_range (cpsFileTgt, &cpsFtable, 0, Nrows, 0, Nrows)) {
+      myAbort("can't write table data");
+    }
+  
+    fclose(cpsFileTgt);
+    fclose(cpsFileSrc);
+
+    gfits_free_header (&cpsHeaderPHU);
+    gfits_free_header (&cpsHeaderTBL);
+    free (secfilt);
+  }
+
+  return (TRUE);
+}
+

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

@@ -2 +2 @@
 
 # define myAbort(MSG) { fprintf (stderr, "%s\n", MSG); abort(); }
-# define myAssert(LOGIC,MSG) { if (!LOGIC) { fprintf (stderr, "%s\n", MSG); abort(); } }
+# 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: