Index: /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/Makefile =================================================================== --- /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/Makefile (revision 29778) +++ /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/Makefile (revision 29779) @@ -72,7 +72,13 @@ DVOREPAIR = \ $(SRC)/dvorepair.$(ARCH).o \ +$(SRC)/dvorepairFixCPT.$(ARCH).o \ +$(SRC)/dvorepairImagesVsMeasures.$(ARCH).o \ +$(SRC)/dvorepairDeleteImageList.$(ARCH).o \ +$(SRC)/dvorepairFixImages.$(ARCH).o \ $(SRC)/psps_ids.$(ARCH).o \ $(SRC)/LoadImages.$(ARCH).o \ +$(SRC)/ReadDeleteList.$(ARCH).o \ $(SRC)/match_image.$(ARCH).o \ +$(SRC)/help.$(ARCH).o \ $(SRC)/ImageOps.$(ARCH).o Index: /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/include/dvomerge.h =================================================================== --- /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/include/dvomerge.h (revision 29778) +++ /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/include/dvomerge.h (revision 29779) @@ -16,4 +16,7 @@ # include # include + +# define myAbort(MSG) { fprintf (stderr, "%s\n", MSG); abort(); } +# define myAssert(LOGIC,MSG) { if (!(LOGIC)) { fprintf (stderr, "%s\n", MSG); abort(); } } int VERBOSE; @@ -87,5 +90,17 @@ int dvo_update_image_IDs PROTO((IDmapType *IDmap, Catalog *catalog)); -Image *LoadImages (char *filename, off_t *Nimage); -Image *MatchImage (Image *image, off_t Nimage, unsigned int time, short int source, unsigned int imageID); -off_t match_image (Image *image, off_t Nimage, unsigned int T, short int S); +// dvorepair prototypes +Image *LoadImages PROTO((FITS_DB *db, char *filename, off_t *Nimage)); +Image *MatchImage PROTO((Image *image, off_t Nimage, unsigned int time, short int source, unsigned int imageID)); +off_t match_image PROTO((Image *image, off_t Nimage, unsigned int T, short int S)); + +int SaveImages PROTO((FITS_DB *oldDB, char *filename, Image *imagesOut, off_t Nout)); + +int dvorepairFixCPT PROTO((int argc, char **argv)); +int dvorepairImagesVsMeasures PROTO((int argc, char **argv)); +int dvorepairDeleteImageList PROTO((int argc, char **argv)); +int dvorepairFixImages PROTO((int argc, char **argv)); + +int *ReadDeleteList PROTO((char *filename, int *nindex)); +int RepairTableCPT PROTO((char *cptFilenameSrc, char *cptFilenameTgt, char *cpsFilenameSrc, char *cpsFilenameTgt, Measure *measure, off_t Nmeasure, Image *image, off_t Nimage, char catformat)); +void dvorepair_help PROTO((int argc, char **argv)); Index: /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/src/LoadImages.c =================================================================== --- /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/src/LoadImages.c (revision 29778) +++ /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/src/LoadImages.c (revision 29779) @@ -1,46 +1,23 @@ # include "dvomerge.h" -int dvoUseImageCache = 1; - -static char *lastFilename = NULL; -static time_t lastModified = 0; -static Image *imageCache = NULL; -static off_t cacheNimage = 0; - -static time_t getLastModified(char *filename); - -Image *LoadImages (char *filename, off_t *Nimage) { +Image *LoadImages (FITS_DB *db, char *filename, off_t *Nimage) { int status; Image *image; - FITS_DB db; + double ZeroPoint; - if (lastFilename) { - if (dvoUseImageCache && !strcmp(lastFilename, filename)) { - // Make sure the file hasn't changed since we loaded it - if (getLastModified(filename) == lastModified) { - *Nimage = cacheNimage; - return imageCache; - } - } - free(lastFilename); - lastFilename = NULL; - free(imageCache); - imageCache = NULL; - cacheNimage = 0; - lastModified = 0; - } + myAssert(filename, "programming errror"); + myAssert(Nimage, "programming errror"); + + gfits_db_init (db); + db->lockstate = LCK_SOFT; + db->timeout = 120.0; - - gfits_db_init (&db); - db.lockstate = LCK_SOFT; - db.timeout = 120.0; - - if (!gfits_db_lock (&db, filename)) { + if (!gfits_db_lock (db, filename)) { fprintf (stderr, "error opening image catalog %s (1)\n", filename); return (NULL); } - if (db.dbstate == LCK_EMPTY) { + if (db->dbstate == LCK_EMPTY) { fprintf (stderr, "note: image catalog is empty\n"); ALLOCATE (image, Image, 1); @@ -49,6 +26,6 @@ } - status = dvo_image_load (&db, TRUE, FALSE); - gfits_db_close (&db); + status = dvo_image_load (db, TRUE, FALSE); + gfits_db_close (db); if (!status) { @@ -57,34 +34,13 @@ } - image = gfits_table_get_Image (&db.ftable, Nimage, &db.swapped); + image = gfits_table_get_Image (&db->ftable, Nimage, &db->swapped); if (!image) { fprintf (stderr, "ERROR: failed to read images\n"); return (NULL); } - if (dvoUseImageCache && image) { - cacheNimage = *Nimage; - imageCache = image; - lastFilename = strdup(filename); - lastModified = getLastModified(filename); - } + + gfits_scan (&db->header, "ZERO_PT", "%lf", 1, &ZeroPoint); + SetZeroPoint (ZeroPoint); return (image); } - -static time_t getLastModified(char *filename) { - struct stat statbuf; - if (!stat(filename, &statbuf)) { - return statbuf.st_mtime; - } else { - return 0; - } -} - -void FreeImages(Image *images) { - if (!dvoUseImageCache && (images != NULL)) { - free(images); - } else { - // defer free until next LoadImages with a different or modified Images.dat - } -} - Index: /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/src/ReadDeleteList.c =================================================================== --- /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/src/ReadDeleteList.c (revision 29779) +++ /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/src/ReadDeleteList.c (revision 29779) @@ -0,0 +1,99 @@ +# include "dvomerge.h" + +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 +*/ + Index: /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/src/dvorepair.c =================================================================== --- /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/src/dvorepair.c (revision 29778) +++ /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/src/dvorepair.c (revision 29779) @@ -1,663 +1,13 @@ # include "dvomerge.h" - -# define myAbort(MSG) { fprintf (stderr, "%s\n", MSG); abort(); } -# 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]) -// * 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) -// * 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, 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 *cpmFilenameSrc = NULL; - char *cptFilenameSrc = NULL; - char *cpsFilenameSrc = NULL; - char *cpmFilenameTgt = NULL; - char *cptFilenameTgt = NULL; - char *cpsFilenameTgt = NULL; - - char *imageFilename = NULL; - - Header cpmHeaderPHU; - Header cpmHeaderTBL; - FTable cpmFtable; - - FILE *cpmFile = NULL; - - char catformat; - - if (argc != 3) { - fprintf (stderr, "USAGE: dvorepair (catdir) (deleteList)\n"); - fprintf (stderr, " catdir : database of interest\n"); - fprintf (stderr, " deleteList : table from dvorepair giving images to be deleted\n"); - exit (2); - } - - char *catdir = argv[1]; - char *delList = 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); - - // 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; - 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; - } - - 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 - insky = SkyTableLoadOptimal (catdir, NULL, NULL, FALSE, SKY_DEPTH_HST, VERBOSE); - myAssert(insky, "can't read SkyTable"); - SkyTableSetFilenames (insky, catdir, "cpt"); + dvorepair_help(argc, argv); - // XXX apply this...generate the subset matching the user-selected region - SkyRegion UserPatch; - - 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(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; - - /*** 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 ++; - } - } - - 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); - } - gfits_free_header (&cpmHeaderPHU); - gfits_free_header (&cpmHeaderTBL); - free (measure); - free (measureNew); - - Ncheck ++; - if (Ncheck % 1000 == 0) { - fprintf (stderr, "%s...", inlist[0].regions[i][0].name); - } - } - 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); + if (!strcmp(argv[1], "-fix-cpt")) dvorepairFixCPT(argc, argv); + if (!strcmp(argv[1], "-images-vs-measures")) dvorepairImagesVsMeasures(argc, argv); + if (!strcmp(argv[1], "-delete-image-list")) dvorepairDeleteImageList(argc, argv); + if (!strcmp(argv[1], "-fix-images")) dvorepairFixImages(argc, argv); + dvorepair_help(0, NULL); + exit (2); } - -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); -} - Index: /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/src/dvorepairDeleteImageList.c =================================================================== --- /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/src/dvorepairDeleteImageList.c (revision 29779) +++ /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/src/dvorepairDeleteImageList.c (revision 29779) @@ -0,0 +1,566 @@ +# include "dvomerge.h" + +// 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]) +// * 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) +// * 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 dvorepairDeleteImageList (int argc, char **argv) { + + FITS_DB db; // database handle pointing to input image table + + off_t i, j, Nmeasure, NmeasureNew, Ndelete, Nvalid, Nimage, Nindex, *imageIdx, index; + int N, 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 *cpmFilenameSrc = NULL; + char *cptFilenameSrc = NULL; + char *cpsFilenameSrc = NULL; + char *cpmFilenameTgt = NULL; + char *cptFilenameTgt = NULL; + char *cpsFilenameTgt = NULL; + + char *imageFilename = NULL; + + Header cpmHeaderPHU; + Header cpmHeaderTBL; + FTable cpmFtable; + + FILE *cpmFile = NULL; + + char catformat; + + N = get_argument (argc, argv, "-delete-image-list"); + myAssert(N == 1, "programming error: -delete-image-list must be first from main"); + remove_argument (N, &argc, argv); + + if (argc != 3) { + fprintf (stderr, "USAGE: dvorepair -delete-image-list (catdir) (deleteList)\n"); + fprintf (stderr, " catdir : database of interest\n"); + fprintf (stderr, " deleteList : table from dvorepair giving images to be deleted\n"); + exit (2); + } + + char *catdir = argv[1]; + char *delList = argv[2]; + + // load the image data + ALLOCATE(imageFilename, char, strlen(catdir) + 12); + sprintf (imageFilename, "%s/Images.dat", catdir); + if ((image = LoadImages (&db, imageFilename, &Nimage)) == NULL) return (FALSE); + BuildChipMatch (image, Nimage); + + // 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; + 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; + } + + 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 + 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; + + 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(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; + + /*** 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 ++; + } + } + + 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); + } + gfits_free_header (&cpmHeaderPHU); + gfits_free_header (&cpmHeaderTBL); + free (measure); + free (measureNew); + + Ncheck ++; + if (Ncheck % 1000 == 0) { + fprintf (stderr, "%s...", inlist[0].regions[i][0].name); + } + } + 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 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); +} + Index: /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/src/dvorepairFixCPT.c =================================================================== --- /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/src/dvorepairFixCPT.c (revision 29779) +++ /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/src/dvorepairFixCPT.c (revision 29779) @@ -0,0 +1,190 @@ +# include "dvomerge.h" + +// 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 + +int dvorepairFixCPT (int argc, char **argv) { + + FITS_DB db; // database handle pointing to input image table + + off_t Nmeasure, Nimage; + int i, N, Nbytes, NaveMax, Naverage, NAVERAGE, Nave, Nold; + int *found; + + Image *image, *thisImage; + Average *average; + Measure *measure; + + Matrix matrix; + + char *cpmFilename, *cptFilenameSrc, *cptFilenameTgt, *imageFilename; + + Header cptHeaderPHU, cptHeaderTBL; + Header cpmHeaderPHU, cpmHeaderTBL; + FTable cpmFtable, cptFtable; + + FILE *cptFileSrc = NULL; + FILE *cptFileTgt = NULL; + FILE *cpmFile = NULL; + + char catformat; + + N = get_argument (argc, argv, "-fix-cpt"); + myAssert(N == 1, "programming error: -fix-cpt must be first from main"); + remove_argument (N, &argc, argv); + + if (argc != 5) { + fprintf (stderr, "USAGE: dvorepair -fix-cpt (images) (cpm) (cptInput) (cptOutput)\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.. + + // load the image data + if ((image = LoadImages (&db, imageFilename, &Nimage)) == NULL) return (FALSE); + BuildChipMatch (image, Nimage); + + // open cpm file + cpmFile = fopen(cpmFilename, "r"); + 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, FALSE)) { + myAbort("can't read data for cpm table"); + } + + 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); + } + 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); +} Index: /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/src/dvorepairFixImages.c =================================================================== --- /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/src/dvorepairFixImages.c (revision 29779) +++ /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/src/dvorepairFixImages.c (revision 29779) @@ -0,0 +1,289 @@ +# include "dvomerge.h" + +// 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]) +// * create an index of the imageIDs to be deleted (delete (T/F) = imageDELindex[i]) +// * XXX need to find the DIS images for the WRP images to be deleted +// ** BuildChipMatch +// ** copy ChipMatch, create index +// ** sort index, ChipMatch by ChipMatch +// ** loop over ChipMatch, counting delete/no delete for uniq values + +void SortChipMatch (off_t *M, off_t *I, off_t N); +int MarkMosaicsToDelete(Image *image, off_t Nimage, int *deleteIndex); + +int dvorepairFixImages (int argc, char **argv) { + + FITS_DB db; // database handle pointing to input image table + + off_t i, Nimage, Nout, Nindex, *imageIdx, index; + + int N; + + Image *image, *imageOut; + + char *imageFilenameOld = NULL; + char *imageFilenameNew = NULL; + + N = get_argument (argc, argv, "-fix-images"); + myAssert(N == 1, "programming error: -fix-images must be first from main"); + remove_argument (N, &argc, argv); + + if (argc != 3) { + fprintf (stderr, "USAGE: dvorepair -fix-images (catdir) (deleteList)\n"); + fprintf (stderr, " catdir : database of interest\n"); + fprintf (stderr, " deleteList : table from dvorepair giving images to be deleted\n"); + exit (2); + } + + char *catdir = argv[1]; + char *delList = argv[2]; + + // load the image data + ALLOCATE(imageFilenameOld, char, strlen(catdir) + 12); + sprintf (imageFilenameOld, "%s/Images.dat", catdir); + + ALLOCATE(imageFilenameNew, char, strlen(catdir) + 24); + sprintf (imageFilenameNew, "%s/Images.dat.fixed", catdir); + + if ((image = LoadImages (&db, imageFilenameOld, &Nimage)) == NULL) { + fprintf (stderr, "error loading images\n"); + exit (1); + } + BuildChipMatch (image, Nimage); + + // 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; + 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; + } + + MarkMosaicsToDelete(image, Nimage, deleteIndex); + + // now create a new output Image table with only the non-deletions + ALLOCATE (imageOut, Image, Nimage); + + Nout = 0; + for (i = 0; i < Nimage; i++) { + index = image[i].imageID; + if (deleteIndex[index]) continue; + imageOut[Nout] = image[i]; + Nout++; + } + + SaveImages(&db, imageFilenameNew, imageOut, Nout); + + free (imageFilenameOld); + free (imageFilenameNew); + free (imageIdx); + free (deleteIndex); + free (deleteIDs); + + exit (0); +} + +// sort two times vectors and an index by first time vector +void SortChipMatch (off_t *M, off_t *I, off_t N) { + +# define SWAPFUNC(A,B){ off_t tmp_m; off_t tmp_i; \ + tmp_m = M[A]; M[A] = M[B]; M[B] = tmp_m; \ + tmp_i = I[A]; I[A] = I[B]; I[B] = tmp_i; \ + } +# define COMPARE(A,B)(M[A] < M[B]) + OHANA_SORT (N, COMPARE, SWAPFUNC); +# undef SWAPFUNC +# undef COMPARE +} + +int MarkMosaicsToDelete(Image *image, off_t Nimage, int *deleteIndex) { + + off_t i; + + off_t *ChipByMosaicID, *ChipIndex, *mosaicIDs; + off_t mosaic, Nmosaic, NMOSAIC, chipSeq, chipID, mosaicID; + int *nChildren, *nToDelete; + int Na, Nb, Nc; + + /* We now have a table of images to delete by image ID : but only WRP images have been + marked. Now find all WRP images for which none/some/all images are to be deleted. + */ + + // make a copy of the ChipMatch table (generated by BuildChipMatch, see libdvo/src/mosaic_astrom.c) + // generate an associated image index and sort by the mosaic IDs (value of ChipMatch table) + ALLOCATE(ChipByMosaicID, off_t, Nimage); + off_t *ChipMatchRaw = GetChipMatch(); + myAssert(ChipMatchRaw, "need to run BuildChipMatch first"); + + memcpy (ChipByMosaicID, ChipMatchRaw, sizeof(off_t)*Nimage); + ALLOCATE (ChipIndex, off_t, Nimage); + for (i = 0; i < Nimage; i++) { + ChipIndex[i] = i; + } + SortChipMatch(ChipByMosaicID, ChipIndex, Nimage); + + // now count the number of children images with none/some/all marked for deletion + // we don't know the number of mosaic images, so allocate as we go + + Nmosaic = 0; + NMOSAIC = 1000; + ALLOCATE (nChildren, int, NMOSAIC); memset(&nChildren[Nmosaic], 0, (NMOSAIC-Nmosaic)*sizeof(int)); + ALLOCATE (nToDelete, int, NMOSAIC); memset(&nToDelete[Nmosaic], 0, (NMOSAIC-Nmosaic)*sizeof(int)); + ALLOCATE (mosaicIDs, off_t, NMOSAIC); memset(&mosaicIDs[Nmosaic], 0, (NMOSAIC-Nmosaic)*sizeof(off_t)); + + // scan past the mosaic (-3), unassigned (-2) and non-chip (-1) entries + for (i = 0; (i < Nimage) && (ChipByMosaicID[i] == -3); i++); + Na = i; + fprintf (stderr, "skipping %d DIS entries\n", Na); + + for (; (i < Nimage) && (ChipByMosaicID[i] == -2); i++); + Nb = i - Na; + fprintf (stderr, "skipping %d failed entries\n", Nb); + + for (; (i < Nimage) && (ChipByMosaicID[i] == -1); i++); + Nc = i - Na - Nb; + fprintf (stderr, "skipping %d single-chip entries\n", Nc); + + if (i == Nimage) { + fprintf (stderr, "there are no mosaic images in this database\n"); + free (nChildren); + free (nToDelete); + free (mosaicIDs); + + free (ChipIndex); + free (ChipByMosaicID); + + return TRUE; + } + + // mosaicIDs[] actually contains the sequence in images[] of the mosaic images + // do not confuse the sequence numbers and the imageID values + + mosaic = ChipByMosaicID[i]; + mosaicIDs[Nmosaic] = mosaic; + for (; i < Nimage; i++) { + if (ChipByMosaicID[i] != mosaic) { + // time for the next entry + Nmosaic ++; + if (Nmosaic == NMOSAIC - 1) { + NMOSAIC += 1000; + REALLOCATE (nChildren, int, NMOSAIC); memset(&nChildren[Nmosaic], 0, (NMOSAIC-Nmosaic)*sizeof(int)); + REALLOCATE (nToDelete, int, NMOSAIC); memset(&nToDelete[Nmosaic], 0, (NMOSAIC-Nmosaic)*sizeof(int)); + REALLOCATE (mosaicIDs, off_t, NMOSAIC); memset(&mosaicIDs[Nmosaic], 0, (NMOSAIC-Nmosaic)*sizeof(off_t)); + } + mosaic = ChipByMosaicID[i]; + mosaicIDs[Nmosaic] = mosaic; + } + nChildren[Nmosaic] ++; + chipSeq = ChipIndex[i]; + chipID = image[chipSeq].imageID; + if (deleteIndex[chipID]) { + nToDelete[Nmosaic] ++; + } + } + Nmosaic ++; + + // XXX test: + for (i = 0; i < Nmosaic; i++) { + if (nToDelete[i]) { + fprintf (stderr, "mosaic %lld (%lld) : %s : %d children : %d to delete\n", + (long long) i, (long long) mosaicIDs[i], image[mosaicIDs[i]].name, nChildren[i], nToDelete[i]); + } + + // mark the mosaic for deletion IFF all children are to be deleted + if (nToDelete[i] == nChildren[i]) { + mosaic = mosaicIDs[i]; + mosaicID = image[mosaic].imageID; + deleteIndex[mosaicID] = TRUE; + } + } + + free (nChildren); + free (nToDelete); + free (mosaicIDs); + + free (ChipIndex); + free (ChipByMosaicID); + + return TRUE; +} + +int SaveImages(FITS_DB *oldDB, char *filename, Image *imageOut, off_t Nout) { + + int status, Nx, nbytes; + off_t IDstart; + FITS_DB db; + + /* setup image table format and lock */ + db.mode = oldDB->mode; + db.format = oldDB->format; + status = dvo_image_lock (&db, filename, 3600.0, LCK_XCLD); // shorter timeout? + if (!status) { + fprintf (stderr, "ERROR: failure to lock image catalog %s\n", db.filename); + exit (1); + } + + /* load or create the image table */ + myAssert (db.dbstate == LCK_EMPTY, "do not overwrite exiting image table"); + myAssert (db.format == DVO_FORMAT_PS1_V2, "format mismatch"); + + dvo_image_create (&db, GetZeroPoint()); + + // replace the existing buffer with the image array + gfits_scan (db.ftable.header, "NAXIS1", "%d", 1, &Nx); + gfits_modify (db.ftable.header, "NAXIS2", OFF_T_FMT, 1, Nout); + db.ftable.header[0].Naxis[1] = Nout; + + nbytes = gfits_data_size (db.ftable.header); + db.ftable.datasize = nbytes; + db.ftable.buffer = (char *) imageOut; + REALLOCATE (db.ftable.buffer, char, MAX (nbytes, 1)); + memset (&db.ftable.buffer[Nx*Nout], ' ', nbytes - Nx*Nout); + + db.swapped = TRUE; // internal representation + + gfits_modify (&db.header, "NIMAGES", OFF_T_FMT, 1, Nout); + + status = gfits_scan (&oldDB->header, "IMAGEID", OFF_T_FMT, 1, &IDstart); + myAssert(status, "db image table lacks IMAGEID"); + + gfits_modify (&db.header, "IMAGEID", OFF_T_FMT, 1, IDstart); + + dvo_image_save (&db, TRUE); + + dvo_image_unlock (&db); + + return TRUE; +} + Index: /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/src/dvorepairImagesVsMeasures.c =================================================================== --- /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/src/dvorepairImagesVsMeasures.c (revision 29779) +++ /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/src/dvorepairImagesVsMeasures.c (revision 29779) @@ -0,0 +1,181 @@ +# include "dvomerge.h" + +// 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 dvorepairImagesVsMeasures (int argc, char **argv) { + + FITS_DB db; // database handle pointing to input image table + + off_t i, j, Nmeasure, Nvalid, Nimage, Nindex, *imageIdx, index, imageSeq; + int N, Nbytes, NbytesPerRow, Nbad, Ncheck, Ntol; + int *detCounts; + + Image *image; + Measure *measure; + + SkyTable *insky; + SkyList *inlist; + + char *catdir = NULL; + char *cpmFilename = NULL; + char *imageFilename = NULL; + + Header cpmHeaderPHU; + Header cpmHeaderTBL; + FTable cpmFtable; + + FILE *cpmFile = NULL; + + char catformat; + + N = get_argument (argc, argv, "-images-vs-measures"); + myAssert(N == 1, "programming error: -images-vs-measures must be first from main"); + remove_argument (N, &argc, argv); + + if (argc != 3) { + fprintf (stderr, "USAGE: dvorepair -images-vs-measures (catdir) (Ntol)\n"); + fprintf (stderr, " catdir : database of interest\n"); + fprintf (stderr, " Ntol : allow Ntol missing detections\n"); + exit (2); + } + + 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 (&db, imageFilename, &Nimage)) == NULL) return (FALSE); + BuildChipMatch (image, Nimage); + + // 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)); + + 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; + } + + // generate a list of the detection counts: + ALLOCATE (detCounts, int, Nimage); + memset (detCounts, 0, Nimage*sizeof(int)); + + // 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); + } + } + } + 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"); + } + + exit (0); +} Index: /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/src/help.c =================================================================== --- /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/src/help.c (revision 29778) +++ /branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/src/help.c (revision 29779) @@ -86,2 +86,26 @@ } +void dvorepair_help (int argc, char **argv) { + + /* check for help request */ + if (!argv) goto show_help; + if (get_argument (argc, argv, "-help")) goto show_help; + if (get_argument (argc, argv, "-h")) goto show_help; + if (argc < 2) goto show_help; + return; + +show_help: + + fprintf (stderr, "USAGE\n"); + fprintf (stderr, " dvorepair -fix-cpt (images) (cpm) (cptInput) (cptOutput) - regenerate a specific cpt file from the cpm file\n"); + fprintf (stderr, " dvorepair -images-vs-measures (catdir) (Ntol) - find images with too many missing detections\n"); + fprintf (stderr, " dvorepair -delete-image-list (catdir) (deleteList) - delete a set of images based on image IDs (output from -images-vs-measures)\n"); + fprintf (stderr, " dvorepair -fix-images (catdir) (deleteList) - delete a set of images based on image IDs (output from -images-vs-measures)\n"); + fprintf (stderr, "\n"); + + fprintf (stderr, " optional flags:\n"); + fprintf (stderr, " -help : this list\n"); + fprintf (stderr, " -h : this list\n\n"); + exit (2); +} + Index: /branches/eam_branches/ipp-20101103/Ohana/src/libdvo/include/dvo.h =================================================================== --- /branches/eam_branches/ipp-20101103/Ohana/src/libdvo/include/dvo.h (revision 29778) +++ /branches/eam_branches/ipp-20101103/Ohana/src/libdvo/include/dvo.h (revision 29779) @@ -293,4 +293,5 @@ int isRegisteredMosaic (void); off_t GetRegisteredMosaic (void); +off_t *GetChipMatch (); int GetMosaicCoords (Coords *coords); int FindMosaicForImage (Image *images, off_t Nimages, off_t entry); Index: /branches/eam_branches/ipp-20101103/Ohana/src/libdvo/src/mosaic_astrom.c =================================================================== --- /branches/eam_branches/ipp-20101103/Ohana/src/libdvo/src/mosaic_astrom.c (revision 29778) +++ /branches/eam_branches/ipp-20101103/Ohana/src/libdvo/src/mosaic_astrom.c (revision 29779) @@ -17,4 +17,9 @@ } +/* what is the currently registered mosaic image? */ +off_t *GetChipMatch () { + return (ChipMatch); +} + /* given an image array and a current entry of type WRP, find matching DIS & Register it */ int FindMosaicForImage (Image *images, off_t Nimages, off_t entry) { @@ -40,5 +45,5 @@ if (strcmp(&images[entry].coords.ctype[4], "-WRP")) { /* not a wrp image, do nothing */ - return (TRUE); /* error or not */ + return (entry + 1); /* error or not */ } @@ -52,5 +57,5 @@ RegisterMosaic (&images[i].coords); iDIS = i; - return (TRUE); + return (i + 1); } @@ -62,5 +67,5 @@ RegisterMosaic (&images[i].coords); iDIS = i; - return (TRUE); + return (i + 1); } return (FALSE); @@ -74,5 +79,5 @@ if (strcmp(&images[entry].coords.ctype[4], "-WRP")) { /* not a wrp image, do nothing */ - return (TRUE); /* error or not? */ + return (entry + 1); /* error or not? */ } @@ -85,5 +90,5 @@ RegisterMosaic (&images[N].coords); iDIS = N; - return (TRUE); + return (N + 1); } @@ -165,8 +170,18 @@ SortDISindex (DIStzero, DISentry, Ndis); + // ChipMatch has a few possible values: + // -3 : image is a mosaic (DIS) + // -2 : image is an unassigned WRP image + // -1 : image is not a WRP or DIS images + // >= 0 : value is the mosaic seq number for this WRP image + /* find all matched WRP images */ ALLOCATE (ChipMatch, off_t, Nimages); for (i = 0; i < Nimages; i++) { ChipMatch[i] = -1; + if (!strcmp(&images[i].coords.ctype[4], "-DIS")) { + ChipMatch[i] = -3; + continue; + } if (strcmp(&images[i].coords.ctype[4], "-WRP")) continue; Index: /branches/eam_branches/ipp-20101103/Ohana/src/opihi/dvo/images.c =================================================================== --- /branches/eam_branches/ipp-20101103/Ohana/src/opihi/dvo/images.c (revision 29778) +++ /branches/eam_branches/ipp-20101103/Ohana/src/opihi/dvo/images.c (revision 29779) @@ -11,9 +11,9 @@ int images (int argc, char **argv) { - off_t i, Nimage; + off_t i, Nimage, Nmosaic; int j, status, InPic, leftside, *plist, TimeSelect, ByName; int WITH_MOSAIC, SOLO_MOSAIC, HIDDEN; time_t tzero, tend; - int N, NPTS, n, npts, Npts, kapa; + int N, NPTS, n, npts, Npts, kapa, *foundMosaic; Vector Xvec, Yvec; double r[8], d[8], x[8], y[8], Rmin, Rmax, Rmid, trange, Radius; @@ -132,4 +132,9 @@ BuildChipMatch (image, Nimage); + if (SOLO_MOSAIC && photcode) { + ALLOCATE(foundMosaic, int, Nimage); + memset(foundMosaic, 0, Nimage*sizeof(int)); + } + Rmin = graphmode.coords.crval1 - 180.0; Rmax = graphmode.coords.crval1 + 180.0; @@ -137,4 +142,5 @@ int DistortImage = wordhash ("-DIS"); + int ChipImage = wordhash ("-WRP"); int TriangleUp = wordhash ("TRP-"); int TriangleDn = wordhash ("TRM-"); @@ -150,5 +156,6 @@ if (ByName && strncmp (image[i].name, name, strlen(name))) continue; if (TimeSelect && ((image[i].tzero < tzero) || (image[i].tzero+image[i].trate*image[i].NY > tzero + trange))) continue; - if (!FindMosaicForImage (image, Nimage, i)) continue; + if (!(Nmosaic = FindMosaicForImage (image, Nimage, i))) continue; + Nmosaic --; // XXX kind of a hack: FindMosaicForImage returns 0 or the mosaic seq number + 1 if (photcode) { if ( photcodeEquiv && (photcode[0].code != GetPhotcodeEquivCodebyCode(image[i].photcode))) continue; @@ -161,4 +168,31 @@ typehash = wordhash (&image[i].coords.ctype[4]); + + if (photcode && SOLO_MOSAIC) { + // mosaic (DIS) images are not currently given a photcode : plot these via the WRP entries + /* DIS images represent a field, not a chip */ + if (typehash != ChipImage) continue; + if (foundMosaic[Nmosaic]) continue; + x[0] = -0.5*image[Nmosaic].NX; y[0] = -0.5*image[Nmosaic].NY; + x[1] = +0.5*image[Nmosaic].NX; y[1] = -0.5*image[Nmosaic].NY; + x[2] = +0.5*image[Nmosaic].NX; y[2] = +0.5*image[Nmosaic].NY; + x[3] = -0.5*image[Nmosaic].NX; y[3] = +0.5*image[Nmosaic].NY; + for (j = 0; j < Npts; j++) { + status = XY_to_RD (&r[j], &d[j], x[j], y[j], &image[Nmosaic].coords); + if (!status) break; + r[j] = ohana_normalize_angle (r[j]); + while (r[j] < Rmin) { r[j] += 360.0; } + while (r[j] > Rmax) { r[j] -= 360.0; } + if (j == 0) { + leftside = (r[j] < Rmid); + } + if (j > 0) { + if ( leftside && (r[j] > Rmid + 90)) { r[j] -= 360.0; } + if (!leftside && (r[j] < Rmid - 90)) { r[j] += 360.0; } + } + } + foundMosaic[Nmosaic] = TRUE; + goto plot_points; + } /* DIS images represent a field, not a chip */ @@ -270,4 +304,6 @@ if (Npts == 0) continue; + plot_points: + status = FALSE; for (j = 0; j < Npts; j++) { @@ -320,4 +356,5 @@ free (Yvec.elements.Flt); FreeImages (image); + free (foundMosaic); return (TRUE);