IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jan 11, 2015, 2:14:39 PM (12 years ago)
Author:
eugene
Message:

merge changes from eam branch ipp-20140904

Location:
trunk/Ohana
Files:
14 edited
1 copied

Legend:

Unmodified
Added
Removed

  • trunk/Ohana

  • trunk/Ohana/src/relphot/include/relphot.h

    r37116 r37807  
    452452int setMrel_catalog_alt (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt);
    453453int setMrelAverageExposure (off_t meas, int cat, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT, off_t *found);
    454 int setMrelAverageStack (off_t meas, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure);
     454int setMrelAverageStack (off_t measureOffset, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t Nmeasure);
     455int setMrelAverageForcedWarp (off_t measureOffset, int cat, int pass, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t *found, off_t Nmeasure);
     456
    455457int setGlobalObjStats (Average *average, Measure *measure);
    456458

  • trunk/Ohana/src/relphot/src/BoundaryTreeOps.c

    r36491 r37807  
    5151  if (!tess) return FALSE;
    5252
    53   status = TessellationPrimaryCellIDs(tess, Ntess, tessID, projID, skycellID, ra, dec);
     53  double R = ohana_normalize_angle(ra);
     54  status = TessellationPrimaryCellIDs(tess, Ntess, tessID, projID, skycellID, R, dec);
    5455  return status;
    5556}

  • trunk/Ohana/src/relphot/src/ImageOps.c

    r37116 r37807  
    3535static off_t        *imageIdx; // list of index for image IDs
    3636
     37// QQ // as an alternative, we generate imageSeq, which directly maps imageID -> seq
     38// QQ static off_t        *imageSeq; // list of index for image IDs
     39// QQ
     40// QQ // MAX_ID requires 512M to store the image index
     41// QQ # define MAX_ID 0x8000000
     42// QQ static off_t         minImageID = MAX_ID;
     43// QQ static off_t         maxImageID = 0;
     44
    3745// elsewhere, we have loaded a set of catalogs with measures (catalog[cat].measure[meas])
    3846// each image has N_onImage[ImageIndex] measurements
    … …  
    95103
    96104    TessellationIDsByImageName (&tessID[i], &projectID[i], &skycellID[i], image[i].name);
     105
     106    // QQ minImageID = MIN(minImageID, image[i].imageID);
     107    // QQ maxImageID = MAX(maxImageID, image[i].imageID);
     108    // QQ myAssert (image[i].imageID < MAX_ID, "image IDs too large for index memory");
    97109  }
    98110
    … …  
    100112  // XXX does this break the imageID <-> projectID, etc match?
    101113  llsortpair (imageIDs, imageIdx, Nimage);
     114
     115  // QQ ALLOCATE (imageSeq, off_t, maxImageID + 1);
     116  // QQ for (i = 0; i < maxImageID + 1; i++)  {
     117  // QQ   imageSeq[i] = -1; // not yet assigned
     118  // QQ }
     119  // QQ
     120  // QQ for (i = 0; i < Nimage; i++) {
     121  // QQ   off_t N = image[i].imageID;
     122  // QQ   myAssert (imageSeq[N] == -1, "previously assigned");
     123  // QQ   imageSeq[N] = i;
     124  // QQ }
    102125}
    103126
    … …  
    167190
    168191off_t getImageByID (off_t ID) {
     192
     193  // QQ  if (imageSeq) {
     194  // QQ    myAssert (ID >= minImageID, "oops");
     195  // QQ    myAssert (ID <= maxImageID, "oops");
     196  // QQ    off_t N = imageSeq[ID];
     197  // QQ    return N;
     198  // QQ  }
    169199
    170200  // we have a pair of vectors (imageIDs, imageIdx) sorted by imageIDs

  • trunk/Ohana/src/relphot/src/ImageTable.c

    r36630 r37807  
    3232  }
    3333
     34  // XXX do not make the chip match -- we have not loaded the PHU entries (and do not need them here)
     35  // BuildChipMatch (image, Nimage);
     36
    3437  *nimage = Nimage;
    3538  return image;

  • trunk/Ohana/src/relphot/src/MosaicOps.c

    r37037 r37807  
    520520      m = MosaicToImage[i][j];
    521521
    522       if (!FindMosaicForImage (image, Nimage, m)) {
    523         if (VERBOSE2) fprintf (stderr, "cannot find mosaic for "OFF_T_FMT"\n", i);
    524         continue;
    525       }
    526 
    527522      NX = image[m].NX;
    528523      NY = image[m].NY;
    529       XY_to_RD (&R, &D, 0.5*NX, 0.5*NY, &image[m].coords);
     524      if (!XY_to_RD (&R, &D, 0.5*NX, 0.5*NY, &image[m].coords)) continue;
    530525      R = ohana_normalize_angle_to_midpoint (R, 180.0);
    531526
    … …  
    587582    }
    588583
    589     strcpy (mosaic[i].coords.ctype, "DEC--TAN");
     584    InitCoords (&mosaic[i].coords, "DEC--TAN");
    590585    mosaic[i].coords.crval1 = Rmid;
    591586    mosaic[i].coords.crval2 = Dmid;
    592     mosaic[i].coords.crpix1 = 0.0;
    593     mosaic[i].coords.crpix2 = 0.0;
    594587    mosaic[i].coords.cdelt1 = 1.0 / 3600.0;
    595588    mosaic[i].coords.cdelt2 = 1.0 / 3600.0;
    596     mosaic[i].coords.pc1_1  = 1.0;
    597     mosaic[i].coords.pc2_2  = 1.0;
    598     mosaic[i].coords.pc1_2  = 0.0;
    599     mosaic[i].coords.pc2_1  = 0.0;
    600589
    601590    mosaic[i].Mcal   = 0.0;
    … …  
    631620      NY = image[m].NY;
    632621      dS += hypot (image[m].coords.cdelt1*image[m].coords.pc1_1, image[m].coords.cdelt1*image[m].coords.pc2_1);
    633       XY_to_RD (&R, &D, 0.0, 0.0, &image[m].coords);
     622
     623      OhanaProjection proj = GetProjection (image[m].coords.ctype);
     624      if ((proj == PROJ_WRP) && !image[m].coords.mosaic) {
     625        XY_to_LM (&R, &D, 0.0, 0.0, &image[m].coords);
     626      } else {
     627        XY_to_RD (&R, &D, 0.0, 0.0, &image[m].coords);
     628      }
    634629      Rmin = MIN (Rmin, R);
    635630      Rmax = MAX (Rmax, R);
    636631      Dmin = MIN (Dmin, D);
    637632      Dmax = MAX (Dmax, D);
    638       XY_to_RD (&R, &D, (double) NX, 0.0, &image[m].coords);
     633
     634      if ((proj == PROJ_WRP) && !image[m].coords.mosaic) {
     635        XY_to_LM (&R, &D, (double) NX, 0.0, &image[m].coords);
     636      } else {
     637        XY_to_RD (&R, &D, (double) NX, 0.0, &image[m].coords);
     638      }
    639639      Rmin = MIN (Rmin, R);
    640640      Rmax = MAX (Rmax, R);
    641641      Dmin = MIN (Dmin, D);
    642642      Dmax = MAX (Dmax, D);
    643       XY_to_RD (&R, &D, (double) NX, (double) NY, &image[m].coords);
     643
     644      if ((proj == PROJ_WRP) && !image[m].coords.mosaic) {
     645        XY_to_LM (&R, &D, (double) NX, (double) NY, &image[m].coords);
     646      } else {
     647        XY_to_RD (&R, &D, (double) NX, (double) NY, &image[m].coords);
     648      }
    644649      Rmin = MIN (Rmin, R);
    645650      Rmax = MAX (Rmax, R);
    646651      Dmin = MIN (Dmin, D);
    647652      Dmax = MAX (Dmax, D);
    648       XY_to_RD (&R, &D, 0.0, (double) NY, &image[m].coords);
     653
     654      if ((proj == PROJ_WRP) && !image[m].coords.mosaic) {
     655        XY_to_LM (&R, &D, 0.0, (double) NY, &image[m].coords);
     656      } else {
     657        XY_to_RD (&R, &D, 0.0, (double) NY, &image[m].coords);
     658      }
    649659      Rmin = MIN (Rmin, R);
    650660      Rmax = MAX (Rmax, R);
    … …  
    652662      Dmax = MAX (Dmax, D);
    653663
     664      // XXX : this probably does not handle mosaics at RA = 0,360 well
     665
    654666      /* we are using mosaic.Mcal, not image.Mcal. reset image.Mcal */
    655667
    … …  
    668680    }
    669681    dS /= MosaicN_Image[i];
    670     strcpy (mosaic[i].coords.ctype, "DEC--TAN");
     682
     683    InitCoords (&mosaic[i].coords, "DEC--TAN");
    671684    mosaic[i].coords.crval1 = Rmin;
    672685    mosaic[i].coords.crval2 = Dmin;
    673     mosaic[i].coords.crpix1 = 0.0;
    674     mosaic[i].coords.crpix2 = 0.0;
    675686    mosaic[i].coords.cdelt1 = dS;
    676687    mosaic[i].coords.cdelt2 = dS;
    677     mosaic[i].coords.pc1_1  = 1.0;
    678     mosaic[i].coords.pc2_2  = 1.0;
    679     mosaic[i].coords.pc1_2  = 0.0;
    680     mosaic[i].coords.pc2_1  = 0.0;
     688
    681689    RD_to_XY (&dX, &dY, Rmax, Dmax, &mosaic[i].coords);
    682690

  • trunk/Ohana/src/relphot/src/assign_images.c

    r37037 r37807  
    8080    if (!strcmp(&image[j].coords.ctype[4], "-DIS")) continue;
    8181
    82     // match the image to its corresponding ASTROMETRIC mosaic (not the same as the Mosaic above)
    83     if (!FindMosaicForImage (image, Nimage, j)) {
    84       if (VERBOSE2) fprintf (stderr, "cannot find mosaic for "OFF_T_FMT"\n", i);
    85       continue;
    86     }
    87    
    8882    // Exclude images with crazy astrometry
    8983    // XXX NOTE : this is gpc1-specific
    … …  
    123117
    124118    if (MOSAIC_ZEROPT) {
    125       // use the coords of the associated mosaic to select
     119      // use the coords of the associated mosaic to select (only for chips; stacks use their own center)
    126120      Mosaic *mosaic = getMosaicForImage (j);
    127       Rc = mosaic->coords.crval1;
    128       Dc = mosaic->coords.crval2;
    129       // NOTE : have defined mosaic Rc,Dc to choose the side of 0,360 on which most of the
    130       // chips are located.  but, for host assignment, we rationalize to 0.0 - 360.0
    131       Rc = ohana_normalize_angle_to_midpoint (Rc, 180.0);
     121      if (mosaic) {
     122        Rc = mosaic->coords.crval1;
     123        Dc = mosaic->coords.crval2;
     124        // NOTE : have defined mosaic Rc,Dc to choose the side of 0,360 on which most of the
     125        // chips are located.  but, for host assignment, we rationalize to 0.0 - 360.0
     126        Rc = ohana_normalize_angle_to_midpoint (Rc, 180.0);
     127      }
    132128    }
    133129
    … …  
    206202}
    207203
    208 int calculate_host_image_bounds (RegionHostTable *regionHosts, double Rmid) {
    209 
    210   int i, n;
    211   off_t j;
    212 
    213   for (i = 0; i < regionHosts->Nhosts; i++) {
    214 
    215     RegionHostInfo *host = &regionHosts->hosts[i];
    216 
    217     // XXX clear the chip match?
    218 
    219     BuildChipMatch (host->image, host->Nimage);
    220 
    221     double Rmin =  720.0;
    222     double Rmax = -360.0;
    223     double Dmin =  +90.0;
    224     double Dmax = -90.0;
    225 
    226     for (j = 0; j < host->Nimage; j++) {
    227 
    228       Image *image = &host->image[j];
    229 
    230       if (!FindMosaicForImage (host->image, host->Nimage, j)) {
    231         fprintf (stderr, "missing astrometry? programming error?\n");
    232         abort ();
    233       }
    234      
    235       // define image corners
    236       for (n = 0; n < 4; n++) {
    237         double Xc, Yc, Rc, Dc;
    238         Xc = Xf[n]*image->NX;
    239         Yc = Yf[n]*image->NY;
    240         XY_to_RD (&Rc, &Dc, Xc, Yc, &image->coords);
    241         Rc = ohana_normalize_angle_to_midpoint (Rc, Rmid);
    242      
    243         Rmin = MIN (Rmin, Rc);
    244         Dmin = MIN (Dmin, Rc);
    245        
    246         Rmax = MAX (Rmax, Rc);
    247         Dmax = MAX (Dmax, Rc);
    248       }
    249 
    250     }
    251 
    252     // RminCat, RmaxCat may extended beyond 0.0 - 360.0
    253     host->RminCat = Rmin;
    254     host->DminCat = Dmin;
    255     host->RmaxCat = Rmax;
    256     host->DmaxCat = Dmax;
    257 
    258     // regionHosts needs to have the full outer boundary
    259     // (so reload_catalogs covers the correct region)
    260     regionHosts->Rmin = MIN(Rmin, regionHosts->Rmin);
    261     regionHosts->Rmax = MAX(Rmax, regionHosts->Rmax);
    262     regionHosts->Dmin = MIN(Dmin, regionHosts->Dmin);
    263     regionHosts->Dmax = MAX(Dmax, regionHosts->Dmax);
    264   }
    265   return TRUE;
    266 }
    267 
    268 
    269204// Rc is in range 0.0 - 360.0, hosts Rmin,Rmax in range 0.0 - 360.0
    270205int find_host_for_coords (RegionHostTable *regionHosts, double Rc, double Dc) {

  • trunk/Ohana/src/relphot/src/load_images.c

    r35759 r37807  
    2727  }
    2828  MARKTIME("read image table: %f sec\n", dtime);
     29
     30  BuildChipMatch (image, Nimage);
     31  MARKTIME("build chip match for %d images: %f sec\n", (int) Nimage, dtime);
     32
     33  char mapfile[DVO_MAX_PATH];
     34  snprintf (mapfile, DVO_MAX_PATH, "%s/AstroMap.fits", CATDIR);
     35  AstromOffsetTable *table = AstromOffsetMapLoad (mapfile, VERBOSE);
     36
     37  // assign images.coords.offsetMap -> table->map[i]
     38  if (table) {
     39    AstromOffsetTableMatchChips (image, Nimage, table);
     40  }
    2941
    3042  // allocate and init an array to identify the images included in the subset

  • trunk/Ohana/src/relphot/src/reload_catalogs.c

    r36630 r37807  
    128128}
    129129
     130int reload_catalog_parallel_group (HostTableGroup *group, SkyList *sky, char *imageFile);
     131
    130132// XXX Image to Image Subset
    131133int reload_catalogs_parallel (SkyList *sky) {
    … …  
    165167    free (BOUNDARY_TREE);
    166168    BOUNDARY_TREE = tmppath;
    167    
    168   }
     169  }
     170
     171  int Ngroups;
     172  HostTableGroup *groups = HostTableGroups (table, &Ngroups);
     173  // split the table into Ngroups, each with a unique set of machines (this avoids overloading the remote host machines)
     174
     175  int i;
     176  for (i = 0; i < Ngroups; i++) {
     177    // update only a group of unique machines at a time
     178    reload_catalog_parallel_group (&groups[i], sky, imageFile);
     179  }
     180  return TRUE;
     181}
     182
     183int reload_catalog_parallel_group (HostTableGroup *group, SkyList *sky, char *imageFile) {
    169184
    170185  int i, j;
    171   for (i = 0; i < table->Nhosts; i++) {
    172 
    173     if (sky->Nregions < table->Nhosts) {
     186  for (i = 0; i < group->Nhosts; i++) {
     187
     188    if (sky->Nregions < group->Nhosts) {
    174189      // do any of the regions want this host?
    175190      int wantThisHost = FALSE;
    176191      for (j = 0; j < sky->Nregions; j++) {
    177         if (HostTableTestHost (sky->regions[j], table->hosts[i].hostID)) {
     192        if (HostTableTestHost (sky->regions[j], group->hosts[i][0].hostID)) {
    178193          wantThisHost = TRUE;
    179194          break;
    … …  
    181196      }
    182197      if (!wantThisHost) {
    183         // fprintf (stderr, "skip host %s\n", table->hosts[i].hostname);
     198        // fprintf (stderr, "skip host %s\n", group->hosts[i][0].hostname);
    184199        continue;
    185200      }
    … …  
    187202
    188203    // ensure that the paths are absolute path names
    189     char *tmppath = abspath (table->hosts[i].pathname, DVO_MAX_PATH);
    190     free (table->hosts[i].pathname);
    191     table->hosts[i].pathname = tmppath;
     204    char *tmppath = abspath (group->hosts[i][0].pathname, DVO_MAX_PATH);
     205    free (group->hosts[i][0].pathname);
     206    group->hosts[i][0].pathname = tmppath;
    192207
    193208    char command[1024];
    194209    snprintf (command, 1024, "relphot_client %s -update-catalogs %s -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -statmode %s -D CAMERA %s -D STAR_TOOFEW %d -minerror %f",
    195               PhotcodeList, imageFile, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, STATMODE, CAMERA, STAR_TOOFEW, MIN_ERROR);
     210              PhotcodeList, imageFile, group->hosts[i][0].hostID, CATDIR, group->hosts[i][0].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, STATMODE, CAMERA, STAR_TOOFEW, MIN_ERROR);
    196211
    197212    // options & configs which affect relphot_client -update-catalogs
    … …  
    238253      // launch the job on the remote machine (no handshake)
    239254      int errorInfo = 0;
    240       int pid = rconnect ("ssh", table->hosts[i].hostname, command, table->hosts[i].stdio, &errorInfo, FALSE);
     255      int pid = rconnect ("ssh", group->hosts[i][0].hostname, command, group->hosts[i][0].stdio, &errorInfo, FALSE);
    241256      if (!pid) {
    242         if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", table->hosts[i].hostname, errorInfo);
     257        if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", group->hosts[i][0].hostname, errorInfo);
    243258        exit (1);
    244259      }
    245       table->hosts[i].pid = pid; // save for future reference
     260      group->hosts[i][0].pid = pid; // save for future reference
    246261    }
    247262  }
    … …  
    252267  }
    253268  if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) {
    254     int status = HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE);
     269    int status = HostTableGroupWaitJobsGetIO (group, __FILE__, __LINE__, VERBOSE);
    255270    if (!status) {
    256271      fprintf (stderr, "at least one remote client job failed to load data, exiting\n");

  • trunk/Ohana/src/relphot/src/relphot.c

    r36630 r37807  
    99  switch (mode) {
    1010    case UPDATE_IMAGES:
     11      // IF CALLED WITH NLOOP == 0, DOES NOT LOAD bcatalog, just loads image table and generates ImageSubset.dat
    1112      relphot_images ();
    1213      exit (0);
    … …  
    1415    case UPDATE_AVERAGES:
    1516      // take the current set of detections and set the mean magnitudes
     17      // DOES NOT LOAD THE IMAGE TABLE
    1618      relphot_objects (0, NULL);
    1719      exit (0);
    … …  
    4345  exit (1);
    4446}
     47
     48/***
     49
     50    I would like to merge the functionality of relphot_images (Nloop == 0), relphot_objects, and apply_offsets
     51   
     52***/

  • trunk/Ohana/src/relphot/src/relphot_objects.c

    r35104 r37807  
    3636    if (!HostTableTestHost(skylist[0].regions[i], hostID)) continue;
    3737
    38     // set up the basic catalog info
    3938    // set up the basic catalog info
    4039    char hostfile[1024];

  • trunk/Ohana/src/relphot/src/relphot_parallel_regions.c

    r37037 r37807  
    4343  MARKTIME("-- load image data: %f sec\n", dtime);
    4444
     45  // save Images.dat using the copied structure
     46  if (UPDATE_CATFORMAT) {
     47    // ensure the db format is updated
     48    db.format = dvo_catalog_catformat (UPDATE_CATFORMAT);
     49  }
     50
    4551  /* assign the images to the different region hosts */
    4652  if (!assign_images (&db, regionHosts)) Shutdown ("error assigning images to region hosts");
    … …  
    4955  /* launch processing on the parallel region hosts */
    5056  if (!launch_region_hosts (regionHosts)) Shutdown ("error launching region hosts");
     57  MARKTIME ("finished relphot -parallel-regions: %f sec total\n", dtime);
    5158
    5259  // retrieve updated image parameters from the remote hosts (also set Image.mcal)

  • trunk/Ohana/src/relphot/src/select_images.c

    r36491 r37807  
    3636
    3737  // the comparison is made in the catalog local projection. below we set crval1,2
    38   tcoords.crpix1 = tcoords.crpix2 = 0.0;
     38  InitCoords (&tcoords, "DEC--TAN");
    3939  tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
    40   tcoords.pc1_1  = tcoords.pc2_2 = 1.0;
    41   tcoords.pc1_2  = tcoords.pc2_1 = 0.0;
    42   strcpy (tcoords.ctype, "DEC--TAN");
    4340
    4441  double RminSkyRegion = region[0].Rmin;
    … …  
    9390
    9491  if (VERBOSE) fprintf (stderr, "finding images\n");
    95   BuildChipMatch (timage, Ntimage);
    96   MARKTIME("build chip match for %d images: %f sec\n", (int) Ntimage, dtime);
    9792
    9893  D_NIMAGE = 6000;
    … …  
    107102     
    108103    if (!(i % 300000)) fprintf (stderr, ".");
     104
     105    // XXX I am just going for force this for the moment:
     106    goto found_it;
    109107
    110108    /* exclude images by photcode */
    … …  
    123121    }
    124122   
    125     // this adds 1.3 sec for 3M images
    126     if (!FindMosaicForImage (timage, Ntimage, i)) {
    127       if (VERBOSE2) fprintf (stderr, "cannot find mosaic for "OFF_T_FMT"\n", i);
    128       continue;
    129     }
    130 
    131123    /* define image corners - note the DIS images (mosaic phu) are special */
    132124    if (!strcmp(&timage[i].coords.ctype[4], "-DIS")) {
    … …  
    405397
    406398  // the comparison is made in the catalog local projection. below we set crval1,2
    407   tcoords.crpix1 = tcoords.crpix2 = 0.0;
     399  InitCoords (&tcoords, "DEC--TAN");
    408400  tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
    409   tcoords.pc1_1  = tcoords.pc2_2 = 1.0;
    410   tcoords.pc1_2  = tcoords.pc2_1 = 0.0;
    411 
    412   tcoords.Npolyterms = 0;
    413   memset (tcoords.polyterms, 0, 14*sizeof(float));
    414   strcpy (tcoords.ctype, "DEC--TAN");
    415401
    416402  /* compare with each region file */

  • trunk/Ohana/src/relphot/src/setMrelCatalog.c

    r37740 r37807  
    3636}
    3737
    38 int setMrelAverageForcedWarp (off_t measureOffset, int cat, int pass, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t *found);
    39 
    4038int setMrel_catalog_alt (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) {
    4139
    … …  
    6058    Average *average = catalog[Nc].average ? &catalog[Nc].average[j] : NULL;
    6159    Measure *measure = catalog[Nc].measure ? &catalog[Nc].measure[m] : NULL;
    62     setMrelAverageExposure (m, Nc, pass, flatcorr, results, average, &catalog[Nc].averageT[j], &catalog[Nc].secfilt[j*Nsecfilt], measure, &catalog[Nc].measureT[m], &catalog[Nc].found_t[Nsecfilt*j]);
     60    SecFilt *secfilt = catalog[Nc].secfilt ? &catalog[Nc].secfilt[j*Nsecfilt] : NULL;
     61    setMrelAverageExposure (m, Nc, pass, flatcorr, results, average, &catalog[Nc].averageT[j], secfilt, measure, &catalog[Nc].measureT[m], &catalog[Nc].found_t[Nsecfilt*j]);
    6362
    6463    // only apply Stack operation on setMrelFinal in first pass
    6564    if (isSetMrelFinal && (pass == 0)) {
    66       setMrelAverageStack (m, Nc, flatcorr, results, &catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt], &catalog[Nc].measure[m]);
     65      setMrelAverageStack (m, Nc, flatcorr, results, average, secfilt, &catalog[Nc].measure[m], catalog[Nc].Nmeasure);
    6766      setGlobalObjStats (&catalog[Nc].average[j], &catalog[Nc].measure[m]);
    6867    }
    … …  
    7069    // only measure force-warp mean values if issetMrelFinal (make it optional?)
    7170    if (isSetMrelFinal) {
    72       setMrelAverageForcedWarp (m, Nc, pass, results, &catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt], &catalog[Nc].measure[m], &catalog[Nc].foundWarp_t[Nsecfilt*j]);
    73       setGlobalObjStats (&catalog[Nc].average[j], &catalog[Nc].measure[m]);
     71      setMrelAverageForcedWarp (m, Nc, pass, results, average, secfilt, measure, &catalog[Nc].foundWarp_t[Nsecfilt*j], catalog[Nc].Nmeasure);
     72      setGlobalObjStats (average, measure);
    7473    }
    7574  }
    … …  
    456455// 2) select the BEST detections per filter (regardless of PRIMARY)
    457456// 3) apply the zero point and AB->Jy transformations
    458 int setMrelAverageStack (off_t measureOffset, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure) {
     457int setMrelAverageStack (off_t measureOffset, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t Nmeasure) {
    459458
    460459  off_t k;
    … …  
    496495    int Nstack = 0; // number for this photcode
    497496    int NstackDet = 0; // number for this photcode
     497
     498    secfilt[Nsec].stackBestOff = -1;
     499    secfilt[Nsec].stackPrmryOff = -1;
    498500
    499501    off_t meas = measureOffset;
    … …  
    526528        measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY;
    527529        secfilt[Nsec].stackPrmryOff = meas;
     530        myAssert (secfilt[Nsec].stackPrmryOff <= Nmeasure, "stackPrmryOff out of range");
    528531      }
    529532
    … …  
    637640
    638641    secfilt[Nsec].stackBestOff = k + measureOffset;
     642    myAssert (secfilt[Nsec].stackBestOff <= Nmeasure, "stackBestOff out of range");
    639643
    640644    secfilt[Nsec].Nstack    = Nstack;
    … …  
    669673// somewhat simplified relative to chip-photometry:
    670674// * no grid, no mosaic, no 2MASS, no SYNTH, no Ubercal, no flatcorr
    671 int setMrelAverageForcedWarp (off_t measureOffset, int cat, int pass, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t *found) {
     675int setMrelAverageForcedWarp (off_t measureOffset, int cat, int pass, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t *found, off_t Nmeasure) {
    672676
    673677  off_t k;
    … …  
    713717    // if we cannot make that association, skip this set of warps
    714718    int stkTessID, stkProjID, stkSkycellID;
     719
     720    myAssert (secfilt[Nsec].stackBestOff <= Nmeasure, "stackBestOff out of range");
    715721    if (!FindImageSkycellID (secfilt[Nsec].stackBestOff, cat, &stkTessID, &stkProjID, &stkSkycellID)) continue;
    716722
Note: See TracChangeset for help on using the changeset viewer.