Index: trunk/Ohana/src/dvomerge/Makefile =================================================================== --- trunk/Ohana/src/dvomerge/Makefile (revision 37762) +++ trunk/Ohana/src/dvomerge/Makefile (revision 37807) @@ -23,5 +23,7 @@ dvoverify_client : $(BIN)/dvoverify_client.$(ARCH) -all: dvomerge dvomerge_client dvoconvert dvosecfilt dvosecfilt_client dvoverify dvoverify_client +dvoutils : $(BIN)/dvoutils.$(ARCH) + +all: dvomerge dvomerge_client dvoconvert dvosecfilt dvosecfilt_client dvorepair dvoverify dvoverify_client dvoutils # $(SRC)/dvomergeContinue.$(ARCH).o @@ -120,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 \ @@ -149,5 +155,14 @@ $(BIN)/dvoverify_client.$(ARCH) : $(DVOVERIFY_CLIENT) -INSTALL = dvomerge dvomerge_client dvoconvert dvosecfilt dvosecfilt_client dvorepair dvoverify dvoverify_client +DVOUTILS = \ +$(SRC)/dvoutils.$(ARCH).o \ +$(SRC)/dvoutils_args.$(ARCH).o \ +$(SRC)/dvoutils_uniq_images.$(ARCH).o \ +$(SRC)/dvoutils_load_image_index.$(ARCH).o + +$(DVOUTILS) : $(INC)/dvoutils.h +$(BIN)/dvoutils.$(ARCH) : $(DVOUTILS) + +INSTALL = dvomerge dvomerge_client dvoconvert dvosecfilt dvosecfilt_client dvorepair dvoverify dvoverify_client dvoverify # dependancy rules for binary code ######################### Index: trunk/Ohana/src/dvomerge/include/dvomerge.h =================================================================== --- trunk/Ohana/src/dvomerge/include/dvomerge.h (revision 37762) +++ trunk/Ohana/src/dvomerge/include/dvomerge.h (revision 37807) @@ -41,7 +41,22 @@ int REPLACE_BY_PHOTCODE; int FORCE_MERGE; +char *UPDATE_CATFORMAT; + +int MATCHED_TABLES; char *SINGLE_CPT; 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 { @@ -138,4 +153,9 @@ int add_lens_link PROTO((Average *average, off_t *next, off_t Nlensing, off_t NLENSING)); Lensing *sort_lensing PROTO((Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing, off_t *next)); + +off_t *build_starpar_links PROTO((Average *average, off_t Naverage, StarPar *starpar, off_t Nstarpar)); +off_t *init_starpar_links PROTO((Average *average, off_t Naverage, StarPar *starpar, off_t Nstarpar)); +int add_star_link PROTO((Average *average, off_t *next, off_t Nstarpar, off_t NSTARPAR)); +StarPar *sort_starpar PROTO((Average *average, off_t Naverage, StarPar *starpar, off_t Nstarpar, off_t *next)); off_t *init_missing_links PROTO((Average *average, off_t Naverage, Missing *missing, off_t Nmissing)); @@ -200,2 +220,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: trunk/Ohana/src/dvomerge/include/dvoutils.h =================================================================== --- trunk/Ohana/src/dvomerge/include/dvoutils.h (revision 37807) +++ trunk/Ohana/src/dvomerge/include/dvoutils.h (revision 37807) @@ -0,0 +1,34 @@ +# include +# include +# include +# include +# include +# include + +/* solaris requires both of these instead of ip.h: + # include + # include +*/ + +/* linux is happy with this, not solaris */ +# include +# include +# include +# include + +typedef enum { + DVOUTILS_NONE, + DVOUTILS_UNIQ_IMAGES, +} DVOUTILS_OP_TYPE; + +int VERBOSE; +int DVOUTILS_OP; + +char *IMAGES_LIST; + +int main (int argc, char **argv); +int dvoutils_args (int *argc, char **argv); +int dvoutils_uniq_images (char *filename); +int *dvoutils_load_image_index (char *filename, int *nindex); + +Image *dvoutils_load_image_table (char *filename, int *nimage); Index: trunk/Ohana/src/dvomerge/src/ImageOps.c =================================================================== --- trunk/Ohana/src/dvomerge/src/ImageOps.c (revision 37762) +++ trunk/Ohana/src/dvomerge/src/ImageOps.c (revision 37807) @@ -18,5 +18,4 @@ } 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 37762) +++ trunk/Ohana/src/dvomerge/src/LoadCatalog.c (revision 37807) @@ -8,5 +8,5 @@ // always load all of the data (if any exists) - catalog[0].catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF | LOAD_LENSING | LOAD_LENSOBJ; + catalog[0].catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF | LOAD_LENSING | LOAD_LENSOBJ | LOAD_STARPAR; catalog[0].catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data Index: trunk/Ohana/src/dvomerge/src/ReadDeleteListExternID.c =================================================================== --- trunk/Ohana/src/dvomerge/src/ReadDeleteListExternID.c (revision 37807) +++ trunk/Ohana/src/dvomerge/src/ReadDeleteListExternID.c (revision 37807) @@ -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: trunk/Ohana/src/dvomerge/src/args.c =================================================================== --- trunk/Ohana/src/dvomerge/src/args.c (revision 37762) +++ trunk/Ohana/src/dvomerge/src/args.c (revision 37807) @@ -37,4 +37,11 @@ } + /* extra error messages */ + MATCHED_TABLES = FALSE; + if ((N = get_argument (*argc, argv, "-matched-tables"))) { + MATCHED_TABLES = TRUE; + remove_argument (N, argc, argv); + } + /* use a different photcode file to define mean values */ if ((N = get_argument (*argc, argv, "-photcode-file"))) { @@ -48,4 +55,11 @@ if ((N = get_argument (*argc, argv, "-force-merge"))) { FORCE_MERGE = TRUE; + remove_argument (N, argc, argv); + } + + UPDATE_CATFORMAT = NULL; + if ((N = get_argument (*argc, argv, "-update-catformat"))) { + remove_argument (N, argc, argv); + UPDATE_CATFORMAT = strcreate (argv[N]); remove_argument (N, argc, argv); } @@ -185,4 +199,11 @@ } + /* extra error messages */ + MATCHED_TABLES = FALSE; + if ((N = get_argument (*argc, argv, "-matched-tables"))) { + MATCHED_TABLES = TRUE; + remove_argument (N, argc, argv); + } + /* replace measurement, don't duplicate */ REPLACE_BY_PHOTCODE = FALSE; Index: trunk/Ohana/src/dvomerge/src/build_links.c =================================================================== --- trunk/Ohana/src/dvomerge/src/build_links.c (revision 37762) +++ trunk/Ohana/src/dvomerge/src/build_links.c (revision 37807) @@ -382,3 +382,153 @@ } - +/******** StarPar ********************************************************************************/ + +/* build the initial links assuming the table is sorted, + not partial, and has a correct set of average[].starparOffset,Nstarpar values */ +off_t *init_starpar_links (Average *average, off_t Naverage, StarPar *starpar, off_t Nstarpar) { + + off_t i, j, N; + off_t *next_star; + + N = 0; + + ALLOCATE (next_star, off_t, Nstarpar); + for (i = 0; i < Naverage; i++) { + if (!average[i].Nstarpar) continue; + for (j = 0; j < average[i].Nstarpar - 1; j++, N++) { + next_star[N] = N + 1; + if (N >= Nstarpar) { + fprintf (stderr, "WARNING: N out of bounds (1)\n"); + } + } + next_star[N] = -1; + if (N >= Nstarpar) { + fprintf (stderr, "WARNING: N out of bounds (2)\n"); + } + + if (N >= Nstarpar) { + fprintf (stderr, "overflow in init_starpar_links\n"); + abort (); + } + N++; + } + return (next_star); +} + +/* construct starpar links which are valid FOR THIS LOAD + * - if we have a full load, we will get links which can + * be used by other programs (eg, relphot, etc) + * - if we have a partial load, the links are only valid + * for that partial load + */ + +off_t *build_starpar_links (Average *average, off_t Naverage, StarPar *starpar, off_t Nstarpar) { + + off_t i, m, k, Nm, averef; + off_t *next_star; + + ALLOCATE (next_star, off_t, Nstarpar); + + /* reset the Nm, offset values for average */ + for (i = 0; i < Naverage; i++) { + average[i].starparOffset = -1; + average[i].Nstarpar = 0; + } + + for (Nm = 0; Nm < Nstarpar; Nm++) { + averef = starpar[Nm].averef; + m = average[averef].starparOffset; + next_star[Nm] = -1; + + if (m == -1) { /* no links yet for source */ + average[averef].starparOffset = Nm; + average[averef].Nstarpar = 1; + continue; + } + + for (k = 0; next_star[m] != -1; k++) { + m = next_star[m]; + if (m >= Nstarpar) { + fprintf (stderr, "WARNING: m out of bounds (1)\n"); + } + } + + average[averef].Nstarpar = k + 2; + next_star[m] = Nm; + if (m >= Nstarpar) { + fprintf (stderr, "WARNING: m out of bounds (2)\n"); + } + } + return (next_star); +} + +/* average[].starparOffset, average[].Nstarpar are valid within an addstar run */ +int add_star_link (Average *average, off_t *next_star, off_t Nstarpar, off_t NSTARPAR) { + + off_t k, m; + + /* if we have trouble, check validity of next_star[m] : m < Nstarpar */ + m = average[0].starparOffset; + + for (k = 0; k < average[0].Nstarpar - 1; k++) { + m = next_star[m]; + if (m >= NSTARPAR) { + fprintf (stderr, "WARNING: m out of bounds (3)\n"); + } + } + + /* set up references */ + next_star[Nstarpar] = -1; + if (Nstarpar >= NSTARPAR) { + fprintf (stderr, "WARNING: Nstarpar out of bounds (1)\n"); + } + + if (m == -1) { + average[0].starparOffset = Nstarpar; + } else { + next_star[m] = Nstarpar; + if (m >= NSTARPAR) { + fprintf (stderr, "WARNING: m out of bounds (4)\n"); + } + } + + return (TRUE); +} + +StarPar *sort_starpar (Average *average, off_t Naverage, StarPar *starpar, off_t Nstarpar, off_t *next_star) { + + off_t i, k, n, np, N; + StarPar *tmpstarpar; + + /* + for (i = 0; i < Naverage; i++) { + if (average[i].Nstarpar != 4) { + fprintf (stderr, "check %d %d %d\n", (int) i, (int) average[i].Nstarpar, (int) average[i].starparOffset); + } + } + */ + + /* fix order of StarPar (memory intensive, but fast) */ + N = 0; + ALLOCATE (tmpstarpar, StarPar, Nstarpar); + for (i = 0; i < Naverage; i++) { + if (!average[i].Nstarpar) continue; + n = average[i].starparOffset; + average[i].starparOffset = N; + for (k = 0; k < average[i].Nstarpar; k++, N++) { + if (n == -1) { + fprintf (stderr, "entry after %d has a problem\n", (int) np); + abort(); + } + tmpstarpar[N] = starpar[n]; + if (starpar[n].averef != i) abort(); + tmpstarpar[N].averef = i; + np = n; + n = next_star[n]; + } + } + free (starpar); + return (tmpstarpar); +} + + Index: trunk/Ohana/src/dvomerge/src/dvo_image_merge_dbs.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvo_image_merge_dbs.c (revision 37762) +++ trunk/Ohana/src/dvomerge/src/dvo_image_merge_dbs.c (revision 37807) @@ -113,9 +113,4 @@ create_IDmap_lookup (IDmap); - if (!out[0].swapped) { - gfits_convert_Image ((Image *) out[0].ftable.buffer, sizeof(Image), Nout); - out[0].swapped = TRUE; - } - Nout += Nimages; IDstart += Nimages; Index: trunk/Ohana/src/dvomerge/src/dvomergeUpdate_catalogs.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvomergeUpdate_catalogs.c (revision 37762) +++ trunk/Ohana/src/dvomerge/src/dvomergeUpdate_catalogs.c (revision 37807) @@ -58,8 +58,20 @@ } +# ifdef OLD_CODE + // NOTE: i was having trouble dropping objects on edges with the option below. Since an object can move outside the catalog grab the neighbors. // 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, -1, 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); +# endif + + if (MATCHED_TABLES) { + // if we know the output will have the same depth / layout as the input, just choose the matched table + outlist = SkyRegionByIndex (outsky, inlist[0].regions[i][0].index); + } else { + // Since an object can move outside the catalog grab the neighbors. + float dPos = 2.0/3600.0; + outlist = SkyListByBounds (outsky, -1, 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); + } OutputStatus *outstat = OutputStatusInit (outlist->Nregions); @@ -163,7 +175,11 @@ LoadCatalog (&outcatalog, outlist[0].regions[j], outcatalog.filename, "w", NsecfiltOutput); - // if no catalog already exists, use the input catalog to define the format + // if no catalog already exists, use the input catalog to define the format or the specified format if (outcatalog.Naves_disk == 0) { - outcatalog.catformat = incatalog.catformat; + if (UPDATE_CATFORMAT) { + outcatalog.catformat = dvo_catalog_catformat (UPDATE_CATFORMAT); + } else { + outcatalog.catformat = incatalog.catformat; + } } @@ -187,9 +203,13 @@ SetProtect (FALSE); + if (!dvo_catalog_unlock (&outcatalog)) { + fprintf (stderr, "ERROR: failed to unlock catalog %s\n", outlist[0].filename[j]); + exit (1); + } + if (!dvo_catalog_unlink_backup (&outcatalog, TRUE)) { fprintf (stderr, "WARNING: failed to remove backup for catalog %s\n", outlist[0].filename[j]); } - dvo_catalog_unlock (&outcatalog); dvo_catalog_free (&outcatalog); @@ -256,4 +276,5 @@ if (PARALLEL_INPUT) { snprintf (tmpline, DVO_MAX_PATH, "%s -parallel-input", command); strcpy (command, tmpline); } if (FORCE_MERGE) { snprintf (tmpline, DVO_MAX_PATH, "%s -force-merge", command); strcpy (command, tmpline); } + if (MATCHED_TABLES) { snprintf (tmpline, DVO_MAX_PATH, "%s -matched-tables", command); strcpy (command, tmpline); } // add some config variables: Index: trunk/Ohana/src/dvomerge/src/dvorepair.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvorepair.c (revision 37762) +++ trunk/Ohana/src/dvomerge/src/dvorepair.c (revision 37807) @@ -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: trunk/Ohana/src/dvomerge/src/dvorepairCPT.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvorepairCPT.c (revision 37762) +++ trunk/Ohana/src/dvomerge/src/dvorepairCPT.c (revision 37807) @@ -119,4 +119,5 @@ 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].R = ohana_normalize_angle (average[Nave].R); average[Nave].Nmeasure = 1; Index: trunk/Ohana/src/dvomerge/src/dvorepairDeleteImageList.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvorepairDeleteImageList.c (revision 37762) +++ trunk/Ohana/src/dvomerge/src/dvorepairDeleteImageList.c (revision 37807) @@ -114,8 +114,4 @@ 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); @@ -404,4 +400,5 @@ 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].R = ohana_normalize_angle (average[Nave].R); average[Nave].Nmeasure = 1; Index: trunk/Ohana/src/dvomerge/src/dvorepairDeleteImagesByExternID.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvorepairDeleteImagesByExternID.c (revision 37807) +++ trunk/Ohana/src/dvomerge/src/dvorepairDeleteImagesByExternID.c (revision 37807) @@ -0,0 +1,770 @@ +# 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 SaveImageTable (Image *image, off_t Nimage, char *catdir, FITS_DB *oldDB); +Image *DeleteSelectedImages (Image *image, off_t Nimage, int *deleteImage, myIndexType *imageIDindex, off_t *NnewImage); + +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++) { + myIndexUpdateLimits (imageIDindex, image[i].imageID); + + if (!image[i].externID) continue; + myAssert (image[i].imageID, "image must have imageID > 0"); + myIndexUpdateLimits (externIDindex, image[i].externID); + } + + 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++) { + myIndexSetEntry (imageIDindex, image[i].imageID, i); + + if (!image[i].externID) continue; + myIndexSetEntry (externIDindex, image[i].externID, i); + } + + // read the list of externIDs to delete + int NdeleteIDs; + int *deleteIDs = ReadDeleteListExternID(delList, &NdeleteIDs); + + int NdeleteImages = 0; + 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; + NdeleteImages ++; + + if (VERBOSE) fprintf (stderr, "delete %s extern ID: %d = image ID %d (seq %d)\n", image[n].name, externID, image[n].imageID, n); + } + + if (NdeleteImages == 0) { + fprintf (stderr, "no images to be deleted in this database, exiting\n"); + exit (0); + } + + // 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 + 0.1) { + 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); + + int NdeleteTotal = 0; + + 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"); + // XXX this case is valid if we have REF detections (no associated image) + + int N = myIndexGetEntry(imageIDindex, imageID); + if (N < 0) { + // this detection comes from a non-existant image, delete + Ndelete ++; + continue; + } + + if (deleteImage[N]) { + Ndelete ++; + continue; + } + measureNew[NmeasureNew] = measure[j]; + NmeasureNew ++; + } + } + + NdeleteTotal += Ndelete; + + // 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); + } + + fprintf (stderr, "Deleted %d detections\n", NdeleteTotal); + + off_t NnewImage; + Image *newImage = DeleteSelectedImages (image, Nimage, deleteImage, imageIDindex, &NnewImage); + + SaveImageTable (newImage, NnewImage, catdir, &db); + + myIndexFree (externIDindex); + myIndexFree (imageIDindex); + + free (imageFilename); + free (deleteImage); + free (deleteIDs); + + free(cpmFilenameSrc); + free(cptFilenameSrc); + free(cpsFilenameSrc); + free(cpmFilenameTgt); + free(cptFilenameTgt); + free(cpsFilenameTgt); + + SkyTableFree(insky); + + gfits_db_free (&db); + + exit (0); +} + +Image *DeleteSelectedImages (Image *image, off_t Nimage, int *deleteImage, myIndexType *imageIDindex, off_t *NnewImage) { + + // we need to remove the (WRP) images, but we also want to remove the DIS [PHU] image entries if they have + // had all of their children deleted + + // first, generate an array of the number of children for each parent + + int *Nchildren = NULL; + ALLOCATE (Nchildren, int, Nimage); + memset (Nchildren, 0, Nimage*sizeof(int)); + + int i; + for (i = 0; i < Nimage; i++) { + if (!image[i].parentID) continue; + if (strcmp(&image[i].coords.ctype[4], "-WRP")) { + fprintf (stderr, "warning: child with non WRP coords : %s\n", image[i].name); + } + + int n = myIndexGetEntry (imageIDindex, image[i].parentID); + myAssert (n != -1, "invalid parentID?"); + myAssert (n < Nimage, "invalid parentID??"); + Nchildren[n] ++; + } + + // now subtract the children which need to be deleted + for (i = 0; i < Nimage; i++) { + if (!deleteImage[i]) continue; + if (!image[i].parentID) continue; // should I warn on this? + + int n = myIndexGetEntry (imageIDindex, image[i].parentID); + myAssert (n != -1, "invalid parentID?"); + myAssert (n < Nimage, "invalid parentID??"); + Nchildren[n] --; + myAssert (Nchildren[n] >= 0, "over deleted???"); + } + + // now check for any parent and delete it + for (i = 0; i < Nimage; i++) { + if (strcmp(&image[i].coords.ctype[4], "-DIS")) continue; + myAssert (!image[i].parentID, "parent ID for DIS image??"); + if (!Nchildren[i]) deleteImage[i] = TRUE; + } + + fprintf (stderr, "deleting the following images:\n"); + + off_t Nnew = 0; + Image *newImage = NULL; + ALLOCATE (newImage, Image, Nimage); + for (i = 0; i < Nimage; i++) { + if (deleteImage[i]) { + fprintf (stderr, "DELETE: %s : %d : %d : %d\n", image[i].name, image[i].imageID, image[i].parentID, image[i].externID); + continue; + } + newImage[Nnew] = image[i]; + Nnew ++; + } + + REALLOCATE (newImage, Image, Nnew); + + *NnewImage = Nnew; + return newImage; +} + +int SaveImageTable (Image *image, off_t Nimage, char *catdir, FITS_DB *oldDB) { + + FITS_DB db; + + char ImageCatSrc[DVO_MAX_PATH]; + char ImageCatTgt[DVO_MAX_PATH]; + + // we are creating a new image table (what about db ID?) + snprintf (ImageCatSrc, DVO_MAX_PATH, "%s/Images.dat", catdir); + snprintf (ImageCatTgt, DVO_MAX_PATH, "%s/Images.dat.broken", catdir); + + // rename the old cpt file: + if (rename (ImageCatSrc, ImageCatTgt)) { + perror ("tried to rename file"); + exit (2); + } + + /* setup image table format and lock */ + db.mode = dvo_catalog_catmode ("SPLIT"); + db.format = dvo_catalog_catformat ("PS1_V5"); + int status = dvo_image_lock (&db, ImageCatSrc, 3600.0, LCK_XCLD); // shorter timeout? + if (!status) { + fprintf (stderr, "ERROR: failure to lock image catalog %s", db.filename); + exit (2); + } + + /* load or create the image table */ + if (db.dbstate != LCK_EMPTY) { + fprintf (stderr, "image table %s exists, exiting", db.filename); + exit (2); + } + dvo_image_create (&db, GetZeroPoint()); + + // I want to keep the old DVO_DBID + char dbID[33]; + int imageIDmax; + if (!gfits_scan (&oldDB->header, "DVO_DBID", "%s", 1, dbID)) { + fprintf (stderr, "source image table is missing DVO_DBID, exiting\n"); + exit (3); + } + if (!gfits_scan (&oldDB->header, "IMAGEID", "%d", 1, &imageIDmax)) { + fprintf (stderr, "max image ID is missing, exiting\n"); + exit (3); + } + gfits_modify (&db.header, "IMAGEID", "%d", 1, imageIDmax); + gfits_modify (&db.header, "DVO_DBID", "%s", 1, dbID); + + /* add the new image and save */ + dvo_image_addrows (&db, image, Nimage); + dvo_image_update (&db, VERBOSE); + dvo_image_unlock (&db); /* unlock? */ + + return TRUE; +} + +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!"); + + dvo_average_init (&average[Nave]); + + // 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].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); + + for (i = 0; i < Nrows; i++) { + dvo_secfilt_init (&secfilt[i]); + } + + /* 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: trunk/Ohana/src/dvomerge/src/dvorepairFixImages.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvorepairFixImages.c (revision 37762) +++ trunk/Ohana/src/dvomerge/src/dvorepairFixImages.c (revision 37807) @@ -130,4 +130,8 @@ int MarkMosaicsToDelete(Image *image, off_t Nimage, int *deleteIndex) { + myAbort ("this function needs to be redone with a concept to match the old BuildChipMatch"); + +# if (0) + off_t i; @@ -140,4 +144,6 @@ 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) @@ -237,4 +243,6 @@ free (ChipByMosaicID); +# endif + return TRUE; } Index: trunk/Ohana/src/dvomerge/src/dvorepairFixTables.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvorepairFixTables.c (revision 37762) +++ trunk/Ohana/src/dvomerge/src/dvorepairFixTables.c (revision 37807) @@ -259,4 +259,5 @@ 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].R = ohana_normalize_angle (average[Nave].R); average[Nave].Nmeasure = 1; Index: trunk/Ohana/src/dvomerge/src/dvoutils.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvoutils.c (revision 37807) +++ trunk/Ohana/src/dvomerge/src/dvoutils.c (revision 37807) @@ -0,0 +1,20 @@ +# include "dvoutils.h" + +int main (int argc, char **argv) { + + // check various options + dvoutils_args (&argc, argv); + + switch (DVOUTILS_OP) { + + case DVOUTILS_UNIQ_IMAGES: + dvoutils_uniq_images(IMAGES_LIST); + exit (0); + + default: + fprintf (stderr, "ERROR: unknown option\n"); + exit (2); + } + + exit (2); +} Index: trunk/Ohana/src/dvomerge/src/dvoutils_args.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvoutils_args.c (revision 37807) +++ trunk/Ohana/src/dvomerge/src/dvoutils_args.c (revision 37807) @@ -0,0 +1,37 @@ +# include "dvoutils.h" + +int dvoutils_args (int *argc, char **argv) { + + int N; + + VERBOSE = FALSE; + if ((N = get_argument (*argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, argc, argv); + } + if ((N = get_argument (*argc, argv, "-verbose"))) { + VERBOSE = TRUE; + remove_argument (N, argc, argv); + } + + DVOUTILS_OP = DVOUTILS_NONE; + if ((N = get_argument (*argc, argv, "-uniq-images"))) { + DVOUTILS_OP = DVOUTILS_UNIQ_IMAGES; + remove_argument (N, argc, argv); + IMAGES_LIST = strcreate (argv[N]); + remove_argument (N, argc, argv); + } + + if (DVOUTILS_OP == DVOUTILS_NONE) { + fprintf (stderr, "ERROR: no valid mode selected\n"); + exit (2); + } + + if (*argc != 1) { + fprintf (stderr, "USAGE: dvoutils -uniq-images (filename) [-v,-verbose]\n\n"); + fprintf (stderr, " -v : VERBOSE\n"); + exit (2); + } + + return TRUE; +} Index: trunk/Ohana/src/dvomerge/src/dvoutils_load_image_index.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvoutils_load_image_index.c (revision 37807) +++ trunk/Ohana/src/dvomerge/src/dvoutils_load_image_index.c (revision 37807) @@ -0,0 +1,97 @@ +# include "dvoutils.h" + +// load the array of EXTERN_ID from the ImageIndex.fits tables +int *dvoutils_load_image_index (char *filename, int *nindex) { + + int Ncol; + off_t Nrow; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + char type[16]; + + FILE *f = fopen (filename, "r"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image index file %s\n", filename); + return NULL; + } + + /* load in PHU segment (ignore) */ + if (!gfits_fread_header (f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image index header\n"); + fclose (f); + return NULL; + } + if (!gfits_fread_matrix (f, &matrix, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image index matrix\n"); + gfits_free_header (&header); + fclose (f); + return NULL; + } + + ftable.header = &theader; + + // load data for this header + if (!gfits_load_header (f, &theader)) { + fclose (f); + return NULL; + } + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) { + fclose (f); + return (NULL); + } + fclose (f); + + int *externID = gfits_get_bintable_column_data (&theader, &ftable, "EXTERN_ID", type, &Nrow, &Ncol); + myAssert (!strcmp(type, "int"), "wrong column type"); + + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + *nindex = Nrow; + + return externID; +} + +Image *dvoutils_load_image_table (char *filename, int *nimage) { + + FITS_DB db; + + 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); + exit (3); + } + + if (db.dbstate == LCK_EMPTY) { + fprintf (stderr, "note: image catalog is empty\n"); + gfits_db_free (&db); + return NULL; + } + + int status = dvo_image_load (&db, TRUE, FALSE); + gfits_db_close (&db); + + if (!status) { + fprintf (stderr, "problem loading image database table\n"); + exit (4); + } + + off_t Nimage; + Image *image = gfits_table_get_Image (&db.ftable, &Nimage, &db.swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (5); + } + + *nimage = Nimage; + return image; +} Index: trunk/Ohana/src/dvomerge/src/dvoutils_uniq_images.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvoutils_uniq_images.c (revision 37807) +++ trunk/Ohana/src/dvomerge/src/dvoutils_uniq_images.c (revision 37807) @@ -0,0 +1,169 @@ +# include "dvoutils.h" + +# ifndef MAX_INT +# define MAX_INT 2147483647 +# endif + +int dvoutils_uniq_images(char *filename) { + + // given a list of Images.dat tables, find the unique subset of extern_id + // format of the input list is just (dirname)/Images.dat. commented lines + // starting with # are allowed + + FILE *f = fopen (filename, "r"); + if (!f) { + fprintf (stderr, "cannot read input list %s\n", filename); + exit (1); + } + + int minID = -1; + int maxID = -1; + int *IDlist = NULL; + int *DBfound = NULL; + int *DBindex = NULL; + + int Ndb = 0; + int NDB = 100; + char **DBlist = NULL; + ALLOCATE (DBlist, char *, NDB); // list of the names of the image database tables. + + char imFile[1024]; + while (scan_line_maxlen (f, imFile, 1024) != EOF) { + + stripwhite (imFile); + + if (imFile[0] == '#') continue; + + if (strchr(imFile, ' ') != NULL) { + fprintf (stderr, "line has extra spaces\n%s\n", imFile); + exit (2); + } + + int Nimage; + // Image *image = dvoutils_load_image_table (imFile, &Nimage); + int *extern_id = dvoutils_load_image_index (imFile, &Nimage); + if (!extern_id) continue; + + DBlist[Ndb] = strcreate (imFile); + + // now loop over all of the images and accumulate the array of extern_id values (skip + // the ones with extern_id == 0) + + // things I need to track: + // minID, maxID (only allocate an array of length maxID - minID + 1 + padding) + + int myMaxID = 0; + int myMinID = MAX_INT; + + int i; + for (i = 0; i < Nimage; i++) { + if (!extern_id[i]) continue; + myMaxID = MAX(extern_id[i], myMaxID); + myMinID = MIN(extern_id[i], myMinID); + } + + fprintf (stderr, "read %s, %d images, %d - %d vs %d - %d\n", imFile, (int) Nimage, myMinID, myMaxID, minID, maxID); + + // if this is new, we treat it a bit differently + if (!IDlist) { + int Nindex = myMaxID - myMinID + 1; + ALLOCATE (IDlist, int, Nindex); // number of times this ID is seen. IDlist[i] corresponds to extern_id = i + minID + ALLOCATE (DBfound, int, Nindex); // index of first DB in which this ID is seen (if n = DBfound[i], the correspond ID was seen in DB n + ALLOCATE (DBindex, int, Nindex); // sequence of this image in the first DB for which this ID was found + + for (i = 0; i < Nindex; i++) { + IDlist[i] = 0; + DBfound[i] = -1; + DBindex[i] = -1; + } + + minID = myMinID; + maxID = myMaxID; + } else { + int newMinID = MIN(myMinID, minID); + int newMaxID = MAX(myMaxID, maxID); + + int Nindex = newMaxID - newMinID + 1; + + int oldNindex = maxID - minID + 1; + myAssert (oldNindex <= Nindex, "impossible!"); + + int *newIDlist = NULL; + int *newDBfound = NULL; + int *newDBindex = NULL; + + ALLOCATE (newIDlist, int, Nindex); + ALLOCATE (newDBfound, int, Nindex); + ALLOCATE (newDBindex, int, Nindex); + + for (i = 0; i < Nindex; i++) { + newIDlist[i] = 0; + newDBfound[i] = -1; + newDBindex[i] = -1; + } + + // if myMinID >= minID, newMinID = minID -> offset = 0 + // if myMinID < minID, newMinID = myMinID -> offset = minID - myMinID + int offset = minID - newMinID; + + // save the old values + memcpy (&newIDlist [offset], IDlist, oldNindex*sizeof(int)); + memcpy (&newDBfound[offset], DBfound, oldNindex*sizeof(int)); + memcpy (&newDBindex[offset], DBindex, oldNindex*sizeof(int)); + + free (IDlist); + free (DBfound); + free (DBindex); + + IDlist = newIDlist; + DBfound = newDBfound; + DBindex = newDBindex; + + minID = newMinID; + maxID = newMaxID; + } + + for (i = 0; i < Nimage; i++) { + if (!extern_id[i]) continue; + int n = extern_id[i] - minID; + if (IDlist[n] == 0) { + DBfound[n] = Ndb; + DBindex[n] = i; + } + IDlist[n] ++; + } + + if (0) { + fprintf (stderr, "------------ %s ----------\n", imFile); + + for (i = 0; i < Nimage; i++) { + if (!extern_id[i]) continue; + fprintf (stderr, "%d => %d\n", extern_id[i], extern_id[i] - minID); + } + + int Nindex = maxID - minID + 1; + + for (i = 0; i < Nindex; i++) { + if (!IDlist[i]) continue; + fprintf (stderr, "%d : %d : %d\n", i, i + minID, IDlist[i]); + } + } + + // we are done with the db file, bump the counter + Ndb ++; + CHECK_REALLOCATE (DBlist, char *, NDB, Ndb, 100); + } + + int Nindex = maxID - minID + 1; + + int i; + for (i = 0; i < Nindex; i++) { + if (IDlist[i] < 2) continue; + int n = DBfound[i]; + myAssert (n > -1, "impossible!"); + fprintf (stdout, "%d : %d -- in %s (%d) : %d \n", i + minID, IDlist[i], DBlist[n], DBfound[i], DBindex[i]); + } + + exit (0); +} + Index: trunk/Ohana/src/dvomerge/src/dvoverify_utils.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvoverify_utils.c (revision 37762) +++ trunk/Ohana/src/dvomerge/src/dvoverify_utils.c (revision 37807) @@ -238,6 +238,6 @@ } // lensingOffset needs to be in range for each object - lensingOffsetOK &= (catalog.average[i].lensingOffset < catalog.Nlensing); - if (VERBOSE && !(catalog.average[i].lensingOffset < catalog.Nlensing)) { + lensingOffsetOK &= (!catalog.Nlensing || (catalog.average[i].lensingOffset < catalog.Nlensing)); + if (VERBOSE && catalog.Nlensing && !(catalog.average[i].lensingOffset < catalog.Nlensing)) { fprintf (stderr, "lensingOffset >= catalog.Nlensing: %d %d %d\n", i, catalog.average[i].lensingOffset, (int) catalog.Nlensing); } Index: trunk/Ohana/src/dvomerge/src/merge_catalogs_new.c =================================================================== --- trunk/Ohana/src/dvomerge/src/merge_catalogs_new.c (revision 37762) +++ trunk/Ohana/src/dvomerge/src/merge_catalogs_new.c (revision 37807) @@ -11,5 +11,5 @@ off_t i, j, offset; - off_t NAVERAGE, NMEASURE, NLENSING, Naverage, Nmeasure, Nlensing, NsecfiltIn, NsecfiltOut, Nm; + off_t NAVERAGE, NMEASURE, NLENSING, NSTARPAR, Naverage, Nmeasure, Nlensing, Nstarpar, NsecfiltIn, NsecfiltOut, Nm; Naverage = output[0].Naverage; @@ -81,4 +81,20 @@ output[0].average[Naverage].Nlensing = Nm; + Nm = 0; + for (j = 0; j < input[0].average[i].Nstarpar; j++) { + offset = input[0].average[i].starparOffset + j; + + output[0].starpar[Nstarpar] = input[0].starpar[offset]; + output[0].starpar[Nstarpar].averef = Naverage; + + Nstarpar ++; + Nm ++; + if (Nstarpar == NSTARPAR) { + NSTARPAR += 1000; + REALLOCATE (output[0].starpar, StarPar, NSTARPAR); + } + } + output[0].average[Naverage].Nstarpar = Nm; + Naverage ++; if (Naverage == NAVERAGE) { @@ -91,7 +107,9 @@ REALLOCATE (output[0].measure, Measure, MAX (Nmeasure, 1)); REALLOCATE (output[0].lensing, Lensing, MAX (Nlensing, 1)); + REALLOCATE (output[0].starpar, StarPar, MAX (Nstarpar, 1)); REALLOCATE (output[0].secfilt, SecFilt, NsecfiltOut*MAX (Naverage, 1)); output[0].Naverage = Naverage; output[0].Nlensing = Nlensing; + output[0].Nstarpar = Nstarpar; output[0].Nsecf_mem = Naverage * NsecfiltOut; Index: trunk/Ohana/src/dvomerge/src/merge_catalogs_old.c =================================================================== --- trunk/Ohana/src/dvomerge/src/merge_catalogs_old.c (revision 37762) +++ trunk/Ohana/src/dvomerge/src/merge_catalogs_old.c (revision 37807) @@ -4,5 +4,5 @@ # define IN_REGION(R,D) ( \ ((D) >= region[0].Dmin) && ((D) < region[0].Dmax) && \ - ((R) >= region[0].Rmin) && ((R) < region[0].Rmax)) + ((R) >= region[0].Rmin) && ((R) < region[0].Rmax)) // merge the input data into the output catalog @@ -14,6 +14,6 @@ double *X1, *Y1, *X2, *Y2; double dX, dY, dR; - off_t *N1, *N2, *next_meas, *next_lens; - off_t Nave, NAVE, Nmeas, NMEAS, Nmatch, Nlens, NLENS; + off_t *N1, *N2, *next_meas, *next_lens, *next_star; + off_t Nave, NAVE, Nmeas, NMEAS, Nmatch, Nlens, NLENS, Nstar, NSTAR; int NsecfiltIn; int NsecfiltOut; @@ -53,4 +53,5 @@ NMEAS = Nmeas = output[0].Nmeasure; NLENS = Nlens = output[0].Nlensing; + NSTAR = Nstar = output[0].Nstarpar; // current max obj ID for this catalog @@ -64,4 +65,6 @@ * coordops.update.c). We use the center of the region (catalog) for crval1,2. */ + + InitCoords (&tcoords, "DEC--ARC"); tcoords.crval1 = 0.5*(region[0].Rmin + region[0].Rmax); if (region[0].Dmax < 90) { @@ -70,11 +73,5 @@ tcoords.crval2 = 90.0; } - tcoords.crpix1 = 0; - tcoords.crpix2 = 0; tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; - tcoords.pc1_1 = tcoords.pc2_2 = 1.0; - tcoords.pc1_2 = tcoords.pc2_1 = 0.0; - tcoords.Npolyterms = 1; - strcpy (tcoords.ctype, "DEC--ARC"); if (VERBOSE) fprintf (stderr, "merging %s into %s\n", input[0].filename, output[0].filename); @@ -120,7 +117,9 @@ next_meas = init_measure_links (output[0].average, Nave, output[0].measure, Nmeas); next_lens = init_lensing_links (output[0].average, Nave, output[0].lensing, Nlens); + next_star = init_starpar_links (output[0].average, Nave, output[0].starpar, Nstar); } else { next_meas = build_measure_links (output[0].average, Nave, output[0].measure, Nmeas); next_lens = build_lensing_links (output[0].average, Nave, output[0].lensing, Nlens); + next_star = build_starpar_links (output[0].average, Nave, output[0].starpar, Nstar); } @@ -184,8 +183,13 @@ REALLOCATE (output[0].measure, Measure, NMEAS); } - if (Nlens + input[0].average[N].Nmeasure >= NLENS) { - NLENS = Nlens + input[0].average[N].Nmeasure + 1000; + if (Nlens + input[0].average[N].Nlensing >= NLENS) { + NLENS = Nlens + input[0].average[N].Nlensing + 1000; REALLOCATE (next_lens, off_t, NLENS); REALLOCATE (output[0].lensing, Lensing, NLENS); + } + if (Nstar + input[0].average[N].Nstarpar >= NSTAR) { + NSTAR = Nstar + input[0].average[N].Nstarpar + 1000; + REALLOCATE (next_star, off_t, NSTAR); + REALLOCATE (output[0].starpar, StarPar, NSTAR); } @@ -274,4 +278,22 @@ } + // if lensing measurements exist, add them too + if (output[0].starpar) { + for (Nin = 0; Nin < input[0].average[N].Nstarpar; Nin++) { + /* add to end of lensing list */ + add_star_link (&output[0].average[n], next_star, Nstar, NSTAR); + + // set the new starpar + off_t staroff = input[0].average[N].starparOffset + Nin; + output[0].starpar[Nstar] = input[0].starpar[staroff]; + + output[0].starpar[Nstar].averef = n; + output[0].starpar[Nstar].objID = output[0].average[n].objID; + output[0].starpar[Nstar].catID = output[0].catID; + output[0].average[n].Nstarpar ++; + Nstar ++; + } + } + // update the average properties to reflect the incoming entries: // if the original value is NAN but the input value is not, accept the input: @@ -319,4 +341,9 @@ REALLOCATE (output[0].lensing, Lensing, NLENS); } + if (Nstar + input[0].average[N].Nstarpar >= NSTAR) { + NSTAR = Nstar + input[0].average[N].Nstarpar + 1000; + REALLOCATE (next_star, off_t, NSTAR); + REALLOCATE (output[0].starpar, StarPar, NSTAR); + } if (Nave >= NAVE) { NAVE = Nave + 1000; @@ -326,5 +353,8 @@ if (input[0].found_t[N] >= 0) continue; - if (!IN_REGION (input[0].average[N].R, input[0].average[N].D)) continue; + + // if we are using MATCHED_TABLES, we are going to leave the edge cases in their + // source catalog, even if they have leaked beyond the edge + if (!MATCHED_TABLES && !IN_REGION (input[0].average[N].R, input[0].average[N].D)) continue; // XXX should we accept the input measurements for these fields? @@ -413,4 +443,30 @@ } + /** add starpar for this input average object **/ + if (output[0].starpar) { + output[0].average[Nave].starparOffset = Nstar; + for (Nin = 0; Nin < input[0].average[N].Nstarpar; Nin ++) { + // supply the starpar values from this detection + off_t staroff = input[0].average[N].starparOffset + Nin; + output[0].starpar[Nstar] = input[0].starpar[staroff]; + + // the following starpar elements cannot be set until here: + output[0].starpar[Nstar].averef = Nave; + output[0].starpar[Nstar].objID = output[0].average[Nave].objID; + output[0].starpar[Nstar].catID = output[0].catID; + + // as we add starpar, update Nstarpar to match + output[0].average[Nave].Nstarpar ++; + + /* we set next[Nstar] to -1 here, and update correctly below */ + next_star[Nstar] = -1; + Nstar ++; + } + int Ngroup = input[0].average[N].Nstarpar; + for (j = 0; j < Ngroup - 1; j++) { + next_star[Nstar - Ngroup + j] = Nstar - Ngroup + j + 1; + } + } + Nave ++; } @@ -421,4 +477,5 @@ REALLOCATE (output[0].measure, Measure, Nmeas); REALLOCATE (output[0].lensing, Lensing, Nlens); + REALLOCATE (output[0].starpar, StarPar, Nstar); # define NOSORT 0 @@ -429,4 +486,5 @@ output[0].measure = sort_measure (output[0].average, Nave, output[0].measure, Nmeas, next_meas); output[0].lensing = sort_lensing (output[0].average, Nave, output[0].lensing, Nlens, next_lens); + output[0].starpar = sort_starpar (output[0].average, Nave, output[0].starpar, Nstar, next_star); } @@ -445,4 +503,5 @@ output[0].Nmeasure = Nmeas; output[0].Nlensing = Nlens; + output[0].Nstarpar = Nstar; output[0].Nsecf_mem = Nave*NsecfiltOut; if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas, Nlens: "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT", ("OFF_T_FMT" matches)\n", Nstars, Nave, Nmeas, Nlens, Nmatch); @@ -450,4 +509,5 @@ free (next_meas); free (next_lens); + free (next_star); free (X2); Index: trunk/Ohana/src/dvomerge/src/myIndex.c =================================================================== --- trunk/Ohana/src/dvomerge/src/myIndex.c (revision 37807) +++ trunk/Ohana/src/dvomerge/src/myIndex.c (revision 37807) @@ -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]; +}