IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 34827


Ignore:
Timestamp:
Dec 16, 2012, 6:28:32 AM (14 years ago)
Author:
eugene
Message:

boundary tree now basically works for skycells (not just projection cells)

Location:
branches/eam_branches/ipp-20121130/Ohana/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r34826 r34827  
    2424int apply_tree (char *treefile, char *datafile);
    2525
     26float SCALE = 1.0;
     27int NX_SUB = 1;
     28int NY_SUB = 1;
     29
    2630int main (int argc, char **argv) {
    2731
    28   int N, NX_SUB, NY_SUB;
     32  int N;
    2933  char *treefile = NULL;
    3034
    … …  
    4751    remove_argument (N, &argc, argv);
    4852    NY_SUB = atof (argv[N]);
     53    remove_argument (N, &argc, argv);
     54  }
     55
     56  /* pixel scale (arcsec/pixel) */
     57  SCALE = 1.0;
     58  if ((N = get_argument (argc, argv, "-scale"))) {
     59    remove_argument (N, &argc, argv);
     60    SCALE = atof (argv[N]);
    4961    remove_argument (N, &argc, argv);
    5062  }
    … …  
    121133
    122134  tree.Nzone = 46;
     135  tree.NX_SUB = NX_SUB;
     136  tree.NY_SUB = NY_SUB;
     137  tree.dPix = SCALE/3600.0;
    123138
    124139  ALLOCATE (tree.Nband, int, tree.Nzone);
    … …  
    134149  ALLOCATE (tree.ra,   double *, tree.Nzone);
    135150  ALLOCATE (tree.dec,  double *, tree.Nzone);
     151  ALLOCATE (tree.Xo,      int *, tree.Nzone);
     152  ALLOCATE (tree.Yo,      int *, tree.Nzone);
     153  ALLOCATE (tree.dX,      int *, tree.Nzone);
     154  ALLOCATE (tree.dY,      int *, tree.Nzone);
    136155  ALLOCATE (tree.cell,    int *, tree.Nzone);
    137156  ALLOCATE (tree.name,  char **, tree.Nzone);
    … …  
    164183    ALLOCATE (tree.ra[zone],   double, tree.NBAND[zone]);
    165184    ALLOCATE (tree.dec[zone],  double, tree.NBAND[zone]);
     185    ALLOCATE (tree.Xo[zone],      int, tree.NBAND[zone]);
     186    ALLOCATE (tree.Yo[zone],      int, tree.NBAND[zone]);
     187    ALLOCATE (tree.dX[zone],      int, tree.NBAND[zone]);
     188    ALLOCATE (tree.dY[zone],      int, tree.NBAND[zone]);
    166189    ALLOCATE (tree.cell[zone], int,    tree.NBAND[zone]);
    167190    ALLOCATE (tree.name[zone], char *, tree.NBAND[zone]);
    … …  
    169192      tree.ra[zone][band] = NAN;
    170193      tree.dec[zone][band] = NAN;
     194      tree.Xo[zone][band] = -1;
     195      tree.Yo[zone][band] = -1;
     196      tree.dX[zone][band] = -1;
     197      tree.dY[zone][band] = -1;
    171198      tree.cell[zone][band] = -1;
    172199      ALLOCATE (tree.name[zone][band], char, BOUNDARY_TREE_NAME_LENGTH);
    … …  
    191218      REALLOCATE (tree.ra[zone],   double, tree.NBAND[zone]);
    192219      REALLOCATE (tree.dec[zone],  double, tree.NBAND[zone]);
    193       REALLOCATE (tree.cell[zone], int,    tree.NBAND[zone]);
     220      REALLOCATE (tree.Xo[zone],      int, tree.NBAND[zone]);
     221      REALLOCATE (tree.Yo[zone],      int, tree.NBAND[zone]);
     222      REALLOCATE (tree.dX[zone],      int, tree.NBAND[zone]);
     223      REALLOCATE (tree.dY[zone],      int, tree.NBAND[zone]);
     224      REALLOCATE (tree.cell[zone],    int, tree.NBAND[zone]);
    194225      REALLOCATE (tree.name[zone], char *, tree.NBAND[zone]);
    195226      for (j = start; j < tree.NBAND[zone]; j++) {
    196227        tree.ra[zone][j] = NAN;
    197228        tree.dec[zone][j] = NAN;
     229        tree.Xo[zone][band] = -1;
     230        tree.Yo[zone][band] = -1;
     231        tree.dX[zone][band] = -1;
     232        tree.dY[zone][band] = -1;
    198233        tree.cell[zone][j] = -1;
    199234        ALLOCATE (tree.name[zone][j], char, BOUNDARY_TREE_NAME_LENGTH);
    … …  
    204239    tree.Xo[zone][band] = x;
    205240    tree.Yo[zone][band] = y;
     241    tree.dX[zone][band] = image[i].NX / NX_SUB;
     242    tree.dY[zone][band] = image[i].NY / NY_SUB;
     243   
    206244    tree.cell[zone][band] = i;
    207245
    … …  
    349387    // convert R,D to X,Y with hard-wired projection and scale, orientation?
    350388   
    351     TreeCellProjection (&x, &y, ra, dec, tree, zone, band);
    352 
    353     xi = x / tree->NX_SUB;
    354     yi = y / tree->NY_SUB;
    355 
     389    double x = 0.0;
     390    double y = 0.0;
     391    BoundaryTreeProjection (&x, &y, ra, dec, tree, zone, band);
     392    fprintf (stdout, "%6.1f %6.1f\n", x, y);
     393
     394    int xi = x / tree->dX[zone][band];
     395    int yi = y / tree->dY[zone][band];
    356396   
    357 
     397    char format[24], skycellname[128];
     398    int Ndigit = (int)(log10(tree->NX_SUB*tree->NY_SUB)) + 1 ;
     399    snprintf (format, 24, "%s.%%0%dd", tree->name[zone][band], Ndigit);
     400
     401    int N = xi + tree->NX_SUB * yi;
     402    snprintf (skycellname, 128, format, N);
     403
     404    fprintf (stderr, "skycell: %s\n", skycellname);
    358405  }
    359406

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

    r34826 r34827  
    327327  char  ***name; // projection cell name
    328328 
    329   int NX_SUB, NY_SUB;
     329  int NX_SUB;
     330  int NY_SUB;
     331  double dPix;
    330332
    331333  int **Xo;
    332334  int **Yo;
     335  int **dX;
     336  int **dY;
    333337} BoundaryTree;
    334338
    … …  
    706710int BoundaryTreeSave(char *filename, BoundaryTree *tree);
    707711BoundaryTree *BoundaryTreeLoad(char *filename);
     712int BoundaryTreeProjection (double *x, double *y, double r, double d, BoundaryTree *tree, int zone, int band);
    708713
    709714void dvo_average_init (Average *average);

  • branches/eam_branches/ipp-20121130/Ohana/src/libdvo/src/BoundaryTree.c

    r34826 r34827  
    4747  gfits_scan (&header, "DEC_ORI", "%lf", 1, &tree->DEC_origin);
    4848  gfits_scan (&header, "DEC_OFF", "%lf", 1, &tree->DEC_offset);
     49
     50  gfits_scan (&header, "NX_SUB", "%d", 1, &tree->NX_SUB);
     51  gfits_scan (&header, "NY_SUB", "%d", 1, &tree->NY_SUB);
     52  gfits_scan (&header, "PIXSCALE", "%lf", 1, &tree->dPix);
    4953
    5054  ftable.header = &theader;
    … …  
    7579    ALLOCATE (tree->ra,   double *, tree->Nzone);
    7680    ALLOCATE (tree->dec,  double *, tree->Nzone);
    77     ALLOCATE (tree->cell, int *, tree->Nzone);
    78     ALLOCATE (tree->name, char **, tree->Nzone);
     81    ALLOCATE (tree->Xo,      int *, tree->Nzone);
     82    ALLOCATE (tree->Yo,      int *, tree->Nzone);
     83    ALLOCATE (tree->dX,      int *, tree->Nzone);
     84    ALLOCATE (tree->dY,      int *, tree->Nzone);
     85    ALLOCATE (tree->cell,    int *, tree->Nzone);
     86    ALLOCATE (tree->name,  char **, tree->Nzone);
    7987    for (i = 0; i < tree->Nzone; i++) {
    8088      ALLOCATE (tree->ra[i],   double, tree->Nband[i]);
    8189      ALLOCATE (tree->dec[i],  double, tree->Nband[i]);
    82       ALLOCATE (tree->cell[i], int,    tree->Nband[i]);
     90      ALLOCATE (tree->Xo[i],      int, tree->Nband[i]);
     91      ALLOCATE (tree->Yo[i],      int, tree->Nband[i]);
     92      ALLOCATE (tree->dX[i],      int, tree->Nband[i]);
     93      ALLOCATE (tree->dY[i],      int, tree->Nband[i]);
     94      ALLOCATE (tree->cell[i],    int, tree->Nband[i]);
    8395      ALLOCATE (tree->name[i], char *, tree->Nband[i]);
    8496      for (j = 0; j < tree->Nband[i]; j++) {
    … …  
    102114    GET_COLUMN_NEW(band,  "BAND",        int);
    103115    GET_COLUMN_NEW(index, "INDEX",       int);
     116    GET_COLUMN_NEW(Xo,    "X_CENT",      int);
     117    GET_COLUMN_NEW(Yo,    "Y_CENT",      int);
     118    GET_COLUMN_NEW(dX,    "X_GRID",      int);
     119    GET_COLUMN_NEW(dY,    "Y_GRID",      int);
    104120    GET_COLUMN_NEW(name,  "NAME",        char); // XXX how is this done?
    105121    gfits_free_header (&theader);
    … …  
    113129      tree->ra[nz][nb] = R[i];
    114130      tree->dec[nz][nb] = D[i];
     131      tree->Xo[nz][nb] = Xo[i];
     132      tree->Yo[nz][nb] = Yo[i];
     133      tree->dX[nz][nb] = dX[i];
     134      tree->dY[nz][nb] = dY[i];
    115135      tree->cell[nz][nb] = i; // XXX ?
    116136      memcpy(tree->name[nz][nb], &name[i*BOUNDARY_TREE_NAME_LENGTH], BOUNDARY_TREE_NAME_LENGTH);
    … …  
    121141    free (zone  );
    122142    free (band  );
     143    free (Xo    );
     144    free (Yo    );
     145    free (dX    );
     146    free (dY    );
    123147    free (index );
    124148    free (name  );
    … …  
    165189  gfits_modify (&header, "DEC_ORI", "%lf", 1, tree->DEC_origin);
    166190  gfits_modify (&header, "DEC_OFF", "%lf", 1, tree->DEC_offset);
     191
     192  gfits_modify (&header, "NX_SUB", "%d", 1, tree->NX_SUB);
     193  gfits_modify (&header, "NY_SUB", "%d", 1, tree->NY_SUB);
     194  gfits_modify (&header, "PIXSCALE", "%lf", 1, tree->dPix);
    167195
    168196  gfits_fwrite_header  (f, &header);
    … …  
    223251    gfits_define_bintable_column (&theader, "J", "BAND", "band sequence number", "none", 1.0, 0.0);
    224252    gfits_define_bintable_column (&theader, "J", "INDEX","cell index", "none", 1.0, 0.0);
     253    gfits_define_bintable_column (&theader, "J", "X_CENT", "projection cell center pixel", "none", 1.0, 0.0);
     254    gfits_define_bintable_column (&theader, "J", "Y_CENT", "projection cell center pixel", "none", 1.0, 0.0);
     255    gfits_define_bintable_column (&theader, "J", "X_GRID", "skycell grid spacing", "none", 1.0, 0.0);
     256    gfits_define_bintable_column (&theader, "J", "Y_GRID", "skycell grid spacing", "none", 1.0, 0.0);
    225257    gfits_define_bintable_column (&theader, fmt, "NAME", "cell name", "none", 1.0, 0.0);
    226258
    … …  
    240272    int    *band          ; ALLOCATE (band,  int,    Ncell);
    241273    int    *index         ; ALLOCATE (index, int,    Ncell);
     274    int    *Xo            ; ALLOCATE (Xo,    int,    Ncell);
     275    int    *Yo            ; ALLOCATE (Yo,    int,    Ncell);
     276    int    *dX            ; ALLOCATE (dX,    int,    Ncell);
     277    int    *dY            ; ALLOCATE (dY,    int,    Ncell);
    242278    char   *name          ; ALLOCATE (name,  char,   Ncell*BOUNDARY_TREE_NAME_LENGTH);
    243279
    … …  
    251287        R[i]     = tree->ra[nz][nb];
    252288        D[i]     = tree->dec[nz][nb];
     289        Xo[i]    = tree->Xo[nz][nb];
     290        Yo[i]    = tree->Yo[nz][nb];
     291        dX[i]    = tree->dX[nz][nb];
     292        dY[i]    = tree->dY[nz][nb];
    253293        zone[i]  = nz;
    254294        band[i]  = nb;
    … …  
    265305    gfits_set_bintable_column (&theader, &ftable, "BAND",  band,  Ncell);
    266306    gfits_set_bintable_column (&theader, &ftable, "INDEX", index, Ncell);
     307    gfits_set_bintable_column (&theader, &ftable, "X_CENT", Xo,   Ncell);
     308    gfits_set_bintable_column (&theader, &ftable, "Y_CENT", Yo,   Ncell);
     309    gfits_set_bintable_column (&theader, &ftable, "X_GRID", dX,   Ncell);
     310    gfits_set_bintable_column (&theader, &ftable, "Y_GRID", dY,   Ncell);
    267311    gfits_set_bintable_column (&theader, &ftable, "NAME",  name,  Ncell);
    268312
    … …  
    272316    free (band  );
    273317    free (index );
     318    free (Xo    );
     319    free (Yo    );
     320    free (dX    );
     321    free (dY    );
    274322    free (name  );
    275323
    … …  
    322370// need Ro, Do, Xo, Yo, dPix
    323371
    324 int TreeCellProjection (double *x, double *y, double r, double d, BoundaryTree *tree, int zone, int band) {
     372int BoundaryTreeProjection (double *x, double *y, double r, double d, BoundaryTree *tree, int zone, int band) {
    325373
    326374    double Xo = tree->Yo[zone][band];
    … …  
    332380    // this block only depends on Ro, Do
    333381
    334     sdp  = sin(RAD_DEG*Do);
    335     cdp  = cos(RAD_DEG*Do);
    336     salp = sin(RAD_DEG*(ra - Ro));
    337     calp = cos(RAD_DEG*(ra - Ro));
    338     sdel = sin(RAD_DEG*dec);
    339     cdel = cos(RAD_DEG*dec);
     382    double sdp  = sin(RAD_DEG*Do);
     383    double cdp  = cos(RAD_DEG*Do);
     384    double salp = sin(RAD_DEG*(r - Ro));
     385    double calp = cos(RAD_DEG*(r - Ro));
     386    double sdel = sin(RAD_DEG*d);
     387    double cdel = cos(RAD_DEG*d);
    340388   
    341     stht = sdel*sdp + cdel*cdp*calp;    /* sin(theta) */
    342     sphi = cdel*salp;                   /* = cos(theta)*sin(phi) */
    343     cphi = cdel*sdp*calp - sdel*cdp;    /* = cos(theta)*cos(phi) */
     389    double stht = sdel*sdp + cdel*cdp*calp;    /* sin(theta) */
     390    double sphi = cdel*salp;                   /* = cos(theta)*sin(phi) */
     391    double cphi = cdel*sdp*calp - sdel*cdp;    /* = cos(theta)*cos(phi) */
    344392
    345393    // defines the TAN projection (one of zenithal projections available, libdvo/src/coordops.c
    346394    // R = cot (theta) = cos(theta) / sin(theta)
     395    double L, M;
    347396    if (stht == 0) {
    348         Rc = hypot(sphi, cphi);
     397        double Rc = hypot(sphi, cphi);
    349398        L = 180.0 * sphi / Rc;
    350399        M = 180.0 * cphi / Rc;
    … …  
    362411    // Yo = (coords[0].pc1_1*M - coords[0].pc2_1*L) / Ro;
    363412
    364     Xo = L;
    365     Yo = M;
     413    double X = L;
     414    double Y = M;
    366415
    367416    // scale is dPix
    368417
    369     *x = Xo / dPix + Xo;
    370     *y = Yo / dPix + Yo;
     418    *x = Xo - X / dPix;
     419    *y = Yo + Y / dPix;
    371420   
    372421    return TRUE;
Note: See TracChangeset for help on using the changeset viewer.