IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 36471


Ignore:
Timestamp:
Feb 3, 2014, 9:25:52 AM (12 years ago)
Author:
eugene
Message:

working on adding local tessellations to the boundary stuff

Location:
branches/eam_branches/ipp-20131211/Ohana/src
Files:
1 added
2 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20131211/Ohana/src/addstar/src/findskycell.c

    r36468 r36471  
    1212// in an even more specific case, RA[i,zone] = RA_origin[zone] + RA_offset[zone]
    1313
    14 enum {TREE_NONE, TREE_MAKE, TREE_USE};
     14enum {TREE_NONE, TREE_MAKE, TREE_LOCAL, TREE_USE};
     15enum {REGION_NONE, REGION_USER, REGION_MIN, REGION_MAX};
    1516
    1617void usage (void) {
    1718  fprintf (stderr, "USAGE: findcell -mktree (tree) (catdir) [-nx Nx] [-ny Ny]\n");
     19  fprintf (stderr, "USAGE: findcell -mklocal (file) (catdir) [-nx Nx] [-ny Ny]\n");
    1820  fprintf (stderr, "USAGE: findcell -tree (tree) (datafile)\n");
    1921  fprintf (stderr, "   (datafile) should contain a list of RA,DEC pairs\n");
    … …  
    2830int NY_SUB = 1;
    2931
     32double R_MIN = NAN;
     33double R_MAX = NAN;
     34double D_MIN = NAN;
     35double D_MAX = NAN;
     36int REGION_OPTION = REGION_NONE;
     37
     38char *BASENAME = NULL;
     39int projectIDoff = -1;
     40int skycellIDoff = -1;
     41
    3042int main (int argc, char **argv) {
    3143
    … …  
    5971    remove_argument (N, &argc, argv);
    6072    SCALE = atof (argv[N]);
     73    remove_argument (N, &argc, argv);
     74  }
     75
     76  // user-specified region
     77  if ((N = get_argument (argc, argv, "-region"))) {
     78    remove_argument (N, &argc, argv);
     79    if (N > argc - 4) {
     80      fprintf (stderr, "USAGE: -region requires 4 arguments (Rmin Rmax Dmin Dmax)\n");
     81      exit (1);
     82    }
     83    R_MIN = atof (argv[N]); remove_argument (N, &argc, argv);
     84    R_MAX = atof (argv[N]); remove_argument (N, &argc, argv);
     85    D_MIN = atof (argv[N]); remove_argument (N, &argc, argv);
     86    D_MAX = atof (argv[N]); remove_argument (N, &argc, argv);
     87    REGION_OPTION = REGION_USER;
     88  }
     89
     90  if ((N = get_argument (argc, argv, "-region-min"))) {
     91    remove_argument (N, &argc, argv);
     92    REGION_OPTION = REGION_MIN;
     93  }
     94
     95  if ((N = get_argument (argc, argv, "-region-max"))) {
     96    remove_argument (N, &argc, argv);
     97    REGION_OPTION = REGION_MAX;
     98  }
     99
     100  if ((N = get_argument (argc, argv, "-basename"))) {
     101    remove_argument (N, &argc, argv);
     102    if (N > argc - 3) {
     103      fprintf (stderr, "USAGE: -basename (name) (proj_offset) (skycell_offset)\n");
     104      exit (1);
     105    }
     106    BASENAME = strcreate  (argv[N]);
     107    remove_argument (N, &argc, argv);
     108    projectIDoff = atoi(argv[N]);
     109    remove_argument (N, &argc, argv);
     110    skycellIDoff = atoi(argv[N]);
    61111    remove_argument (N, &argc, argv);
    62112  }
    … …  
    70120    remove_argument (N, &argc, argv);
    71121  }
     122  if ((N = get_argument (argc, argv, "-mklocal"))) {
     123    MODE = TREE_LOCAL;
     124    remove_argument (N, &argc, argv);
     125    treefile = strcreate (argv[N]);
     126    remove_argument (N, &argc, argv);
     127  }
    72128  if ((N = get_argument (argc, argv, "-tree"))) {
    73129    MODE = TREE_USE;
    … …  
    83139  if (MODE == TREE_MAKE) {
    84140    mktree (treefile, argv[1]);
     141    exit (0);
     142  }
     143
     144  if (MODE == TREE_LOCAL) {
     145    mklocal (treefile, argv[1]);
    85146    exit (0);
    86147  }
    … …  
    408469// this function takes a catdir and generates an extension for a tess file for a local tess.
    409470
     471/*
     472
     473  LOCAL projections are defined by single projection cells.  should I be supplying the info on
     474  the cmd line or figure out the bounds from the catdir? 
     475
     476 */
     477
    410478int mklocal (char *treefile, char *catdir) {
    411479
    … …  
    416484  double x, y, ra, dec;
    417485
     486  if (REGION_OPTION == REGION_NONE) {
     487    fprintf (stderr, "ERROR: need to define bounding region (-region Rmin Rmax Dmin Dmax | -region-min | -region-max)\n");
     488    exit (2);
     489  }
     490
    418491  char imagefile[DVO_MAX_PATH];
    419492  snprintf (imagefile, DVO_MAX_PATH, "%s/Images.dat", catdir);
    … …  
    431504  image = gfits_table_get_Image (&db.ftable, &Nimage, &db.swapped);
    432505  if (!image) {
    433       fprintf (stderr, "ERROR: failed to read images\n");
    434       exit (2);
     506    fprintf (stderr, "ERROR: failed to read images\n");
     507    exit (2);
    435508  }
    436509 
    437510  // generate an empty BoundaryTree
    438   TessallationBoundary tess;
    439  
    440   // user supplied values, do not try to derive from Image.dat
    441   tree.NX_SUB = NX_SUB;
    442   tree.NY_SUB = NY_SUB;
    443   tree.dPix = SCALE/3600.0;
    444 
    445   // derive these or assign these?
    446   tree.Dmin = NAN;
    447   tree.Dmax = NAN;
    448   tree.Rmin = NAN;
    449   tree.Rmax = NAN;
    450 
     511  TessallationBoundary *tess = NULL;
     512  ALLOCATE (tess, TessallationBoundary, Nimage);
     513 
    451514  // find the RA,DEC of the image centers & assign to cells
    452515  for (i = 0; i < Nimage; i++) {
     516    // user supplied values, do not try to derive from Image.dat
     517    tess[i].NX_SUB = NX_SUB;
     518    tess[i].NY_SUB = NY_SUB;
     519    tess[i].dPix = SCALE/3600.0;
     520
    453521    x = 0.5*image[i].NX;
    454522    y = 0.5*image[i].NY;
    455523    XY_to_RD (&ra, &dec, x, y, &image[i].coords);
    456524
    457     if (!BoundaryTreeCellCoords (&tree, &zone, &band, ra, dec)) {
    458       fprintf (stderr, "mismatch!\n");
    459       continue;
    460     }
    461     // fprintf (stderr, "%d  %f %f  %f  %f %f  %d %d\n", i, x, y, tree.RA_offset[zone], ra, dec, zone, band);
    462    
    463     if (band >= tree.NBAND[zone]) {
    464       int start = tree.NBAND[zone];
    465       tree.NBAND[zone] = band + 10;
    466       REALLOCATE (tree.ra[zone],   double, tree.NBAND[zone]);
    467       REALLOCATE (tree.dec[zone],  double, tree.NBAND[zone]);
    468       REALLOCATE (tree.Xo[zone],   double, tree.NBAND[zone]);
    469       REALLOCATE (tree.Yo[zone],   double, tree.NBAND[zone]);
    470       REALLOCATE (tree.dX[zone],      int, tree.NBAND[zone]);
    471       REALLOCATE (tree.dY[zone],      int, tree.NBAND[zone]);
    472       REALLOCATE (tree.cell[zone],    int, tree.NBAND[zone]);
    473       REALLOCATE (tree.name[zone], char *, tree.NBAND[zone]);
    474       for (j = start; j < tree.NBAND[zone]; j++) {
    475         tree.ra[zone][j] = NAN;
    476         tree.dec[zone][j] = NAN;
    477         tree.Xo[zone][band] = NAN;
    478         tree.Yo[zone][band] = NAN;
    479         tree.dX[zone][band] = -1;
    480         tree.dY[zone][band] = -1;
    481         tree.cell[zone][j] = -1;
    482         ALLOCATE (tree.name[zone][j], char, BOUNDARY_TREE_NAME_LENGTH);
    483       }
    484     }
    485     tree.ra[zone][band] = ra;
    486     tree.dec[zone][band] = dec;
    487     tree.Xo[zone][band] = x;
    488     tree.Yo[zone][band] = y;
    489     tree.dX[zone][band] = image[i].NX / NX_SUB;
    490     tree.dY[zone][band] = image[i].NY / NY_SUB;
    491    
    492     tree.cell[zone][band] = i;
    493 
    494     // what are the min and max DEC values for this zone? (test center and corners of top and bottom edge) 
    495 
    496     if (dec > 0.0) {
    497       // min DEC at bottom of the cell at the center
    498       x = 0.5*image[i].NX;
    499       y = 0.0*image[i].NY;
    500       XY_to_RD (&ra, &dec, x, y, &image[i].coords);
    501       tree.DEC_min_raw[zone] = MIN(tree.DEC_min_raw[zone], dec);
    502  
    503       // max DEC : find the intersection between the RA boundary and the top of the cell
    504       double ra_band_min = tree.RA_origin[zone] + tree.RA_offset[zone] * band;
    505  
    506       // does the parity matter?
    507       x = 0.0*image[i].NX;
    508       y = 1.0*image[i].NY;
    509  
    510       int Niter;
    511       for (Niter = 0; Niter < 3; Niter ++) {
    512         XY_to_RD (&ra, &dec, x, y, &image[i].coords);
    513         RD_to_XY (&x, &y, ra_band_min, dec, &image[i].coords);
    514         y = 1.0*image[i].NY;
    515       }
    516       tree.DEC_max_raw[zone] = MAX(tree.DEC_max_raw[zone], dec);
    517     } else {
    518       // max DEC at top of the cell at the center
    519       x = 0.5*image[i].NX;
    520       y = 1.0*image[i].NY;
    521       XY_to_RD (&ra, &dec, x, y, &image[i].coords);
    522       tree.DEC_max_raw[zone] = MAX(tree.DEC_max_raw[zone], dec);
    523  
    524       // max DEC : find the intersection between the RA boundary and the bottom of the cell
    525       double ra_band_min = tree.RA_origin[zone] + tree.RA_offset[zone] * band;
    526  
    527       // does the parity matter? (NO)
    528       x = 0.0*image[i].NX;
    529       y = 0.0*image[i].NY;
    530  
    531       int Niter;
    532       for (Niter = 0; Niter < 3; Niter ++) {
    533         XY_to_RD (&ra, &dec, x, y, &image[i].coords);
    534         RD_to_XY (&x, &y, ra_band_min, dec, &image[i].coords);
    535         y = 0.0*image[i].NY;
    536       }
    537       tree.DEC_min_raw[zone] = MIN(tree.DEC_min_raw[zone], dec);
    538     }
    539     memcpy (tree.name[zone][band], image[i].name, BOUNDARY_TREE_NAME_LENGTH);
    540   }
    541 
    542   // figure out the max band value for each zone?
    543   for (zone = 0; zone < tree.Nzone; zone++) {
    544     int found_last = FALSE;
    545     int last_band = -1;
    546     for (band = 0; band < tree.NBAND[zone]; band++) {
    547       // all cells should be filled
    548       if (tree.cell[zone][band] < 0) {
    549         if (!found_last) {
    550           found_last = TRUE;
    551           last_band = band;
     525    tess[i].ra  = ra;
     526    tess[i].dec = dec;
     527    tess[i].Xo  = x;
     528    tess[i].Yo  = y;
     529    tess[i].dX  = image[i].NX / NX_SUB;
     530    tess[i].dY  = image[i].NY / NY_SUB;
     531
     532    // find the minimum or maximum containing region
     533    if ((REGION_OPTION == REGION_MIN) || (REGION_OPTION == REGION_MAX)) {
     534      // XXX short cut for now : if the projection cell bounds the equator or 0,360 boundary, this will fail:
     535      double R[4], D[4];
     536      XY_to_RD (&R[0], &D[0],           0,           0, &image[i].coords);
     537      XY_to_RD (&R[1], &D[1], image[i].NX,           0, &image[i].coords);
     538      XY_to_RD (&R[2], &D[2],           0, image[i].NY, &image[i].coords);
     539      XY_to_RD (&R[3], &D[3], image[i].NX, image[i].NY, &image[i].coords);
     540
     541      if (REGION_OPTION == REGION_MIN) {
     542        for (i = 0; i < 4; i++) {
     543          R_MIN = (R[i] < ra)  ? (isfinite(R_MIN) ? MAX(R_MIN, R[i]) : R[i]) : R_MIN;
     544          R_MAX = (R[i] > ra)  ? (isfinite(R_MAX) ? MIN(R_MAX, R[i]) : R[i]) : R_MAX;
     545          D_MIN = (D[i] < dec) ? (isfinite(D_MIN) ? MAX(D_MIN, D[i]) : D[i]) : D_MIN;
     546          D_MAX = (D[i] > dec) ? (isfinite(D_MAX) ? MIN(D_MAX, D[i]) : D[i]) : D_MAX;
    552547        }
    553548      } else {
    554         if (found_last) {
    555           fprintf (stderr, "error: empty cell (%d,%d) after last band (%d)\n", zone, band, last_band);
     549        for (i = 0; i < 4; i++) {
     550          R_MIN = (R[i] < ra)  ? (isfinite(R_MIN) ? MIN(R_MIN, R[i]) : R[i]) : R_MIN;
     551          R_MAX = (R[i] > ra)  ? (isfinite(R_MAX) ? MAX(R_MAX, R[i]) : R[i]) : R_MAX;
     552          D_MIN = (D[i] < dec) ? (isfinite(D_MIN) ? MIN(D_MIN, D[i]) : D[i]) : D_MIN;
     553          D_MAX = (D[i] > dec) ? (isfinite(D_MAX) ? MAX(D_MAX, D[i]) : D[i]) : D_MAX;
    556554        }
    557555      }
    558556    }
    559     if (last_band == -1) {
    560       last_band = tree.NBAND[zone];
    561     }
    562     tree.Nband[zone] = last_band;
    563     // fprintf (stderr, "last_band : %d, Nband: %d\n", last_band, tree.Nband[zone]);
    564   }
    565 
    566   // figure out the max band value for each zone?
    567   int Nm = 0;
    568   tree.DEC_min[Nm] = -90.0;
    569   tree.DEC_max[Nm] = 0.5*(tree.DEC_max_raw[Nm] + tree.DEC_min_raw[Nm + 1]);
    570 
    571   int Np = tree.Nzone - 1;
    572   tree.DEC_min[Np] = 0.5*(tree.DEC_min_raw[Np] + tree.DEC_max_raw[Np - 1]);
    573   tree.DEC_max[Np] = +90.0;
    574 
    575   for (zone = 1; zone < tree.Nzone - 1; zone++) {
    576     tree.DEC_min[zone] = 0.5*(tree.DEC_min_raw[zone] + tree.DEC_max_raw[zone - 1]);
    577     tree.DEC_max[zone] = 0.5*(tree.DEC_max_raw[zone] + tree.DEC_min_raw[zone + 1]);
    578   }
    579 
    580   struct timeval start, stop;
    581   gettimeofday (&start, (void *) NULL);
    582 
    583   int Npts = 10000000;
    584 
    585   // test : find skycell for NN random points on the sky
    586   long A = time(NULL);
    587   long B = A + 10000;
    588   srand48(B);
    589   for (i = 0; i < Npts; i++) {
    590     ra  = 360.0 * drand48();
    591     dec = 180.0 * drand48() - 90.0;
    592     if (!BoundaryTreeCellCoords (&tree, &zone, &band, ra, dec)) {
    593       fprintf (stderr, "failure for %f,%f\n", ra, dec);
    594     }
    595   }
    596   MARKTIME("-- test %d pts: %f sec\n", Npts, dtime);
     557    tess[i].Rmin = R_MIN;
     558    tess[i].Rmax = R_MAX;
     559    tess[i].Dmin = D_MIN;
     560    tess[i].Dmax = D_MAX;
     561
     562    tess[i].type = TESS_LOCAL;
     563
     564    // XXX I don't really want to do the work of discovering the rule...
     565    memcpy (tess[i].basename, BASENAME);
     566    tess[i].Nbasename = strlen(BASENAME);
     567
     568    tess[i].projectIDoff = projectIDoff;
     569    tess[i].skycellIDoff = skycellIDoff;
     570  }
    597571
    598572  BoundaryTreeSave (treefile, &tree);

  • branches/eam_branches/ipp-20131211/Ohana/src/libdvo/include/dvo.h

    r36465 r36471  
    346346typedef enum { TESS_NONE, TESS_LOCAL, TESS_RINGS } TessType;
    347347
    348 // TessallationBoundaries is a structure to describe the LOCAL skycell boundaries in terms of lines of constant (RA,DEC)
    349 // the structure describes the boundaries of a SINGLE projection cell with Nx * Ny skycells
     348// TessallationTable is a structure to describe the parameters of a set of "tessellations"
     349// (these are not strictly tessellations but projection sets as only the non-local
     350// versions can cover the full sky).  For LOCAL projection cells, the structure describes
     351// the boundaries of a SINGLE projection cell with Nx * Ny skycells and includes some
     352// basic parameters (not used by the fullsky, eg RINGS, tessellations)
    350353typedef struct {
    351354  double Rmin; // this tessellation is valid only for RA >= Rmin
    … …  
    365368  int NY_SUB;
    366369
    367   TessType type;
     370  TessType type; //
    368371  BoundaryTree *tree;
    369 } TessallationBoundary;
     372} TessallationTable;
    370373
    371374// a reduced-subset structure for relphot
Note: See TracChangeset for help on using the changeset viewer.