Index: trunk/Ohana/src/libdvo/Makefile =================================================================== --- trunk/Ohana/src/libdvo/Makefile (revision 38970) +++ trunk/Ohana/src/libdvo/Makefile (revision 38986) @@ -6,4 +6,5 @@ HOME = $(ROOT)/src/libdvo AUTO = $(ROOT)/src/libautocode +BIN = $(HOME)/bin LIB = $(HOME)/lib SRC = $(HOME)/src @@ -22,5 +23,15 @@ FULL_LDFLAGS = $(BASE_LDFLAGS) -lFITS -lohana -default: install +TEST_CFLAGS = $(BASE_CFLAGS) +TEST_CPPFLAGS = $(BASE_CPPFLAGS) +TEST_LDFLAGS = $(BASE_LDFLAGS) -ldvo -lFITS -lohana -ltap_ohana + +TESTPROG = tap_AstromOffsetMapIO +$(TESTPROG) : % : $(TESTBIN)/% +test: $(TESTPROG) + for i in $(TESTPROG); do $(TESTBIN)/$$i 2>&1 | $(TESTHARNESS); done +test.verbose: $(TESTPROG) + for i in $(TESTPROG); do $(TESTBIN)/$$i 2>&1 |& $(TESTHARNESS) -v; done + install: $(DESTLIB)/libdvo.a $(DESTLIB)/libdvo.$(DLLTYPE) libdvo: $(LIB)/libdvo.$(ARCH).a $(LIB)/libdvo.$(ARCH).$(DLLTYPE) @@ -183,10 +194,3 @@ cd $(AUTO) && $(MAKE) clean -TESTPROG = -$(TESTPROG) : % : $(TESTBIN)/% $(OBJS) -test: $(TESTPROG) - for i in $(TESTPROG); do $(TESTBIN)/$$i 2>&1 | $(TESTHARNESS); done -test.verbose: $(TESTPROG) - for i in $(TESTPROG); do $(TESTBIN)/$$i 2>&1 |& $(TESTHARNESS) -v; done - -.PRECIOUS: $(ASRC)/%.c +# .PRECIOUS: $(ASRC)/%.c Index: trunk/Ohana/src/libdvo/include/dvo.h =================================================================== --- trunk/Ohana/src/libdvo/include/dvo.h (revision 38970) +++ trunk/Ohana/src/libdvo/include/dvo.h (revision 38986) @@ -171,7 +171,9 @@ ID_PROPER = 0x00000010, // star with large proper motion ID_TRANSIENT = 0x00000020, // identified as a non-periodic (stationary) transient - ID_VARIABLE = 0x00000040, // identified as a period variable + ID_VARIABLE = 0x00000040, // identified as a periodic variable ID_ASTEROID = 0x00000080, // identified with a known solar-system object (asteroid or other) - // bits 0x00000100 - 0x00008000 are currently unused + ID_STACK_ASTROM = 0x00000100, // stack position used for astrometry + + // bits 0x00000200 - 0x00008000 are currently unused ID_STAR_FIT_AVE = 0x00010000, // average position fitted ID_STAR_FIT_PM = 0x00020000, // proper motion fitted @@ -813,5 +815,4 @@ int sorted; /* is measure table average-sorted? (NOTE this is an int only because gfits_scan %t requires it) */ - /* pointers for data manipulation */ off_t *found_t; @@ -819,8 +820,4 @@ char *measureRank; -// off_t *image_t; -// off_t *mosaic_t; -// float *X_t; -// float *Y_t; int *nOwn_t; // relastro uses this to count owned detections per object @@ -840,4 +837,6 @@ int gfits_db_free PROTO((FITS_DB *db)); +int gfits_fread_uncompressed (Catalog *catalog, FTable *ftable, char *nativeOrder, char VERBOSE); + char *libdvo_version (void); @@ -903,4 +902,6 @@ float PhotZeroPoint (Measure *measure, Average *average, SecFilt *secfilt); +int PhotSecfiltFlags (PhotCode *code, Average *average, SecFilt *secfilt); +int PhotNcode (PhotCode *code, Average *average, SecFilt *secfilt); int PhotNphot (PhotCode *code, Average *average, SecFilt *secfilt, dvoMagClassType class, dvoMagSourceType source); float PhotMstdev (PhotCode *code, Average *average, SecFilt *secfilt, dvoMagClassType class, dvoMagSourceType source); @@ -969,6 +970,6 @@ void dvo_catalog_test (Catalog *catalog, int halt); -int dvo_catalog_backup (Catalog *catalog, int primary); -int dvo_catalog_unlink_backup (Catalog *catalog, int primary); +int dvo_catalog_backup (Catalog *catalog, char *suffix, int primary); +int dvo_catalog_unlink_backup (Catalog *catalog, char *suffix, int primary); /* catmode-specific APIs */ Index: trunk/Ohana/src/libdvo/include/dvodb.h =================================================================== --- trunk/Ohana/src/libdvo/include/dvodb.h (revision 38970) +++ trunk/Ohana/src/libdvo/include/dvodb.h (revision 38986) @@ -356,4 +356,5 @@ IMAGE_FWHM_MEDIAN, IMAGE_EXPTIME, + IMAGE_EXPNAME_AS_INT, IMAGE_NSTAR, IMAGE_NCAL, Index: trunk/Ohana/src/libdvo/include/get_graphdata.h =================================================================== --- trunk/Ohana/src/libdvo/include/get_graphdata.h (revision 38970) +++ trunk/Ohana/src/libdvo/include/get_graphdata.h (revision 38986) @@ -11,4 +11,6 @@ double labelPadXm, labelPadYm, labelPadXp, labelPadYp; double padXm, padXp, padYm, padYp; + double fLabelRangeXm, fLabelRangeXp, fLabelRangeYm, fLabelRangeYp; + double fMinorXm, fMinorXp, fMinorYm, fMinorYp; Coords coords; int flipeast, flipnorth; Index: trunk/Ohana/src/libdvo/include/libdvo_astro.h =================================================================== --- trunk/Ohana/src/libdvo/include/libdvo_astro.h (revision 38970) +++ trunk/Ohana/src/libdvo/include/libdvo_astro.h (revision 38986) @@ -56,6 +56,6 @@ float dX; float dY; - float **dXv; - float **dYv; + float *dXv; + float *dYv; unsigned int tableID; unsigned int imageID; @@ -65,4 +65,5 @@ typedef struct { int Nmap; + int NMAP; AstromOffsetMap **map; int *imageIDtoTableSeq; @@ -160,12 +161,14 @@ /* in AstromOffsetMapIO.c */ -AstromOffsetTable *AstromOffsetMapLoad (char *filename, int VERBOSE); +AstromOffsetTable *AstromOffsetMapLoad (char *filename, int Nrows, int VERBOSE); int AstromOffsetMapSave (AstromOffsetTable *table, char *filename); AstromOffsetTable *AstromOffsetMapToTable(AstromOffsetMap_Disk_6x6 *map_disk, off_t Nmap); AstromOffsetMap_Disk_6x6 *AstromOffsetTableToMap(AstromOffsetTable *table, off_t *Nmap); +int AstromOffsetTableSetIDs (AstromOffsetTable *table); +AstromOffsetTable *AstromOffsetMapAppendToTable(AstromOffsetTable *table, AstromOffsetMap_Disk_6x6 *map_disk, off_t Nmap); /* in AstromOffsetMapOps.c */ float AstromOffsetMapValue (AstromOffsetMap *map, float x, float y, int xdir); -int AstromOffsetMapFit (AstromOffsetMap *map, float *x, float *y, float *f, int Npts, int xdir); +int AstromOffsetMapFit (AstromOffsetMap *map, float *x, float *y, float *f, float *df, int Npts, int xdir); /* in AstromOffsetMapUtils.c */ Index: trunk/Ohana/src/libdvo/include/ps1_v5_defs.h =================================================================== --- trunk/Ohana/src/libdvo/include/ps1_v5_defs.h (revision 38970) +++ trunk/Ohana/src/libdvo/include/ps1_v5_defs.h (revision 38986) @@ -23,6 +23,8 @@ StarPar_PS1_V5 *StarParInternalTo_PS1_V5 (StarPar *in, off_t Nvalues); -GalPhot *GalPhot_PS1_V5_ToInternal (GalPhot_PS1_V5 *in, off_t Nvalues); -GalPhot_PS1_V5 *GalPhotInternalTo_PS1_V5 (GalPhot *in, off_t Nvalues); +GalPhot *GalPhot_PS1_V5_R1_ToInternal (GalPhot_PS1_V5_R1 *in, off_t Nvalues); +GalPhot *GalPhot_PS1_V5_R0_ToInternal (GalPhot_PS1_V5_R0 *in, off_t Nvalues); + +GalPhot_PS1_V5_R1 *GalPhotInternalTo_PS1_V5_R1 (GalPhot *in, off_t Nvalues); PhotCode *PhotCode_PS1_V5_To_Internal (PhotCode_PS1_V5 *in, off_t Nvalues); Index: trunk/Ohana/src/libdvo/src/AstromOffsetMapIO.c =================================================================== --- trunk/Ohana/src/libdvo/src/AstromOffsetMapIO.c (revision 38970) +++ trunk/Ohana/src/libdvo/src/AstromOffsetMapIO.c (revision 38986) @@ -1,5 +1,5 @@ # include "dvo.h" -AstromOffsetTable *AstromOffsetMapLoad (char *filename, int VERBOSE) { +AstromOffsetTable *AstromOffsetMapLoad (char *filename, int Nrows, int VERBOSE) { off_t Nmap; @@ -7,6 +7,9 @@ Header header; Matrix matrix; + Header theader; FTable ftable; + ftable.header = &theader; + ftable.buffer = NULL; FILE *f = fopen (filename, "r"); @@ -17,5 +20,4 @@ /* load in table data */ - ftable.header = &theader; if (!gfits_fread_header (f, &header)) { if (VERBOSE) fprintf (stderr, "can't read Astrom Offset Map header\n"); @@ -30,26 +32,75 @@ } - // for now, we only have one flavor (6x6) - - if (!gfits_fread_ftable (f, &ftable, "ASTROM_OFFSET_MAP_DISK_6x6")) { - if (VERBOSE) fprintf (stderr, "can't read Astrom Offset Map table\n"); - gfits_free_header (&header); - gfits_free_matrix (&matrix); - fclose (f); - return (NULL); - } - AstromOffsetMap_Disk_6x6 *map_disk = gfits_table_get_AstromOffsetMap_Disk_6x6 (&ftable, &Nmap, NULL, NULL); - if (!map_disk) { - fprintf (stderr, "ERROR: failed to read Astrom Offset Map\n"); - exit (2); - } - - // AstromOffsetMap_Disk_6x6 *map_disk is an external format stored as an array of maps. - // Convert the disk array of maps to then internal format in a rich structure: - AstromOffsetTable *table = AstromOffsetMapToTable (map_disk, Nmap); + // the output table is not what is stored in the FITS file. read and convert + // from the disk file format + AstromOffsetTable *table = NULL; + + // for now, we only have one flavor (6x6) of disk file format + + // loop over the ftable reading blocks of 100k rows at a time + if (Nrows) { + if (!gfits_find_Xheader (f, ftable.header, "ASTROM_OFFSET_MAP_DISK_6x6")) myAbort ("problem 1"); + int NrowsTotal = ftable.header[0].Naxis[1]; + int Nblocks = (NrowsTotal % Nrows) ? (int) (NrowsTotal / Nrows) + 1 : (NrowsTotal / Nrows); + + // allocate the full array up front. passes through AstromOffsetMapAppendToTable + // check the current allocation and will bump the allocation as needed. + ALLOCATE (table, AstromOffsetTable, 1); + ALLOCATE (table->map, AstromOffsetMap *, NrowsTotal); + table->Nmap = 0; + table->NMAP = NrowsTotal; + + // on each call to gfits_fread_ftable_range, we need to reset the FILE pointer to the + // beginning of the data block. again this seems a little wrong. + off_t diskStart = ftello (f); + + int i; + for (i = 0; i < Nblocks; i++) { + if (!gfits_fread_ftable_range (f, FALSE, FALSE, &ftable, i*Nrows, Nrows)) myAbort ("problem 2"); + + AstromOffsetMap_Disk_6x6 *map_disk = gfits_table_get_AstromOffsetMap_Disk_6x6 (&ftable, &Nmap, NULL, NULL); + if (!map_disk) myAbort ("ERROR: failed to read Astrom Offset Map\n"); + + table = AstromOffsetMapAppendToTable (table, map_disk, Nmap); + + // XXX this is kind of hack-ish: gfits_fread_ftable_range modifies Naxis[1] to the + // number of rows actually read, but this makes a subsequent read invalid. I'm + // resetting the value of Naxis[1] before the next read to deal with this problem. + // But is there a more consistent way to address this? + + // Answer: the FTable structure could include information about the table + // representation on disk, including start,Nrows values. I'm not going to mess + // around with that just now. + + ftable.header[0].Naxis[1] = NrowsTotal; + fseeko (f, diskStart, SEEK_SET); + } + AstromOffsetTableSetIDs (table); + gfits_free_header (ftable.header); + gfits_free_table (&ftable); + } else { + + // convert the blocks into the equivalent table entries and append + if (!gfits_fread_ftable (f, &ftable, "ASTROM_OFFSET_MAP_DISK_6x6")) { + if (VERBOSE) fprintf (stderr, "can't read Astrom Offset Map table\n"); + gfits_free_header (&header); + gfits_free_matrix (&matrix); + fclose (f); + return (NULL); + } + AstromOffsetMap_Disk_6x6 *map_disk = gfits_table_get_AstromOffsetMap_Disk_6x6 (&ftable, &Nmap, NULL, NULL); + if (!map_disk) { + fprintf (stderr, "ERROR: failed to read Astrom Offset Map\n"); + exit (2); + } + + // AstromOffsetMap_Disk_6x6 *map_disk is an external format stored as an array of maps. + // Convert the disk array of maps to then internal format in a rich structure: + table = AstromOffsetMapToTable (map_disk, Nmap); + gfits_free_table (&ftable); + } gfits_free_header (&theader); gfits_free_matrix (&matrix); - gfits_free_table (&ftable); return (table); @@ -135,24 +186,24 @@ table->map[i][0].Nx = map_disk[i].Nx; table->map[i][0].Ny = map_disk[i].Ny; - table->map[i][0].tableID = map_disk[i].tableID; - table->map[i][0].imageID = map_disk[i].imageID; + table->map[i][0].tableID = map_disk[i].tableID; + table->map[i][0].imageID = map_disk[i].imageID; // STORE THESE VALUES? - table->map[i][0].dY = map_disk[i].dX; - table->map[i][0].dX = map_disk[i].dY; + table->map[i][0].dX = map_disk[i].dX; + table->map[i][0].dY = map_disk[i].dY; // since this was on disk, we obviously keep it table->map[i][0].keep = TRUE; - ALLOCATE (table->map[i][0].dXv, float *, map_disk[i].Nx); - ALLOCATE (table->map[i][0].dYv, float *, map_disk[i].Nx); - - for (j = 0; j < map_disk[i].Nx; j++) { - ALLOCATE (table->map[i][0].dXv[j], float, map_disk[i].Ny); - ALLOCATE (table->map[i][0].dYv[j], float, map_disk[i].Ny); - - for (k = 0; k < map_disk[i].Ny; k++) { - table->map[i][0].dXv[j][k] = map_disk[i].dXv[j][k]; - table->map[i][0].dYv[j][k] = map_disk[i].dYv[j][k]; + int Nx = table->map[i][0].Nx; + int Ny = table->map[i][0].Ny; + + ALLOCATE (table->map[i][0].dXv, float, Nx*Ny); + ALLOCATE (table->map[i][0].dYv, float, Nx*Ny); + + for (j = 0; j < Nx; j++) { + for (k = 0; k < Ny; k++) { + table->map[i][0].dXv[j + k*Nx] = map_disk[i].dXv[j][k]; + table->map[i][0].dYv[j + k*Nx] = map_disk[i].dYv[j][k]; } } @@ -182,10 +233,11 @@ map_disk[Ndisk].dY = table->map[i][0].dY; - - - for (j = 0; j < map_disk[Ndisk].Nx; j++) { - for (k = 0; k < map_disk[Ndisk].Ny; k++) { - map_disk[Ndisk].dXv[j][k] = table->map[i][0].dXv[j][k]; - map_disk[Ndisk].dYv[j][k] = table->map[i][0].dYv[j][k]; + int Nx = table->map[i][0].Nx; + int Ny = table->map[i][0].Ny; + + for (j = 0; j < Nx; j++) { + for (k = 0; k < Ny; k++) { + map_disk[Ndisk].dXv[j][k] = table->map[i][0].dXv[j + k*Nx]; + map_disk[Ndisk].dYv[j][k] = table->map[i][0].dYv[j + k*Nx]; } for (k = map_disk[Ndisk].Ny; k < 6; k++) { @@ -206,2 +258,82 @@ } +AstromOffsetTable *AstromOffsetMapAppendToTable(AstromOffsetTable *table, AstromOffsetMap_Disk_6x6 *map_disk, off_t Nmap) { + + int i, j, k, Nstart; + + if (table == NULL) { + ALLOCATE (table, AstromOffsetTable, 1); + Nstart = 0; + table->Nmap = Nmap; + ALLOCATE (table->map, AstromOffsetMap *, Nmap); + } else { + Nstart = table->Nmap; + if (table->Nmap + Nmap > table->NMAP) { + table->NMAP = table->Nmap + Nmap; + REALLOCATE (table->map, AstromOffsetMap *, table->NMAP); + } + table->Nmap += Nmap; + } + + // append map values (this allocates just the area needed for each image, not the + // full 6x6, saving some memory while doing the analysis) + for (i = 0; i < Nmap; i++) { + int N = i + Nstart; + ALLOCATE (table->map[N], AstromOffsetMap, 1); + + table->map[N][0].Nx = map_disk[i].Nx; + table->map[N][0].Ny = map_disk[i].Ny; + table->map[N][0].tableID = map_disk[i].tableID; + table->map[N][0].imageID = map_disk[i].imageID; + + // STORE THESE VALUES? + table->map[N][0].dX = map_disk[i].dX; + table->map[N][0].dY = map_disk[i].dY; + + // since this was on disk, we obviously keep it + table->map[N][0].keep = TRUE; + + int Nx = table->map[N][0].Nx; + int Ny = table->map[N][0].Ny; + + ALLOCATE (table->map[N][0].dXv, float, Nx*Ny); + ALLOCATE (table->map[N][0].dYv, float, Nx*Ny); + + for (j = 0; j < Nx; j++) { + for (k = 0; k < Ny; k++) { + table->map[N][0].dXv[j + k*Nx] = map_disk[i].dXv[j][k]; + table->map[N][0].dYv[j + k*Nx] = map_disk[i].dYv[j][k]; + } + } + } + return table; +} + +int AstromOffsetTableSetIDs (AstromOffsetTable *table) { + + int i; + + // find the max value of imageID + int MaxTableID = 0; + int MaxImageID = 0; + for (i = 0; i < table->Nmap; i++) { + MaxTableID = MAX(table->map[i][0].tableID, MaxTableID); + MaxImageID = MAX(table->map[i][0].imageID, MaxImageID); + } + table->MaxTableID = MaxTableID; + table->MaxImageID = MaxImageID; + + // generate the index and init values to -1 + ALLOCATE (table->imageIDtoTableSeq, int, MaxImageID + 1); + for (i = 0; i <= MaxImageID; i++) { + table->imageIDtoTableSeq[i] = -1; + } + + // assign the ID values + for (i = 0; i < table->Nmap; i++) { + int ImageID = table->map[i][0].imageID; + myAssert (table->imageIDtoTableSeq[ImageID] == -1, "oops, duplicate image IDs"); + table->imageIDtoTableSeq[ImageID] = i; + } + return TRUE; +} Index: trunk/Ohana/src/libdvo/src/AstromOffsetMapOps.c =================================================================== --- trunk/Ohana/src/libdvo/src/AstromOffsetMapOps.c (revision 38970) +++ trunk/Ohana/src/libdvo/src/AstromOffsetMapOps.c (revision 38986) @@ -6,5 +6,5 @@ int dump_map_data (float *x, float *y, float *f, int Npts, char *filename); int AstromOffsetMapFit_Chisq (AstromOffsetMap *map, float *x, float *y, float *f, int Npts, int xdir); -int AstromOffsetMapFit_Median (AstromOffsetMap *map, float *x, float *y, float *f, int Npts, int xdir); +int AstromOffsetMapFit_Mean (AstromOffsetMap *map, float *x, float *y, float *f, float *df, int Npts, int xdir); float AstromOffsetMapValue (AstromOffsetMap *map, float x, float y, int xdir) { @@ -12,30 +12,33 @@ // given x,y find the offset (x-direction) at that location from the map - float **V = xdir ? map->dXv : map->dYv; + int Nx = map->Nx; + int Ny = map->Ny; + + float *V = xdir ? map->dXv : map->dYv; // if there is no spatial information, the value at the point is the value in the map - if ((map->Nx == 1) && (map->Ny == 1)) return V[0][0]; + if ((Nx == 1) && (Ny == 1)) return V[0]; // if there is no spatial information in 1D, the interpolation is 1D - if (map->Nx == 1) { + if (Nx == 1) { float ymo = y * map->dY - 0.5; int ymi = floor(ymo); - ymi = MAX(0,MIN(map->Ny - 2, ymi)); // force range of ymi to be 0,Ny-2 (Ny must be > 1) + ymi = MAX(0,MIN(Ny - 2, ymi)); // force range of ymi to be 0,Ny-2 (Ny must be > 1) float ymf = ymo - ymi; - float value = ymf * V[0][ymi+1] + (1.0 - ymf) * V[0][ymi]; + float value = ymf * V[ymi+1] + (1.0 - ymf) * V[ymi]; return value; } // if there is no spatial information in 1D, the interpolation is 1D - if (map->Ny == 1) { + if (Ny == 1) { float xmo = x * map->dX - 0.5; int xmi = floor(xmo); - xmi = MAX(0,MIN(map->Nx - 2, xmi)); // force range of ymi to be 0,Nx-2 (Nx must be > 1) + xmi = MAX(0,MIN(Nx - 2, xmi)); // force range of ymi to be 0,Nx-2 (Nx must be > 1) float xmf = xmo - xmi; - float value = xmf * V[xmi+1][0] + (1.0 - xmf) * V[xmi][0]; + float value = xmf * V[xmi+1] + (1.0 - xmf) * V[xmi]; return value; } @@ -45,16 +48,16 @@ float xmo = x * map->dX - 0.5; int xmi = floor(xmo); - xmi = MAX(0,MIN(map->Nx - 2, xmi)); // force range of ymi to be 0,Nx-1 + xmi = MAX(0,MIN(Nx - 2, xmi)); // force range of ymi to be 0,Nx-1 float xmf = xmo - xmi; float ymo = y * map->dY - 0.5; int ymi = floor(ymo); - ymi = MAX(0,MIN(map->Ny - 2, ymi)); // force range of ymi to be 0,Ny-2 + ymi = MAX(0,MIN(Ny - 2, ymi)); // force range of ymi to be 0,Ny-2 float ymf = ymo - ymi; - float V00 = V[xmi+0][ymi+0]; - float V01 = V[xmi+0][ymi+1]; - float V10 = V[xmi+1][ymi+0]; - float V11 = V[xmi+1][ymi+1]; + float V00 = V[xmi+0 + (ymi+0)*Nx]; + float V01 = V[xmi+0 + (ymi+1)*Nx]; + float V10 = V[xmi+1 + (ymi+0)*Nx]; + float V11 = V[xmi+1 + (ymi+1)*Nx]; float Vx0 = V10*xmf + V00*(1.0 - xmf); @@ -65,7 +68,7 @@ } -int AstromOffsetMapFit (AstromOffsetMap *map, float *x, float *y, float *f, int Npts, int xdir) { - - int status = AstromOffsetMapFit_Median (map, x, y, f, Npts, xdir); +int AstromOffsetMapFit (AstromOffsetMap *map, float *x, float *y, float *f, float *df, int Npts, int xdir) { + + int status = AstromOffsetMapFit_Mean (map, x, y, f, df, Npts, xdir); return status; } @@ -78,5 +81,5 @@ // choose to map direction: - float **V = xdir ? map->dXv : map->dYv; + float *Vptr = xdir ? map->dXv : map->dYv; // special cases: @@ -88,5 +91,5 @@ vstats_getstats_f (f, NULL, NULL, Npts, &stats); - V[0][0] = stats.mean; + Vptr[0] = stats.mean; return TRUE; } @@ -352,5 +355,5 @@ for (iy = 0; iy < Ny; iy++) { int I = ix + Nx * iy; - V[ix][iy] = B[I][0]; + Vptr[I] = B[I][0]; // error[ix][iy] = sqrt(A[I][I]); } @@ -371,10 +374,10 @@ // given (x,y),value vector sets, choose the map that minimizes the difference between // the values and bilinear interpolation of the map -int AstromOffsetMapFit_Median (AstromOffsetMap *map, float *x, float *y, float *f, int Npts, int xdir) { +int AstromOffsetMapFit_Mean (AstromOffsetMap *map, float *x, float *y, float *f, float *df, int Npts, int xdir) { int i, ix, iy; // choose to map direction: - float **V = xdir ? map->dXv : map->dYv; + float *Vptr = xdir ? map->dXv : map->dYv; // measure clipped median in each bin @@ -386,10 +389,12 @@ int Npix = Nx*Ny; - double **values; + double **values, **dvalues; int *Nvalue; ALLOCATE (Nvalue, int, Npix); ALLOCATE (values, double *, Npix); + ALLOCATE (dvalues, double *, Npix); for (i = 0; i < Npix; i++) { ALLOCATE (values[i], double, Npts); + ALLOCATE (dvalues[i], double, Npts); Nvalue[i] = 0; } @@ -400,4 +405,6 @@ // data value & weight for this point if (!isfinite(f[i])) continue; + if (df && !isfinite(df[i])) continue; + if (df && df[i] == 0.0) continue; // if (mask && (mask[i] & maskValue)) continue; @@ -414,4 +421,5 @@ int N = Nvalue[I]; values[I][N] = f[i]; + dvalues[I][N] = df ? df[i] : 1.0; Nvalue[I] ++; } @@ -421,8 +429,8 @@ int I = ix + Nx * iy; if (Nvalue[I] > 5) { - vstats_getstats (values[I], NULL, NULL, Nvalue[I], &stats); - V[ix][iy] = stats.mean; + vstats_getstats (values[I], dvalues[I], NULL, Nvalue[I], &stats); + Vptr[I] = stats.mean; } else { - V[ix][iy] = NAN; + Vptr[I] = NAN; } } @@ -432,6 +440,8 @@ for (i = 0; i < Npix; i++) { free (values[i]); + free (dvalues[i]); } free (values); + free (dvalues); free (Nvalue); @@ -453,10 +463,9 @@ int ix, iy; - float **Vfix = NULL; - ALLOCATE (Vfix, float *, Nx); + float *Vfix = NULL; + ALLOCATE (Vfix, float, Nx*Ny); for (ix = 0; ix < Nx; ix++) { - ALLOCATE (Vfix[ix], float, Ny); - for (iy = 0; iy < Ny; iy++) { - Vfix[ix][iy] = NAN; + for (iy = 0; iy < Ny; iy++) { + Vfix[ix + iy*Nx] = NAN; } } @@ -468,5 +477,6 @@ for (ix = 0; ix < Nx; ix++) { for (iy = 0; iy < Ny; iy++) { - float value = xdir ? map->dXv[ix][iy] : map->dYv[ix][iy]; + int I = ix + iy*Nx; + float value = xdir ? map->dXv[I] : map->dYv[I]; if (isnan(value)) continue; mean += value; @@ -479,8 +489,8 @@ for (ix = 0; ix < Nx; ix++) { for (iy = 0; iy < Ny; iy++) { - - float value = xdir ? map->dXv[ix][iy] : map->dYv[ix][iy]; + int I = ix + iy*Nx; + float value = xdir ? map->dXv[I] : map->dYv[I]; if (!isnan(value)) { - Vfix[ix][iy] = value; + Vfix[I] = value; continue; } @@ -499,5 +509,6 @@ if (ny < 0) continue; if (ny >= Ny) continue; - float value = xdir ? map->dXv[nx][ny] : map->dYv[nx][ny]; + int I = ix + iy*Nx; + float value = xdir ? map->dXv[I] : map->dYv[I]; if (isnan(value)) continue; @@ -508,5 +519,5 @@ // if there are no valid (non-NAN) local pixels, use global mean meanLocal = (npixLocal > 0.0) ? meanLocal / npixLocal : mean; - Vfix[ix][iy] = meanLocal; + Vfix[I] = meanLocal; } } @@ -515,13 +526,11 @@ for (ix = 0; ix < Ny; ix++) { for (iy = 0; iy < Ny; iy++) { + int I = ix + iy*Nx; if (xdir) { - map->dXv[ix][iy] = Vfix[ix][iy]; + map->dXv[I] = Vfix[I]; } else { - map->dYv[ix][iy] = Vfix[ix][iy]; - } - } - } - for (ix = 0; ix < Nx; ix++) { - free (Vfix[ix]); + map->dYv[I] = Vfix[I]; + } + } } free (Vfix); Index: trunk/Ohana/src/libdvo/src/AstromOffsetMapUtils.c =================================================================== --- trunk/Ohana/src/libdvo/src/AstromOffsetMapUtils.c (revision 38970) +++ trunk/Ohana/src/libdvo/src/AstromOffsetMapUtils.c (revision 38986) @@ -36,15 +36,12 @@ map->keep = TRUE; // output map size - ALLOCATE (map->dXv, float *, map->Nx); - ALLOCATE (map->dYv, float *, map->Nx); - - int j, k; - for (j = 0; j < map->Nx; j++) { - ALLOCATE (map->dXv[j], float, map->Ny); - ALLOCATE (map->dYv[j], float, map->Ny); - - for (k = 0; k < map->Ny; k++) { - map->dXv[j][k] = 0.0; - map->dYv[j][k] = 0.0; + ALLOCATE (map->dXv, float, Nx*Ny); + ALLOCATE (map->dYv, float, Nx*Ny); + + int j, k; + for (j = 0; j < Nx; j++) { + for (k = 0; k < Ny; k++) { + map->dXv[j + k*Nx] = 0.0; + map->dYv[j + k*Nx] = 0.0; } } @@ -56,10 +53,4 @@ if (!map) return; - - int j; - for (j = 0; j < map->Nx; j++) { - if (map->dXv) free (map->dXv[j]); - if (map->dYv) free (map->dYv[j]); - } free (map->dXv); @@ -76,8 +67,4 @@ // rather than try to figure out how to resize, and free/allocate, lets just free the old arrays and make new ones - for (j = 0; j < map->Nx; j++) { - if (map->dXv) free (map->dXv[j]); - if (map->dYv) free (map->dYv[j]); - } free (map->dXv); free (map->dYv); @@ -90,14 +77,11 @@ map->keep = TRUE; // output map size - ALLOCATE (map->dXv, float *, map->Nx); - ALLOCATE (map->dYv, float *, map->Nx); - - for (j = 0; j < map->Nx; j++) { - ALLOCATE (map->dXv[j], float, map->Ny); - ALLOCATE (map->dYv[j], float, map->Ny); - - for (k = 0; k < map->Ny; k++) { - map->dXv[j][k] = 0.0; - map->dYv[j][k] = 0.0; + ALLOCATE (map->dXv, float, Nx*Ny); + ALLOCATE (map->dYv, float, Nx*Ny); + + for (j = 0; j < Nx; j++) { + for (k = 0; k < Ny; k++) { + map->dXv[j + k*Nx] = 0.0; + map->dYv[j + k*Nx] = 0.0; } } @@ -116,9 +100,12 @@ tgt->tableID = map->tableID; - int j, k; - for (j = 0; j < tgt->Nx; j++) { - for (k = 0; k < map->Ny; k++) { - tgt->dXv[j][k] = map->dXv[j][k]; - tgt->dYv[j][k] = map->dYv[j][k]; + int Nx = map->Nx; + int Ny = map->Ny; + + int j, k; + for (j = 0; j < Nx; j++) { + for (k = 0; k < Ny; k++) { + tgt->dXv[j + k*Nx] = map->dXv[j + k*Nx]; + tgt->dYv[j + k*Nx] = map->dYv[j + k*Nx]; } } @@ -139,9 +126,12 @@ tgt->keep = src->keep; - int j, k; - for (j = 0; j < tgt->Nx; j++) { - for (k = 0; k < src->Ny; k++) { - tgt->dXv[j][k] = src->dXv[j][k]; - tgt->dYv[j][k] = src->dYv[j][k]; + int Nx = src->dX; + int Ny = src->dY; + + int j, k; + for (j = 0; j < Nx; j++) { + for (k = 0; k < Ny; k++) { + tgt->dXv[j + k*Nx] = src->dXv[j + k*Nx]; + tgt->dYv[j + k*Nx] = src->dYv[j + k*Nx]; } } @@ -153,10 +143,6 @@ if (!table) return; - int i, j; + int i; for (i = 0; i < table->Nmap; i++) { - for (j = 0; j < table->map[i][0].Nx; j++) { - free (table->map[i][0].dXv[j]); - free (table->map[i][0].dYv[j]); - } if (!table->map[i]) continue; FREE (table->map[i][0].dXv); @@ -260,17 +246,20 @@ } + int Nx = map->Nx; + int Ny = map->Ny; + fprintf (f, "imageID: %d, tableID: %d (dX: %f, dY: %f), keep: %d\n", map->imageID, map->tableID, map->dX, map->dY, map->keep); int ix, iy; fprintf (f, "dXv map:\n"); - for (ix = 0; ix < map->Nx; ix++) { - for (iy = 0; iy < map->Ny; iy++) { - fprintf (f, "%9.5f ", map->dXv[ix][iy]); + for (ix = 0; ix < Nx; ix++) { + for (iy = 0; iy < Ny; iy++) { + fprintf (f, "%9.5f ", map->dXv[ix + iy*Nx]); } fprintf (f, "\n"); } fprintf (f, "dYv map:\n"); - for (ix = 0; ix < map->Nx; ix++) { - for (iy = 0; iy < map->Ny; iy++) { - fprintf (f, "%9.5f ", map->dYv[ix][iy]); + for (ix = 0; ix < Nx; ix++) { + for (iy = 0; iy < Ny; iy++) { + fprintf (f, "%9.5f ", map->dYv[ix + iy*Nx]); } fprintf (f, "\n"); Index: trunk/Ohana/src/libdvo/src/LoadImages.c =================================================================== --- trunk/Ohana/src/libdvo/src/LoadImages.c (revision 38970) +++ trunk/Ohana/src/libdvo/src/LoadImages.c (revision 38986) @@ -85,5 +85,5 @@ if (table) AstromOffsetTableFree(table); - table = AstromOffsetMapLoad (mapfile, FALSE); + table = AstromOffsetMapLoad (mapfile, 100000, FALSE); // assign images.coords.offsetMap -> table->map[i] Index: trunk/Ohana/src/libdvo/src/dbExtractImages.c =================================================================== --- trunk/Ohana/src/libdvo/src/dbExtractImages.c (revision 38970) +++ trunk/Ohana/src/libdvo/src/dbExtractImages.c (revision 38986) @@ -233,4 +233,14 @@ value.Flt = image[N].exptime; break; + + case IMAGE_EXPNAME_AS_INT: + // XXX NOTE this does not handle gpc2 images + if ((image[N].name[0] == 'o') && (image[N].name[5] == 'g') && (image[N].name[10] == 'o')) { + int mjd = atoi(&image[N].name[1]); + int Nexp = atoi(&image[N].name[6]); + value.Int = mjd * 10000 + Nexp; + } + break; + case IMAGE_SIDTIME: value.Flt = image[N].sidtime; Index: trunk/Ohana/src/libdvo/src/dbExtractMeasures.c =================================================================== --- trunk/Ohana/src/libdvo/src/dbExtractMeasures.c (revision 38970) +++ trunk/Ohana/src/libdvo/src/dbExtractMeasures.c (revision 38986) @@ -241,12 +241,15 @@ break; - case MAG_OPTION_STDEV: - case MAG_OPTION_CHISQ: - case MAG_OPTION_MIN: - case MAG_OPTION_MAX: - case MAG_OPTION_NCODE: - case MAG_OPTION_NPHOT: - case MAG_OPTION_UC_DIST: - case MAG_OPTION_FLAGS: + // the following are all values which are come from the secfit table + case MAG_OPTION_STDEV: { value.Flt = PhotMstdev (equiv, average, secfilt, field->magClass, field->magSource); break; } + case MAG_OPTION_CHISQ: { value.Flt = PhotXm (equiv, average, secfilt); break; } + case MAG_OPTION_MIN: { value.Flt = PhotMmin (equiv, average, secfilt); break;} + case MAG_OPTION_MAX: { value.Flt = PhotMmax (equiv, average, secfilt); break;} + case MAG_OPTION_NCODE: { value.Int = PhotNcode (equiv, average, secfilt); break; } + case MAG_OPTION_NPHOT: { value.Int = PhotNphot (equiv, average, secfilt, field->magClass, field->magSource); break; } + case MAG_OPTION_UC_DIST: { value.Flt = PhotUCdist (equiv, average, secfilt); break; } + case MAG_OPTION_FLAGS: { value.Int = PhotSecfiltFlags (equiv, average, secfilt); break; } + + // the following values come from the mean lensobj table case MAG_OPTION_X11_SM_OBJ: case MAG_OPTION_X12_SM_OBJ: Index: trunk/Ohana/src/libdvo/src/dbFields.c =================================================================== --- trunk/Ohana/src/libdvo/src/dbFields.c (revision 38970) +++ trunk/Ohana/src/libdvo/src/dbFields.c (revision 38986) @@ -356,4 +356,5 @@ dvoMagSourceType source = GetMagSource (word); if (source != MAG_SRC_NONE) { + // if mySource is already set, then we have two 'source' elements (an error) if (mySource != MAG_SRC_NONE) { fprintf (stderr, "invalid selection %s in %s\n", word, fieldName); free (word); return FALSE; } mySource = source; @@ -364,4 +365,5 @@ dvoMagLevelType level = GetMagLevel (word); if (level != MAG_LEVEL_NONE) { + // if myLevel is already set, then we have two 'level' elements (an error) if (myLevel != MAG_LEVEL_NONE) { fprintf (stderr, "invalid selection %s in %s\n", word, fieldName); free (word); return FALSE; } myLevel = level; @@ -372,4 +374,5 @@ dvoMagClassType class = GetMagClass (word); if (class != MAG_CLASS_NONE) { + // if myClass is already set, then we have two 'class' elements (an error) if (myClass != MAG_CLASS_NONE) { fprintf (stderr, "invalid selection %s in %s\n", word, fieldName); free (word); return FALSE; } myClass = class; @@ -380,36 +383,22 @@ dvoMagOptionType option = GetMagOption (word); if (option != MAG_OPTION_NONE) { - // some combinations are OK: mag + err -> err - if ((myOption == MAG_OPTION_MAG) && (option == MAG_OPTION_ERR)) { + // if myClass is already set, then we have two 'class' elements + // this is an error unless the request was of the form mag:value + // if so, then code->type will be PHOT_MAG + + if (code->type == PHOT_MAG) { + // mag:err -> magerr + if ((myOption == MAG_OPTION_MAG) && (option == MAG_OPTION_ERR)) { + myOption = MAG_OPTION_ERR; + ptr = skipword (ptr); + continue; + } + // flux:err -> fluxerr + if ((myOption == MAG_OPTION_FLUX) && (option == MAG_OPTION_ERR)) { + myOption = MAG_OPTION_FLUX_ERR; + ptr = skipword (ptr); + continue; + } myOption = option; - ptr = skipword (ptr); - continue; - } - // some combinations are OK: err + mag -> err - if ((myOption == MAG_OPTION_ERR) && (option == MAG_OPTION_MAG)) { - ptr = skipword (ptr); - continue; - } - // some combinations are OK: flux + err -> FLUX_ERR - if ((myOption == MAG_OPTION_FLUX) && (option == MAG_OPTION_ERR)) { - myOption = MAG_OPTION_FLUX_ERR; - ptr = skipword (ptr); - continue; - } - // some combinations are OK: err + flux -> FLUX_ERR - if ((myOption == MAG_OPTION_ERR) && (option == MAG_OPTION_FLUX)) { - myOption = MAG_OPTION_FLUX_ERR; - ptr = skipword (ptr); - continue; - } - // some combinations are OK: flux + fluxerr -> FLUX_ERR - if ((myOption == MAG_OPTION_FLUX) && (option == MAG_OPTION_FLUX_ERR)) { - myOption = MAG_OPTION_FLUX_ERR; - ptr = skipword (ptr); - continue; - } - // some combinations are OK: fluxerr + flux -> FLUX_ERR - if ((myOption == MAG_OPTION_FLUX_ERR) && (option == MAG_OPTION_FLUX)) { - myOption = MAG_OPTION_FLUX_ERR; ptr = skipword (ptr); continue; @@ -815,4 +804,5 @@ if (!strcasecmp (fieldName, "photcode" )) ESCAPE (IMAGE_PHOTCODE, OPIHI_INT); if (!strcasecmp (fieldName, "exptime" )) ESCAPE (IMAGE_EXPTIME, OPIHI_FLT); + if (!strcasecmp (fieldName, "expname" )) ESCAPE (IMAGE_EXPNAME_AS_INT, OPIHI_INT); if (!strcasecmp (fieldName, "sidtime" )) ESCAPE (IMAGE_SIDTIME, OPIHI_FLT); Index: trunk/Ohana/src/libdvo/src/dvo_catalog.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_catalog.c (revision 38970) +++ trunk/Ohana/src/libdvo/src/dvo_catalog.c (revision 38986) @@ -200,5 +200,7 @@ average->catID = 0; average->extID = 0; - average->extIDgc = 0; + + average->uRgal = NAN; + average->uDgal = NAN; } @@ -600,11 +602,8 @@ /* pointers for data manipulation */ - // catalog[0].X_t = NULL; - // catalog[0].Y_t = NULL; catalog[0].nOwn_t = NULL; catalog[0].found_t = NULL; - // catalog[0].image_t = NULL; - // catalog[0].mosaic_t = NULL; catalog[0].foundWarp_t = NULL; + catalog[0].measureRank = NULL; } @@ -922,4 +921,6 @@ void dvo_catalog_free (Catalog *catalog) { + if (!catalog) return; + if (catalog[0].catmode == DVO_MODE_SPLIT) { @@ -1044,5 +1045,5 @@ // make a backup of this catalog (including the measure, secfilt, and missing tables as needed) -int dvo_catalog_backup (Catalog *catalog, int primary) { +int dvo_catalog_backup (Catalog *catalog, char *suffix, int primary) { // skip empty cpt files @@ -1052,5 +1053,5 @@ char tmpfilename[DVO_MAX_PATH]; - int status = snprintf (tmpfilename, DVO_MAX_PATH, "%s~", catalog->filename); + int status = snprintf (tmpfilename, DVO_MAX_PATH, "%s%s", catalog->filename, suffix); if (status >= DVO_MAX_PATH) { fprintf (stderr, "path name too long: %s\n", catalog->filename); @@ -1067,4 +1068,18 @@ } + // play it safe: do not overwrite an existing backup file + struct stat fileStats; + status = stat (tmpfilename, &fileStats); + if (!status) { + fprintf (stderr, "ERROR: backup file %s already exists, exiting\n", tmpfilename); + return FALSE; + } + + // some error accessing the file. there is only one acceptable error: file not found + if (status && (errno != ENOENT)) { + perror ("problem with output target"); + return FALSE; + } + status = rename (catalog->filename, tmpfilename); if (status) { @@ -1085,15 +1100,15 @@ if (catalog[0].catmode == DVO_MODE_SPLIT) { if (catalog[0].measure_catalog != NULL) { - if (!dvo_catalog_backup (catalog[0].measure_catalog, FALSE)) { + if (!dvo_catalog_backup (catalog[0].measure_catalog, suffix, FALSE)) { return FALSE; } } if (catalog[0].missing_catalog != NULL) { - if (!dvo_catalog_backup (catalog[0].missing_catalog, FALSE)) { + if (!dvo_catalog_backup (catalog[0].missing_catalog, suffix, FALSE)) { return FALSE; } } if (catalog[0].secfilt_catalog != NULL) { - if (!dvo_catalog_backup (catalog[0].secfilt_catalog, FALSE)) { + if (!dvo_catalog_backup (catalog[0].secfilt_catalog, suffix, FALSE)) { return FALSE; } @@ -1104,5 +1119,5 @@ status = dvo_catalog_lock (catalog[0].lensing_catalog, lockmode); } else { - if (!dvo_catalog_backup (catalog[0].lensing_catalog, FALSE)) { + if (!dvo_catalog_backup (catalog[0].lensing_catalog, suffix, FALSE)) { return FALSE; } @@ -1114,5 +1129,5 @@ status = dvo_catalog_lock (catalog[0].lensobj_catalog, lockmode); } else { - if (!dvo_catalog_backup (catalog[0].lensobj_catalog, FALSE)) { + if (!dvo_catalog_backup (catalog[0].lensobj_catalog, suffix, FALSE)) { return FALSE; } @@ -1124,5 +1139,5 @@ status = dvo_catalog_lock (catalog[0].starpar_catalog, lockmode); } else { - if (!dvo_catalog_backup (catalog[0].starpar_catalog, FALSE)) { + if (!dvo_catalog_backup (catalog[0].starpar_catalog, suffix, FALSE)) { return FALSE; } @@ -1134,5 +1149,5 @@ status = dvo_catalog_lock (catalog[0].galphot_catalog, lockmode); } else { - if (!dvo_catalog_backup (catalog[0].galphot_catalog, FALSE)) { + if (!dvo_catalog_backup (catalog[0].galphot_catalog, suffix, FALSE)) { return FALSE; } @@ -1144,5 +1159,5 @@ // make a backup of this catalog (including the measure, secfilt, and missing tables as needed) -int dvo_catalog_unlink_backup (Catalog *catalog, int primary) { +int dvo_catalog_unlink_backup (Catalog *catalog, char *suffix, int primary) { if (primary && !catalog->Naverage_disk) { @@ -1152,5 +1167,5 @@ char tmpfilename[DVO_MAX_PATH]; - int status = snprintf (tmpfilename, DVO_MAX_PATH, "%s~", catalog->filename); + int status = snprintf (tmpfilename, DVO_MAX_PATH, "%s%s", catalog->filename, suffix); if (status >= DVO_MAX_PATH) { fprintf (stderr, "path name too long: %s\n", catalog->filename); @@ -1168,35 +1183,35 @@ if (catalog[0].catmode == DVO_MODE_SPLIT) { if (catalog[0].measure_catalog != NULL) { - if (!dvo_catalog_unlink_backup (catalog[0].measure_catalog, FALSE)) { + if (!dvo_catalog_unlink_backup (catalog[0].measure_catalog, suffix, FALSE)) { outStatus = FALSE; } } if (catalog[0].missing_catalog != NULL) { - if (!dvo_catalog_unlink_backup (catalog[0].missing_catalog, FALSE)) { + if (!dvo_catalog_unlink_backup (catalog[0].missing_catalog, suffix, FALSE)) { outStatus = FALSE; } } if (catalog[0].secfilt_catalog != NULL) { - if (!dvo_catalog_unlink_backup (catalog[0].secfilt_catalog, FALSE)) { + if (!dvo_catalog_unlink_backup (catalog[0].secfilt_catalog, suffix, FALSE)) { outStatus = FALSE; } } if ((catalog[0].lensing_catalog != NULL) && (catalog[0].Nlensing_disk > 0)) { - if (!dvo_catalog_unlink_backup (catalog[0].lensing_catalog, FALSE)) { + if (!dvo_catalog_unlink_backup (catalog[0].lensing_catalog, suffix, FALSE)) { outStatus = FALSE; } } if ((catalog[0].lensobj_catalog != NULL) && (catalog[0].Nlensobj_disk > 0)) { - if (!dvo_catalog_unlink_backup (catalog[0].lensobj_catalog, FALSE)) { + if (!dvo_catalog_unlink_backup (catalog[0].lensobj_catalog, suffix, FALSE)) { outStatus = FALSE; } } if ((catalog[0].starpar_catalog != NULL) && (catalog[0].Nstarpar_disk > 0)) { - if (!dvo_catalog_unlink_backup (catalog[0].starpar_catalog, FALSE)) { + if (!dvo_catalog_unlink_backup (catalog[0].starpar_catalog, suffix, FALSE)) { outStatus = FALSE; } } if ((catalog[0].galphot_catalog != NULL) && (catalog[0].Ngalphot_disk > 0)) { - if (!dvo_catalog_unlink_backup (catalog[0].galphot_catalog, FALSE)) { + if (!dvo_catalog_unlink_backup (catalog[0].galphot_catalog, suffix, FALSE)) { outStatus = FALSE; } Index: trunk/Ohana/src/libdvo/src/dvo_catalog_split.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_catalog_split.c (revision 38970) +++ trunk/Ohana/src/libdvo/src/dvo_catalog_split.c (revision 38986) @@ -92,4 +92,6 @@ if (OHANA_MEMCHECK) ohana_memcheck (TRUE); + if (!catalog->f) return TRUE; + /* rewind file pointers and truncate (file is still open) */ if (fseeko (catalog->f, 0LL, SEEK_SET)) { @@ -98,4 +100,6 @@ return FALSE; } + + if (OHANA_MEMCHECK) ohana_memcheck_func (TRUE); // write PHU header @@ -118,4 +122,6 @@ gfits_free_matrix (&matrix); # endif + + if (OHANA_MEMCHECK) ohana_memcheck_func (TRUE); FTable *outtable = ftable; @@ -135,5 +141,5 @@ char *compressMode = dvo_catalog_compress_string (catalog->catcompress); // I should test how Ntile affects fpack/funpack and speed - if (!gfits_compress_table (ftable, &cmptable, 1000, compressMode)) { + if (!gfits_compress_table (ftable, &cmptable, 10000, compressMode)) { fprintf (stderr, "compression failure\n"); free (compressMode); @@ -548,5 +554,5 @@ } /* read Average table data : format is irrelevant here */ - if (!gfits_fread_ftable_range (catalog[0].f, FALSE, &ftable, start, Nrows)) { + if (!gfits_fread_ftable_range (catalog[0].f, FALSE, FALSE, &ftable, start, Nrows)) { if (VERBOSE) fprintf (stderr, "can't read table average data"); return (FALSE); @@ -578,5 +584,5 @@ } /* read Secfilt table data : format is irrelevant here */ - if (!gfits_fread_ftable_range (subcat[0].f, FALSE, &ftable, start*catalog[0].Nsecfilt, catalog[0].Naverage*catalog[0].Nsecfilt)) { + if (!gfits_fread_ftable_range (subcat[0].f, FALSE, FALSE, &ftable, start*catalog[0].Nsecfilt, catalog[0].Naverage*catalog[0].Nsecfilt)) { if (VERBOSE) fprintf (stderr, "can't read table measure data"); return (FALSE); @@ -614,5 +620,5 @@ } /* read Measure table data : format is irrelevant here */ - if (!gfits_fread_ftable_range (subcat[0].f, FALSE, &ftable, start, Nrows)) { + if (!gfits_fread_ftable_range (subcat[0].f, FALSE, FALSE, &ftable, start, Nrows)) { if (VERBOSE) fprintf (stderr, "can't read table measure data"); return (FALSE); @@ -645,5 +651,5 @@ } /* read Missing table data : format is irrelevant here */ - if (!gfits_fread_ftable_range (subcat[0].f, FALSE, &ftable, start, Nrows)) { + if (!gfits_fread_ftable_range (subcat[0].f, FALSE, FALSE, &ftable, start, Nrows)) { if (VERBOSE) fprintf (stderr, "can't read table missing data"); return (FALSE); @@ -679,5 +685,5 @@ } /* read Lensing table data : format is irrelevant here */ - if (!gfits_fread_ftable_range (subcat[0].f, FALSE, &ftable, start, Nrows)) { + if (!gfits_fread_ftable_range (subcat[0].f, FALSE, FALSE, &ftable, start, Nrows)) { if (VERBOSE) fprintf (stderr, "can't read table lensing data"); return (FALSE); @@ -710,5 +716,5 @@ } /* read Lensobj table data : format is irrelevant here */ - if (!gfits_fread_ftable_range (subcat[0].f, FALSE, &ftable, start, Nrows)) { + if (!gfits_fread_ftable_range (subcat[0].f, FALSE, FALSE, &ftable, start, Nrows)) { if (VERBOSE) fprintf (stderr, "can't read table lensobj data"); return (FALSE); @@ -741,5 +747,5 @@ } /* read StarPar table data : format is irrelevant here */ - if (!gfits_fread_ftable_range (subcat[0].f, FALSE, &ftable, start, Nrows)) { + if (!gfits_fread_ftable_range (subcat[0].f, FALSE, FALSE, &ftable, start, Nrows)) { if (VERBOSE) fprintf (stderr, "can't read table starpar data"); return (FALSE); @@ -772,5 +778,5 @@ } /* read GalPhot table data : format is irrelevant here */ - if (!gfits_fread_ftable_range (subcat[0].f, FALSE, &ftable, start, Nrows)) { + if (!gfits_fread_ftable_range (subcat[0].f, FALSE, FALSE, &ftable, start, Nrows)) { if (VERBOSE) fprintf (stderr, "can't read table galphot data"); return (FALSE); @@ -842,4 +848,6 @@ } + if (OHANA_MEMCHECK) ohana_memcheck_func (TRUE); + Nmeasure_disk_new = MAX (catalog[0].Nmeasure_disk, catalog[0].Nmeasure + catalog[0].Nmeasure_off); Nmissing_disk_new = MAX (catalog[0].Nmissing_disk, catalog[0].Nmissing + catalog[0].Nmissing_off); @@ -865,4 +873,6 @@ /* in split mode, we can save only part of the data */ + if (OHANA_MEMCHECK) ohana_memcheck_func (TRUE); + /*** Average Table ***/ if ((catalog[0].catflags & DVO_LOAD_AVERAGE) && (catalog[0].average != NULL)) { @@ -890,8 +900,12 @@ } + if (OHANA_MEMCHECK) ohana_memcheck_func (TRUE); + if (!dvo_catalog_save_subcat (catalog, &ftable, start, Nrows, catalog[0].Naverage_disk, Naverage_disk_new)) { fprintf (stderr, "failure writing Average table\n"); goto failure; } + + if (OHANA_MEMCHECK) ohana_memcheck_func (TRUE); gfits_free_header (&header); @@ -922,4 +936,6 @@ assert (catalog[0].Nmeasure_disk >= catalog[0].Nmeasure_off); + if (OHANA_MEMCHECK) ohana_memcheck_func (TRUE); + // if we are going to compress, we need to receive unswapped data -- int swapFromNative = !output_is_compressed (start, Nrows, Nmeasure_disk_new, catalog->catcompress); @@ -930,4 +946,6 @@ goto failure; } + + if (OHANA_MEMCHECK) ohana_memcheck_func (TRUE); // write out Measure table @@ -937,4 +955,6 @@ goto failure; } + if (OHANA_MEMCHECK) ohana_memcheck_func (TRUE); + gfits_free_header (&header); @@ -986,4 +1006,6 @@ int swapFromNative = !output_is_compressed (start, Nrows, Nsecfilt_disk_new, catalog->catcompress); + if (OHANA_MEMCHECK) ohana_memcheck_func (TRUE); + // convert to external table format if (!SecFiltToFtable (&ftable, &secfilt[first], Nrows, catalog[0].catformat, swapFromNative)) { @@ -991,4 +1013,6 @@ goto failure; } + + if (OHANA_MEMCHECK) ohana_memcheck_func (TRUE); // write out SecFilt table @@ -1003,5 +1027,5 @@ /*** Lensing Table ***/ - if ((catalog[0].catflags & DVO_LOAD_LENSING) && (catalog[0].lensing != NULL)) { + if ((catalog[0].catflags & DVO_LOAD_LENSING) && catalog[0].lensing && Nlensing_disk_new) { off_t first = 0; // first row in memory to write @@ -1034,5 +1058,5 @@ /*** Lensobj Table ***/ - if ((catalog[0].catflags & DVO_LOAD_LENSOBJ) && (catalog[0].lensobj != NULL)) { + if ((catalog[0].catflags & DVO_LOAD_LENSOBJ) && catalog[0].lensobj && Nlensobj_disk_new) { off_t first = 0; // first row in memory to write @@ -1065,5 +1089,5 @@ /*** StarPar Table ***/ - if ((catalog[0].catflags & DVO_LOAD_STARPAR) && (catalog[0].starpar != NULL)) { + if ((catalog[0].catflags & DVO_LOAD_STARPAR) && catalog[0].starpar && Nstarpar_disk_new) { off_t first = 0; // first row in memory to write @@ -1096,5 +1120,5 @@ /*** GalPhot Table ***/ - if ((catalog[0].catflags & DVO_LOAD_GALPHOT) && (catalog[0].galphot != NULL)) { + if ((catalog[0].catflags & DVO_LOAD_GALPHOT) && catalog[0].galphot && Ngalphot_disk_new) { off_t first = 0; // first row in memory to write Index: trunk/Ohana/src/libdvo/src/dvo_convert.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert.c (revision 38970) +++ trunk/Ohana/src/libdvo/src/dvo_convert.c (revision 38986) @@ -77,4 +77,5 @@ average = Average_PS1_V4alt_ToInternal (tmpAverage, *Naverage); free (tmpAverage); + ftable[0].buffer = NULL; *format = DVO_FORMAT_PS1_V4; return (average); @@ -91,19 +92,21 @@ average = Average_PS1_V5alt_ToInternal (tmpAverage, *Naverage); free (tmpAverage); + ftable[0].buffer = NULL; *format = DVO_FORMAT_PS1_V5; return (average); } -# define CONVERT_FORMAT(NAME, FORMAT, TYPE) \ - if (!strcmp (extname, NAME)) { \ - Average_##TYPE *tmpAverage; \ - tmpAverage = gfits_table_get_Average_##TYPE (ftable, Naverage, NULL, &nativeBytes); \ - if (!tmpAverage) { \ - fprintf (stderr, "ERROR: failed to read averages\n"); \ - exit (2); \ - } \ +# define CONVERT_FORMAT(NAME, FORMAT, TYPE) \ + if (!strcmp (extname, NAME)) { \ + Average_##TYPE *tmpAverage; \ + tmpAverage = gfits_table_get_Average_##TYPE (ftable, Naverage, NULL, &nativeBytes); \ + if (!tmpAverage) { \ + fprintf (stderr, "ERROR: failed to read averages\n"); \ + exit (2); \ + } \ average = Average_##TYPE##_ToInternal (tmpAverage, *Naverage, primary); \ - free (tmpAverage); \ - *format = DVO_FORMAT_##FORMAT; \ + free (tmpAverage); \ + ftable[0].buffer = NULL; \ + *format = DVO_FORMAT_##FORMAT; \ return (average); } @@ -141,11 +144,11 @@ int AverageToFtable (FTable *ftable, Average *average, off_t Naverage, DVOCatFormat format, SecFilt *primary, int swapFromNative) { -# define FORMAT_CASE(FORMAT, TYPE) \ +# define FORMAT_CASE(FORMAT, TYPE) \ case DVO_FORMAT_##FORMAT: { \ Average_##TYPE *tmpAverage; \ - tmpAverage = AverageInternalTo_##TYPE (average, Naverage, primary); \ - gfits_table_set_Average_##TYPE (ftable, tmpAverage, Naverage, swapFromNative); \ - free (tmpAverage); \ - break; } + tmpAverage = AverageInternalTo_##TYPE (average, Naverage, primary); \ + gfits_table_set_Average_##TYPE (ftable, tmpAverage, Naverage, swapFromNative); \ + free (tmpAverage); \ + break; } /* convert from the internal format */ @@ -195,4 +198,5 @@ if (!strcmp (extname, "DVO_MEASURE_PS1_V4") && (ftable[0].header[0].Naxis[0] == 176)) { fprintf (stderr, "reading alt PS1_V4 format\n"); + myAssert (!nativeBytes, "need to implement optional swap"); Measure_PS1_V4alt *tmpMeasure; tmpMeasure = gfits_table_get_Measure_PS1_V4alt (ftable, Nmeasure, NULL); @@ -204,4 +208,5 @@ measure = Measure_PS1_V4alt_ToInternal (average, tmpMeasure, *Nmeasure); free (tmpMeasure); + ftable[0].buffer = NULL; *format = DVO_FORMAT_PS1_V4; return (measure); @@ -211,4 +216,5 @@ if (!strcmp (extname, "DVO_MEASURE_PS1_V5") && (ftable[0].header[0].Naxis[0] == 232)) { fprintf (stderr, "reading alt PS1_V5 format\n"); + myAssert (!nativeBytes, "need to implement optional swap"); Measure_PS1_V5alt *tmpMeasure; tmpMeasure = gfits_table_get_Measure_PS1_V5alt (ftable, Nmeasure, NULL); @@ -220,19 +226,21 @@ measure = Measure_PS1_V5alt_ToInternal (average, tmpMeasure, *Nmeasure); free (tmpMeasure); + ftable[0].buffer = NULL; *format = DVO_FORMAT_PS1_V5; return (measure); } -# define CONVERT_FORMAT(NAME, FORMAT, TYPE, ABS_COORDS) \ +# define CONVERT_FORMAT(NAME, FORMAT, TYPE, ABS_COORDS) \ if (!strcmp (extname, NAME)) { \ Measure_##TYPE *tmpMeasure; \ - tmpMeasure = gfits_table_get_Measure_##TYPE (ftable, Nmeasure, NULL, &nativeBytes); \ - if (!tmpMeasure) { \ - fprintf (stderr, "ERROR: failed to read measures\n"); \ - exit (2); \ - } \ - myAssert (ABS_COORDS || average, "conversion to internal needs average table"); \ + tmpMeasure = gfits_table_get_Measure_##TYPE (ftable, Nmeasure, NULL, &nativeBytes); \ + if (!tmpMeasure) { \ + fprintf (stderr, "ERROR: failed to read measures\n"); \ + exit (2); \ + } \ + myAssert (ABS_COORDS || average, "conversion to internal needs average table"); \ measure = Measure_##TYPE##_ToInternal (average, tmpMeasure, *Nmeasure); \ free (tmpMeasure); \ + ftable[0].buffer = NULL; \ *format = DVO_FORMAT_##FORMAT; \ return (measure); } @@ -272,12 +280,12 @@ int MeasureToFtable (FTable *ftable, Average *average, Measure *measure, off_t Nmeasure, DVOCatFormat format, int swapFromNative) { -# define FORMAT_CASE(FORMAT, TYPE, ABS_COORDS) \ - case DVO_FORMAT_##FORMAT: { \ - Measure_##TYPE *tmpMeasure; \ - myAssert (ABS_COORDS || average, "conversion from internal needs average table"); \ - tmpMeasure = MeasureInternalTo_##TYPE (average, measure, Nmeasure); \ - gfits_table_set_Measure_##TYPE (ftable, tmpMeasure, Nmeasure, swapFromNative); \ - free (tmpMeasure); \ - break; } +# define FORMAT_CASE(FORMAT, TYPE, ABS_COORDS) \ + case DVO_FORMAT_##FORMAT: { \ + Measure_##TYPE *tmpMeasure; \ + myAssert (ABS_COORDS || average, "conversion from internal needs average table"); \ + tmpMeasure = MeasureInternalTo_##TYPE (average, measure, Nmeasure); \ + gfits_table_set_Measure_##TYPE (ftable, tmpMeasure, Nmeasure, swapFromNative); \ + free (tmpMeasure); \ + break; } /* convert from the internal format */ @@ -322,12 +330,12 @@ } -# define SKIPPING_FORMAT(NAME, FORMAT, TYPE) \ +# define SKIPPING_FORMAT(NAME, FORMAT, TYPE) \ if (!strcmp (extname, NAME)) { \ fprintf (stderr, "ERROR: format %s not defined for missing, skipping\n", NAME); \ - *Nmissing = 0; \ - return NULL; \ - } - -# define CONVERT_FORMAT(NAME, FORMAT, TYPE) \ + *Nmissing = 0; \ + return NULL; \ + } + +# define CONVERT_FORMAT(NAME, FORMAT, TYPE) \ if (!strcmp (extname, NAME)) { \ Missing_##TYPE *tmpMissing; \ @@ -337,6 +345,7 @@ exit (2); \ } \ - missing = Missing_##TYPE##_ToInternal (tmpMissing, *Nmissing); \ + missing = Missing_##TYPE##_ToInternal (tmpMissing, *Nmissing); \ free (tmpMissing); \ + ftable[0].buffer = NULL; \ *format = DVO_FORMAT_##FORMAT; \ return (missing); } @@ -408,4 +417,5 @@ secfilt = SecFilt_PS1_V5alt_ToInternal (tmpSecFilt, *Nsecfilt); free (tmpSecFilt); + ftable[0].buffer = NULL; *format = DVO_FORMAT_PS1_V5; return (secfilt); @@ -413,15 +423,16 @@ -# define CONVERT_FORMAT(NAME, FORMAT, TYPE) \ - if (!strcmp (extname, NAME)) { \ - SecFilt_##TYPE *tmpSecFilt; \ +# define CONVERT_FORMAT(NAME, FORMAT, TYPE) \ + if (!strcmp (extname, NAME)) { \ + SecFilt_##TYPE *tmpSecFilt; \ tmpSecFilt = gfits_table_get_SecFilt_##TYPE (ftable, Nsecfilt, NULL, &nativeBytes); \ - if (!tmpSecFilt) { \ - fprintf (stderr, "ERROR: failed to read secfilts\n"); \ - exit (2); \ - } \ - secfilt = SecFilt_##TYPE##_ToInternal (tmpSecFilt, *Nsecfilt); \ - free (tmpSecFilt); \ - *format = DVO_FORMAT_##FORMAT; \ + if (!tmpSecFilt) { \ + fprintf (stderr, "ERROR: failed to read secfilts\n"); \ + exit (2); \ + } \ + secfilt = SecFilt_##TYPE##_ToInternal (tmpSecFilt, *Nsecfilt); \ + free (tmpSecFilt); \ + ftable[0].buffer = NULL; \ + *format = DVO_FORMAT_##FORMAT; \ return (secfilt); } @@ -459,11 +470,11 @@ int SecFiltToFtable (FTable *ftable, SecFilt *secfilt, off_t Nsecfilt, DVOCatFormat format, int swapFromNative) { -# define FORMAT_CASE(FORMAT, TYPE) \ - case DVO_FORMAT_##FORMAT: { \ - SecFilt_##TYPE *tmpSecFilt; \ - tmpSecFilt = SecFiltInternalTo_##TYPE (secfilt, Nsecfilt); \ - gfits_table_set_SecFilt_##TYPE (ftable, tmpSecFilt, Nsecfilt, swapFromNative); \ - free (tmpSecFilt); \ - break; } +# define FORMAT_CASE(FORMAT, TYPE) \ + case DVO_FORMAT_##FORMAT: { \ + SecFilt_##TYPE *tmpSecFilt; \ + tmpSecFilt = SecFiltInternalTo_##TYPE (secfilt, Nsecfilt); \ + gfits_table_set_SecFilt_##TYPE (ftable, tmpSecFilt, Nsecfilt, swapFromNative); \ + free (tmpSecFilt); \ + break; } /* convert from the internal format */ @@ -517,4 +528,5 @@ lensing = Lensing_PS1_V5_R0_ToInternal (tmpLensing, *Nlensing); free (tmpLensing); + ftable[0].buffer = NULL; *format = DVO_FORMAT_PS1_V5; return (lensing); } @@ -529,4 +541,5 @@ lensing = Lensing_PS1_V5_R1_ToInternal (tmpLensing, *Nlensing); free (tmpLensing); + ftable[0].buffer = NULL; *format = DVO_FORMAT_PS1_V5; return (lensing); } @@ -541,15 +554,16 @@ lensing = Lensing_PS1_V5_R2_ToInternal (tmpLensing, *Nlensing); free (tmpLensing); + ftable[0].buffer = NULL; *format = DVO_FORMAT_PS1_V5; return (lensing); } -# define SKIPPING_FORMAT(NAME, FORMAT, TYPE) \ +# define SKIPPING_FORMAT(NAME, FORMAT, TYPE) \ if (!strcmp (extname, NAME)) { \ fprintf (stderr, "ERROR: format %s not defined for lensing, skipping\n", NAME); \ - *Nlensing = 0; \ - return NULL; \ - } - -# define CONVERT_FORMAT(NAME, FORMAT, TYPE) \ + *Nlensing = 0; \ + return NULL; \ + } + +# define CONVERT_FORMAT(NAME, FORMAT, TYPE) \ if (!strcmp (extname, NAME)) { \ Lensing_##TYPE *tmpLensing; \ @@ -559,6 +573,7 @@ exit (2); \ } \ - lensing = Lensing_##TYPE##_ToInternal (tmpLensing, *Nlensing); \ + lensing = Lensing_##TYPE##_ToInternal (tmpLensing, *Nlensing); \ free (tmpLensing); \ + ftable[0].buffer = NULL; \ *format = DVO_FORMAT_##FORMAT; \ return (lensing); } @@ -599,11 +614,11 @@ int LensingToFtable (FTable *ftable, Lensing *lensing, off_t Nlensing, DVOCatFormat format, int swapFromNative) { -# define FORMAT_CASE(FORMAT, TYPE) \ - case DVO_FORMAT_##FORMAT: { \ - Lensing_##TYPE *tmpLensing; \ - tmpLensing = LensingInternalTo_##TYPE (lensing, Nlensing); \ - gfits_table_set_Lensing_##TYPE (ftable, tmpLensing, Nlensing, swapFromNative); \ - free (tmpLensing); \ - break; } +# define FORMAT_CASE(FORMAT, TYPE) \ + case DVO_FORMAT_##FORMAT: { \ + Lensing_##TYPE *tmpLensing; \ + tmpLensing = LensingInternalTo_##TYPE (lensing, Nlensing); \ + gfits_table_set_Lensing_##TYPE (ftable, tmpLensing, Nlensing, swapFromNative); \ + free (tmpLensing); \ + break; } /* convert from the internal format */ @@ -656,15 +671,16 @@ lensobj = Lensobj_PS1_V5_R0_ToInternal (tmpLensobj, *Nlensobj); free (tmpLensobj); + ftable[0].buffer = NULL; *format = DVO_FORMAT_PS1_V5; return (lensobj); } -# define SKIPPING_FORMAT(NAME, FORMAT, TYPE) \ +# define SKIPPING_FORMAT(NAME, FORMAT, TYPE) \ if (!strcmp (extname, NAME)) { \ fprintf (stderr, "ERROR: format %s not defined for lensobj, skipping\n", NAME); \ - *Nlensobj = 0; \ - return NULL; \ - } - -# define CONVERT_FORMAT(NAME, FORMAT, TYPE) \ + *Nlensobj = 0; \ + return NULL; \ + } + +# define CONVERT_FORMAT(NAME, FORMAT, TYPE) \ if (!strcmp (extname, NAME)) { \ Lensobj_##TYPE *tmpLensobj; \ @@ -674,6 +690,7 @@ exit (2); \ } \ - lensobj = Lensobj_##TYPE##_ToInternal (tmpLensobj, *Nlensobj); \ + lensobj = Lensobj_##TYPE##_ToInternal (tmpLensobj, *Nlensobj); \ free (tmpLensobj); \ + ftable[0].buffer = NULL; \ *format = DVO_FORMAT_##FORMAT; \ return (lensobj); } @@ -714,11 +731,11 @@ int LensobjToFtable (FTable *ftable, Lensobj *lensobj, off_t Nlensobj, DVOCatFormat format, int swapFromNative) { -# define FORMAT_CASE(FORMAT, TYPE) \ - case DVO_FORMAT_##FORMAT: { \ - Lensobj_##TYPE *tmpLensobj; \ - tmpLensobj = LensobjInternalTo_##TYPE (lensobj, Nlensobj); \ - gfits_table_set_Lensobj_##TYPE (ftable, tmpLensobj, Nlensobj, swapFromNative); \ - free (tmpLensobj); \ - break; } +# define FORMAT_CASE(FORMAT, TYPE) \ + case DVO_FORMAT_##FORMAT: { \ + Lensobj_##TYPE *tmpLensobj; \ + tmpLensobj = LensobjInternalTo_##TYPE (lensobj, Nlensobj); \ + gfits_table_set_Lensobj_##TYPE (ftable, tmpLensobj, Nlensobj, swapFromNative); \ + free (tmpLensobj); \ + break; } /* convert from the internal format */ @@ -762,12 +779,12 @@ } -# define SKIPPING_FORMAT(NAME, FORMAT, TYPE) \ +# define SKIPPING_FORMAT(NAME, FORMAT, TYPE) \ if (!strcmp (extname, NAME)) { \ fprintf (stderr, "ERROR: format %s not defined for starpar, skipping\n", NAME); \ - *Nstarpar = 0; \ - return NULL; \ - } - -# define CONVERT_FORMAT(NAME, FORMAT, TYPE) \ + *Nstarpar = 0; \ + return NULL; \ + } + +# define CONVERT_FORMAT(NAME, FORMAT, TYPE) \ if (!strcmp (extname, NAME)) { \ StarPar_##TYPE *tmpStarPar; \ @@ -777,6 +794,7 @@ exit (2); \ } \ - starpar = StarPar_##TYPE##_ToInternal (tmpStarPar, *Nstarpar); \ + starpar = StarPar_##TYPE##_ToInternal (tmpStarPar, *Nstarpar); \ free (tmpStarPar); \ + ftable[0].buffer = NULL; \ *format = DVO_FORMAT_##FORMAT; \ return (starpar); } @@ -817,11 +835,11 @@ int StarParToFtable (FTable *ftable, StarPar *starpar, off_t Nstarpar, DVOCatFormat format, int swapFromNative) { -# define FORMAT_CASE(FORMAT, TYPE) \ - case DVO_FORMAT_##FORMAT: { \ - StarPar_##TYPE *tmpStarPar; \ - tmpStarPar = StarParInternalTo_##TYPE (starpar, Nstarpar); \ - gfits_table_set_StarPar_##TYPE (ftable, tmpStarPar, Nstarpar, swapFromNative); \ - free (tmpStarPar); \ - break; } +# define FORMAT_CASE(FORMAT, TYPE) \ + case DVO_FORMAT_##FORMAT: { \ + StarPar_##TYPE *tmpStarPar; \ + tmpStarPar = StarParInternalTo_##TYPE (starpar, Nstarpar); \ + gfits_table_set_StarPar_##TYPE (ftable, tmpStarPar, Nstarpar, swapFromNative); \ + free (tmpStarPar); \ + break; } /* convert from the internal format */ @@ -866,12 +884,26 @@ } -# define SKIPPING_FORMAT(NAME, FORMAT, TYPE) \ + if (!strcmp (extname, "DVO_GALPHOT_PS1_V5") && (ftable[0].header[0].Naxis[0] == 72)) { + GalPhot_PS1_V5_R0 *tmpGalPhot; + tmpGalPhot = gfits_table_get_GalPhot_PS1_V5_R0 (ftable, Ngalphot, NULL, &nativeBytes); + if (!tmpGalPhot) { + fprintf (stderr, "ERROR: failed to read galphots\n"); + exit (2); + } + galphot = GalPhot_PS1_V5_R0_ToInternal (tmpGalPhot, *Ngalphot); + free (tmpGalPhot); + ftable[0].buffer = NULL; + *format = DVO_FORMAT_PS1_V5; + return (galphot); + } + +# define SKIPPING_FORMAT(NAME, FORMAT, TYPE) \ if (!strcmp (extname, NAME)) { \ fprintf (stderr, "ERROR: format %s not defined for galphot, skipping\n", NAME); \ - *Ngalphot = 0; \ - return NULL; \ - } - -# define CONVERT_FORMAT(NAME, FORMAT, TYPE) \ + *Ngalphot = 0; \ + return NULL; \ + } + +# define CONVERT_FORMAT(NAME, FORMAT, TYPE) \ if (!strcmp (extname, NAME)) { \ GalPhot_##TYPE *tmpGalPhot; \ @@ -881,6 +913,7 @@ exit (2); \ } \ - galphot = GalPhot_##TYPE##_ToInternal (tmpGalPhot, *Ngalphot); \ + galphot = GalPhot_##TYPE##_ToInternal (tmpGalPhot, *Ngalphot); \ free (tmpGalPhot); \ + ftable[0].buffer = NULL; \ *format = DVO_FORMAT_##FORMAT; \ return (galphot); } @@ -907,5 +940,5 @@ SKIPPING_FORMAT ("DVO_GALPHOT_PS1_V3", PS1_V3, PS1_V3); SKIPPING_FORMAT ("DVO_GALPHOT_PS1_V4", PS1_V4, PS1_V4); - CONVERT_FORMAT ("DVO_GALPHOT_PS1_V5", PS1_V5, PS1_V5); + CONVERT_FORMAT ("DVO_GALPHOT_PS1_V5", PS1_V5, PS1_V5_R1); SKIPPING_FORMAT ("DVO_GALPHOT_PS1_SIM", PS1_SIM, PS1_SIM); # undef CONVERT_FORMAT @@ -921,11 +954,11 @@ int GalPhotToFtable (FTable *ftable, GalPhot *galphot, off_t Ngalphot, DVOCatFormat format, int swapFromNative) { -# define FORMAT_CASE(FORMAT, TYPE) \ - case DVO_FORMAT_##FORMAT: { \ - GalPhot_##TYPE *tmpGalPhot; \ - tmpGalPhot = GalPhotInternalTo_##TYPE (galphot, Ngalphot); \ - gfits_table_set_GalPhot_##TYPE (ftable, tmpGalPhot, Ngalphot, swapFromNative); \ - free (tmpGalPhot); \ - break; } +# define FORMAT_CASE(FORMAT, TYPE) \ + case DVO_FORMAT_##FORMAT: { \ + GalPhot_##TYPE *tmpGalPhot; \ + tmpGalPhot = GalPhotInternalTo_##TYPE (galphot, Ngalphot); \ + gfits_table_set_GalPhot_##TYPE (ftable, tmpGalPhot, Ngalphot, swapFromNative); \ + free (tmpGalPhot); \ + break; } /* convert from the internal format */ @@ -946,5 +979,5 @@ // FORMAT_CASE (PS1_V3, PS1_V3); // FORMAT_CASE (PS1_V4, PS1_V4); - FORMAT_CASE (PS1_V5, PS1_V5); + FORMAT_CASE (PS1_V5, PS1_V5_R1); # undef FORMAT_CASE @@ -990,21 +1023,21 @@ } -# define CONVERT_FORMAT(NAME, FORMAT, TYPE) \ - if (!strcmp (extname, NAME)) { \ - Image_##TYPE *tmpimage; \ - *format = DVO_FORMAT_##FORMAT; \ +# define CONVERT_FORMAT(NAME, FORMAT, TYPE) \ + if (!strcmp (extname, NAME)) { \ + Image_##TYPE *tmpimage; \ + *format = DVO_FORMAT_##FORMAT; \ tmpimage = gfits_table_get_Image_##TYPE (ftable, &Nimage, NULL, NULL); \ - if (!tmpimage) { \ - fprintf (stderr, "ERROR: failed to read images\n"); \ - exit (2); \ - } \ + if (!tmpimage) { \ + fprintf (stderr, "ERROR: failed to read images\n"); \ + exit (2); \ + } \ off_t Nalloc = gfits_data_pad_size(Nimage*sizeof(Image)); \ ftable[0].buffer = (char *) Image_##TYPE##_ToInternal (tmpimage, Nimage, Nalloc); \ - free (tmpimage); \ - gfits_free_header (theader); \ - gfits_table_mkheader_Image (theader); \ - gfits_modify (theader, "NAXIS2", OFF_T_FMT, 1, Nimage); \ - theader[0].Naxis[1] = Nimage; \ - ftable[0].datasize = gfits_data_size (theader); \ + free (tmpimage); \ + gfits_free_header (theader); \ + gfits_table_mkheader_Image (theader); \ + gfits_modify (theader, "NAXIS2", OFF_T_FMT, 1, Nimage); \ + theader[0].Naxis[1] = Nimage; \ + ftable[0].datasize = gfits_data_size (theader); \ return (TRUE); } @@ -1035,29 +1068,29 @@ Nimage = theader[0].Naxis[1]; -# define FORMAT_CASE(FORMAT, TYPE) \ - case DVO_FORMAT_##FORMAT: { \ +# define FORMAT_CASE(FORMAT, TYPE) \ + case DVO_FORMAT_##FORMAT: { \ Image_##TYPE *tmpImage; \ - tmpImage = ImageInternalTo_##TYPE ((Image *) ftable[0].buffer, Nimage); \ - free (ftable[0].buffer); \ - gfits_free_header (ftable->header); \ - gfits_table_set_Image_##TYPE (ftable, tmpImage, Nimage, TRUE); \ - free (tmpImage); \ - break; } + tmpImage = ImageInternalTo_##TYPE ((Image *) ftable[0].buffer, Nimage); \ + free (ftable[0].buffer); \ + gfits_free_header (ftable->header); \ + gfits_table_set_Image_##TYPE (ftable, tmpImage, Nimage, TRUE); \ + free (tmpImage); \ + break; } /* convert from the internal format */ switch (format) { - FORMAT_CASE (ELIXIR, Elixir); - FORMAT_CASE (LONEOS, Loneos); - FORMAT_CASE (PANSTARRS_DEV_0, Panstarrs_DEV_0); - FORMAT_CASE (PANSTARRS_DEV_1, Panstarrs_DEV_1); - FORMAT_CASE (PS1_DEV_1, PS1_DEV_1); - FORMAT_CASE (PS1_DEV_2, PS1_DEV_2); - FORMAT_CASE (PS1_V1, PS1_V1); - FORMAT_CASE (PS1_V2, PS1_V2); - FORMAT_CASE (PS1_V3, PS1_V3); - FORMAT_CASE (PS1_V4, PS1_V4); - FORMAT_CASE (PS1_V5, PS1_V5); - FORMAT_CASE (PS1_REF, PS1_REF); - FORMAT_CASE (PS1_SIM, PS1_SIM); + FORMAT_CASE (ELIXIR, Elixir); + FORMAT_CASE (LONEOS, Loneos); + FORMAT_CASE (PANSTARRS_DEV_0, Panstarrs_DEV_0); + FORMAT_CASE (PANSTARRS_DEV_1, Panstarrs_DEV_1); + FORMAT_CASE (PS1_DEV_1, PS1_DEV_1); + FORMAT_CASE (PS1_DEV_2, PS1_DEV_2); + FORMAT_CASE (PS1_V1, PS1_V1); + FORMAT_CASE (PS1_V2, PS1_V2); + FORMAT_CASE (PS1_V3, PS1_V3); + FORMAT_CASE (PS1_V4, PS1_V4); + FORMAT_CASE (PS1_V5, PS1_V5); + FORMAT_CASE (PS1_REF, PS1_REF); + FORMAT_CASE (PS1_SIM, PS1_SIM); # undef FORMAT_CASE @@ -1076,39 +1109,39 @@ Nrow = vtable[0].Nrow; -# define FORMAT_CASE(FORMAT, TYPE) \ - case DVO_FORMAT_##FORMAT: { \ - Image_##TYPE *tmpImage; \ - /* convert table rows from internal to external format */ \ - for (i = 0; i < Nrow; i++) { \ - tmpImage = ImageInternalTo_##TYPE ((Image *) vtable[0].buffer[i], 1); \ - gfits_convert_Image_##TYPE (tmpImage, sizeof(Image_##TYPE), 1); \ - free (vtable[0].buffer[i]); \ - vtable[0].buffer[i] = (char *) tmpImage; \ - } \ - /* convert header from old format to new format */ \ - gfits_scan (theader, "NAXIS2", OFF_T_FMT, 1, &Nimage); \ - gfits_free_header (theader); \ - gfits_table_mkheader_Image_##TYPE (theader); \ - gfits_modify (theader, "NAXIS2", OFF_T_FMT, 1, Nimage); \ - theader[0].Naxis[1] = Nimage; \ - vtable[0].datasize = gfits_data_size (theader); \ - return (TRUE); } +# define FORMAT_CASE(FORMAT, TYPE) \ + case DVO_FORMAT_##FORMAT: { \ + Image_##TYPE *tmpImage; \ + /* convert table rows from internal to external format */ \ + for (i = 0; i < Nrow; i++) { \ + tmpImage = ImageInternalTo_##TYPE ((Image *) vtable[0].buffer[i], 1); \ + gfits_convert_Image_##TYPE (tmpImage, sizeof(Image_##TYPE), 1); \ + free (vtable[0].buffer[i]); \ + vtable[0].buffer[i] = (char *) tmpImage; \ + } \ + /* convert header from old format to new format */ \ + gfits_scan (theader, "NAXIS2", OFF_T_FMT, 1, &Nimage); \ + gfits_free_header (theader); \ + gfits_table_mkheader_Image_##TYPE (theader); \ + gfits_modify (theader, "NAXIS2", OFF_T_FMT, 1, Nimage); \ + theader[0].Naxis[1] = Nimage; \ + vtable[0].datasize = gfits_data_size (theader); \ + return (TRUE); } /* convert from the internal format */ switch (format) { - FORMAT_CASE (ELIXIR, Elixir); - FORMAT_CASE (LONEOS, Loneos); - FORMAT_CASE (PANSTARRS_DEV_0, Panstarrs_DEV_0); - FORMAT_CASE (PANSTARRS_DEV_1, Panstarrs_DEV_1); - FORMAT_CASE (PS1_DEV_1, PS1_DEV_1); - FORMAT_CASE (PS1_DEV_2, PS1_DEV_2); - FORMAT_CASE (PS1_V1, PS1_V1); - FORMAT_CASE (PS1_V2, PS1_V2); - FORMAT_CASE (PS1_V3, PS1_V3); - FORMAT_CASE (PS1_V4, PS1_V4); - FORMAT_CASE (PS1_V5, PS1_V5); - FORMAT_CASE (PS1_REF, PS1_REF); - FORMAT_CASE (PS1_SIM, PS1_SIM); + FORMAT_CASE (ELIXIR, Elixir); + FORMAT_CASE (LONEOS, Loneos); + FORMAT_CASE (PANSTARRS_DEV_0, Panstarrs_DEV_0); + FORMAT_CASE (PANSTARRS_DEV_1, Panstarrs_DEV_1); + FORMAT_CASE (PS1_DEV_1, PS1_DEV_1); + FORMAT_CASE (PS1_DEV_2, PS1_DEV_2); + FORMAT_CASE (PS1_V1, PS1_V1); + FORMAT_CASE (PS1_V2, PS1_V2); + FORMAT_CASE (PS1_V3, PS1_V3); + FORMAT_CASE (PS1_V4, PS1_V4); + FORMAT_CASE (PS1_V5, PS1_V5); + FORMAT_CASE (PS1_REF, PS1_REF); + FORMAT_CASE (PS1_SIM, PS1_SIM); # undef FORMAT_CASE Index: trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V5.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V5.c (revision 38970) +++ trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V5.c (revision 38986) @@ -218,5 +218,7 @@ out[i].catID = in[i].catID; out[i].extID = in[i].extID; - out[i].extIDgc = in[i].extIDgc; + + out[i].uRgal = in[i].uRgal; + out[i].uDgal = in[i].uDgal; } return (out); @@ -286,5 +288,7 @@ out[i].catID = in[i].catID; out[i].extID = in[i].extID; - out[i].extIDgc = in[i].extIDgc; + + out[i].uRgal = in[i].uRgal; + out[i].uDgal = in[i].uDgal; } return (out); @@ -979,5 +983,5 @@ } -GalPhot *GalPhot_PS1_V5_ToInternal (GalPhot_PS1_V5 *in, off_t Nvalues) { +GalPhot *GalPhot_PS1_V5_R1_ToInternal (GalPhot_PS1_V5_R1 *in, off_t Nvalues) { off_t i; @@ -1013,10 +1017,44 @@ } -GalPhot_PS1_V5 *GalPhotInternalTo_PS1_V5 (GalPhot *in, off_t Nvalues) { - - off_t i; - GalPhot_PS1_V5 *out; - - ALLOCATE_ZERO (out, GalPhot_PS1_V5, Nvalues); +GalPhot *GalPhot_PS1_V5_R0_ToInternal (GalPhot_PS1_V5_R0 *in, off_t Nvalues) { + + off_t i; + GalPhot *out; + + ALLOCATE_ZERO (out, GalPhot, Nvalues); + + for (i = 0; i < Nvalues; i++) { + dvo_galphot_init (&out[i]); + + out[i].Xfit = in[i].Xfit; + out[i].Yfit = in[i].Yfit; + out[i].mag = in[i].mag; + out[i].magErr = in[i].magErr; + out[i].majorAxis = in[i].majorAxis; + out[i].minorAxis = in[i].minorAxis; + out[i].majorAxisErr = in[i].majorAxisErr; + out[i].minorAxisErr = in[i].minorAxisErr; + out[i].theta = in[i].theta; + out[i].thetaErr = in[i].thetaErr; + out[i].index = in[i].index; + out[i].chisq = in[i].chisq; + out[i].Npix = in[i].Npix; + out[i].objID = in[i].objID; + out[i].catID = in[i].catID; + out[i].detID = in[i].detID; + out[i].imageID = in[i].imageID; +// out[i].averef = in[i].averef; + out[i].photcode = in[i].photcode; + out[i].modelType = in[i].modelType; + } + return (out); +} + +GalPhot_PS1_V5_R1 *GalPhotInternalTo_PS1_V5_R1 (GalPhot *in, off_t Nvalues) { + + off_t i; + GalPhot_PS1_V5_R1 *out; + + ALLOCATE_ZERO (out, GalPhot_PS1_V5_R1, Nvalues); for (i = 0; i < Nvalues; i++) { @@ -1318,5 +1356,4 @@ out[i].catID = in[i].catID; out[i].extID = in[i].extID; - out[i].extIDgc = in[i].extIDgc; } return (out); Index: trunk/Ohana/src/libdvo/src/dvo_photcode_ops.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_photcode_ops.c (revision 38970) +++ trunk/Ohana/src/libdvo/src/dvo_photcode_ops.c (revision 38986) @@ -205,4 +205,5 @@ /* return equivalent photcode.code for given code */ int GetPhotcodeEquivCodebyCode (int code) { + myAssert (photcodes, "photcodes not loaded"); int entry; @@ -218,4 +219,5 @@ // returns Nsec if code is PRI/SEC, else -1 int GetPhotcodeNsec (int code) { + myAssert (photcodes, "photcodes not loaded"); int Nsec; @@ -230,4 +232,5 @@ /* Nsec of 0 is PRI */ PhotCode *GetPhotcodebyNsec (int Nsec) { + myAssert (photcodes, "photcodes not loaded"); int Ncode; @@ -245,4 +248,5 @@ int GetPhotcodeNsecfilt () { + myAssert (photcodes, "photcodes not loaded"); return (photcodes[0].Nsecfilt); } @@ -745,4 +749,26 @@ } return (Mstdev); +} + +// return the number of detections in this filter (gpc1) +int PhotNcode (PhotCode *code, Average *average, SecFilt *secfilt) { + + if (code == NULL) return 0; + + int Ns = photcodes[0].hashNsec[code[0].code]; + if (Ns == -1) return 0; + + return secfilt[Ns].Ncode; +} + +// return the number of detections in this filter (gpc1) +int PhotSecfiltFlags (PhotCode *code, Average *average, SecFilt *secfilt) { + + if (code == NULL) return 0; + + int Ns = photcodes[0].hashNsec[code[0].code]; + if (Ns == -1) return 0; + + return secfilt[Ns].flags; } Index: trunk/Ohana/src/libdvo/test/tap_AstromOffsetMapIO.c =================================================================== --- trunk/Ohana/src/libdvo/test/tap_AstromOffsetMapIO.c (revision 38986) +++ trunk/Ohana/src/libdvo/test/tap_AstromOffsetMapIO.c (revision 38986) @@ -0,0 +1,95 @@ +# include "dvo.h" +# include "tap_ohana.h" + +int compare_tables (AstromOffsetTable *table_src, AstromOffsetTable *table_tgt); + +int main (void) { + + int i, j, k, status; + + plan_tests (181); + + diag ("libdvo AstromOffsetMapIO.c tests"); + + /*** create a small, simple reference AstroMap ***/ + AstromOffsetTable *table = NULL; + ALLOCATE (table, AstromOffsetTable, 1); + + table->Nmap = 4; + ALLOCATE (table->map, AstromOffsetMap *, table->Nmap); + + // assign the map values (this allocates just the area needed for each image, not the + // full 6x6, saving some memory while doing the analysis) + for (i = 0; i < table->Nmap; i++) { + ALLOCATE (table->map[i], AstromOffsetMap, 1); + + table->map[i][0].Nx = i + 1; + table->map[i][0].Ny = i + 1; + table->map[i][0].tableID = i + 10; + table->map[i][0].imageID = i + 20; + + // STORE THESE VALUES? + table->map[i][0].dX = 2*i + 0.5; + table->map[i][0].dY = 3*i + 0.5; + + // since this was on disk, we obviously keep it + table->map[i][0].keep = TRUE; + + ALLOCATE (table->map[i][0].dXv, float *, table->map[i][0].Nx); + ALLOCATE (table->map[i][0].dYv, float *, table->map[i][0].Nx); + + for (j = 0; j < table->map[i][0].Nx; j++) { + ALLOCATE (table->map[i][0].dXv[j], float, table->map[i][0].Ny); + ALLOCATE (table->map[i][0].dYv[j], float, table->map[i][0].Ny); + + for (k = 0; k < table->map[i][0].Ny; k++) { + table->map[i][0].dXv[j][k] = i+j+k+1; + table->map[i][0].dYv[j][k] = i+j+k+2; + } + } + } + + status = AstromOffsetMapSave (table, "test.v0.fits"); + ok (status, "wrote test table to file"); + + AstromOffsetTable *table_full = AstromOffsetMapLoad ("test.v0.fits", 0, TRUE); + ok (table_full, "read test table from file"); + + /*** compare table_full and table above ***/ + diag ("compare table_full and table"); + compare_tables (table_full, table); + + /*** compare AstromOffsetMapLoad in blocks vs single read ***/ + AstromOffsetTable *table_rows = AstromOffsetMapLoad ("test.v0.fits", 100, TRUE); + ok (table_rows, "read in table in 100 rows per read"); + + diag ("compare table_rows and table"); + compare_tables (table_rows, table); + + return exit_status(); +} + +int compare_tables (AstromOffsetTable *table_src, AstromOffsetTable *table_tgt) { + + int i, j, k; + + ok (table_tgt->Nmap == table_src->Nmap, "number of maps matches"); + + for (i = 0; i < table_src->Nmap; i++) { + ok (table_tgt->map[i][0].Nx == table_src->map[i][0].Nx, "Nx sizes match"); + ok (table_tgt->map[i][0].Ny == table_src->map[i][0].Ny, "Ny sizes match"); + ok (table_tgt->map[i][0].tableID == table_src->map[i][0].tableID, "tableIDs match"); + ok (table_tgt->map[i][0].imageID == table_src->map[i][0].imageID, "imageIDs match"); + ok (table_tgt->map[i][0].dX == table_src->map[i][0].dX, "dX matches"); + ok (table_tgt->map[i][0].dY == table_src->map[i][0].dY, "dY matches"); + ok (table_tgt->map[i][0].keep == table_src->map[i][0].keep, "keep matches"); + + for (j = 0; j < table_src->map[i][0].Nx; j++) { + for (k = 0; k < table_src->map[i][0].Ny; k++) { + ok (table_tgt->map[i][0].dXv[j][k] == table_src->map[i][0].dXv[j][k], "dX values match"); + ok (table_tgt->map[i][0].dYv[j][k] == table_src->map[i][0].dYv[j][k], "dY values match"); + } + } + } + return TRUE; +}