IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 41625


Ignore:
Timestamp:
May 31, 2021, 11:45:18 AM (5 years ago)
Author:
eugene
Message:

major revision of the flagging and resets. cleanup of old cruft and deprecated flag values

Location:
branches/eam_branches/relphot.20210521
Files:
15 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/relphot.20210521/include/relphot.h

    r41624 r41625  
    88
    99# define NBOOTSTRAP 100
     10
     11// MEAS_BAD is used in ImageOps.c, MosaicOps.c, TGroupOps.c to skip specific measurements.
     12// ID_MEAS_AREA is raised for measurements outside the user-specified region of the chips
     13// ID_MEAS_NOCAL is raised for measurements outside the valid time range or for inactive photcodes
     14// Note: these bits are only raised on the temporary MeasureTiny elements and are not saved
     15// to the database files.
     16# define MEAS_BAD (ID_MEAS_AREA | ID_MEAS_NOCAL)
    1017
    1118# ifndef MAX_INT
    … …  
    453460char   MOSAICNAME[256];
    454461char   STATMODE[32];
    455 int    STAR_BAD;
    456 int    MEAS_BAD;
    457462int    STAR_TOOFEW;
    458463int    GRID_TOOFEW;
    459464int    IMAGE_TOOFEW;
    460465double IMAGE_GOOD_FRACTION;
    461 int    IMAGE_BAD;
    462466int    CALIBRATE_STACKS_AND_WARPS;
    463467
    … …  
    848852int GetZptIteration (void);
    849853
    850 void SetZeroPointModes (void);
     854void SetZeroPointModes (Catalog *catalog, int Ncatalog);
    851855int UseStandardOLS (ZptFitModeType mode);
    852856int GetActivePhotcodeIndex (int photcode);

  • branches/eam_branches/relphot.20210521/src/GridOps.c

    r41620 r41625  
    192192      for (int k = 0; k < catalog[nc].averageT[na].Nmeasure; k++, nm++) {
    193193
    194         // XXX I need to skip bad measurements, but I am not certain this
    195         // flag is set correctly.  see note below
     194        // skip measurements marked by AREA or TIME
    196195        if (catalog[nc].measureT[nm].dbFlags & MEAS_BAD) continue;
    197196

  • branches/eam_branches/relphot.20210521/src/ImageOps.c

    r41608 r41625  
    385385  if (i == -1) return NAN;
    386386
    387   if (image[i].flags & IMAGE_BAD) return NAN; 
     387  // allow context dependence?: in relphot_images, we want to exclude detections from poor images
     388  // but in the final pass, we want to allow even bad detections
     389  if (image[i].flags & (ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW)) return NAN; 
    388390
    389391  switch (class) {
    … …  
    409411  if (i == -1) return (1000);
    410412
    411   if (image[i].flags & IMAGE_BAD) return (1000); 
    412413  distance = image[i].ubercalDist; // was dummy3 in structure
    413414  return (distance);
    … …  
    689690    if (mark) {
    690691      image[i].flags |= ID_IMAGE_PHOTOM_FEW;
     692      image[i].McalPSF    = 0.0;
     693      image[i].McalAPER   = 0.0;
     694      image[i].McalChiSq  = NAN;
    691695      image[i].dMcal      = NAN;
    692       image[i].McalChiSq  = NAN;
     696      image[i].dMagSys    = NAN;
    693697      image[i].nFitPhotom = 0;
    694698      Nfew ++;
    … …  
    894898  // measure stats for Mcal and dMcal
    895899  for (i = N = 0; i < Nimage; i++) {
    896     if (image[i].flags & IMAGE_BAD) continue;
     900    if (image[i].flags & (ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
    897901
    898902    if (FREEZE_IMAGES && isMosaicChip(image[i].photcode)) continue;
    … …  
    922926
    923927    mark = FALSE;
    924     image[i].flags &= ~ID_IMAGE_PHOTOM_POOR;
    925     mark = (image[i].dMcal > MaxScatter) || (fabs(image[i].McalPSF - MedOffset) > MaxOffset);
     928
     929    if (image[i].flags & ID_IMAGE_PHOTOM_FEW) {
     930      mark = TRUE;
     931    }
     932    if (image[i].dMcal > MaxScatter) {
     933      mark = TRUE;
     934    }
     935    if (fabs(image[i].McalPSF - MedOffset) > MaxOffset) {
     936      mark = TRUE;
     937    }
     938
    926939    if (mark) {
    927940      Nmark ++;
    … …  
    10891102  n = 0;
    10901103  for (i = 0; i < Nimage; i++) {
    1091     if (image[i].flags & IMAGE_BAD) continue;
     1104    if (image[i].flags & (ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
    10921105
    10931106    if (FREEZE_IMAGES && isMosaicChip(image[i].photcode)) continue;
    … …  
    11381151  for (i = 0; i < Nimage; i++) {
    11391152
    1140     if (image[i].flags & IMAGE_BAD) continue;
     1153    if (image[i].flags & (ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
    11411154
    11421155    if (FREEZE_IMAGES && isMosaicChip(image[i].photcode)) continue;
    … …  
    11711184  for (i = 0; i < Nimage; i++) {
    11721185
    1173     if (image[i].flags & IMAGE_BAD) continue;
     1186    if (image[i].flags & (ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
    11741187
    11751188    if (FREEZE_IMAGES && isMosaicChip(image[i].photcode)) continue;
    … …  
    12041217  for (i = 0; i < Nimage; i++) {
    12051218
    1206     if (image[i].flags & IMAGE_BAD) continue;
     1219    if (image[i].flags & (ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
    12071220
    12081221    if (FREEZE_IMAGES && isMosaicChip(image[i].photcode)) continue;

  • branches/eam_branches/relphot.20210521/src/MagResidSave.c

    r41608 r41625  
    215215    off_t c = ImageToCatalog[j];
    216216     
    217     // XXX drop this?
    218217    if (catalog[c].measureT[m].dbFlags & MEAS_BAD) continue;
    219218

  • branches/eam_branches/relphot.20210521/src/MosaicOps.c

    r41608 r41625  
    780780  if (i == -1) return (0.0);
    781781
    782   if (mosaic[i].flags & IMAGE_BAD) return (NAN); 
     782  // if the mosaic cannot be calibrated (too few measurements, skip it)
     783  // XXX: if the mosaic is bad (few) but the chips are OK, I can used them.
     784  // XXX: if (mosaic[i].flags & ID_IMAGE_PHOTOM_FEW) return (NAN); 
    783785  value = mosaic[i].McalPSF;
    784786  return (value);
    … …  
    13301332}
    13311333
    1332 // When we rationalize the images/mosaics, we are driving the negative cloud images back
    1333 // to 0.0.  At the same time, we make a guess to the effective impact on all other images,
    1334 // driven by the coupling of common stars.
    1335 // XXX this function is no longer used because we force significantly negative clouds to 0.0
    1336 int rationalize_mosaics (Catalog *catalog, int Ncatalog) {
    1337 
    1338   double *imageOffset, **starOffset;
    1339   int **starNcount, *seclist;
    1340   int **Slist, *NSlist, *NSLIST;
    1341   int i, j, k, m, nMos, Ns, found;
    1342 
    1343   off_t Nimage;
    1344   Image *image;
    1345 
    1346   // set a test value for now
    1347   float CLOUD_TOLERANCE = 0.025;
    1348 
    1349   if (!MOSAIC_ZEROPT) return (FALSE);
    1350   if (FREEZE_MOSAICS) return (FALSE);
    1351 
    1352   image = getimages (&Nimage, NULL);
    1353 
    1354   ALLOCATE (imageOffset, double, Nmosaic);
    1355   ALLOCATE ( Slist, int *, Nmosaic); // array of calibrated star indexes on this mosaic
    1356   ALLOCATE (NSlist, int,   Nmosaic); // number of stars on mosaic
    1357   ALLOCATE (NSLIST, int,   Nmosaic); // number of Slist entries allocated
    1358   memset (Slist, 0, Nmosaic*sizeof(int *));
    1359 
    1360   // find the images / mosaics with negative clouds and save their offset
    1361   for (i = 0; i < Nmosaic; i++) {
    1362  
    1363     NSlist[i] =   0;
    1364     NSLIST[i] = 100;
    1365     ALLOCATE (Slist[i], int, NSLIST[i]);
    1366 
    1367     imageOffset[i] = 0.0;
    1368 
    1369     if (VERBOSE2 && (fabs(mosaic[i].McalPSF) < CLOUD_TOLERANCE)) {
    1370       fprintf (stderr, "cloud-free: %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].McalPSF);
    1371     }
    1372     if (VERBOSE2 && (mosaic[i].McalPSF < -CLOUD_TOLERANCE)) {
    1373       imageOffset[i] = -mosaic[i].McalPSF;
    1374       // NOTE the negative sign: down below, we are going to add in the negative of McalPSF
    1375       // to this image, and the propagated mean values for other images
    1376       fprintf (stderr, "anti-clouds: %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].McalPSF);
    1377     }
    1378     if (VERBOSE2 && (mosaic[i].McalPSF > CLOUD_TOLERANCE)) {
    1379       fprintf (stderr, "cloudy    : %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].McalPSF);
    1380     }
    1381   }
    1382 
    1383   int Nsecfilt = GetPhotcodeNsecfilt ();
    1384   ALLOCATE (seclist, int, Nphotcodes);
    1385   for (Ns = 0; Ns < Nphotcodes; Ns ++) {
    1386     int thisCode = photcodes[Ns][0].code;
    1387     seclist[Ns] = GetPhotcodeNsec(thisCode);
    1388   }
    1389 
    1390   // allocate an array for star offsets
    1391   int Nstars = 0;
    1392   int NSTARS = 1000;
    1393   ALLOCATE (starOffset, double *, NSTARS);
    1394   ALLOCATE (starNcount, int *,    NSTARS);
    1395   memset (starOffset, 0, NSTARS*sizeof(double *));
    1396   memset (starNcount, 0, NSTARS*sizeof(int *));
    1397 
    1398   // find the mean offset for each star
    1399   for (i = 0; i < Ncatalog; i++) {
    1400     for (j = 0; j < catalog[i].Naverage; j++) {
    1401       ALLOCATE (starOffset[Nstars], double, Nphotcodes);
    1402       ALLOCATE (starNcount[Nstars], int,    Nphotcodes);
    1403       memset (starOffset[Nstars], 0, Nphotcodes*sizeof(double));
    1404       memset (starNcount[Nstars], 0, Nphotcodes*sizeof(int));
    1405      
    1406       m = catalog[i].averageT[j].measureOffset;
    1407 
    1408       // determine the mosaic for each measurement
    1409       for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
    1410 
    1411         // skip unused measurements
    1412         if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue;
    1413 
    1414         // skip unused measurements
    1415 
    1416         // XXX replace with:
    1417         // Ns = GetActivePhotcodeIndex (photcode);
    1418         // if (Ns < 0) continue;
    1419 
    1420         int Nsec;
    1421         int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
    1422         found = FALSE;
    1423         for (Ns = 0; !found && (Ns < Nphotcodes); Ns ++) {
    1424           if (ecode == photcodes[Ns][0].code) {
    1425             found = TRUE;
    1426             break;
    1427           }
    1428         }
    1429         if (!found) continue;
    1430         Nsec = seclist[Ns];
    1431 
    1432         // bad stars for this secfilt
    1433         if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; 
    1434 
    1435         // skip REF measurements (not tied to an image)
    1436         if (getImageEntry (m, i) < 0) continue;
    1437 
    1438         // find the source of this measurement (skip unassigned measurements)
    1439         nMos = MeasureToMosaic[i][m];
    1440         if (nMos == -1) continue;
    1441 
    1442         // XXX how can this not be true? 
    1443         if (mosaic[nMos].photcode != ecode) {
    1444           fprintf (stderr, "*");
    1445         }
    1446 
    1447         assert (Ns >= 0);
    1448         assert (Ns < Nphotcodes);
    1449 
    1450         // accumulate the offsets from the negative cloud images (others have 0.0 value)
    1451         starOffset[Nstars][Ns] += imageOffset[nMos];
    1452         starNcount[Nstars][Ns] ++;
    1453        
    1454         // record the mosaic->star reference
    1455         Slist[nMos][NSlist[nMos]] = Nstars;
    1456         NSlist[nMos] ++;
    1457         if (NSlist[nMos] == NSLIST[nMos]) {
    1458           NSLIST[nMos] += 100;
    1459           REALLOCATE (Slist[nMos], int, NSLIST[nMos]);
    1460         }         
    1461       }
    1462       Nstars ++;
    1463       if (Nstars == NSTARS) {
    1464         NSTARS += 1000;
    1465         REALLOCATE (starOffset, double *, NSTARS);
    1466         REALLOCATE (starNcount, int *,    NSTARS);
    1467         memset (&starOffset[NSTARS-1000], 0, 1000*sizeof(double *));
    1468         memset (&starNcount[NSTARS-1000], 0, 1000*sizeof(int *));
    1469       }
    1470     }
    1471   }
    1472 
    1473   // find the mean offset of the images without negative clouds
    1474   for (i = 0; i < Nmosaic; i++) {
    1475 
    1476     found = FALSE;
    1477     for (Ns = 0; !found && (Ns < Nphotcodes); Ns ++) {
    1478       if (mosaic[i].photcode == photcodes[Ns][0].code) {
    1479         found = TRUE;
    1480         break;
    1481       }
    1482     }
    1483     if (!found) {
    1484       fprintf (stderr, "invalid photcode for mosaic?\n");
    1485       abort();
    1486     }
    1487 
    1488     // a negative cloud image (cloud: McalPSF > 0; anti-clouds: McalPSF < 0; imageOffset = -McalPSF)
    1489     if (imageOffset[i] > 0.0) continue;
    1490 
    1491     // we need to actually have cross-references to count
    1492     if (NSlist[i] < 2) continue;
    1493 
    1494     float dM = 0.0;
    1495     for (j = 0; j < NSlist[i]; j++) {
    1496       k = Slist[i][j];
    1497       if (starNcount[k][Ns] > 1) {
    1498         dM += (starOffset[k][Ns] / starNcount[k][Ns]);
    1499       }
    1500     }
    1501     imageOffset[i] = dM / NSlist[i];
    1502     // fprintf (stderr, "adjust image: %s : (%f %d) : %f\n", image[MosaicToImage[i][0]].name, dM, NSlist[i], imageOffset[i]);
    1503   }
    1504 
    1505   // for (i = 0; i < Nmosaic; i++) {
    1506   //   fprintf (stderr, "correction: %s : %f\n", image[MosaicToImage[i][0]].name, imageOffset[i]);
    1507   // }
    1508 
    1509   // apply all offset values to the mosaics
    1510   // find the images / mosaics with negative clouds and save their offset
    1511   for (i = 0; i < Nmosaic; i++) {
    1512     mosaic[i].McalPSF += imageOffset[i];
    1513   }
    1514 
    1515   for (i = 0; i < Nstars; i++) {
    1516     free (starOffset[i]);
    1517     free (starNcount[i]);
    1518   }
    1519   free (starOffset);
    1520   free (starNcount);
    1521 
    1522   free (seclist);
    1523   free (NSlist);
    1524   free (NSLIST);
    1525   for (i = 0; i < Nmosaic; i++){
    1526     free (Slist[i]);
    1527   }
    1528   free (Slist);
    1529   free (imageOffset);
    1530   return (TRUE);
    1531 }
    1532 
    15331334StatType statsMosaicM (Catalog *catalog) {
    15341335  OHANA_UNUSED_PARAM(catalog);
    … …  
    15471348  n = 0;
    15481349  for (i = 0; i < Nmosaic; i++) {
    1549     if (mosaic[i].flags & IMAGE_BAD) continue;
     1350    if (mosaic[i].flags & (ID_IMAGE_MOSAIC_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
    15501351    if (mosaic[i].skipCal) continue;
    15511352    if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
    … …  
    15801381  n = 0;
    15811382  for (i = 0; i < Nmosaic; i++) {
    1582     if (mosaic[i].flags & IMAGE_BAD) continue;
     1383    if (mosaic[i].flags & (ID_IMAGE_MOSAIC_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
    15831384    if (mosaic[i].skipCal) continue;
    15841385    if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
    … …  
    16131414  n = 0;
    16141415  for (i = 0; i < Nmosaic; i++) {
    1615     if (mosaic[i].flags & IMAGE_BAD)  continue;
     1416    if (mosaic[i].flags & (ID_IMAGE_MOSAIC_POOR | ID_IMAGE_PHOTOM_FEW))  continue;
    16161417    if (mosaic[i].skipCal) continue;
    16171418    if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
    … …  
    16631464  n = 0;
    16641465  for (i = 0; i < Nmosaic; i++) {
    1665     if (mosaic[i].flags & IMAGE_BAD) continue;
     1466    if (mosaic[i].flags & (ID_IMAGE_MOSAIC_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
    16661467    if (mosaic[i].skipCal) continue;
    16671468    if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
    … …  
    17951596      c = MosaicToCatalog[i][j];
    17961597     
    1797       if (catalog[c].measureT[m].dbFlags & (ID_MEAS_AREA | ID_MEAS_NOCAL)) continue;
     1598      if (catalog[c].measureT[m].dbFlags & MEAS_BAD) continue;
    17981599
    17991600      // ave = catalog[c].measureT[m].averef;

  • branches/eam_branches/relphot.20210521/src/ResetOps.c

    r41624 r41625  
    7575    }
    7676  }
     77  if (RESET_FLATCORR) {
     78    measure->Mflat = 0.0;
     79  }
    7780  measure->dbFlags &= ~photomBitsMeasure;
    7881}
    … …  
    8891  for (off_t i = 0; i < Nsubset; i++) {
    8992    // reset these bits regardless (we will re-determine)
    90     subset[i].flags &= ~(ID_IMAGE_PHOTOM_NOCAL | ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_SKIP | ID_IMAGE_PHOTOM_FEW);
     93    subset[i].flags &= ~(ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW);
    9194    if (RESET) {
    9295      if (RESET_ZEROPTS) {
    … …  
    117120
    118121  for (off_t i = 0; i < catalog->Naverage; i++) {
    119     catalog[0].average[i].psfQF     = NAN;      // force recalculation in setMrelCatalog
    120     catalog[0].average[i].psfQFperf = NAN;      // force recalculation in setMrelCatalog
    121     catalog[0].average[i].stargal   = NAN;      // force recalculation in setMrelCatalog
    122     catalog[0].average[i].photFlagsUpper = 0;   // reset (will be re-calculated)
    123     catalog[0].average[i].photFlagsLower = 0;   // reset (will be re-calculated)
    124     catalog[0].average[i].NwarpOK        = 0;   // reset (will be re-calculated)
     122    catalog->average[i].psfQF     = NAN;        // force recalculation in setMrelCatalog
     123    catalog->average[i].psfQFperf = NAN;        // force recalculation in setMrelCatalog
     124    catalog->average[i].stargal   = NAN;        // force recalculation in setMrelCatalog
     125    catalog->average[i].photFlagsUpper = 0;     // reset (will be re-calculated)
     126    catalog->average[i].photFlagsLower = 0;     // reset (will be re-calculated)
     127    catalog->average[i].NwarpOK        = 0;     // reset (will be re-calculated)
     128
     129    // RESET_FLATCORR independent of RESET?
    125130
    126131    if (!RESET) continue;
    … …  
    128133    catalog->average[i].flags    &= ~photomBitsAverage; // reset all photometry bits (but not astrom)
    129134
    130     for (int Ns = 0; Ns < Nsecfilt; Ns++) {
     135    // Reset secfilt values only for the active photcodes
     136    for (int Ns = 0; Ns < Nphotcodes; Ns++) {
     137     
     138      int thisCode = photcodes[Ns][0].code;
     139      int Nsec = GetPhotcodeNsec(thisCode);
     140      off_t N = Nsecfilt*i+Nsec;
     141     
     142      dvo_secfilt_init (&catalog->secfilt[N], SECFILT_RESET_ALL);
     143    }
     144
     145    // Reset measure values only for the active photcodes
     146   
     147    off_t m = catalog->average[i].measureOffset;
     148   
     149    for (off_t j = 0; j < catalog->average[i].Nmeasure; j++, m++) {
     150     
     151      /* select measurements for active photcodes */
     152      PhotCode *code = GetPhotcodebyCode (catalog->measure[m].photcode);
     153      if (!code) continue;
     154      // photcode exists, but is it one of the active ones?
     155      int found = FALSE;
     156      for (int Ns = 0; !found && (Ns < Nphotcodes); Ns++) {
     157        if (code->equiv != photcodes[Ns][0].code) continue;
     158        found = TRUE;
     159      }
     160      if (!found) continue;
     161
     162      // reset the Mflat value if requested
     163      if (RESET_FLATCORR) {
     164        catalog->measure[m].Mflat = 0.0;
     165      }
     166     
     167      /* select measurements by time */
     168      if (TimeSelect) {
     169        if (catalog->measure[m].t < TSTART) continue;
     170        if (catalog->measure[m].t > TSTOP) continue;
     171      }
    131172       
    132       off_t N = Nsecfilt*i+Ns;
    133       dvo_secfilt_init (&catalog->secfilt[N], SECFILT_RESET_ALL);
     173      // only reset Mcal for measures with a matching image
     174      // do not reset Mcal for ubercal images unless explicitly requested
    134175
    135       off_t m = catalog->average[i].measureOffset;
    136 
    137       for (off_t j = 0; j < catalog->average[i].Nmeasure; j++, m++) {
    138        
    139         /* select measurements by time */
    140         if (TimeSelect) {
    141           if (catalog->measure[m].t < TSTART) continue;
    142           if (catalog->measure[m].t > TSTOP) continue;
     176      if (catalog->measure[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
     177        if (!KEEP_UBERCAL) {
     178          catalog->measure[m].McalPSF  = 0.0;
     179          catalog->measure[m].McalAPER = 0.0;
     180          catalog->measure[m].dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL;
     181        }
     182      } else {
     183        if (RESET_ZEROPTS && (getImageEntry (m, 0) >= 0)) {
     184          catalog->measure[m].McalPSF  = 0.0;
     185          catalog->measure[m].McalAPER = 0.0;
    143186        }
    144        
    145         // only reset Mcal for measures with a matching image
    146         // do not reset Mcal for ubercal images unless explicitly requested
    147 
    148         if (catalog->measure[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
    149           if (!KEEP_UBERCAL) {
    150             catalog->measure[m].McalPSF  = 0.0;
    151             catalog->measure[m].McalAPER = 0.0;
    152             catalog->measure[m].dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL;
    153           }
    154         } else {
    155           if (RESET_ZEROPTS && (getImageEntry (m, 0) >= 0)) {
    156             catalog->measure[m].McalPSF  = 0.0;
    157             catalog->measure[m].McalAPER = 0.0;
    158           }
    159         }
    160         catalog->measure[m].dbFlags &= ~photomBitsMeasure;
    161187      }
     188      catalog->measure[m].dbFlags &= ~photomBitsMeasure;
    162189    }
    163190  }
    164191}
    165192
     193

  • branches/eam_branches/relphot.20210521/src/StarOps.c

    r41624 r41625  
    11# include "relphot.h"
     2
     3# define STAR_BAD (ID_SECF_STAR_FEW | ID_SECF_STAR_POOR)
    24
    35enum {THREAD_RUN, THREAD_DONE};
    … …  
    290292
    291293// setMrel and setMrelOutput are extremely similar, but have slightly different implications:
     294// * both call setMrelCatalog, but setMrel set isSetMrelFinal to FALSE
    292295// * setMrel only uses the internal Tiny structures
    293 // * setMrel skips stars for which there are too few good measurements (setMrelOutput assigns an average mag with only 1)
     296// * setMrel skips stars for which there are too few good measurements
     297// -> (setMrelOutput assigns an average mag for all objects, even those with only 1 measurement)
    294298// * setMrelOutput sets psfQF, psfQFperf, extNsigma, rank info, etc
    295299// * setMrelOutput sets average Map
    296300// * setMrelOutput updates 2MASS average flags
    297301// * setMrelOutput updates average EXT flags (PS1 and 2MASS)
    298 // * setMrelOutput sets average stack & warp photometry (if requested)
     302// * setMrelOutput sets average stack & warp photometry (if requested), setMrel always skips those
     303// * setMrelOutput is NOT threaded it reads/writes catalogs one at a time
    299304
    300305int setMrel (Catalog *catalog, int Ncatalog) {
    … …  
    325330}
    326331
    327 // setMrelOutput is NOT threaded because we read/write catalogs one at a time (for now)
    328332int setMrelOutput (Catalog *catalog, int Ncatalog) {
    329333
    … …  
    443447}
    444448
    445 # define SKIP_THIS_MEAS(REASON) {                               \
    446     catalog[Nc].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM;     \
    447     if (catalog[Nc].measure) {                                  \
    448       catalog[Nc].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM;   \
    449     }                                                           \
    450     results->REASON ++;                                         \
    451     continue; }
    452 
    453 
    454449int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure) {
    455450
    … …  
    471466}
    472467
    473 /* set measure.Mcal for all measures except ID_MEAS_NOCAL and ID_IMAGE_PHOTOM_NOCAL */
     468/* set measure.Mcal for all measures except ID_MEAS_NOCAL */
    474469int setMcalOutput (Catalog *catalog, int Ncatalog) {
    475470
    … …  
    477472  off_t j, k, m;
    478473
    479   MEAS_BAD = ID_MEAS_NOCAL;
    480   IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL;
    481 
    482474  for (i = 0; i < Ncatalog; i++) {
    483475    for (j = 0; j < catalog[i].Naverage; j++) {
    … …  
    485477      m = catalog[i].averageT[j].measureOffset;
    486478      for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
    487         if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue;
     479        if (catalog[i].measureT[m].dbFlags & ID_MEAS_NOCAL) continue;
    488480        float McalPSF  = getMcal  (m, i, MAG_CLASS_PSF);
    489481        float McalAPER = getMcal  (m, i, MAG_CLASS_APER);
    490482        if (isnan(McalPSF)) continue;
    491         float Mmos  = getMmos  (m, i);
    492         if (isnan(Mmos)) continue;
    493         float Mgrp  = getMgrp  (m, i, catalog[i].measureT[m].airmass, NULL);
    494         if (isnan(Mgrp)) continue;
    495         float Mgrid = getMgridTiny (&catalog[i].measureT[m]);
    496         if (isnan(Mgrid)) continue;
     483
     484        // Mmos, Mgrp have already been transfered to Mcal and should be zero here
     485        // float Mmos  = getMmos  (m, i);
     486        // if (isnan(Mmos)) continue;
     487        // float Mgrp  = getMgrp  (m, i, catalog[i].measureT[m].airmass, NULL);
     488        // if (isnan(Mgrp)) continue;
     489
     490        // Mgrid is applied to Mflat and should be ignored here
     491        // float Mgrid = getMgridTiny (&catalog[i].measureT[m]);
     492        // if (isnan(Mgrid)) continue;
    497493
    498494        // Note that this operation is setting measure->McalAPER to image->McalAPER only
    … …  
    507503
    508504        // set the output calibration
    509         // float Moff = Mmos + Mgrp + Mgrid;
    510         float Moff = Mmos + Mgrid; // Mgrp has already been applied to the image.Mcal
    511         catalog[i].measure[m].McalPSF  = McalPSF + Moff;
    512         catalog[i].measure[m].McalAPER = useStackAper ? McalAPER + Moff : McalPSF + Moff;
     505        catalog[i].measure[m].McalPSF  = McalPSF;
     506        catalog[i].measure[m].McalAPER = useStackAper ? McalAPER : McalPSF;
    513507
    514508        if (catalog[i].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
    … …  
    612606    for (i = Ntot = 0; i < Ncatalog; i++) {
    613607      for (j = 0; j < catalog[i].Naverage; j++) {
    614         if ( catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD ) continue;
     608        if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;
    615609        float Mchisq = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq;
    616610        if (isnan(Mchisq)) continue;
    … …  
    621615      }
    622616    }
    623  
     617
     618    // XXX the limits for MaxChisq and MaxScatter should be
     619    // user-defined.  This uses
    624620    liststats (xlist, dlist, NULL, Ntot, &stats);
    625     MaxChisq = MAX (STAR_CHISQ, 2*stats.median);
     621    float ChiSqUpper90 = stats.Upper90;
     622    MaxChisq = MAX (STAR_CHISQ, ChiSqUpper90);
    626623
    627624    liststats (slist, dlist, NULL, Ntot, &stats);
    628     MaxScatter = MAX (STAR_SCATTER, 2*stats.median);
     625    float ScatterUpper90 = stats.Upper90;
     626    MaxScatter = MAX (STAR_SCATTER, ScatterUpper90);
    629627    fprintf (stderr, "Max Scatter: %f, Max Chisq: %f\n", MaxScatter, MaxChisq);
    630628
    … …  
    636634        mark = (dM > MaxScatter) || (isnan(Mchisq)) || (Mchisq > MaxChisq);
    637635        if (mark) {
    638           catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_OBJ_POOR;
     636          catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_STAR_POOR;
    639637          Ndel ++;
    640638          if (dM > MaxScatter)   { Nscat ++; }
    … …  
    642640          if (Mchisq > MaxChisq) { Nchi ++; }
    643641        } else {
    644           catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_OBJ_POOR;
     642          catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_SECF_STAR_POOR;
    645643        }
    646644        Nave ++;
    … …  
    652650  free (slist);
    653651  free (dlist);
    654 }
    655 
    656 // clean_measures examines the stats for a single star.  It measures the INNER 50% mean
    657 // and sigma, it then re-measures the mean and sigma using all stars with NSIGMA_CLIP (3)
    658 // sigma of the INNER 50% mean.  it then flags any measurements which are more than
    659 // NSIGMA_REJECT (5) sigma of the mean
    660 
    661 // this function only operations on the PSF magnitudes
    662 
    663 // XXX EAM clean_measures is no longer used -- we rely on IRLS to down-weight outliers
    664 
    665 # define NSIGMA_CLIP 3.0
    666 # define NSIGMA_REJECT 5.0
    667 void clean_measures (Catalog *catalog, int Ncatalog, int final) {
    668 
    669   off_t j, k, m, Nmax, Ndel, Nave;
    670   int i, N, image_bad, TOOFEW;
    671   off_t *ilist;
    672   double *tlist, *list, *dlist;
    673   float Msys, Mcal, Mmos, Mgrp, Mgrid;
    674   int Ncal, Nmos, Ngrp, Ngrid, Nfew;
    675 
    676   int Nsecfilt = GetPhotcodeNsecfilt ();
    677 
    678   if (VERBOSE) fprintf (stderr, "marking poor measures\n");
    679   /* Nmeasure is now different, need to reallocate */
    680   Nmax = 0;
    681   for (i = 0; i < Ncatalog; i++) {
    682     for (j = 0; j < catalog[i].Naverage; j++) {
    683       Nmax = MAX (Nmax, catalog[i].averageT[j].Nmeasure);
    684     }
    685   }
    686   ALLOCATE (list, double, MAX (1, Nmax));
    687   ALLOCATE (dlist, double, MAX (1, Nmax));
    688   ALLOCATE (ilist, off_t, Nmax);
    689   ALLOCATE (tlist, double, Nmax);
    690  
    691   /* it makes no sense to mark 3-sigma outliers with <5 measurements */
    692   TOOFEW = MAX (5, STAR_TOOFEW);
    693 
    694   // stats structures for inner and full stats
    695   StatType instats, stats;
    696   liststats_setmode (&instats, "INNER_MEAN");
    697   liststats_setmode (&stats, "MEAN");
    698 
    699   Ndel = Nave = 0;
    700   Ncal = Nmos = Ngrp = Ngrid = Nfew = 0;
    701   for (i = 0; i < Ncatalog; i++) {
    702     for (j = 0; j < catalog[i].Naverage; j++) {
    703 
    704       int Ns;
    705       for (Ns = 0; Ns < Nphotcodes; Ns++) {
    706        
    707         int thisCode = photcodes[Ns][0].code;
    708         int Nsec = GetPhotcodeNsec(thisCode);
    709        
    710         /* on final processing, skip stars already measured */
    711         if (final && catalog[i].found_t[Nsecfilt*j + Nsec]) continue; 
    712 
    713         /* skip bad stars to prevent them from becoming good (on inner sample) */
    714         if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; 
    715 
    716         /* accumulate list of valid measurements */
    717         m = catalog[i].averageT[j].measureOffset;
    718         N = 0;
    719         for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
    720           // skip measurements that do not match the current photcode
    721           int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
    722           if (ecode != thisCode) { continue; }
    723 
    724           // SKIP gpc1 stack data
    725           if (isGPC1stack(catalog[i].measureT[m].photcode)) continue;
    726          
    727           // SKIP gpc1 forced-warp data
    728           if (isGPC1warp(catalog[i].measureT[m].photcode)) continue;
    729          
    730           // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers
    731           Mcal  = getMcal  (m, i, MAG_CLASS_PSF);
    732           if (isnan(Mcal)) { Ncal ++; continue; }
    733           Mmos  = getMmos  (m, i);
    734           if (isnan(Mmos)) { Nmos ++; continue; }
    735           Mgrp  = getMgrp  (m, i, catalog[i].measureT[m].airmass, NULL);
    736           if (isnan(Mgrp)) { Ngrp ++; continue; }
    737           Mgrid = getMgridTiny (&catalog[i].measureT[m]);
    738           if (isnan(Mgrid)) { Ngrid ++; continue; }
    739 
    740           Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt], MAG_CLASS_PSF);
    741           list[N] = Msys - Mcal - Mmos - Mgrp - Mgrid;
    742           dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR);
    743           N++;
    744         }
    745         if (N <= TOOFEW) {
    746           Nfew ++;
    747           continue;
    748         }
    749 
    750         /* 3-sigma clip based on stats of inner 50% */
    751 
    752         // calculated mean of inner 50%
    753         liststats (list, dlist, NULL, N, &instats);
    754         instats.sigma = MAX (MIN_ERROR, instats.sigma); /* if measurements agree too well, sigma -> 0.0 */
    755 
    756         // ignore entries > 3sigma from inner mean
    757         for (k = m = 0; k < N; k++) {
    758           if (fabs (list[k] - instats.median) < NSIGMA_CLIP*instats.sigma) {
    759             list[m] = list[k];
    760             m++;
    761           }
    762         }
    763         // recalculate the mean & sigma of the accepted measurements
    764         liststats (list, dlist, NULL, m, &stats);
    765         stats.sigma = MAX (MIN_ERROR, stats.sigma);
    766 
    767         /* apply to list of all relevant measurements, including IMAGE_POOR & IMAGE_FEW */
    768         image_bad = IMAGE_BAD;
    769         IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL;
    770         m = catalog[i].averageT[j].measureOffset;
    771         N = 0;
    772         for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
    773           // skip measurements that do not match the current photcode
    774           int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
    775           if (ecode != thisCode) { continue; }
    776 
    777           // SKIP gpc1 stack data
    778           if (isGPC1stack(catalog[i].measureT[m].photcode)) continue;
    779          
    780           // SKIP gpc1 forced-warp data
    781           if (isGPC1warp(catalog[i].measureT[m].photcode)) continue;
    782          
    783           // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers
    784           Mcal  = getMcal  (m, i, MAG_CLASS_PSF);
    785           if (isnan(Mcal)) continue;
    786           Mmos  = getMmos  (m, i);
    787           if (isnan(Mmos)) continue;
    788           Mgrid = getMgridTiny (&catalog[i].measureT[m]);
    789           if (isnan(Mgrid)) continue;
    790 
    791           Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt], MAG_CLASS_PSF);
    792           list[N] = Msys - Mcal - Mmos - Mgrid;
    793           dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR);
    794           ilist[N] = m;
    795           N++;
    796           Nave ++;
    797         }
    798         if (N < TOOFEW) continue;
    799 
    800         /* mark bad measures (> 5 sigma deviant) */
    801         for (k = 0; k < N; k++) {
    802           // treat the scatter of the star as a systematic term.  this is a bit of an
    803           // over-estimage (a perfect Gauss distribution with perfect errors would have
    804           // mySigma = sqrt(2) too large)
    805           float mySigma = hypot (stats.sigma, dlist[k]);
    806           if (fabs (list[k] - stats.median) > NSIGMA_REJECT*mySigma) {
    807             catalog[i].measureT[ilist[k]].dbFlags |= ID_MEAS_POOR_PHOTOM;
    808             if (final) {
    809               // for the final pass, we have a duplicate set of values in measure and measureT
    810               catalog[i].measure[ilist[k]].dbFlags |= ID_MEAS_POOR_PHOTOM;
    811             }
    812             Ndel ++;
    813           } else {
    814             catalog[i].measureT[ilist[k]].dbFlags &= ~ID_MEAS_POOR_PHOTOM;
    815             if (final) {
    816               // for the final pass, we have a duplicate set of values in measure and measureT
    817               catalog[i].measure[ilist[k]].dbFlags &= ~ID_MEAS_POOR_PHOTOM;
    818             }
    819           }
    820         }
    821         IMAGE_BAD = image_bad;
    822       }
    823     }
    824   }
    825   if (VERBOSE) fprintf (stderr, OFF_T_FMT" measures marked poor, "OFF_T_FMT" total\n", Ndel, Nave);
    826 
    827   free (list);
    828   free (dlist);
    829   free (ilist);
    830   free (tlist);
    831652}
    832653

  • branches/eam_branches/relphot.20210521/src/TGroupOps.c

    r41608 r41625  
    11721172  ALLOCATE (dlist, double, NtgroupTimes*Nphotcodes);
    11731173
    1174   int NIGHT_BAD = ID_IMAGE_NIGHT_POOR | ID_IMAGE_PHOTOM_FEW;
    1175 
    11761174  int n = 0;
    11771175  for (int i = 0; i < NtgroupTimes; i++) {
    11781176    TGroup *tgroup = tgroupTimes[i][0].byCode;
    11791177    for (int j = 0; j < tgroupTimes[i][0].nCode; j++) {
    1180       if (tgroup[j].flags & NIGHT_BAD) continue;
     1178      if (tgroup[j].flags & (ID_IMAGE_NIGHT_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
    11811179      list[n] = tgroup[j].McalPSF;
    11821180      dlist[n] = 1;
    … …  
    12091207    TGroup *tgroup = tgroupTimes[i][0].byCode;
    12101208    for (int j = 0; j < tgroupTimes[i][0].nCode; j++) {
    1211       if (tgroup[j].flags & IMAGE_BAD) continue;
     1209      if (tgroup[j].flags & (ID_IMAGE_NIGHT_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
    12121210      list[n] = tgroup[j].dMcal;
    12131211      dlist[n] = 1;
    … …  
    12401238    TGroup *tgroup = tgroupTimes[i][0].byCode;
    12411239    for (int j = 0; j < tgroupTimes[i][0].nCode; j++) {
    1242       if (tgroup[j].flags & IMAGE_BAD) continue;
     1240      if (tgroup[j].flags & (ID_IMAGE_NIGHT_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
    12431241      list[n] = tgroup[j].McalChiSq;
    12441242      dlist[n] = 1;
    … …  
    13961394        c = tgroup[i].catalog[j];
    13971395     
    1398         if (catalog[c].measureT[m].dbFlags & (ID_MEAS_AREA | ID_MEAS_NOCAL)) continue;
     1396        if (catalog[c].measureT[m].dbFlags & MEAS_BAD) continue;
    13991397
    14001398        // ave = catalog[c].measureT[m].averef;

  • branches/eam_branches/relphot.20210521/src/ZeroPointModes.c

    r41603 r41625  
    2222}
    2323
    24 void SetZeroPointModes (void) {
     24void SetZeroPointModes (Catalog *catalog, int Ncatalog) {
    2525
    2626  if (TGROUP_ZEROPT && MOSAIC_ZEROPT) {
    … …  
    5959      clean_mosaics(); // do this on every pass
    6060      clean_images(); // do this on every pass
     61      clean_stars(catalog, Ncatalog); // do this on every pass
    6162    }
    6263    if ((CurrentLoop >= 21) && (CurrentLoop <= 999)) {
    … …  
    7172      clean_mosaics(); // do this on every pass
    7273      clean_images(); // do this on every pass
     74      clean_stars(catalog, Ncatalog); // do this on every pass
    7375    }
    7476  }

  • branches/eam_branches/relphot.20210521/src/initialize.c

    r41603 r41625  
    2727  }   
    2828   
    29   IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP;
    30   STAR_BAD  = ID_OBJ_POOR | ID_OBJ_FEW;
    31   MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA;
    32 
    3329  if (SHOW_PARAMS) {
    3430    int Ns;

  • branches/eam_branches/relphot.20210521/src/plot_scatter.c

    r41606 r41625  
    2828    N = 0;
    2929    for (i = 0; i < Ncatalog; i++) {
    30         for (j = 0; j < catalog[i].Naverage; j++) {
     30      for (j = 0; j < catalog[i].Naverage; j++) {
    3131
    32             /* calculate the average value for a single star */
    33             if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; 
    34             m = catalog[i].averageT[j].measureOffset;
     32        /* calculate the average value for a single star */
     33        // if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; 
     34        m = catalog[i].averageT[j].measureOffset;
    3535
    36             for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
    37                 // skip measurements that do not match the current photcode
    38                 int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
    39                 if (ecode != thisCode) { continue; }
     36        for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
     37          // skip measurements that do not match the current photcode
     38          int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
     39          if (ecode != thisCode) { continue; }
    4040
    41                 if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue;
    42                 Mcal = getMcal  (m, i, MAG_CLASS_PSF);
    43                 if (isnan(Mcal)) continue;
    44                 Mmos = getMmos  (m, i);
    45                 if (isnan(Mmos)) continue;
    46                 Mgrid = getMgridTiny (&catalog[i].measureT[m]);
    47                 if (isnan(Mgrid)) continue;
     41          if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue;
     42          Mcal = getMcal  (m, i, MAG_CLASS_PSF);
     43          if (isnan(Mcal)) continue;
     44          Mmos = getMmos  (m, i);
     45          if (isnan(Mmos)) continue;
     46          Mgrid = getMgridTiny (&catalog[i].measureT[m]);
     47          if (isnan(Mgrid)) continue;
    4848
    49                 Mrel = catalog[i].secfilt[Nsecfilt*j+Nsec].MpsfChp;
    50                 if (isnan(Mrel)) continue;
     49          Mrel = catalog[i].secfilt[Nsecfilt*j+Nsec].MpsfChp;
     50          if (isnan(Mrel)) continue;
    5151
    52                 xlist[N] = Mrel;
    53                 ylist[N] = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt], MAG_CLASS_PSF) - Mcal - Mmos - Mgrid - Mrel;
    54                 ilist[N] = PhotInstTiny (&catalog[i].measureT[m], MAG_CLASS_PSF);
    55                 N++;
    56             }
     52          xlist[N] = Mrel;
     53          ylist[N] = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt], MAG_CLASS_PSF) - Mcal - Mmos - Mgrid - Mrel;
     54          ilist[N] = PhotInstTiny (&catalog[i].measureT[m], MAG_CLASS_PSF);
     55          N++;
    5756        }
     57      }
    5858    }
    5959

  • branches/eam_branches/relphot.20210521/src/relphot_images.c

    r41624 r41625  
    8383    SAVEPLOT = FALSE;
    8484
     85    // mark ID_MEAS_AREA (user-selected chip region) and ID_MEAS_NOCAL (time range & active photcode)
    8586    setExclusions (catalog, Ncatalog, TRUE);
    8687
    … …  
    114115
    115116      global_stats (catalog, Ncatalog, i);
    116       SetZeroPointModes ();
     117      SetZeroPointModes (catalog, Ncatalog);
    117118
    118119      MARKTIME("-- finished loop %d: %f sec\n", i, dtime);

  • branches/eam_branches/relphot.20210521/src/relphot_parallel_images.c

    r41606 r41625  
    88 * update the unowned detections for owned objects from neighbor regions
    99 * update the unowned objects for owned detections
     10
     11 * TGROUP and GRID corrections are not calculated.
    1012 */
    1113
    … …  
    105107    setMcal  (catalog);
    106108    setMmos  (catalog);
    107     // setMgrp  (catalog, FALSE); XXX think this through: this may not make sense for tgroups
    108109    MARKTIME("-- set Mrel, Mcal, Mmos, Mgrid : %f sec\n", dtime);
    109110   
    … …  
    117118    global_stats (catalog, Ncatalog, i);
    118119
    119     SetZeroPointModes ();
     120    SetZeroPointModes (catalog, Ncatalog);
    120121    MARKTIME("-- finished loop %d: %f sec\n", i, dtime);
    121122  }
    … …  
    138139
    139140  setMcalFromMosaics (); // copy per-mosaic calibrations to the images
    140   // XXX keep? setMcalFromTGroups (); // copy per-tgroup calibrations to the images
    141141
    142142  share_image_mags (regionHosts, -1);
    … …  
    152152  freeImageBins(Ncatalog, TRUE);
    153153  freeMosaicBins (Ncatalog, TRUE);
    154   freeGridBins ();
    155154  freeImages((char *)image);
    156155  free (image);

  • branches/eam_branches/relphot.20210521/src/setMrelCatalog.c

    r41624 r41625  
    225225    }
    226226
     227    /* some things to note here:
     228
     229       Mcal, Mmos, Mgrid are only relevant (non-zero) during the relphot_images and
     230       relphot_parallel_images analysis steps. For the final output steps (setMrelFinal),
     231       these are turned off
     232
     233       During the image analysis, Mflat is the prior stored value while Mgrid is the new
     234       measurement.  In the final output step, Mgrid is transferred to Mflat and zero
     235
     236       On the final calculation of Mrel, the value of Mgrid has been applied to Mflat, but
     237       it has not been removed from the GridOps.c structures.  In order to avoid
     238       double-counting, we need to skip Mgrid on the final calculation.
     239    */
     240
    227241    // ** Choose the calibration (depends on the mode : do I have an image reference or not?)
    228242    float Mcal = 0.0, Mmos = 0.0, Mgrid = 0.0, Mflat = 0.0, Mgrp = 0.0;
    … …  
    235249      Mcal = measureT[k].McalPSF; // check that this is zero for loaded REF value
    236250    } else {
    237       // getMcal returns image[].Mcal; note the flat-field correction is stored in measure.Mflat
     251      // getMcal returns image[].Mcal
    238252      Mcal  = getMcal (meas, cat, MAG_CLASS_PSF);
    239253      if (isnan(Mcal))  SKIP_THIS_MEAS(Ncal);
    240       Mmos  = getMmos  (meas, cat);
    241       if (isnan(Mmos))  SKIP_THIS_MEAS(Nmos);
    242       Mgrp  = getMgrp  (meas, cat, measureT[k].airmass, &dMgrp);
    243       if (isnan(Mgrp))  SKIP_THIS_MEAS(Ngrp);
    244       Mgrid = isSetMrelFinal ? 0.0 : getMgridTiny (&measureT[k]);
    245       if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid);
     254
     255      // the flat-field correction is stored in measure.Mflat
    246256      Mflat = isnan (measureT[k].Mflat) ? 0.0 : measureT[k].Mflat;
    247     }
    248     // on the final calculation of Mrel the value of Mgrid has been applied to Mflat, but
    249     // it has not been removed from the GridOps.c structures.  In order to avoid double-counting,
    250     // we need to skip Mgrid on the final calculation.
     257
     258      // see note above re: final output vs image analysis
     259      if (!isSetMrelFinal) {
     260        Mmos  = getMmos  (meas, cat);
     261        if (isnan(Mmos))  SKIP_THIS_MEAS(Nmos);
     262        Mgrp  = getMgrp  (meas, cat, measureT[k].airmass, &dMgrp);
     263        if (isnan(Mgrp))  SKIP_THIS_MEAS(Ngrp);
     264        Mgrid = getMgridTiny (&measureT[k]);
     265        if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid);
     266      }
     267    }
    251268
    252269    int myUbercalDist = getUbercalDist(meas, cat);
    … …  
    367384    }
    368385
     386    // if -preserve-ps1 is set and this object has PS1 data, skip the rest of the steps:
     387    // NOTE: This test only applies if this bit was set by an earlier analysis
    369388    if (PRESERVE_PS1 && (secfilt[Nsec].flags & ID_SECF_HAS_PS1)) continue;
    370     // if -preserve-ps1 is set and this object has PS1 data, skip the rest of the steps:
    371389
    372390    // XXX hardwired grizy = (01234) JHK = (567) w = (8)
    … …  
    419437    int NrankingPSF = useOLS ? magStatsByRanking (&results->psfData[Nsec], psfstats) : magStatsByRankingIRLS (&results->psfData[Nsec], psfstats);
    420438    if (NrankingPSF < Nminmeas) {
    421       secfilt[Nsec].flags |= ID_OBJ_FEW;
     439      secfilt[Nsec].flags |= ID_SECF_STAR_FEW;
    422440    } else {
    423441      secfilt[Nsec].MpsfChp  = psfstats->mean;
    424442      secfilt[Nsec].dMpsfChp = psfstats->error;
    425443      secfilt[Nsec].Mchisq   = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN;
     444      secfilt[Nsec].flags &= ~ID_SECF_STAR_FEW;
    426445    }
    427446
    … …  
    629648      haveStackObject = TRUE;
    630649
    631       // if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS_STACK(Nbad);
    632 
    633650      int isPrimary = FALSE;
    634651     

  • branches/eam_branches/relphot.20210521/src/setMrelFinal.c

    r41624 r41625  
    1313  ResetAverageAndMeasure (catalog);
    1414 
    15   // this sets flags in the measureT element, not the measure element
     15  // This sets ID_MEAS flags in the measureT element, not the measure element
     16  // These bits are thus ephemeral and not saved. 
    1617  setExclusions (catalog, 1, VERBOSE);  /* mark by area */
    17   // XXXXX THis probably does not make sense
    1818
    1919  /* set catalog[0].found[i] = FALSE */
    … …  
    3131  setMrelOutput (catalog, 1); // sets the values secfilt.MpsfChp = <measure.M - image.Mcal - measure.Mflat>
    3232  setMcalOutput (catalog, 1); // sets measure.Mcal = image.Mcal
    33 
    34   // clear ID_OBJ_POOR, ID_OBJ_FEW, ID_MEAS_NOCAL values before writing
    35   // These are internal bits, so they should be cleared
    36   // XXX are these still used?
    37   // XXX eventually: move these internal flags to an internal-only catalog->FOO element (like found)
    38   for (i = 0; i < catalog[0].Naverage; i++) {
    39     catalog[0].average[i].flags &= ~ID_OBJ_FEW;
    40     catalog[0].average[i].flags &= ~ID_OBJ_POOR;
    41     off_t j;
    42     for (j = 0; j < Nsecfilt; j++) {
    43         catalog[0].secfilt[i*Nsecfilt+j].flags &= ~ID_OBJ_FEW;
    44         catalog[0].secfilt[i*Nsecfilt+j].flags &= ~ID_OBJ_POOR;
    45     }
    46     off_t m = catalog[0].average[i].measureOffset;
    47     for (j = 0; j < catalog[0].average[i].Nmeasure; j++, m++) {
    48       catalog[0].measure[m].dbFlags &= ~ID_MEAS_NOCAL;
    49     }
    50   }
    5133}
    52 
    53 /* ID_MEAS_SKIP marks measurements which were not used to calculate Mrel */
    54 // For each star & photcode, we are selecting the measurements to use in the average
    55 // photometry value.  This function is called several times, lowering the bar to
    56 // acceptance on each pass. The rejections at the different passes are:
    57 
    58 // 0 : only use measurements thought to be GOOD (photflags not POOR or BAD)
    59 // 1 : accept measurements thought to be POOR based on photflags
    60 // 2 : accept the measurements marked as outliers; accept images marked as outliers
    61 // 3 : accept measurements thought to be BAD based on photflags
    62 // 4 : accept the measurements outside of the instrumental magnitude limits (eg, SAT)
    63 
    64 // XXX if I'm reading this correctly, passes 0, 1, 2 are identical & accept all
    65 // non-outliers with Imag range
    6634
    6735// int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure) {
    … …  
    8452// }
    8553
    86 int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure);
    87 
    88 // XXX this function is no longer used:
    89 void skip_measurements (Catalog *catalog, int pass) {
    90 
    91   off_t i, k, m;
    92   off_t Ntot, Ntry, Nkeep, Nskip;
    93   float mag;
    94 
    95   int Nsecfilt = GetPhotcodeNsecfilt ();
    96 
    97   Ntot = Ntry = Nskip = Nkeep = 0;
    98 
    99   /* allow measures from images marked POOR and FEW */
    100   if (pass >= 3) IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL;
    101  
    102   /* allow measures marked as outliers (POOR) and off image region (AREA) */
    103   if (pass >= 3) {
    104     MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_SKIP_PHOTOM;
    105   } else {
    106     MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA;
    107   }
    108 
    109   /* mark measures which should be ignored on this pass */
    110   for (i = 0; i < catalog[0].Naverage; i++) {
    111     Ntot += catalog[0].average[i].Nmeasure;
    112 
    113     if (FALSE && (catalog[0].average[i].objID == 0x46a4) && (catalog[0].average[i].catID == 0xf40e)) {
    114       fprintf (stderr, "test obj\n");
    115       print_measure_set (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], catalog[0].measure);
    116     }
    117 
    118     // mark measurements for each secfilt separately
    119     int Ns;
    120     for (Ns = 0; Ns < Nphotcodes; Ns++) {
    121      
    122       int thisCode = photcodes[Ns][0].code;
    123       int Nsec = GetPhotcodeNsec(thisCode);
    124      
    125       /* star/photcodes already calibrated */
    126       if (catalog[0].found_t[Nsecfilt*i+Nsec]) continue; 
    127 
    128       m = catalog[0].average[i].measureOffset;
    129       for (k = 0; k < catalog[0].average[i].Nmeasure; k++, m++) {
    130         Ntry++;
    131 
    132         // skip measurements not related to this photcode
    133         PhotCode *code = GetPhotcodebyCode (catalog[0].measure[m].photcode);
    134         if (!code) continue;
    135         if (code->equiv != thisCode) continue;
    136 
    137         // skip measurements by time range (mark as skipped for this secfilt)
    138         // XXX note that this is a bit dangerous : some objects may never get calibrated
    139         if (TimeSelect) {
    140           if (catalog[0].measure[m].t < TSTART) goto skip;
    141           if (catalog[0].measure[m].t > TSTOP) goto skip;
    142         }
    143        
    144         // we now have a measurement relavant to this photcode (and time range)
    145 
    146         // clear SKIP for all measures at first
    147         catalog[0].measureT[m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
    148         catalog[0].measure [m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
    149 
    150         // skip measurements from BAD images and mosaics (not REF mags)
    151         // do NOT skip measurements without a matching image (REF mags)
    152         off_t Nim = getImageEntry (m, 0);
    153         if (Nim > -1) {
    154           if (isnan(getMcal (m, 0, MAG_CLASS_PSF))) goto skip;
    155           if (isnan(getMmos (m, 0))) goto skip;
    156         }
    157 
    158         // PASS 4 : skip measurements by inst mag limit (not REF mags)
    159         if ((pass < 4) && ImagSelect) {
    160           if (Nim > -1) {
    161             mag = PhotInst (&catalog[0].measure[m], MAG_CLASS_PSF);
    162             if (mag < ImagMin) goto skip;
    163             if (mag > ImagMax) goto skip;
    164           }
    165         }
    166 
    167         // PASS 3 : accept bad measurements (eg, SAT, CR), internal outliers
    168         if (pass < 3) {
    169           if (catalog[0].measure[m].photFlags & code->photomBadMask) goto skip;
    170           if ((catalog[0].measure[m].photcode > 10000) && (catalog[0].measure[m].photcode < 10500)) {
    171             if (catalog[0].measure[m].psfQF < 0.85) goto skip;
    172           }
    173         }
    174        
    175         // PASS 2 : (no additional cut)
    176 
    177         // PASS 1 : accept poor measurements as well (eg, POOR FIT, etc)
    178         if ((pass < 1) && (catalog[0].measure[m].photFlags & code->photomPoorMask)) goto skip;
    179 
    180         // PASS 0 : skip poor, outlier, bad measurements as well (eg, POOR FIT, etc)
    181 
    182         Nkeep ++;
    183        
    184         continue;
    185 
    186       skip:
    187         catalog[0].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM;
    188         catalog[0].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM;
    189         Nskip ++;
    190       }
    191     }
    192   }
    193   if (VERBOSE) fprintf (stderr, "pass %d, Ntot: "OFF_T_FMT", Ntry: "OFF_T_FMT", Nskip: "OFF_T_FMT", Nkeep: "OFF_T_FMT"\n",
    194                         pass, Ntot, Ntry, Nskip, Nkeep);
    195 }
    196 
    197 // XXX review this, especially the good/bad images
    198 // by this point, we have propagated the mosaic and tgroup flags to each image
     54// This function is only called for the final output step.  By this point, we have
     55// propagated the mosaic and tgroup flags to each image.
    19956void setMeasureRank (Catalog *catalog) {
    20057
    20158  int i;
    20259
    203   // make these global
    204   int IMAGE_BAD  = ID_IMAGE_PHOTOM_NOCAL;
    205   int IMAGE_POOR = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP;
    206   int MEAS_BAD   = ID_MEAS_NOCAL        | ID_MEAS_SKIP_PHOTOM;
    207   int MEAS_POOR  = ID_MEAS_POOR_PHOTOM  | ID_MEAS_AREA;
    208 
    20960  Measure     *measure     = catalog[0].measure;
    210   MeasureTiny *measureT    = catalog[0].measureT;
    21161  char        *measureRank = catalog[0].measureRank;
    21262
    … …  
    21565    measureRank[i] = 11; // start at a low rank
    21666
    217     // clear SKIP for all measures at first
    218     measureT[i].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
    219     measure [i].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
    220 
    221     // skip measurements without a valid photcode
     67    // measurements without a valid photcode have lowest rank (should not be used anyway)
    22268    PhotCode *code = GetPhotcodebyCode (measure[i].photcode);
    22369    if (!code) continue;
    … …  
    22571    // measurements outside time range have poor rank
    22672    if (TimeSelect) {
    227       if (measure[i].t < TSTART)                     { measureRank[i] = 10; continue; }
    228       if (measure[i].t > TSTOP)                      { measureRank[i] = 10; continue; }
     73      if (measure[i].t < TSTART)                                        { measureRank[i] = 10; continue; }
     74      if (measure[i].t > TSTOP)                                         { measureRank[i] = 10; continue; }
    22975    }
    23076   
    231     // measurements from BAD images and mosaics (not REF mags)
    232     off_t Nim = getImageEntry (i, 0);
     77    // Nim < 0 for REF mags, imageFlags have bits for IMAGE, MOSAIC, NIGHT
     78    off_t Nim      = getImageEntry (i, 0);
    23379    int imageFlags = getImageFlags (i, 0);
    234     int mosaicFlags = getMosaicFlags (i, 0);
    23580
    23681    if (Nim > -1) {
    237       if (imageFlags  & IMAGE_BAD)                   { measureRank[i] = 9; continue; }
    238 
    239       if (mosaicFlags & IMAGE_BAD)                   { measureRank[i] = 9; continue; }
    240    
    24182      // measurements ranked by inst mag limit (not REF mags)
    24283      if (ImagSelect) {
    24384        float mag = PhotInst (&measure[i], MAG_CLASS_PSF);
    244         if (mag < ImagMin)                           { measureRank[i] = 8; continue; }
    245         if (mag > ImagMax)                           { measureRank[i] = 8; continue; }
     85        if (mag < ImagMin)                                              { measureRank[i] = 9; continue; }
     86        if (mag > ImagMax)                                              { measureRank[i] = 9; continue; }
    24687      }
    24788    }
    24889
    249     // RANK 1 : BAD photFlags
    250     if (measure[i].dbFlags & MEAS_BAD)               { measureRank[i] = 7; continue; }
    251 
    252     // RANK 3 : BAD photFlags (eg, SAT, CR), internal outliers
    253     if (measure[i].photFlags & code->photomBadMask)  { measureRank[i] = 6; continue; }
    254 
     90    // RANK 8 : Poor image
    25591    if (Nim > -1) {
    256       if (imageFlags  & IMAGE_POOR)                  { measureRank[i] = 5; continue; }
    257       if (mosaicFlags & IMAGE_POOR)                  { measureRank[i] = 5; continue; }
     92      if (imageFlags & ID_IMAGE_PHOTOM_POOR)                            { measureRank[i] = 8; continue; }
    25893    }
    25994
    260     // RANK 2 : psfQF value
    261     if (!isfinite(measure[i].psfQF) || measure[i].psfQF < 0.85) { measureRank[i] = 4; continue; }
     95    // RANK 7 : BAD photFlags (eg, SAT, CR), internal outliers
     96    if (measure[i].photFlags & code->photomBadMask)                     { measureRank[i] = 7; continue; }
     97
     98    // RANK 6 : bad psfQF value
     99    if (!isfinite(measure[i].psfQF) || measure[i].psfQF < 0.85)         { measureRank[i] = 6; continue; }
    262100       
    263     // RANK 1 : POOR photFlags
    264     if (measure[i].dbFlags & MEAS_POOR)              { measureRank[i] = 3; continue; }
     101    // RANK 5 : not in valid chip region
     102    if (measure[i].dbFlags & ID_MEAS_AREA)                              { measureRank[i] = 5; continue; }
    265103
    266     // RANK 1 : POOR photFlags
    267     if (measure[i].photFlags & code->photomPoorMask) { measureRank[i] = 2; continue; }
     104    // RANK 4 : POOR photFlags
     105    if (measure[i].photFlags & code->photomPoorMask)                    { measureRank[i] = 4; continue; }
    268106   
    269     // RANK 1 : psfQFperfect value
    270     if (!isfinite(measure[i].psfQFperf) || measure[i].psfQFperf < 0.85) { measureRank[i] = 1; continue; }
     107    // RANK 3 : bad psfQFperfect value
     108    if (!isfinite(measure[i].psfQFperf) || measure[i].psfQFperf < 0.85) { measureRank[i] = 3; continue; }
    271109       
     110    // RANK 2 : Poor mosaic
     111    if (Nim > -1) {
     112      if (imageFlags & ID_IMAGE_MOSAIC_POOR)                            { measureRank[i] = 2; continue; }
     113    }
     114
     115    // RANK 1 : Poor night
     116    if (Nim > -1) {
     117      if (imageFlags & ID_IMAGE_NIGHT_POOR)                             { measureRank[i] = 1; continue; }
     118    }
    272119    // RANK 0 : perfect measurement:
    273120    measureRank[i] = 0;
Note: See TracChangeset for help on using the changeset viewer.