Changeset 12775

Timestamp:

Apr 6, 2007, 5:47:41 PM (19 years ago)

Author:

eugene

Message:

improvements to skycells

Location:

trunk/Ohana/src/addstar

Files:

• include/skycells.h (modified) (2 diffs)
• src/args_skycells.c (modified) (2 diffs)
• src/gstars.c (modified) (1 diff)
• src/sky_tessalation.c (modified) (15 diffs)
• src/skycells.c (modified) (2 diffs)

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

@@ -52 +52 @@
 int    VERBOSE;
 int    MODE;
+double SCALE;
 
 void       SetProtect             PROTO((int mode));
@@ -65 +66 @@
 Point sky_divide_edge (Point v1, Point v2);
 SkyTriangle *sky_base_triangles (int *ntriangles);
-int sky_tessalation_init ();
+int sky_tessalation_init (double scale);
 int sky_triangle_to_rectangle (Image *image, SkyTriangle *triangle);
 
-Image *sky_tessalation (int level, int *nimages, int mode);
-Image *sky_tessalation_triangles (int level, int *nimages);
-Image *sky_tessalation_squares (int level, int *nimages);
+Image *sky_tessalation (int level, int *nimages, int mode, double scale);
+Image *sky_tessalation_triangles (int level, int *nimages, double scale);
+Image *sky_tessalation_squares (int level, int *nimages, double scale);

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

@@ -5 +5 @@
   
   int i, N;
+  char *ptr;
 
   /* check for help request */
@@ -26 +27 @@
   }
 
-  if (argc == 2) return (TRUE);
+  /* what type of output files? */
+  SCALE = 1.0;
+  if ((N = get_argument (argc, argv, "-scale"))) {
+    remove_argument (N, &argc, argv);
+    SCALE = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  if (argc == 2) {
+    strtol (argv[1], &ptr, 10);
+    if (*ptr == 0) return (TRUE);
+  }
 
   fprintf (stderr, "USAGE: skycells (level) [-triangles] [-scale arcsec/pix]\n");

trunk/Ohana/src/addstar/src/gstars.c

@@ -1 +1 @@
 # include "addstar.h"
-
-// XXX where and when should I call RegisterMosaic??
 
 Stars *gstars (char *filename, int *Nstars, Image **images, int *Nimages, int photcode) {

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

@@ -1 +1 @@
 # include "skycells.h"
 # include "assert.h"
+# define iSWAP(X,Y) {int tmp=(X); (X) = (Y); (Y) = tmp;}
 
 static Coords *refcoords = NULL;
-
-Image *sky_tessalation (int level, int *nimages, int mode) {
+static int warned = FALSE;
+
+Image *sky_tessalation (int level, int *nimages, int mode, double scale) {
 
   Image *images;
 
   if (mode == SQUARES) {
-    images = sky_tessalation_squares (level, nimages);
+    images = sky_tessalation_squares (level, nimages, scale);
     return images;
   }
 
   if (mode == TRIANGLES) {
-    images = sky_tessalation_triangles (level, nimages);
+    images = sky_tessalation_triangles (level, nimages, scale);
     return images;
   }
@@ -21 +23 @@
 }
 
-Image *sky_tessalation_triangles (int level, int *nimages) {
+Image *sky_tessalation_triangles (int level, int *nimages, double scale) {
 
   int i, Ndigit, Ntriangles;
@@ -31 +33 @@
   snprintf (format, 16, "skytri.%%0%dd", Ndigit);
 
-  sky_tessalation_init ();
+  sky_tessalation_init (scale);
 
   tri = sky_base_triangles (&Ntriangles);
@@ -51 +53 @@
 }
 
-Image *sky_tessalation_squares (int level, int *nimages) {
+Image *sky_tessalation_squares (int level, int *nimages, double scale) {
 
   int i, N, Ndigit, Ntriangles;
@@ -62 +64 @@
   snprintf (format, 16, "skycell.%%0%dd", Ndigit);
 
-  sky_tessalation_init ();
+  sky_tessalation_init (scale);
 
   tri = sky_base_triangles (&Ntriangles);
@@ -77 +79 @@
   }  
 
+  // FILE *f1 = fopen ("center.keep.dat", "w");
+  // FILE *f2 = fopen ("center.drop.dat", "w");
+
   ALLOCATE (outimage, Image, Ntriangles);
   for (i = N = 0; i < Ntriangles; i++) {
-    if (!strcmp(image[i].coords.ctype, "DROP")) continue;
+    if (!strcmp(image[i].coords.ctype, "DROP")) {
+      // fprintf (f2, "%f %f\n", tri[i].r, tri[i].d);
+      continue;
+    } else {
+      // fprintf (f1, "%f %f\n", tri[i].r, tri[i].d);
+    }
     memcpy (&outimage[N], &image[i], sizeof(Image));
     snprintf (outimage[N].name, 32, format, i);
@@ -85 +95 @@
   }  
   free (image);
+
+  // fclose (f1);
+  // fclose (f2);
 
   *nimages = N;
@@ -159 +172 @@
 int sky_triangle_to_rectangle (Image *image, SkyTriangle *triangle) {
 
-  int i, n, parity, peak, b1, b2, NX, NY;
+  int i, n, 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 dB, dP, s1, s2, r1, r2, dr, dx, dy;
+  double angle_b1, angle_b2, slope;
 
   // calculate the triangle coordinates in r,d
@@ -195 +209 @@
   b1 = (peak + 1) % 3;
   b2 = (peak + 2) % 3;
+
+  // angle is from the center to the peak corner
+  angle_b1 = DEG_RAD*atan(yv[b1] / xv[b1]);
+  angle_b2 = DEG_RAD*atan(yv[b2] / xv[b2]);
+
+  // if one of the base-center angles is very small, the parity is marginal.  Use additional
+  // information to choose the parity.  note that both angle_b1 and angle_b2 cannot be close to
+  // zero.
+  if (fabs(angle_b1) < 10.0) {
+    right = (xv[b1] > 0);     // pointing left or right?
+    slope = (xv[peak] -  xv[b2]) / (yv[peak] - yv[b2]);
+    if ( right && (slope >= 0.0)) parity = +1;
+    if ( right && (slope <  0.0)) parity = -1;
+    if (!right && (slope <= 0.0)) parity = +1;
+    if (!right && (slope >  0.0)) parity = -1;
+    if (parity > 0) {
+      if (yv[peak] < yv[b2]) iSWAP(peak, b2); // require peak to be top (bottom) point
+    } else {
+      if (yv[peak] > yv[b2]) iSWAP(peak, b2); // require peak to be top (bottom) point
+    }
+    angle = atan2(parity*xv[peak], parity*yv[peak]); // note that this is x/y not y/x (and in radians)
+  }
+  if (fabs(angle_b2) < 10.0) {
+    right = (xv[b2] > 0);     // pointing left or right?
+    slope = (xv[peak] - xv[b1]) / (yv[peak] - yv[b1]); // tilt of opposite line
+    if ( right && (slope >= 0.0)) parity = +1;
+    if ( right && (slope <  0.0)) parity = -1;
+    if (!right && (slope <= 0.0)) parity = +1;
+    if (!right && (slope >  0.0)) parity = -1;
+    if (parity > 0) {
+      if (yv[peak] < yv[b1]) iSWAP(peak, b1); // require peak to be top (bottom) point
+    } else {
+      if (yv[peak] > yv[b1]) iSWAP(peak, b1); // require peak to be top (bottom) point
+    }
+    angle = atan2(parity*xv[peak], parity*yv[peak]); // note that this is x/y not y/x (and in radians)
+  }
 
   // xo, yo is the center of the baseline
@@ -229 +279 @@
   NY = hypot((xv[peak]-xo),(yv[peak]-yo));
 
-
   memset (image, 0, sizeof(Image));
   image[0].coords = *refcoords;
@@ -237 +286 @@
   image[0].coords.pc2_2 = +cos(angle);
 
-  // XXX crpix1,2 should be one of the corners for RA---TAN!!!
-  // 0.5NX + xc, 0.5NY + yc or something like that 
-  // is this the right sign?
+  // 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
-  if (parity == +1) {
+  if (((triangle[0].d >= 0) && (parity == +1)) || ((triangle[0].d < 0) && (parity == -1))) {
     strcpy (image[0].coords.ctype, "RA---TAN");
   } else {
@@ -250 +297 @@
   }
   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);
@@ -317 +368 @@
       out[4*i + 3].vertex[j] = out[4*i + j].vertex[1];
     }
-# if (0)
-    fprintf (stdout, "%f %f %f  %f %f %f  %f %f %f\n", 
-             in[i].vertex[0].x, in[i].vertex[0].y, in[i].vertex[0].z, 
-             in[i].vertex[1].x, in[i].vertex[1].y, in[i].vertex[1].z, 
-             in[i].vertex[2].x, in[i].vertex[2].y, in[i].vertex[2].z);
-    for (j = 0; j < 4; j++) {
-      fprintf (stdout, "%f %f %f  %f %f %f  %f %f %f\n", 
-               out[4*i+j].vertex[0].x, out[4*i+j].vertex[0].y, out[4*i+j].vertex[0].z, 
-               out[4*i+j].vertex[1].x, out[4*i+j].vertex[1].y, out[4*i+j].vertex[1].z, 
-               out[4*i+j].vertex[2].x, out[4*i+j].vertex[2].y, out[4*i+j].vertex[2].z);
-    }
-# endif  
     Nt += 4;
   }
@@ -440 +479 @@
 }
 
-int sky_tessalation_init () {
+int sky_tessalation_init (double scale) {
 
   ALLOCATE (refcoords, Coords, 1);
   refcoords[0].crval1 = refcoords[0].crval2 = 0.0;
   refcoords[0].crpix1 = refcoords[0].crpix2 = 0.0;
-  refcoords[0].cdelt1 = refcoords[0].cdelt2 = 1.0 / 3600;
+  refcoords[0].cdelt1 = refcoords[0].cdelt2 = scale / 3600;
   refcoords[0].pc1_1 = refcoords[0].pc2_2 = 1.0;
   refcoords[0].pc1_2 = refcoords[0].pc2_1 = 0.0;
@@ -453 +492 @@
   return (TRUE);
 }
-
-
-
-/*** this code was an attempt to fit the triangle boundaries into NX, NY ***/
-# if (0)
-  for (n = 0; n < 3; n++) {
-    // we will project to the triangle center position
-    refcoords[0].crval1 = triangle[0].r;
-    refcoords[0].crval2 = triangle[0].d;
-
-    // find the size, rotation, and parity of the image
-    // project the vertices and find the image parity
-    parity = 1;
-    for (i = 0; i < 3; i++) {
-      RD_to_XY (&xv[i], &yv[i], triangle[0].rv[i], triangle[0].dv[i], refcoords);
-      parity *= SIGN(yv[i]);
-    }
-
-    // choose the peak vertex
-    peak = -1;
-    for (i = 0; (peak == -1) && (i < 3); i++) {
-      if (parity == SIGN(yv[i])) {
-        peak = i;
-      }
-    }
-    assert (peak != -1);
-
-    // find the base and height
-    b1 = (peak + 1) % 3;
-    b2 = (peak + 2) % 3;
-
-    // xo, yo is the center of the baseline
-    xo = 0.5*(xv[b2] + xv[b1]);
-    yo = 0.5*(yv[b2] + yv[b1]);
-
-    if (n != 2) {
-      // xc, yc is the true image center
-      xc = 0.5*(xv[peak] + xo);
-      yc = 0.5*(yv[peak] + yo);
-      XY_to_RD (&triangle[0].r, &triangle[0].d, xc, yc, refcoords);
-    }
-  }  
-
-  NX = hypot((xv[b2]-xv[b1]),(yv[b2]-yv[b1]));
-  NY = hypot((xv[peak]-xo),(yv[peak]-yo));
-  angle = atan2(parity*xv[peak], parity*yv[peak]); // note that this is x/y not y/x (and in radians)
-
-  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);
-
-  if (parity == 1) {
-    strcpy (image[0].coords.ctype, "TRP--TAN");
-  } else {
-    strcpy (image[0].coords.ctype, "TRM--TAN");
-  }
-  sprintf (image[0].name, "test.%03d", i);
-  if ((NX > 60000) || (NY > 60000)) {
-    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].NX = NX;
-  image[0].NY = NY;
-  image[0].code = 1;
-# endif

trunk/Ohana/src/addstar/src/skycells.c

@@ -8 +8 @@
   FITS_DB db;
 
-  // XXX for test data
-  Coords coords;
-  double r, d;
-
   SetSignals ();
   ConfigInit_skycells (&argc, argv);
   args_skycells (argc, argv);
   level = atoi (argv[1]);
-
-  images = sky_tessalation (level, &Nimages, MODE);
+  
+  images = sky_tessalation (level, &Nimages, MODE, SCALE);
 
   /*** update the image table ***/
@@ -45 +41 @@
   exit (0);
 }
-
-# if (0)
-  coords.crval1 = coords.crval2 = 0.0;
-  coords.crpix1 = coords.crpix2 = 0.0;
-  coords.cdelt1 = coords.cdelt2 = 1.0 / 3600;
-  coords.pc1_1 = coords.pc2_2 = 1.0;
-  coords.pc1_2 = coords.pc2_1 = 0.0;
-  coords.Npolyterms = 0;
-  memset (coords.polyterms, 0, 14*sizeof(float));
-  strcpy (coords.ctype, "TRP--TAN");
-
-  /** quick demo to create a set of triangular images **/
-  Nimages = 100;
-  ALLOCATE (images, Image, Nimages);
-  for (r = 0, i = 0; r < 10; r++) {
-      for (d = 0; d < 10; d++, i++) {
-          memset (&images[i], 0, sizeof(Image));
-          images[i].coords = coords;
-          images[i].coords.crval1 = r;
-          images[i].coords.crval2 = d;
-          if (i%2) strcpy (images[i].coords.ctype, "TRM--TAN");
-          sprintf (images[i].name, "test.%03d", i);
-          images[i].NX = 1000;
-          images[i].NY = 1000;
-          images[i].code = 1;
-      }
-  }
-# endif
-