IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 33797


Ignore:
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:
7 edited

Legend:

Unmodified
Added
Removed

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

    r33651 r33797  
    33# include <kapa.h>
    44# include <signal.h>
     5# include <pthread.h>
    56
    67/* # define GRID_V1 */
    … …  
    111112int    PARALLEL_MANUAL;
    112113int    PARALLEL_SERIAL;
     114
     115int    NTHREADS;
    113116
    114117int    VERBOSE;
    … …  
    236239void          initMosaicBins      PROTO((Catalog *catalog, int Ncatalog, int doMosaicList));
    237240void          initMosaicGrid      PROTO((Image *image, off_t Nimage));
    238 void          initMosaics         PROTO((Image *image, off_t Nimage));
     241void          initMosaics         PROTO((Image *subset, off_t Nsubset, Image *image, char *inSubset, off_t Nimage));
    239242void          initMrel            PROTO((Catalog *catalog, int Ncatalog));
    240243void          initialize          PROTO((int argc, char **argv));
    … …  
    246249
    247250SkyList      *load_images         PROTO((FITS_DB *db, char *regionName, SkyRegion *region));
    248 Image        *select_images       PROTO((SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage));
     251Image        *select_images       PROTO((SkyList *skylist, Image *timage, off_t Ntimage, char *inSubset, off_t **LineNumber, off_t *Nimage));
    249252
    250253int           main                PROTO((int argc, char **argv));

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

    r33736 r33797  
    5050  ALLOCATE (mosaic, Mosaic, NMOSAIC);
    5151
    52   ALLOCATE (MosaicToImage,  off_t *, NMOSAIC);
     52  ALLOCATE (MosaicToImage, off_t *, NMOSAIC);
    5353  ALLOCATE (MosaicN_Image, off_t,   NMOSAIC);
    5454  ALLOCATE (MosaicN_IMAGE, off_t,   NMOSAIC);
    … …  
    9595    if (found) continue;
    9696   
    97     /* a new mosaic, define ranges */
     97    /* a new mosaic, define ranges -- preserve the original values incase this image is not used */
    9898    mosaic[Nmosaic].start = start;
    9999    mosaic[Nmosaic].stop  = stop;
    100     mosaic[Nmosaic].Mcal  = 0.0;
    101     mosaic[Nmosaic].dMcal = 0.0;
     100    mosaic[Nmosaic].Mcal  = 0.0; // note : mosaic stores only offsets relative to the original image values
     101    mosaic[Nmosaic].dMcal = 0.0; // note : at the end, mosaic.Mcal is added back to the input images
    102102    mosaic[Nmosaic].dMsys = 0.0;
    103103    mosaic[Nmosaic].Xm    = 0.0;
    … …  
    152152}
    153153
    154 /* find mosaic frames (unique time periods & photcode name matches mosaic) */
    155 void initMosaics (Image *image, off_t Nimage) {
     154/* find mosaic frames (unique time periods) (NOTE : we do NOT require matching photcodes...)
     155   this function will also identify the images NOT in the subset which belong to a selected mosaic
     156 */
     157void initMosaics (Image *subset, off_t Nsubset, Image *image, char *inSubset, off_t Nimage) {
    156158
    157159  off_t i, j, status, found, NMOSAIC, *MosaicN_IMAGE;
    … …  
    166168   */
    167169
    168   // generate a list of unique start times (these define the mosaics)
    169   ALLOCATE (startTimes, unsigned int, Nimage);
    170   for (i = 0; i < Nimage; i++) {
    171     startTimes[i] = image[i].tzero;
    172   }
    173   sort_times (startTimes, Nimage);
     170  // generate a list of all subset image start times
     171  ALLOCATE (startTimes, unsigned int, Nsubset);
     172  for (i = 0; i < Nsubset; i++) {
     173    startTimes[i] = subset[i].tzero;
     174  }
     175  sort_times (startTimes, Nsubset);
    174176 
    175177  Nmosaic = 0;
    … …  
    178180  startTimesMosaic[0] = startTimes[0];
    179181
    180   for (i = 0; i < Nimage; i++) {
     182  // generate a list of the unique start times (these define the mosaics)
     183  for (i = 0; i < Nsubset; i++) {
    181184    if (startTimes[i] < startTimesMosaic[Nmosaic]) {
    182185      fprintf (stderr, "error?\n");
    … …  
    192195  }
    193196  Nmosaic ++;
     197
     198  int Nskip, Nmiss, Nmark;
     199  Nskip = Nmiss = Nmark = 0;
     200  // find any start times which match unselected images
     201  for (i = 0; i < Nimage; i++) {
     202    if (inSubset[i]) {
     203      Nskip ++;
     204      continue;
     205    }
     206   
     207    /* select valid mosaic images by photcode */
     208    pname = GetPhotcodeNamebyCode (image[i].photcode);
     209    if (!pname) continue;
     210    status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME));
     211    if (status) continue;
     212
     213    /* set image time range */
     214    start = subset[i].tzero;
     215    stop  = subset[i].tzero + MAX(1.01*subset[i].trate*subset[i].NY, 1);
     216
     217    /* find a matching mosaic */
     218    j = findMosaic(startTimesMosaic, Nmosaic, start);
     219    if (j == -1) {
     220      // no matching mosaic
     221      Nmiss ++;
     222    } else {
     223      // mark this mosaic as bad
     224      Nmark ++;
     225    }
     226  }
     227  fprintf (stderr, "%d images, %d skip, %d miss, %d mark\n", (int) Nimage, (int) Nskip, (int) Nmiss, (int) Nmark);
    194228
    195229  // now I have a list of uniq start times, and they are in order
    … …  
    219253  }
    220254
    221   ALLOCATE (ImageToMosaic, off_t, Nimage); // mosaic to which image belongs
     255  ALLOCATE (ImageToMosaic, off_t, Nsubset); // mosaic to which image belongs
    222256
    223257  // assign each image to a mosaic
    224   for (i = 0; i < Nimage; i++) {
     258  for (i = 0; i < Nsubset; i++) {
    225259    ImageToMosaic[i] = -1;
    226260
    227261    /* select valid mosaic images by photcode */
    228     pname = GetPhotcodeNamebyCode (image[i].photcode);
     262    pname = GetPhotcodeNamebyCode (subset[i].photcode);
    229263    status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME));
    230264    if (status) continue;
    231265
    232266    /* set image time range */
    233     // start = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1);
    234     // stop  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
    235 
    236     start = image[i].tzero;
    237     stop  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
     267    // start = subset[i].tzero - MAX(0.01*subset[i].trate*subset[i].NY, 1);
     268    // stop  = subset[i].tzero + MAX(1.01*subset[i].trate*subset[i].NY, 1);
     269
     270    start = subset[i].tzero;
     271    stop  = subset[i].tzero + MAX(1.01*subset[i].trate*subset[i].NY, 1);
    238272
    239273    j = findMosaic(startTimesMosaic, Nmosaic, start);
    … …  
    268302    mosaic[j].dMsys = 0.0;
    269303    mosaic[j].Xm    = 0.0;
    270     mosaic[j].flags  = image[i].flags;
    271     mosaic[j].secz  = image[i].secz;
    272     mosaic[j].photcode = GetPhotcodeEquivCodebyCode (image[i].photcode);
     304    mosaic[j].flags  = subset[i].flags;
     305    mosaic[j].secz  = subset[i].secz;
     306    mosaic[j].photcode = GetPhotcodeEquivCodebyCode (subset[i].photcode);
    273307  }
    274308
    … …  
    278312  free (startTimesMosaic);
    279313
    280   initMosaicGrid (image, Nimage);
    281 
    282   fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nimage, (int) Nmosaic);
     314  initMosaicGrid (subset, Nsubset);
     315
     316  fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nsubset, (int) Nmosaic);
    283317  return;
    284318}
    … …  
    569603}
    570604
     605typedef struct {
     606  int Nfew;
     607  int Nbad;
     608  int Ncal;
     609  int Ngrid;
     610  int Nrel;
     611  int Nsys;
     612  off_t Nmax;
     613  int PoorImages;
     614  double *list;
     615  double *dlist;
     616  double *Mlist;
     617  double *dMlist;
     618} SetMmosInfo;
     619
     620enum {THREAD_RUN, THREAD_DONE};
     621
     622typedef struct {
     623  int entry;
     624  int state;
     625  Catalog *catalog;
     626  Image *image;
     627  FlatCorrectionTable *flatcorr;
     628  SetMmosInfo info;
     629} ThreadInfo;
     630
     631int setMmos_mosaic (Mosaic *mosaic, off_t Nmos, Image *image, Catalog *catalog, SetMmosInfo *info, FlatCorrectionTable *flatcorr);
     632void *setMmos_worker (void *data);
     633int setMmos_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr);
     634
     635void SetMmosInfoInit (SetMmosInfo *info, off_t Nmax, int allocLists, int PoorImages) {
     636  info->Nfew = 0;
     637  info->Nbad = 0;
     638  info->Ncal = 0;
     639  info->Nrel = 0;
     640  info->Ngrid = 0;
     641  info->Nsys = 0;
     642
     643  info->Nmax = Nmax;
     644  info->PoorImages = PoorImages;
     645
     646  if (allocLists) {
     647    ALLOCATE (info->list, double, Nmax);
     648    ALLOCATE (info->dlist, double, Nmax);
     649    ALLOCATE (info->Mlist, double, Nmax);
     650    ALLOCATE (info->dMlist, double, Nmax);
     651  }
     652}
     653
     654void SetMmosInfoFree (SetMmosInfo *info) {
     655  free (info->list);
     656  free (info->dlist);
     657  free (info->Mlist);
     658  free (info->dMlist);
     659}
     660
     661void SetMmosInfoAccum (SetMmosInfo *summary, SetMmosInfo *results) {
     662  summary->Nfew  += results->Nfew ;
     663  summary->Nsys  += results->Nsys ;
     664  summary->Nbad  += results->Nbad ;
     665  summary->Ncal  += results->Ncal ;
     666  summary->Nrel  += results->Nrel ;
     667  summary->Ngrid += results->Ngrid;
     668}
     669
    571670static int npass_output = 0;
    572671
    573 int setMmos (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
     672// mutex to lock setMmos_worker operations
     673static pthread_mutex_t setMmos_mutex = PTHREAD_MUTEX_INITIALIZER;
     674static int nextMosaic = 0;
     675
     676// we have an array of mosaics (mosaic, Nmosaic).  we need to hand out mosaics one at a time to
     677// the worker threads as they need
     678off_t getNextMosaicForThread () {
     679
     680  pthread_mutex_lock (&setMmos_mutex);
     681  if (nextMosaic >= Nmosaic) {
     682    pthread_mutex_unlock (&setMmos_mutex);
     683    return (-1);
     684  }
     685  int thisMosaic = nextMosaic;
     686  nextMosaic ++;
     687
     688  pthread_mutex_unlock (&setMmos_mutex);
     689  return (thisMosaic);
     690}
     691
     692int setMmos_old (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
    574693
    575694  off_t i, j, m, c, n, N, Nmax;
    … …  
    798917}
    799918 
     919int setMmos (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
     920
     921  off_t i, N, Nmax;
     922  Image *image;
     923
     924  if (!MOSAIC_ZEROPT) return (FALSE);
     925  if (FREEZE_MOSAICS) return (FALSE);
     926
     927  if (NTHREADS) {
     928    int status = setMmos_threaded (catalog, PoorImages, flatcorr);
     929    return status;
     930  }
     931
     932  image = getimages (&N, NULL);
     933
     934  fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE);
     935
     936  if (PoorImages) {
     937    // XXX use bad stars and measurements for PoorImages? or not?
     938    // IMAGE_BAD = STAR_BAD = MEAS_BAD = 0;
     939    IMAGE_BAD = 0;
     940  }
     941
     942  Nmax = 0;
     943  for (i = 0; i < Nmosaic; i++) {
     944    Nmax = MAX (Nmax, N_onMosaic[i]);
     945  }
     946
     947  SetMmosInfo info;
     948  SetMmosInfoInit (&info, Nmax, TRUE, PoorImages);
     949
     950  for (i = 0; i < Nmosaic; i++) {
     951    setMmos_mosaic (mosaic, i, image, catalog, &info, flatcorr);
     952  }
     953  SetMmosInfoFree (&info);
     954
     955  npass_output ++;
     956
     957  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);
     958
     959  if (PoorImages) {
     960    IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP;
     961    STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
     962    MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA;
     963  }
     964  return (TRUE);
     965}
     966 
     967int setMmos_mosaic (Mosaic *mosaic, off_t Nmos, Image *image, Catalog *catalog, SetMmosInfo *info, FlatCorrectionTable *flatcorr) {
     968
     969  off_t j;
     970  StatType stats;
     971
     972  double *list   = info->list;
     973  double *dlist  = info->dlist;
     974  double *Mlist  = info->Mlist;
     975  double *dMlist = info->dMlist;
     976
     977  /* on PoorImages run, skip good images */
     978  if (info->PoorImages) {
     979    int bad = mosaic[0].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_SKIP);
     980    if (!bad) return TRUE;
     981  }     
     982
     983  // UBERCAL image: if this is an ubercal image, set minUbercalDist to 0:
     984  // we optionally do not recalibrate images with UBERCAL zero points
     985  if (mosaic[0].flags & ID_IMAGE_PHOTOM_UBERCAL) {
     986    mosaic[0].ubercalDist = 0;
     987    // propagate ubercalDist to the images
     988    for (j = 0; j < MosaicN_Image[Nmos]; j++) {
     989      off_t im = MosaicToImage[Nmos][j];
     990      image[im].ubercalDist = mosaic[0].ubercalDist;
     991      // fprintf (stderr, "%d %d %d\n", (int) i, (int) im, image[im].ubercalDist);
     992    }
     993    if (KEEP_UBERCAL) return TRUE;
     994  }
     995
     996  int Nsecfilt = GetPhotcodeNsecfilt ();
     997
     998  int minUbercalDist = 1000;
     999
     1000  int testImage = FALSE;
     1001  // testImage |= (abs(mosaic[0].start - 1323003245) < 10);
     1002  // testImage |= (abs(mosaic[0].start - 1323003069) < 10);
     1003  // testImage |= (abs(mosaic[0].start - 1323003125) < 10);
     1004  // testImage |= (abs(mosaic[0].start - 1323003300) < 10);
     1005  // testImage |= (abs(mosaic[0].start - 1323003365) < 10);
     1006  // testImage |= (abs(mosaic[0].start - 1323003191) < 10);
     1007  // testImage |= (abs(mosaic[0].start - 1323003014) < 10);
     1008  // testImage |= (abs(mosaic[0].start - 1323003484) < 10);
     1009  // testImage |= (abs(mosaic[0].start - 1323003419) < 10);
     1010  // testImage |= (abs(mosaic[0].start - 1323002949) < 10);
     1011
     1012  FILE *fout = NULL;
     1013  if (testImage) {
     1014    char filename[64];
     1015    snprintf (filename, 64, "test.%05d.%02d.dat", (int) Nmos, npass_output);
     1016    fout = fopen (filename, "w");
     1017  }
     1018
     1019  // number of stars to measure the bright-end scatter
     1020  int Nbright = 0;
     1021
     1022  int N = 0;
     1023  for (j = 0; j < N_onMosaic[Nmos]; j++) {
     1024    float Msys, Mrel, Mcal, Mgrid, Mflat;
     1025     
     1026    off_t m = MosaicToMeasure[Nmos][j];
     1027    off_t c = MosaicToCatalog[Nmos][j];
     1028     
     1029    if (fout) {
     1030      Mcal  = getMcal  (m, c, flatcorr, catalog);
     1031      Mgrid = getMgrid (m, c);
     1032      Mrel  = getMrel  (catalog, m, c);
     1033      Mflat = getMflat (m, c, flatcorr, catalog);
     1034
     1035      off_t n = catalog[c].measureT[m].averef;
     1036      Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]);
     1037
     1038      float delta = Msys - Mrel - Mcal - Mgrid + Mflat;
     1039
     1040      int isBad = (catalog[c].measureT[m].dbFlags & MEAS_BAD);
     1041
     1042      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);
     1043    }
     1044
     1045    if (catalog[c].measureT[m].dbFlags & MEAS_BAD) {
     1046      info->Nbad ++;
     1047      return TRUE;
     1048    }
     1049    Mcal  = getMcal  (m, c, flatcorr, catalog);
     1050    if (isnan(Mcal)) {
     1051      info->Ncal++;
     1052      return TRUE;
     1053    }
     1054    Mgrid = getMgrid (m, c);
     1055    if (isnan(Mgrid)) {
     1056      info->Ngrid ++;
     1057      return TRUE;
     1058    }
     1059    Mrel  = getMrel  (catalog, m, c);
     1060    if (isnan(Mrel)) {
     1061      info->Nrel ++;
     1062      return TRUE;
     1063    }
     1064     
     1065    // image.Mcal is not supposed to include the flat-field correction, so we need to
     1066    // apply that offset as well here for this image (in other words, each detection is
     1067    // being compared to the model, excluding the zero point, Mcal.  The model includes
     1068    // the flat-correction.  NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat)
     1069
     1070    Mflat = getMflat (m, c, flatcorr, catalog);
     1071
     1072    off_t n = catalog[c].measureT[m].averef;
     1073    Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]);
     1074    if (isnan(Msys)) {
     1075      info->Nsys++;
     1076      return TRUE;
     1077    }
     1078
     1079    PhotCode *code = GetPhotcodebyCode (catalog[c].measureT[m].photcode);
     1080    if (!code) goto skip;
     1081    if (code->equiv < 1) goto skip;
     1082    int Nsec = GetPhotcodeNsec (code->equiv);
     1083    if (Nsec == -1) goto skip;
     1084    minUbercalDist = MIN (catalog[c].secfilt[n*Nsecfilt + Nsec].ubercalDist, minUbercalDist);
     1085
     1086  skip:
     1087    list[N]  = Msys - Mrel - Mcal - Mgrid + Mflat;
     1088    dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
     1089    if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) {
     1090      Mlist[Nbright] = list[N];
     1091      dMlist[Nbright] = dlist[N];
     1092      Nbright ++;
     1093    }
     1094    N++;
     1095  }
     1096  /* N_onMosaic[Nmos] is all measurements, N is good measurements */
     1097
     1098  if (fout) {
     1099    fclose (fout);
     1100  }
     1101
     1102  /* too few good measurements or too many bad measurements (skip in PoorImages run) */
     1103  if (!info->PoorImages) {
     1104    int mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*N_onMosaic[Nmos]);
     1105    if (mark) {
     1106      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]); }
     1107      mosaic[0].flags |= ID_IMAGE_PHOTOM_FEW;
     1108      info->Nfew ++;
     1109      if (testImage) {
     1110        fprintf (stderr, "NOTE: *** marked test image poor : %d %d %d***\n", (int) N, (int) IMAGE_TOOFEW, (int) (IMAGE_GOOD_FRACTION*N_onMosaic[Nmos]));
     1111      }
     1112    } else {
     1113      mosaic[0].flags &= ~ID_IMAGE_PHOTOM_FEW;
     1114    }
     1115  }
     1116  liststats (list, dlist, N, &stats);
     1117  if (VERBOSE2 && info->PoorImages) fprintf (stderr, "Mmos: %f %f %d %d\n", stats.mean, stats.sigma, stats.Nmeas, N);
     1118
     1119  mosaic[0].Mcal  = stats.mean;
     1120  mosaic[0].dMcal = stats.error;
     1121  mosaic[0].nFitPhotom = N;
     1122  mosaic[0].Xm    = 100.0*log10(stats.chisq);
     1123
     1124  if (testImage) {
     1125    fprintf (stderr, "test image %d (%d) %f %f %d ... ", (int) Nmos, mosaic[0].start, stats.mean, stats.error, mosaic[0].nFitPhotom);
     1126  }
     1127
     1128  plot_setMcal (list, N, &stats, CLOUD_TOLERANCE);
     1129
     1130  // bright end scatter
     1131  liststats (Mlist, dMlist, Nbright, &stats);
     1132  mosaic[0].dMsys = stats.sigma;
     1133
     1134  if (mosaic[0].Mcal < -CLOUD_TOLERANCE) {
     1135    mosaic[0].Mcal = 0.0;
     1136  }
     1137
     1138  if (testImage) {
     1139    fprintf (stderr, "%f %f\n", mosaic[0].Mcal, mosaic[0].dMsys);
     1140  }
     1141
     1142  // minUbercalDist calculated here is the min value for any star owned by this image
     1143  // since this particular image is tied to that star, bump its distance by 1
     1144  mosaic[0].ubercalDist = minUbercalDist + 1;
     1145
     1146  // propagate ubercalDist to the images
     1147  for (j = 0; j < MosaicN_Image[Nmos]; j++) {
     1148    off_t im = MosaicToImage[Nmos][j];
     1149    image[im].ubercalDist = mosaic[0].ubercalDist;
     1150    // fprintf (stderr, "%d %d %d\n", (int) i, (int) im, image[im].ubercalDist);
     1151  }
     1152  return TRUE;
     1153}
     1154 
     1155int setMmos_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
     1156
     1157  int i;
     1158  off_t N;
     1159
     1160  Image *image = getimages (&N, NULL);
     1161
     1162  fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE);
     1163
     1164  if (PoorImages) {
     1165    // XXX use bad stars and measurements for PoorImages? or not?
     1166    // IMAGE_BAD = STAR_BAD = MEAS_BAD = 0;
     1167    IMAGE_BAD = 0;
     1168  }
     1169
     1170  off_t Nmax = 0;
     1171  for (i = 0; i < Nmosaic; i++) {
     1172    Nmax = MAX (Nmax, N_onMosaic[i]);
     1173  }
     1174
     1175  SetMmosInfo summary;
     1176  SetMmosInfoInit (&summary, Nmax, FALSE, PoorImages);
     1177
     1178  pthread_attr_t attr;
     1179  pthread_attr_init (&attr);
     1180  pthread_attr_setdetachstate (&attr, PTHREAD_CREATE_DETACHED);
     1181 
     1182  pthread_t *threads;
     1183  ALLOCATE (threads, pthread_t, NTHREADS);
     1184
     1185  ThreadInfo *threadinfo;
     1186  ALLOCATE (threadinfo, ThreadInfo, NTHREADS);
     1187
     1188  // launch N worker threads
     1189  for (i = 0; i < NTHREADS; i++) {
     1190    threadinfo[i].entry = i;
     1191    threadinfo[i].state = THREAD_RUN;
     1192    threadinfo[i].catalog  =  catalog;
     1193    threadinfo[i].image    =    image;
     1194    threadinfo[i].flatcorr = flatcorr;
     1195    SetMmosInfoInit (&threadinfo[i].info, Nmax, FALSE, PoorImages);
     1196    pthread_create (&threads[i], NULL, setMmos_worker, &threadinfo[i]);
     1197  }
     1198  pthread_attr_destroy (&attr);
     1199
     1200  // wait until all threads have finished
     1201  while (1) {
     1202    int allDone = TRUE;
     1203    for (i = 0; i < NTHREADS; i++) {
     1204      if (threadinfo[i].state == THREAD_RUN) allDone = FALSE;
     1205    }
     1206    if (allDone) {
     1207      break;
     1208    }
     1209    usleep (500000);
     1210  }
     1211
     1212  // all threads are done, free the threads array and grab the info
     1213  free (threads);
     1214 
     1215  // report stats & summary from the threads
     1216  for (i = 0; i < NTHREADS; i++) {
     1217    fprintf (stderr, "setMmos thread %d : %d mosaics marked having too few measurements (Nbad: %d, Ncal: %d, Ngrid: %d, Nrel: %d, Nsys: %d)\n",
     1218             i,
     1219             threadinfo[i].info.Nfew,
     1220             threadinfo[i].info.Nbad,
     1221             threadinfo[i].info.Ncal,
     1222             threadinfo[i].info.Ngrid,
     1223             threadinfo[i].info.Nrel,
     1224             threadinfo[i].info.Nsys);
     1225    SetMmosInfoAccum (&summary, &threadinfo[i].info);
     1226  }
     1227  fprintf (stderr, "total : %d mosaics marked having too few measurements (Nbad: %d, Ncal: %d, Ngrid: %d, Nrel: %d, Nsys: %d)\n",
     1228           summary.Nfew,
     1229           summary.Nbad,
     1230           summary.Ncal,
     1231           summary.Ngrid,
     1232           summary.Nrel,
     1233           summary.Nsys);
     1234  free (threadinfo);
     1235
     1236  return TRUE;
     1237}
     1238
     1239void *setMmos_worker (void *data) {
     1240
     1241  ThreadInfo *threadinfo = data;
     1242
     1243  SetMmosInfo results;
     1244  SetMmosInfoInit (&results, threadinfo->info.Nmax, TRUE, threadinfo->info.PoorImages); // allocate list, dlist arrays here
     1245
     1246  while (1) {
     1247
     1248    off_t i = getNextMosaicForThread();
     1249    if (i == -1) {
     1250      threadinfo->state = THREAD_DONE;
     1251      return NULL;
     1252    }
     1253
     1254    Catalog *catalog = threadinfo->catalog;
     1255    FlatCorrectionTable *flatcorr = threadinfo->flatcorr;
     1256    Image *image = threadinfo->image;
     1257
     1258    setMmos_mosaic (mosaic, i, image, catalog, &results, flatcorr);
     1259    SetMmosInfoAccum (&threadinfo->info, &results);
     1260  }
     1261
     1262  SetMmosInfoFree (&results);
     1263  return NULL;
     1264}
     1265
    8001266// When we rationalize the images/mosaics, we are driving the negative cloud images back
    8011267// 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

    r33736 r33797  
    22
    33static int Nmax;
    4 static double *list;
    5 static double *dlist;
    6 
    7 // When we rationalize the images/mosaics, we are driving the negative cloud images back
    8 // to 0.0.  At the same time, we make a guess to the effective impact on all other images,
    9 // driven by the coupling of common stars.  This array carries the impact of those offsets
    10 // on each star
    11 static double *Moffset;
    12 
     4
     5typedef struct {
     6  int Nfew;
     7  int Nsys;
     8  int Nbad;
     9  int Ncal;
     10  int Nmos;
     11  int Ngrid;
     12  double *list;
     13  double *dlist;
     14} SetMrelInfo;
     15
     16enum {THREAD_RUN, THREAD_DONE};
     17
     18typedef struct {
     19  int entry;
     20  int state;
     21  Catalog *catalog;
     22  int Ncatalog;
     23  FlatCorrectionTable *flatcorr;
     24  SetMrelInfo summary;
     25} ThreadInfo;
     26
     27void *setMrel_worker (void *data);
     28int setMrel_threaded (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr);
     29int setMrel_catalog (Catalog *catalog, int Nc, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt);
     30
     31// we want to allocate the stats list,dlist arrays only once (or once per thread).
     32// this function finds the largest array so we can allocate that max size when needed
    1333void initMrel (Catalog *catalog, int Ncatalog) {
    1434
    … …  
    2141    }
    2242  }
    23 
    24   ALLOCATE (list,    double, MAX (1, Nmax));
    25   ALLOCATE (dlist,   double, MAX (1, Nmax));
    26   ALLOCATE (Moffset, double, MAX (1, Nmax));
    2743} 
    2844
    … …  
    4763}
    4864
     65void SetMrelInfoInit (SetMrelInfo *results, int allocLists) {
     66  results->Nfew  = 0;
     67  results->Nsys  = 0;
     68  results->Nbad  = 0;
     69  results->Ncal  = 0;
     70  results->Nmos  = 0;
     71  results->Ngrid = 0;
     72  if (allocLists) {
     73    ALLOCATE (results->list,  double, Nmax);
     74    ALLOCATE (results->dlist, double, Nmax);
     75  }
     76}
     77
     78void SetMrelInfoFree (SetMrelInfo *results) {
     79  free (results->list);
     80  free (results->dlist);
     81}
     82
     83void SetMrelInfoAccum (SetMrelInfo *summary, SetMrelInfo *results) {
     84  summary->Nfew  += results->Nfew ;
     85  summary->Nsys  += results->Nsys ;
     86  summary->Nbad  += results->Nbad ;
     87  summary->Ncal  += results->Ncal ;
     88  summary->Nmos  += results->Nmos ;
     89  summary->Ngrid += results->Ngrid;
     90}
     91
     92// mutex to lock setMrel_worker operations
     93static pthread_mutex_t setMrel_mutex = PTHREAD_MUTEX_INITIALIZER;
     94static int nextCatalog = 0;
     95
     96// we have an array of catalogs (catalog, Ncatalog).  we need to hand out catalogs one at a time to
     97// the worker threads as they need
     98off_t getNextCatalogForThread (int Ncatalog) {
     99
     100  pthread_mutex_lock (&setMrel_mutex);
     101  if (nextCatalog >= Ncatalog) {
     102    pthread_mutex_unlock (&setMrel_mutex);
     103    return (-1);
     104  }
     105  int thisCatalog = nextCatalog;
     106  nextCatalog ++;
     107
     108  pthread_mutex_unlock (&setMrel_mutex);
     109  return (thisCatalog);
     110}
     111
    49112int setMrel (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
    50113
     114  int i;
     115
     116  if (NTHREADS) {
     117    int status = setMrel_threaded (catalog, Ncatalog, flatcorr);
     118    return status;
     119  }
     120
     121  int Nsecfilt = GetPhotcodeNsecfilt ();
     122
     123  SetMrelInfo summary, results;
     124  SetMrelInfoInit (&summary, FALSE);
     125  SetMrelInfoInit (&results, TRUE);
     126
     127  for (i = 0; i < Ncatalog; i++) {
     128    setMrel_catalog (catalog, i, flatcorr, &results, Nsecfilt);
     129    SetMrelInfoAccum (&summary, &results);
     130  }
     131  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);
     132
     133  SetMrelInfoFree (&results);
     134
     135  return (TRUE);
     136}
     137
     138// set the Mrel values for the specified catalog
     139// NOTE: here 'catalog' is a pointer to a specific catalog, not the root of the array
     140int setMrel_catalog (Catalog *catalog, int Nc, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) {
     141
    51142  off_t j, k, m;
    52   int i, N, Nfew, Nsys, Nbad, Ncal, Nmos, Ngrid;
     143  int N;
    53144  float Msys, Mcal, Mmos, Mgrid;
    54145  StatType stats;
    55146
     147  double *list  = results->list;
     148  double *dlist = results->dlist;
     149
     150  SetMrelInfoInit (results, FALSE); // do not allocate list,dlist arrays
     151
     152  for (j = 0; j < catalog[Nc].Naverage; j++) {
     153    // XXX accumulate all secfilt values in a single pass?
     154
     155    int minUbercalDist = 1000;
     156
     157    int Ns;
     158    for (Ns = 0; Ns < Nphotcodes; Ns++) {
     159
     160      int thisCode = photcodes[Ns][0].code;
     161      int Nsec = GetPhotcodeNsec(thisCode);
     162
     163      /* calculate the average mag in this SEC photcode for a single star */
     164
     165      // skip bad stars
     166      if (catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;
     167
     168      N = 0;
     169      m = catalog[Nc].averageT[j].measureOffset;
     170      for (k = 0; k < catalog[Nc].averageT[j].Nmeasure; k++, m++) {
     171
     172        // skip measurements that do not match the current photcode
     173        PhotCode *code = GetPhotcodebyCode (catalog[Nc].measureT[m].photcode);
     174        if (!code) continue;
     175        if (code->equiv != thisCode) { continue; }
     176
     177        if (catalog[Nc].measureT[m].dbFlags & MEAS_BAD) { results->Nbad ++; continue; }
     178
     179        if (getImageEntry (m, Nc) < 0) {
     180          Mcal = Mmos = Mgrid = 0;
     181        } else {
     182          Mcal  = getMcal  (m, Nc, flatcorr, catalog);
     183          if (isnan(Mcal))  { results->Ncal ++; continue; }
     184          Mmos  = getMmos  (m, Nc);
     185          if (isnan(Mmos))  { results->Nmos ++; continue; }
     186          Mgrid = getMgrid (m, Nc);
     187          if (isnan(Mgrid)) { results->Ngrid++; continue; }
     188        }
     189
     190        Msys = PhotSysTiny (&catalog[Nc].measureT[m], &catalog[Nc].averageT[j], &catalog[Nc].secfilt[j*Nsecfilt]);
     191        if (isnan(Msys)) { results->Nsys++; continue; }
     192        list[N] = Msys - Mcal - Mmos - Mgrid;
     193
     194        int myUbercalDist = getUbercalDist(m,Nc);
     195        minUbercalDist = MIN(minUbercalDist, myUbercalDist);
     196
     197        // dlist gives the error, which is used as the weight in WT_MEAN.
     198        // we can modify the resulting weight in a few ways:
     199        // 1) MIN_ERROR guarantees a floor
     200        // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode
     201        // 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages
     202        // 4) some reference photcode of some kind can be specified as fixed and have a high weight
     203        dlist[N] = MAX (hypot(catalog[Nc].measureT[m].dM, code->photomErrSys), MIN_ERROR);
     204
     205        // up-weight the ubercal values (or convergence can take a long time...)
     206        if (catalog[Nc].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
     207          dlist[N] = MAX (0.1*catalog[Nc].measureT[m].dM, MIN_ERROR);
     208        }
     209
     210        // tie down reference photometry if the -refcode (code) option is selected
     211        // eg, -refcode g_SDSS
     212        // this probably makes no sense in the context of multifilter analysis
     213        if (refPhotcode) {
     214          if (code->code == refPhotcode->code) {
     215            // tiny error -> large weight
     216            // dlist[N] = MAX (0.01*catalog[Nc].measureT[m].dM, MIN_ERROR);
     217            dlist[N] = 0.0001;
     218          }
     219        }
     220        N++;
     221      }
     222
     223      // when performing the grid analysis, STAR_TOOFEW will be set to 1;
     224      if (N <= STAR_TOOFEW) { /* too few measurements */
     225        // 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);
     226        catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_STAR_FEW;
     227        results->Nfew ++;
     228      } else {
     229        catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_STAR_FEW;
     230      }
     231
     232      liststats (list, dlist, N, &stats);
     233       
     234      catalog[Nc].secfilt[Nsecfilt*j+Nsec].M  = stats.mean;
     235      catalog[Nc].secfilt[Nsecfilt*j+Nsec].dM = stats.sigma;
     236      catalog[Nc].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq + 1e-4) : NAN_S_SHORT;
     237       
     238      catalog[Nc].secfilt[Nsecfilt*j+Nsec].ubercalDist = minUbercalDist;
     239    }
     240  }
     241
     242  // values which need to be passed out: Nfew, Nbad, Ncal, Nmos, Ngrid, Nsys
     243  return (TRUE);
     244}
     245
     246int setMrel_threaded (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
     247
     248  int i;
     249
     250  SetMrelInfo summary;
     251  SetMrelInfoInit (&summary, FALSE);
     252
     253  pthread_attr_t attr;
     254  pthread_attr_init (&attr);
     255  pthread_attr_setdetachstate (&attr, PTHREAD_CREATE_DETACHED);
     256 
     257  pthread_t *threads;
     258  ALLOCATE (threads, pthread_t, NTHREADS);
     259
     260  ThreadInfo *threadinfo;
     261  ALLOCATE (threadinfo, ThreadInfo, NTHREADS);
     262
     263  // launch N worker threads
     264  for (i = 0; i < NTHREADS; i++) {
     265    threadinfo[i].entry = i;
     266    threadinfo[i].state = THREAD_RUN;
     267    threadinfo[i].catalog  =  catalog;
     268    threadinfo[i].Ncatalog = Ncatalog;
     269    threadinfo[i].flatcorr = flatcorr;
     270    SetMrelInfoInit (&threadinfo[i].summary, FALSE);
     271    pthread_create (&threads[i], NULL, setMrel_worker, &threadinfo[i]);
     272  }
     273  pthread_attr_destroy (&attr);
     274
     275  // wait until all threads have finished
     276  while (1) {
     277    int allDone = TRUE;
     278    for (i = 0; i < NTHREADS; i++) {
     279      if (threadinfo[i].state == THREAD_RUN) allDone = FALSE;
     280    }
     281    if (allDone) {
     282      break;
     283    }
     284    usleep (500000);
     285  }
     286
     287  // all threads are done, free the threads array and grab the info
     288  free (threads);
     289 
     290  // report stats & summary from the threads
     291  for (i = 0; i < NTHREADS; i++) {
     292    fprintf (stderr, "setMrel thread %d : %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n",
     293             i,
     294             threadinfo[i].summary.Nfew,
     295             threadinfo[i].summary.Nbad,
     296             threadinfo[i].summary.Ncal,
     297             threadinfo[i].summary.Nmos,
     298             threadinfo[i].summary.Ngrid,
     299             threadinfo[i].summary.Nsys);
     300    SetMrelInfoAccum (&summary, &threadinfo[i].summary);
     301  }
     302  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);
     303  free (threadinfo);
     304
     305  return TRUE;
     306}
     307
     308void *setMrel_worker (void *data) {
     309
     310  ThreadInfo *threadinfo = data;
     311
    56312  int Nsecfilt = GetPhotcodeNsecfilt ();
    57   Nfew = Nsys = Nbad = Ncal = Nmos = Ngrid = 0;
    58 
    59   for (i = 0; i < Ncatalog; i++) {
    60     for (j = 0; j < catalog[i].Naverage; j++) {
    61       // XXX accumulate all secfilt values in a single pass?
    62 
    63       int minUbercalDist = 1000;
    64 
    65       int Ns;
    66       for (Ns = 0; Ns < Nphotcodes; Ns++) {
    67 
    68         int thisCode = photcodes[Ns][0].code;
    69         int Nsec = GetPhotcodeNsec(thisCode);
    70 
    71         /* calculate the average mag in this SEC photcode for a single star */
    72 
    73         // skip bad stars
    74         if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;
    75 
    76         N = 0;
    77         m = catalog[i].averageT[j].measureOffset;
    78         for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
    79 
    80           // skip measurements that do not match the current photcode
    81           PhotCode *code = GetPhotcodebyCode (catalog[i].measureT[m].photcode);
    82           if (!code) continue;
    83           if (code->equiv != thisCode) { continue; }
    84 
    85           if (catalog[i].measureT[m].dbFlags & MEAS_BAD) { Nbad ++; continue; }
    86 
    87           if (getImageEntry (m, i) < 0) {
    88             Mcal = Mmos = Mgrid = 0;
    89           } else {
    90             Mcal  = getMcal  (m, i, flatcorr, catalog);
    91             if (isnan(Mcal)) { Ncal ++; continue; }
    92             Mmos  = getMmos  (m, i);
    93             if (isnan(Mmos)) { Nmos ++; continue; }
    94             Mgrid = getMgrid (m, i);
    95             if (isnan(Mgrid)) { Ngrid++; continue; }
    96           }
    97 
    98           Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt]);
    99           if (isnan(Msys)) { Nsys++; continue; }
    100           list[N] = Msys - Mcal - Mmos - Mgrid;
    101 
    102           int myUbercalDist = getUbercalDist(m,i);
    103           minUbercalDist = MIN(minUbercalDist, myUbercalDist);
    104 
    105           // dlist gives the error, which is used as the weight in WT_MEAN.
    106           // we can modify the resulting weight in a few ways:
    107           // 1) MIN_ERROR guarantees a floor
    108           // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode
    109           // 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages
    110           // 4) some reference photcode of some kind can be specified as fixed and have a high weight
    111           dlist[N] = MAX (hypot(catalog[i].measureT[m].dM, code->photomErrSys), MIN_ERROR);
    112 
    113           // up-weight the ubercal values (or convergence can take a long time...)
    114           if (catalog[i].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
    115             dlist[N] = MAX (0.1*catalog[i].measureT[m].dM, MIN_ERROR);
    116           }
    117 
    118           // tie down reference photometry if the -refcode (code) option is selected
    119           // eg, -refcode g_SDSS
    120           // this probably makes no sense in the context of multifilter analysis
    121           if (refPhotcode) {
    122             if (code->code == refPhotcode->code) {
    123               // tiny error -> large weight
    124               // dlist[N] = MAX (0.01*catalog[i].measureT[m].dM, MIN_ERROR);
    125               dlist[N] = 0.0001;
    126             }
    127           }
    128           N++;
    129         }
    130 
    131         // when performing the grid analysis, STAR_TOOFEW will be set to 1;
    132         if (N <= STAR_TOOFEW) { /* too few measurements */
    133           // 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);
    134           catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_STAR_FEW;
    135           Nfew ++;
    136         } else {
    137           catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_STAR_FEW;
    138         }       
    139 
    140         liststats (list, dlist, N, &stats);
    141        
    142         catalog[i].secfilt[Nsecfilt*j+Nsec].M  = stats.mean;
    143         catalog[i].secfilt[Nsecfilt*j+Nsec].dM = stats.sigma;
    144         catalog[i].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq + 1e-4) : NAN_S_SHORT;
    145        
    146         catalog[i].secfilt[Nsecfilt*j+Nsec].ubercalDist = minUbercalDist;
    147       }
    148     }
    149   }
    150   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);
    151 
    152   return (TRUE);
     313
     314  SetMrelInfo results;
     315  SetMrelInfoInit (&results, TRUE); // allocate list, dlist arrays here
     316
     317  while (1) {
     318
     319    off_t i = getNextCatalogForThread(threadinfo->Ncatalog);
     320    if (i == -1) {
     321      threadinfo->state = THREAD_DONE;
     322      return NULL;
     323    }
     324
     325    Catalog *catalog = threadinfo->catalog;
     326    FlatCorrectionTable *flatcorr = threadinfo->flatcorr;
     327
     328    setMrel_catalog (catalog, i, flatcorr, &results, Nsecfilt);
     329    SetMrelInfoAccum (&threadinfo->summary, &results);
     330  }
     331
     332  SetMrelInfoFree (&results);
     333  return NULL;
    153334}
    154335

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

    r33651 r33797  
    7474  }
    7575
     76  NTHREADS = 0;
     77  if ((N = get_argument (argc, argv, "-threads"))) {
     78    remove_argument (N, &argc, argv);
     79    NTHREADS = atof (argv[N]);
     80    remove_argument (N, &argc, argv);
     81  }
     82
    7683  // XXX for the moment, make this selection manual.  it needs to be automatic
    7784  // based on the state of the SkyTable

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

    r33651 r33797  
    1515  off_t     *LineNumber;
    1616  struct timeval start, stop;
     17  char *inSubset;
    1718
    1819  SkyTable *sky = NULL;
    … …  
    3435  MARKTIME("read image table: %f sec\n", dtime);
    3536
     37  // allocate and init an array to identify the images included in the subset
     38  ALLOCATE (inSubset, char, Nimage);
     39  memset (inSubset, 0, Nimage);
     40
    3641  // determine the populated SkyRegions overlapping the requested area
    3742
    … …  
    4550    // select the images which overlap the selected sky regions
    4651    // 'subset' points to a new copy of the data (different from 'image')
    47     subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset);
     52    subset = select_images (skylist, image, Nimage, inSubset, &LineNumber, &Nsubset);
    4853    MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime);
    4954  } else {
    … …  
    5964      // select the images which overlap the selected sky regions
    6065      // 'subset' points to a new copy of the data (different from 'image')
    61       subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset);
     66      subset = select_images (skylist, image, Nimage, inSubset, &LineNumber, &Nsubset);
    6267      MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime);
    6368
    … …  
    6671      off_t *LineNumberExtra, i;
    6772
    68       subsetExtra = select_images (extraList, image, Nimage, &LineNumberExtra, &NsubsetExtra);
     73      subsetExtra = select_images (extraList, image, Nimage, inSubset, &LineNumberExtra, &NsubsetExtra);
    6974      MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime);
    7075
    … …  
    8388      // select the images which overlap the selected sky regions
    8489      // 'subset' points to a new copy of the data (different from 'image')
    85       subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset);
     90      subset = select_images (skylist, image, Nimage, inSubset, &LineNumber, &Nsubset);
    8691      MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime);
    8792    }
    8893  }
    8994
     95  // match chips to mosaics (if applicable)
     96  initMosaics (subset, Nsubset, image, inSubset, Nimage);
     97  MARKTIME("init mosaics: %f sec\n", dtime);
     98 
    9099  // save the subset of images in the static reference in ImageOps, set up indexes
    91100  initImages (subset, LineNumber, Nsubset);
    92101  MARKTIME("init images: %f sec\n", dtime);
    93102
    94   // match chips to mosaics (if applicable)
    95   initMosaics (subset, Nsubset);
    96   MARKTIME("init mosaics: %f sec\n", dtime);
    97  
    98103  return (skylist);
    99104}

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

    r33651 r33797  
    121121        plot_scatter (catalog, Ncatalog, flatcorr);
    122122      }
    123       setMrel  (catalog, Ncatalog, flatcorr);
     123      setMrel  (catalog, Ncatalog, flatcorr); // threaded
    124124      if (PLOTSTUFF) {
    125125        plot_scatter (catalog, Ncatalog, flatcorr);

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

    r33788 r33797  
    2222  fprintf (stderr, MSG, __VA_ARGS__); }
    2323
    24 Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage) {
     24Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, char *inSubset, off_t **LineNumber, off_t *Nimage) {
    2525 
    2626  Image *image;
    … …  
    125125  double RmidSkyRegion = 0.5*(RminSkyRegion + RmaxSkyRegion);
    126126  MARKTIME("create sky region coords: %f sec\n", dtime);
     127
     128  // if we overlap 0,360, RminSkyRegion will be 180 (or so)
    127129
    128130  dsortindex (RmaxSky, index, skylist[0].Nregions);
    … …  
    252254  found_it:
    253255    image[nimage] = timage[i];
     256    inSubset[i] = TRUE;
    254257    /* always allow 'few' images to succeed, if possible (new images / detections may have
    255258     * been added) */
Note: See TracChangeset for help on using the changeset viewer.