Index: trunk/Ohana/src/dvomerge/Makefile =================================================================== --- trunk/Ohana/src/dvomerge/Makefile (revision 29181) +++ trunk/Ohana/src/dvomerge/Makefile (revision 29938) @@ -22,8 +22,13 @@ all: dvomerge dvoconvert dvosecfilt +# $(SRC)/dvomergeContinue.$(ARCH).o + DVOMERGE = \ $(SRC)/dvomerge.$(ARCH).o \ $(SRC)/dvomergeUpdate.$(ARCH).o \ $(SRC)/dvomergeCreate.$(ARCH).o \ +$(SRC)/dvomergeUpdate_threaded.$(ARCH).o \ +$(SRC)/dvomergeFromList.$(ARCH).o \ +$(SRC)/dvomergeImageIDs.$(ARCH).o \ $(SRC)/dvo_image_merge_dbs.$(ARCH).o \ $(SRC)/SetSignals.$(ARCH).o \ @@ -70,5 +75,29 @@ $(BIN)/dvosecfilt.$(ARCH) : $(DVOSECFILT) -INSTALL = dvomerge dvoconvert dvosecfilt +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 + +$(DVOREPAIR) : $(INC)/dvomerge.h + +$(BIN)/dvorepair.$(ARCH) : $(DVOREPAIR) + +DVOVERIFY = \ +$(SRC)/dvoverify.$(ARCH).o + +$(DVOVERIFY) : $(INC)/dvomerge.h + +$(BIN)/dvoverify.$(ARCH) : $(DVOVERIFY) + +INSTALL = dvomerge dvoconvert dvosecfilt dvorepair dvoverify # dependancy rules for binary code ######################### Index: trunk/Ohana/src/dvomerge/doc/notes.txt =================================================================== --- trunk/Ohana/src/dvomerge/doc/notes.txt (revision 29938) +++ trunk/Ohana/src/dvomerge/doc/notes.txt (revision 29938) @@ -0,0 +1,24 @@ + +2010.11.08 : DVO repairs + + * I have made a dvorepair program that reconstructs cpt files from a valid cpm file + + * the ThreePi.V1 database got busted; I think we have the following issues: + + * cpt files which are short (and can be fixed with the command above) + + * Image.dat has many image that have no matching detections (repeated dvomerge runs after the failure) + + * a number of cpm files which are short : in this case, I believe + these files have all detections before the failure, but are + missing some of the detections from the failure. + + * repair step 1 : scan image IDs on the database : count the number + of detections for each image ID and find those image IDs which are + missing detections + + * remove images that are missing their detections + + * create repaired cpm files excluding the now-orphaned detections + + Index: trunk/Ohana/src/dvomerge/include/dvomerge.h =================================================================== --- trunk/Ohana/src/dvomerge/include/dvomerge.h (revision 29181) +++ trunk/Ohana/src/dvomerge/include/dvomerge.h (revision 29938) @@ -17,4 +17,7 @@ # include +# define myAbort(MSG) { fprintf (stderr, "%s\n", MSG); abort(); } +# define myAssert(LOGIC,MSG) { if (!(LOGIC)) { fprintf (stderr, "%s\n", MSG); abort(); } } + int VERBOSE; char CATDIR[256]; @@ -25,4 +28,5 @@ float RADIUS; int SKY_DEPTH; +int NTHREADS; char *ALTERNATE_PHOTCODE_FILE; @@ -65,5 +69,5 @@ SkyList *SkyTablePopulatedList_old PROTO((SkyTable *sky, off_t Ns, off_t Ne)); -int LoadCatalog PROTO((Catalog *catalog, SkyRegion *region, char *filename, char *mode)); +int LoadCatalog PROTO((Catalog *catalog, SkyRegion *region, char *filename, char *mode, int Nsecfilt)); int merge_catalogs_new PROTO((SkyRegion *region, Catalog *output, Catalog *input, int *secflitMap)); @@ -86,2 +90,25 @@ off_t dvo_map_image_ID PROTO((IDmapType *IDmap, off_t oldID)); int dvo_update_image_IDs PROTO((IDmapType *IDmap, Catalog *catalog)); + +// 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)); + +off_t getTgtIndex PROTO((e_time start, e_time stop, short photcode, off_t *TgtIndex, e_time *TgtTimes, short *TgtCodes, off_t NimagesTgt)); +void SortTgtByTimes PROTO((e_time *S, off_t *I, short *C, off_t N)); +int dvo_image_match_dbs PROTO((IDmapType *IDmap, FITS_DB *tgt, FITS_DB *src)); +int dvomergeImagesGetMap PROTO((IDmapType *IDmap, char *input, char *output)); +int dvomergeFromList PROTO((int argc, char **argv)); +int dvomergeUpdate_threaded PROTO((int argc, char **argv)); Index: trunk/Ohana/src/dvomerge/src/ImageOps.c =================================================================== --- trunk/Ohana/src/dvomerge/src/ImageOps.c (revision 29938) +++ trunk/Ohana/src/dvomerge/src/ImageOps.c (revision 29938) @@ -0,0 +1,22 @@ +# include "dvomerge.h" + +Image *MatchImage (Image *image, off_t Nimage, unsigned int time, short int source, unsigned int imageID) { + + int m = -1; + + if ((imageID != 0) && (imageID < Nimage)) { + // imageID is in range for the array of images. If the table is still in order and + // no images have been deleted the index of the image we are looking for will be imageID - 1 + // If this is the case, we have it. Otherwise we'll have to go search for it below + int guess = (int) imageID - 1; + if (image[guess].imageID == imageID) { + m = guess; + } + } + if (m == -1) { + m = match_image (image, Nimage, time, source); + } + if (m == -1) return (NULL); + if (!FindMosaicForImage (image, Nimage, m)) return (NULL); + return (&image[m]); +} Index: trunk/Ohana/src/dvomerge/src/LoadCatalog.c =================================================================== --- trunk/Ohana/src/dvomerge/src/LoadCatalog.c (revision 29181) +++ trunk/Ohana/src/dvomerge/src/LoadCatalog.c (revision 29938) @@ -1,9 +1,9 @@ # include "dvomerge.h" -int LoadCatalog (Catalog *catalog, SkyRegion *region, char *filename, char *mode) { +int LoadCatalog (Catalog *catalog, SkyRegion *region, char *filename, char *mode, int Nsecfilt) { // set the parameters which guide catalog open/load/create catalog[0].filename = filename; - catalog[0].Nsecfilt = GetPhotcodeNsecfilt (); + catalog[0].Nsecfilt = Nsecfilt; // always load all of the data (if any exists) Index: trunk/Ohana/src/dvomerge/src/LoadImages.c =================================================================== --- trunk/Ohana/src/dvomerge/src/LoadImages.c (revision 29938) +++ trunk/Ohana/src/dvomerge/src/LoadImages.c (revision 29938) @@ -0,0 +1,46 @@ +# include "dvomerge.h" + +Image *LoadImages (FITS_DB *db, char *filename, off_t *Nimage) { + + int status; + Image *image; + double ZeroPoint; + + myAssert(filename, "programming errror"); + myAssert(Nimage, "programming errror"); + + gfits_db_init (db); + db->lockstate = LCK_SOFT; + db->timeout = 120.0; + + if (!gfits_db_lock (db, filename)) { + fprintf (stderr, "error opening image catalog %s (1)\n", filename); + return (NULL); + } + + if (db->dbstate == LCK_EMPTY) { + fprintf (stderr, "note: image catalog is empty\n"); + ALLOCATE (image, Image, 1); + *Nimage = 1; + return (image); + } + + status = dvo_image_load (db, TRUE, FALSE); + gfits_db_close (db); + + if (!status) { + fprintf (stderr, "problem loading image database table\n"); + return (NULL); + } + + image = gfits_table_get_Image (&db->ftable, Nimage, &db->swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + return (NULL); + } + + gfits_scan (&db->header, "ZERO_PT", "%lf", 1, &ZeroPoint); + SetZeroPoint (ZeroPoint); + + return (image); +} Index: trunk/Ohana/src/dvomerge/src/ReadDeleteList.c =================================================================== --- trunk/Ohana/src/dvomerge/src/ReadDeleteList.c (revision 29938) +++ trunk/Ohana/src/dvomerge/src/ReadDeleteList.c (revision 29938) @@ -0,0 +1,100 @@ +# 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, *buffer; + + char *indexPoint = NULL; + 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: trunk/Ohana/src/dvomerge/src/args.c =================================================================== --- trunk/Ohana/src/dvomerge/src/args.c (revision 29181) +++ trunk/Ohana/src/dvomerge/src/args.c (revision 29938) @@ -15,4 +15,11 @@ remove_argument (N, argc, argv); ALTERNATE_PHOTCODE_FILE = strdup(argv[N]); + remove_argument (N, argc, argv); + } + + NTHREADS = 0; + if ((N = get_argument (*argc, argv, "-threads"))) { + remove_argument (N, argc, argv); + NTHREADS = MAX(0, atoi(argv[N])); remove_argument (N, argc, argv); } @@ -36,5 +43,5 @@ } - if ((*argc != 6) && (*argc != 4)) dvomerge_usage(); + if ((*argc < 4) || (*argc > 6)) dvomerge_usage(); return TRUE; } Index: trunk/Ohana/src/dvomerge/src/dvo_image_merge_dbs.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvo_image_merge_dbs.c (revision 29181) +++ trunk/Ohana/src/dvomerge/src/dvo_image_merge_dbs.c (revision 29938) @@ -2,4 +2,108 @@ void sort_IDmap (IDmapType *IDmap); + +off_t getTgtIndex (e_time start, e_time stop, short photcode, off_t *TgtIndex, e_time *TgtTimes, short *TgtCodes, off_t NimagesTgt) { + + // use bisection to find the starting entry by time + + off_t Nlo, Nhi, N; + + // find the last TGT before start + Nlo = 0; Nhi = NimagesTgt; + while (Nhi - Nlo > 10) { + N = 0.5*(Nlo + Nhi); + if (TgtTimes[N] < start) { + Nlo = MAX(N, 0); + } else { + Nhi = MIN(N + 1, NimagesTgt); + } + } + + // check for the matched mosaic starting from Nlo + // we may have to go much beyond Nlo since stop is not sorted + // can we use a sorted version of stop to check when we are beyond the valid range?? + for (N = Nlo; N < NimagesTgt; N++) { + if (TgtTimes[N] < start) continue; + if (TgtTimes[N] > stop) return (-1); + if (TgtCodes[N] != photcode) continue; + return (TgtIndex[N]); + } + return (-1); +} + +// sort two times vectors and an index by first time vector +void SortTgtByTimes (e_time *S, off_t *I, short *C, off_t N) { + +# define SWAPFUNC(A,B){ e_time tmp_t; off_t tmp_i; short tmp_c; \ + tmp_t = S[A]; S[A] = S[B]; S[B] = tmp_t; \ + tmp_i = I[A]; I[A] = I[B]; I[B] = tmp_i; \ + tmp_c = C[A]; C[A] = C[B]; C[B] = tmp_c; \ +} +# define COMPARE(A,B)(S[A] < S[B]) + + OHANA_SORT (N, COMPARE, SWAPFUNC); + +# undef SWAPFUNC +# undef COMPARE + +} + +// we have two tables; 'tgt' contains some exposures from 'src' : find them and match a map +int dvo_image_match_dbs (IDmapType *IDmap, FITS_DB *tgt, FITS_DB *src) { + + Image *imagesSrc, *imagesTgt; + off_t NimagesSrc, NimagesTgt; + off_t iSrc, iTgt; + off_t *TgtIndex; + e_time *TgtTimes; + short *TgtCodes; + + imagesSrc = gfits_table_get_Image (&src[0].ftable, &NimagesSrc, &src[0].swapped); + if (!imagesSrc) { + fprintf (stderr, "ERROR: failed to read images from src\n"); + exit (2); + } + + imagesTgt = gfits_table_get_Image (&tgt[0].ftable, &NimagesTgt, &tgt[0].swapped); + if (!imagesTgt) { + fprintf (stderr, "ERROR: failed to read images from tgt\n"); + exit (2); + } + + ALLOCATE (IDmap->old, off_t, NimagesSrc); + ALLOCATE (IDmap->new, off_t, NimagesSrc); + IDmap->Nmap = NimagesSrc; + + // match by time and photcode + ALLOCATE (TgtIndex, off_t, NimagesTgt); + ALLOCATE (TgtTimes, e_time, NimagesTgt); + ALLOCATE (TgtCodes, short, NimagesTgt); + + // save the Index, Times, Codes for all TGT images + for (iTgt = 0; iTgt < NimagesTgt; iTgt++) { + TgtIndex[iTgt] = iTgt; + TgtTimes[iTgt] = imagesTgt[iTgt].tzero; + TgtCodes[iTgt] = imagesTgt[iTgt].photcode; + } + + // sort the index, start, and stop by the start times: + SortTgtByTimes (TgtTimes, TgtIndex, TgtCodes, NimagesTgt); + + for (iSrc = 0; iSrc < NimagesSrc; iSrc++) { + iTgt = getTgtIndex (imagesSrc[iSrc].tzero, imagesSrc[iSrc].tzero + (int) imagesSrc[iSrc].exptime, imagesSrc[iSrc].photcode, TgtIndex, TgtTimes, TgtCodes, NimagesTgt); + if (iTgt < 0) Shutdown ("failure to match images: %s\n", imagesSrc[iSrc].name); + IDmap[0].old[iSrc] = imagesSrc[iSrc].imageID; + IDmap[0].new[iSrc] = imagesTgt[iTgt].imageID; + } + + FREE(TgtIndex); + FREE(TgtTimes); + FREE(TgtCodes); + + // sort IDmap->old,new on the basis of IDmap->old: + sort_IDmap (IDmap); + + return TRUE; +} // merge db2 into db1 Index: trunk/Ohana/src/dvomerge/src/dvoconvert.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvoconvert.c (revision 29181) +++ trunk/Ohana/src/dvomerge/src/dvoconvert.c (revision 29938) @@ -6,5 +6,5 @@ char filename[256], *input, *output; - int depth; + int depth, Nsecfilt; off_t i, j, Ns, Ne; SkyTable *outsky, *insky; @@ -32,4 +32,6 @@ exit (1); } + Nsecfilt = GetPhotcodeNsecfilt(); + // save the photcodes in the output catdir sprintf (filename, "%s/Photcodes.dat", output); @@ -67,5 +69,5 @@ // load / create output catalog - LoadCatalog (&outcatalog, &outsky[0].regions[i], outsky[0].filename[i], "w"); + LoadCatalog (&outcatalog, &outsky[0].regions[i], outsky[0].filename[i], "w", Nsecfilt); // combine only tables at equal or larger depth @@ -77,5 +79,5 @@ // load input catalog - LoadCatalog (&incatalog, inlist[0].regions[j], inlist[0].filename[j], "r"); + LoadCatalog (&incatalog, inlist[0].regions[j], inlist[0].filename[j], "r", Nsecfilt); // skip empty input catalogs Index: trunk/Ohana/src/dvomerge/src/dvomerge.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvomerge.c (revision 29181) +++ trunk/Ohana/src/dvomerge/src/dvomerge.c (revision 29938) @@ -8,8 +8,18 @@ dvomerge_args (&argc, argv); - // XXX require both inputs to be sorted? - - if (argc == 6) dvomergeCreate (argc, argv); - if (argc == 4) dvomergeUpdate (argc, argv); + if (argc == 6) { + if (!strcasecmp (argv[2], "and")) { + dvomergeCreate (argc, argv); + } else { + dvomergeFromList (argc, argv); + } + } + if ((argc == 4) || (argc == 5)) { + if (NTHREADS) { + dvomergeUpdate_threaded (argc, argv); + } else { + dvomergeUpdate (argc, argv); + } + } dvomerge_usage(); exit (2); // cannot reach here. @@ -18,6 +28,7 @@ /* we have two possible modes of operation: - Create : dvomerge (in1) and (in2) to (out) -- create a new db from two input dbs - Update : dvomerge (in) into (out) -- merge a new db into an existing db + Create : dvomerge (in1) and (in2) to (out) -- create a new db from two input dbs + Update : dvomerge (in) into (out) -- merge a new db into an existing db + Continue : dvomerge (in) into (out) continue -- merge a new db into an existing db */ Index: trunk/Ohana/src/dvomerge/src/dvomergeContinue.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvomergeContinue.c (revision 29938) +++ trunk/Ohana/src/dvomerge/src/dvomergeContinue.c (revision 29938) @@ -0,0 +1,164 @@ +# include "dvomerge.h" + +int dvomergeContinue (int argc, char **argv) { + + int depth; + off_t i, j, Ns, Ne; + SkyTable *outsky, *insky; + SkyList *outlist, *inlist; + Catalog incatalog, outcatalog; + char filename[256], *input, *output; + IDmapType IDmap; + PhotCodeData *inputPhotcodes; + PhotCodeData *outputPhotcodes; + int *secfiltMap = NULL; + int NsecfiltInput, NsecfiltOutput; + + double dtime; + struct timeval start, stop; + gettimeofday (&start, NULL); + + if (strcasecmp (argv[2], "into")) dvomerge_usage(); + if (strcasecmp (argv[4], "continue")) dvomerge_usage(); + + input = argv[1]; + output = argv[3]; + + if (ALTERNATE_PHOTCODE_FILE) { + fprintf (stderr, "cannot specify photcodes when merging into an existing catdir\n"); + exit (1); + } + + SetPhotcodeTable(NULL); + sprintf (filename, "%s/Photcodes.dat", input); + if (!LoadPhotcodes (filename, NULL, FALSE)) { + fprintf (stderr, "error reading input database directory %s\n", input); + exit (1); + } + inputPhotcodes = GetPhotcodeTable(); + NsecfiltInput = GetPhotcodeNsecfilt(); + + // since we are merging the input db into the output db, the output defines the photcode + // table & db layout but, this requires the output to exist. if it does not, instead use the + // input. + + SetPhotcodeTable(NULL); + sprintf (filename, "%s/Photcodes.dat", output); + if (LoadPhotcodes (filename, NULL, FALSE)) { + + outputPhotcodes = GetPhotcodeTable(); + NsecfiltOutput = GetPhotcodeNsecfilt(); + + secfiltMap = GetSecFiltMap(outputPhotcodes, inputPhotcodes); + if (!secfiltMap) { + fprintf (stderr, "failed to map input secfilt photcodes to output photcodes table\n"); + exit (1); + } + } else { + fprintf(stderr, "%s not found using photcodes from %s/Photcodes.dat\n", filename, input); + outputPhotcodes = inputPhotcodes; + NsecfiltOutput = NsecfiltInput; + + if (!check_dir_access (output, VERBOSE)) { + fprintf (stderr, "error creating output database directory %s\n", output); + exit (1); + } + SetPhotcodeTable(inputPhotcodes); + if (!SavePhotcodesFITS (filename)) { + fprintf (stderr, "error saving photcode table in %s/Photcodes.dat\n", output); + exit (1); + } + } + + // need to determine the mapping from the input to the output images + dvomergeImagesGetMap (&IDmap, input, output); + + // load the sky table for the existing database + insky = SkyTableLoadOptimal (input, NULL, NULL, FALSE, SKY_DEPTH_HST, VERBOSE); + if (!insky) { + Shutdown ("can't read SkyTable for %s", input); + } + SkyTableSetFilenames (insky, input, "cpt"); + + // XXX apply this...generate the subset matching the user-selected region + inlist = SkyListByPatch (insky, -1, &UserPatch); + + // generate an output table populated at the desired depth + outsky = SkyTableLoadOptimal (output, NULL, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); + if (!outsky) { + Shutdown ("can't read or create SkyTable for %s", output); + } + SkyTableSetFilenames (outsky, output, "cpt"); + + // loop over the populatable output tables; check for data in input in the corresponding regions + + SkyListPopulatedRange (&Ns, &Ne, inlist, 0); + depth = inlist[0].regions[Ns][0].depth; + + SetPhotcodeTable(NULL); + + // loop over the populated input regions + for (i = 0; i < inlist[0].Nregions; i++) { + if (!inlist[0].regions[i][0].table) continue; + if (VERBOSE) fprintf (stderr, "input: %s\n", inlist[0].regions[i][0].name); + + // load / create output catalog (if catalog does not exist, it will be created) + LoadCatalog (&incatalog, &inlist[0].regions[i][0], inlist[0].filename[i], "r", NsecfiltInput); + // skip empty input catalogs + if (!incatalog.Naves_disk) { + dvo_catalog_unlock (&incatalog); + dvo_catalog_free (&incatalog); + continue; + } + + // combine only tables at equal or larger depth + + // load in all of the tables from input for this region + // SkyListByBounds will return neighbor catalogs if the boundaries exactly match (due to rounding). Since the regions are not infinitely small, + // compare to a slightly reduced footprint + float dPos = 2.0/3600.0; + outlist = SkyListByBounds (outsky, depth, inlist[0].regions[i][0].Rmin + dPos, inlist[0].regions[i][0].Rmax - dPos, inlist[0].regions[i][0].Dmin + dPos, inlist[0].regions[i][0].Dmax - dPos); + for (j = 0; j < outlist[0].Nregions; j++) { + if (VERBOSE) fprintf (stderr, "output : %s\n", outlist[0].regions[j][0].name); + + // load input catalog + LoadCatalog (&outcatalog, outlist[0].regions[j], outlist[0].filename[j], "w", NsecfiltOutput); + + dvo_update_image_IDs (&IDmap, &incatalog); + merge_catalogs_old (&outsky[0].regions[j], &outcatalog, &incatalog, RADIUS, secfiltMap); + + outcatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; + + // if we receive a signal which would cause us to exit, wait until the full catalog is written + SetProtect (TRUE); + if (!dvo_catalog_save (&outcatalog, VERBOSE)) { + fprintf (stderr, "ERROR: failed to save catalog %s\n", outlist[0].filename[j]); + exit (1); + } + SetProtect (FALSE); + + dvo_catalog_unlock (&outcatalog); + dvo_catalog_free (&outcatalog); + + fprintf (stderr, "merged %s into %s\n", inlist[0].regions[i][0].name, outlist[0].regions[j][0].name); + } + SkyListFree (outlist); + + dvo_catalog_unlock (&incatalog); + dvo_catalog_free (&incatalog); + } + + // save the output sky table copy + sprintf (filename, "%s/SkyTable.fits", output); + check_file_access (filename, TRUE, TRUE, VERBOSE); + if (!SkyTableSave (outsky, filename)) { + fprintf (stderr, "ERROR: failed to save sky table for %s\n", output); + exit (1); + } + + gettimeofday (&stop, NULL); + dtime = DTIME (stop, start); + fprintf (stderr, "SUCCESS: elapsed time %9.4f sec\n", dtime); + + exit (0); +} Index: trunk/Ohana/src/dvomerge/src/dvomergeContinue_threaded.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvomergeContinue_threaded.c (revision 29938) +++ trunk/Ohana/src/dvomerge/src/dvomergeContinue_threaded.c (revision 29938) @@ -0,0 +1,313 @@ +# include "dvomerge.h" +# include + +enum { + TS_WAIT, + TS_RUN, + TS_DONE, + TS_FAIL, + TS_EXIT, +}; + +typedef struct { + int iThread; + int state; + // elements used to carry data to the threadjob + char *name; + char *filename; + SkyRegion *region; + int NsecfiltInput; + int NsecfiltOutput; + IDmapType *IDmap; + SkyTable *outsky; + int *secfiltMap; + int depth; +} ThreadData; + +void *dvomergeUpdate_threadjob (void *inputData) { + + off_t j; + Catalog incatalog, outcatalog; + SkyList *outlist; + + ThreadData *threadData = (ThreadData *) inputData; + + while (TRUE) { + + while (threadData->state != TS_RUN) { + usleep (100); + } + + if (VERBOSE) fprintf (stderr, "input: %s\n", threadData->name); + + // load / create output catalog (if catalog does not exist, it will be created) + LoadCatalog (&incatalog, threadData->region, threadData->filename, "r", threadData->NsecfiltInput); + // skip empty input catalogs + if (!incatalog.Naves_disk) { + dvo_catalog_unlock (&incatalog); + dvo_catalog_free (&incatalog); + threadData->state = TS_DONE; + continue; + } + + // combine only tables at equal or larger depth. + + // NOTE: this requirement is especially important for the threaded version: if the + // output catalog is guaranteed to be the same size or finer than the input, then two + // input catalogs cannot collide on their target output catalog. + + // load in all of the tables from input for this region + // SkyListByBounds will return neighbor catalogs if the boundaries exactly match (due to rounding). Since the regions are not infinitely small, + // compare to a slightly reduced footprint + float dPos = 2.0/3600.0; + outlist = SkyListByBounds (threadData->outsky, threadData->depth, threadData->region->Rmin + dPos, threadData->region->Rmax - dPos, threadData->region->Dmin + dPos, threadData->region->Dmax - dPos); + for (j = 0; (j < outlist[0].Nregions) && (threadData->state != TS_FAIL); j++) { + if (VERBOSE) fprintf (stderr, "output : %s\n", outlist[0].regions[j][0].name); + + // load input catalog + LoadCatalog (&outcatalog, outlist[0].regions[j], outlist[0].filename[j], "w", threadData->NsecfiltOutput); + + dvo_update_image_IDs (threadData->IDmap, &incatalog); + merge_catalogs_old (&threadData->outsky->regions[j], &outcatalog, &incatalog, RADIUS, threadData->secfiltMap); + + outcatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; + + // if we receive a signal which would cause us to exit, wait until the full catalog is written + SetProtect (TRUE); + if (!dvo_catalog_save (&outcatalog, VERBOSE)) { + fprintf (stderr, "ERROR: failed to save catalog %s\n", outlist[0].filename[j]); + threadData->state = TS_FAIL; + continue; + } + SetProtect (FALSE); + + dvo_catalog_unlock (&outcatalog); + dvo_catalog_free (&outcatalog); + + fprintf (stderr, "merged %s into %s\n", threadData->name, outlist[0].regions[j][0].name); + } + SkyListFree (outlist); + + dvo_catalog_unlock (&incatalog); + dvo_catalog_free (&incatalog); + threadData->state = TS_DONE; + } + return NULL; +} + +int dvomergeUpdate_threaded (int argc, char **argv) { + + int j, launched, done, Nwait, Nrun; + + ThreadData *threadData; + pthread_t *threads; + + int depth; + off_t i, Ns, Ne; + SkyTable *outsky, *insky; + SkyList *inlist; + char filename[256], *input, *output; + IDmapType IDmap; + PhotCodeData *inputPhotcodes; + PhotCodeData *outputPhotcodes; + int *secfiltMap = NULL; + int NsecfiltInput, NsecfiltOutput; + + double dtime; + struct timeval start, stop; + gettimeofday (&start, NULL); + + CONTINUE = FALSE; + if (strcasecmp (argv[2], "into")) dvomerge_usage(); + if (argc == 5) { + if (strcasecmp (argv[4], "continue")) dvomerge_usage(); + CONTINUE = TRUE; + } + + input = argv[1]; + output = argv[3]; + + if (ALTERNATE_PHOTCODE_FILE) { + fprintf (stderr, "cannot specify photcodes when merging into an existing catdir\n"); + exit (1); + } + + SetPhotcodeTable(NULL); + sprintf (filename, "%s/Photcodes.dat", input); + if (!LoadPhotcodes (filename, NULL, FALSE)) { + fprintf (stderr, "error reading input database directory %s\n", input); + exit (1); + } + inputPhotcodes = GetPhotcodeTable(); + NsecfiltInput = GetPhotcodeNsecfilt(); + + // since we are merging the input db into the output db, the output defines the photcode + // table & db layout but, this requires the output to exist. if it does not, instead use the + // input. + + SetPhotcodeTable(NULL); + sprintf (filename, "%s/Photcodes.dat", output); + if (LoadPhotcodes (filename, NULL, FALSE)) { + + outputPhotcodes = GetPhotcodeTable(); + NsecfiltOutput = GetPhotcodeNsecfilt(); + + secfiltMap = GetSecFiltMap(outputPhotcodes, inputPhotcodes); + if (!secfiltMap) { + fprintf (stderr, "failed to map input secfilt photcodes to output photcodes table\n"); + exit (1); + } + } else { + fprintf(stderr, "%s not found using photcodes from %s/Photcodes.dat\n", filename, input); + outputPhotcodes = inputPhotcodes; + NsecfiltOutput = NsecfiltInput; + + if (!check_dir_access (output, VERBOSE)) { + fprintf (stderr, "error creating output database directory %s\n", output); + exit (1); + } + SetPhotcodeTable(inputPhotcodes); + if (!SavePhotcodesFITS (filename)) { + fprintf (stderr, "error saving photcode table in %s/Photcodes.dat\n", output); + exit (1); + } + } + + if (CONTINUE) { + // need to determine the mapping from the input to the output images + dvomergeImagesGetMap (&IDmap, input, output); + } else { + dvomergeImagesUpdate (&IDmap, input, output); + } + + // load the sky table for the existing database + insky = SkyTableLoadOptimal (input, NULL, NULL, FALSE, SKY_DEPTH_HST, VERBOSE); + if (!insky) { + Shutdown ("can't read SkyTable for %s", input); + } + SkyTableSetFilenames (insky, input, "cpt"); + + // XXX apply this...generate the subset matching the user-selected region + inlist = SkyListByPatch (insky, -1, &UserPatch); + + // generate an output table populated at the desired depth + outsky = SkyTableLoadOptimal (output, NULL, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); + if (!outsky) { + Shutdown ("can't read or create SkyTable for %s", output); + } + SkyTableSetFilenames (outsky, output, "cpt"); + + // loop over the populatable output tables; check for data in input in the corresponding regions + + SkyListPopulatedRange (&Ns, &Ne, inlist, 0); + depth = inlist[0].regions[Ns][0].depth; + + SetPhotcodeTable(NULL); + + ALLOCATE(threads, pthread_t, NTHREADS); + ALLOCATE(threadData, ThreadData, NTHREADS); + for (i = 0; i < NTHREADS; i++) { + threadData[i].state = TS_WAIT; + threadData[i].iThread = i; + } + for (i = 0; i < NTHREADS; i++) { + pthread_create (&threads[i], NULL, &dvomergeUpdate_threadjob, (void *) &threadData[i]); + } + + // loop over the populated input regions + for (i = 0; i < inlist[0].Nregions; i++) { + if (!inlist[0].regions[i][0].table) continue; + + launched = FALSE; + while (!launched) { + + // are any threads done? can we harvest them? + for (j = 0; j < NTHREADS; j++) { + if (threadData[j].state == TS_FAIL) goto failure; + if (threadData[j].state != TS_DONE) continue; + // fprintf (stderr, "harvested thread %d\n", threadData[j].iThread); + threadData[j].state = TS_WAIT; + } + + // are any threads available? can we launch this job? + for (j = 0; !launched && (j < NTHREADS); j++) { + if (threadData[j].state == TS_FAIL) goto failure; + if (threadData[j].state != TS_WAIT) continue; + + threadData[j].name = inlist[0].regions[i][0].name; + threadData[j].filename = inlist[0].filename[i]; + threadData[j].region = inlist[0].regions[i]; + threadData[j].NsecfiltInput = NsecfiltInput; + threadData[j].NsecfiltOutput = NsecfiltOutput; + threadData[j].IDmap = &IDmap; + threadData[j].outsky = outsky; + threadData[j].secfiltMap = secfiltMap; + threadData[j].depth = depth; + threadData[j].state = TS_RUN; + launched = TRUE; + } + if (!launched) usleep (100); + } + } + + // wait for all threads to be done + done = FALSE; + while (!done) { + // are any threads done? can we harvest them? + for (j = 0; j < NTHREADS; j++) { + if (threadData[j].state == TS_FAIL) goto failure; + if (threadData[j].state != TS_DONE) continue; + threadData[j].state = TS_WAIT; + } + + // how many are waiting now? + Nwait = 0; + for (j = 0; j < NTHREADS; j++) { + if (threadData[j].state == TS_FAIL) goto failure; + if (threadData[j].state != TS_WAIT) continue; + Nwait ++; + } + + if (Nwait < NTHREADS) { + usleep (100); + } else { + done = TRUE; + } + } + + // save the output sky table copy + sprintf (filename, "%s/SkyTable.fits", output); + check_file_access (filename, TRUE, TRUE, VERBOSE); + if (!SkyTableSave (outsky, filename)) { + fprintf (stderr, "ERROR: failed to save sky table for %s\n", output); + exit (1); + } + + gettimeofday (&stop, NULL); + dtime = DTIME (stop, start); + fprintf (stderr, "SUCCESS: elapsed time %9.4f sec\n", dtime); + + exit (0); + +failure: + + // wait for all threads to be done + done = FALSE; + while (!done) { + + // how many are still running? + Nrun = 0; + for (j = 0; j < NTHREADS; j++) { + if (threadData[j].state != TS_RUN) continue; + Nrun ++; + } + + if (Nrun) { + usleep (100); + } else { + done = TRUE; + } + } + exit (1); +} + Index: trunk/Ohana/src/dvomerge/src/dvomergeCreate.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvomergeCreate.c (revision 29181) +++ trunk/Ohana/src/dvomerge/src/dvomergeCreate.c (revision 29938) @@ -15,4 +15,5 @@ char *input1, *input2, *output; IDmapType IDmap1, IDmap2; + int NsecfiltInput1, NsecfiltInput2, NsecfiltOutput; if (strcasecmp (argv[2], "and")) dvomerge_usage(); @@ -39,4 +40,5 @@ } input1Photcodes = GetPhotcodeTable(); + NsecfiltInput1 = GetPhotcodeNsecfilt(); // Read the input2 photcodes @@ -48,4 +50,5 @@ } input2Photcodes = GetPhotcodeTable(); + NsecfiltInput2 = GetPhotcodeNsecfilt(); if (ALTERNATE_PHOTCODE_FILE) { @@ -60,4 +63,5 @@ } outputPhotcodes = GetPhotcodeTable(); + NsecfiltOutput = GetPhotcodeNsecfilt(); secfiltMap1 = GetSecFiltMap(outputPhotcodes, input1Photcodes); @@ -74,4 +78,5 @@ outputPhotcodes = input1Photcodes; codesFilename = filename1; + NsecfiltOutput = NsecfiltInput1; // secfitMap1 can remain NULL since input1 defines the set of codes } @@ -84,6 +89,9 @@ } + // trigger any dependencies, if any + SetPhotcodeTable(NULL); + // save the output photcodes in the output catdir - SetPhotcodeTable(outputPhotcodes); + // SetPhotcodeTable(outputPhotcodes); sprintf (filename, "%s/Photcodes.dat", output); if (!check_file_access (filename, TRUE, TRUE, VERBOSE)) { @@ -121,13 +129,13 @@ if (VERBOSE) fprintf (stderr, "output: %s\n", outsky[0].regions[i].name); - if (outputPhotcodes) { - SetPhotcodeTable(outputPhotcodes); - } + // if (outputPhotcodes) { + // SetPhotcodeTable(outputPhotcodes); + // } // load / create output catalog - LoadCatalog (&outcatalog, &outsky[0].regions[i], outsky[0].filename[i], "w"); - - if (input1Photcodes) { - SetPhotcodeTable(input1Photcodes); - } + LoadCatalog (&outcatalog, &outsky[0].regions[i], outsky[0].filename[i], "w", NsecfiltOutput); + + // if (input1Photcodes) { + // SetPhotcodeTable(input1Photcodes); + // } // combine only tables at equal or larger depth @@ -138,5 +146,5 @@ // load input catalog (1) - LoadCatalog (&incatalog, inlist[0].regions[j], inlist[0].filename[j], "r"); + LoadCatalog (&incatalog, inlist[0].regions[j], inlist[0].filename[j], "r", NsecfiltInput1); // skip empty input catalogs @@ -153,7 +161,7 @@ SkyListFree (inlist); - if (input2Photcodes) { - SetPhotcodeTable(input2Photcodes); - } + // if (input2Photcodes) { + // SetPhotcodeTable(input2Photcodes); + // } // load in all of the tables from input2 for this region @@ -163,5 +171,5 @@ // load input catalog (2) - LoadCatalog (&incatalog, inlist[0].regions[j], inlist[0].filename[j], "r"); + LoadCatalog (&incatalog, inlist[0].regions[j], inlist[0].filename[j], "r", NsecfiltInput2); // skip empty input catalogs Index: trunk/Ohana/src/dvomerge/src/dvomergeFromList.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvomergeFromList.c (revision 29938) +++ trunk/Ohana/src/dvomerge/src/dvomergeFromList.c (revision 29938) @@ -0,0 +1,158 @@ +# include "dvomerge.h" + +int dvomergeFromList (int argc, char **argv) { + + off_t i, j; + Catalog incatalog, outcatalog; + char filename[256], *input, *output; + IDmapType IDmap; + PhotCodeData *inputPhotcodes; + PhotCodeData *outputPhotcodes; + int *secfiltMap = NULL; + char *listname, **list, inputfile[256], outputfile[256]; + int Nlist, NLIST; + int NsecfiltInput, NsecfiltOutput; + + if (strcasecmp (argv[2], "into")) dvomerge_usage(); + if (strcasecmp (argv[4], "from")) dvomerge_usage(); + + input = argv[1]; + output = argv[3]; + listname = argv[5]; + + Nlist = 0; + NLIST = 100; + ALLOCATE(list, char *, NLIST); + for (i = 0; i < NLIST; i++) { + ALLOCATE(list[i], char, 256); + memset(list[i], 0, 256); + } + + FILE *listfile = fopen(listname, "r"); + myAssert(listfile, "error opening list"); + + for (i = 0; fscanf (listfile, "%s", list[i]) != EOF; i++) { + if (i == NLIST - 1) { + NLIST += 100; + REALLOCATE(list, char *, NLIST); + for (j = i + 1; j < NLIST; j++) { + ALLOCATE(list[j], char, 256); + memset(list[j], 0, 256); + } + } + } + Nlist = i; + fclose(listfile); + + if (ALTERNATE_PHOTCODE_FILE) { + fprintf (stderr, "cannot specify photcodes when merging into an existing catdir\n"); + exit (1); + } + + SetPhotcodeTable(NULL); + sprintf (filename, "%s/Photcodes.dat", input); + if (!LoadPhotcodes (filename, NULL, FALSE)) { + fprintf (stderr, "error reading input database directory %s\n", input); + exit (1); + } + inputPhotcodes = GetPhotcodeTable(); + NsecfiltInput = GetPhotcodeNsecfilt(); + + // since we are merging the input db into the output db, the output defines the photcode + // table & db layout but, this requires the output to exist. if it does not, instead use the + // input. + + SetPhotcodeTable(NULL); + sprintf (filename, "%s/Photcodes.dat", output); + if (LoadPhotcodes (filename, NULL, FALSE)) { + + outputPhotcodes = GetPhotcodeTable(); + NsecfiltOutput = GetPhotcodeNsecfilt(); + + secfiltMap = GetSecFiltMap(outputPhotcodes, inputPhotcodes); + if (!secfiltMap) { + fprintf (stderr, "failed to map input secfilt photcodes to output photcodes table\n"); + exit (1); + } + } else { + fprintf(stderr, "%s not found using photcodes from %s/Photcodes.dat\n", filename, input); + outputPhotcodes = inputPhotcodes; + NsecfiltOutput = NsecfiltInput; + + if (!check_dir_access (output, VERBOSE)) { + fprintf (stderr, "error creating output database directory %s\n", output); + exit (1); + } + SetPhotcodeTable(inputPhotcodes); + if (!SavePhotcodesFITS (filename)) { + fprintf (stderr, "error saving photcode table in %s/Photcodes.dat\n", output); + exit (1); + } + } + + // XXX need to determine the mapping from the input to the output images + dvomergeImagesGetMap (&IDmap, input, output); + + SetPhotcodeTable(NULL); + + // loop over the populated input regions + for (i = 0; i < Nlist; i++) { + if (VERBOSE) fprintf (stderr, "input: %s\n", list[i]); + + sprintf (inputfile, "%s/%s", input, list[i]); + sprintf (outputfile, "%s/%s", output, list[i]); + + LoadCatalog (&incatalog, NULL, inputfile, "r", NsecfiltInput); + + // skip empty input catalogs + if (!incatalog.Naves_disk) { + dvo_catalog_unlock (&incatalog); + dvo_catalog_free (&incatalog); + continue; + } + + // combine only tables at equal depth + + SkyRegion skyregion; + gfits_scan (&incatalog.header, "RA0", "%f", 1, &skyregion.Rmin); + gfits_scan (&incatalog.header, "RA1", "%f", 1, &skyregion.Rmax); + gfits_scan (&incatalog.header, "DEC0", "%f", 1, &skyregion.Dmin); + gfits_scan (&incatalog.header, "DEC1", "%f", 1, &skyregion.Dmax); + + // load input catalog + // SetPhotcodeTable(outputPhotcodes); + LoadCatalog (&outcatalog, NULL, outputfile, "w", NsecfiltOutput); + + if (outcatalog.Naverage == 0) { + gfits_modify (&outcatalog.header, "RA0", "%f", 1, skyregion.Rmin); + gfits_modify (&outcatalog.header, "DEC0", "%f", 1, skyregion.Dmin); + gfits_modify (&outcatalog.header, "RA1", "%f", 1, skyregion.Rmax); + gfits_modify (&outcatalog.header, "DEC1", "%f", 1, skyregion.Dmax); + gfits_modify (&outcatalog.header, "CATID", "%d", 1, incatalog.catID); + outcatalog.catID = incatalog.catID; + } + + dvo_update_image_IDs (&IDmap, &incatalog); + merge_catalogs_old (&skyregion, &outcatalog, &incatalog, RADIUS, secfiltMap); + + outcatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; + + // if we receive a signal which would cause us to exit, wait until the full catalog is written + SetProtect (TRUE); + if (!dvo_catalog_save (&outcatalog, VERBOSE)) { + fprintf (stderr, "ERROR: failed to save catalog %s\n", outputfile); + exit (1); + } + SetProtect (FALSE); + + fprintf (stderr, "merged %s into %s\n", inputfile, outputfile); + + dvo_catalog_unlock (&outcatalog); + dvo_catalog_free (&outcatalog); + + dvo_catalog_unlock (&incatalog); + dvo_catalog_free (&incatalog); + } + + exit (0); +} Index: trunk/Ohana/src/dvomerge/src/dvomergeImageIDs.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvomergeImageIDs.c (revision 29938) +++ trunk/Ohana/src/dvomerge/src/dvomergeImageIDs.c (revision 29938) @@ -0,0 +1,100 @@ +# include "dvomerge.h" + +/*** update the image table ***/ +int dvomergeImagesUpdate (IDmapType *IDmap, char *input, char *output) { + + FITS_DB inDB; + FITS_DB outDB; + int status; + + /*** load input1/Images.dat ***/ + sprintf (ImageCat, "%s/Images.dat", input); + inDB.mode = dvo_catalog_catmode (CATMODE); + inDB.format = dvo_catalog_catformat (CATFORMAT); + status = dvo_image_lock (&inDB, ImageCat, 3600.0, LCK_XCLD); // shorter timeout? + if (!status) Shutdown ("ERROR: failure to lock image catalog %s", inDB.filename); + + // load the image table + if (inDB.dbstate == LCK_EMPTY) { + IDmap->Nmap = 0; + return TRUE; + } + if (!dvo_image_load (&inDB, VERBOSE, TRUE)) { + Shutdown ("can't read input image catalog %s", inDB.filename); + } + + /*** load output/Images.dat ***/ + sprintf (ImageCat, "%s/Images.dat", output); + outDB.mode = dvo_catalog_catmode (CATMODE); + outDB.format = dvo_catalog_catformat (CATFORMAT); + status = dvo_image_lock (&outDB, ImageCat, 3600.0, LCK_XCLD); // shorter timeout? + if (!status) Shutdown ("ERROR: failure to lock image catalog %s", outDB.filename); + + /* load the image table */ + if (outDB.dbstate == LCK_EMPTY) { + dvo_image_create (&outDB, GetZeroPoint()); + } else { + if (!dvo_image_load (&outDB, VERBOSE, TRUE)) { + Shutdown ("can't read output image catalog %s", outDB.filename); + } + } + + // convert database table to internal structure & add to output image db + dvo_image_merge_dbs(IDmap, &outDB, &inDB); + dvo_image_unlock (&inDB); // unlock input + + SetProtect (TRUE); + dvo_image_save (&outDB, VERBOSE); + SetProtect (FALSE); + dvo_image_unlock (&outDB); // unlock output + + return TRUE; +} + +/*** update the image table ***/ +int dvomergeImagesGetMap (IDmapType *IDmap, char *input, char *output) { + + FITS_DB inDB; + FITS_DB outDB; + int status; + + /*** load input1/Images.dat ***/ + sprintf (ImageCat, "%s/Images.dat", input); + inDB.mode = dvo_catalog_catmode (CATMODE); + inDB.format = dvo_catalog_catformat (CATFORMAT); + status = dvo_image_lock (&inDB, ImageCat, 3600.0, LCK_XCLD); // shorter timeout? + if (!status) Shutdown ("ERROR: failure to lock image catalog %s", inDB.filename); + + // load the image table + if (inDB.dbstate == LCK_EMPTY) { + IDmap->Nmap = 0; + return TRUE; + } + + // the input database is allowed to have no images + if (!dvo_image_load (&inDB, VERBOSE, TRUE)) { + Shutdown ("can't read input image catalog %s", inDB.filename); + } + + /*** load output/Images.dat ***/ + sprintf (ImageCat, "%s/Images.dat", output); + outDB.mode = dvo_catalog_catmode (CATMODE); + outDB.format = dvo_catalog_catformat (CATFORMAT); + status = dvo_image_lock (&outDB, ImageCat, 3600.0, LCK_XCLD); // shorter timeout? + if (!status) Shutdown ("ERROR: failure to lock image catalog %s", outDB.filename); + + /* load the image table */ + if (outDB.dbstate == LCK_EMPTY) { + Shutdown ("only use -continue for an existing, partially merged database"); + } + if (!dvo_image_load (&outDB, VERBOSE, TRUE)) { + Shutdown ("can't read output image catalog %s", outDB.filename); + } + + // convert database table to internal structure & add to output image db + dvo_image_match_dbs(IDmap, &outDB, &inDB); + dvo_image_unlock (&inDB); // unlock input + dvo_image_unlock (&outDB); // unlock output + + return TRUE; +} Index: trunk/Ohana/src/dvomerge/src/dvomergeUpdate.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvomergeUpdate.c (revision 29181) +++ trunk/Ohana/src/dvomerge/src/dvomergeUpdate.c (revision 29938) @@ -3,5 +3,5 @@ int dvomergeUpdate (int argc, char **argv) { - int depth; + int depth, CONTINUE; off_t i, j, Ns, Ne; SkyTable *outsky, *insky; @@ -13,6 +13,16 @@ PhotCodeData *outputPhotcodes; int *secfiltMap = NULL; + int NsecfiltInput, NsecfiltOutput; + double dtime; + struct timeval start, stop; + gettimeofday (&start, NULL); + + CONTINUE = FALSE; if (strcasecmp (argv[2], "into")) dvomerge_usage(); + if (argc == 5) { + if (strcasecmp (argv[4], "continue")) dvomerge_usage(); + CONTINUE = TRUE; + } input = argv[1]; @@ -31,4 +41,5 @@ } inputPhotcodes = GetPhotcodeTable(); + NsecfiltInput = GetPhotcodeNsecfilt(); // since we are merging the input db into the output db, the output defines the photcode @@ -41,4 +52,5 @@ outputPhotcodes = GetPhotcodeTable(); + NsecfiltOutput = GetPhotcodeNsecfilt(); secfiltMap = GetSecFiltMap(outputPhotcodes, inputPhotcodes); @@ -50,4 +62,6 @@ fprintf(stderr, "%s not found using photcodes from %s/Photcodes.dat\n", filename, input); outputPhotcodes = inputPhotcodes; + NsecfiltOutput = NsecfiltInput; + if (!check_dir_access (output, VERBOSE)) { fprintf (stderr, "error creating output database directory %s\n", output); @@ -61,5 +75,10 @@ } - dvomergeImagesUpdate (&IDmap, input, output); + if (CONTINUE) { + // need to determine the mapping from the input to the output images + dvomergeImagesGetMap (&IDmap, input, output); + } else { + dvomergeImagesUpdate (&IDmap, input, output); + } // load the sky table for the existing database @@ -85,4 +104,6 @@ depth = inlist[0].regions[Ns][0].depth; + SetPhotcodeTable(NULL); + // loop over the populated input regions for (i = 0; i < inlist[0].Nregions; i++) { @@ -90,7 +111,6 @@ if (VERBOSE) fprintf (stderr, "input: %s\n", inlist[0].regions[i][0].name); - SetPhotcodeTable(inputPhotcodes); // load / create output catalog (if catalog does not exist, it will be created) - LoadCatalog (&incatalog, &inlist[0].regions[i][0], inlist[0].filename[i], "r"); + LoadCatalog (&incatalog, &inlist[0].regions[i][0], inlist[0].filename[i], "r", NsecfiltInput); // skip empty input catalogs if (!incatalog.Naves_disk) { @@ -111,6 +131,5 @@ // load input catalog - SetPhotcodeTable(outputPhotcodes); - LoadCatalog (&outcatalog, outlist[0].regions[j], outlist[0].filename[j], "w"); + LoadCatalog (&outcatalog, outlist[0].regions[j], outlist[0].filename[j], "w", NsecfiltOutput); dvo_update_image_IDs (&IDmap, &incatalog); @@ -121,5 +140,8 @@ // if we receive a signal which would cause us to exit, wait until the full catalog is written SetProtect (TRUE); - dvo_catalog_save (&outcatalog, VERBOSE); + if (!dvo_catalog_save (&outcatalog, VERBOSE)) { + fprintf (stderr, "ERROR: failed to save catalog %s\n", outlist[0].filename[j]); + exit (1); + } SetProtect (FALSE); @@ -143,56 +165,8 @@ } + gettimeofday (&stop, NULL); + dtime = DTIME (stop, start); + fprintf (stderr, "SUCCESS: elapsed time %9.4f sec\n", dtime); + exit (0); } - -/*** update the image table ***/ -int dvomergeImagesUpdate (IDmapType *IDmap, char *input, char *output) { - - FITS_DB inDB; - FITS_DB outDB; - int status; - - /*** load input1/Images.dat ***/ - sprintf (ImageCat, "%s/Images.dat", input); - inDB.mode = dvo_catalog_catmode (CATMODE); - inDB.format = dvo_catalog_catformat (CATFORMAT); - status = dvo_image_lock (&inDB, ImageCat, 3600.0, LCK_XCLD); // shorter timeout? - if (!status) Shutdown ("ERROR: failure to lock image catalog %s", inDB.filename); - - // load the image table - if (inDB.dbstate == LCK_EMPTY) { - // Shutdown ("can't find input image catalog %s", inDB.filename); - IDmap->Nmap = 0; - return TRUE; - } - if (!dvo_image_load (&inDB, VERBOSE, TRUE)) { - Shutdown ("can't read input image catalog %s", inDB.filename); - } - - /*** load output/Images.dat ***/ - sprintf (ImageCat, "%s/Images.dat", output); - outDB.mode = dvo_catalog_catmode (CATMODE); - outDB.format = dvo_catalog_catformat (CATFORMAT); - status = dvo_image_lock (&outDB, ImageCat, 3600.0, LCK_XCLD); // shorter timeout? - if (!status) Shutdown ("ERROR: failure to lock image catalog %s", outDB.filename); - - /* load the image table */ - if (outDB.dbstate == LCK_EMPTY) { - dvo_image_create (&outDB, GetZeroPoint()); - } else { - if (!dvo_image_load (&outDB, VERBOSE, TRUE)) { - Shutdown ("can't read output image catalog %s", outDB.filename); - } - } - - // convert database table to internal structure & add to output image db - dvo_image_merge_dbs(IDmap, &outDB, &inDB); - dvo_image_unlock (&inDB); // unlock input - - SetProtect (TRUE); - dvo_image_save (&outDB, VERBOSE); - SetProtect (FALSE); - dvo_image_unlock (&outDB); // unlock output - - return TRUE; -} Index: trunk/Ohana/src/dvomerge/src/dvomergeUpdate_threaded.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvomergeUpdate_threaded.c (revision 29938) +++ trunk/Ohana/src/dvomerge/src/dvomergeUpdate_threaded.c (revision 29938) @@ -0,0 +1,313 @@ +# include "dvomerge.h" +# include + +enum { + TS_WAIT, + TS_RUN, + TS_DONE, + TS_FAIL, + TS_EXIT, +}; + +typedef struct { + int iThread; + int state; + // elements used to carry data to the threadjob + char *name; + char *filename; + SkyRegion *region; + int NsecfiltInput; + int NsecfiltOutput; + IDmapType *IDmap; + SkyTable *outsky; + int *secfiltMap; + int depth; +} ThreadData; + +void *dvomergeUpdate_threadjob (void *inputData) { + + off_t j; + Catalog incatalog, outcatalog; + SkyList *outlist; + + ThreadData *threadData = (ThreadData *) inputData; + + while (TRUE) { + + while (threadData->state != TS_RUN) { + usleep (100); + } + + if (VERBOSE) fprintf (stderr, "input: %s\n", threadData->name); + + // load / create output catalog (if catalog does not exist, it will be created) + LoadCatalog (&incatalog, threadData->region, threadData->filename, "r", threadData->NsecfiltInput); + // skip empty input catalogs + if (!incatalog.Naves_disk) { + dvo_catalog_unlock (&incatalog); + dvo_catalog_free (&incatalog); + threadData->state = TS_DONE; + continue; + } + + // combine only tables at equal or larger depth. + + // NOTE: this requirement is especially important for the threaded version: if the + // output catalog is guaranteed to be the same size or finer than the input, then two + // input catalogs cannot collide on their target output catalog. + + // load in all of the tables from input for this region + // SkyListByBounds will return neighbor catalogs if the boundaries exactly match (due to rounding). Since the regions are not infinitely small, + // compare to a slightly reduced footprint + float dPos = 2.0/3600.0; + outlist = SkyListByBounds (threadData->outsky, threadData->depth, threadData->region->Rmin + dPos, threadData->region->Rmax - dPos, threadData->region->Dmin + dPos, threadData->region->Dmax - dPos); + for (j = 0; (j < outlist[0].Nregions) && (threadData->state != TS_FAIL); j++) { + if (VERBOSE) fprintf (stderr, "output : %s\n", outlist[0].regions[j][0].name); + + // load input catalog + LoadCatalog (&outcatalog, outlist[0].regions[j], outlist[0].filename[j], "w", threadData->NsecfiltOutput); + + dvo_update_image_IDs (threadData->IDmap, &incatalog); + merge_catalogs_old (&threadData->outsky->regions[j], &outcatalog, &incatalog, RADIUS, threadData->secfiltMap); + + outcatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; + + // if we receive a signal which would cause us to exit, wait until the full catalog is written + SetProtect (TRUE); + if (!dvo_catalog_save (&outcatalog, VERBOSE)) { + fprintf (stderr, "ERROR: failed to save catalog %s\n", outlist[0].filename[j]); + threadData->state = TS_FAIL; + continue; + } + SetProtect (FALSE); + + dvo_catalog_unlock (&outcatalog); + dvo_catalog_free (&outcatalog); + + fprintf (stderr, "merged %s into %s\n", threadData->name, outlist[0].regions[j][0].name); + } + SkyListFree (outlist); + + dvo_catalog_unlock (&incatalog); + dvo_catalog_free (&incatalog); + threadData->state = TS_DONE; + } + return NULL; +} + +int dvomergeUpdate_threaded (int argc, char **argv) { + + int j, launched, done, Nwait, Nrun; + + ThreadData *threadData; + pthread_t *threads; + + int depth, CONTINUE; + off_t i, Ns, Ne; + SkyTable *outsky, *insky; + SkyList *inlist; + char filename[256], *input, *output; + IDmapType IDmap; + PhotCodeData *inputPhotcodes; + PhotCodeData *outputPhotcodes; + int *secfiltMap = NULL; + int NsecfiltInput, NsecfiltOutput; + + double dtime; + struct timeval start, stop; + gettimeofday (&start, NULL); + + CONTINUE = FALSE; + if (strcasecmp (argv[2], "into")) dvomerge_usage(); + if (argc == 5) { + if (strcasecmp (argv[4], "continue")) dvomerge_usage(); + CONTINUE = TRUE; + } + + input = argv[1]; + output = argv[3]; + + if (ALTERNATE_PHOTCODE_FILE) { + fprintf (stderr, "cannot specify photcodes when merging into an existing catdir\n"); + exit (1); + } + + SetPhotcodeTable(NULL); + sprintf (filename, "%s/Photcodes.dat", input); + if (!LoadPhotcodes (filename, NULL, FALSE)) { + fprintf (stderr, "error reading input database directory %s\n", input); + exit (1); + } + inputPhotcodes = GetPhotcodeTable(); + NsecfiltInput = GetPhotcodeNsecfilt(); + + // since we are merging the input db into the output db, the output defines the photcode + // table & db layout but, this requires the output to exist. if it does not, instead use the + // input. + + SetPhotcodeTable(NULL); + sprintf (filename, "%s/Photcodes.dat", output); + if (LoadPhotcodes (filename, NULL, FALSE)) { + + outputPhotcodes = GetPhotcodeTable(); + NsecfiltOutput = GetPhotcodeNsecfilt(); + + secfiltMap = GetSecFiltMap(outputPhotcodes, inputPhotcodes); + if (!secfiltMap) { + fprintf (stderr, "failed to map input secfilt photcodes to output photcodes table\n"); + exit (1); + } + } else { + fprintf(stderr, "%s not found using photcodes from %s/Photcodes.dat\n", filename, input); + outputPhotcodes = inputPhotcodes; + NsecfiltOutput = NsecfiltInput; + + if (!check_dir_access (output, VERBOSE)) { + fprintf (stderr, "error creating output database directory %s\n", output); + exit (1); + } + SetPhotcodeTable(inputPhotcodes); + if (!SavePhotcodesFITS (filename)) { + fprintf (stderr, "error saving photcode table in %s/Photcodes.dat\n", output); + exit (1); + } + } + + if (CONTINUE) { + // need to determine the mapping from the input to the output images + dvomergeImagesGetMap (&IDmap, input, output); + } else { + dvomergeImagesUpdate (&IDmap, input, output); + } + + // load the sky table for the existing database + insky = SkyTableLoadOptimal (input, NULL, NULL, FALSE, SKY_DEPTH_HST, VERBOSE); + if (!insky) { + Shutdown ("can't read SkyTable for %s", input); + } + SkyTableSetFilenames (insky, input, "cpt"); + + // XXX apply this...generate the subset matching the user-selected region + inlist = SkyListByPatch (insky, -1, &UserPatch); + + // generate an output table populated at the desired depth + outsky = SkyTableLoadOptimal (output, NULL, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); + if (!outsky) { + Shutdown ("can't read or create SkyTable for %s", output); + } + SkyTableSetFilenames (outsky, output, "cpt"); + + // loop over the populatable output tables; check for data in input in the corresponding regions + + SkyListPopulatedRange (&Ns, &Ne, inlist, 0); + depth = inlist[0].regions[Ns][0].depth; + + SetPhotcodeTable(NULL); + + ALLOCATE(threads, pthread_t, NTHREADS); + ALLOCATE(threadData, ThreadData, NTHREADS); + for (i = 0; i < NTHREADS; i++) { + threadData[i].state = TS_WAIT; + threadData[i].iThread = i; + } + for (i = 0; i < NTHREADS; i++) { + pthread_create (&threads[i], NULL, &dvomergeUpdate_threadjob, (void *) &threadData[i]); + } + + // loop over the populated input regions + for (i = 0; i < inlist[0].Nregions; i++) { + if (!inlist[0].regions[i][0].table) continue; + + launched = FALSE; + while (!launched) { + + // are any threads done? can we harvest them? + for (j = 0; j < NTHREADS; j++) { + if (threadData[j].state == TS_FAIL) goto failure; + if (threadData[j].state != TS_DONE) continue; + // fprintf (stderr, "harvested thread %d\n", threadData[j].iThread); + threadData[j].state = TS_WAIT; + } + + // are any threads available? can we launch this job? + for (j = 0; !launched && (j < NTHREADS); j++) { + if (threadData[j].state == TS_FAIL) goto failure; + if (threadData[j].state != TS_WAIT) continue; + + threadData[j].name = inlist[0].regions[i][0].name; + threadData[j].filename = inlist[0].filename[i]; + threadData[j].region = inlist[0].regions[i]; + threadData[j].NsecfiltInput = NsecfiltInput; + threadData[j].NsecfiltOutput = NsecfiltOutput; + threadData[j].IDmap = &IDmap; + threadData[j].outsky = outsky; + threadData[j].secfiltMap = secfiltMap; + threadData[j].depth = depth; + threadData[j].state = TS_RUN; + launched = TRUE; + } + if (!launched) usleep (100); + } + } + + // wait for all threads to be done + done = FALSE; + while (!done) { + // are any threads done? can we harvest them? + for (j = 0; j < NTHREADS; j++) { + if (threadData[j].state == TS_FAIL) goto failure; + if (threadData[j].state != TS_DONE) continue; + threadData[j].state = TS_WAIT; + } + + // how many are waiting now? + Nwait = 0; + for (j = 0; j < NTHREADS; j++) { + if (threadData[j].state == TS_FAIL) goto failure; + if (threadData[j].state != TS_WAIT) continue; + Nwait ++; + } + + if (Nwait < NTHREADS) { + usleep (100); + } else { + done = TRUE; + } + } + + // save the output sky table copy + sprintf (filename, "%s/SkyTable.fits", output); + check_file_access (filename, TRUE, TRUE, VERBOSE); + if (!SkyTableSave (outsky, filename)) { + fprintf (stderr, "ERROR: failed to save sky table for %s\n", output); + exit (1); + } + + gettimeofday (&stop, NULL); + dtime = DTIME (stop, start); + fprintf (stderr, "SUCCESS: elapsed time %9.4f sec\n", dtime); + + exit (0); + +failure: + + // wait for all threads to be done + done = FALSE; + while (!done) { + + // how many are still running? + Nrun = 0; + for (j = 0; j < NTHREADS; j++) { + if (threadData[j].state != TS_RUN) continue; + Nrun ++; + } + + if (Nrun) { + usleep (100); + } else { + done = TRUE; + } + } + exit (1); +} + Index: trunk/Ohana/src/dvomerge/src/dvorepair.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvorepair.c (revision 29938) +++ trunk/Ohana/src/dvomerge/src/dvorepair.c (revision 29938) @@ -0,0 +1,14 @@ +# include "dvomerge.h" + +int main (int argc, char **argv) { + + dvorepair_help(argc, argv); + + if (!strcmp(argv[1], "-fix-cpt")) dvorepairFixCPT(argc, argv); + // if (!strcmp(argv[1], "-fix-tables")) dvorepairFixTables(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); +} Index: trunk/Ohana/src/dvomerge/src/dvorepairCPT.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvorepairCPT.c (revision 29938) +++ trunk/Ohana/src/dvomerge/src/dvorepairCPT.c (revision 29938) @@ -0,0 +1,190 @@ +# include "dvomerge.h" + +# define myAbort(MSG) { 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: +// * 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 +// * + +// XXX absorb this code into a unified dvorepair program + +int main (int argc, char **argv) { + + off_t Nmeasure, Nimage; + int i, 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; + + if (argc != 5) { + fprintf (stderr, "USAGE: dvorepair (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 (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: trunk/Ohana/src/dvomerge/src/dvorepairDeleteImageList.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvorepairDeleteImageList.c (revision 29938) +++ trunk/Ohana/src/dvomerge/src/dvorepairDeleteImageList.c (revision 29938) @@ -0,0 +1,568 @@ +# 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; + Qmin = 360.0; + Qmax = 0.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: trunk/Ohana/src/dvomerge/src/dvorepairFixCPT.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvorepairFixCPT.c (revision 29938) +++ trunk/Ohana/src/dvomerge/src/dvorepairFixCPT.c (revision 29938) @@ -0,0 +1,160 @@ +# 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 FixCPTfile (char *rootname, Image *image, off_t Nimage); + +int dvorepairFixCPT (int argc, char **argv) { + + int i, j, N, Nfiles, NFILES; + off_t Nimage; + FITS_DB db; // database handle pointing to input image table + Image *image; + char *listFilename, *imageFilename; + char **files; + + 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 != 3) { + fprintf (stderr, "USAGE: dvorepair -fix-cpt (images) (filelist)\n"); + exit (2); + } + + imageFilename = argv[1]; + listFilename = argv[2]; + + // 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); + + FILE *file = fopen (listFilename, "r"); + myAssert(file, "failed to open list"); + + Nfiles = 0; + NFILES = 100; + ALLOCATE(files, char *, NFILES); + for (i = 0; i < NFILES; i++) { + ALLOCATE (files[i], char, 256); + memset (files[i], 0, 256); + } + + for (i = 0; TRUE; i++) { + if (fscanf (file, "%s", files[i]) == EOF) { + break; + } + if (i == NFILES - 1) { + NFILES += 100; + REALLOCATE(files, char *, NFILES); + for (j = i + 1; j < NFILES; j++) { + ALLOCATE (files[j], char, 256); + memset (files[j], 0, 256); + } + } + } + Nfiles = i; + fclose(file); + + for (i = 0; i < Nfiles; i++) { + FixCPTfile(files[i], image, Nimage); + } + + exit (0); +} + +// fix the CPT file. in the process, we rewrite the CPM file so that the averef entries are correctly set. +int FixCPTfile (char *rootFilename, Image *image, off_t Nimage) { + + char cptFilenameSrc[256], cptFilenameTgt[256], cpsFilenameSrc[256], cpsFilenameTgt[256], cpmFilenameSrc[256], cpmFilenameTgt[256]; + + Header cpmHeaderPHU, cpmHeaderTBL; + FTable cpmFtable; + FILE *cpmFile = NULL; + + off_t Nbytes, Nmeasure; + + Measure *measure; + + Matrix matrix; + + char catformat; + + cpmFtable.header = &cpmHeaderTBL; + + sprintf (cpmFilenameSrc, "%s.cpm", rootFilename); + sprintf (cpmFilenameTgt, "%s.cpm.fixed", rootFilename); + sprintf (cptFilenameSrc, "%s.cpt", rootFilename); + sprintf (cptFilenameTgt, "%s.cpt.fixed", rootFilename); + sprintf (cpsFilenameSrc, "%s.cps", rootFilename); + sprintf (cpsFilenameTgt, "%s.cps.fixed", rootFilename); + + // open cpm file + cpmFile = fopen(cpmFilenameSrc, "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 */ + // XXX this should not be running on broken CPM files + 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"); + + // 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, measure, Nmeasure, image, Nimage, catformat); + + // convert internal to external format + if (!MeasureToFtable (&cpmFtable, measure, Nmeasure, 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, Nmeasure, 0, Nmeasure)) { + myAbort("can't write table data"); + } + + gfits_free_header (&cpmHeaderPHU); + gfits_free_header (&cpmHeaderTBL); + free (measure); + fclose(cpmFile); + + return TRUE; +} Index: trunk/Ohana/src/dvomerge/src/dvorepairFixImages.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvorepairFixImages.c (revision 29938) +++ trunk/Ohana/src/dvomerge/src/dvorepairFixImages.c (revision 29938) @@ -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: trunk/Ohana/src/dvomerge/src/dvorepairFixTables.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvorepairFixTables.c (revision 29938) +++ trunk/Ohana/src/dvomerge/src/dvorepairFixTables.c (revision 29938) @@ -0,0 +1,423 @@ +# include "dvomerge.h" + +// repair short cpm files, fix inconsistent cpt / cps / cpm files +// * 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; + + Measure *measure; + Measure *measureNew; + + SkyTable *insky; + SkyList *inlist; + + char *cpmFilenameSrc = NULL; + char *cptFilenameSrc = NULL; + char *cpsFilenameSrc = NULL; + char *cpmFilenameTgt = NULL; + char *cptFilenameTgt = NULL; + char *cpsFilenameTgt = NULL; + + Header cpmHeaderPHU; + Header cpmHeaderTBL; + FTable cpmFtable; + + FILE *cpmFile = NULL; + + char catformat; + + N = get_argument (argc, argv, "-fix-tables"); + myAssert(N == 1, "programming error: -fix-tables must be first from main"); + remove_argument (N, &argc, argv); + + // restrict to a portion of the sky + UserPatch.Rmin = 0; + UserPatch.Rmax = 360; + UserPatch.Dmin = -90; + UserPatch.Dmax = +90; + if ((N = get_argument (argc, argv, "-region"))) { + remove_argument (N, &argc, argv); + UserPatch.Rmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Rmax = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmax = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + if (argc != 2) { + fprintf (stderr, "USAGE: dvorepair -fix-tables (catdir) [-region Rmin Rmax Dmin Dmax]\n"); + fprintf (stderr, " catdir : database of interest\n"); + exit (2); + } + + char *catdir = argv[1]; + + // 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"); + inlist = SkyListByPatch (insky, -1, &UserPatch); + + 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); + + // loop over the tables + 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; + + measure = LoadTableCPM(cpmFilenameSrc, &Nmeasure, &Nexpect); + + // 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); + + // if the file is short, create + if (Nmeasure != Nexpect) { + fprintf (stderr, "deleting %d of %d (keep %d) detections from %s -> %s\n", (int) Ndelete, (int) Nvalid, (int) NmeasureNew, cpmFilenameSrc, cpmFilenameTgt); + + // 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); + } + } + SkyListFree(inlist); + SkyTableFree(insky); + + free(cpmFilenameSrc); + free(cptFilenameSrc); + free(cpsFilenameSrc); + free(cpmFilenameTgt); + free(cptFilenameTgt); + free(cpsFilenameTgt); + + exit (0); +} + +// load this CPM file : is it short? +Measure *LoadTableCPM (char *cpmFilenameSrc, off_t *nmeasure, off_t *nexpect) { + + fprintf (stderr, "check %s\n", cpmFilenameSrc); + + // open cpm file + cpmFile = fopen(cpmFilenameSrc, "r"); + if (!cpmFile) continue; + + // 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); + gfits_scan(&cpmHeaderTBL, "NAXIS2", "%d", 1, &Nrows); + + 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); + + myAssert(Nrows >= Nvalid, "how can we find more entries than expected??"); + + *nmeasure = Nvalid; + *nexpect = Nrows; + + return measure; +} + +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: trunk/Ohana/src/dvomerge/src/dvorepairImageVsMeasure.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvorepairImageVsMeasure.c (revision 29938) +++ trunk/Ohana/src/dvomerge/src/dvorepairImageVsMeasure.c (revision 29938) @@ -0,0 +1,178 @@ +# include "dvomerge.h" + +# 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: +// * load the Images.dat file +// * create an index for imageID (seq = imageIDindex[i]) +// * scan over all catalog files +// * load the cpm file (pad if short, identify the padded section) +// * loop over detections: increment detection count for each imageID + +int main (int argc, char **argv) { + + off_t i, j, Nmeasure, Nvalid, Nimage, Nindex, *imageIdx, index, imageSeq; + int 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; + + if (argc != 3) { + fprintf (stderr, "USAGE: dvorepair (catdir) (Ntol)\n"); + fprintf (stderr, " catdir : database of interest\n"); + fprintf (stderr, " Ntol : allow Ntol missing detections\n"); + exit (2); + } + + catdir = argv[1]; + Ntol = atoi(argv[2]); + + // load the image data + ALLOCATE(imageFilename, char, strlen(catdir) + 12); + sprintf (imageFilename, "%s/Images.dat", catdir); + if ((image = LoadImages (imageFilename, &Nimage)) == NULL) return (FALSE); + BuildChipMatch (image, Nimage); + + // 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: trunk/Ohana/src/dvomerge/src/dvorepairImagesVsMeasures.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvorepairImagesVsMeasures.c (revision 29938) +++ trunk/Ohana/src/dvomerge/src/dvorepairImagesVsMeasures.c (revision 29938) @@ -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: trunk/Ohana/src/dvomerge/src/dvoverify.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvoverify.c (revision 29181) +++ trunk/Ohana/src/dvomerge/src/dvoverify.c (revision 29938) @@ -8,170 +8,210 @@ */ +int VerifyTableFile (char *filename); + +# define DEBUG 0 + int main (int argc, char **argv) { - char filename[256], *input, *output; - - int depth; - off_t i, j, Ns, Ne; - SkyTable *outsky, *insky; + char filename[1024]; + + int N, Nbad; + off_t i; + SkyTable *insky; SkyList *inlist; - Catalog incatalog, outcatalog; - - SetSignals (); - dvoconvert_help (argc, argv); - ConfigInit (&argc, argv); - dvoconvert_args (&argc, argv); - - if (strcasecmp (argv[2], "to")) dvoconvert_usage(); - input = argv[1]; - output = argv[3]; - - // load input images, save to output images - dvoConvert_copy_images (input, output); - - // copy photcode table - { - // the first input defines the photcode table & db layout - sprintf (filename, "%s/Photcodes.dat", input); - if (!LoadPhotcodes (filename, NULL, FALSE)) { - fprintf (stderr, "error loading photcode table %s\n", filename); - exit (1); - } - // save the photcodes in the output catdir - sprintf (filename, "%s/Photcodes.dat", output); - if (!check_file_access (filename, TRUE, TRUE, VERBOSE)) { - fprintf (stderr, "error creating output catdir %s\n", output); - exit (1); - } - if (!SavePhotcodesFITS (filename)) { - fprintf (stderr, "error saving photcode table %s\n", filename); - exit (1); - } - } - - // copy skytable - { - // load the sky table for the existing database - insky = SkyTableLoadOptimal (input, NULL, NULL, FALSE, SKY_DEPTH_HST, VERBOSE); - SkyTableSetFilenames (insky, input, "cpt"); - - // generate an output table populated at the desired depth - // XXX force this to match the input? - outsky = SkyTableLoadOptimal (output, NULL, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); - SkyTableSetFilenames (outsky, output, "cpt"); - - SkyTablePopulatedRange (&Ns, &Ne, insky, 0); - depth = insky[0].regions[Ns].depth; - // XXX this seems to imply that insky is a uniform depth... - } - - // loop over all input catalogs, save to output catalogs - // loop over the populatable output regions - for (i = 0; i < outsky[0].Nregions; i++) { - if (!outsky[0].regions[i].table) continue; - if (VERBOSE) fprintf (stderr, "output: %s\n", outsky[0].regions[i].name); - - // load / create output catalog - LoadCatalog (&outcatalog, &outsky[0].regions[i], outsky[0].filename[i], "w"); - - // combine only tables at equal or larger depth - - // load in all of the tables from input for this region - inlist = SkyListByBounds (insky, depth, outsky[0].regions[i].Rmin, outsky[0].regions[i].Rmax, outsky[0].regions[i].Dmin, outsky[0].regions[i].Dmax); - for (j = 0; j < inlist[0].Nregions; j++) { - if (VERBOSE) fprintf (stderr, "input : %s\n", inlist[0].regions[j][0].name); - - // load input catalog - LoadCatalog (&incatalog, inlist[0].regions[j], inlist[0].filename[j], "r"); - - // skip empty input catalogs - if (!incatalog.Naves_disk) { - dvo_catalog_unlock (&incatalog); - dvo_catalog_free (&incatalog); - continue; - } - merge_catalogs_new (&outsky[0].regions[i], &outcatalog, &incatalog); - dvo_catalog_unlock (&incatalog); - dvo_catalog_free (&incatalog); - } - SkyListFree (inlist); - - outcatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; - dvo_catalog_save (&outcatalog, VERBOSE); - dvo_catalog_unlock (&outcatalog); - dvo_catalog_free (&outcatalog); - } - - // save the output sky table copy - sprintf (filename, "%s/SkyTable.fits", output); - check_file_access (filename, TRUE, TRUE, VERBOSE); - if (!SkyTableSave (outsky, filename)) { - fprintf (stderr, "ERROR: failed to save sky table for %s\n", output); + SkyRegion UserPatch; + // Catalog catalog; + + // restrict to a portion of the sky + UserPatch.Rmin = 0; + UserPatch.Rmax = 360; + UserPatch.Dmin = -90; + UserPatch.Dmax = +90; + if ((N = get_argument (argc, argv, "-region"))) { + remove_argument (N, &argc, argv); + UserPatch.Rmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Rmax = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmax = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + if (argc != 2) { + fprintf (stderr, "USAGE: dvoverify (catdir) [-region Rmin Rmax Dmin Dmax]\n"); + fprintf (stderr, " catdir : database of interest\n"); + exit (2); + } + + char *catdir = argv[1]; + + Nbad = 0; + + // XXX make this step optional + if (1) { + // check the photcode table + sprintf (filename, "%s/Photcodes.dat", catdir); + if (!VerifyTableFile (filename)) { + Nbad ++; + } + + // check the skytable + sprintf (filename, "%s/SkyTable.fits", catdir); + if (!VerifyTableFile (filename)) { + Nbad ++; + } + + // check the image table + sprintf (filename, "%s/Images.dat", catdir); + if (!VerifyTableFile (filename)) { + Nbad ++; + } + } + + // 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"); + inlist = SkyListByPatch (insky, -1, &UserPatch); + + // loop over all catalogs, save to output catalogs + for (i = 0; i < inlist[0].Nregions; i++) { + if (!inlist[0].regions[i][0].table) continue; + + if (i % 1000 == 0) fprintf (stderr, "."); + + sprintf (filename, "%s/%s.cpt", catdir, inlist[0].regions[i][0].name); + if (!VerifyTableFile (filename)) { + Nbad ++; + } + + sprintf (filename, "%s/%s.cps", catdir, inlist[0].regions[i][0].name); + if (!VerifyTableFile (filename)) { + Nbad ++; + } + + sprintf (filename, "%s/%s.cpm", catdir, inlist[0].regions[i][0].name); + if (!VerifyTableFile (filename)) { + Nbad ++; + } + } + + if (Nbad > 0) { + fprintf (stderr, "ERROR: %d files are bad\n", Nbad); exit (1); } + fprintf (stderr, "SUCCESS: no files are bad\n"); exit (0); } -int dvoConvert_copy_images (char *input, char *output) { - - FITS_DB inDB; - FITS_DB outDB; - int status; - - Image *images; - off_t Nimages; - off_t ID; - - /*** load output/Images.dat ***/ - sprintf (ImageCat, "%s/Images.dat", output); - outDB.mode = dvo_catalog_catmode (CATMODE); - outDB.format = dvo_catalog_catformat (CATFORMAT); - status = dvo_image_lock (&outDB, ImageCat, 3600.0, LCK_XCLD); // shorter timeout? - if (!status) Shutdown ("ERROR: failure to lock image catalog %s", outDB.filename); - - // output image table should not already exist - if (outDB.dbstate != LCK_EMPTY) { - Shutdown ("ERROR: image table %s already exists", outDB.filename); - } - dvo_image_create (&outDB, GetZeroPoint()); - - /*** load input/Images.dat ***/ - sprintf (ImageCat, "%s/Images.dat", input); - // inDB.mode = dvo_catalog_catmode (CATMODE); - // inDB.format = dvo_catalog_catformat (CATFORMAT); - status = dvo_image_lock (&inDB, ImageCat, 3600.0, LCK_XCLD); // shorter timeout? - if (!status) Shutdown ("ERROR: failure to lock image catalog %s", inDB.filename); - - // load the image table - if (inDB.dbstate == LCK_EMPTY) { - Shutdown ("can't find input image catalog %s", inDB.filename); - } - if (!dvo_image_load (&inDB, VERBOSE, TRUE)) { - Shutdown ("can't read input image catalog %s", inDB.filename); - } - - // convert the raw image table to Image type (byteswap if needed) - images = gfits_table_get_Image (&inDB.ftable, &Nimages, &inDB.swapped); - if (!images) { - fprintf (stderr, "ERROR: failed to read images\n"); - exit (2); - } - - // update additional metadata - gfits_scan (&inDB.header, "IMAGEID", OFF_T_FMT, 1, &ID); - gfits_modify (&outDB.header, "NIMAGES", OFF_T_FMT, 1, Nimages); - gfits_modify (&outDB.header, "IMAGEID", OFF_T_FMT, 1, ID); - - // copy input rows to output table - gfits_add_rows (&outDB.ftable, (char *) images, Nimages, sizeof(Image)); - - // write out the image table to disk - SetProtect (TRUE); - dvo_image_save (&outDB, VERBOSE); - SetProtect (FALSE); - dvo_image_unlock (&outDB); // unlock output - dvo_image_unlock (&inDB); // unlock input1 - +// is this file a consistent FITS file? +int VerifyTableFile (char *filename) { + + int status; + off_t Nbytes; + Header header; + + struct stat fileStats; + FILE *file; + + // does the file exist? + status = stat (filename, &fileStats); + if (status) { + // some error accessing the file. there is only one acceptable error: file not found + switch (errno) { + case ENOENT: + if (DEBUG) fprintf (stderr, "file does not exist, skipping %s\n", filename); + return TRUE; + case ENOMEM: + fprintf (stderr, "Out of memory: %s\n", filename); + return TRUE; + case EACCES: + fprintf (stderr, "Permission error on %s\n", filename); + return FALSE; + case EFAULT: + fprintf (stderr, "Bad address: %s\n", filename); + return FALSE; + case ELOOP: + fprintf (stderr, "Too many symbolic links encountered while traversing the path: %s\n", filename); + return FALSE; + case ENAMETOOLONG: + fprintf (stderr, "File name too long: %s\n", filename); + return FALSE; + case ENOTDIR: + fprintf (stderr, "A component of the path is not a directory: %s\n", filename); + return FALSE; + case EOVERFLOW: + fprintf (stderr, "file too large for program version: %s\n", filename); + return FALSE; + default: + fprintf (stderr, "unknown error: %s\n", filename); + return FALSE; + } + } + + // does it have any data? + if (fileStats.st_size == 0) { + fprintf (stderr, "file is empty: %s\n", filename); + return FALSE; + } + + // can we open it? + file = fopen(filename, "r"); + if (!file) { + fprintf (stderr, "unable to open valid file: %s\n", filename); + return FALSE; + } + + // scan all extentions + Nbytes = 0; + if (DEBUG) fprintf (stderr, "sizes: ("OFF_T_FMT" vs "OFF_T_FMT")\n", Nbytes, fileStats.st_size); + while (Nbytes < fileStats.st_size) { + + // Check on the PHU + if (!gfits_fread_header (file, &header)) { + fprintf (stderr, "unable to read PHU header for %s\n", filename); + fclose(file); + return (FALSE); + } + + // move to TBL header + Nbytes += header.datasize + gfits_data_size (&header); + if (DEBUG) fprintf (stderr, "sizes: ("OFF_T_FMT" vs "OFF_T_FMT")\n", Nbytes, fileStats.st_size); + if (Nbytes > fileStats.st_size) { + fprintf (stderr, "file is short ("OFF_T_FMT" vs "OFF_T_FMT"): %s\n", Nbytes, fileStats.st_size, filename); + gfits_free_header(&header); + fclose (file); + return FALSE; + } + gfits_free_header(&header); + + status = fseeko (file, Nbytes, SEEK_SET); + if (status) { + switch (errno) { + case EBADF: + fprintf (stderr, "something wrong with file handle: %s\n", filename); + fclose (file); + return FALSE; + case EINVAL: + fprintf (stderr, "invalid offset: %s\n", filename); + fclose (file); + return FALSE; + default: + fprintf (stderr, "other error in fseeko: %s\n", filename); + fclose (file); + return FALSE; + } + } + } + if (DEBUG) fprintf (stderr, "file is good: %s\n", filename); + fclose (file); return TRUE; } + + // gfits_scan(&cpmHeaderTBL, "NAXIS1", "%d", 1, &NbytesPerRow); + // gfits_scan(&cpmHeaderTBL, "NAXIS2", "%d", 1, &Nrows); + + Index: trunk/Ohana/src/dvomerge/src/help.c =================================================================== --- trunk/Ohana/src/dvomerge/src/help.c (revision 29181) +++ trunk/Ohana/src/dvomerge/src/help.c (revision 29938) @@ -4,4 +4,6 @@ fprintf (stderr, "USAGE: dvomerge (input1) and (input2) to (output)\n"); fprintf (stderr, " OR: dvomerge (input) into (output)\n"); + fprintf (stderr, " OR: dvomerge (input) into (output) continue\n"); + fprintf (stderr, " OR: dvomerge (input) into (output) from (list)\n"); fprintf (stderr, " [-region Rmin Rmax Dmin Dmax]\n"); exit (2); @@ -17,5 +19,5 @@ void dvosecfilt_usage(void) { - fprintf (stderr, "USAGE: dvosecfilt (input) -photcodes (photcodes.txt)\n"); + fprintf (stderr, "USAGE: dvosecfilt (catdir) (Nsecfilt)\n"); exit (2); @@ -33,5 +35,7 @@ fprintf (stderr, "USAGE\n"); fprintf (stderr, " dvomerge (input1) and (input2) to (output)\n"); - fprintf (stderr, " dvomerge (input) into (output)\n\n"); + fprintf (stderr, " dvomerge (input) into (output)\n"); + fprintf (stderr, " dvomerge (input) into (output) continue\n\n"); + fprintf (stderr, " dvomerge (input) into (output) from (list)\n\n"); fprintf (stderr, " merge DVO databases\n"); fprintf (stderr, " optional flags:\n"); @@ -39,4 +43,5 @@ fprintf (stderr, " -help : this list\n"); fprintf (stderr, " -h : this list\n\n"); + fprintf (stderr, " -region Rmin Rmax Dmin Dmax : region to merge\n\n"); exit (2); } @@ -75,7 +80,8 @@ fprintf (stderr, "USAGE\n"); - fprintf (stderr, " dvosecfilt (input) (Nsecfilt)\n\n"); + fprintf (stderr, " dvosecfilt (catdir) (Nsecfilt)\n\n"); - fprintf (stderr, " change number of secfilt entries in photcode table (updates secfilt tables only)\n"); + fprintf (stderr, " change number of secfilt entries in photcode table (updates secfilt tables (cps), NSECFILT in cpt files)\n"); + fprintf (stderr, " NOTE: the user must change the photcode table to reflect the change (use photcode-table -export / -import)\n"); fprintf (stderr, " optional flags:\n"); @@ -86,2 +92,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) (rootlist) - regenerate cpt & cps files from the cpm files\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: trunk/Ohana/src/dvomerge/src/match_image.c =================================================================== --- trunk/Ohana/src/dvomerge/src/match_image.c (revision 29938) +++ trunk/Ohana/src/dvomerge/src/match_image.c (revision 29938) @@ -0,0 +1,37 @@ +# include "dvomerge.h" + +off_t match_image (Image *image, off_t Nimage, unsigned int T, short int S) { + + off_t N, Nlo, Nhi, N1, N2; + + /* bracket first value of interest */ + Nlo = 0; Nhi = Nimage; + while (Nhi - Nlo > 10) { + N = 0.5*(Nlo + Nhi); + if (image[N].tzero < T) { + Nlo = N; + } else { + Nhi = N + 1; + } + } + N1 = Nlo; + + /* bracket last value of interest */ + Nlo = 0; Nhi = Nimage; + while (Nhi - Nlo > 10) { + N = 0.5*(Nlo + Nhi); + if (image[N].tzero > T) { + Nhi = N; + } else { + Nlo = N - 1; + } + } + N2 = Nhi; + + for (N = N1; N < N2; N++) { + if ((image[N].tzero == T) && (image[N].photcode == S)) { + return (N); + } + } + return (-1); +} Index: trunk/Ohana/src/dvomerge/src/merge_catalogs_old.c =================================================================== --- trunk/Ohana/src/dvomerge/src/merge_catalogs_old.c (revision 29181) +++ trunk/Ohana/src/dvomerge/src/merge_catalogs_old.c (revision 29938) @@ -28,5 +28,5 @@ unsigned int objID, catID; Coords tcoords; - + // struct timeval start, stop; // gettimeofday (&start, (void *) NULL); @@ -240,4 +240,6 @@ } else { for (k = 0; k < NsecfiltIn; k++) { + if (secfiltMap[k] < 0) continue; // skip secfilt entries from input that are not in the output + // index for this entry in output's secfilt list int outputIndex = n * NsecfiltOut + secfiltMap[k]; @@ -254,5 +256,5 @@ i++; } - // MARKTIME("find matched stars: %f sec for "OFF_T_FMT","OFF_T_FMT" stars\n", dtime, Nstars, Nave); + // MARKTIME("find matched stars: %f sec for "OFF_T_FMT","OFF_T_FMT" stars ("OFF_T_FMT" meas)\n", dtime, Nstars, Nave, Nmeas); /** incorporate unmatched image stars, if this star is in field of this catalog **/