Index: /branches/eam_branches/ipp-20111122/Ohana/src/libautocode/def/FlatCorrection.d =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/libautocode/def/FlatCorrection.d (revision 33103) +++ /branches/eam_branches/ipp-20111122/Ohana/src/libautocode/def/FlatCorrection.d (revision 33104) @@ -2,5 +2,5 @@ EXTNAME FLAT_CORRECTION TYPE BINTABLE -SIZE 20 +SIZE 16 FIELD ID, ID, int Index: /branches/eam_branches/ipp-20111122/Ohana/src/libdvo/Makefile =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/libdvo/Makefile (revision 33103) +++ /branches/eam_branches/ipp-20111122/Ohana/src/libdvo/Makefile (revision 33104) @@ -70,4 +70,5 @@ $(SRC)/dvo_convert_PS1_V2.$(ARCH).o \ $(SRC)/dvo_convert_PS1_REF.$(ARCH).o \ +$(SRC)/flatcorr_io.$(ARCH).o \ $(SRC)/skyregion_io.$(ARCH).o \ $(SRC)/skyregion_gsc.$(ARCH).o \ Index: /branches/eam_branches/ipp-20111122/Ohana/src/libdvo/include/dvo.h =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/libdvo/include/dvo.h (revision 33103) +++ /branches/eam_branches/ipp-20111122/Ohana/src/libdvo/include/dvo.h (revision 33104) @@ -205,4 +205,13 @@ } SkyList; +typedef struct { + off_t Nimage; + FlatCorrectionImage *image; + off_t Ncorr; + int *IDtoSeq; + FlatCorrection *corr; + float ***offset; // the correction images represented as a set of arrays (same sequence as *image) +} FlatCorrectionTable; + // special-case function: CMF_PS1_V2 *gfits_table_get_CMF_PS1_V1_Alt (FTable *ftable, off_t *Ndata, char *swapped); @@ -484,4 +493,10 @@ void dvo_image_create (FITS_DB *db, double ZeroPoint); +/* flatcorr APIs */ +FlatCorrectionTable *FlatCorrectionLoad (char *filename, int VERBOSE); +int FlatCorrectionInternal(FlatCorrectionTable *flatcorrTable); +int FlatCorrectionSave (FlatCorrectionTable *flatcorrTable, char *filename); +float FlatCorrectionOffset (FlatCorrectionTable *flatcorr, int ID, int X, int Y); + /* skyregion APIs */ int SkyTableSave PROTO((SkyTable *table, char *filename)); Index: /branches/eam_branches/ipp-20111122/Ohana/src/libdvo/src/flatcorr_io.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/libdvo/src/flatcorr_io.c (revision 33103) +++ /branches/eam_branches/ipp-20111122/Ohana/src/libdvo/src/flatcorr_io.c (revision 33104) @@ -10,12 +10,4 @@ */ -typedef struct { - int Nimage; - FlatCorrectionImage *image; - int Ncorr; - FlatCorrection *corr; - float ***offset; // the correction images represented as a set of arrays (same sequence as *image) -} FlatCorrectionTable; - FlatCorrectionTable *FlatCorrectionLoad (char *filename, int VERBOSE) { @@ -25,5 +17,4 @@ Header theader; FTable ftable; - int i; FlatCorrectionTable *flatcorrTable; @@ -85,38 +76,58 @@ // convert the table format to map format: + FlatCorrectionInternal (flatcorrTable); - // generate the arrays to hold the corrections + return (flatcorrTable); +} + +int FlatCorrectionInternal(FlatCorrectionTable *flatcorrTable) { + + int i, j; + + // assert that the internal arrays are not yet allocated? + + // generate the arrays to hold the corrections & initialize ALLOCATE (flatcorrTable->offset, float **, flatcorrTable->Nimage); for (i = 0; i < flatcorrTable->Nimage; i++) { - ALLOCATE (flatcorrTable->offset[i], float *, flatcorrTable->image[i]->Nx); - for (j = 0; j < flatcorrTable->image[i]->Nx; j++) { - ALLOCATE (flatcorrTable->offset[i][j], float, flatcorrTable->image[i]->Ny); - memset (flatcorrTable->offset[i][j], 0.0, flatcorrTable->image[i]->Ny*sizeof(float)); - } + ALLOCATE (flatcorrTable->offset[i], float *, flatcorrTable->image[i].Nx); + for (j = 0; j < flatcorrTable->image[i].Nx; j++) { + ALLOCATE (flatcorrTable->offset[i][j], float, flatcorrTable->image[i].Ny); + memset (flatcorrTable->offset[i][j], 0.0, flatcorrTable->image[i].Ny*sizeof(float)); + } } // create the lookup table (ID -> Seq) + + // find the max value of ID + int MaxID = 0; + for (i = 0; i < flatcorrTable->Nimage; i++) { + MaxID = MAX(flatcorrTable->image[i].ID, MaxID); + } + + // generate the index and init values to -1 ALLOCATE (flatcorrTable->IDtoSeq, int, MaxID); + for (i = 0; i < MaxID; i++) { + flatcorrTable->IDtoSeq[i] = -1; + } + + // assign the ID values for each image for (i = 0; i < flatcorrTable->Nimage; i++) { - flatcorrTable->IDtoSeq[i] = -1; - } - for (i = 0; i < flatcorrTable->Nimage; i++) { - ID = flatcorrTable->image[i].ID; - assert (flatcorrTable->IDtoSeq[ID] == -1); - flatcorrTable->IDtoSeq[ID] = i; + int ID = flatcorrTable->image[i].ID; + assert (flatcorrTable->IDtoSeq[ID] == -1); + flatcorrTable->IDtoSeq[ID] = i; } // assign the known specific correction values for (i = 0; i < flatcorrTable->Ncorr; i++) { - ID = flatcorrTable->corr[i].ID; - x = flatcorrTable->corr[i].x; - y = flatcorrTable->corr[i].y; - seq = flatcorrTable->IDtoSeq[ID]; - assert (seq != -1); + int ID = flatcorrTable->corr[i].ID; + int x = flatcorrTable->corr[i].x; + int y = flatcorrTable->corr[i].y; + int seq = flatcorrTable->IDtoSeq[ID]; + assert (seq != -1); - flatcorrTable->offset[i][x][y] = flatcorrTable->corr[i].offset; + flatcorrTable->offset[i][x][y] = flatcorrTable->corr[i].offset; } - return (flatcorrTable); + return TRUE; } @@ -142,5 +153,5 @@ ftableCorr.header = &theaderCorr; - gfits_table_set_FlatCorrectionCorr (&ftableCorr, flatcorrTable->corr, flatcorrTable->Ncorr); + gfits_table_set_FlatCorrection (&ftableCorr, flatcorrTable->corr, flatcorrTable->Ncorr); f = fopen (filename, "w"); @@ -163,15 +174,15 @@ float FlatCorrectionOffset (FlatCorrectionTable *flatcorr, int ID, int X, int Y) { - // validate the flat_id (not out of range?) - seq = flatcorr->IDtoSeq[ID]; + // validate the flat_id (not out of range?) + int seq = flatcorr->IDtoSeq[ID]; - // convert X,Y to Xbin, Ybin: - Xbin = MAX(MIN(X / flatcorr->image[seq].DX, flatcorr->image[seq].NX - 1), 0); - Ybin = MAX(MIN(Y / flatcorr->image[seq].DY, flatcorr->image[seq].NY - 1), 0); + // convert X,Y to Xbin, Ybin: + int Xbin = MAX(MIN(X * flatcorr->image[seq].Nx / flatcorr->image[seq].DX, flatcorr->image[seq].Nx - 1), 0); + int Ybin = MAX(MIN(Y * flatcorr->image[seq].Ny / flatcorr->image[seq].DY, flatcorr->image[seq].Ny - 1), 0); - // XXX warn if X,Y are out of range? not super important unless *way* out of range + // XXX warn if X,Y are out of range? not super important unless *way* out of range - offset = flatcorr->offset[seq][Xbin][Ybin]; - return offset; + float offset = flatcorr->offset[seq][Xbin][Ybin]; + return offset; } Index: /branches/eam_branches/ipp-20111122/Ohana/src/uniphot/Makefile =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/uniphot/Makefile (revision 33103) +++ /branches/eam_branches/ipp-20111122/Ohana/src/uniphot/Makefile (revision 33104) @@ -52,4 +52,5 @@ $(SRC)/load_images.$(ARCH).o \ $(SRC)/match_zpts_to_images.$(ARCH).o \ +$(SRC)/match_flatcorr_to_images.$(ARCH).o \ $(SRC)/update_catalog_setphot.$(ARCH).o \ $(SRC)/SetSignals.$(ARCH).o \ Index: /branches/eam_branches/ipp-20111122/Ohana/src/uniphot/include/uniphot.h =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/uniphot/include/uniphot.h (revision 33103) +++ /branches/eam_branches/ipp-20111122/Ohana/src/uniphot/include/uniphot.h (revision 33104) @@ -66,4 +66,5 @@ char STATMODE[64]; int VERBOSE; +int UBERCAL; // load the supplied ubercal zero point fits table (with flat-field corrections) int NLOOP; int TimeSelect; @@ -132,6 +133,9 @@ Image *load_images_setphot PROTO((FITS_DB *db, off_t *Nimage)); int match_zpts_to_images PROTO((Image *image, off_t Nimage, ZptTable *zpts, int Nzpts)); -int update_setphot PROTO((Image *image, off_t Nimage)); -void update_catalog_setphot PROTO((Catalog *catalog, Image *image, off_t *index, off_t Nimage)); +int update_setphot PROTO((Image *image, off_t Nimage, FlatCorrectionTable *flatcorr)); +void update_catalog_setphot PROTO((Catalog *catalog, Image *image, off_t *index, off_t Nimage, FlatCorrectionTable *flatcorr)); + +ZptTable *load_zpt_ubercal (char *filename, int *nzpts, FlatCorrectionTable *flatcorrTable); +int match_flatcorr_to_images (Image *image, off_t Nimage, FlatCorrectionTable *flatcorrTable); /*** time/coord conversion functions not supplied by libohana ***/ Index: /branches/eam_branches/ipp-20111122/Ohana/src/uniphot/src/load_zpt_table.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/uniphot/src/load_zpt_table.c (revision 33103) +++ /branches/eam_branches/ipp-20111122/Ohana/src/uniphot/src/load_zpt_table.c (revision 33104) @@ -42,9 +42,27 @@ } -ZptTable *load_zpt_ubercal (char *filename, int *nzpts) { +# define NFILTER 5 +# define NSEASON 4 +char filter[NFILTER] = {'g', 'r', 'i', 'z', 'y'}; +double tstart[NSEASON] = {55000.0, 55296.0, 55327.0, 55662.0}; +double tstop[NSEASON] = {55296.0, 55327.0, 55662.0, 60000.0}; + +# define NCHIP_X 8 +# define NCHIP_Y 8 +# define NCELL_X 2 +# define NCELL_Y 2 + +ZptTable *load_zpt_ubercal (char *filename, int *nzpts, FlatCorrectionTable *flatcorrTable) { + + int i, nfilter, nseason, ix, iy, ixc, iyc, Ncol, seq; + off_t Nrow; + char type[16]; int Nzpts, NZPTS; ZptTable *zpts; - double zpt, mjd, zpt_err; + Header header; + Header theader; + Matrix matrix; + FTable ftable; FILE *f; @@ -57,5 +75,4 @@ /* load in PHU segment (ignore) */ - ftable.header = &theader; if (!gfits_fread_header (f, &header)) { if (VERBOSE) fprintf (stderr, "can't read Flat Correction header\n"); @@ -74,49 +91,92 @@ ALLOCATE (zpts, ZptTable, NZPTS); - ALLOCATE (flatcorr, FlatCorrection, 4*5*X); + // we have 5 filters, and 4 flat-field correction sets for each + flatcorrTable->Ncorr = NFILTER*NSEASON*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y; + flatcorrTable->Nimage = NFILTER*NSEASON*NCHIP_X*NCHIP_Y; - for (i = 0; i < 5; i++) { - // load data for this header - if (!gfits_load_header (f, &theader)) return (FALSE); + ALLOCATE (flatcorrTable->corr, FlatCorrection, flatcorrTable->Ncorr); + ALLOCATE (flatcorrTable->image, FlatCorrectionImage, flatcorrTable->Nimage); - // read the fits table bytes - if (gfits_fread_ftable_data (f, &table, FALSE)) return (TRUE); + int corrID = 1; + int Nimage = 0; + ftable.header = &theader; + for (nfilter = 0; nfilter < NFILTER; nfilter++) { + // load data for this header + if (!gfits_load_header (f, &theader)) return (NULL); + + // read the fits table bytes + if (gfits_fread_ftable_data (f, &ftable, FALSE)) return (NULL); - // skip over remaining bytes in data segment - fseeko (f, table[0].datasize - table[0].validsize, SEEK_CUR); + // skip over remaining bytes in data segment + fseeko (f, ftable.datasize - ftable.validsize, SEEK_CUR); - // load data for this header - if (!gfits_load_header (f, &header)) return (FALSE); + // need to create and assign to flat-field correction + double *mjd = gfits_get_bintable_column_data (&theader, &ftable, "mjd_obs", type, &Nrow, &Ncol); + assert (!strcmp(type, "double")); - // read the fits table bytes - if (gfits_fread_matrix (f, &matrix, FALSE)) return (TRUE); + double *zp = gfits_get_bintable_column_data (&theader, &ftable, "zp", type, &Nrow, &Ncol); + assert (!strcmp(type, "double")); + + double *zperr = gfits_get_bintable_column_data (&theader, &ftable, "resid", type, &Nrow, &Ncol); + assert (!strcmp(type, "double")); + + NZPTS += Nrow; + REALLOCATE (zpts, ZptTable, NZPTS); + for (i = 0; i < Nrow; i++) { + zpts[i+Nzpts].time = ohana_mjd_to_sec(mjd[i]); + zpts[i+Nzpts].zpt = zp[i]; + zpts[i+Nzpts].zpt_err = zperr[i]; + } + Nzpts += Nrow; - // need to create and assign to flat-field correction - double *mjd = gfits_get_bintable_column_data (header, table, "mjd_obs", &type, &Nrow, &Ncol); - assert (!strcmp(type, "double")); + // the image contains the flat-field corrections - double *zpts = gfits_get_bintable_column_data (header, table, "zp", &type, &Nrow, &Ncol); - assert (!strcmp(type, "double")); - - double *zperr = gfits_get_bintable_column_data (header, table, "resid", &type, &Nrow, &Ncol); - assert (!strcmp(type, "double")); - - NZPTS += Nrow; - REALLOCATE (zpts, ZptTable, NZPTS); - for (i = 0; i < Nrow; i++) { - zpts[i+Nzpts].time = ohana_mjd_to_sec(mjd[i]); - zpts[i+Nzpts].zpt = zp[i]; - zpts[i+Nzpts].zpt_err = zperr[i]; + // load data for this header + if (!gfits_load_header (f, &header)) return (NULL); + + // read the fits table bytes + if (gfits_fread_matrix (f, &matrix, FALSE)) return (NULL); + + for (nseason = 0; nseason < NSEASON; nseason++) { // seasons + for (iy = 0; iy < NCHIP_Y; iy++) { // y-chip + for (ix = 0; ix < NCHIP_X; ix++) { // x-chip + flatcorrTable->image[Nimage].Nx = NCELL_X; + flatcorrTable->image[Nimage].Ny = NCELL_Y; + flatcorrTable->image[Nimage].ID = corrID; + flatcorrTable->image[Nimage].DX = 4880; + flatcorrTable->image[Nimage].DY = 4864; + flatcorrTable->image[Nimage].tstart = ohana_mjd_to_sec(tstart[nseason]); + flatcorrTable->image[Nimage].tstop = ohana_mjd_to_sec(tstop[nseason]); + + // photcode name + char photname[64]; + snprintf (photname, 64, "GPC1.%c.XY%d%d", filter[nfilter], ix, iy); + // note that the XY00, XY07, etc, chips will have photcode values of 0 + flatcorrTable->image[Nimage].photcode = GetPhotcodeCodebyName(photname); + flatcorrTable->image[Nimage].tstart = tstart[nseason]; + flatcorrTable->image[Nimage].tstop = tstop[nseason]; + + for (iyc = 0; iyc < NCELL_Y; iyc++) { + for (ixc = 0; ixc < NCELL_X; ixc++) { + if (ix > 3) { + seq = nfilter*NSEASON*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + nseason*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + (iy*NCELL_Y + NCELL_Y - 1 - iyc)*NCHIP_X + (ix*NCELL_X + NCELL_X - 1 - ixc); + } else { + seq = nfilter*NSEASON*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + nseason*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + (iy*NCELL_Y + iyc)*NCHIP_X + (ix*NCELL_X + ixc); + } + flatcorrTable->corr[seq].x = ixc; + flatcorrTable->corr[seq].y = iyc; + flatcorrTable->corr[seq].offset = matrix.buffer[seq]; + flatcorrTable->corr[seq].ID = corrID; + } + } + corrID ++; + Nimage ++; + } } - Nzpts += Nrow; - - // the image is the flat-field correction. - // need to byteswap or whatever + } } - if (status != EOF) { - fprintf (stderr, "unexpected formatting on line %d\n", Nzpts); - exit (2); - } + /*** convert from corr,image format to offsets ***/ + FlatCorrectionInternal (flatcorrTable); fprintf (stderr, "loaded %d zero points\n", Nzpts); Index: /branches/eam_branches/ipp-20111122/Ohana/src/uniphot/src/match_flatcorr_to_images.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/uniphot/src/match_flatcorr_to_images.c (revision 33104) +++ /branches/eam_branches/ipp-20111122/Ohana/src/uniphot/src/match_flatcorr_to_images.c (revision 33104) @@ -0,0 +1,23 @@ +# include "uniphot.h" + +int match_flatcorr_to_images (Image *image, off_t Nimage, FlatCorrectionTable *flatcorrTable) { + + int i, j; + + for (i = 0; i < Nimage; i++) { + + int found = FALSE; + for (j = 0; !found && (j < flatcorrTable->Nimage); j++) { + + if (image[i].tzero < flatcorrTable->image[j].tstart) continue; + if (image[i].tzero > flatcorrTable->image[j].tstop) continue; + if (image[i].photcode != flatcorrTable->image[j].photcode) continue; + image[i].photom_map_id = flatcorrTable->image[j].ID; + found = TRUE; + } + if (!found) { + fprintf (stderr, "image does not match flatcorr Table ranges: %s\n", image[i].name); + } + } + return TRUE; +} Index: /branches/eam_branches/ipp-20111122/Ohana/src/uniphot/src/setphot.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/uniphot/src/setphot.c (revision 33103) +++ /branches/eam_branches/ipp-20111122/Ohana/src/uniphot/src/setphot.c (revision 33104) @@ -8,4 +8,5 @@ ZptTable *zpts; Image *image; + FlatCorrectionTable flatcorrTable; /* get configuration info, args, lockfile */ @@ -18,5 +19,8 @@ if (UBERCAL) { - zpts = load_zpt_ubercal (argv[1], &Nzpts); + zpts = load_zpt_ubercal (argv[1], &Nzpts, &flatcorrTable); + char flatcorrfile[64]; + snprintf (flatcorrfile, 64, "%s/flatcorr.fits", CATDIR); + FlatCorrectionSave(&flatcorrTable, flatcorrfile); } else { zpts = load_zpt_table (argv[1], &Nzpts); @@ -29,5 +33,9 @@ match_zpts_to_images (image, Nimage, zpts, Nzpts); - update_setphot(image, Nimage); + if (UBERCAL) { + match_flatcorr_to_images (image, Nimage, &flatcorrTable); + } + + update_setphot(image, Nimage, &flatcorrTable); // write image table Index: /branches/eam_branches/ipp-20111122/Ohana/src/uniphot/src/update_catalog_setphot.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/uniphot/src/update_catalog_setphot.c (revision 33103) +++ /branches/eam_branches/ipp-20111122/Ohana/src/uniphot/src/update_catalog_setphot.c (revision 33104) @@ -2,5 +2,5 @@ # include -int update_setphot (Image *image, off_t Nimage) { +int update_setphot (Image *image, off_t Nimage, FlatCorrectionTable *flatcorr) { SkyRegion UserPatch; @@ -54,5 +54,5 @@ } - update_catalog_setphot (&catalog, image, index, Nimage); + update_catalog_setphot (&catalog, image, index, Nimage, flatcorr); if (!UPDATE) { @@ -73,22 +73,21 @@ void update_catalog_setphot (Catalog *catalog, Image *image, off_t *index, off_t Nimage, FlatCorrectionTable *flatcorr) { - off_t i, j, m, id, idx, found; - float Mcal, dMcal; + off_t i, j, found; found = 0; for (i = 0; i < catalog[0].Naverage; i++) { - m = catalog[0].average[i].measureOffset; + off_t m = catalog[0].average[i].measureOffset; for (j = 0; j < catalog[0].average[i].Nmeasure; j++, m++) { - idx = catalog[0].measure[m].imageID; + off_t idx = catalog[0].measure[m].imageID; if (idx <= 0) continue; // detections with imageID == 0 do not have a valid image (eg, ref photcode) - id = index[idx]; - Mcal = image[id].Mcal; - dMcal = image[id].dMcal; - Mcal_offset = 0.0; + off_t id = index[idx]; + float Mcal = image[id].Mcal; + float dMcal = image[id].dMcal; + float Mcal_offset = 0.0; // if we know about a flat-field correction, then we need to apply the sub-chip correction - flat_id = image[id].photom_map_id; + int flat_id = image[id].photom_map_id; if (flat_id > 0) { Mcal_offset = FlatCorrectionOffset (flatcorr, flat_id, catalog[0].measure[m].Xccd, catalog[0].measure[m].Yccd);