Changeset 27839 for branches/simmosaic_branches/Ohana/src/relastro/src/MosaicOps.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/relastro (modified) (1 prop)
• Ohana/src/relastro/src/MosaicOps.c (modified) (10 diffs)

branches/simmosaic_branches/Ohana/src/relastro/src/MosaicOps.c

@@ -2 +2 @@
 
 // array of mosaic definition structures
-static int    Nmosaic;
+static off_t   Nmosaic;
 static Mosaic *mosaic;
 
 // list of all images associated with a mosaic
-static int   *Nmosaic_own_images; // number of images for this mosaic
-static int   *Amosaic_own_images; // size of allocated array
-static int   **mosaic_own_images; // array of arrays: mosaic -> images
+static off_t   *Nmosaic_own_images; // number of images for this mosaic
+static off_t   *Amosaic_own_images; // size of allocated array
+static off_t   **mosaic_own_images; // array of arrays: mosaic -> images
 
 // list of mosaic associated with each image  
-static int    Nmosaic_for_images; // number of images (for internal checks)
-static int    *mosaic_for_images; // array of: image -> mosaic
-
-Mosaic *getmosaics (int *N) {
+static off_t    Nmosaic_for_images; // number of images (for off_ternal checks)
+static off_t    *mosaic_for_images; // array of: image -> mosaic
+
+Mosaic *getmosaics (off_t *N) {
   *N = Nmosaic;
   return (mosaic);
 }
 
+off_t getMosaicByTimes (unsigned int start, unsigned int stop, unsigned int *startMos, unsigned int *stopMos, off_t *indexMos) {
+
+  // use bisection to find the overlapping mosaic
+
+  off_t Nlo, Nhi, N;
+
+  // find the last mosaic before start
+  Nlo = 0; Nhi = Nmosaic;
+  while (Nhi - Nlo > 10) {
+    N = 0.5*(Nlo + Nhi);
+    if (startMos[N] < start) {
+      Nlo = MAX(N, 0);
+    } else {
+      Nhi = MIN(N + 1, Nmosaic);
+    }
+  }
+
+  // check for the matched mosaic starting from Nlo 
+  // we may have to go much beyond Nlo since stop is not sorted
+  // can we use a sorted version of stop to check when we are beyond the valid range??
+  for (N = Nlo; N < Nmosaic; N++) { 
+    if (stop  < stopMos[N]) continue;
+    if (start > startMos[N])  continue;
+    if (stop  < startMos[N]) return (-1);
+    return (indexMos[N]);
+  }
+
+  return (-1);
+}
+
+// sort two times vectors and an index by first time vector
+void sort_mosaic_times (unsigned int *S, unsigned int *E, off_t *I, off_t N) {
+
+# define SWAPFUNC(A,B){ unsigned int tmp_t; off_t tmp_i; \
+  tmp_t = S[A]; S[A] = S[B]; S[B] = tmp_t; \
+  tmp_t = E[A]; E[A] = E[B]; E[B] = tmp_t; \
+  tmp_i = I[A]; I[A] = I[B]; I[B] = tmp_i; \
+}
+# define COMPARE(A,B)(S[A] < S[B])
+
+  OHANA_SORT (N, COMPARE, SWAPFUNC);
+
+# undef SWAPFUNC
+# undef COMPARE
+
+}
+
+// first, let's continue to use the time to make the match, but use bracketing to make it faster:
+
 // find mosaic frames (unique time periods & photcode name matches mosaic) 
-void initMosaics (Image *image, int Nimage) {
-
-  int i, j, found, NMOSAIC;
+void initMosaics (Image *image, off_t Nimage) {
+
+  off_t i, Nmos, NMOSAIC;
   unsigned int start, stop;
+
+  unsigned int *startMos, *stopMos;
+  off_t *indexMos;
 
   Nmosaic = 0;
@@ -29 +81 @@
   ALLOCATE (mosaic, Mosaic, NMOSAIC);
 
-  ALLOCATE (Nmosaic_own_images, int, NMOSAIC);
-  ALLOCATE (Amosaic_own_images, int, NMOSAIC);
-  ALLOCATE (mosaic_own_images, int *, NMOSAIC);
+  ALLOCATE (Nmosaic_own_images, off_t, NMOSAIC);
+  ALLOCATE (Amosaic_own_images, off_t, NMOSAIC);
+  ALLOCATE (mosaic_own_images, off_t *, NMOSAIC);
+  ALLOCATE (startMos, unsigned int, NMOSAIC);
+  ALLOCATE (stopMos, unsigned int, NMOSAIC);
+  ALLOCATE (indexMos, off_t, NMOSAIC);
 
   /* find the mosaic images (coords.ctype = DIS); generate list of unique mosaics */
@@ -54 +109 @@
     Nmosaic_own_images[Nmosaic] = 0;
     Amosaic_own_images[Nmosaic] = 10;
-    ALLOCATE (mosaic_own_images[Nmosaic], int, Amosaic_own_images[Nmosaic]);
+    ALLOCATE (mosaic_own_images[Nmosaic], off_t, Amosaic_own_images[Nmosaic]);
+
+    startMos[Nmosaic] = start;
+    stopMos[Nmosaic] = stop;
+    indexMos[Nmosaic] = Nmosaic;
 
     Nmosaic ++;
@@ -60 +119 @@
       NMOSAIC += 10;
       REALLOCATE (mosaic, Mosaic, NMOSAIC);
-      REALLOCATE (mosaic_own_images, int *, NMOSAIC);
-      REALLOCATE (Nmosaic_own_images, int, NMOSAIC);
-      REALLOCATE (Amosaic_own_images, int, NMOSAIC);
-    }
-  }
-
+      REALLOCATE (mosaic_own_images, off_t *, NMOSAIC);
+      REALLOCATE (Nmosaic_own_images, off_t, NMOSAIC);
+      REALLOCATE (Amosaic_own_images, off_t, NMOSAIC);
+      REALLOCATE (startMos, unsigned int, NMOSAIC);
+      REALLOCATE (stopMos, unsigned int, NMOSAIC);
+      REALLOCATE (indexMos, off_t, NMOSAIC);
+    }
+  }
+
+  // sort the index, start, and stop by the start times:
+  sort_mosaic_times (startMos, stopMos, indexMos, Nmosaic);
+  
   // array to store image->mosaic index
   Nmosaic_for_images = Nimage;
-  ALLOCATE (mosaic_for_images, int, Nmosaic_for_images);
+  ALLOCATE (mosaic_for_images, off_t, Nmosaic_for_images);
 
   /* now assign the WRP images to these mosaics */
@@ -80 +145 @@
     stop  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
 
-    /* find existing mosaic with this time range */
-    found = FALSE;
-    for (j = 0; !found && (j < Nmosaic); j++) { 
-      if (stop  < mosaic[j].start) continue;
-      if (start > mosaic[j].stop)  continue;
-      found = TRUE;
-      break;
-    }
-    /* if no matching mosaic exists, skip this image */
-    if (!found) continue;
+    Nmos = getMosaicByTimes (start, stop, startMos, stopMos, indexMos);
+    if (Nmos == -1) {
+      fprintf (stderr, "cannot match mosaic for %s\n", image[i].name);
+      continue;
+    }
 
     // mosaic corresponding to this image
-    mosaic_for_images[i] = j;
+    mosaic_for_images[i] = Nmos;
 
     // add image to mosaic_own_image list 
-    mosaic_own_images[j][Nmosaic_own_images[j]] = i;
-    Nmosaic_own_images[j] ++;
-    if (Nmosaic_own_images[j] == Amosaic_own_images[j]) {
-      Amosaic_own_images[j] += 10;
-      REALLOCATE (mosaic_own_images[j], int, Amosaic_own_images[j]);
-    }
-    assert (Nmosaic_own_images[j] < Amosaic_own_images[j]);
-  }
-
+    mosaic_own_images[Nmos][Nmosaic_own_images[Nmos]] = i;
+    Nmosaic_own_images[Nmos] ++;
+    if (Nmosaic_own_images[Nmos] == Amosaic_own_images[Nmos]) {
+      Amosaic_own_images[Nmos] += 10;
+      REALLOCATE (mosaic_own_images[Nmos], off_t, Amosaic_own_images[Nmos]);
+    }
+    assert (Nmosaic_own_images[Nmos] < Amosaic_own_images[Nmos]);
+  }
+
+  free (startMos);
+  free (stopMos);
+  free (indexMos);
   return;
 }
@@ -109 +172 @@
 // return StarData values for detections in the specified image, converting coordinates from the
 // chip positions: X,Y -> L,M -> P,Q -> R,D
-StarData *getMosaicRaw (Catalog *catalog, int Ncatalog, int mos, int *Nstars) {
-
-  int i, j, im, Nraw, Nnew;
+StarData *getMosaicRaw (Catalog *catalog, int Ncatalog, off_t mos, off_t *Nstars) {
+
+  off_t i, j, im, Nraw, Nnew;
   StarData *raw, *new;
 
@@ -125 +188 @@
     // this function does the reverse-lookup for the mosaic corresponding to this image
     new = getImageRaw (catalog, Ncatalog, im, &Nnew, MODE_MOSAIC);
+    if (!new) {
+      fprintf (stderr, "inconsistent: missing mosaic for image already associated with a mosaic? (1)\n");
+      abort();
+    }
     
     // merge new and raw
@@ -142 +209 @@
 // return StarData values for averages positions in the specified image, converting coordinates from
 // the sky positions: R,D -> P,Q -> L,M -> X,Y
-StarData *getMosaicRef (Catalog *catalog, int Ncatalog, int mos, int *Nstars) {
-
-  int i, j, im, Nref, Nnew;
+StarData *getMosaicRef (Catalog *catalog, int Ncatalog, off_t mos, off_t *Nstars) {
+
+  off_t i, j, im, Nref, Nnew;
   StarData *ref, *new;
   
@@ -157 +224 @@
     // this function does the reverse-lookup for the mosaic corresponding to this image
     new = getImageRef (catalog, Ncatalog, im, &Nnew, MODE_MOSAIC);
+    if (!new) {
+      fprintf (stderr, "inconsistent: missing mosaic for image already associated with a mosaic? (2)\n");
+      abort();
+    }
     
     // merge new and ref
@@ -172 +243 @@
 }
 
-Mosaic *getMosaicForImage (int im) {
-
-  int mos;
+Mosaic *getMosaicForImage (off_t im) {
+
+  off_t mos;
 
   if (im < 0) abort();