Changeset 12753

Timestamp:

Apr 5, 2007, 1:36:09 PM (19 years ago)

Author:

eugene

Message:

use TRI images to represent skycell triangles

Location:

trunk/Ohana/src

Files:

• addstar/src/ConfigInit_skycells.c (modified) (1 diff)
• addstar/src/args_skycells.c (modified) (1 diff)
• addstar/src/sky_tessalation.c (modified) (3 diffs)
• addstar/src/skycells.c (modified) (2 diffs)
• opihi/dvo/gimages.c (modified) (1 diff)
• opihi/dvo/images.c (modified) (5 diffs)

trunk/Ohana/src/addstar/src/ConfigInit_skycells.c

@@ -1 @@
-# include "addstar.h"
+# include "skycells.h"
 
 void GetConfig (char *config, char *field, char *format, int N, void *ptr);

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

@@ -1 @@
-# include "addstar.h"
+# include "skycells.h"
 static void help (void);
 

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

@@ -23 +23 @@
   for (i = 0; i < Ntriangles; i++) {
     sky_triangle_to_image (&image[i], &tri[i]);
+    sprintf (image[i].name, "test.%03d", i);
   }  
 
@@ -35 +36 @@
   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);
 
-  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)
+  // we will project to the triangle center position
+  refcoords[0].crval1 = triangle[0].r;
+  refcoords[0].crval2 = triangle[0].d;
+
+  // project the vertices to this projection, find bounds
+  Xmin = Xmax = Ymin = Ymax = 0.0; // 0,0 is center of triangle
+  for (i = 0; i < 3; i++) {
+    RD_to_XY (&xv[i], &yv[i], triangle[0].rv[i], triangle[0].dv[i], refcoords);
+    Xmin = MIN (xv[i], Xmin);
+    Xmax = MAX (xv[i], Xmax);
+    Ymin = MIN (yv[i], Ymin);
+    Ymax = MAX (yv[i], Ymax);
+  }
+
+  // set NX, NY to the roughly full-width box (centered at 0,0)
+  NX = Xmax - Xmin;
+  NY = Ymax - Ymin;
 
   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);
+  image[0].coords.pc1_1 = image[0].coords.pc2_2 = 1.0;
+  image[0].coords.pc1_2 = image[0].coords.pc2_1 = 0.0;
+
+  // We cannot use the correction below if we want to set cdelt1,2 to our desired pixel scale
+  // use this test to raise an error (60000 x 60000 is a very large image...)
+  strcpy (image[0].coords.ctype, "TRI--TAN");
+  scale = 0;
+  for (i = 0; i < 3; i++) {
+    scale = MAX (abs(xv[i]), scale);
+    scale = MAX (abs(yv[i]), scale);
+  }
+  if (scale > 32000) {
+    scale /= 30000.0;
     NX /= scale;
     NY /= scale;
     image[0].coords.cdelt1 *= scale;
     image[0].coords.cdelt2 *= scale;
+    for (i = 0; i < 3; i++) {
+      xv[i] /= scale;
+      yv[i] /= scale;
+    }
   }
   image[0].NX = NX;
   image[0].NY = NY;
-  image[0].code = 1;
-
-  // fprintf (stderr, "%d  %f  %d %d %d  %d %d\n", parity, angle*DEG_RAD, peak, b1, b2, NX, NY);
+
+  image[0].code = 1; // this needs to be set more sensibly
+
+  image[0].Mx   = xv[0];  image[0].My   = yv[0];
+  image[0].Mxxx = xv[1];  image[0].Mxyy = yv[1];
+  image[0].Mxxy = xv[2];  image[0].Myyy = yv[2];
 
   return (TRUE);
@@ -298 +282 @@
 }
 
+
+
+/*** 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

@@ -1 @@
-# include "addstar.h"
+# include "skycells.h"
 
 int main (int argc, char **argv) {
@@ -29 +29 @@
   if (db.dbstate == LCK_EMPTY) {
     if (VERBOSE) fprintf (stderr, "can't find %s, creating a new one\n", ImageCat);
-    dvo_image_create (&db, ZeroPt);
+    dvo_image_create (&db, 25.0);
   } else {
-    if (!dvo_image_load (&db, VERBOSE, FORCE_READ)) {
+    if (!dvo_image_load (&db, VERBOSE, FALSE)) {
       Shutdown ("can't read image catalog %s", db.filename);
     }

trunk/Ohana/src/opihi/dvo/gimages.c

@@ -80 +80 @@
   int TriangleUp   = wordhash ("TRP-");
   int TriangleDn   = wordhash ("TRM-");
+  int TrianglePts  = wordhash ("TRI-");
 
   Nfound = 0;

trunk/Ohana/src/opihi/dvo/images.c

@@ -114 +114 @@
   int TriangleUp   = wordhash ("TRP-");
   int TriangleDn   = wordhash ("TRM-");
+  int TrianglePts  = wordhash ("TRI-");
 
   npts = NPTS = 200;
@@ -155 +156 @@
       x[1] = +0.5*image[i].NX; y[1] = +0.5*image[i].NY;
       x[2] = -0.5*image[i].NX; y[2] = +0.5*image[i].NY;
+      goto got_type;
+    }
+    // For 'TrianglePts' (TRI-), we are using the Mx,My, etc terms to save the coordinates
+    // this means triangular images cannot carry photometric zero-point variations
+    if (typehash == TrianglePts) {
+      Npts = 3;
+      x[0] = image[i].Mx;   y[0] = image[i].My;
+      x[1] = image[i].Mxxx; y[1] = image[i].Mxyy;
+      x[2] = image[i].Mxxy; y[2] = image[i].Myyy;
       goto got_type;
     }
@@ -228 +238 @@
     if (Npts == 0) continue;
 
-    // if any of the points is on the screen, plot the image
     status = FALSE;
     for (j = 0; j < Npts; j++) {
@@ -239 +248 @@
     Xvec.elements[N+2*Npts-1] = Xvec.elements[N];
     Yvec.elements[N+2*Npts-1] = Yvec.elements[N];
+    if (!status) continue;
+    // if none of the points are on the visible side of the projection, do not plot the image
+
     InPic = FALSE;
     for (j = 0; j < 2*Npts; j+=2) {
@@ -247 +259 @@
         InPic = TRUE;
     }
-    if (!status && HIDDEN) status = TRUE;
-    if (InPic && status) {
-      plist[n] = i;
-      n++;
-      if (n > npts - 1) {
-        npts += 200;
-        REALLOCATE (plist, int, npts);
-      }
-      N+=2*Npts;
-      if (N + 16 >= NPTS) {  /* need to leave room for 8 point image */
-        NPTS += 400;
-        REALLOCATE (Xvec.elements, float, NPTS);
-        REALLOCATE (Yvec.elements, float, NPTS);
-      }
+    // if (!status && HIDDEN) status = TRUE;
+    if (!InPic) continue;
+
+    plist[n] = i;
+    n++;
+    if (n > npts - 1) {
+      npts += 200;
+      REALLOCATE (plist, int, npts);
+    }
+    N+=2*Npts;
+    if (N + 16 >= NPTS) {  /* need to leave room for 8 point image */
+      NPTS += 400;
+      REALLOCATE (Xvec.elements, float, NPTS);
+      REALLOCATE (Yvec.elements, float, NPTS);
     }
   }