Changeset 33481

Timestamp:

Mar 11, 2012, 10:56:40 PM (14 years ago)

Author:

eugene

Message:

speed up select_images and initMosaic

Location:

branches/eam_branches/ipp-20111122/Ohana/src/relphot/src

Files:

• MosaicOps.c (modified) (2 diffs)
• select_images.c (modified) (5 diffs)

branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/MosaicOps.c

@@ -36 +36 @@
 
 /* find mosaic frames (unique time periods & photcode name matches mosaic) */
-void initMosaics (Image *image, off_t Nimage) {
+void initMosaics_old (Image *image, off_t Nimage) {
 
   off_t i, j, status, found, NMOSAIC, *MosaicN_IMAGE;
@@ -129 +129 @@
 
   initMosaicGrid (image, Nimage);
+
+  fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nimage, (int) Nmosaic);
   return;
+}
+
+off_t findMosaic (unsigned int *startTimes, off_t Nmosaic, unsigned int start);
+
+void sort_times (unsigned int *T, int N) {
+
+# define SWAPFUNC(A,B){ unsigned int tmp; \
+  tmp = T[A]; T[A] = T[B]; T[B] = tmp; \
+}
+# define COMPARE(A,B)(T[A] < T[B])
+
+  OHANA_SORT (N, COMPARE, SWAPFUNC);
+
+# undef SWAPFUNC
+# undef COMPARE
+
+}
+
+/* find mosaic frames (unique time periods & photcode name matches mosaic) */
+void initMosaics (Image *image, off_t Nimage) {
+
+  off_t i, j, status, found, NMOSAIC, *MosaicN_IMAGE;
+  unsigned int start, stop, *startTimes, *startTimesMosaic;
+  char *pname;
+
+  if (!MOSAIC_ZEROPT) return;
+
+  /* a 'mosaic' in relphot is (unlike relastro) a virtual concept: there is no 
+   * entry in the image table that represents this mosaic.  Instead, it is an
+   * internal construct that defines a group of related images 
+   */
+
+  // generate a list of unique start times (these define the mosaics)
+  ALLOCATE (startTimes, unsigned int, Nimage);
+  for (i = 0; i < Nimage; i++) {
+    startTimes[i] = image[i].tzero;
+  }
+  sort_times (startTimes, Nimage);
+  
+  Nmosaic = 0;
+  NMOSAIC = 1000;
+  ALLOCATE (startTimesMosaic, unsigned int, NMOSAIC);
+  startTimesMosaic[0] = startTimes[0];
+
+  for (i = 0; i < Nimage; i++) {
+    if (startTimes[i] < startTimesMosaic[Nmosaic]) {
+      fprintf (stderr, "error?\n");
+      abort();
+    }
+    if (startTimes[i] == startTimesMosaic[Nmosaic]) continue;
+    Nmosaic ++;
+    if (Nmosaic >= NMOSAIC) {
+      NMOSAIC += 1000;
+      REALLOCATE (startTimesMosaic, unsigned int, NMOSAIC);
+    }
+    startTimesMosaic[Nmosaic] = startTimes[i];
+  }
+  Nmosaic ++;
+
+  // now I have a list of uniq start times, and they are in order
+  // create the mosaic arrays for these times
+  ALLOCATE (mosaic, Mosaic, Nmosaic);
+
+  ALLOCATE (MosaicToImage, off_t *, Nmosaic);
+  ALLOCATE (MosaicN_Image, off_t,   Nmosaic);
+  ALLOCATE (MosaicN_IMAGE, off_t,   Nmosaic);
+
+  for (i = 0; i < Nmosaic; i++) {
+    /* a new mosaic, define ranges */
+    mosaic[i].start = startTimesMosaic[i];
+    mosaic[i].stop  = 0;
+    mosaic[i].Mcal  = 0.0;
+    mosaic[i].dMcal = 0.0;
+    mosaic[i].dMsys = 0.0;
+    mosaic[i].Xm    = 0.0;
+    mosaic[i].flags = 0;
+    mosaic[i].secz  = NAN;
+    mosaic[i].photcode = 0;
+
+    MosaicN_IMAGE[i] = 10;
+    MosaicN_Image[i] = 0;
+    ALLOCATE (MosaicToImage[i], off_t, MosaicN_IMAGE[i]);
+    MosaicToImage[i][0] = -1;
+  }
+
+  ALLOCATE (ImageToMosaic, off_t, Nimage); // mosaic to which image belongs
+
+  // assign each image to a mosaic
+  for (i = 0; i < Nimage; i++) {
+    ImageToMosaic[i] = -1;
+
+    /* select valid mosaic images by photcode */
+    pname = GetPhotcodeNamebyCode (image[i].photcode);
+    status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME));
+    if (status) continue;
+
+    /* set image time range */
+    // start = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1);
+    // stop  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
+
+    start = image[i].tzero;
+    stop  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
+
+    j = findMosaic(startTimesMosaic, Nmosaic, start);
+    if (j == -1) {
+      fprintf (stderr, "error?\n");
+      abort();
+    }
+
+    // add reference from image to mosaic
+    ImageToMosaic[i] = j;
+
+    // have we already found this mosaic?
+    found = (MosaicN_Image[j] > 0);
+
+    /* add image to mosaic image list */
+    MosaicToImage[j][MosaicN_Image[j]] = i;
+    MosaicN_Image[j] ++;
+    if (MosaicN_Image[j] == MosaicN_IMAGE[j]) {
+      MosaicN_IMAGE[j] += 10;
+      REALLOCATE (MosaicToImage[j], off_t, MosaicN_IMAGE[j]);
+    }
+    if (found) continue;
+    
+    /* a new mosaic, define ranges */
+    if (mosaic[j].start != start) { 
+      fprintf (stderr, "error?\n");
+      abort();
+    }
+    mosaic[j].stop  = stop;
+    mosaic[j].Mcal  = 0.0;
+    mosaic[j].dMcal = 0.0;
+    mosaic[j].dMsys = 0.0;
+    mosaic[j].Xm    = 0.0;
+    mosaic[j].flags  = image[i].flags;
+    mosaic[j].secz  = image[i].secz;
+    mosaic[j].photcode = GetPhotcodeEquivCodebyCode (image[i].photcode);
+  }
+
+  // free this or not?
+  free (MosaicN_IMAGE);
+  free (startTimes);
+  free (startTimesMosaic);
+
+  initMosaicGrid (image, Nimage);
+
+  fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nimage, (int) Nmosaic);
+  return;
+}
+
+// use bisection to find the overlapping mosaic
+off_t findMosaic (unsigned int *startTimes, off_t Nmosaic, unsigned int start) {
+
+  off_t Nlo, Nhi, N;
+
+  // find the last mosaic before start
+  Nlo = 0; // startTimes[Nlo] guaranteed to be <= start
+  Nhi = Nmosaic - 1; // startTimes[Nhi] guaranteed to be >= start
+  while (Nhi - Nlo > 10) {
+    N = 0.5*(Nlo + Nhi);
+    if (startTimes[N] < start) {
+      Nlo = MAX(N, 0);
+    } else {
+      Nhi = MIN(N, Nmosaic - 1);
+    }
+  }
+
+  // check for the matched mosaic starting from Nlo 
+  for (N = Nlo; N < Nmosaic; N++) { 
+    if (start > startTimes[N])  continue;
+    return (N);
+  }
+  return (-1);
 }
 

branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/select_images.c

@@ -25 +25 @@
   
   Image *image;
-  off_t i, j, k, m, nStart, iSky, nimage, NIMAGE;
+  off_t i, j, k, m, nStart, iSky, nimage, NIMAGE, D_NIMAGE;
   off_t *line_number;
   int InRange, ecode, found;
@@ -63 +63 @@
   DminSkyRegion = +90.0;
   DmaxSkyRegion = -90.0;
+
+  D_NIMAGE = 6000;
 
   // FILE *ftest = fopen ("relphot.dump.dat", "w");
@@ -110 +112 @@
 
   nimage = 0;
-  NIMAGE = 100;
+  NIMAGE = D_NIMAGE;
   ALLOCATE (image, Image, NIMAGE);
   ALLOCATE (line_number, off_t, NIMAGE);
@@ -117 +119 @@
   for (i = 0; i < Ntimage; i++) {
       
+    if (!(i % 300000)) fprintf (stderr, ".");
+
     /* exclude images by photcode */
     ecode = GetPhotcodeEquivCodebyCode (timage[i].photcode);
@@ -248 +252 @@
     nimage ++;
     if (nimage == NIMAGE) {
-      NIMAGE += 100;
+      NIMAGE += D_NIMAGE;
+      D_NIMAGE = MAX (100000, D_NIMAGE * 1.5);
       REALLOCATE (image, Image, NIMAGE);
       REALLOCATE (line_number, off_t, NIMAGE);