Changeset 27480 for trunk/Ohana/src/relphot/src/select_images.c

Timestamp:

Mar 26, 2010, 1:54:17 PM (16 years ago)

Author:

eugene

Message:

add detection->image, image->mosaic optimizations from relastro

File:

• trunk/Ohana/src/relphot/src/select_images.c (modified) (8 diffs)

trunk/Ohana/src/relphot/src/select_images.c

@@ -13 +13 @@
 } SkyRegionCoords;
 
+void dsortindex (double *X, off_t *Y, int N);
+off_t getRegionStartByRA (double R, double *Rref, off_t Nregions);
+
 Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage) {
   
   Image *image;
-  off_t i, j, k, m;
-  off_t *line_number, nimage, NIMAGE;
+  off_t i, j, k, m, nStart, iSky, nimage, NIMAGE;
+  off_t *line_number;
   int InRange, ecode, found;
   double Ri[5], Di[5], Xi[5], Yi[5], dx, dy;
@@ -23 +26 @@
   SkyRegionCoords *skycoords;
   
+  double *RmaxSky;
+  off_t *index;
+
   if (skylist[0].Nregions < 1) {
     *Nimage = 0;
@@ -38 +44 @@
 
   ALLOCATE (skycoords, SkyRegionCoords, skylist[0].Nregions);
+
+  ALLOCATE (RmaxSky, double, skylist[0].Nregions);
+  ALLOCATE (index, off_t, skylist[0].Nregions);
 
   /* compare with each region file */
@@ -56 +65 @@
     skycoords[i].Yc[4] = skycoords[i].Yc[0];    
 
+    RmaxSky[i] = skylist[0].regions[i][0].Rmax;
+    index[i] = i;
+
     dx = 0.02*(skycoords[i].Xc[2] - skycoords[i].Xc[0]);
     dy = 0.02*(skycoords[i].Yc[2] - skycoords[i].Yc[0]);
@@ -65 +77 @@
   }
 
+  dsortindex (RmaxSky, index, skylist[0].Nregions);
+
   if (VERBOSE) fprintf (stderr, "finding images\n");
   BuildChipMatch (timage, Ntimage);
@@ -86 +100 @@
     }
     
-    if (!FindMosaicForImage (timage, Ntimage, i)) continue;
+    if (!FindMosaicForImage (timage, Ntimage, i)) {
+      fprintf (stderr, "cannot find mosaic for %lld\n", (long long) i);
+      continue;
+    }
 
     /* define image corners */
@@ -96 +113 @@
     found = FALSE;
 
-    /* transform to ra,dec */
+    /* transform corners to ra,dec */
+    double RminImage = 360.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:
+    nStart = getRegionStartByRA (RminImage, RmaxSky, skylist[0].Nregions);
 
     /* compare with each region file */
-    for (m = 0; (m < skylist[0].Nregions) && !found; m++) { 
+    for (iSky = 0; (iSky < skylist[0].Nregions) && !found; iSky++) { 
+
+      m = index[iSky];
 
       /* we make positional comparisons in the projection of catalog */
@@ -238 +262 @@
 }
 
+void dsortindex (double *X, off_t *Y, int N) {
+
+# define SWAPFUNC(A,B){ double tmpf; off_t tmpi; \
+  tmpf = X[A]; X[A] = X[B]; X[B] = tmpf; \
+  tmpi = Y[A]; Y[A] = Y[B]; Y[B] = tmpi; \
+}
+# define COMPARE(A,B)(X[A] < X[B])
+
+  OHANA_SORT (N, COMPARE, SWAPFUNC);
+
+# undef SWAPFUNC
+# undef COMPARE
+
+}
+
+off_t getRegionStartByRA (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);
+}