Index: /branches/eam_branches/ipp-20131211/Ohana/src/addstar/src/findskycell.c =================================================================== --- /branches/eam_branches/ipp-20131211/Ohana/src/addstar/src/findskycell.c (revision 36471) +++ /branches/eam_branches/ipp-20131211/Ohana/src/addstar/src/findskycell.c (revision 36472) @@ -24,4 +24,5 @@ int mktree (char *treefile, char *catdir); +int mklocal (char *treefile, char *catdir); int apply_tree (char *treefile, char *datafile); @@ -478,5 +479,5 @@ int mklocal (char *treefile, char *catdir) { - int i, j, zone, band, status; + int i, status; FITS_DB db; Image *image; @@ -509,9 +510,10 @@ // generate an empty BoundaryTree - TessallationBoundary *tess = NULL; - ALLOCATE (tess, TessallationBoundary, Nimage); + int Ntess = Nimage; + TessellationTable *tess = NULL; + ALLOCATE (tess, TessellationTable, Ntess); // find the RA,DEC of the image centers & assign to cells - for (i = 0; i < Nimage; i++) { + for (i = 0; i < Ntess; i++) { // user supplied values, do not try to derive from Image.dat tess[i].NX_SUB = NX_SUB; @@ -523,6 +525,6 @@ XY_to_RD (&ra, &dec, x, y, &image[i].coords); - tess[i].ra = ra; - tess[i].dec = dec; + tess[i].Ro = ra; + tess[i].Do = dec; tess[i].Xo = x; tess[i].Yo = y; @@ -563,6 +565,6 @@ // XXX I don't really want to do the work of discovering the rule... - memcpy (tess[i].basename, BASENAME); tess[i].Nbasename = strlen(BASENAME); + tess[i].basename = strcreate(BASENAME); tess[i].projectIDoff = projectIDoff; @@ -570,5 +572,5 @@ } - BoundaryTreeSave (treefile, &tree); + TessellationTableSave (treefile, tess, Ntess); return TRUE; Index: /branches/eam_branches/ipp-20131211/Ohana/src/libdvo/Makefile =================================================================== --- /branches/eam_branches/ipp-20131211/Ohana/src/libdvo/Makefile (revision 36471) +++ /branches/eam_branches/ipp-20131211/Ohana/src/libdvo/Makefile (revision 36472) @@ -107,6 +107,6 @@ $(SRC)/convert.$(ARCH).o \ $(SRC)/HostTable.$(ARCH).o \ -$(SRC)/BoundaryTree.$(ARCH).o - +$(SRC)/BoundaryTree.$(ARCH).o \ +$(SRC)/TessellationTable.$(ARCH).o # $(SRC)/dvo_convert_panstarrs.$(ARCH).o Index: /branches/eam_branches/ipp-20131211/Ohana/src/libdvo/include/dvo.h =================================================================== --- /branches/eam_branches/ipp-20131211/Ohana/src/libdvo/include/dvo.h (revision 36471) +++ /branches/eam_branches/ipp-20131211/Ohana/src/libdvo/include/dvo.h (revision 36472) @@ -346,5 +346,5 @@ typedef enum { TESS_NONE, TESS_LOCAL, TESS_RINGS } TessType; -// TessallationTable is a structure to describe the parameters of a set of "tessellations" +// TessellationTable is a structure to describe the parameters of a set of "tessellations" // (these are not strictly tessellations but projection sets as only the non-local // versions can cover the full sky). For LOCAL projection cells, the structure describes @@ -368,7 +368,12 @@ int NY_SUB; + char *basename; + int Nbasename; + int projectIDoff; + int skycellIDoff; + TessType type; // BoundaryTree *tree; -} TessallationTable; +} TessellationTable; // a reduced-subset structure for relphot @@ -794,8 +799,15 @@ int free_tiny_values (Catalog *catalog); +BoundaryTree *BoundaryTreeLoad(char *filename); +BoundaryTree *BoundaryTreeRead(Header *headerPHU, Header *headerZone, FILE *f); + +int BoundaryTreeSave(char *filename, BoundaryTree *tree); +int BoundaryTreeWrite(FILE *f, BoundaryTree *tree); + int BoundaryTreeCellCoords (BoundaryTree *tree, int *zone, int *band, double ra, double dec); -int BoundaryTreeSave(char *filename, BoundaryTree *tree); -BoundaryTree *BoundaryTreeLoad(char *filename); int BoundaryTreeProjection (double *x, double *y, double r, double d, BoundaryTree *tree, int zone, int band); + +TessellationTable *TessellationTableLoad(char *filename, int *Ntess); +int TessellationTableSave(char *filename, TessellationTable *tess, int Ntess); void dvo_average_init (Average *average); Index: /branches/eam_branches/ipp-20131211/Ohana/src/libdvo/src/BoundaryTree.c =================================================================== --- /branches/eam_branches/ipp-20131211/Ohana/src/libdvo/src/BoundaryTree.c (revision 36471) +++ /branches/eam_branches/ipp-20131211/Ohana/src/libdvo/src/BoundaryTree.c (revision 36472) @@ -1,10 +1,10 @@ # include "dvo.h" -# define GET_COLUMN_NEW(OUT,NAME,TYPE) \ - TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \ +# define GET_COLUMN_NEW(HEADER,FTABLE,OUT,NAME,TYPE) \ + TYPE *OUT = gfits_get_bintable_column_data (HEADER, FTABLE, NAME, type, &Nrow, &Ncol); \ myAssert (!strcmp(type, #TYPE), "wrong column type"); -# define GET_COLUMN_RAW(OUT,NAME,TYPE) \ - OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \ +# define GET_COLUMN_RAW(HEADER,FTABLE,OUT,NAME,TYPE) \ + OUT = gfits_get_bintable_column_data (HEADER, FTABLE, NAME, type, &Nrow, &Ncol); \ myAssert (!strcmp(type, #TYPE), "wrong column type"); @@ -12,159 +12,182 @@ BoundaryTree *BoundaryTreeLoad(char *filename) { + + Header header; + Header theader; + Matrix matrix; + + header.buffer = NULL; + matrix.buffer = NULL; + theader.buffer = NULL; + BoundaryTree *tree = NULL; + + FILE *f = fopen (filename, "r"); + if (!f) { + fprintf (stderr, "ERROR: cannot open boundary tree file %s\n", filename); + return NULL; + } + + /* load in PHU segment (ignore) */ + if (!gfits_fread_header (f, &header)) { + if (DEBUG) fprintf (stderr, "can't read boundary tree header\n"); + goto escape; + } + if (!gfits_fread_matrix (f, &matrix, &header)) { + if (DEBUG) fprintf (stderr, "can't read boundary tree matrix\n"); + goto escape; + } + + // load data for this header + if (!gfits_load_header (f, &theader)) { + if (DEBUG) fprintf (stderr, "can't read boundary tree zone table header\n"); + goto escape; + } + + tree = BoundaryTreeRead (&header, &theader, f); + +escape: + + gfits_free_header (&theader); + gfits_free_header (&header); + gfits_free_matrix (&matrix); + fclose (f); + + return tree; +} + +// assume we are pointing at the relevant table portion +BoundaryTree *BoundaryTreeRead(Header *headerPHU, Header *headerZone, FILE *f) { int i, j, nz, nb, Ncol; off_t Nrow; char type[16]; - Header header; - Header theader; - Matrix matrix; - FTable ftable; - - header.buffer = NULL; - matrix.buffer = NULL; - ftable.buffer = NULL; - theader.buffer = NULL; + + FTable ftableZone; + + Header headerCell; + FTable ftableCell; + + ftableZone.buffer = NULL; + ftableCell.buffer = NULL; + headerCell.buffer = NULL; + + // we must have already read in the Zone table header section + ftableZone.header = headerZone; + ftableCell.header = &headerCell; + BoundaryTree *tree = NULL; - FILE *f = fopen (filename, "r"); - if (!f) { - fprintf (stderr, "ERROR: cannot open image subset file %s\n", filename); - return NULL; - } - - /* load in PHU segment (ignore) */ - if (!gfits_fread_header (f, &header)) { - if (DEBUG) fprintf (stderr, "can't read image subset header\n"); - goto escape; - } - if (!gfits_fread_matrix (f, &matrix, &header)) { - if (DEBUG) fprintf (stderr, "can't read image subset matrix\n"); - goto escape; - } - ALLOCATE (tree, BoundaryTree, 1); - gfits_scan (&header, "DEC_ORI", "%lf", 1, &tree->DEC_origin); - gfits_scan (&header, "DEC_OFF", "%lf", 1, &tree->DEC_offset); - - gfits_scan (&header, "NX_SUB", "%d", 1, &tree->NX_SUB); - gfits_scan (&header, "NY_SUB", "%d", 1, &tree->NY_SUB); - gfits_scan (&header, "PIXSCALE", "%lf", 1, &tree->dPix); - - ftable.header = &theader; + // we need to read the boundary tree parameters from the correct header + // put them in the PHU header in any case? + gfits_scan (headerPHU, "DEC_ORI", "%lf", 1, &tree->DEC_origin); + gfits_scan (headerPHU, "DEC_OFF", "%lf", 1, &tree->DEC_offset); + gfits_scan (headerPHU, "NX_SUB", "%d", 1, &tree->NX_SUB); + gfits_scan (headerPHU, "NY_SUB", "%d", 1, &tree->NY_SUB); + gfits_scan (headerPHU, "PIXSCALE", "%lf", 1, &tree->dPix); /*** zone information table ***/ - { - // load data for this header - if (!gfits_load_header (f, &theader)) goto escape; - - // read the fits table bytes - if (!gfits_fread_ftable_data (f, &ftable, FALSE)) goto escape; + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftableZone, FALSE)) goto escape; - // need to create and assign to flat-field correction - GET_COLUMN_RAW(tree->Nband, "NBAND", int); - GET_COLUMN_RAW(tree->RA_origin, "RA_ORIGIN", double); - GET_COLUMN_RAW(tree->RA_offset, "RA_OFFSET", double); - GET_COLUMN_RAW(tree->DEC_min , "DEC_MIN", double); - GET_COLUMN_RAW(tree->DEC_max , "DEC_MAX", double); - GET_COLUMN_RAW(tree->DEC_min_raw, "DEC_MIN_RAW", double); - GET_COLUMN_RAW(tree->DEC_max_raw, "DEC_MAX_RAW", double); - gfits_free_header (&theader); - gfits_free_table (&ftable); - - fprintf (stderr, "loaded data for %lld zones\n", (long long) Nrow); - tree->Nzone = Nrow; - - // allocate the storage arrays - ALLOCATE (tree->ra, double *, tree->Nzone); - ALLOCATE (tree->dec, double *, tree->Nzone); - ALLOCATE (tree->Xo, double *, tree->Nzone); - ALLOCATE (tree->Yo, double *, tree->Nzone); - ALLOCATE (tree->dX, int *, tree->Nzone); - ALLOCATE (tree->dY, int *, tree->Nzone); - ALLOCATE (tree->cell, int *, tree->Nzone); - ALLOCATE (tree->projID, int *, tree->Nzone); - ALLOCATE (tree->name, char **, tree->Nzone); - for (i = 0; i < tree->Nzone; i++) { - ALLOCATE (tree->ra[i], double, tree->Nband[i]); - ALLOCATE (tree->dec[i], double, tree->Nband[i]); - ALLOCATE (tree->Xo[i], double, tree->Nband[i]); - ALLOCATE (tree->Yo[i], double, tree->Nband[i]); - ALLOCATE (tree->dX[i], int, tree->Nband[i]); - ALLOCATE (tree->dY[i], int, tree->Nband[i]); - ALLOCATE (tree->cell[i], int, tree->Nband[i]); - ALLOCATE (tree->projID[i], int, tree->Nband[i]); - ALLOCATE (tree->name[i], char *, tree->Nband[i]); - for (j = 0; j < tree->Nband[i]; j++) { - ALLOCATE (tree->name[i][j], char, BOUNDARY_TREE_NAME_LENGTH); - } + // need to create and assign to flat-field correction + GET_COLUMN_RAW(headerZone, &ftableZone, tree->Nband, "NBAND", int); + GET_COLUMN_RAW(headerZone, &ftableZone, tree->RA_origin, "RA_ORIGIN", double); + GET_COLUMN_RAW(headerZone, &ftableZone, tree->RA_offset, "RA_OFFSET", double); + GET_COLUMN_RAW(headerZone, &ftableZone, tree->DEC_min , "DEC_MIN", double); + GET_COLUMN_RAW(headerZone, &ftableZone, tree->DEC_max , "DEC_MAX", double); + GET_COLUMN_RAW(headerZone, &ftableZone, tree->DEC_min_raw, "DEC_MIN_RAW", double); + GET_COLUMN_RAW(headerZone, &ftableZone, tree->DEC_max_raw, "DEC_MAX_RAW", double); + gfits_free_table (&ftableZone); + + fprintf (stderr, "loaded data for %lld zones\n", (long long) Nrow); + tree->Nzone = Nrow; + + // allocate the storage arrays + ALLOCATE (tree->ra, double *, tree->Nzone); + ALLOCATE (tree->dec, double *, tree->Nzone); + ALLOCATE (tree->Xo, double *, tree->Nzone); + ALLOCATE (tree->Yo, double *, tree->Nzone); + ALLOCATE (tree->dX, int *, tree->Nzone); + ALLOCATE (tree->dY, int *, tree->Nzone); + ALLOCATE (tree->cell, int *, tree->Nzone); + ALLOCATE (tree->projID, int *, tree->Nzone); + ALLOCATE (tree->name, char **, tree->Nzone); + for (i = 0; i < tree->Nzone; i++) { + ALLOCATE (tree->ra[i], double, tree->Nband[i]); + ALLOCATE (tree->dec[i], double, tree->Nband[i]); + ALLOCATE (tree->Xo[i], double, tree->Nband[i]); + ALLOCATE (tree->Yo[i], double, tree->Nband[i]); + ALLOCATE (tree->dX[i], int, tree->Nband[i]); + ALLOCATE (tree->dY[i], int, tree->Nband[i]); + ALLOCATE (tree->cell[i], int, tree->Nband[i]); + ALLOCATE (tree->projID[i], int, tree->Nband[i]); + ALLOCATE (tree->name[i], char *, tree->Nband[i]); + for (j = 0; j < tree->Nband[i]; j++) { + ALLOCATE (tree->name[i][j], char, BOUNDARY_TREE_NAME_LENGTH); } } /*** cell information table ***/ - { - // load data for this header - if (!gfits_load_header (f, &theader)) goto escape; - - // read the fits table bytes - if (!gfits_fread_ftable_data (f, &ftable, FALSE)) goto escape; + + // load data for this header + if (!gfits_load_header (f, &headerCell)) goto escape; + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftableCell, FALSE)) goto escape; - // need to create and assign to flat-field correction - GET_COLUMN_NEW(R, "RA", double); - GET_COLUMN_NEW(D, "DEC", double); - GET_COLUMN_NEW(zone, "ZONE", int); - GET_COLUMN_NEW(band, "BAND", int); - GET_COLUMN_NEW(index, "INDEX", int); - GET_COLUMN_NEW(Xo, "X_CENT", double); - GET_COLUMN_NEW(Yo, "Y_CENT", double); - GET_COLUMN_NEW(dX, "X_GRID", int); - GET_COLUMN_NEW(dY, "Y_GRID", int); - GET_COLUMN_NEW(name, "NAME", char); // XXX how is this done? - gfits_free_header (&theader); - gfits_free_table (&ftable); - fprintf (stderr, "loaded data for %lld cells\n", (long long) Nrow); - - // assign the storage arrays - for (i = 0; i < Nrow; i++) { - nz = zone[i]; - nb = band[i]; - tree->ra[nz][nb] = R[i]; - tree->dec[nz][nb] = D[i]; - tree->Xo[nz][nb] = Xo[i]; - tree->Yo[nz][nb] = Yo[i]; - tree->dX[nz][nb] = dX[i]; - tree->dY[nz][nb] = dY[i]; - tree->cell[nz][nb] = i; // XXX ? - memcpy(tree->name[nz][nb], &name[i*BOUNDARY_TREE_NAME_LENGTH], BOUNDARY_TREE_NAME_LENGTH); - // XXX parse out the ID from the name (skycell.NNNN) - tree->projID[nz][nb] = atoi(&tree->name[nz][nb][8]); - } - - free (R ); - free (D ); - free (zone ); - free (band ); - free (Xo ); - free (Yo ); - free (dX ); - free (dY ); - free (index ); - free (name ); - } - - gfits_free_header (&header); - gfits_free_matrix (&matrix); - fclose (f); + // need to create and assign to flat-field correction + GET_COLUMN_NEW(&headerCell, &ftableCell, R, "RA", double); + GET_COLUMN_NEW(&headerCell, &ftableCell, D, "DEC", double); + GET_COLUMN_NEW(&headerCell, &ftableCell, zone, "ZONE", int); + GET_COLUMN_NEW(&headerCell, &ftableCell, band, "BAND", int); + GET_COLUMN_NEW(&headerCell, &ftableCell, index, "INDEX", int); + GET_COLUMN_NEW(&headerCell, &ftableCell, Xo, "X_CENT", double); + GET_COLUMN_NEW(&headerCell, &ftableCell, Yo, "Y_CENT", double); + GET_COLUMN_NEW(&headerCell, &ftableCell, dX, "X_GRID", int); + GET_COLUMN_NEW(&headerCell, &ftableCell, dY, "Y_GRID", int); + GET_COLUMN_NEW(&headerCell, &ftableCell, name, "NAME", char); // XXX how is this done? + gfits_free_header (&headerCell); + gfits_free_table (&ftableCell); + + fprintf (stderr, "loaded data for %lld cells\n", (long long) Nrow); + + // assign the storage arrays + for (i = 0; i < Nrow; i++) { + nz = zone[i]; + nb = band[i]; + tree->ra[nz][nb] = R[i]; + tree->dec[nz][nb] = D[i]; + tree->Xo[nz][nb] = Xo[i]; + tree->Yo[nz][nb] = Yo[i]; + tree->dX[nz][nb] = dX[i]; + tree->dY[nz][nb] = dY[i]; + tree->cell[nz][nb] = i; // XXX ? + memcpy(tree->name[nz][nb], &name[i*BOUNDARY_TREE_NAME_LENGTH], BOUNDARY_TREE_NAME_LENGTH); + // XXX parse out the ID from the name (skycell.NNNN) + tree->projID[nz][nb] = atoi(&tree->name[nz][nb][8]); + } + + free (R ); + free (D ); + free (zone ); + free (band ); + free (Xo ); + free (Yo ); + free (dX ); + free (dY ); + free (index ); + free (name ); return tree; escape: - gfits_free_header (&header); - gfits_free_matrix (&matrix); - gfits_free_header (&theader); - gfits_free_table (&ftable); + gfits_free_header (&headerCell); + gfits_free_table (&ftableCell); + gfits_free_table (&ftableZone); if (tree) free (tree); - fclose (f); return NULL; } @@ -173,9 +196,6 @@ int BoundaryTreeSave(char *filename, BoundaryTree *tree) { - int i, nz, nb; Header header; - Header theader; Matrix matrix; - FTable ftable; gfits_init_header (&header); @@ -202,4 +222,16 @@ gfits_free_header (&header); gfits_free_matrix (&matrix); + + BoundaryTreeWrite (f, tree); + fclose (f); + + return TRUE; +} + +int BoundaryTreeWrite(FILE *f, BoundaryTree *tree) { + + int i, nz, nb; + Header theader; + FTable ftable; /*** zone information table ***/ Index: /branches/eam_branches/ipp-20131211/Ohana/src/libdvo/src/TessellationTable.c =================================================================== --- /branches/eam_branches/ipp-20131211/Ohana/src/libdvo/src/TessellationTable.c (revision 36471) +++ /branches/eam_branches/ipp-20131211/Ohana/src/libdvo/src/TessellationTable.c (revision 36472) @@ -11,7 +11,8 @@ # define DEBUG 0 -TessellationTable *TessellationTableLoad(char *filename) { - - int i, j, nz, nb, Ncol; +// backwards compatible load : if the first table is a BoundaryTree, not a TessellationTable, just load that. +TessellationTable *TessellationTableLoad(char *filename, int *Ntess) { + + int i, Ncol; off_t Nrow; char type[16]; @@ -26,4 +27,5 @@ theader.buffer = NULL; TessellationTable *tess = NULL; + *Ntess = 0; FILE *f = fopen (filename, "r"); @@ -43,107 +45,118 @@ } - ALLOCATE (tess, TessellationTable, 1); - ftable.header = &theader; - /*** zone information table ***/ - { - // load data for this header - if (!gfits_load_header (f, &theader)) goto escape; - - // read the fits table bytes - if (!gfits_fread_ftable_data (f, &ftable, FALSE)) goto escape; + // load data for this header + if (!gfits_load_header (f, &theader)) goto escape; + + // what kind of table have I just read (if any!) + char extname[64]; + gfits_scan (&theader, "EXTNAME", "%s", 1, extname); + + // is this a BoundaryTree? if so, it must be stand-alone: load just that and generate a + // containing TessellationTable. + if (!strcmp(extname, "ZONE_DATA")) { + BoundaryTree *tree = BoundaryTreeRead (&header, &theader, f); + gfits_free_header (&theader); + gfits_free_header (&header); + gfits_free_matrix (&matrix); + fclose (f); + + ALLOCATE (tess, TessellationTable, 1); + tess[0].tree = tree; + tess[0].type = TESS_RINGS; + *Ntess = 1; + return tess; + } + + if (strcmp(extname, "TESS_DATA")) { + fprintf (stderr, "header is neither ZONE_DATA nor TESS_DATA, problem with file\n"); + goto escape; + } + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) goto escape; - // need to create and assign to flat-field correction - GET_COLUMN_RAW(tess->Nband, "NBAND", int); - GET_COLUMN_RAW(tess->RA_origin, "RA_ORIGIN", double); - GET_COLUMN_RAW(tess->RA_offset, "RA_OFFSET", double); - GET_COLUMN_RAW(tess->DEC_min , "DEC_MIN", double); - GET_COLUMN_RAW(tess->DEC_max , "DEC_MAX", double); - GET_COLUMN_RAW(tess->DEC_min_raw, "DEC_MIN_RAW", double); - GET_COLUMN_RAW(tess->DEC_max_raw, "DEC_MAX_RAW", double); + // need to create and assign to flat-field correction + GET_COLUMN_NEW(Ro, "RA", double); + GET_COLUMN_NEW(Do, "DEC", double); + GET_COLUMN_NEW(Xo, "X_CENT", double); + GET_COLUMN_NEW(Yo, "Y_CENT", double); + GET_COLUMN_NEW(dX, "X_GRID", int); + GET_COLUMN_NEW(dY, "Y_GRID", int); + GET_COLUMN_NEW(Rmin, "R_MIN", double); + GET_COLUMN_NEW(Rmax, "R_MIN", double); + GET_COLUMN_NEW(Dmin, "D_MIN", double); + GET_COLUMN_NEW(Dmax, "D_MAX", double); + GET_COLUMN_NEW(dPix, "SCALE", double); + GET_COLUMN_NEW(NX_SUB, "NX_SUB", int); + GET_COLUMN_NEW(NY_SUB, "NY_SUB", int); + + *Ntess = Nrow; + ALLOCATE (tess, TessellationTable, *Ntess); + for (i = 0; i < *Ntess; i++) { + tess[i].Ro = Ro[i] ; + tess[i].Do = Do[i] ; + tess[i].Xo = Xo[i] ; + tess[i].Yo = Yo[i] ; + tess[i].dX = dX[i] ; + tess[i].dY = dY[i] ; + tess[i].Rmin = Rmin[i] ; + tess[i].Rmax = Rmax[i] ; + tess[i].Dmin = Dmin[i] ; + tess[i].Dmax = Dmax[i] ; + tess[i].dPix = dPix[i] ; + tess[i].NX_SUB = NX_SUB[i]; + tess[i].NY_SUB = NY_SUB[i]; + tess[i].type = TESS_LOCAL; + tess[i].tree = NULL; + } + + free(Ro ); + free(Do ); + free(Xo ); + free(Yo ); + free(dX ); + free(dY ); + free(Rmin ); + free(Rmax ); + free(Dmin ); + free(Dmax ); + free(dPix ); + free(NX_SUB); + free(NY_SUB); + + gfits_free_header (&theader); + gfits_free_table (&ftable); + + // fprintf (stderr, "loaded data for %lld projection cell\n", (long long) *Ntess); + + /*** check for a RINGS entry ***/ + + // load data for this header, or exit + if (!gfits_load_header (f, &theader)) { + return tess; + } + + // what kind of table have I just read (if any!) + gfits_scan (&theader, "EXTNAME", "%s", 1, extname); + + // is this a BoundaryTree? if so, it must be stand-alone: load just that and generate a + // containing TessellationTable. + if (!strcmp(extname, "ZONE_DATA")) { + BoundaryTree *tree = BoundaryTreeRead (&header, &theader, f); gfits_free_header (&theader); - gfits_free_table (&ftable); - - fprintf (stderr, "loaded data for %lld zones\n", (long long) Nrow); - tess->Nzone = Nrow; - - // allocate the storage arrays - ALLOCATE (tess->ra, double *, tess->Nzone); - ALLOCATE (tess->dec, double *, tess->Nzone); - ALLOCATE (tess->Xo, double *, tess->Nzone); - ALLOCATE (tess->Yo, double *, tess->Nzone); - ALLOCATE (tess->dX, int *, tess->Nzone); - ALLOCATE (tess->dY, int *, tess->Nzone); - ALLOCATE (tess->cell, int *, tess->Nzone); - ALLOCATE (tess->projID, int *, tess->Nzone); - ALLOCATE (tess->name, char **, tess->Nzone); - for (i = 0; i < tess->Nzone; i++) { - ALLOCATE (tess->ra[i], double, tess->Nband[i]); - ALLOCATE (tess->dec[i], double, tess->Nband[i]); - ALLOCATE (tess->Xo[i], double, tess->Nband[i]); - ALLOCATE (tess->Yo[i], double, tess->Nband[i]); - ALLOCATE (tess->dX[i], int, tess->Nband[i]); - ALLOCATE (tess->dY[i], int, tess->Nband[i]); - ALLOCATE (tess->cell[i], int, tess->Nband[i]); - ALLOCATE (tess->projID[i], int, tess->Nband[i]); - ALLOCATE (tess->name[i], char *, tess->Nband[i]); - for (j = 0; j < tess->Nband[i]; j++) { - ALLOCATE (tess->name[i][j], char, BOUNDARY_TESS_NAME_LENGTH); - } - } - } - - /*** cell information table ***/ - { - // load data for this header - if (!gfits_load_header (f, &theader)) goto escape; - - // read the fits table bytes - if (!gfits_fread_ftable_data (f, &ftable, FALSE)) goto escape; - - // need to create and assign to flat-field correction - GET_COLUMN_NEW(R, "RA", double); - GET_COLUMN_NEW(D, "DEC", double); - GET_COLUMN_NEW(zone, "ZONE", int); - GET_COLUMN_NEW(band, "BAND", int); - GET_COLUMN_NEW(index, "INDEX", int); - GET_COLUMN_NEW(Xo, "X_CENT", double); - GET_COLUMN_NEW(Yo, "Y_CENT", double); - GET_COLUMN_NEW(dX, "X_GRID", int); - GET_COLUMN_NEW(dY, "Y_GRID", int); - GET_COLUMN_NEW(name, "NAME", char); // XXX how is this done? - gfits_free_header (&theader); - gfits_free_table (&ftable); - fprintf (stderr, "loaded data for %lld cells\n", (long long) Nrow); - - // assign the storage arrays - for (i = 0; i < Nrow; i++) { - nz = zone[i]; - nb = band[i]; - tess->ra[nz][nb] = R[i]; - tess->dec[nz][nb] = D[i]; - tess->Xo[nz][nb] = Xo[i]; - tess->Yo[nz][nb] = Yo[i]; - tess->dX[nz][nb] = dX[i]; - tess->dY[nz][nb] = dY[i]; - tess->cell[nz][nb] = i; // XXX ? - memcpy(tess->name[nz][nb], &name[i*BOUNDARY_TESS_NAME_LENGTH], BOUNDARY_TESS_NAME_LENGTH); - // XXX parse out the ID from the name (skycell.NNNN) - tess->projID[nz][nb] = atoi(&tess->name[nz][nb][8]); - } - - free (R ); - free (D ); - free (zone ); - free (band ); - free (Xo ); - free (Yo ); - free (dX ); - free (dY ); - free (index ); - free (name ); - } - + gfits_free_header (&header); + gfits_free_matrix (&matrix); + fclose (f); + + (*Ntess) ++; + REALLOCATE (tess, TessellationTable, *Ntess); + tess[(*Ntess - 1)].tree = tree; + tess[(*Ntess - 1)].type = TESS_RINGS; + return tess; + } + + gfits_free_header (&theader); gfits_free_header (&header); gfits_free_matrix (&matrix); @@ -163,12 +176,11 @@ } -// we are passed a TessellationTable structure, write it to a FITS table (3 ext) -int TessellationTableSave(char *filename, TessellationTable *tess) { - - int i, nz, nb; +// we are passed a TessellationTable structure, write it to a FITS table. +// if one of the tess entries is a RINGS tessellation, write it as a boundary tree +int TessellationTableSave(char *filename, TessellationTable *tess, int Ntess) { + + int i; Header header; - Header theader; Matrix matrix; - FTable ftable; gfits_init_header (&header); @@ -183,11 +195,21 @@ } - // we need some information in the header to define the layout - gfits_modify (&header, "DEC_ORI", "%lf", 1, tess->DEC_origin); - gfits_modify (&header, "DEC_OFF", "%lf", 1, tess->DEC_offset); - - gfits_modify (&header, "NX_SUB", "%d", 1, tess->NX_SUB); - gfits_modify (&header, "NY_SUB", "%d", 1, tess->NY_SUB); - gfits_modify (&header, "PIXSCALE", "%lf", 1, tess->dPix); + // find the RINGS tessellation (if present) + BoundaryTree *rings = NULL; + for (i = 0; i < Ntess; i++) { + if (tess[i].type == TESS_LOCAL) continue; + if (tess[i].type != TESS_RINGS) continue; + rings = tess[i].tree; + } + + // if we have a RINGS entry, we need to save some information in the header to define the layout + if (rings) { + gfits_modify (&header, "DEC_ORI", "%lf", 1, rings->DEC_origin); + gfits_modify (&header, "DEC_OFF", "%lf", 1, rings->DEC_offset); + + gfits_modify (&header, "NX_SUB", "%d", 1, rings->NX_SUB); + gfits_modify (&header, "NY_SUB", "%d", 1, rings->NY_SUB); + gfits_modify (&header, "PIXSCALE", "%lf", 1, rings->dPix); + } gfits_fwrite_header (f, &header); @@ -196,134 +218,105 @@ gfits_free_matrix (&matrix); - /*** zone information table ***/ - { - gfits_create_table_header (&theader, "BINTABLE", "ZONE_DATA"); - - gfits_define_bintable_column (&theader, "J", "ZONE", "zone sequence number", "none", 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "NBAND", "number of cells in each zone", "none", 1.0, 0.0); - gfits_define_bintable_column (&theader, "D", "RA_ORIGIN", "origin of ra cell sequence", "degree", 1.0, 0.0); - gfits_define_bintable_column (&theader, "D", "RA_OFFSET", "offset per cell of ra cell sequence", "degree/cell", 1.0, 0.0); - gfits_define_bintable_column (&theader, "D", "DEC_MIN", "min dec for zone", "degree", 1.0, 0.0); - gfits_define_bintable_column (&theader, "D", "DEC_MAX", "max dec for zone", "degree", 1.0, 0.0); - gfits_define_bintable_column (&theader, "D", "DEC_MIN_RAW", "min dec for zone", "degree", 1.0, 0.0); - gfits_define_bintable_column (&theader, "D", "DEC_MAX_RAW", "max dec for zone", "degree", 1.0, 0.0); - - // generate the output array that carries the data - gfits_create_table (&theader, &ftable); - - // create intermediate storage arrays - int *zone = NULL; ALLOCATE (zone, int, tess->Nzone); - - // assign the storage arrays - for (i = 0; i < tess->Nzone; i++) { - zone[i] = i; - } - - // add the columns to the output array - gfits_set_bintable_column (&theader, &ftable, "ZONE", zone, tess->Nzone); - gfits_set_bintable_column (&theader, &ftable, "NBAND", tess->Nband, tess->Nzone); - gfits_set_bintable_column (&theader, &ftable, "RA_ORIGIN", tess->RA_origin, tess->Nzone); - gfits_set_bintable_column (&theader, &ftable, "RA_OFFSET", tess->RA_offset, tess->Nzone); - gfits_set_bintable_column (&theader, &ftable, "DEC_MIN", tess->DEC_min, tess->Nzone); - gfits_set_bintable_column (&theader, &ftable, "DEC_MAX", tess->DEC_max, tess->Nzone); - gfits_set_bintable_column (&theader, &ftable, "DEC_MIN_RAW", tess->DEC_min_raw, tess->Nzone); - gfits_set_bintable_column (&theader, &ftable, "DEC_MAX_RAW", tess->DEC_max_raw, tess->Nzone); - free (zone); - - gfits_fwrite_Theader (f, &theader); - gfits_fwrite_table (f, &ftable); - gfits_free_header (&theader); - gfits_free_table (&ftable); - } - - /*** cell information table ***/ - { - gfits_create_table_header (&theader, "BINTABLE", "CELL_DATA"); - - char fmt[16]; - snprintf (fmt, 16, "%dA", BOUNDARY_TESS_NAME_LENGTH); - gfits_define_bintable_column (&theader, "D", "RA", "ra (J2000) of cell center", "degree", 1.0, 0.0); - gfits_define_bintable_column (&theader, "D", "DEC", "dec (J2000) of cell center", "degree", 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "ZONE", "zone sequence number", "none", 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "BAND", "band sequence number", "none", 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "INDEX","cell index", "none", 1.0, 0.0); - gfits_define_bintable_column (&theader, "D", "X_CENT", "projection cell center pixel", "none", 1.0, 0.0); - gfits_define_bintable_column (&theader, "D", "Y_CENT", "projection cell center pixel", "none", 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "X_GRID", "skycell grid spacing", "none", 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "Y_GRID", "skycell grid spacing", "none", 1.0, 0.0); - gfits_define_bintable_column (&theader, fmt, "NAME", "cell name", "none", 1.0, 0.0); - - // generate the output array that carries the data - gfits_create_table (&theader, &ftable); - - int Ncell = 0; - for (i = 0; i < tess->Nzone; i++) { - Ncell += tess->Nband[i]; - } - - // create intermediate storage arrays - // NOTE: we have to unroll the 2D arrays in tess into 1D arrays - double *R ; ALLOCATE (R, double, Ncell); - double *D ; ALLOCATE (D, double, Ncell); - int *zone ; ALLOCATE (zone, int, Ncell); - int *band ; ALLOCATE (band, int, Ncell); - int *index ; ALLOCATE (index, int, Ncell); - double *Xo ; ALLOCATE (Xo, double, Ncell); - double *Yo ; ALLOCATE (Yo, double, Ncell); - int *dX ; ALLOCATE (dX, int, Ncell); - int *dY ; ALLOCATE (dY, int, Ncell); - char *name ; ALLOCATE (name, char, Ncell*BOUNDARY_TESS_NAME_LENGTH); - - // NOTE: a table column of characters must be fixed width, and is passed as a - // contiguous array of Nchar * Nrow values - - // assign the storage arrays - i = 0; - for (nz = 0; nz < tess->Nzone; nz++) { - for (nb = 0; nb < tess->Nband[nz]; nb++) { - R[i] = tess->ra[nz][nb]; - D[i] = tess->dec[nz][nb]; - Xo[i] = tess->Xo[nz][nb]; - Yo[i] = tess->Yo[nz][nb]; - dX[i] = tess->dX[nz][nb]; - dY[i] = tess->dY[nz][nb]; - zone[i] = nz; - band[i] = nb; - index[i] = i; // or tess->cells[nz][nb] ? - memcpy(&name[i*BOUNDARY_TESS_NAME_LENGTH], tess->name[nz][nb], BOUNDARY_TESS_NAME_LENGTH); - i++; - } - } - - // add the columns to the output array - gfits_set_bintable_column (&theader, &ftable, "RA", R, Ncell); - gfits_set_bintable_column (&theader, &ftable, "DEC", D, Ncell); - gfits_set_bintable_column (&theader, &ftable, "ZONE", zone, Ncell); - gfits_set_bintable_column (&theader, &ftable, "BAND", band, Ncell); - gfits_set_bintable_column (&theader, &ftable, "INDEX", index, Ncell); - gfits_set_bintable_column (&theader, &ftable, "X_CENT", Xo, Ncell); - gfits_set_bintable_column (&theader, &ftable, "Y_CENT", Yo, Ncell); - gfits_set_bintable_column (&theader, &ftable, "X_GRID", dX, Ncell); - gfits_set_bintable_column (&theader, &ftable, "Y_GRID", dY, Ncell); - gfits_set_bintable_column (&theader, &ftable, "NAME", name, Ncell); - - free (R ); - free (D ); - free (zone ); - free (band ); - free (index ); - free (Xo ); - free (Yo ); - free (dX ); - free (dY ); - free (name ); - - gfits_fwrite_Theader (f, &theader); - gfits_fwrite_table (f, &ftable); - gfits_free_header (&theader); - gfits_free_table (&ftable); - } + /** LOCAL tessellation table info **/ + Header theader; + FTable ftable; + + gfits_create_table_header (&theader, "BINTABLE", "TESS_DATA"); + + gfits_define_bintable_column (&theader, "D", "RA", "ra (J2000) of cell center", "degree", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "DEC", "dec (J2000) of cell center", "degree", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "X_CENT", "projection cell center pixel", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "Y_CENT", "projection cell center pixel", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "X_GRID", "skycell grid spacing", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "Y_GRID", "skycell grid spacing", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "R_MIN", "RA limit (lower)", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "R_MAX", "RA limit (upper)", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "D_MIN", "DEC limit (lower)", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "D_MAX", "DEC limit (upper)", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "SCALE", "pixel scale for projection cell", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "NX_SUB", "skycell subdivision in x", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "NY_SUB", "skycell subdivision in y", "none", 1.0, 0.0); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + int Nout = (rings == NULL) ? Ntess : Ntess - 1; + + // create intermediate storage arrays + // NOTE: we have to unroll the 2D arrays in tess into 1D arrays + double *Ro ; ALLOCATE (Ro , double, Nout); + double *Do ; ALLOCATE (Do , double, Nout); + double *Xo ; ALLOCATE (Xo , double, Nout); + double *Yo ; ALLOCATE (Yo , double, Nout); + int *dX ; ALLOCATE (dX , int, Nout); + int *dY ; ALLOCATE (dY , int, Nout); + double *Rmin ; ALLOCATE (Rmin , double, Nout); + double *Rmax ; ALLOCATE (Rmax , double, Nout); + double *Dmin ; ALLOCATE (Dmin , double, Nout); + double *Dmax ; ALLOCATE (Dmax , double, Nout); + double *dPix ; ALLOCATE (dPix , double, Nout); + int *NX_SUB ; ALLOCATE (NX_SUB, int, Nout); + int *NY_SUB ; ALLOCATE (NY_SUB, int, Nout); + + // assign the storage arrays + int N = 0; + for (i = 0; i < Ntess; i++) { + if (tess->type != TESS_LOCAL) continue; + Ro[N] = tess[i].Ro; + Do[N] = tess[i].Do; + Xo[N] = tess[i].Xo; + Yo[N] = tess[i].Yo; + dX[N] = tess[i].dX; + dY[N] = tess[i].dY; + Rmin[N] = tess[i].Rmin; + Rmax[N] = tess[i].Rmax; + Dmin[N] = tess[i].Dmin; + Dmax[N] = tess[i].Dmax; + dPix[N] = tess[i].dPix; + NX_SUB[N] = tess[i].NX_SUB; + NY_SUB[N] = tess[i].NY_SUB; + N++; + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "RA", Ro, Nout); + gfits_set_bintable_column (&theader, &ftable, "DEC", Do, Nout); + gfits_set_bintable_column (&theader, &ftable, "X_CENT", Xo, Nout); + gfits_set_bintable_column (&theader, &ftable, "Y_CENT", Yo, Nout); + gfits_set_bintable_column (&theader, &ftable, "X_GRID", dX, Nout); + gfits_set_bintable_column (&theader, &ftable, "Y_GRID", dY, Nout); + gfits_set_bintable_column (&theader, &ftable, "R_MIN", Rmin, Nout); + gfits_set_bintable_column (&theader, &ftable, "R_MAX", Rmax, Nout); + gfits_set_bintable_column (&theader, &ftable, "D_MIN", Dmin, Nout); + gfits_set_bintable_column (&theader, &ftable, "D_MAX", Dmax, Nout); + gfits_set_bintable_column (&theader, &ftable, "SCALE", dPix, Nout); + gfits_set_bintable_column (&theader, &ftable, "NX_SUB", NX_SUB, Nout); + gfits_set_bintable_column (&theader, &ftable, "NY_SUB", NY_SUB, Nout); + + free (Ro ); + free (Do ); + free (Xo ); + free (Yo ); + free (dX ); + free (dY ); + free (Rmin ); + free (Rmax ); + free (Dmin ); + free (Dmax ); + free (dPix ); + free (NX_SUB); + free (NY_SUB); + + gfits_fwrite_Theader (f, &theader); + gfits_fwrite_table (f, &ftable); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + if (rings) { + BoundaryTreeWrite(f, rings); + } + + fclose (f); + return TRUE; } - -