Changeset 30474

Timestamp:

Feb 3, 2011, 9:21:49 AM (15 years ago)

Author:

eugene

Message:

speed up select_images with better range checking; correctly handle bounaries of DIS images

File:

• branches/eam_branches/ipp-20101205/Ohana/src/relastro/src/select_images.c (modified) (7 diffs)

branches/eam_branches/ipp-20101205/Ohana/src/relastro/src/select_images.c

@@ -15 +15 @@
 void dsortindex (double *X, off_t *Y, int N);
 off_t getRegionStartByRA (double R, double *Rref, off_t Nregions);
+
+# define MARKTIME(MSG,...) { \
+  float dtime; \
+  gettimeofday (&stop, (void *) NULL); \
+  dtime = DTIME (stop, start); \
+  fprintf (stderr, MSG, __VA_ARGS__); }
 
 Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage) {
@@ -25 +31 @@
   Coords tcoords;
   SkyRegionCoords *skycoords;
-  
+  struct timeval start, stop;
+  
+  double RmaxSkyRegion, RminSkyRegion, DminSkyRegion, DmaxSkyRegion;
+
   double *RmaxSky;
   off_t *index;
@@ -40 +49 @@
     return NULL;
   }
+
+  gettimeofday (&start, (void *) NULL);
 
   // the comparison is made in the catalog local projection. below we set crval1,2
@@ -52 +63 @@
   ALLOCATE (RmaxSky, double, skylist[0].Nregions);
   ALLOCATE (index, off_t, skylist[0].Nregions);
+
+  RminSkyRegion = +360.0;
+  RmaxSkyRegion = -360.0;
+  DminSkyRegion = +90.0;
+  DmaxSkyRegion = -90.0;
 
   /* compare with each region file */
@@ -80 +96 @@
     skycoords[i].Xc[3] -= dx; skycoords[i].Yc[3] += dy;
     skycoords[i].Xc[4] -= dx; skycoords[i].Yc[4] -= dy;
-  }
+
+    RminSkyRegion = MIN(RminSkyRegion, skylist[0].regions[i][0].Rmin);
+    RmaxSkyRegion = MAX(RmaxSkyRegion, skylist[0].regions[i][0].Rmax);
+    DminSkyRegion = MIN(DminSkyRegion, skylist[0].regions[i][0].Dmin);
+    DmaxSkyRegion = MAX(DmaxSkyRegion, skylist[0].regions[i][0].Dmax);
+  }
+  MARKTIME("create sky region coords: %f sec\n", dtime);
 
   dsortindex (RmaxSky, index, skylist[0].Nregions);
+  MARKTIME("sort sky coords: %f sec\n", dtime);
 
   if (VERBOSE) fprintf (stderr, "finding images\n");
   BuildChipMatch (timage, Ntimage);
+  MARKTIME("build chip match: %f sec\n", dtime);
 
   nimage = 0;
@@ -124 +148 @@
     }
 
-    /* define image corners */
-    Xi[0] = 0;            Yi[0] = 0;
-    Xi[1] = timage[i].NX; Yi[1] = 0;
-    Xi[2] = timage[i].NX; Yi[2] = timage[i].NY;
-    Xi[3] = 0;            Yi[3] = timage[i].NY;
-    Xi[4] = 0;            Yi[4] = 0;
+    /* define image corners - note the DIS images (mosaic phu) are special */
+    if (!strcmp(&timage[i].coords.ctype[4], "-DIS")) {
+      Xi[0] = -0.5*timage[i].NX; Yi[0] = -0.5*timage[i].NY;
+      Xi[1] = +0.5*timage[i].NX; Yi[1] = -0.5*timage[i].NY;
+      Xi[2] = +0.5*timage[i].NX; Yi[2] = +0.5*timage[i].NY;
+      Xi[3] = -0.5*timage[i].NX; Yi[3] = +0.5*timage[i].NY;
+      Xi[4] = -0.5*timage[i].NX; Yi[4] = -0.5*timage[i].NY;
+    } else {
+      Xi[0] = 0;            Yi[0] = 0;
+      Xi[1] = timage[i].NX; Yi[1] = 0;
+      Xi[2] = timage[i].NX; Yi[2] = timage[i].NY;
+      Xi[3] = 0;            Yi[3] = timage[i].NY;
+      Xi[4] = 0;            Yi[4] = 0;
+    }
     found = FALSE;
 
     /* transform corners to ra,dec */
     double RminImage = 360.0;
+    double RmaxImage =   0.0;
+    double DminImage = +90.0;
+    double DmaxImage = -90.0;
     for (j = 0; j < 5; j++) {
       XY_to_RD (&Ri[j], &Di[j], Xi[j], Yi[j], &timage[i].coords);
       RminImage = MIN(RminImage, Ri[j]);
-    }
+      RmaxImage = MAX(RmaxImage, Ri[j]);
+      DminImage = MIN(DminImage, Di[j]);
+      DmaxImage = MAX(DmaxImage, Di[j]);
+    }
+    if (RmaxImage - RminImage > 180.0) {
+        double tmp = RminImage;
+        RmaxImage = RminImage;
+        RminImage = tmp - 360.0;
+    }
+
+    // check that this image is even in range of the searched region
+    if (DminImage > DmaxSkyRegion) continue;
+    if (DmaxImage < DminSkyRegion) continue;
+    
+    // the sky region RA is defined to be 0 - 360.0
+    if (RminImage > RmaxSkyRegion) continue;
+    if (RmaxImage < RminSkyRegion) continue;
 
     // RA(nStart) is guaranteed to be < RminImage:
@@ -194 +245 @@
     }
   }
-      
+  MARKTIME("finish image selection: %f sec\n", dtime);
+    
   if (VERBOSE) fprintf (stderr, "found "OFF_T_FMT" images\n",  nimage);