IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 15997


Ignore:
Timestamp:
Jan 4, 2008, 11:15:56 AM (18 years ago)
Author:
eugene
Message:

added option to align y-axis with north; added euler rotations; added padding

Location:
trunk/Ohana/src/addstar
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/addstar/include/skycells.h

    r13186 r15997  
    2424  double rv[3], dv[3];        // triangle center (2d)
    2525} SkyTriangle;
     26
     27typedef struct {
     28  Coords coords;
     29  int NX;
     30  int NY;
     31  int code;
     32  char name[64];
     33} SkyRectangle;
    2634
    2735/* globals which define database info / data sources (KEEP) */
    … …  
    5260int    VERBOSE;
    5361int    MODE;
     62int    FIX_NS;
    5463int    NMAX;
    5564int    NX_SUB, NY_SUB;
    5665double SCALE;
     66double PADDING;
    5767
    58 void       SetProtect             PROTO((int mode));
    59 int        SetSignals             PROTO(());
    60 int        Shutdown               PROTO((char *message, ...); )
    61 void       TrapSignal             PROTO((int sig));
    62 int args_skycells (int argc, char **argv);
    63 int ConfigInit_skycells (int *argc, char **argv);
     68double EULER_A;
     69double EULER_B;
    6470
    65 int sky_triangle_to_image (Image *image, SkyTriangle *triangle);
    66 int sky_triangle_coords (SkyTriangle *triangle);
    67 SkyTriangle *sky_divide_triangles (SkyTriangle *in, int *ntriangles);
    68 Point sky_divide_edge (Point v1, Point v2);
    69 SkyTriangle *sky_base_triangles (int *ntriangles);
    70 int sky_tessalation_init (double scale);
    71 int sky_triangle_to_rectangle (Image *image, SkyTriangle *triangle);
     71void         SetProtect                PROTO((int mode));
     72int          SetSignals                PROTO(());
     73int          Shutdown                  PROTO((char *message, ...); )
     74void         TrapSignal                PROTO((int sig));
    7275
    73 int sky_tessalation (FITS_DB *db, int level, int Nmax, int mode, double scale);
    74 int sky_tessalation_triangles (FITS_DB *db, int level, int Nmax);
    75 int sky_tessalation_squares (FITS_DB *db, int level, int Nmax);
     76int          args_skycells             PROTO((int argc, char **argv));
     77int          ConfigInit_skycells       PROTO((int *argc, char **argv));
     78int          sky_tessalation           PROTO((FITS_DB *db, int level, int Nmax, int mode, double scale));
     79int          sky_tessalation_init      PROTO((double scale));
     80int          sky_tessalation_triangles PROTO((FITS_DB *db, int level, int Nmax));
     81int          sky_tessalation_squares   PROTO((FITS_DB *db, int level, int Nmax));
     82int          sky_triangle_to_image     PROTO((Image *image, SkyTriangle *triangle));
     83int          sky_triangle_to_rectangle PROTO((SkyRectangle *image, SkyTriangle *triangle));
     84int          sky_subdivide_image       PROTO((Image *output, SkyRectangle *input, int Nx, int Ny));
     85int          sky_triangle_coords       PROTO((SkyTriangle *triangle));
     86SkyTriangle *sky_divide_triangles      PROTO((SkyTriangle *in, int *ntriangles));
     87SkyTriangle *sky_base_triangles        PROTO((int *ntriangles));
     88int          sky_base_rotation         PROTO((SkyTriangle *base, int Nbase));
     89Point        sky_divide_edge           PROTO((Point v1, Point v2));
    7690
    77 // migrate to libdvo eventually
    78 int dvo_image_clear_vtable (FITS_DB *db);
     91// XXX migrate to libdvo eventually
     92int dvo_image_clear_vtable             PROTO((FITS_DB *db));

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

    r13186 r15997  
    44int args_skycells (int argc, char **argv) {
    55 
    6   int i, N;
     6  int N;
    77  char *ptr;
    88
    … …  
    2727  }
    2828
     29  /* what type of output files? */
     30  FIX_NS = FALSE;
     31  if ((N = get_argument (argc, argv, "-fix-ns"))) {
     32    FIX_NS = TRUE;
     33    remove_argument (N, &argc, argv);
     34  }
     35
    2936  /* pixel scale (arcsec/pixel) */
    3037  SCALE = 1.0;
    … …  
    3239    remove_argument (N, &argc, argv);
    3340    SCALE = atof (argv[N]);
     41    remove_argument (N, &argc, argv);
     42  }
     43
     44  /* pixel scale (arcsec/pixel) */
     45  EULER_A = EULER_B = 0.0;
     46  if ((N = get_argument (argc, argv, "-euler"))) {
     47    remove_argument (N, &argc, argv);
     48    EULER_A = RAD_DEG*atof (argv[N]);
     49    remove_argument (N, &argc, argv);
     50    EULER_B = RAD_DEG*atof (argv[N]);
     51    remove_argument (N, &argc, argv);
     52  }
     53
     54  /* pixel scale (arcsec/pixel) */
     55  PADDING = 0.0;
     56  if ((N = get_argument (argc, argv, "-padding"))) {
     57    remove_argument (N, &argc, argv);
     58    PADDING = atof (argv[N]);
    3459    remove_argument (N, &argc, argv);
    3560  }
    … …  
    6186  }
    6287
    63   fprintf (stderr, "USAGE: skycells (level) [-triangles] [-scale arcsec/pix] [-nmax Max Cells in memory] [-nx (Nx cells)] [-ny (Ny cells)]\n");
     88  fprintf (stderr, "USAGE: skycells (level) [-scale arcsec/pix] [-nx (Nx cells)] [-ny (Ny cells)]\n");
     89  fprintf (stderr, "  [-h for details and other options]\n");
    6490  exit (2);
    6591}
    … …  
    7197
    7298  fprintf (stderr, "  -v                          : verbose mode\n");
     99  fprintf (stderr, "  -triangles                  : save base triangles instead of skycells\n");
     100  fprintf (stderr, "  -fix-ns                     : orient skycells with y-axis aligned with Dec\n");
     101  fprintf (stderr, "  -scale                      : set pixel scale (default 1.0 arcsec / pixel)\n");
     102  fprintf (stderr, "  -padding                    : pad skycells by this fraction in each dimension\n");
     103  fprintf (stderr, "  -nmax                       : only keep nmax skycells in memory\n");
     104  fprintf (stderr, "  -nx                         : subdivide skycell projection in x\n");
     105  fprintf (stderr, "  -ny                         : subdivide skycell projection in y\n");
    73106  fprintf (stderr, "  -help                       : this list\n");
    74107  fprintf (stderr, "  -h                          : this list\n\n");

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

    r13886 r15997  
    33# define iSWAP(X,Y) {int tmp=(X); (X) = (Y); (Y) = tmp;}
    44
     5// we use a static refcoords structure to avoid multiple alloc / init steps
    56static Coords *refcoords = NULL;
    6 static int warned = FALSE;
    77
    88int sky_tessalation (FITS_DB *db, int level, int Nmax, int mode, double scale) {
    … …  
    3333  // generate the initial base set
    3434  base = sky_base_triangles (&Nbase);
     35
     36  sky_base_rotation (base, Nbase);
    3537
    3638  // how many triangles total for this level?
    … …  
    9597int sky_tessalation_squares (FITS_DB *db, int level, int Nmax) {
    9698
    97   int i, j, Nname, Ndigit, Ntriangles, Nbase, Ntotal, Ntop, Ltop, Nout, Nx, Ny;
     99  int i, j, Nname, Ndigit, Ntriangles, Nbase, Nimage, Ntotal, Ntop, Ltop, Nsubset, Nx, Ny;
    98100  double fLtop;
    99101  SkyTriangle *base, *tri, *new;
    100   Image *image, *out;
     102  SkyRectangle *rectangle, *subset;
     103  Image *image;
    101104  char format[16];
    102105
    … …  
    106109  // generate the initial base set
    107110  base = sky_base_triangles (&Nbase);
     111
     112  sky_base_rotation (base, Nbase);
    108113
    109114  // how many total cells for this level (multiply by subdivisions, if used)?
    … …  
    144149    }
    145150
    146     // convert the SkyTriangles to Image
    147     ALLOCATE (image, Image, Ntriangles);
     151    // convert the SkyTriangles to SkyRectangles
     152    ALLOCATE (rectangle, SkyRectangle, Ntriangles);
    148153    for (j = 0; j < Ntriangles; j++) {
    149       sky_triangle_to_rectangle (&image[j], &tri[j]);
     154      sky_triangle_to_rectangle (&rectangle[j], &tri[j]);
    150155    } 
    151156
    152157    // drop the appropriate subset
    153     ALLOCATE (out, Image, Ntriangles);
    154     for (j = Nout = 0; j < Ntriangles; j++) {
    155       if (!strcmp(image[j].coords.ctype, "DROP")) continue;
    156       memcpy (&out[Nout], &image[j], sizeof(Image));
    157       snprintf (out[Nout].name, 64, format, Nname);
     158    ALLOCATE (subset, SkyRectangle, Ntriangles);
     159    for (j = Nsubset = 0; j < Ntriangles; j++) {
     160      if (!strcmp(rectangle[j].coords.ctype, "DROP")) continue;
     161      memcpy (&subset[Nsubset], &rectangle[j], sizeof(SkyRectangle));
     162      snprintf (subset[Nsubset].name, 64, format, Nname);
    158163      Nname++;
    159       Nout++;
     164      Nsubset++;
    160165    } 
    161     free (image);
     166    free (rectangle);
    162167
    163168    // subdivide each image (Nx x Ny subcells)
    164     if (Nx*Ny > 1) {
    165       ALLOCATE (image, Image, Nx*Ny*Nout);
    166       for (j = 0; j < Nout; j++) {
    167         sky_subdivide_image (&image[j*Nx*Ny], &out[j], Nx, Ny);
    168       }
    169       free (out);
    170       out = image;
    171       Nout *= Nx*Ny;
     169    Nimage = Nx*Ny*Nsubset;
     170    ALLOCATE (image, Image, Nimage);
     171    for (j = 0; j < Nsubset; j++) {
     172      // convert the SkyRectangles to Images for output
     173      sky_subdivide_image (&image[j*Nx*Ny], &subset[j], Nx, Ny);
    172174    }
    173175
    174176    /* add the new images and save */
    175     dvo_image_addrows (db, out, Nout);
     177    dvo_image_addrows (db, image, Nimage);
    176178    SetProtect (TRUE);
    177179    dvo_image_update (db, VERBOSE);
    … …  
    179181    dvo_image_clear_vtable (db);
    180182
    181     free (out);
     183    free (subset);
     184    free (image);
    182185    free (tri);
    183186  }
    … …  
    188191int sky_triangle_to_image (Image *image, SkyTriangle *triangle) {
    189192
    190   int i, n, parity, peak, b1, b2, NX, NY;
     193  int i, NX, NY;
    191194  double xv[3], yv[3];        // coordinates of the vertex in the reference projection
    192   double xo, yo, xc, yc, angle, scale;
     195  double scale;
    193196  double Xmin, Xmax, Ymin, Ymax;
    194197
    … …  
    250253}
    251254
    252 // an allocated image set is supplied, we fill in the values
    253 int sky_subdivide_image (Image *output, Image *input, int Nx, int Ny) {
    254 
    255   int i, j, N, NX, NY, Ndigit;
    256   char format[24];
    257 
    258   NX = input[0].NX/(double)Nx + 0.5;
    259   NY = input[0].NY/(double)Ny + 0.5;
    260 
    261   Ndigit = (int)(log10(Nx*Ny)) + 1 ;
    262   snprintf (format, 24, "%s.%%0%dd", input[0].name, Ndigit);
    263 
    264   N = 0;
    265   for (j = 0; j < Ny; j++) {
    266     for (i = 0; i < Nx; i++) {
    267 
    268       memcpy (&output[N], &input[0], sizeof(Image));
    269       snprintf (output[N].name, 64, format, N);
    270 
    271       // output[0].coords = input[0].coords;
    272       // strcpy (output[0].coords.ctype, input[0].coords.ctype);
    273 
    274       output[N].NX = NX;
    275       output[N].NY = NY;
    276 
    277       output[N].coords.crpix1 = input[0].coords.crpix1 - i*NX;
    278       output[N].coords.crpix2 = input[0].coords.crpix2 - j*NY;
    279       N++;
    280     }
    281   }
    282   return (TRUE);
    283 }
    284 
    285255// an allocated image is supplied, we fill in the values
    286256// we are only keeping ~half of the images
    287 int sky_triangle_to_rectangle (Image *image, SkyTriangle *triangle) {
    288 
    289   int i, n, parity, peak, b1, b2, NX, NY, right;
     257int sky_triangle_to_rectangle (SkyRectangle *rectangle, SkyTriangle *triangle) {
     258
     259  int i, parity, peak, b1, b2, NX, NY, right;
    290260  double xv[3], yv[3];        // coordinates of the vertex in the reference projection
    291   double xo, yo, xc, yc, angle, scale;
    292   double Xmin, Xmax, Ymin, Ymax;
     261  double xo, yo, xc, yc, xcr, ycr, angle;
    293262  double dB, dP, s1, s2, r1, r2, dr, dx, dy;
    294263  double angle_b1, angle_b2, slope;
    … …  
    390359  xc = 0.5*(xv[peak] + xo);
    391360  yc = 0.5*(yv[peak] + yo);
    392    
    393   NX = hypot((xv[b2]-xv[b1]),(yv[b2]-yv[b1]));
    394   NY = hypot((xv[peak]-xo),(yv[peak]-yo));
    395 
    396   memset (image, 0, sizeof(Image));
    397   image[0].coords = *refcoords;
    398   image[0].coords.pc1_1 = +cos(angle);
    399   image[0].coords.pc1_2 = +sin(angle);
    400   image[0].coords.pc2_1 = -sin(angle);
    401   image[0].coords.pc2_2 = +cos(angle);
    402 
     361
     362  // NX,NY are the size of the circumscribed square, expanded by PADDING
     363  NX = hypot((xv[b2]-xv[b1]),(yv[b2]-yv[b1])) * (1 + PADDING);
     364  NY = hypot((xv[peak]-xo),(yv[peak]-yo)) * (1 + PADDING);
     365
     366  memset (rectangle, 0, sizeof(SkyRectangle));
     367  rectangle[0].coords = *refcoords;
     368
     369  if (FIX_NS) {
     370    rectangle[0].coords.pc1_1 = +1.0;
     371    rectangle[0].coords.pc1_2 = +0.0;
     372    rectangle[0].coords.pc2_1 = -0.0;
     373    rectangle[0].coords.pc2_2 = +1.0;
     374    xcr = xc*cos(angle) - yc*sin(angle);
     375    ycr = yc*cos(angle) + xc*sin(angle);
     376  } else {
     377    rectangle[0].coords.pc1_1 = +cos(angle);
     378    rectangle[0].coords.pc1_2 = +sin(angle);
     379    rectangle[0].coords.pc2_1 = -sin(angle);
     380    rectangle[0].coords.pc2_2 = +cos(angle);
     381    xcr = xc;
     382    ycr = yc;
     383  }
     384 
    403385  // crpix1,crpix2 is the projection center
    404   image[0].coords.crpix1 = 0.5*NX - xc;
    405   image[0].coords.crpix2 = 0.5*NY - yc;
    406 
    407   // only keep one of the parity images
     386  rectangle[0].coords.crpix1 = 0.5*NX - xcr;
     387  rectangle[0].coords.crpix2 = 0.5*NY - ycr;
     388
     389  // only keep one of the parity rectangles
    408390  if (((triangle[0].d >= 0) && (parity == +1)) || ((triangle[0].d < 0) && (parity == -1))) {
    409     strcpy (image[0].coords.ctype, "DEC--TAN");
     391    strcpy (rectangle[0].coords.ctype, "DEC--TAN");
    410392  } else {
    411     strcpy (image[0].coords.ctype, "DROP");
    412   }
    413   if ((NX > 60000) || (NY > 60000)) {
    414     if (!warned) {
    415       fprintf (stderr, "warning: NX,NY too big for DVO limits; modifying pixel scale\n");
    416       warned = TRUE;
    417     }
    418     scale = NX / 60000.0;
    419     scale = MAX(scale, NY / 60000.0);
    420     NX /= scale;
    421     NY /= scale;
    422     image[0].coords.cdelt1 *= scale;
    423     image[0].coords.cdelt2 *= scale;
    424     image[0].coords.crpix1 /= scale;
    425     image[0].coords.crpix2 /= scale;
    426   }
    427 
    428   image[0].NX = NX;
    429   image[0].NY = NY;
    430   image[0].code = 1; // this needs to be set more sensibly
    431 
     393    strcpy (rectangle[0].coords.ctype, "DROP");
     394  }
     395
     396  rectangle[0].NX = NX;
     397  rectangle[0].NY = NY;
     398  rectangle[0].code = 1; // this needs to be set more sensibly
     399
     400  return (TRUE);
     401}
     402
     403// an allocated image set is supplied, we fill in the values
     404int sky_subdivide_image (Image *output, SkyRectangle *input, int Nx, int Ny) {
     405
     406  int i, j, N, NX, NY, Ndigit;
     407  char format[24];
     408
     409  NX = input[0].NX/(double)Nx + 0.5;
     410  NY = input[0].NY/(double)Ny + 0.5;
     411
     412  // image[0].NX,NY are unsigned short: abort is we overflow
     413  if ((NX > 0xffff) || (NY > 0xffff)) {
     414    fprintf (stderr, "error: NX,NY too big for DVO limits; modifying pixel scale\n");
     415    exit (1);
     416  }
     417
     418  Ndigit = (int)(log10(Nx*Ny)) + 1 ;
     419  snprintf (format, 24, "%s.%%0%dd", input[0].name, Ndigit);
     420
     421  N = 0;
     422  for (j = 0; j < Ny; j++) {
     423    for (i = 0; i < Nx; i++) {
     424
     425      memset (&output[N], 0, sizeof(Image));
     426      memcpy (&output[N].coords, &input[0].coords, sizeof(Coords));
     427
     428      snprintf (output[N].name, 64, format, N);
     429      output[N].NX = NX;
     430      output[N].NY = NY;
     431      output[N].code = input[0].code;
     432
     433      output[N].coords.crpix1 = input[0].coords.crpix1 - i*NX;
     434      output[N].coords.crpix2 = input[0].coords.crpix2 - j*NY;
     435      N++;
     436    }
     437  }
    432438  return (TRUE);
    433439}
    … …  
    520526  // generate 0-level triangles
    521527  ALLOCATE (tri, SkyTriangle, 20);
     528
     529  for (i = 0; i < 20; i++) {
     530    memset (&tri[i], 0, sizeof(SkyTriangle));
     531  }
    522532
    523533  ctht = cos(THETA);
    … …  
    594604}
    595605
     606int sky_base_rotation (SkyTriangle *base, int Nbase) {
     607
     608  // apply the three euler angles (A, B, C)
     609  // XXX for now, just apply A and B
     610
     611  int i, j, ix;
     612  float rot[3][3], v[3];
     613
     614  rot[0][0] = +cos(EULER_A)*cos(EULER_B);
     615  rot[1][0] = +sin(EULER_A)*cos(EULER_B);
     616  rot[2][0] = +sin(EULER_B);
     617
     618  rot[0][1] = -sin(EULER_A);
     619  rot[1][1] = +cos(EULER_A);
     620  rot[2][1] = +0.0;
     621
     622  rot[0][2] = -cos(EULER_A)*sin(EULER_B);
     623  rot[1][2] = -sin(EULER_A)*sin(EULER_B);
     624  rot[2][2] = +cos(EULER_B);
     625
     626  for (i = 0; i < Nbase; i++) {
     627    for (j = 0; j < 3; j++) {
     628      for (ix = 0; ix < 3; ix++) {
     629        v[ix] = 0.0;
     630        v[ix] += base[i].vertex[j].x * rot[0][ix];
     631        v[ix] += base[i].vertex[j].y * rot[1][ix];
     632        v[ix] += base[i].vertex[j].z * rot[2][ix];
     633      }
     634      base[i].vertex[j].x = v[0];
     635      base[i].vertex[j].y = v[1];
     636      base[i].vertex[j].z = v[2];
     637    }
     638  }
     639  return (TRUE);
     640}
     641
    596642int sky_tessalation_init (double scale) {
    597643
    … …  
    611657int dvo_image_clear_vtable (FITS_DB *db) {
    612658
    613   int Nimages;
    614   int i, N, Nx, Ny;
     659  int i;
    615660
    616661  // free memory used by the current vtable rows
Note: See TracChangeset for help on using the changeset viewer.