IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 40621


Ignore:
Timestamp:
Feb 8, 2019, 8:19:26 PM (7 years ago)
Author:
eugene
Message:

add CFIS-style skycell to findskycell.c

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/addstar/src/findskycell.c

    r38441 r40621  
    1212// in an even more specific case, RA[i,zone] = RA_origin[zone] + RA_offset[zone]
    1313
    14 enum {TREE_NONE, TREE_MAKE, TREE_LOCAL, TREE_USE};
     14enum {TREE_NONE, TREE_MAKE, TREE_CFIS, TREE_LOCAL, TREE_USE};
    1515enum {REGION_NONE, REGION_USER, REGION_MIN, REGION_MAX};
    1616
    … …  
    2424
    2525int mktree (char *treefile, char *catdir);
     26int mkcfis (char *treefile, char *catdir);
    2627int mklocal (char *treefile, char *catdir);
    2728int apply_tree (char *treefile, char *datafile);
    … …  
    131132    remove_argument (N, &argc, argv);
    132133  }
     134  if ((N = get_argument (argc, argv, "-mkcfis"))) {
     135    MODE = TREE_CFIS;
     136    remove_argument (N, &argc, argv);
     137    treefile = strcreate (argv[N]);
     138    remove_argument (N, &argc, argv);
     139  }
    133140  if ((N = get_argument (argc, argv, "-mklocal"))) {
    134141    MODE = TREE_LOCAL;
    … …  
    153160  }
    154161
     162  if (MODE == TREE_CFIS) {
     163    mkcfis (treefile, argv[1]);
     164    exit (0);
     165  }
     166
    155167  if (MODE == TREE_LOCAL) {
    156168    mklocal (treefile, argv[1]);
    … …  
    162174}
    163175
    164 int mktree (char *treefile, char *catdir) {
     176int mkcfis (char *treefile, char *catdir) {
    165177
    166178  int i, j, zone, band, status;
    … …  
    195207  tree.FixedGridRA = TRUE;
    196208 
     209  float CELLSIZE = 0.5;
     210
    197211  // for the moment, I'm going to hardwire the DEC bands to match RINGS.V3
    198   tree.DEC_origin = -92.0;
    199   tree.DEC_offset =   4.0;
    200 
    201   tree.Nzone = 46;
     212  tree.DEC_origin = -90.0 - 0.5*CELLSIZE;
     213  tree.DEC_offset = CELLSIZE;
     214
     215  tree.Nzone = (180 / CELLSIZE) + 1;
    202216  tree.NX_SUB = NX_SUB;
    203217  tree.NY_SUB = NY_SUB;
    204   tree.dPix = SCALE/3600.0;
     218  tree.dPix = SCALE/3600.0; // user-supplied (0.185768447409 for CFIS)
    205219
    206220  ALLOCATE (tree.Nband, int, tree.Nzone);
    … …  
    238252    tree.Nband[zone] = 1000000;
    239253    tree.NBAND[zone] = 10;
    240     tree.RA_origin[zone] = -0.5*RA_offset_RINGS_V3[zone];
    241     tree.RA_offset[zone] = RA_offset_RINGS_V3[zone];
     254
     255    float zoneDec = -90.0 + zone * CELLSIZE;
     256    float dRA = CELLSIZE / cos(RAD_DEG*zoneDec); // dRA is a size in RA degrees == \alpha_n
     257
     258    tree.RA_origin[zone] = 0.0 - 0.5*dRA;
     259    tree.RA_offset[zone] = dRA;
    242260
    243261    // set the starting values here.  In fact, the MIN and MAX values need to be tweaked
    … …  
    414432  MARKTIME("-- test %d pts: %f sec\n", Npts, dtime);
    415433
     434  int Ntess = 0;
     435  TessellationTable *tess = NULL;
     436
     437  if (APPEND) {
     438    tess = TessellationTableLoad (treefile, &Ntess);
     439    REALLOCATE (tess, TessellationTable, Ntess + 1);
     440  } else {
     441    ALLOCATE (tess, TessellationTable, Ntess + 1);
     442  }   
     443
     444  TessellationTableInit (&tess[Ntess], 1);
     445  tess[Ntess].tree = &tree;
     446  tess[Ntess].type = TESS_RINGS;
     447  tess[Ntess].Rmin =   0;
     448  tess[Ntess].Rmax = 360;
     449  tess[Ntess].Dmin = -90;
     450  tess[Ntess].Dmax = +90;
     451
     452  if (BASENAME) {
     453    tess[Ntess].Nbasename = strlen(BASENAME);
     454    tess[Ntess].basename = strcreate(BASENAME);
     455    tess[Ntess].projectIDoff = projectIDoff;
     456    tess[Ntess].skycellIDoff = skycellIDoff;
     457  }
     458     
     459  Ntess ++;
     460  TessellationTableSave (treefile, tess, Ntess);
     461
     462  return TRUE;
     463}
     464
     465int mktree (char *treefile, char *catdir) {
     466
     467  int i, j, zone, band, status;
     468  FITS_DB db;
     469  Image *image;
     470  off_t Nimage;
     471  double x, y, ra, dec;
     472
     473  char imagefile[DVO_MAX_PATH];
     474  snprintf (imagefile, DVO_MAX_PATH, "%s/Images.dat", catdir);
     475
     476  gfits_db_init (&db);
     477  status = dvo_image_lock (&db, imagefile, 2.0, LCK_XCLD);
     478  if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db.filename);
     479
     480  /* load or create the image table */
     481  if (db.dbstate == LCK_EMPTY) Shutdown ("can't read image catalog %s", db.filename);
     482
     483  if (!dvo_image_load (&db, TRUE, FALSE)) Shutdown ("can't read image catalog %s", db.filename);
     484
     485  // convert database table to internal structure (binary to Image)
     486  // 'image' points to the same memory as db->ftable->buffer
     487  image = gfits_table_get_Image (&db.ftable, &Nimage, &db.scaledValue, &db.nativeOrder);
     488  if (!image) {
     489      fprintf (stderr, "ERROR: failed to read images\n");
     490      exit (2);
     491  }
     492 
     493  // generate an empty BoundaryTree
     494  BoundaryTree tree;
     495  tree.FixedGridDEC = TRUE;
     496  tree.FixedGridRA = TRUE;
     497 
     498  // for the moment, I'm going to hardwire the DEC bands to match RINGS.V3
     499  tree.DEC_origin = -92.0;
     500  tree.DEC_offset =   4.0;
     501
     502  tree.Nzone = 46;
     503  tree.NX_SUB = NX_SUB;
     504  tree.NY_SUB = NY_SUB;
     505  tree.dPix = SCALE/3600.0;
     506
     507  ALLOCATE (tree.Nband, int, tree.Nzone);
     508  ALLOCATE (tree.NBAND, int, tree.Nzone);
     509
     510  ALLOCATE (tree.RA_origin,   double, tree.Nzone);
     511  ALLOCATE (tree.RA_offset,   double, tree.Nzone);
     512  ALLOCATE (tree.DEC_min,     double, tree.Nzone);
     513  ALLOCATE (tree.DEC_max,     double, tree.Nzone);
     514  ALLOCATE (tree.DEC_min_raw, double, tree.Nzone);
     515  ALLOCATE (tree.DEC_max_raw, double, tree.Nzone);
     516
     517  ALLOCATE (tree.ra,   double *, tree.Nzone);
     518  ALLOCATE (tree.dec,  double *, tree.Nzone);
     519  ALLOCATE (tree.Xo,   double *, tree.Nzone);
     520  ALLOCATE (tree.Yo,   double *, tree.Nzone);
     521  ALLOCATE (tree.dX,    float *, tree.Nzone);
     522  ALLOCATE (tree.dY,    float *, tree.Nzone);
     523  ALLOCATE (tree.cell,    int *, tree.Nzone);
     524  ALLOCATE (tree.name,  char **, tree.Nzone);
     525
     526  // NOTE 1: the RA_origin, RA_offsets for RINGS.V3 are defined so that the first projection cell
     527  // overlaps the RA = 0,360 boundary.  This make the split a bit of a hack.  We end up
     528  // with Nbands, but the max boundary of the last band only goes to 360.0 -
     529  // 0.5*RA_offset.  To get the right band number for the boundary region, if the
     530  // calculation for the band number lands beyond the Nbands (ie, band >= Nbands), then
     531  // we need to loop back to the first band.
     532
     533  // NOTE 2: when we generate the zone & bands initially, we do not know the number of
     534  // bands in the end.  we cannot use the test of band >= Nband unless we set an absurdly
     535  // large default value.  Thus the value of 1000000 below.
     536
     537  // assign the bands for RINGS.V3
     538  for (zone = 0; zone < tree.Nzone; zone++) {
     539    tree.Nband[zone] = 1000000;
     540    tree.NBAND[zone] = 10;
     541    tree.RA_origin[zone] = -0.5*RA_offset_RINGS_V3[zone];
     542    tree.RA_offset[zone] = RA_offset_RINGS_V3[zone];
     543
     544    // set the starting values here.  In fact, the MIN and MAX values need to be tweaked
     545    // on the basis of the real images...
     546    tree.DEC_min_raw[zone] = +90.0;
     547    tree.DEC_max_raw[zone] = -90.0;
     548    tree.DEC_min[zone] = NAN;
     549    tree.DEC_max[zone] = NAN;
     550
     551    ALLOCATE (tree.ra[zone],   double, tree.NBAND[zone]);
     552    ALLOCATE (tree.dec[zone],  double, tree.NBAND[zone]);
     553    ALLOCATE (tree.Xo[zone],   double, tree.NBAND[zone]);
     554    ALLOCATE (tree.Yo[zone],   double, tree.NBAND[zone]);
     555    ALLOCATE (tree.dX[zone],    float, tree.NBAND[zone]);
     556    ALLOCATE (tree.dY[zone],    float, tree.NBAND[zone]);
     557    ALLOCATE (tree.cell[zone], int,    tree.NBAND[zone]);
     558    ALLOCATE (tree.name[zone], char *, tree.NBAND[zone]);
     559    for (band = 0; band < tree.NBAND[zone]; band++) {
     560      tree.ra[zone][band] = NAN;
     561      tree.dec[zone][band] = NAN;
     562      tree.Xo[zone][band] = NAN;
     563      tree.Yo[zone][band] = NAN;
     564      tree.dX[zone][band] = NAN;
     565      tree.dY[zone][band] = NAN;
     566      tree.cell[zone][band] = -1;
     567      ALLOCATE (tree.name[zone][band], char, BOUNDARY_TREE_NAME_LENGTH);
     568    }
     569  }
     570
     571  // find the RA,DEC of the image centers & assign to cells
     572  for (i = 0; i < Nimage; i++) {
     573    x = 0.5*image[i].NX;
     574    y = 0.5*image[i].NY;
     575    XY_to_RD (&ra, &dec, x, y, &image[i].coords);
     576
     577    if (!BoundaryTreeCellCoords (&tree, &zone, &band, ra, dec)) {
     578      fprintf (stderr, "mismatch!\n");
     579      continue;
     580    }
     581    // fprintf (stderr, "%d  %f %f  %f  %f %f  %d %d\n", i, x, y, tree.RA_offset[zone], ra, dec, zone, band);
     582   
     583    if (band >= tree.NBAND[zone]) {
     584      int start = tree.NBAND[zone];
     585      tree.NBAND[zone] = band + 10;
     586      REALLOCATE (tree.ra[zone],   double, tree.NBAND[zone]);
     587      REALLOCATE (tree.dec[zone],  double, tree.NBAND[zone]);
     588      REALLOCATE (tree.Xo[zone],   double, tree.NBAND[zone]);
     589      REALLOCATE (tree.Yo[zone],   double, tree.NBAND[zone]);
     590      REALLOCATE (tree.dX[zone],    float, tree.NBAND[zone]);
     591      REALLOCATE (tree.dY[zone],    float, tree.NBAND[zone]);
     592      REALLOCATE (tree.cell[zone],    int, tree.NBAND[zone]);
     593      REALLOCATE (tree.name[zone], char *, tree.NBAND[zone]);
     594      for (j = start; j < tree.NBAND[zone]; j++) {
     595        tree.ra[zone][j] = NAN;
     596        tree.dec[zone][j] = NAN;
     597        tree.Xo[zone][band] = NAN;
     598        tree.Yo[zone][band] = NAN;
     599        tree.dX[zone][band] = NAN;
     600        tree.dY[zone][band] = NAN;
     601        tree.cell[zone][j] = -1;
     602        ALLOCATE (tree.name[zone][j], char, BOUNDARY_TREE_NAME_LENGTH);
     603      }
     604    }
     605    tree.ra[zone][band] = ra;
     606    tree.dec[zone][band] = dec;
     607    tree.Xo[zone][band] = x;
     608    tree.Yo[zone][band] = y;
     609    tree.dX[zone][band] = image[i].NX / NX_SUB;
     610    tree.dY[zone][band] = image[i].NY / NY_SUB;
     611   
     612    tree.cell[zone][band] = i;
     613
     614    // what are the min and max DEC values for this zone? (test center and corners of top and bottom edge) 
     615
     616    if (dec > 0.0) {
     617      // min DEC at bottom of the cell at the center
     618      x = 0.5*image[i].NX;
     619      y = 0.0*image[i].NY;
     620      XY_to_RD (&ra, &dec, x, y, &image[i].coords);
     621      tree.DEC_min_raw[zone] = MIN(tree.DEC_min_raw[zone], dec);
     622 
     623      // max DEC : find the intersection between the RA boundary and the top of the cell
     624      double ra_band_min = tree.RA_origin[zone] + tree.RA_offset[zone] * band;
     625 
     626      // does the parity matter?
     627      x = 0.0*image[i].NX;
     628      y = 1.0*image[i].NY;
     629 
     630      int Niter;
     631      for (Niter = 0; Niter < 3; Niter ++) {
     632        XY_to_RD (&ra, &dec, x, y, &image[i].coords);
     633        RD_to_XY (&x, &y, ra_band_min, dec, &image[i].coords);
     634        y = 1.0*image[i].NY;
     635      }
     636      tree.DEC_max_raw[zone] = MAX(tree.DEC_max_raw[zone], dec);
     637    } else {
     638      // max DEC at top of the cell at the center
     639      x = 0.5*image[i].NX;
     640      y = 1.0*image[i].NY;
     641      XY_to_RD (&ra, &dec, x, y, &image[i].coords);
     642      tree.DEC_max_raw[zone] = MAX(tree.DEC_max_raw[zone], dec);
     643 
     644      // max DEC : find the intersection between the RA boundary and the bottom of the cell
     645      double ra_band_min = tree.RA_origin[zone] + tree.RA_offset[zone] * band;
     646 
     647      // does the parity matter? (NO)
     648      x = 0.0*image[i].NX;
     649      y = 0.0*image[i].NY;
     650 
     651      int Niter;
     652      for (Niter = 0; Niter < 3; Niter ++) {
     653        XY_to_RD (&ra, &dec, x, y, &image[i].coords);
     654        RD_to_XY (&x, &y, ra_band_min, dec, &image[i].coords);
     655        y = 0.0*image[i].NY;
     656      }
     657      tree.DEC_min_raw[zone] = MIN(tree.DEC_min_raw[zone], dec);
     658    }
     659    memcpy (tree.name[zone][band], image[i].name, BOUNDARY_TREE_NAME_LENGTH);
     660  }
     661
     662  // figure out the max band value for each zone?
     663  for (zone = 0; zone < tree.Nzone; zone++) {
     664    int found_last = FALSE;
     665    int last_band = -1;
     666    for (band = 0; band < tree.NBAND[zone]; band++) {
     667      // all cells should be filled
     668      if (tree.cell[zone][band] < 0) {
     669        if (!found_last) {
     670          found_last = TRUE;
     671          last_band = band;
     672        }
     673      } else {
     674        if (found_last) {
     675          fprintf (stderr, "error: empty cell (%d,%d) after last band (%d)\n", zone, band, last_band);
     676        }
     677      }
     678    }
     679    if (last_band == -1) {
     680      last_band = tree.NBAND[zone];
     681    }
     682    tree.Nband[zone] = last_band;
     683    // fprintf (stderr, "last_band : %d, Nband: %d\n", last_band, tree.Nband[zone]);
     684  }
     685
     686  // figure out the max band value for each zone?
     687  int Nm = 0;
     688  tree.DEC_min[Nm] = -90.0;
     689  tree.DEC_max[Nm] = 0.5*(tree.DEC_max_raw[Nm] + tree.DEC_min_raw[Nm + 1]);
     690
     691  int Np = tree.Nzone - 1;
     692  tree.DEC_min[Np] = 0.5*(tree.DEC_min_raw[Np] + tree.DEC_max_raw[Np - 1]);
     693  tree.DEC_max[Np] = +90.0;
     694
     695  for (zone = 1; zone < tree.Nzone - 1; zone++) {
     696    tree.DEC_min[zone] = 0.5*(tree.DEC_min_raw[zone] + tree.DEC_max_raw[zone - 1]);
     697    tree.DEC_max[zone] = 0.5*(tree.DEC_max_raw[zone] + tree.DEC_min_raw[zone + 1]);
     698  }
     699
     700  INITTIME;
     701
     702  int Npts = 10000000;
     703
     704  // test : find skycell for NN random points on the sky
     705  long A = time(NULL);
     706  long B = A + 10000;
     707  srand48(B);
     708  for (i = 0; i < Npts; i++) {
     709    ra  = 360.0 * drand48();
     710    dec = 180.0 * drand48() - 90.0;
     711    if (!BoundaryTreeCellCoords (&tree, &zone, &band, ra, dec)) {
     712      fprintf (stderr, "failure for %f,%f\n", ra, dec);
     713    }
     714  }
     715  MARKTIME("-- test %d pts: %f sec\n", Npts, dtime);
     716
    416717  if (APPEND) {
    417718    int Ntess = 0;
Note: See TracChangeset for help on using the changeset viewer.