Changeset 12775 for trunk/Ohana/src/addstar/src/sky_tessalation.c

Timestamp:

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

Author:

eugene

Message:

improvements to skycells

File:

• trunk/Ohana/src/addstar/src/sky_tessalation.c (modified) (15 diffs)

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