Changeset 33797 for branches/eam_branches/ipp-20120405/Ohana

Timestamp:

Apr 18, 2012, 8:16:11 AM (14 years ago)

Author:

eugene

Message:

threaded setMrel and setMmos; working on excluding mosaics with partial image selection (not finished)

Location:

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

Files:

• include/relphot.h (modified) (4 diffs)
• src/MosaicOps.c (modified) (11 diffs)
• src/StarOps.c (modified) (3 diffs)
• src/args.c (modified) (1 diff)
• src/load_images.c (modified) (6 diffs)
• src/relphot.c (modified) (1 diff)
• src/select_images.c (modified) (3 diffs)

branches/eam_branches/ipp-20120405/Ohana/src/relphot/include/relphot.h

@@ -3 +3 @@
 # include <kapa.h>
 # include <signal.h>
+# include <pthread.h>
 
 /* # define GRID_V1 */
@@ -111 +112 @@
 int    PARALLEL_MANUAL;
 int    PARALLEL_SERIAL;
+
+int    NTHREADS;
 
 int    VERBOSE;
@@ -236 +239 @@
 void          initMosaicBins      PROTO((Catalog *catalog, int Ncatalog, int doMosaicList));
 void          initMosaicGrid      PROTO((Image *image, off_t Nimage));
-void          initMosaics         PROTO((Image *image, off_t Nimage));
+void          initMosaics         PROTO((Image *subset, off_t Nsubset, Image *image, char *inSubset, off_t Nimage));
 void          initMrel            PROTO((Catalog *catalog, int Ncatalog));
 void          initialize          PROTO((int argc, char **argv));
@@ -246 +249 @@
 
 SkyList      *load_images         PROTO((FITS_DB *db, char *regionName, SkyRegion *region));
-Image        *select_images       PROTO((SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage));
+Image        *select_images       PROTO((SkyList *skylist, Image *timage, off_t Ntimage, char *inSubset, off_t **LineNumber, off_t *Nimage));
 
 int           main                PROTO((int argc, char **argv));

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

@@ -50 +50 @@
   ALLOCATE (mosaic, Mosaic, NMOSAIC);
 
-  ALLOCATE (MosaicToImage,  off_t *, NMOSAIC);
+  ALLOCATE (MosaicToImage, off_t *, NMOSAIC);
   ALLOCATE (MosaicN_Image, off_t,   NMOSAIC);
   ALLOCATE (MosaicN_IMAGE, off_t,   NMOSAIC);
@@ -95 +95 @@
     if (found) continue;
     
-    /* a new mosaic, define ranges */
+    /* a new mosaic, define ranges -- preserve the original values incase this image is not used */
     mosaic[Nmosaic].start = start;
     mosaic[Nmosaic].stop  = stop;
-    mosaic[Nmosaic].Mcal  = 0.0;
-    mosaic[Nmosaic].dMcal = 0.0;
+    mosaic[Nmosaic].Mcal  = 0.0; // note : mosaic stores only offsets relative to the original image values
+    mosaic[Nmosaic].dMcal = 0.0; // note : at the end, mosaic.Mcal is added back to the input images
     mosaic[Nmosaic].dMsys = 0.0;
     mosaic[Nmosaic].Xm    = 0.0;
@@ -152 +152 @@
 }
 
-/* find mosaic frames (unique time periods & photcode name matches mosaic) */
-void initMosaics (Image *image, off_t Nimage) {
+/* find mosaic frames (unique time periods) (NOTE : we do NOT require matching photcodes...)
+   this function will also identify the images NOT in the subset which belong to a selected mosaic
+ */
+void initMosaics (Image *subset, off_t Nsubset, Image *image, char *inSubset, off_t Nimage) {
 
   off_t i, j, status, found, NMOSAIC, *MosaicN_IMAGE;
@@ -166 +168 @@
    */
 
-  // 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);
+  // generate a list of all subset image start times
+  ALLOCATE (startTimes, unsigned int, Nsubset);
+  for (i = 0; i < Nsubset; i++) {
+    startTimes[i] = subset[i].tzero;
+  }
+  sort_times (startTimes, Nsubset);
   
   Nmosaic = 0;
@@ -178 +180 @@
   startTimesMosaic[0] = startTimes[0];
 
-  for (i = 0; i < Nimage; i++) {
+  // generate a list of the unique start times (these define the mosaics)
+  for (i = 0; i < Nsubset; i++) {
     if (startTimes[i] < startTimesMosaic[Nmosaic]) {
       fprintf (stderr, "error?\n");
@@ -192 +195 @@
   }
   Nmosaic ++;
+
+  int Nskip, Nmiss, Nmark;
+  Nskip = Nmiss = Nmark = 0;
+  // find any start times which match unselected images 
+  for (i = 0; i < Nimage; i++) {
+    if (inSubset[i]) {
+      Nskip ++;
+      continue;
+    }
+    
+    /* select valid mosaic images by photcode */
+    pname = GetPhotcodeNamebyCode (image[i].photcode);
+    if (!pname) continue;
+    status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME));
+    if (status) continue;
+
+    /* set image time range */
+    start = subset[i].tzero;
+    stop  = subset[i].tzero + MAX(1.01*subset[i].trate*subset[i].NY, 1);
+
+    /* find a matching mosaic */
+    j = findMosaic(startTimesMosaic, Nmosaic, start);
+    if (j == -1) {
+      // no matching mosaic
+      Nmiss ++;
+    } else {
+      // mark this mosaic as bad
+      Nmark ++;
+    }
+  }
+  fprintf (stderr, "%d images, %d skip, %d miss, %d mark\n", (int) Nimage, (int) Nskip, (int) Nmiss, (int) Nmark);
 
   // now I have a list of uniq start times, and they are in order
@@ -219 +253 @@
   }
 
-  ALLOCATE (ImageToMosaic, off_t, Nimage); // mosaic to which image belongs
+  ALLOCATE (ImageToMosaic, off_t, Nsubset); // mosaic to which image belongs
 
   // assign each image to a mosaic
-  for (i = 0; i < Nimage; i++) {
+  for (i = 0; i < Nsubset; i++) {
     ImageToMosaic[i] = -1;
 
     /* select valid mosaic images by photcode */
-    pname = GetPhotcodeNamebyCode (image[i].photcode);
+    pname = GetPhotcodeNamebyCode (subset[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);
+    // start = subset[i].tzero - MAX(0.01*subset[i].trate*subset[i].NY, 1);
+    // stop  = subset[i].tzero + MAX(1.01*subset[i].trate*subset[i].NY, 1);
+
+    start = subset[i].tzero;
+    stop  = subset[i].tzero + MAX(1.01*subset[i].trate*subset[i].NY, 1);
 
     j = findMosaic(startTimesMosaic, Nmosaic, start);
@@ -268 +302 @@
     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);
+    mosaic[j].flags  = subset[i].flags;
+    mosaic[j].secz  = subset[i].secz;
+    mosaic[j].photcode = GetPhotcodeEquivCodebyCode (subset[i].photcode);
   }
 
@@ -278 +312 @@
   free (startTimesMosaic);
 
-  initMosaicGrid (image, Nimage);
-
-  fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nimage, (int) Nmosaic);
+  initMosaicGrid (subset, Nsubset);
+
+  fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nsubset, (int) Nmosaic);
   return;
 }
@@ -569 +603 @@
 }
 
+typedef struct {
+  int Nfew;
+  int Nbad;
+  int Ncal;
+  int Ngrid;
+  int Nrel;
+  int Nsys;
+  off_t Nmax;
+  int PoorImages;
+  double *list;
+  double *dlist;
+  double *Mlist;
+  double *dMlist;
+} SetMmosInfo;
+
+enum {THREAD_RUN, THREAD_DONE};
+
+typedef struct {
+  int entry;
+  int state;
+  Catalog *catalog;
+  Image *image;
+  FlatCorrectionTable *flatcorr;
+  SetMmosInfo info;
+} ThreadInfo;
+
+int setMmos_mosaic (Mosaic *mosaic, off_t Nmos, Image *image, Catalog *catalog, SetMmosInfo *info, FlatCorrectionTable *flatcorr);
+void *setMmos_worker (void *data);
+int setMmos_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr);
+
+void SetMmosInfoInit (SetMmosInfo *info, off_t Nmax, int allocLists, int PoorImages) {
+  info->Nfew = 0;
+  info->Nbad = 0;
+  info->Ncal = 0;
+  info->Nrel = 0;
+  info->Ngrid = 0;
+  info->Nsys = 0;
+
+  info->Nmax = Nmax;
+  info->PoorImages = PoorImages;
+
+  if (allocLists) {
+    ALLOCATE (info->list, double, Nmax);
+    ALLOCATE (info->dlist, double, Nmax);
+    ALLOCATE (info->Mlist, double, Nmax);
+    ALLOCATE (info->dMlist, double, Nmax);
+  }
+}
+
+void SetMmosInfoFree (SetMmosInfo *info) {
+  free (info->list);
+  free (info->dlist);
+  free (info->Mlist);
+  free (info->dMlist);
+}
+
+void SetMmosInfoAccum (SetMmosInfo *summary, SetMmosInfo *results) {
+  summary->Nfew  += results->Nfew ;
+  summary->Nsys  += results->Nsys ;
+  summary->Nbad  += results->Nbad ;
+  summary->Ncal  += results->Ncal ;
+  summary->Nrel  += results->Nrel ;
+  summary->Ngrid += results->Ngrid;
+}
+
 static int npass_output = 0;
 
-int setMmos (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
+// mutex to lock setMmos_worker operations 
+static pthread_mutex_t setMmos_mutex = PTHREAD_MUTEX_INITIALIZER;
+static int nextMosaic = 0;
+
+// we have an array of mosaics (mosaic, Nmosaic).  we need to hand out mosaics one at a time to
+// the worker threads as they need
+off_t getNextMosaicForThread () {
+
+  pthread_mutex_lock (&setMmos_mutex);
+  if (nextMosaic >= Nmosaic) {
+    pthread_mutex_unlock (&setMmos_mutex);
+    return (-1);
+  }
+  int thisMosaic = nextMosaic;
+  nextMosaic ++;
+
+  pthread_mutex_unlock (&setMmos_mutex);
+  return (thisMosaic);
+}
+
+int setMmos_old (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
 
   off_t i, j, m, c, n, N, Nmax;
@@ -798 +917 @@
 }
   
+int setMmos (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
+
+  off_t i, N, Nmax;
+  Image *image;
+
+  if (!MOSAIC_ZEROPT) return (FALSE);
+  if (FREEZE_MOSAICS) return (FALSE);
+
+  if (NTHREADS) {
+    int status = setMmos_threaded (catalog, PoorImages, flatcorr);
+    return status;
+  }
+
+  image = getimages (&N, NULL);
+
+  fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE);
+
+  if (PoorImages) {
+    // XXX use bad stars and measurements for PoorImages? or not?
+    // IMAGE_BAD = STAR_BAD = MEAS_BAD = 0;
+    IMAGE_BAD = 0;
+  }
+
+  Nmax = 0;
+  for (i = 0; i < Nmosaic; i++) {
+    Nmax = MAX (Nmax, N_onMosaic[i]);
+  }
+
+  SetMmosInfo info;
+  SetMmosInfoInit (&info, Nmax, TRUE, PoorImages);
+
+  for (i = 0; i < Nmosaic; i++) {
+    setMmos_mosaic (mosaic, i, image, catalog, &info, flatcorr);
+  }
+  SetMmosInfoFree (&info);
+
+  npass_output ++;
+
+  fprintf (stderr, "%d mosaics marked having too few measurements (Nbad: %d, Ncal: %d, Ngrid: %d, Nrel: %d, Nsys: %d)\n", info.Nfew, info.Nbad, info.Ncal, info.Ngrid, info.Nrel, info.Nsys);
+
+  if (PoorImages) {
+    IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP;
+    STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
+    MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA;
+  }
+  return (TRUE);
+}
+  
+int setMmos_mosaic (Mosaic *mosaic, off_t Nmos, Image *image, Catalog *catalog, SetMmosInfo *info, FlatCorrectionTable *flatcorr) {
+
+  off_t j;
+  StatType stats;
+
+  double *list   = info->list;
+  double *dlist  = info->dlist;
+  double *Mlist  = info->Mlist;
+  double *dMlist = info->dMlist;
+
+  /* on PoorImages run, skip good images */
+  if (info->PoorImages) {
+    int bad = mosaic[0].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_SKIP);
+    if (!bad) return TRUE;
+  }      
+
+  // UBERCAL image: if this is an ubercal image, set minUbercalDist to 0:
+  // we optionally do not recalibrate images with UBERCAL zero points 
+  if (mosaic[0].flags & ID_IMAGE_PHOTOM_UBERCAL) {
+    mosaic[0].ubercalDist = 0;
+    // propagate ubercalDist to the images
+    for (j = 0; j < MosaicN_Image[Nmos]; j++) {
+      off_t im = MosaicToImage[Nmos][j];
+      image[im].ubercalDist = mosaic[0].ubercalDist;
+      // fprintf (stderr, "%d %d %d\n", (int) i, (int) im, image[im].ubercalDist);
+    }
+    if (KEEP_UBERCAL) return TRUE;
+  }
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  int minUbercalDist = 1000;
+
+  int testImage = FALSE;
+  // testImage |= (abs(mosaic[0].start - 1323003245) < 10);
+  // testImage |= (abs(mosaic[0].start - 1323003069) < 10);
+  // testImage |= (abs(mosaic[0].start - 1323003125) < 10);
+  // testImage |= (abs(mosaic[0].start - 1323003300) < 10);
+  // testImage |= (abs(mosaic[0].start - 1323003365) < 10);
+  // testImage |= (abs(mosaic[0].start - 1323003191) < 10);
+  // testImage |= (abs(mosaic[0].start - 1323003014) < 10);
+  // testImage |= (abs(mosaic[0].start - 1323003484) < 10);
+  // testImage |= (abs(mosaic[0].start - 1323003419) < 10);
+  // testImage |= (abs(mosaic[0].start - 1323002949) < 10);
+
+  FILE *fout = NULL;
+  if (testImage) {
+    char filename[64];
+    snprintf (filename, 64, "test.%05d.%02d.dat", (int) Nmos, npass_output);
+    fout = fopen (filename, "w");
+  }
+
+  // number of stars to measure the bright-end scatter
+  int Nbright = 0;
+
+  int N = 0;
+  for (j = 0; j < N_onMosaic[Nmos]; j++) {
+    float Msys, Mrel, Mcal, Mgrid, Mflat;
+      
+    off_t m = MosaicToMeasure[Nmos][j];
+    off_t c = MosaicToCatalog[Nmos][j];
+      
+    if (fout) {
+      Mcal  = getMcal  (m, c, flatcorr, catalog);
+      Mgrid = getMgrid (m, c);
+      Mrel  = getMrel  (catalog, m, c);
+      Mflat = getMflat (m, c, flatcorr, catalog);
+
+      off_t n = catalog[c].measureT[m].averef;
+      Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]);
+
+      float delta = Msys - Mrel - Mcal - Mgrid + Mflat;
+
+      int isBad = (catalog[c].measureT[m].dbFlags & MEAS_BAD);
+
+      fprintf (fout, "%f %f : %f %f %f %f %f  : %f %d\n", catalog[c].averageT[n].R, catalog[c].averageT[n].D, Msys, Mrel, Mcal, Mgrid, Mflat, delta, isBad);
+    }
+
+    if (catalog[c].measureT[m].dbFlags & MEAS_BAD) {
+      info->Nbad ++;
+      return TRUE;
+    }
+    Mcal  = getMcal  (m, c, flatcorr, catalog);
+    if (isnan(Mcal)) {
+      info->Ncal++;
+      return TRUE;
+    }
+    Mgrid = getMgrid (m, c);
+    if (isnan(Mgrid)) {
+      info->Ngrid ++;
+      return TRUE;
+    }
+    Mrel  = getMrel  (catalog, m, c);
+    if (isnan(Mrel)) {
+      info->Nrel ++;
+      return TRUE;
+    }
+      
+    // image.Mcal is not supposed to include the flat-field correction, so we need to
+    // apply that offset as well here for this image (in other words, each detection is
+    // being compared to the model, excluding the zero point, Mcal.  The model includes
+    // the flat-correction.  NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat)
+
+    Mflat = getMflat (m, c, flatcorr, catalog);
+
+    off_t n = catalog[c].measureT[m].averef;
+    Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]);
+    if (isnan(Msys)) {
+      info->Nsys++;
+      return TRUE;
+    }
+
+    PhotCode *code = GetPhotcodebyCode (catalog[c].measureT[m].photcode);
+    if (!code) goto skip;
+    if (code->equiv < 1) goto skip;
+    int Nsec = GetPhotcodeNsec (code->equiv);
+    if (Nsec == -1) goto skip;
+    minUbercalDist = MIN (catalog[c].secfilt[n*Nsecfilt + Nsec].ubercalDist, minUbercalDist);
+
+  skip:
+    list[N]  = Msys - Mrel - Mcal - Mgrid + Mflat;
+    dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
+    if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) {
+      Mlist[Nbright] = list[N];
+      dMlist[Nbright] = dlist[N];
+      Nbright ++;
+    }
+    N++;
+  }
+  /* N_onMosaic[Nmos] is all measurements, N is good measurements */
+
+  if (fout) {
+    fclose (fout);
+  }
+
+  /* too few good measurements or too many bad measurements (skip in PoorImages run) */
+  if (!info->PoorImages) {
+    int mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*N_onMosaic[Nmos]);
+    if (mark) {
+      if (VERBOSE2) { fprintf (stderr, "marked mosaic %s ("OFF_T_FMT"), (%d < %d) || (%d < %f*"OFF_T_FMT")\n", image[MosaicToImage[Nmos][0]].name,  Nmos,  N, IMAGE_TOOFEW,  N, IMAGE_GOOD_FRACTION,  N_onMosaic[Nmos]); }
+      mosaic[0].flags |= ID_IMAGE_PHOTOM_FEW;
+      info->Nfew ++;
+      if (testImage) {
+        fprintf (stderr, "NOTE: *** marked test image poor : %d %d %d***\n", (int) N, (int) IMAGE_TOOFEW, (int) (IMAGE_GOOD_FRACTION*N_onMosaic[Nmos]));
+      }
+    } else {
+      mosaic[0].flags &= ~ID_IMAGE_PHOTOM_FEW;
+    }
+  }
+  liststats (list, dlist, N, &stats);
+  if (VERBOSE2 && info->PoorImages) fprintf (stderr, "Mmos: %f %f %d %d\n", stats.mean, stats.sigma, stats.Nmeas, N);
+
+  mosaic[0].Mcal  = stats.mean;
+  mosaic[0].dMcal = stats.error;
+  mosaic[0].nFitPhotom = N;
+  mosaic[0].Xm    = 100.0*log10(stats.chisq);
+
+  if (testImage) {
+    fprintf (stderr, "test image %d (%d) %f %f %d ... ", (int) Nmos, mosaic[0].start, stats.mean, stats.error, mosaic[0].nFitPhotom);
+  }
+
+  plot_setMcal (list, N, &stats, CLOUD_TOLERANCE);
+
+  // bright end scatter
+  liststats (Mlist, dMlist, Nbright, &stats);
+  mosaic[0].dMsys = stats.sigma;
+
+  if (mosaic[0].Mcal < -CLOUD_TOLERANCE) {
+    mosaic[0].Mcal = 0.0;
+  }
+
+  if (testImage) {
+    fprintf (stderr, "%f %f\n", mosaic[0].Mcal, mosaic[0].dMsys);
+  }
+
+  // minUbercalDist calculated here is the min value for any star owned by this image
+  // since this particular image is tied to that star, bump its distance by 1
+  mosaic[0].ubercalDist = minUbercalDist + 1;
+
+  // propagate ubercalDist to the images
+  for (j = 0; j < MosaicN_Image[Nmos]; j++) {
+    off_t im = MosaicToImage[Nmos][j];
+    image[im].ubercalDist = mosaic[0].ubercalDist;
+    // fprintf (stderr, "%d %d %d\n", (int) i, (int) im, image[im].ubercalDist);
+  }
+  return TRUE;
+}
+  
+int setMmos_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
+
+  int i;
+  off_t N;
+
+  Image *image = getimages (&N, NULL);
+
+  fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE);
+
+  if (PoorImages) {
+    // XXX use bad stars and measurements for PoorImages? or not?
+    // IMAGE_BAD = STAR_BAD = MEAS_BAD = 0;
+    IMAGE_BAD = 0;
+  }
+
+  off_t Nmax = 0;
+  for (i = 0; i < Nmosaic; i++) {
+    Nmax = MAX (Nmax, N_onMosaic[i]);
+  }
+
+  SetMmosInfo summary;
+  SetMmosInfoInit (&summary, Nmax, FALSE, PoorImages);
+
+  pthread_attr_t attr;
+  pthread_attr_init (&attr);
+  pthread_attr_setdetachstate (&attr, PTHREAD_CREATE_DETACHED);
+  
+  pthread_t *threads;
+  ALLOCATE (threads, pthread_t, NTHREADS);
+
+  ThreadInfo *threadinfo;
+  ALLOCATE (threadinfo, ThreadInfo, NTHREADS);
+
+  // launch N worker threads
+  for (i = 0; i < NTHREADS; i++) {
+    threadinfo[i].entry = i;
+    threadinfo[i].state = THREAD_RUN;
+    threadinfo[i].catalog  =  catalog;
+    threadinfo[i].image    =    image;
+    threadinfo[i].flatcorr = flatcorr;
+    SetMmosInfoInit (&threadinfo[i].info, Nmax, FALSE, PoorImages);
+    pthread_create (&threads[i], NULL, setMmos_worker, &threadinfo[i]);
+  }
+  pthread_attr_destroy (&attr);
+
+  // wait until all threads have finished
+  while (1) {
+    int allDone = TRUE;
+    for (i = 0; i < NTHREADS; i++) {
+      if (threadinfo[i].state == THREAD_RUN) allDone = FALSE;
+    }
+    if (allDone) {
+      break;
+    }
+    usleep (500000);
+  }
+
+  // all threads are done, free the threads array and grab the info
+  free (threads);
+  
+  // report stats & summary from the threads
+  for (i = 0; i < NTHREADS; i++) {
+    fprintf (stderr, "setMmos thread %d : %d mosaics marked having too few measurements (Nbad: %d, Ncal: %d, Ngrid: %d, Nrel: %d, Nsys: %d)\n", 
+             i, 
+             threadinfo[i].info.Nfew, 
+             threadinfo[i].info.Nbad, 
+             threadinfo[i].info.Ncal, 
+             threadinfo[i].info.Ngrid, 
+             threadinfo[i].info.Nrel, 
+             threadinfo[i].info.Nsys);
+    SetMmosInfoAccum (&summary, &threadinfo[i].info);
+  }
+  fprintf (stderr, "total : %d mosaics marked having too few measurements (Nbad: %d, Ncal: %d, Ngrid: %d, Nrel: %d, Nsys: %d)\n", 
+           summary.Nfew, 
+           summary.Nbad, 
+           summary.Ncal, 
+           summary.Ngrid, 
+           summary.Nrel, 
+           summary.Nsys);
+  free (threadinfo);
+
+  return TRUE;
+}
+
+void *setMmos_worker (void *data) {
+
+  ThreadInfo *threadinfo = data;
+
+  SetMmosInfo results;
+  SetMmosInfoInit (&results, threadinfo->info.Nmax, TRUE, threadinfo->info.PoorImages); // allocate list, dlist arrays here
+
+  while (1) {
+
+    off_t i = getNextMosaicForThread();
+    if (i == -1) {
+      threadinfo->state = THREAD_DONE;
+      return NULL;
+    }
+
+    Catalog *catalog = threadinfo->catalog;
+    FlatCorrectionTable *flatcorr = threadinfo->flatcorr;
+    Image *image = threadinfo->image;
+
+    setMmos_mosaic (mosaic, i, image, catalog, &results, flatcorr);
+    SetMmosInfoAccum (&threadinfo->info, &results);
+  }
+
+  SetMmosInfoFree (&results);
+  return NULL;
+}
+
 // When we rationalize the images/mosaics, we are driving the negative cloud images back
 // to 0.0.  At the same time, we make a guess to the effective impact on all other images,

branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/StarOps.c

@@ -2 +2 @@
 
 static int Nmax;
-static double *list;
-static double *dlist;
-
-// When we rationalize the images/mosaics, we are driving the negative cloud images back
-// to 0.0.  At the same time, we make a guess to the effective impact on all other images,
-// driven by the coupling of common stars.  This array carries the impact of those offsets
-// on each star
-static double *Moffset;
-
+
+typedef struct {
+  int Nfew;
+  int Nsys;
+  int Nbad;
+  int Ncal;
+  int Nmos;
+  int Ngrid;
+  double *list;
+  double *dlist;
+} SetMrelInfo;
+
+enum {THREAD_RUN, THREAD_DONE};
+
+typedef struct {
+  int entry;
+  int state;
+  Catalog *catalog;
+  int Ncatalog;
+  FlatCorrectionTable *flatcorr;
+  SetMrelInfo summary;
+} ThreadInfo;
+
+void *setMrel_worker (void *data);
+int setMrel_threaded (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr);
+int setMrel_catalog (Catalog *catalog, int Nc, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt);
+
+// we want to allocate the stats list,dlist arrays only once (or once per thread).
+// this function finds the largest array so we can allocate that max size when needed
 void initMrel (Catalog *catalog, int Ncatalog) {
 
@@ -21 +41 @@
     }
   }
-
-  ALLOCATE (list,    double, MAX (1, Nmax));
-  ALLOCATE (dlist,   double, MAX (1, Nmax));
-  ALLOCATE (Moffset, double, MAX (1, Nmax));
 }  
 
@@ -47 +63 @@
 }
 
+void SetMrelInfoInit (SetMrelInfo *results, int allocLists) {
+  results->Nfew  = 0;
+  results->Nsys  = 0;
+  results->Nbad  = 0;
+  results->Ncal  = 0;
+  results->Nmos  = 0;
+  results->Ngrid = 0;
+  if (allocLists) {
+    ALLOCATE (results->list,  double, Nmax);
+    ALLOCATE (results->dlist, double, Nmax);
+  }
+}
+
+void SetMrelInfoFree (SetMrelInfo *results) {
+  free (results->list);
+  free (results->dlist);
+}
+
+void SetMrelInfoAccum (SetMrelInfo *summary, SetMrelInfo *results) {
+  summary->Nfew  += results->Nfew ;
+  summary->Nsys  += results->Nsys ;
+  summary->Nbad  += results->Nbad ;
+  summary->Ncal  += results->Ncal ;
+  summary->Nmos  += results->Nmos ;
+  summary->Ngrid += results->Ngrid;
+}
+
+// mutex to lock setMrel_worker operations 
+static pthread_mutex_t setMrel_mutex = PTHREAD_MUTEX_INITIALIZER;
+static int nextCatalog = 0;
+
+// we have an array of catalogs (catalog, Ncatalog).  we need to hand out catalogs one at a time to
+// the worker threads as they need
+off_t getNextCatalogForThread (int Ncatalog) {
+
+  pthread_mutex_lock (&setMrel_mutex);
+  if (nextCatalog >= Ncatalog) {
+    pthread_mutex_unlock (&setMrel_mutex);
+    return (-1);
+  }
+  int thisCatalog = nextCatalog;
+  nextCatalog ++;
+
+  pthread_mutex_unlock (&setMrel_mutex);
+  return (thisCatalog);
+}
+
 int setMrel (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
 
+  int i;
+
+  if (NTHREADS) {
+    int status = setMrel_threaded (catalog, Ncatalog, flatcorr);
+    return status;
+  }
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  SetMrelInfo summary, results;
+  SetMrelInfoInit (&summary, FALSE);
+  SetMrelInfoInit (&results, TRUE);
+
+  for (i = 0; i < Ncatalog; i++) {
+    setMrel_catalog (catalog, i, flatcorr, &results, Nsecfilt);
+    SetMrelInfoAccum (&summary, &results);
+  }
+  fprintf (stderr, "%d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys);
+
+  SetMrelInfoFree (&results);
+
+  return (TRUE);
+}
+
+// set the Mrel values for the specified catalog
+// NOTE: here 'catalog' is a pointer to a specific catalog, not the root of the array
+int setMrel_catalog (Catalog *catalog, int Nc, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) {
+
   off_t j, k, m;
-  int i, N, Nfew, Nsys, Nbad, Ncal, Nmos, Ngrid;
+  int N;
   float Msys, Mcal, Mmos, Mgrid;
   StatType stats;
 
+  double *list  = results->list;
+  double *dlist = results->dlist;
+
+  SetMrelInfoInit (results, FALSE); // do not allocate list,dlist arrays
+
+  for (j = 0; j < catalog[Nc].Naverage; j++) {
+    // XXX accumulate all secfilt values in a single pass?
+
+    int minUbercalDist = 1000;
+
+    int Ns;
+    for (Ns = 0; Ns < Nphotcodes; Ns++) {
+
+      int thisCode = photcodes[Ns][0].code;
+      int Nsec = GetPhotcodeNsec(thisCode);
+
+      /* calculate the average mag in this SEC photcode for a single star */
+
+      // skip bad stars
+      if (catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;
+
+      N = 0;
+      m = catalog[Nc].averageT[j].measureOffset;
+      for (k = 0; k < catalog[Nc].averageT[j].Nmeasure; k++, m++) {
+
+        // skip measurements that do not match the current photcode
+        PhotCode *code = GetPhotcodebyCode (catalog[Nc].measureT[m].photcode);
+        if (!code) continue;
+        if (code->equiv != thisCode) { continue; }
+
+        if (catalog[Nc].measureT[m].dbFlags & MEAS_BAD) { results->Nbad ++; continue; }
+
+        if (getImageEntry (m, Nc) < 0) {
+          Mcal = Mmos = Mgrid = 0;
+        } else {
+          Mcal  = getMcal  (m, Nc, flatcorr, catalog);
+          if (isnan(Mcal))  { results->Ncal ++; continue; }
+          Mmos  = getMmos  (m, Nc);
+          if (isnan(Mmos))  { results->Nmos ++; continue; }
+          Mgrid = getMgrid (m, Nc);
+          if (isnan(Mgrid)) { results->Ngrid++; continue; }
+        }
+
+        Msys = PhotSysTiny (&catalog[Nc].measureT[m], &catalog[Nc].averageT[j], &catalog[Nc].secfilt[j*Nsecfilt]);
+        if (isnan(Msys)) { results->Nsys++; continue; }
+        list[N] = Msys - Mcal - Mmos - Mgrid;
+
+        int myUbercalDist = getUbercalDist(m,Nc);
+        minUbercalDist = MIN(minUbercalDist, myUbercalDist);
+
+        // dlist gives the error, which is used as the weight in WT_MEAN.
+        // we can modify the resulting weight in a few ways:
+        // 1) MIN_ERROR guarantees a floor
+        // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode
+        // 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages
+        // 4) some reference photcode of some kind can be specified as fixed and have a high weight
+        dlist[N] = MAX (hypot(catalog[Nc].measureT[m].dM, code->photomErrSys), MIN_ERROR);
+
+        // up-weight the ubercal values (or convergence can take a long time...)
+        if (catalog[Nc].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
+          dlist[N] = MAX (0.1*catalog[Nc].measureT[m].dM, MIN_ERROR);
+        }
+
+        // tie down reference photometry if the -refcode (code) option is selected
+        // eg, -refcode g_SDSS
+        // this probably makes no sense in the context of multifilter analysis
+        if (refPhotcode) {
+          if (code->code == refPhotcode->code) {
+            // tiny error -> large weight
+            // dlist[N] = MAX (0.01*catalog[Nc].measureT[m].dM, MIN_ERROR);
+            dlist[N] = 0.0001;
+          }
+        }
+        N++;
+      }
+
+      // when performing the grid analysis, STAR_TOOFEW will be set to 1;
+      if (N <= STAR_TOOFEW) { /* too few measurements */
+        // fprintf (f, "%10.6f %10.6f %d %d %d\n", catalog[Nc].averageT[j].R, catalog[Nc].averageT[j].D, catalog[Nc].measureT[catalog[Nc].averageT[j].measureOffset].imageID, N, STAR_TOOFEW); 
+        catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_STAR_FEW;
+        results->Nfew ++;
+      } else {
+        catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_STAR_FEW;
+      } 
+
+      liststats (list, dlist, N, &stats);
+        
+      catalog[Nc].secfilt[Nsecfilt*j+Nsec].M  = stats.mean;
+      catalog[Nc].secfilt[Nsecfilt*j+Nsec].dM = stats.sigma;
+      catalog[Nc].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq + 1e-4) : NAN_S_SHORT;
+        
+      catalog[Nc].secfilt[Nsecfilt*j+Nsec].ubercalDist = minUbercalDist;
+    }
+  }
+
+  // values which need to be passed out: Nfew, Nbad, Ncal, Nmos, Ngrid, Nsys
+  return (TRUE);
+}
+
+int setMrel_threaded (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
+
+  int i;
+
+  SetMrelInfo summary;
+  SetMrelInfoInit (&summary, FALSE);
+
+  pthread_attr_t attr;
+  pthread_attr_init (&attr);
+  pthread_attr_setdetachstate (&attr, PTHREAD_CREATE_DETACHED);
+  
+  pthread_t *threads;
+  ALLOCATE (threads, pthread_t, NTHREADS);
+
+  ThreadInfo *threadinfo;
+  ALLOCATE (threadinfo, ThreadInfo, NTHREADS);
+
+  // launch N worker threads
+  for (i = 0; i < NTHREADS; i++) {
+    threadinfo[i].entry = i;
+    threadinfo[i].state = THREAD_RUN;
+    threadinfo[i].catalog  =  catalog;
+    threadinfo[i].Ncatalog = Ncatalog;
+    threadinfo[i].flatcorr = flatcorr;
+    SetMrelInfoInit (&threadinfo[i].summary, FALSE);
+    pthread_create (&threads[i], NULL, setMrel_worker, &threadinfo[i]);
+  }
+  pthread_attr_destroy (&attr);
+
+  // wait until all threads have finished
+  while (1) {
+    int allDone = TRUE;
+    for (i = 0; i < NTHREADS; i++) {
+      if (threadinfo[i].state == THREAD_RUN) allDone = FALSE;
+    }
+    if (allDone) {
+      break;
+    }
+    usleep (500000);
+  }
+
+  // all threads are done, free the threads array and grab the info
+  free (threads);
+  
+  // report stats & summary from the threads
+  for (i = 0; i < NTHREADS; i++) {
+    fprintf (stderr, "setMrel thread %d : %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", 
+             i, 
+             threadinfo[i].summary.Nfew, 
+             threadinfo[i].summary.Nbad, 
+             threadinfo[i].summary.Ncal, 
+             threadinfo[i].summary.Nmos, 
+             threadinfo[i].summary.Ngrid, 
+             threadinfo[i].summary.Nsys);
+    SetMrelInfoAccum (&summary, &threadinfo[i].summary);
+  }
+  fprintf (stderr, "total : %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys);
+  free (threadinfo);
+
+  return TRUE;
+}
+
+void *setMrel_worker (void *data) {
+
+  ThreadInfo *threadinfo = data;
+
   int Nsecfilt = GetPhotcodeNsecfilt ();
-  Nfew = Nsys = Nbad = Ncal = Nmos = Ngrid = 0;
-
-  for (i = 0; i < Ncatalog; i++) {
-    for (j = 0; j < catalog[i].Naverage; j++) {
-      // XXX accumulate all secfilt values in a single pass?
-
-      int minUbercalDist = 1000;
-
-      int Ns;
-      for (Ns = 0; Ns < Nphotcodes; Ns++) {
-
-        int thisCode = photcodes[Ns][0].code;
-        int Nsec = GetPhotcodeNsec(thisCode);
-
-        /* calculate the average mag in this SEC photcode for a single star */
-
-        // skip bad stars
-        if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;
-
-        N = 0;
-        m = catalog[i].averageT[j].measureOffset;
-        for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
-
-          // skip measurements that do not match the current photcode
-          PhotCode *code = GetPhotcodebyCode (catalog[i].measureT[m].photcode);
-          if (!code) continue;
-          if (code->equiv != thisCode) { continue; }
-
-          if (catalog[i].measureT[m].dbFlags & MEAS_BAD) { Nbad ++; continue; }
-
-          if (getImageEntry (m, i) < 0) {
-            Mcal = Mmos = Mgrid = 0;
-          } else {
-            Mcal  = getMcal  (m, i, flatcorr, catalog);
-            if (isnan(Mcal)) { Ncal ++; continue; }
-            Mmos  = getMmos  (m, i);
-            if (isnan(Mmos)) { Nmos ++; continue; }
-            Mgrid = getMgrid (m, i);
-            if (isnan(Mgrid)) { Ngrid++; continue; }
-          }
-
-          Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt]);
-          if (isnan(Msys)) { Nsys++; continue; }
-          list[N] = Msys - Mcal - Mmos - Mgrid;
-
-          int myUbercalDist = getUbercalDist(m,i);
-          minUbercalDist = MIN(minUbercalDist, myUbercalDist);
-
-          // dlist gives the error, which is used as the weight in WT_MEAN.
-          // we can modify the resulting weight in a few ways:
-          // 1) MIN_ERROR guarantees a floor
-          // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode
-          // 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages
-          // 4) some reference photcode of some kind can be specified as fixed and have a high weight
-          dlist[N] = MAX (hypot(catalog[i].measureT[m].dM, code->photomErrSys), MIN_ERROR);
-
-          // up-weight the ubercal values (or convergence can take a long time...)
-          if (catalog[i].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
-            dlist[N] = MAX (0.1*catalog[i].measureT[m].dM, MIN_ERROR);
-          }
-
-          // tie down reference photometry if the -refcode (code) option is selected
-          // eg, -refcode g_SDSS
-          // this probably makes no sense in the context of multifilter analysis
-          if (refPhotcode) {
-            if (code->code == refPhotcode->code) {
-              // tiny error -> large weight
-              // dlist[N] = MAX (0.01*catalog[i].measureT[m].dM, MIN_ERROR);
-              dlist[N] = 0.0001;
-            }
-          }
-          N++;
-        }
-
-        // when performing the grid analysis, STAR_TOOFEW will be set to 1;
-        if (N <= STAR_TOOFEW) { /* too few measurements */
-          // fprintf (f, "%10.6f %10.6f %d %d %d\n", catalog[i].averageT[j].R, catalog[i].averageT[j].D, catalog[i].measureT[catalog[i].averageT[j].measureOffset].imageID, N, STAR_TOOFEW); 
-          catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_STAR_FEW;
-          Nfew ++;
-        } else {
-          catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_STAR_FEW;
-        }       
-
-        liststats (list, dlist, N, &stats);
-        
-        catalog[i].secfilt[Nsecfilt*j+Nsec].M  = stats.mean;
-        catalog[i].secfilt[Nsecfilt*j+Nsec].dM = stats.sigma;
-        catalog[i].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq + 1e-4) : NAN_S_SHORT;
-        
-        catalog[i].secfilt[Nsecfilt*j+Nsec].ubercalDist = minUbercalDist;
-      }
-    }
-  }
-  fprintf (stderr, "%d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", Nfew, Nbad, Ncal, Nmos, Ngrid, Nsys);
-
-  return (TRUE);
+
+  SetMrelInfo results;
+  SetMrelInfoInit (&results, TRUE); // allocate list, dlist arrays here
+
+  while (1) {
+
+    off_t i = getNextCatalogForThread(threadinfo->Ncatalog);
+    if (i == -1) {
+      threadinfo->state = THREAD_DONE;
+      return NULL;
+    }
+
+    Catalog *catalog = threadinfo->catalog;
+    FlatCorrectionTable *flatcorr = threadinfo->flatcorr;
+
+    setMrel_catalog (catalog, i, flatcorr, &results, Nsecfilt);
+    SetMrelInfoAccum (&threadinfo->summary, &results);
+  }
+
+  SetMrelInfoFree (&results);
+  return NULL;
 }
 

branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/args.c

@@ -74 +74 @@
   }
 
+  NTHREADS = 0;
+  if ((N = get_argument (argc, argv, "-threads"))) {
+    remove_argument (N, &argc, argv);
+    NTHREADS = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
   // XXX for the moment, make this selection manual.  it needs to be automatic 
   // based on the state of the SkyTable

branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/load_images.c

@@ -15 +15 @@
   off_t     *LineNumber;
   struct timeval start, stop;
+  char *inSubset;
 
   SkyTable *sky = NULL;
@@ -34 +35 @@
   MARKTIME("read image table: %f sec\n", dtime);
 
+  // allocate and init an array to identify the images included in the subset
+  ALLOCATE (inSubset, char, Nimage);
+  memset (inSubset, 0, Nimage);
+
   // determine the populated SkyRegions overlapping the requested area
 
@@ -45 +50 @@
     // select the images which overlap the selected sky regions
     // 'subset' points to a new copy of the data (different from 'image')
-    subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset);
+    subset = select_images (skylist, image, Nimage, inSubset, &LineNumber, &Nsubset);
     MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime);
   } else {
@@ -59 +64 @@
       // select the images which overlap the selected sky regions
       // 'subset' points to a new copy of the data (different from 'image')
-      subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset);
+      subset = select_images (skylist, image, Nimage, inSubset, &LineNumber, &Nsubset);
       MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime);
 
@@ -66 +71 @@
       off_t *LineNumberExtra, i;
 
-      subsetExtra = select_images (extraList, image, Nimage, &LineNumberExtra, &NsubsetExtra);
+      subsetExtra = select_images (extraList, image, Nimage, inSubset, &LineNumberExtra, &NsubsetExtra);
       MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime);
 
@@ -83 +88 @@
       // select the images which overlap the selected sky regions
       // 'subset' points to a new copy of the data (different from 'image')
-      subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset);
+      subset = select_images (skylist, image, Nimage, inSubset, &LineNumber, &Nsubset);
       MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime);
     }
   } 
 
+  // match chips to mosaics (if applicable)
+  initMosaics (subset, Nsubset, image, inSubset, Nimage);
+  MARKTIME("init mosaics: %f sec\n", dtime);
+  
   // save the subset of images in the static reference in ImageOps, set up indexes
   initImages (subset, LineNumber, Nsubset);
   MARKTIME("init images: %f sec\n", dtime);
 
-  // match chips to mosaics (if applicable)
-  initMosaics (subset, Nsubset);
-  MARKTIME("init mosaics: %f sec\n", dtime);
-  
   return (skylist);
 }

branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/relphot.c

@@ -121 +121 @@
         plot_scatter (catalog, Ncatalog, flatcorr); 
       }
-      setMrel  (catalog, Ncatalog, flatcorr);
+      setMrel  (catalog, Ncatalog, flatcorr); // threaded
       if (PLOTSTUFF) {
         plot_scatter (catalog, Ncatalog, flatcorr); 

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

@@ -22 +22 @@
   fprintf (stderr, MSG, __VA_ARGS__); }
 
-Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage) {
+Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, char *inSubset, off_t **LineNumber, off_t *Nimage) {
   
   Image *image;
@@ -125 +125 @@
   double RmidSkyRegion = 0.5*(RminSkyRegion + RmaxSkyRegion);
   MARKTIME("create sky region coords: %f sec\n", dtime);
+
+  // if we overlap 0,360, RminSkyRegion will be 180 (or so)
 
   dsortindex (RmaxSky, index, skylist[0].Nregions);
@@ -252 +254 @@
   found_it:
     image[nimage] = timage[i]; 
+    inSubset[i] = TRUE;
     /* always allow 'few' images to succeed, if possible (new images / detections may have
      * been added) */