Changeset 12772

Timestamp:

Apr 6, 2007, 4:42:36 AM (19 years ago)

Author:

eugene

Message:

make rectangles from the triangles

Location:

trunk/Ohana/src/addstar

Files:

• include/skycells.h (modified) (1 diff)
• src/sky_tessalation.c (modified) (3 diffs)

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

@@ -64 +64 @@
 SkyTriangle *sky_base_triangles (int *ntriangles);
 int sky_tessalation_init ();
+int sky_triangle_to_rectangle (Image *image, SkyTriangle *triangle);

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

@@ -6 +6 @@
 Image *sky_tessalation (int level, int *nimages) {
 
-  int i, Ntriangles;
+  int i, N, Ntriangles;
   SkyTriangle *tri, *new;
-  Image *image;
+  Image *image, *outimage;
 
   sky_tessalation_init ();
@@ -22 +22 @@
   ALLOCATE (image, Image, Ntriangles);
   for (i = 0; i < Ntriangles; i++) {
-    sky_triangle_to_image (&image[i], &tri[i]);
-    sprintf (image[i].name, "test.%03d", i);
+    // XXX add a user-option to choose betwen these
+    // sky_triangle_to_image (&image[i], &tri[i]);
+    sky_triangle_to_rectangle (&image[i], &tri[i]);
   }  
 
-  *nimages = Ntriangles;
-  return (image);
+  FILE *f1 = fopen ("par.1.dat", "w");
+  FILE *f2 = fopen ("par.2.dat", "w");
+
+  ALLOCATE (outimage, Image, Ntriangles);
+  for (i = N = 0; i < Ntriangles; i++) {
+    if (!strcmp(image[i].coords.ctype, "DROP")) { 
+      fprintf (f2, "%f %f\n", image[i].coords.crval1, image[i].coords.crval2);
+      continue;
+    } else {
+      fprintf (f1, "%f %f\n", image[i].coords.crval1, image[i].coords.crval2);
+    }
+    memcpy (&outimage[N], &image[i], sizeof(Image));
+    sprintf (outimage[N].name, "test.%03d", i);
+    N++;
+  }  
+  free (image);
+
+  fclose (f1);
+  fclose (f2);
+
+  *nimages = N;
+  return (outimage);
 }
 
@@ -91 +112 @@
   image[0].Mxxx = xv[1];  image[0].Mxyy = yv[1];
   image[0].Mxxy = xv[2];  image[0].Myyy = yv[2];
+
+  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;
+  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;
+
+  // calculate the triangle coordinates in r,d
+  sky_triangle_coords (triangle);
+
+  // 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]);
+
+  // xc, yc is the true image center
+  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));
+  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);
+
+  // 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?
+  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) {
+    strcpy (image[0].coords.ctype, "RA---TAN");
+  } else {
+    strcpy (image[0].coords.ctype, "DROP");
+  }
+  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].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
 
   return (TRUE);