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

Timestamp:

Apr 6, 2007, 6:07:30 AM (19 years ago)

Author:

eugene

Message:

adding triangle vs square options, expanding square to outer points

File:

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

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

@@ -4 +4 @@
 static Coords *refcoords = NULL;
 
-Image *sky_tessalation (int level, int *nimages) {
-
-  int i, N, Ntriangles;
+Image *sky_tessalation (int level, int *nimages, int mode) {
+
+  Image *images;
+
+  if (mode == SQUARES) {
+    images = sky_tessalation_squares (level, nimages);
+    return images;
+  }
+
+  if (mode == TRIANGLES) {
+    images = sky_tessalation_triangles (level, nimages);
+    return images;
+  }
+
+  return NULL;
+}
+
+Image *sky_tessalation_triangles (int level, int *nimages) {
+
+  int i, Ndigit, Ntriangles;
   SkyTriangle *tri, *new;
-  Image *image, *outimage;
+  Image *image;
+  char format[16];
+
+  Ndigit = (int)(log10(20*pow(4.0, level))) + 1 ;
+  snprintf (format, 16, "skytri.%%0%dd", Ndigit);
 
   sky_tessalation_init ();
@@ -22 +43 @@
   ALLOCATE (image, Image, Ntriangles);
   for (i = 0; i < Ntriangles; i++) {
-    // XXX add a user-option to choose betwen these
-    // sky_triangle_to_image (&image[i], &tri[i]);
+    sky_triangle_to_image (&image[i], &tri[i]);
+    snprintf (image[i].name, 32, format, i);
+  }  
+
+  *nimages = Ntriangles;
+  return (image);
+}
+
+Image *sky_tessalation_squares (int level, int *nimages) {
+
+  int i, N, Ndigit, Ntriangles;
+  SkyTriangle *tri, *new;
+  Image *image, *outimage;
+  char format[16];
+
+  // XXX I need to adjust for the number of cell subdivisions
+  Ndigit = (int)(log10(20*pow(4.0, level))) + 1 ;
+  snprintf (format, 16, "skycell.%%0%dd", Ndigit);
+
+  sky_tessalation_init ();
+
+  tri = sky_base_triangles (&Ntriangles);
+
+  for (i = 0; i < level; i++) {
+    new = sky_divide_triangles (tri, &Ntriangles);
+    free (tri);
+    tri = new;
+  }
+
+  ALLOCATE (image, Image, Ntriangles);
+  for (i = 0; i < Ntriangles; i++) {
     sky_triangle_to_rectangle (&image[i], &tri[i]);
   }  
 
-  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);
-    }
+    if (!strcmp(image[i].coords.ctype, "DROP")) continue;
     memcpy (&outimage[N], &image[i], sizeof(Image));
-    sprintf (outimage[N].name, "test.%03d", i);
+    snprintf (outimage[N].name, 32, format, i);
     N++;
   }  
   free (image);
-
-  fclose (f1);
-  fclose (f2);
 
   *nimages = N;
@@ -124 +163 @@
   double xo, yo, xc, yc, angle, scale;
   double Xmin, Xmax, Ymin, Ymax;
+  double dB, dP, s1, s2, r1, r2, dr, dx, dy;
 
   // calculate the triangle coordinates in r,d
@@ -149 +189 @@
   assert (peak != -1);
 
+  // angle is from the center to the peak corner
+  angle = atan2(parity*xv[peak], parity*yv[peak]); // note that this is x/y not y/x (and in radians)
+
   // find the base and height
   b1 = (peak + 1) % 3;
@@ -156 +199 @@
   xo = 0.5*(xv[b2] + xv[b1]);
   yo = 0.5*(yv[b2] + yv[b1]);
+
+  // find the max perpendicular distance from the peak
+
+  // dB[b1] == dB[b2] (since xo,yo is the midpoint of [b1] to [b2]
+  dB = hypot(xo      -xv[b1], yo      -yv[b1]);
+  dP = hypot(xv[peak]-xo,     yv[peak] -yo);
+             
+  // XXX we could just choose the point based on s1 vs s2...
+  s1 = hypot(xv[peak]-xv[b1], yv[peak]-yv[b1]);
+  s2 = hypot(xv[peak]-xv[b2], yv[peak]-yv[b2]);
+
+  r1 = (SQ(s1) - SQ(dB) - SQ(dP)) / (2*dP);
+  r2 = (SQ(s2) - SQ(dB) - SQ(dP)) / (2*dP);
+
+  // dr >= 0
+  dr = MAX (r1, r2);
+
+  dx = -parity*dr*sin(angle);
+  dy = -parity*dr*cos(angle);
+
+  xo += dx;
+  yo += dy;
 
   // xc, yc is the true image center
@@ -163 +228 @@
   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));