Index: /branches/eam_branches/ipp-20140904/Ohana/src/dvomerge/Makefile =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/dvomerge/Makefile (revision 37494) +++ /branches/eam_branches/ipp-20140904/Ohana/src/dvomerge/Makefile (revision 37495) @@ -122,8 +122,12 @@ $(SRC)/dvorepairImagesVsMeasures.$(ARCH).o \ $(SRC)/dvorepairDeleteImageList.$(ARCH).o \ +$(SRC)/dvorepairDeleteImagesByExternID.$(ARCH).o \ +$(SRC)/dvorepairDeleteImagesByExternID.$(ARCH).o \ $(SRC)/dvorepairFixImages.$(ARCH).o \ $(SRC)/psps_ids.$(ARCH).o \ $(SRC)/LoadImages.$(ARCH).o \ $(SRC)/ReadDeleteList.$(ARCH).o \ +$(SRC)/ReadDeleteListExternID.$(ARCH).o \ +$(SRC)/myIndex.$(ARCH).o \ $(SRC)/match_image.$(ARCH).o \ $(SRC)/help.$(ARCH).o \ Index: /branches/eam_branches/ipp-20140904/Ohana/src/dvomerge/include/dvomerge.h =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/dvomerge/include/dvomerge.h (revision 37494) +++ /branches/eam_branches/ipp-20140904/Ohana/src/dvomerge/include/dvomerge.h (revision 37495) @@ -45,4 +45,16 @@ SkyRegion UserPatch; // used by MODE CAT +# ifndef MAX_INT +# define MAX_INT 2147483647 +# endif + +typedef struct { + int minID; + int maxID; + int *index; + int Nindex; + int NINDEX; +} myIndexType; + typedef struct { unsigned int Nmap; @@ -200,2 +212,12 @@ int dvoCreateID (char *catdir); +myIndexType *myIndexInit (); +int myIndexFree (myIndexType *myIndex); +int myIndexUpdateLimits (myIndexType *myIndex, int value); +int myIndexSetRange (myIndexType *myIndex); +int myIndexSetEntry (myIndexType *myIndex, int value, int entry); +int myIndexGetEntry (myIndexType *myIndex, int value); + +int dvorepairDeleteImagesByExternID (int argc, char **argv); +int RepairTableCPT_V1(char *cptFilenameSrc, char *cptFilenameTgt, char *cpsFilenameSrc, char *cpsFilenameTgt, Measure *measure, off_t Nmeasure, Image *image, off_t Nimage, myIndexType *imageIDindex, char catformat); +int *ReadDeleteListExternID(char *filename, int *nindex); Index: /branches/eam_branches/ipp-20140904/Ohana/src/dvomerge/src/ReadDeleteListExternID.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/dvomerge/src/ReadDeleteListExternID.c (revision 37495) +++ /branches/eam_branches/ipp-20140904/Ohana/src/dvomerge/src/ReadDeleteListExternID.c (revision 37495) @@ -0,0 +1,93 @@ +# include "dvomerge.h" + +int *ReadDeleteListExternID(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 we have 9 fields broken by 8 spaces: + space = c0; // pointer to track the space-separated words + indexPoint = c0; // pointer to the ID on this line + + for (j = 0; j < 8; j++) { + space = strchr(space, ' '); + if (!space) { + fprintf (stderr, "line only has %d word(s)\n", j + 1); + goto next_line; + } + space ++; + } + + // save the value we have found: + index = atoi(indexPoint); + indexList[Nindex] = index; + + // fprintf (stderr, "index: %d, line: %s\n", index, c0); + + Nindex ++; + if (Nindex == NINDEX) { + NINDEX += 1000; + REALLOCATE (indexList, int, NINDEX); + } + + next_line: + Nline ++; + c0 = c1 + 1; + } + } + + *nindex = Nindex; + return (indexList); +} + +/* Delete List Format: + image o5252g0035o.140669.cm.51017.smf[XY11.hdr] (219974) : 368 vs 375 + image o5252g0051o.140685.cm.51123.smf[XY11.hdr] (226417) : 350 vs 356 + image o5399g0207o.194950.cm.98351.smf[XY46.hdr] (1683277) : 10341 vs 10347 + image o5464g0284o.230438.cm.123774.smf[XY21.hdr] (2634913) : 214 vs 242 + image o5464g0305o.230459.cm.123795.smf[XY21.hdr] (2635758) : 232 vs 257 + image o5465g0306o.231171.cm.124201.smf[XY01.hdr] (2654941) : 0 vs 494 + image o5465g0306o.231171.cm.124201.smf[XY02.hdr] (2654942) : 10 vs 813 + image o5465g0306o.231171.cm.124201.smf[XY10.hdr] (2654947) : 0 vs 534 + image o5465g0306o.231171.cm.124201.smf[XY11.hdr] (2654948) : 0 vs 447 +*/ + Index: /branches/eam_branches/ipp-20140904/Ohana/src/dvomerge/src/dvorepair.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/dvomerge/src/dvorepair.c (revision 37494) +++ /branches/eam_branches/ipp-20140904/Ohana/src/dvomerge/src/dvorepair.c (revision 37495) @@ -3,7 +3,7 @@ int main (int argc, char **argv) { - fprintf (stderr, "this program needs to be updated to load old format Measure tables (pre PV1_V5) in which dR,dD are saved, not R,D\n"); - fprintf (stderr, "reminder: relastro can re-construct R,D from X,Y; FtableToMeasure and vice versa could just NAN those values\n"); - exit (2); + // fprintf (stderr, "this program needs to be updated to load old format Measure tables (pre PV1_V5) in which dR,dD are saved, not R,D\n"); + // fprintf (stderr, "reminder: relastro can re-construct R,D from X,Y; FtableToMeasure and vice versa could just NAN those values\n"); + // exit (2); dvorepair_help(argc, argv); @@ -13,4 +13,5 @@ if (!strcmp(argv[1], "-images-vs-measures")) dvorepairImagesVsMeasures(argc, argv); if (!strcmp(argv[1], "-delete-image-list")) dvorepairDeleteImageList(argc, argv); + if (!strcmp(argv[1], "-delete-images-by-extern-id")) dvorepairDeleteImagesByExternID(argc, argv); if (!strcmp(argv[1], "-fix-images")) dvorepairFixImages(argc, argv); dvorepair_help(0, NULL); Index: /branches/eam_branches/ipp-20140904/Ohana/src/dvomerge/src/dvorepairDeleteImagesByExternID.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/dvomerge/src/dvorepairDeleteImagesByExternID.c (revision 37495) +++ /branches/eam_branches/ipp-20140904/Ohana/src/dvomerge/src/dvorepairDeleteImagesByExternID.c (revision 37495) @@ -0,0 +1,622 @@ +# 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 dvorepairDeleteImagesByExternID (int argc, char **argv) { + + FITS_DB db; // database handle pointing to input image table + + int i, j, N, NmeasureNew, Ndelete, Nvalid; + off_t Nimage, Nmeasure; + int nPass, raPass; + + Image *image; + Measure *measure; + Measure *measureNew; + + char *cpmFilenameSrc = NULL; + char *cptFilenameSrc = NULL; + char *cpsFilenameSrc = NULL; + char *cpmFilenameTgt = NULL; + char *cptFilenameTgt = NULL; + char *cpsFilenameTgt = NULL; + + char *imageFilename = NULL; + + char catformat; + + N = get_argument (argc, argv, "-delete-images-by-extern-id"); + myAssert(N == 1, "programming error: -delete-images-by-extern-id must be first argument"); + remove_argument (N, &argc, argv); + + VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + + if (argc != 3) { + fprintf (stderr, "USAGE: dvorepair -delete-images-by-extern-id (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); // needed for astrometric tranformations + + // generate an index for imageIDs + myIndexType *externIDindex = myIndexInit(); + myIndexType *imageIDindex = myIndexInit(); + + // set the min and max ID values for the two strucures: + for (i = 0; i < Nimage; i++) { + if (!image[i].externID) continue; + myAssert (image[i].imageID, "image must have imageID > 0"); + myIndexUpdateLimits (externIDindex, image[i].externID); + myIndexUpdateLimits (imageIDindex, image[i].imageID); + } + + myIndexSetRange (externIDindex); + myIndexSetRange (imageIDindex); + + fprintf (stderr, "loaded %d images, image ID range: %d - %d, extern ID range: %d - %d\n", (int) Nimage, imageIDindex->minID, imageIDindex->maxID, externIDindex->minID, externIDindex->maxID); + + // deleteImage[i] will be TRUE if we want to delete that image + int *deleteImage = NULL; + ALLOCATE (deleteImage, int, Nimage); + for (i = 0; i < Nimage; i++) { + deleteImage[i] = FALSE; + } + + // assign the externIDs and imageIDs to their index structures + for (i = 0; i < Nimage; i++) { + if (!image[i].externID) continue; + + myIndexSetEntry (externIDindex, image[i].externID, i); + myIndexSetEntry (imageIDindex, image[i].imageID, i); + } + + // read the list of externIDs to delete + int NdeleteIDs; + int *deleteIDs = ReadDeleteListExternID(delList, &NdeleteIDs); + + for (i = 0; i < NdeleteIDs; i++) { + int externID = deleteIDs[i]; + myAssert(externID > 0, "extern ID to delete cannot be <= 0"); + + // we do not require all externIDs to be deleted to actually be in this db, + // but skip any externIDs not in this db (myIndexGetEntry returns -1) + int n = myIndexGetEntry(externIDindex, externID); + if (n < 0) continue; + + if (deleteImage[n]) { + fprintf (stderr, "multiple entries for extern ID %d == imageID %d, re-run\n", externID, image[n].imageID); + } + deleteImage[n] = TRUE; + + if (VERBOSE) fprintf (stderr, "delete extern ID: %d = image ID %d (seq %d)\n", externID, image[n].imageID, n); + } + + // find the RA & DEC range of the images we want to delete + double Rmin = 360.0; + double Rmax = 0.0; + double Dmin = +90.0; + double Dmax = -90.0; + double Qmin = 360.0; + double Qmax = 0.0; + for (i = 0; i < Nimage; i++) { + if (!deleteImage[i]) continue; + + double Rthis, Dthis, Qthis; + XY_to_RD(&Rthis, &Dthis, 0, 0, &image[i].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[i].NX, 0, &image[i].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[i].NY, &image[i].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[i].NX, image[i].NY, &image[i].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); + } + + double dQ = Qmax - Qmin; + double dR = Rmax - Rmin; + + // load the sky table for the existing database + SkyTable *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 + 0.1) { + UserPatch.Rmin = Rmin - 0.1; + UserPatch.Rmax = Rmax + 0.1; + } else { + if (raPass == 0) { + UserPatch.Rmin = 0.0; + UserPatch.Rmax = Qmax + 0.1; + } else { + UserPatch.Rmin = Qmin + 360.0 - 0.1; + UserPatch.Rmax = 360.0; + } + } + + SkyList *inlist = SkyListByPatch (insky, -1, &UserPatch); + + fprintf (stderr, "%d cpt regions affected\n", (int) inlist->Nregions); + + // SkyListPopulatedRange (&Ns, &Ne, inlist, 0); + // depth = inlist[0].regions[Ns][0].depth; + + // loop over the populated input regions + int 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.broken", catdir, inlist[0].regions[i][0].name); + sprintf (cptFilenameTgt, "%s/%s.cpt.broken", catdir, inlist[0].regions[i][0].name); + sprintf (cpsFilenameTgt, "%s/%s.cps.broken", catdir, inlist[0].regions[i][0].name); + + Header cpmHeaderPHU; + Header cpmHeaderTBL; + FTable cpmFtable; + + /*** read and examine the CPM file ***/ + { + if (VERBOSE) fprintf (stderr, "check %s\n", cpmFilenameSrc); + + cpmFtable.header = &cpmHeaderTBL; + + // open cpm file + 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 + int 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"); + } + + int NbytesPerRow; + gfits_scan(&cpmHeaderTBL, "NAXIS1", "%d", 1, &NbytesPerRow); + + // this function can only handle PS1_V5 and later formats + // NOTE: FtableToMeasure (and equivalent) free the data associated with cpmFtable. + // We are left with the data in measure. + measure = FtableToMeasure (&cpmFtable, NULL, &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++) { + int imageID = measure[j].imageID; + myAssert(imageID, "measure is missing an image ID"); + + int N = myIndexGetEntry(imageIDindex, imageID); + myAssert (N > -1, "impossible"); + + if (deleteImage[N]) { + Ndelete ++; + continue; + } + measureNew[NmeasureNew] = measure[j]; + NmeasureNew ++; + } + } + + // if we actually want to delete any measurements, write out a new cpm file + if (Ndelete > 0) { + + fprintf (stderr, "deleting %d of %d (keep %d) detections from %s -> %s\n", (int) Ndelete, (int) Nvalid, (int) NmeasureNew, cpmFilenameSrc, cpmFilenameTgt); + + // the CPT and CPS tables need to be regenerated. This must happen first because, in the process, we also update measure->averef + // what about other tables? + RepairTableCPT_V1(cptFilenameSrc, cptFilenameTgt, cpsFilenameSrc, cpsFilenameTgt, measureNew, NmeasureNew, image, Nimage, imageIDindex, catformat); + + // convert internal to external format + if (!MeasureToFtable (&cpmFtable, NULL, measureNew, NmeasureNew, catformat)) { + myAbort("trouble converting format"); + } + + // rename the old cpm file: + if (rename (cpmFilenameSrc, cpmFilenameTgt)) { + perror ("tried to rename file"); + exit (2); + } + + // create and write the output file + FILE *cpmFile = fopen(cpmFilenameSrc, "w"); + myAssert(cpmFile, "failed to open cpt file"); + + // write PHU header + if (!gfits_fwrite_header (cpmFile, &cpmHeaderPHU)) { + myAbort("can't write primary header"); + } + + Matrix matrix; + + // 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_table (&cpmFtable); + } else { + if (VERBOSE) fprintf (stderr, "nothing to delete in %s\n", cpmFilenameSrc); + } + + 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); + } + + myIndexFree (externIDindex); + myIndexFree (imageIDindex); + + free (imageFilename); + free (deleteImage); + free (deleteIDs); + + free(cpmFilenameSrc); + free(cptFilenameSrc); + free(cpsFilenameSrc); + free(cpmFilenameTgt); + free(cptFilenameTgt); + free(cpsFilenameTgt); + + SkyTableFree(insky); + + // XXX: NOTE: I am not deleting the images... + gfits_db_free (&db); + + exit (0); +} + +int RepairTableCPT_V1(char *cptFilenameSrc, char *cptFilenameTgt, char *cpsFilenameSrc, char *cpsFilenameTgt, Measure *measure, off_t Nmeasure, Image *image, off_t Nimage, myIndexType *imageIDindex, char catformat) { + + off_t *averefMatch; + off_t i, NaveMax, Naverage, NAVERAGE, NaverageOut, Nave, Nout, Nold; + int *found, Nsecfilt; + + Average *average, *averageOut; + + Matrix matrix; + + Header cptHeaderPHU; + Header cptHeaderTBL; + FTable cptFtable; + + 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 create new objects as needed + // we start with a valid, sorted dvo database, so we use averef to link objects + for (i = 0; i < Nmeasure; i++) { + Nave = measure[i].averef; + + // we only allocate as many as we need + 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)); + } + + // this measure matches an existing average, just check it is OK and bump Nmeasure + 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 + int n = myIndexGetEntry (imageIDindex, measure[i].imageID); + myAssert (n > -1, "impossible!"); + + // I actually have measure->R,D, I could just use those... + XY_to_RD (&average[Nave].R, &average[Nave].D, measure[i].Xccd, measure[i].Yccd, &image[n].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].refColorBlue = NAN; + average[Nave].refColorRed = NAN; + + 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; + + // modify measure.averef to match the new sequence + 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 + FILE *cptFile = fopen(cptFilenameSrc, "r"); + myAssert(cptFile, "failed to open cpt file"); + + // load the cpt header (use for CATID, RA,DEC range, filenames) + if (!gfits_fread_header (cptFile, &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"); + } + fclose(cptFile); + + // rename the old cpt file: + if (rename (cptFilenameSrc, cptFilenameTgt)) { + perror ("tried to rename file"); + exit (2); + } + + // create and write the output file + cptFile = fopen(cptFilenameSrc, "w"); + myAssert(cptFile, "failed to open cpt file"); + + // write PHU header + if (!gfits_fwrite_header (cptFile, &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 (cptFile, &matrix)) { + myAbort("can't write primary matrix"); + } + gfits_free_matrix (&matrix); + + // write the table data + if (!gfits_fwrite_ftable_range (cptFile, &cptFtable, 0, NaverageOut, 0, NaverageOut)) { + myAbort("can't write table data"); + } + fclose(cptFile); + + gfits_free_table (&cptFtable); + gfits_free_header (&cptHeaderPHU); + gfits_free_header (&cptHeaderTBL); + free (average); + free (averageOut); + + free (found); + free (averefMatch); + + { + Header cpsHeaderPHU; + Header cpsHeaderTBL; + FTable cpsFtable; + + SecFilt *secfilt = NULL; + + cpsFtable.header = &cpsHeaderTBL; + + // open source cpt file + FILE *cpsFile = fopen(cpsFilenameSrc, "r"); + myAssert(cpsFile, "failed to open cps file"); + + // load the cps header (use for CATID, RA,DEC range, filenames) + if (!gfits_fread_header (cpsFile, &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"); + } + fclose(cpsFile); + + // rename the old cpt file: + if (rename (cpsFilenameSrc, cpsFilenameTgt)) { + perror ("tried to rename file"); + exit (2); + } + + // create and write the output file + cpsFile = fopen(cpsFilenameSrc, "w"); + myAssert(cpsFile, "failed to open cps file"); + + // write PHU header + if (!gfits_fwrite_header (cpsFile, &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 (cpsFile, &matrix)) { + myAbort("can't write primary matrix"); + } + gfits_free_matrix (&matrix); + + // write the table data + if (!gfits_fwrite_ftable_range (cpsFile, &cpsFtable, 0, Nrows, 0, Nrows)) { + myAbort("can't write table data"); + } + fclose(cpsFile); + + gfits_free_table (&cpsFtable); + gfits_free_header (&cpsHeaderPHU); + gfits_free_header (&cpsHeaderTBL); + free (secfilt); + } + + return (TRUE); +} + Index: /branches/eam_branches/ipp-20140904/Ohana/src/dvomerge/src/myIndex.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/dvomerge/src/myIndex.c (revision 37495) +++ /branches/eam_branches/ipp-20140904/Ohana/src/dvomerge/src/myIndex.c (revision 37495) @@ -0,0 +1,72 @@ +# include "dvomerge.h" + +myIndexType *myIndexInit () { + + myIndexType *myIndex; + ALLOCATE (myIndex, myIndexType, 1); + + myIndex->minID = MAX_INT; + myIndex->maxID = 0; + + myIndex->NINDEX = 1000; + myIndex->Nindex = 0; + + ALLOCATE (myIndex->index, int, myIndex->NINDEX); + + return myIndex; +} + +int myIndexFree (myIndexType *myIndex) { + + free (myIndex->index); + free (myIndex); + return TRUE; +} + +// set minID and maxID externally +int myIndexUpdateLimits (myIndexType *myIndex, int value) { + + myIndex->maxID = MAX(value, myIndex->maxID); + myIndex->minID = MIN(value, myIndex->minID); + + return TRUE; +} + +// once minID and maxID are set, this function defines the range and inits +int myIndexSetRange (myIndexType *myIndex) { + + myIndex->Nindex = myIndex->maxID - myIndex->minID + 1; + myIndex->NINDEX = myIndex->Nindex; + + REALLOCATE (myIndex->index, int, myIndex->NINDEX); + + int i; + for (i = 0; i < myIndex->Nindex; i++) { + myIndex->index[i] = -1; + } + + return TRUE; +} + +// once minID and maxID are set, this function defines the range and inits +int myIndexSetEntry (myIndexType *myIndex, int value, int entry) { + + int n = value - myIndex->minID; + + myAssert (n >= 0, "impossible!"); + myAssert(myIndex->index[n] == -1, "stomping on an earlier index"); + + myIndex->index[n] = entry; + return n; +} + +// once minID and maxID are set, this function defines the range and inits +int myIndexGetEntry (myIndexType *myIndex, int value) { + + if (value < myIndex->minID) return -1; // not in this index + if (value > myIndex->maxID) return -1; // not in this index + + int n = value - myIndex->minID; + + return myIndex->index[n]; +}