Changeset 30443 for branches/eam_branches/ipp-20101205/Ohana/src/relphot/src/select_images.c

Timestamp:

Jan 31, 2011, 2:47:28 PM (15 years ago)

Author:

eugene

Message:

speed up mosaic match to detections by limiting the search region and sorting

File:

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

branches/eam_branches/ipp-20101205/Ohana/src/relphot/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;
@@ -35 +44 @@
     return NULL;
   }
+
+  gettimeofday (&start, (void *) NULL);
 
   // the comparison is made in the catalog local projection. below we set crval1,2
@@ -47 +58 @@
   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 */
@@ -75 +91 @@
     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;
@@ -100 +124 @@
     }
     
+    // this adds 1.3 sec for 3M images
     if (!FindMosaicForImage (timage, Ntimage, i)) {
       fprintf (stderr, "cannot find mosaic for "OFF_T_FMT"\n", i);
@@ -113 +138 @@
     found = FALSE;
 
-    /* transform corners to ra,dec */
+    /* transform corners to ra,dec -- costs ~3sec for 3M images */
     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]);
-    }
-
-    // RA(nStart) is guaranteed to be < RminImage:
+      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: -- costs 0.5sec for 3M images
     nStart = getRegionStartByRA (RminImage, RmaxSky, skylist[0].Nregions);
 
     /* compare with each region file */
-    for (iSky = 0; (iSky < skylist[0].Nregions) && !found; iSky++) { 
+    for (iSky = nStart; (iSky < skylist[0].Nregions) && !found; iSky++) { 
 
       m = index[iSky];
@@ -174 +218 @@
     }
   }
-      
+  MARKTIME("finish image selection: %f sec\n", dtime);
+
   if (VERBOSE) fprintf (stderr, "found "OFF_T_FMT" images\n", nimage);
 
@@ -295 +340 @@
   return (Nlo);
 }
+
+off_t getRegionStopByRA (double R, double *Rref, off_t Nregions) {
+
+  // use bisection to find the overlapping mosaic
+
+  off_t Nlo, Nhi, N;
+
+  // find the last mosaic before start
+  Nlo = 0; Nhi = Nregions;
+  while (Nhi - Nlo > 10) {
+    N = 0.5*(Nlo + Nhi);
+    if (Rref[N] < R) {
+      Nlo = MAX(N, 0);
+    } else {
+      Nhi = MIN(N, Nregions);
+    }
+  }
+  return (Nlo);
+}