Changeset 36484 for trunk/Ohana/src/addstar/src/findskycell.c

Timestamp:

Feb 6, 2014, 12:22:33 PM (12 years ago)

Author:

eugene

Message:

updates to findskycell for general tess table

File:

• trunk/Ohana/src/addstar/src/findskycell.c (modified) (10 diffs)

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

@@ -12 +12 @@
 // in an even more specific case, RA[i,zone] = RA_origin[zone] + RA_offset[zone]
 
-enum {TREE_NONE, TREE_MAKE, TREE_USE};
+enum {TREE_NONE, TREE_MAKE, TREE_LOCAL, TREE_USE};
+enum {REGION_NONE, REGION_USER, REGION_MIN, REGION_MAX};
 
 void usage (void) {
   fprintf (stderr, "USAGE: findcell -mktree (tree) (catdir) [-nx Nx] [-ny Ny]\n");
+  fprintf (stderr, "USAGE: findcell -mklocal (file) (catdir) [-nx Nx] [-ny Ny]\n");
   fprintf (stderr, "USAGE: findcell -tree (tree) (datafile)\n");
   fprintf (stderr, "   (datafile) should contain a list of RA,DEC pairs\n");
@@ -22 +24 @@
 
 int mktree (char *treefile, char *catdir);
+int mklocal (char *treefile, char *catdir);
 int apply_tree (char *treefile, char *datafile);
 
@@ -28 +31 @@
 int NY_SUB = 1;
 
+double R_MIN = NAN;
+double R_MAX = NAN;
+double D_MIN = NAN;
+double D_MAX = NAN;
+int REGION_OPTION = REGION_NONE;
+
+char *BASENAME = NULL;
+int projectIDoff = -1;
+int skycellIDoff = -1;
+
+int APPEND = FALSE;
+
 int main (int argc, char **argv) {
 
@@ -42 +57 @@
   if (get_argument (argc, argv, "-h")) usage ();
 
-  NX_SUB = NY_SUB = 1;
   if ((N = get_argument (argc, argv, "-nx"))) {
     remove_argument (N, &argc, argv);
@@ -55 +69 @@
 
   /* pixel scale (arcsec/pixel) */
-  SCALE = 1.0;
   if ((N = get_argument (argc, argv, "-scale"))) {
     remove_argument (N, &argc, argv);
     SCALE = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  /* pixel scale (arcsec/pixel) */
+  if ((N = get_argument (argc, argv, "-append"))) {
+    remove_argument (N, &argc, argv);
+    APPEND = TRUE;
+  }
+
+  // user-specified region
+  if ((N = get_argument (argc, argv, "-region"))) {
+    remove_argument (N, &argc, argv);
+    if (N > argc - 4) {
+      fprintf (stderr, "USAGE: -region requires 4 arguments (Rmin Rmax Dmin Dmax)\n");
+      exit (1);
+    }
+    R_MIN = atof (argv[N]); remove_argument (N, &argc, argv);
+    R_MAX = atof (argv[N]); remove_argument (N, &argc, argv);
+    D_MIN = atof (argv[N]); remove_argument (N, &argc, argv);
+    D_MAX = atof (argv[N]); remove_argument (N, &argc, argv);
+    REGION_OPTION = REGION_USER;
+  }
+
+  if ((N = get_argument (argc, argv, "-region-min"))) {
+    remove_argument (N, &argc, argv);
+    REGION_OPTION = REGION_MIN;
+  }
+
+  if ((N = get_argument (argc, argv, "-region-max"))) {
+    remove_argument (N, &argc, argv);
+    REGION_OPTION = REGION_MAX;
+  }
+
+  if ((N = get_argument (argc, argv, "-basename"))) {
+    remove_argument (N, &argc, argv);
+    if (N > argc - 3) {
+      fprintf (stderr, "USAGE: -basename (name) (proj_offset) (skycell_offset)\n");
+      exit (1);
+    }
+    BASENAME = strcreate  (argv[N]);
+    remove_argument (N, &argc, argv);
+    projectIDoff = atoi(argv[N]);
+    remove_argument (N, &argc, argv);
+    skycellIDoff = atoi(argv[N]);
     remove_argument (N, &argc, argv);
   }
@@ -70 +127 @@
     remove_argument (N, &argc, argv);
   }
+  if ((N = get_argument (argc, argv, "-mklocal"))) {
+    MODE = TREE_LOCAL;
+    remove_argument (N, &argc, argv);
+    treefile = strcreate (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
   if ((N = get_argument (argc, argv, "-tree"))) {
     MODE = TREE_USE;
@@ -83 +146 @@
   if (MODE == TREE_MAKE) {
     mktree (treefile, argv[1]);
+    exit (0);
+  }
+
+  if (MODE == TREE_LOCAL) {
+    mklocal (treefile, argv[1]);
     exit (0);
   }
@@ -348 +416 @@
   MARKTIME("-- test %d pts: %f sec\n", Npts, dtime);
 
-  BoundaryTreeSave (treefile, &tree);
+  if (APPEND) {
+    int Ntess = 0;
+    TessellationTable *tess = TessellationTableLoad (treefile, &Ntess);
+    REALLOCATE (tess, TessellationTable, Ntess + 1);
+    TessellationTableInit (&tess[Ntess], 1);
+    tess[Ntess].tree = &tree;
+    tess[Ntess].type = TESS_RINGS;
+    tess[Ntess].Rmin =   0;
+    tess[Ntess].Rmax = 360;
+    tess[Ntess].Dmin = -90;
+    tess[Ntess].Dmax = +90;
+
+    if (BASENAME) {
+      tess[Ntess].Nbasename = strlen(BASENAME);
+      tess[Ntess].basename = strcreate(BASENAME);
+      tess[Ntess].projectIDoff = projectIDoff;
+      tess[Ntess].skycellIDoff = skycellIDoff;
+    }
+      
+    // add basename an related here...
+    Ntess ++;
+    TessellationTableSave (treefile, tess, Ntess);
+  } else {
+    BoundaryTreeSave (treefile, &tree);
+  }
 
   return TRUE;
@@ -354 +446 @@
 
 int apply_tree (char *treefile, char *datafile) {
+
+  int Ntess = 0;
+  TessellationTable *tess = TessellationTableLoad (treefile, &Ntess);
+  if (!tess) {
+    fprintf (stderr, "error loading tessellation table file %s\n", treefile);
+    exit (2);
+  }
+
+  FILE *f = fopen (datafile, "r");
+  if (!f) {
+    fprintf (stderr, "error opening data file %s\n", datafile);
+    exit (3);
+  }
+
+  double ra, dec;
+  int Nvalue = 0;
+  while ((Nvalue = fscanf (f, "%lf %lf", &ra, &dec)) != EOF) {
+
+    int tessID, projID, skycellID;
+
+    if (!TessellationPrimaryCellIDs (tess, Ntess, &tessID, &projID, &skycellID, ra, dec)) {
+      fprintf (stderr, "error finding cell for %f,%f\n", ra, dec);
+      continue;
+    }
+
+    fprintf (stdout, "%10.6f %10.6f : %2d  %04d %03d : %s\n", ra, dec, tessID, projID, skycellID, tess[tessID].basename);
+  }
+
+  exit (0);
+}
+
+int apply_tree_old (char *treefile, char *datafile) {
 
   BoundaryTree *tree = BoundaryTreeLoad (treefile);
@@ -405 +529 @@
 }
 
+// a given tess is defined by a catdir (more than one per catdir?)
+// this function takes a catdir and generates an extension for a tess file for a local tess.
+
+/*
+
+  LOCAL projections are defined by single projection cells.  should I be supplying the info on
+  the cmd line or figure out the bounds from the catdir?  
+
+ */
+
+int mklocal (char *treefile, char *catdir) {
+
+  int i, j, status;
+  FITS_DB db;
+  Image *image;
+  off_t Nimage;
+  double x, y, ra, dec;
+
+  if (REGION_OPTION == REGION_NONE) {
+    fprintf (stderr, "ERROR: need to define bounding region (-region Rmin Rmax Dmin Dmax | -region-min | -region-max)\n");
+    exit (2);
+  }
+
+  if (!BASENAME) {
+    fprintf (stderr, "ERROR: need to define BASENAME -basename (name) (proj_offset) (skycell_offset)\n");
+    exit (2);
+  }
+
+  char imagefile[DVO_MAX_PATH];
+  snprintf (imagefile, DVO_MAX_PATH, "%s/Images.dat", catdir);
+
+  status = dvo_image_lock (&db, imagefile, 2.0, LCK_XCLD);
+  if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db.filename);
+
+  /* load or create the image table */
+  if (db.dbstate == LCK_EMPTY) Shutdown ("can't read image catalog %s", db.filename);
+
+  if (!dvo_image_load (&db, TRUE, FALSE)) Shutdown ("can't read image catalog %s", db.filename);
+
+  // convert database table to internal structure (binary to Image)
+  // 'image' points to the same memory as db->ftable->buffer
+  image = gfits_table_get_Image (&db.ftable, &Nimage, &db.swapped);
+  if (!image) {
+    fprintf (stderr, "ERROR: failed to read images\n");
+    exit (2);
+  }
+  
+  // generate an empty BoundaryTree
+  TessellationTable *tess = NULL;
+
+  int NtessDisk = 0;
+  if (APPEND) {
+    tess = TessellationTableLoad (treefile, &NtessDisk);
+    REALLOCATE (tess, TessellationTable, Nimage + NtessDisk);
+  } else {
+    ALLOCATE (tess, TessellationTable, Nimage + NtessDisk);
+  }
+
+  TessellationTableInit (&tess[NtessDisk], Nimage);
+
+  // find the RA,DEC of the image centers & assign to cells
+  int Ntess = NtessDisk;
+  for (i = 0; i < Nimage; i++) {
+    // user supplied values, do not try to derive from Image.dat
+    tess[Ntess].NX_SUB = NX_SUB;
+    tess[Ntess].NY_SUB = NY_SUB;
+    tess[Ntess].dPix = SCALE/3600.0;
+
+    x = 0.5*image[i].NX;
+    y = 0.5*image[i].NY;
+    XY_to_RD (&ra, &dec, x, y, &image[i].coords);
+
+    tess[Ntess].Ro  = ra;
+    tess[Ntess].Do  = dec;
+    tess[Ntess].Xo  = x;
+    tess[Ntess].Yo  = y;
+    tess[Ntess].dX  = image[i].NX / NX_SUB;
+    tess[Ntess].dY  = image[i].NY / NY_SUB;
+
+    // find the minimum or maximum containing region
+    if ((REGION_OPTION == REGION_MIN) || (REGION_OPTION == REGION_MAX)) {
+      // XXX short cut for now : if the projection cell bounds the equator or 0,360 boundary, this will fail:
+      double R[4], D[4];
+      XY_to_RD (&R[0], &D[0],           0,           0, &image[i].coords);
+      XY_to_RD (&R[1], &D[1], image[i].NX,           0, &image[i].coords);
+      XY_to_RD (&R[2], &D[2],           0, image[i].NY, &image[i].coords);
+      XY_to_RD (&R[3], &D[3], image[i].NX, image[i].NY, &image[i].coords);
+
+      if (REGION_OPTION == REGION_MIN) { 
+        for (j = 0; j < 4; j++) {
+          R_MIN = (R[j] < ra)  ? (isfinite(R_MIN) ? MAX(R_MIN, R[j]) : R[j]) : R_MIN;
+          R_MAX = (R[j] > ra)  ? (isfinite(R_MAX) ? MIN(R_MAX, R[j]) : R[j]) : R_MAX;
+          D_MIN = (D[j] < dec) ? (isfinite(D_MIN) ? MAX(D_MIN, D[j]) : D[j]) : D_MIN;
+          D_MAX = (D[j] > dec) ? (isfinite(D_MAX) ? MIN(D_MAX, D[j]) : D[j]) : D_MAX;
+        }
+      } else {
+        for (j = 0; j < 4; j++) {
+          R_MIN = (R[j] < ra)  ? (isfinite(R_MIN) ? MIN(R_MIN, R[j]) : R[j]) : R_MIN;
+          R_MAX = (R[j] > ra)  ? (isfinite(R_MAX) ? MAX(R_MAX, R[j]) : R[j]) : R_MAX;
+          D_MIN = (D[j] < dec) ? (isfinite(D_MIN) ? MIN(D_MIN, D[j]) : D[j]) : D_MIN;
+          D_MAX = (D[j] > dec) ? (isfinite(D_MAX) ? MAX(D_MAX, D[j]) : D[j]) : D_MAX;
+        }
+      }
+    }
+    tess[Ntess].Rmin = R_MIN;
+    tess[Ntess].Rmax = R_MAX;
+    tess[Ntess].Dmin = D_MIN;
+    tess[Ntess].Dmax = D_MAX;
+
+    tess[Ntess].type = TESS_LOCAL;
+
+    // XXX I don't really want to do the work of discovering the rule...
+    tess[Ntess].Nbasename = strlen(BASENAME);
+    tess[Ntess].basename = strcreate(BASENAME);
+
+    tess[Ntess].projectIDoff = projectIDoff;
+    tess[Ntess].skycellIDoff = skycellIDoff;
+
+    Ntess ++;
+  }
+
+  TessellationTableSave (treefile, tess, Ntess);
+
+  return TRUE;
+}
+