Changeset 10505

Timestamp:

Dec 6, 2006, 12:43:52 PM (20 years ago)

Author:

eugene

Message:

added sane plotting for very large images

File:

• trunk/Ohana/src/opihi/dvo/images.c (modified) (7 diffs)

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

@@ -1 +1 @@
 # include "dvoshell.h"
+# define BETA 0.41421
 
 int images (int argc, char **argv) {
@@ -5 +6 @@
   int i, j, Nimage, status, InPic, leftside, *plist, TimeSelect, ByName;
   time_t tzero, tend;
-  int N, NPTS, n, npts;
+  int N, NPTS, n, npts, Npts;
   Vector Xvec, Yvec;
-  double r, d, x[4], y[4], Rmin, Rmax, Rmid, trange;
+  double r[8], d[8], x[8], y[8], Rmin, Rmax, Rmid, trange, Radius;
   Image *image;
   Graphdata graphmode;
@@ -20 +21 @@
     remove_argument (N, &argc, argv);
     ByName = TRUE;
+  }
+
+  Radius = 45;
+  if ((N = get_argument (argc, argv, "-radius"))) {
+    remove_argument (N, &argc, argv);
+    Radius = atof (argv[N]);
+    remove_argument (N, &argc, argv);
   }
 
@@ -97 +105 @@
       x[3] = 0;           y[3] = image[i].NY;
     }
-    status = FALSE;
-
+
+    Npts = 4;
+    status = TRUE;
     leftside = FALSE;
     for (j = 0; j < 4; j++) {
-      XY_to_RD (&r, &d, x[j], y[j], &image[i].coords);
-      while (r < Rmin) { r += 360.0; }
-      while (r > Rmax) { r -= 360.0; }
+      status = XY_to_RD (&r[j], &d[j], x[j], y[j], &image[i].coords);
+      if (!status) break;
+      while (r[j] < Rmin) { r[j] += 360.0; }
+      while (r[j] > Rmax) { r[j] -= 360.0; }
       if (j == 0) {
-        leftside = (r < Rmid);
+        leftside = (r[j] < Rmid);
       } 
       if (j > 0) { 
-        if ( leftside && (r > Rmid + 90)) { r -= 360.0; }
-        if (!leftside && (r < Rmid - 90)) { r += 360.0; }
-      }
-      status |= fRD_to_XY (&Xvec.elements[N+2*j], &Yvec.elements[N+2*j], r, d, &graphmode.coords);
+        if ( leftside && (r[j] > Rmid + 90)) { r[j] -= 360.0; }
+        if (!leftside && (r[j] < Rmid - 90)) { r[j] += 360.0; }
+      }
+    }
+
+    // extremely large-scale images with certain projections will have odd boundaries.
+    // eg, the ASCA images are essentially circles of radius ~60 degrees.  plot these as 
+    // octagons with some dummy size.
+    if (!status) {
+      int jp, xo, yo;
+      double rc, dc;
+      xo = 0.5*image[i].NX;
+      yo = 0.5*image[i].NY;
+      // is the image center on the screen?
+      status = XY_to_RD (&rc, &dc, 0.5*x[2], 0.5*y[2], &image[i].coords);
+      if (status && (image[i].NX * image[i].coords.cdelt2 > 90)) {
+        // draw an octagon with radius 45 degrees
+        double dX = Radius / image[i].coords.cdelt2;
+        x[0] = xo+dX;      y[0] = yo+BETA*dX;
+        x[1] = xo+BETA*dX; y[1] = yo+dX;
+        x[2] = xo-BETA*dX; y[2] = yo+dX;
+        x[3] = xo-dX;      y[3] = yo+BETA*dX;
+        x[4] = xo-dX;      y[4] = yo-BETA*dX;
+        x[5] = xo-BETA*dX; y[5] = yo-dX;
+        x[6] = xo+BETA*dX; y[6] = yo-dX;
+        x[7] = xo+dX;      y[7] = yo-BETA*dX;
+        Npts = 8;
+        j = 0;
+        for (jp = 0; jp < 8; jp++) {
+          status = XY_to_RD (&r[j], &d[j], x[jp], y[jp], &image[i].coords);
+          if (!status) continue;
+          while (r[j] < Rmin) { r[j] += 360.0; }
+          while (r[j] > Rmax) { r[j] -= 360.0; }
+          if (j == 0) {
+            leftside = (r[j] < Rmid);
+          } 
+          if (j > 0) { 
+            if ( leftside && (r[j] > Rmid + 90)) { r[j] -= 360.0; }
+            if (!leftside && (r[j] < Rmid - 90)) { r[j] += 360.0; }
+          }
+          j++;
+        }
+        Npts = j;
+      } else {
+        continue;
+      }
+    }
+    if (Npts == 0) continue;
+
+    // if any of the points is on the screen, plot the image
+    status = FALSE;
+    for (j = 0; j < Npts; j++) {
+      status |= fRD_to_XY (&Xvec.elements[N+2*j], &Yvec.elements[N+2*j], r[j], d[j], &graphmode.coords);
       if (j > 0) {
         Xvec.elements[N+2*j - 1] = Xvec.elements[N+2*j];
@@ -117 +176 @@
       }
     }
-    Xvec.elements[N+7] = Xvec.elements[N];
-    Yvec.elements[N+7] = Yvec.elements[N];
+    Xvec.elements[N+2*Npts-1] = Xvec.elements[N];
+    Yvec.elements[N+2*Npts-1] = Yvec.elements[N];
     InPic = FALSE;
-    for (j = 0; j < 8; j+=2) {
+    for (j = 0; j < 2*Npts; j+=2) {
       if ((Xvec.elements[N+j] >= graphmode.xmin) && 
           (Xvec.elements[N+j] <= graphmode.xmax) && 
@@ -134 +193 @@
         REALLOCATE (plist, int, npts);
       }
-      N+=8;
-      if (N > NPTS - 1) {  /* this is OK because NPTS is made always a multiple of 8 */
-        NPTS += 200;
+      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);
@@ -143 +202 @@
   }
 
-  gprint (GP_ERR, "plotting %d images\n", N/8);
+  gprint (GP_ERR, "plotting %d images\n", n);
   Xvec.Nelements = Xvec.Nelements = N;
   if (N > 0) {