IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Apr 1, 2012, 2:48:41 PM (14 years ago)
Author:
eugene
Message:

updates to relphot to support ubercal, distributed dvo; some optmimization (eg, image matching to area); cleanup names of static helper arrays in ImageOps, MosaicOps; cleanup rules for -reset options (relphot does not reset astrom bits); add flat-field correction; determine ubercal distance for images; better plot ranges; bcatalog limit density selecting objects by number of measurements; allow relphot to update the output catalog schema (catformat); make relphot -averages and relphot -update-objects equivalent in how low-quality objects are handled; record stage at which average is determined

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r31668 r33651  
    11# include "relphot.h"
    22
    3 extern double drand48();
    4 
    5 int CopyAverageTiny (AverageTiny *averageT, Average *average) {
    6 
    7   averageT[0].R             = average[0].R;
    8   averageT[0].D             = average[0].D;
    9   averageT[0].flags         = average[0].flags;
    10   averageT[0].Nmeasure      = average[0].Nmeasure;
    11   averageT[0].measureOffset = average[0].measureOffset;
    12 
    13   // make Nmeasure & measureOffset optional?
    14 
    15   return (TRUE);
    16 }
    17 
    18 int CopyMeasureTiny (MeasureTiny *measureT, Measure *measure) {
    19 
    20   measureT[0].dR       = measure[0].dR;
    21   measureT[0].dD       = measure[0].dD;
    22   measureT[0].M        = measure[0].M;
    23   measureT[0].Mcal     = measure[0].Mcal;
    24   measureT[0].dM       = measure[0].dM;
    25   measureT[0].airmass  = measure[0].airmass;
    26   measureT[0].Xccd     = measure[0].Xccd;
    27   measureT[0].Yccd     = measure[0].Yccd;
    28   measureT[0].t        = measure[0].t;
    29   measureT[0].dt       = measure[0].dt;
    30   measureT[0].averef   = measure[0].averef;
    31   measureT[0].imageID  = measure[0].imageID;
    32   measureT[0].dbFlags  = measure[0].dbFlags;
    33   measureT[0].photcode = measure[0].photcode;
    34 
    35   return (TRUE);
    36 }
     3int LimitDensityCatalog_ByNmeasure (Catalog *subcatalog, Catalog *oldcatalog);
    374
    385int bcatalog (Catalog *subcatalog, Catalog *catalog) {
    396 
    407  off_t i, j, offset;
    41   int ecode, found, Ns;
     8  int found, Ns, *Nvalid;
    429  off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm;
    4310  float mag;
    44   int Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Ngalaxy, Npsfqf;
     11  int Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Ngalaxy, Npsfqf, Nnan, Nbad, Npoor;
    4512
    4613  int Nsecfilt = GetPhotcodeNsecfilt ();
    … …  
    5320  ALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE);
    5421  ALLOCATE (subcatalog[0].averageT, AverageTiny, NAVERAGE);
    55   ALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt);
     22  ALLOCATE (subcatalog[0].secfilt,  SecFilt,    NAVERAGE*Nsecfilt);
    5623  Nmeasure = Naverage = 0;
    5724
    58   Ncode = Ntime = Ndophot = Nmag = Nsigma = Nimag = Nfew = Npsfqf = Ngalaxy = 0;
     25  Ncode = Ntime = Ndophot = Nmag = Nsigma = Nimag = Nfew = Npsfqf = Ngalaxy = Nnan = Nbad = Npoor = 0;
     26
     27  // copy the following fields to the subcatalog:
     28  subcatalog[0].catID = catalog[0].catID;
     29  subcatalog[0].filename = catalog[0].filename;
     30
     31  // array to count measures with photcodes->equiv to each of the secfilts
     32  ALLOCATE (Nvalid, int, Nsecfilt);
     33
     34  // flags used by the photometry analysis (excluding UBERCAL)
     35  unsigned int PHOTOM_FLAGS =
     36    ID_MEAS_NOCAL |       // detection ignored for this analysis (photcode, time range)
     37    ID_MEAS_POOR_PHOTOM | // detection is photometry outlier
     38    ID_MEAS_SKIP_PHOTOM | // detection was ignored for photometry measurement
     39    ID_MEAS_AREA;         // detetion was outside acceptable area of device
    5940
    6041  /* exclude stars not in range or with too few measurements */
    … …  
    6344
    6445    /* start with all stars good */
    65     CopyAverageTiny (&subcatalog[0].averageT[Naverage], &catalog[0].average[i]);
     46    CopyAverageToTiny (&subcatalog[0].averageT[Naverage], &catalog[0].average[i]);
    6647    subcatalog[0].averageT[Naverage].measureOffset = Nmeasure;
    6748
    6849    for (j = 0; j < Nsecfilt; j++) {
    6950      subcatalog[0].secfilt[Nsecfilt*Naverage+j] = catalog[0].secfilt[Nsecfilt*i+j];
    70     }
    71 
     51      Nvalid[j] = 0; // reset the Nvalid array for this star
     52    }
     53
     54    // reset the calculated average magnitudes (does not affect ubercal-tied measurements or images)
    7255    if (RESET) {
    7356      int Ns;
     57      DVOSecfiltFlags secfiltBits =
     58        ID_SECF_STAR_FEW |
     59        ID_SECF_STAR_POOR |
     60        ID_PHOTOM_PASS_0 |
     61        ID_PHOTOM_PASS_1 |
     62        ID_PHOTOM_PASS_2 |
     63        ID_PHOTOM_PASS_3 |
     64        ID_PHOTOM_PASS_4 |
     65        ID_SECF_USE_SYNTH |
     66        ID_SECF_USE_UBERCAL |
     67        ID_SECF_OBJ_EXT;
    7468      for (Ns = 0; Ns < Nphotcodes; Ns++) {
    7569
    … …  
    7771        int Nsec = GetPhotcodeNsec(thisCode);
    7872
    79         subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M  = NAN;
    80         subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].dM = NAN;
    81         subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].flags &= ~ID_STAR_FEW;
    82         subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].flags &= ~ID_STAR_POOR;
     73        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M      = NAN;
     74        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Map    = NAN;
     75        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].dM     = NAN;
     76        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Xm     = NAN;
     77        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M_20   = NAN_S_SHORT;
     78        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M_80   = NAN_S_SHORT;
     79        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Ncode  = 0;
     80        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Nused  = 0;
     81        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].flags &= ~secfiltBits;
     82        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].ubercalDist = 1000;
    8383      }
    8484    }
    … …  
    9292
    9393      /* select measurements by photcode */
    94       ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[offset].photcode);
     94      PhotCode *code = GetPhotcodebyCode (catalog[0].measure[offset].photcode);
     95      if (!code) continue;
    9596      found = FALSE;
     97      int Nsec = -1; // Nsec equivalent for the measurement
    9698      for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) {
    97         if (ecode == photcodes[Ns][0].code) found = TRUE;
     99        if (code->equiv != photcodes[Ns][0].code) continue;
     100        found = TRUE;
     101        Nsec = GetPhotcodeNsec(code->equiv);
    98102      }
    99103      if (!found) {
    … …  
    113117      // skip garbage measurements
    114118      if (catalog[0].measure[offset].psfQual < 0.85) { Npsfqf ++; continue; }
    115 
    116       // check for galaxies
    117       if (!isnan(catalog[0].measure[offset].Map)) {
     119      if (isnan(catalog[0].measure[offset].M)) { Nnan ++; continue; }
     120
     121      if (catalog[0].measure[offset].photFlags & code->photomBadMask) {
     122        Nbad++;
     123        continue;
     124      }
     125      if (catalog[0].measure[offset].photFlags & code->photomPoorMask) { Npoor++; continue;}
     126
     127      // check for galaxies (XXX skip for now)
     128      if (FALSE && !isnan(catalog[0].measure[offset].Map)) {
    118129        if (catalog[0].measure[offset].M - catalog[0].measure[offset].Map > 0.15) {
    119130          nEXT ++;
    … …  
    137148      }
    138149
    139       CopyMeasureTiny (&subcatalog[0].measureT[Nmeasure], &catalog[0].measure[offset]);
     150      // count this measurement as valid for this secfilt entry
     151      if (Nsec > -1) {
     152        assert (Nsec < Nsecfilt);
     153        Nvalid[Nsec] ++;
     154      }
     155
     156      CopyMeasureToTiny (&subcatalog[0].measureT[Nmeasure], &catalog[0].measure[offset]);
    140157      subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
    141158      subcatalog[0].measureT[Nmeasure].averef = Naverage;
    142159      if (RESET) {
    143         subcatalog[0].measureT[Nmeasure].Mcal = 0;
    144         subcatalog[0].measureT[Nmeasure].dbFlags &= 0xff00;
    145         subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_POOR_PHOTOM;
    146         subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_AREA;
    147         subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_NOCAL;
     160        if (!KEEP_UBERCAL) {
     161          subcatalog[0].measureT[Nmeasure].Mcal = 0;
     162          subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL;
     163        }
     164        subcatalog[0].measureT[Nmeasure].dbFlags &= ~PHOTOM_FLAGS;
    148165      }
    149166      Nmeasure ++;
    … …  
    156173
    157174    // skip object if it is likely to be a galaxy
    158     if (nEXT >= nPSF) {
     175    if (FALSE && (nEXT >= nPSF)) {
    159176      Nmeasure -= Nm;
    160177      Ngalaxy ++;
    … …  
    163180
    164181    // XXXX test : what checks do I need to make elsewhere to avoid problems here?
    165     if (Nm <= STAR_TOOFEW) { /* enough measurements in band? */
     182    if (Nm <= STAR_TOOFEW) { /* enough measurements total? */
    166183      Nmeasure -= Nm;
    167184      Nfew ++;
    168185      continue;
    169186    }
     187    int anySecfiltGood = FALSE;
     188    for (Ns = 0; Ns < Nsecfilt; Ns++) { /* enough measurements in at least one band? */
     189      if (Nvalid[Ns] <= STAR_TOOFEW) continue;
     190      anySecfiltGood = TRUE;
     191    }
     192    if (!anySecfiltGood) {
     193      Nmeasure -= Nm;
     194      Nfew ++;
     195      continue;
     196    }
     197
    170198    subcatalog[0].averageT[Naverage].Nmeasure = Nm;
    171199    Naverage ++;
    … …  
    186214  // limit the total number of stars in the catalog
    187215  if (MaxDensityUse) {
    188     LimitDensityCatalog (subcatalog, catalog);
     216    LimitDensityCatalog_ByNmeasure (subcatalog, catalog);
    189217  }
    190218
    … …  
    192220    fprintf (stderr, "using "OFF_T_FMT" stars ("OFF_T_FMT" measures) of "OFF_T_FMT" for catalog %s\n",
    193221             subcatalog[0].Naverage,  subcatalog[0].Nmeasure,  i, catalog[0].filename);
    194     fprintf (stderr, "rejections: %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d few, %d psfqf, %d galaxies\n",
    195              Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Npsfqf, Ngalaxy);
     222    fprintf (stderr, "rejections: %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d few, %d psfqf, %d Nnan, %d galaxies, %d bad, %d poor\n",
     223             Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Npsfqf, Nnan, Ngalaxy, Nbad, Npoor);
    196224  }
    197225  return (TRUE);
    198226}
    199227
    200 int LimitDensityCatalog (Catalog *subcatalog, Catalog *catalog) {
     228// sort by decreasing Nmeasure (X)
     229void sort_by_Nmeasure (int *X, off_t *Y, off_t N) {
     230
     231# define SWAPFUNC(A,B){ int tmpI; off_t tmpT;   \
     232  tmpI = X[A]; X[A] = X[B]; X[B] = tmpI; \
     233  tmpT = Y[A]; Y[A] = Y[B]; Y[B] = tmpT; \
     234}
     235# define COMPARE(A,B)(X[A] > X[B])
     236
     237  OHANA_SORT (N, COMPARE, SWAPFUNC);
     238
     239# undef SWAPFUNC
     240# undef COMPARE
     241
     242}
     243
     244int LimitDensityCatalog_ByNmeasure (Catalog *subcatalog, Catalog *oldcatalog) {
    201245
    202246  Catalog tmpcatalog;
    … …  
    206250  int Nsecfilt = GetPhotcodeNsecfilt ();
    207251
    208   gfits_scan (&catalog[0].header, "RA0",  "%lf", 1, &Rmin);
    209   gfits_scan (&catalog[0].header, "DEC0", "%lf", 1, &Dmin);
    210   gfits_scan (&catalog[0].header, "RA1",  "%lf", 1, &Rmax);
    211   gfits_scan (&catalog[0].header, "DEC1", "%lf", 1, &Dmax);
     252  gfits_scan (&oldcatalog[0].header, "RA0",  "%lf", 1, &Rmin);
     253  gfits_scan (&oldcatalog[0].header, "DEC0", "%lf", 1, &Dmin);
     254  gfits_scan (&oldcatalog[0].header, "RA1",  "%lf", 1, &Rmax);
     255  gfits_scan (&oldcatalog[0].header, "DEC1", "%lf", 1, &Dmax);
    212256
    213257  if (VERBOSE2) fprintf (stderr, "extracting from catalog covering region %f,%f to %f,%f\n", Rmin, Dmin, Rmax, Dmax);
    … …  
    229273
    230274  // we are going to select Nmax entries by generating a random-sorted index list
    231   off_t *index, tmp, i, j, ave;
     275  int *value;
     276  off_t *index, i, j, ave;
    232277  ALLOCATE (index, off_t, Naverage);
     278  ALLOCATE (value, int, Naverage);
    233279  for (i = 0; i < Naverage; i++) {
    234280    index[i] = i;
    235   }
    236   for (i = 0; i < Naverage; i++) {
    237     j = (Naverage - i) * drand48() + i; // a number between i and Naverage
    238     tmp = index[j];
    239     index[j] = index[i];
    240     index[i] = tmp;
    241   }
     281    value[i] = subcatalog[0].averageT[i].Nmeasure;
     282  }
     283  sort_by_Nmeasure (value, index, Naverage);
    242284
    243285  // count the number of measurements this selection will yield
    … …  
    266308      Nmeasure ++;
    267309    }
     310    for (j = 0; j < Nsecfilt; j++) {
     311      tmpcatalog.secfilt[i*Nsecfilt + j] = subcatalog[0].secfilt[ave*Nsecfilt + j];
     312    }
     313  }
     314
     315  if (VERBOSE) {
     316    char *basename = filebasename (oldcatalog[0].filename);
     317    fprintf (stderr, "limited to "OFF_T_FMT" ("OFF_T_FMT" subset, "OFF_T_FMT" total) stars, "OFF_T_FMT" ("OFF_T_FMT" subset, "OFF_T_FMT" total) measures for catalog %s\n",
     318             Nmax, subcatalog[0].Naverage, oldcatalog[0].Naverage, Nmeasure, subcatalog[0].Nmeasure,  oldcatalog[0].Nmeasure, basename);
     319    free (basename);
     320  }
     321
     322  free (index);
     323  free (value);
     324  free (subcatalog[0].averageT);
     325  free (subcatalog[0].measureT);
     326  free (subcatalog[0].secfilt);
     327
     328  subcatalog[0].averageT = tmpcatalog.averageT;
     329  subcatalog[0].measureT = tmpcatalog.measureT;
     330  subcatalog[0].secfilt = tmpcatalog.secfilt;
     331  subcatalog[0].Naverage = Nmax;
     332  subcatalog[0].Nmeasure = Nmeasure;
     333  subcatalog[0].Nsecfilt = oldcatalog[0].Nsecfilt;
     334  subcatalog[0].Nsecf_mem = Naverage * oldcatalog[0].Nsecfilt;
     335
     336  return (TRUE);
     337}
     338
     339int LimitDensityCatalog (Catalog *subcatalog, Catalog *catalog) {
     340
     341  Catalog tmpcatalog;
     342
     343  double Rmin, Rmax, Dmin, Dmax;
     344
     345  int Nsecfilt = GetPhotcodeNsecfilt ();
     346
     347  gfits_scan (&catalog[0].header, "RA0",  "%lf", 1, &Rmin);
     348  gfits_scan (&catalog[0].header, "DEC0", "%lf", 1, &Dmin);
     349  gfits_scan (&catalog[0].header, "RA1",  "%lf", 1, &Rmax);
     350  gfits_scan (&catalog[0].header, "DEC1", "%lf", 1, &Dmax);
     351
     352  if (VERBOSE2) fprintf (stderr, "extracting from catalog covering region %f,%f to %f,%f\n", Rmin, Dmin, Rmax, Dmax);
     353
     354  float AREA = fabs(Dmax - Dmin) * fabs(Rmax - Rmin) * cos (0.5*RAD_DEG*(Dmax + Dmin));
     355  assert (AREA > 0);
     356
     357  off_t Nmax = MaxDensityValue * AREA;
     358  if (subcatalog[0].Naverage <= Nmax) {
     359    if (VERBOSE) {
     360      fprintf (stderr, "subcatalog has less than the max density\n");
     361    }
     362    return (TRUE);
     363  }
     364
     365  off_t Naverage = subcatalog[0].Naverage;
     366
     367  // select a random subset of Nmax stars from subcatalog using Fisher-Yates
     368
     369  // we are going to select Nmax entries by generating a random-sorted index list
     370  off_t *index, tmp, i, j, ave;
     371  ALLOCATE (index, off_t, Naverage);
     372  for (i = 0; i < Naverage; i++) {
     373    index[i] = i;
     374  }
     375  for (i = 0; i < Naverage; i++) {
     376    j = (Naverage - i) * drand48() + i; // a number between i and Naverage
     377    tmp = index[j];
     378    index[j] = index[i];
     379    index[i] = tmp;
     380  }
     381
     382  // count the number of measurements this selection will yield
     383  off_t NMEASURE = 0;
     384  for (i = 0; i < Nmax; i++) {
     385    ave = index[i];
     386    NMEASURE += subcatalog[0].averageT[ave].Nmeasure;
     387  }
     388
     389  // allocate the output data
     390  ALLOCATE (tmpcatalog.averageT, AverageTiny, Nmax);
     391  ALLOCATE (tmpcatalog.measureT, MeasureTiny, NMEASURE);
     392  ALLOCATE (tmpcatalog.secfilt,  SecFilt, Nmax * Nsecfilt);
     393
     394  off_t Nmeasure = 0;
     395
     396  // copy the Nmax selected entries from subcatalog to tmpcatalog (adjusting links)
     397  for (i = 0; i < Nmax; i++) {
     398    ave = index[i];
     399    tmpcatalog.averageT[i] = subcatalog[0].averageT[ave];
     400    tmpcatalog.averageT[i].measureOffset = Nmeasure;
     401    for (j = 0; j < tmpcatalog.averageT[i].Nmeasure; j++) {
     402      off_t offset = subcatalog[0].averageT[ave].measureOffset + j;
     403      tmpcatalog.measureT[Nmeasure] = subcatalog[0].measureT[offset];
     404      tmpcatalog.measureT[Nmeasure].averef = i;
     405      Nmeasure ++;
     406    }
     407    for (j = 0; j < Nsecfilt; j++) {
     408      tmpcatalog.secfilt[i*Nsecfilt + j] = subcatalog[0].secfilt[ave*Nsecfilt + j];
     409    }
    268410  }
    269411
    … …  
    287429  return (TRUE);
    288430}
    289 
    290 // for the cases where we are not using a subset of the data, we still need to have a copy of these fields
    291 int populate_tiny_values (Catalog *catalog) {
    292 
    293   off_t i;
    294 
    295   AverageTiny *averageT;
    296   Average *average;
    297 
    298   MeasureTiny *measureT;
    299   Measure *measure;
    300 
    301   ALLOCATE (catalog[0].measureT, MeasureTiny, catalog[0].Nmeasure);
    302   ALLOCATE (catalog[0].averageT, AverageTiny, catalog[0].Naverage);
    303 
    304   average  = catalog[0].average;
    305   averageT = catalog[0].averageT;
    306 
    307   measure  = catalog[0].measure;
    308   measureT = catalog[0].measureT;
    309 
    310   for (i = 0; i < catalog[0].Naverage; i++) {
    311     // CopyAverageTiny (&catalog[0].averageT[i], &catalog[0].average[i]);
    312     averageT[i].R             = average[i].R;
    313     averageT[i].D             = average[i].D;
    314     averageT[i].flags         = average[i].flags;
    315     averageT[i].Nmeasure      = average[i].Nmeasure;
    316     averageT[i].measureOffset = average[i].measureOffset;
    317   }
    318 
    319   for (i = 0; i < catalog[0].Nmeasure; i++) {
    320     // CopyMeasureTiny (&catalog[0].measureT[i], &catalog[0].measure[i]);
    321     measureT[i].dR       = measure[i].dR;
    322     measureT[i].dD       = measure[i].dD;
    323     measureT[i].M        = measure[i].M;
    324     measureT[i].Mcal     = measure[i].Mcal;
    325     measureT[i].dM       = measure[i].dM;
    326     measureT[i].airmass  = measure[i].airmass;
    327     measureT[i].Xccd     = measure[i].Xccd;
    328     measureT[i].Yccd     = measure[i].Yccd;
    329     measureT[i].t        = measure[i].t;
    330     measureT[i].dt       = measure[i].dt;
    331     measureT[i].averef   = measure[i].averef;
    332     measureT[i].imageID  = measure[i].imageID;
    333     measureT[i].dbFlags  = measure[i].dbFlags;
    334     measureT[i].photcode = measure[i].photcode;
    335   }
    336 
    337   return (TRUE);
    338 }
    339 
    340 int free_tiny_values (Catalog *catalog) {
    341 
    342   if (catalog[0].averageT) free (catalog[0].averageT);
    343   if (catalog[0].measureT) free (catalog[0].measureT);
    344   return (TRUE);
    345 }
    346 
Note: See TracChangeset for help on using the changeset viewer.