Changeset 27839 for branches/simmosaic_branches/Ohana/src/relphot/src/select_images.c

Timestamp:

May 3, 2010, 8:45:22 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simmosaic_branches

Files:

• . (modified) (1 prop)
• Ohana (modified) (1 prop)
• Ohana/src/relphot/src/select_images.c (modified) (9 diffs)

branches/simmosaic_branches/Ohana/src/relphot/src/select_images.c

@@ -13 +13 @@
 } SkyRegionCoords;
 
-Image *select_images (SkyList *skylist, Image *timage, int Ntimage, int **LineNumber, int *Nimage) {
+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;
-  int i, j, k, m, found, nimage, NIMAGE;
-  int InRange, ecode;
+  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;
-  int *line_number;
   Coords tcoords;
   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);
@@ -71 +85 @@
   NIMAGE = 100;
   ALLOCATE (image, Image, NIMAGE);
-  ALLOCATE (line_number, int, NIMAGE);
+  ALLOCATE (line_number, off_t, NIMAGE);
   
   // go through the complete list of images, selecting ones which overlap any region
@@ -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 */
@@ -146 +170 @@
         NIMAGE += 100;
         REALLOCATE (image, Image, NIMAGE);
-        REALLOCATE (line_number, int, NIMAGE);
+        REALLOCATE (line_number, off_t, NIMAGE);
       }
     }
   }
       
-  if (VERBOSE) fprintf (stderr, "found %d images\n", nimage);
+  if (VERBOSE) fprintf (stderr, "found %lld images\n", (long long) nimage);
 
   REALLOCATE (image, Image, MAX (nimage, 1));
-  REALLOCATE (line_number, int, MAX (nimage, 1));
+  REALLOCATE (line_number, off_t, MAX (nimage, 1));
   free (skycoords);
 
@@ -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);
+}