Index: /trunk/Ohana/src/addstar/include/skycells.h =================================================================== --- /trunk/Ohana/src/addstar/include/skycells.h (revision 15996) +++ /trunk/Ohana/src/addstar/include/skycells.h (revision 15997) @@ -24,4 +24,12 @@ double rv[3], dv[3]; // triangle center (2d) } SkyTriangle; + +typedef struct { + Coords coords; + int NX; + int NY; + int code; + char name[64]; +} SkyRectangle; /* globals which define database info / data sources (KEEP) */ @@ -52,27 +60,33 @@ int VERBOSE; int MODE; +int FIX_NS; int NMAX; int NX_SUB, NY_SUB; double SCALE; +double PADDING; -void SetProtect PROTO((int mode)); -int SetSignals PROTO(()); -int Shutdown PROTO((char *message, ...); ) -void TrapSignal PROTO((int sig)); -int args_skycells (int argc, char **argv); -int ConfigInit_skycells (int *argc, char **argv); +double EULER_A; +double EULER_B; -int sky_triangle_to_image (Image *image, SkyTriangle *triangle); -int sky_triangle_coords (SkyTriangle *triangle); -SkyTriangle *sky_divide_triangles (SkyTriangle *in, int *ntriangles); -Point sky_divide_edge (Point v1, Point v2); -SkyTriangle *sky_base_triangles (int *ntriangles); -int sky_tessalation_init (double scale); -int sky_triangle_to_rectangle (Image *image, SkyTriangle *triangle); +void SetProtect PROTO((int mode)); +int SetSignals PROTO(()); +int Shutdown PROTO((char *message, ...); ) +void TrapSignal PROTO((int sig)); -int sky_tessalation (FITS_DB *db, int level, int Nmax, int mode, double scale); -int sky_tessalation_triangles (FITS_DB *db, int level, int Nmax); -int sky_tessalation_squares (FITS_DB *db, int level, int Nmax); +int args_skycells PROTO((int argc, char **argv)); +int ConfigInit_skycells PROTO((int *argc, char **argv)); +int sky_tessalation PROTO((FITS_DB *db, int level, int Nmax, int mode, double scale)); +int sky_tessalation_init PROTO((double scale)); +int sky_tessalation_triangles PROTO((FITS_DB *db, int level, int Nmax)); +int sky_tessalation_squares PROTO((FITS_DB *db, int level, int Nmax)); +int sky_triangle_to_image PROTO((Image *image, SkyTriangle *triangle)); +int sky_triangle_to_rectangle PROTO((SkyRectangle *image, SkyTriangle *triangle)); +int sky_subdivide_image PROTO((Image *output, SkyRectangle *input, int Nx, int Ny)); +int sky_triangle_coords PROTO((SkyTriangle *triangle)); +SkyTriangle *sky_divide_triangles PROTO((SkyTriangle *in, int *ntriangles)); +SkyTriangle *sky_base_triangles PROTO((int *ntriangles)); +int sky_base_rotation PROTO((SkyTriangle *base, int Nbase)); +Point sky_divide_edge PROTO((Point v1, Point v2)); -// migrate to libdvo eventually -int dvo_image_clear_vtable (FITS_DB *db); +// XXX migrate to libdvo eventually +int dvo_image_clear_vtable PROTO((FITS_DB *db)); Index: /trunk/Ohana/src/addstar/src/args_skycells.c =================================================================== --- /trunk/Ohana/src/addstar/src/args_skycells.c (revision 15996) +++ /trunk/Ohana/src/addstar/src/args_skycells.c (revision 15997) @@ -4,5 +4,5 @@ int args_skycells (int argc, char **argv) { - int i, N; + int N; char *ptr; @@ -27,4 +27,11 @@ } + /* what type of output files? */ + FIX_NS = FALSE; + if ((N = get_argument (argc, argv, "-fix-ns"))) { + FIX_NS = TRUE; + remove_argument (N, &argc, argv); + } + /* pixel scale (arcsec/pixel) */ SCALE = 1.0; @@ -32,4 +39,22 @@ remove_argument (N, &argc, argv); SCALE = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + /* pixel scale (arcsec/pixel) */ + EULER_A = EULER_B = 0.0; + if ((N = get_argument (argc, argv, "-euler"))) { + remove_argument (N, &argc, argv); + EULER_A = RAD_DEG*atof (argv[N]); + remove_argument (N, &argc, argv); + EULER_B = RAD_DEG*atof (argv[N]); + remove_argument (N, &argc, argv); + } + + /* pixel scale (arcsec/pixel) */ + PADDING = 0.0; + if ((N = get_argument (argc, argv, "-padding"))) { + remove_argument (N, &argc, argv); + PADDING = atof (argv[N]); remove_argument (N, &argc, argv); } @@ -61,5 +86,6 @@ } - fprintf (stderr, "USAGE: skycells (level) [-triangles] [-scale arcsec/pix] [-nmax Max Cells in memory] [-nx (Nx cells)] [-ny (Ny cells)]\n"); + fprintf (stderr, "USAGE: skycells (level) [-scale arcsec/pix] [-nx (Nx cells)] [-ny (Ny cells)]\n"); + fprintf (stderr, " [-h for details and other options]\n"); exit (2); } @@ -71,4 +97,11 @@ fprintf (stderr, " -v : verbose mode\n"); + fprintf (stderr, " -triangles : save base triangles instead of skycells\n"); + fprintf (stderr, " -fix-ns : orient skycells with y-axis aligned with Dec\n"); + fprintf (stderr, " -scale : set pixel scale (default 1.0 arcsec / pixel)\n"); + fprintf (stderr, " -padding : pad skycells by this fraction in each dimension\n"); + fprintf (stderr, " -nmax : only keep nmax skycells in memory\n"); + fprintf (stderr, " -nx : subdivide skycell projection in x\n"); + fprintf (stderr, " -ny : subdivide skycell projection in y\n"); fprintf (stderr, " -help : this list\n"); fprintf (stderr, " -h : this list\n\n"); Index: /trunk/Ohana/src/addstar/src/sky_tessalation.c =================================================================== --- /trunk/Ohana/src/addstar/src/sky_tessalation.c (revision 15996) +++ /trunk/Ohana/src/addstar/src/sky_tessalation.c (revision 15997) @@ -3,6 +3,6 @@ # define iSWAP(X,Y) {int tmp=(X); (X) = (Y); (Y) = tmp;} +// we use a static refcoords structure to avoid multiple alloc / init steps static Coords *refcoords = NULL; -static int warned = FALSE; int sky_tessalation (FITS_DB *db, int level, int Nmax, int mode, double scale) { @@ -33,4 +33,6 @@ // generate the initial base set base = sky_base_triangles (&Nbase); + + sky_base_rotation (base, Nbase); // how many triangles total for this level? @@ -95,8 +97,9 @@ int sky_tessalation_squares (FITS_DB *db, int level, int Nmax) { - int i, j, Nname, Ndigit, Ntriangles, Nbase, Ntotal, Ntop, Ltop, Nout, Nx, Ny; + int i, j, Nname, Ndigit, Ntriangles, Nbase, Nimage, Ntotal, Ntop, Ltop, Nsubset, Nx, Ny; double fLtop; SkyTriangle *base, *tri, *new; - Image *image, *out; + SkyRectangle *rectangle, *subset; + Image *image; char format[16]; @@ -106,4 +109,6 @@ // generate the initial base set base = sky_base_triangles (&Nbase); + + sky_base_rotation (base, Nbase); // how many total cells for this level (multiply by subdivisions, if used)? @@ -144,34 +149,31 @@ } - // convert the SkyTriangles to Image - ALLOCATE (image, Image, Ntriangles); + // convert the SkyTriangles to SkyRectangles + ALLOCATE (rectangle, SkyRectangle, Ntriangles); for (j = 0; j < Ntriangles; j++) { - sky_triangle_to_rectangle (&image[j], &tri[j]); + sky_triangle_to_rectangle (&rectangle[j], &tri[j]); } // drop the appropriate subset - ALLOCATE (out, Image, Ntriangles); - for (j = Nout = 0; j < Ntriangles; j++) { - if (!strcmp(image[j].coords.ctype, "DROP")) continue; - memcpy (&out[Nout], &image[j], sizeof(Image)); - snprintf (out[Nout].name, 64, format, Nname); + ALLOCATE (subset, SkyRectangle, Ntriangles); + for (j = Nsubset = 0; j < Ntriangles; j++) { + if (!strcmp(rectangle[j].coords.ctype, "DROP")) continue; + memcpy (&subset[Nsubset], &rectangle[j], sizeof(SkyRectangle)); + snprintf (subset[Nsubset].name, 64, format, Nname); Nname++; - Nout++; + Nsubset++; } - free (image); + free (rectangle); // subdivide each image (Nx x Ny subcells) - if (Nx*Ny > 1) { - ALLOCATE (image, Image, Nx*Ny*Nout); - for (j = 0; j < Nout; j++) { - sky_subdivide_image (&image[j*Nx*Ny], &out[j], Nx, Ny); - } - free (out); - out = image; - Nout *= Nx*Ny; + Nimage = Nx*Ny*Nsubset; + ALLOCATE (image, Image, Nimage); + for (j = 0; j < Nsubset; j++) { + // convert the SkyRectangles to Images for output + sky_subdivide_image (&image[j*Nx*Ny], &subset[j], Nx, Ny); } /* add the new images and save */ - dvo_image_addrows (db, out, Nout); + dvo_image_addrows (db, image, Nimage); SetProtect (TRUE); dvo_image_update (db, VERBOSE); @@ -179,5 +181,6 @@ dvo_image_clear_vtable (db); - free (out); + free (subset); + free (image); free (tri); } @@ -188,7 +191,7 @@ int sky_triangle_to_image (Image *image, SkyTriangle *triangle) { - int i, n, parity, peak, b1, b2, NX, NY; + int i, NX, NY; double xv[3], yv[3]; // coordinates of the vertex in the reference projection - double xo, yo, xc, yc, angle, scale; + double scale; double Xmin, Xmax, Ymin, Ymax; @@ -250,45 +253,11 @@ } -// an allocated image set is supplied, we fill in the values -int sky_subdivide_image (Image *output, Image *input, int Nx, int Ny) { - - int i, j, N, NX, NY, Ndigit; - char format[24]; - - NX = input[0].NX/(double)Nx + 0.5; - NY = input[0].NY/(double)Ny + 0.5; - - Ndigit = (int)(log10(Nx*Ny)) + 1 ; - snprintf (format, 24, "%s.%%0%dd", input[0].name, Ndigit); - - N = 0; - for (j = 0; j < Ny; j++) { - for (i = 0; i < Nx; i++) { - - memcpy (&output[N], &input[0], sizeof(Image)); - snprintf (output[N].name, 64, format, N); - - // output[0].coords = input[0].coords; - // strcpy (output[0].coords.ctype, input[0].coords.ctype); - - output[N].NX = NX; - output[N].NY = NY; - - output[N].coords.crpix1 = input[0].coords.crpix1 - i*NX; - output[N].coords.crpix2 = input[0].coords.crpix2 - j*NY; - N++; - } - } - return (TRUE); -} - // an allocated image is supplied, we fill in the values // we are only keeping ~half of the images -int sky_triangle_to_rectangle (Image *image, SkyTriangle *triangle) { - - int i, n, parity, peak, b1, b2, NX, NY, right; +int sky_triangle_to_rectangle (SkyRectangle *rectangle, SkyTriangle *triangle) { + + int i, parity, peak, b1, b2, NX, NY, right; double xv[3], yv[3]; // coordinates of the vertex in the reference projection - double xo, yo, xc, yc, angle, scale; - double Xmin, Xmax, Ymin, Ymax; + double xo, yo, xc, yc, xcr, ycr, angle; double dB, dP, s1, s2, r1, r2, dr, dx, dy; double angle_b1, angle_b2, slope; @@ -390,44 +359,81 @@ xc = 0.5*(xv[peak] + xo); yc = 0.5*(yv[peak] + yo); - - NX = hypot((xv[b2]-xv[b1]),(yv[b2]-yv[b1])); - NY = hypot((xv[peak]-xo),(yv[peak]-yo)); - - memset (image, 0, sizeof(Image)); - image[0].coords = *refcoords; - image[0].coords.pc1_1 = +cos(angle); - image[0].coords.pc1_2 = +sin(angle); - image[0].coords.pc2_1 = -sin(angle); - image[0].coords.pc2_2 = +cos(angle); - + + // NX,NY are the size of the circumscribed square, expanded by PADDING + NX = hypot((xv[b2]-xv[b1]),(yv[b2]-yv[b1])) * (1 + PADDING); + NY = hypot((xv[peak]-xo),(yv[peak]-yo)) * (1 + PADDING); + + memset (rectangle, 0, sizeof(SkyRectangle)); + rectangle[0].coords = *refcoords; + + if (FIX_NS) { + rectangle[0].coords.pc1_1 = +1.0; + rectangle[0].coords.pc1_2 = +0.0; + rectangle[0].coords.pc2_1 = -0.0; + rectangle[0].coords.pc2_2 = +1.0; + xcr = xc*cos(angle) - yc*sin(angle); + ycr = yc*cos(angle) + xc*sin(angle); + } else { + rectangle[0].coords.pc1_1 = +cos(angle); + rectangle[0].coords.pc1_2 = +sin(angle); + rectangle[0].coords.pc2_1 = -sin(angle); + rectangle[0].coords.pc2_2 = +cos(angle); + xcr = xc; + ycr = yc; + } + // crpix1,crpix2 is the projection center - image[0].coords.crpix1 = 0.5*NX - xc; - image[0].coords.crpix2 = 0.5*NY - yc; - - // only keep one of the parity images + rectangle[0].coords.crpix1 = 0.5*NX - xcr; + rectangle[0].coords.crpix2 = 0.5*NY - ycr; + + // only keep one of the parity rectangles if (((triangle[0].d >= 0) && (parity == +1)) || ((triangle[0].d < 0) && (parity == -1))) { - strcpy (image[0].coords.ctype, "DEC--TAN"); + strcpy (rectangle[0].coords.ctype, "DEC--TAN"); } else { - strcpy (image[0].coords.ctype, "DROP"); - } - if ((NX > 60000) || (NY > 60000)) { - if (!warned) { - fprintf (stderr, "warning: NX,NY too big for DVO limits; modifying pixel scale\n"); - warned = TRUE; - } - scale = NX / 60000.0; - scale = MAX(scale, NY / 60000.0); - NX /= scale; - NY /= scale; - image[0].coords.cdelt1 *= scale; - image[0].coords.cdelt2 *= scale; - image[0].coords.crpix1 /= scale; - image[0].coords.crpix2 /= scale; - } - - image[0].NX = NX; - image[0].NY = NY; - image[0].code = 1; // this needs to be set more sensibly - + strcpy (rectangle[0].coords.ctype, "DROP"); + } + + rectangle[0].NX = NX; + rectangle[0].NY = NY; + rectangle[0].code = 1; // this needs to be set more sensibly + + return (TRUE); +} + +// an allocated image set is supplied, we fill in the values +int sky_subdivide_image (Image *output, SkyRectangle *input, int Nx, int Ny) { + + int i, j, N, NX, NY, Ndigit; + char format[24]; + + NX = input[0].NX/(double)Nx + 0.5; + NY = input[0].NY/(double)Ny + 0.5; + + // image[0].NX,NY are unsigned short: abort is we overflow + if ((NX > 0xffff) || (NY > 0xffff)) { + fprintf (stderr, "error: NX,NY too big for DVO limits; modifying pixel scale\n"); + exit (1); + } + + Ndigit = (int)(log10(Nx*Ny)) + 1 ; + snprintf (format, 24, "%s.%%0%dd", input[0].name, Ndigit); + + N = 0; + for (j = 0; j < Ny; j++) { + for (i = 0; i < Nx; i++) { + + memset (&output[N], 0, sizeof(Image)); + memcpy (&output[N].coords, &input[0].coords, sizeof(Coords)); + + snprintf (output[N].name, 64, format, N); + output[N].NX = NX; + output[N].NY = NY; + output[N].code = input[0].code; + + output[N].coords.crpix1 = input[0].coords.crpix1 - i*NX; + output[N].coords.crpix2 = input[0].coords.crpix2 - j*NY; + N++; + } + } return (TRUE); } @@ -520,4 +526,8 @@ // generate 0-level triangles ALLOCATE (tri, SkyTriangle, 20); + + for (i = 0; i < 20; i++) { + memset (&tri[i], 0, sizeof(SkyTriangle)); + } ctht = cos(THETA); @@ -594,4 +604,40 @@ } +int sky_base_rotation (SkyTriangle *base, int Nbase) { + + // apply the three euler angles (A, B, C) + // XXX for now, just apply A and B + + int i, j, ix; + float rot[3][3], v[3]; + + rot[0][0] = +cos(EULER_A)*cos(EULER_B); + rot[1][0] = +sin(EULER_A)*cos(EULER_B); + rot[2][0] = +sin(EULER_B); + + rot[0][1] = -sin(EULER_A); + rot[1][1] = +cos(EULER_A); + rot[2][1] = +0.0; + + rot[0][2] = -cos(EULER_A)*sin(EULER_B); + rot[1][2] = -sin(EULER_A)*sin(EULER_B); + rot[2][2] = +cos(EULER_B); + + for (i = 0; i < Nbase; i++) { + for (j = 0; j < 3; j++) { + for (ix = 0; ix < 3; ix++) { + v[ix] = 0.0; + v[ix] += base[i].vertex[j].x * rot[0][ix]; + v[ix] += base[i].vertex[j].y * rot[1][ix]; + v[ix] += base[i].vertex[j].z * rot[2][ix]; + } + base[i].vertex[j].x = v[0]; + base[i].vertex[j].y = v[1]; + base[i].vertex[j].z = v[2]; + } + } + return (TRUE); +} + int sky_tessalation_init (double scale) { @@ -611,6 +657,5 @@ int dvo_image_clear_vtable (FITS_DB *db) { - int Nimages; - int i, N, Nx, Ny; + int i; // free memory used by the current vtable rows