IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 33374


Ignore:
Timestamp:
Feb 26, 2012, 4:19:37 PM (14 years ago)
Author:
eugene
Message:

big cleanup of outlier rejections and flagging; measure ubercal distance

Location:
branches/eam_branches/ipp-20111122/Ohana/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20111122/Ohana/src/libdvo/include/dvo.h

    r33303 r33374  
    142142*/
    143143
    144 /* Average.code values -- these values are 32 bit (as of PS1_V1) */
     144/* Average.flags values -- these values are 32 bit (as of PS1_V1) */
    145145typedef enum {
    146146  ID_STAR_FEW     = 0x00000001, // used within relphot: skip star
    … …  
    169169} DVOAverageFlags;
    170170
     171/* Secfilt.flags values -- these values are 32 bit (as of PS1_V1) */
     172typedef enum {
     173  ID_SECF_STAR_FEW  = 0x00000001, // used within relphot: skip star
     174  ID_SECF_STAR_POOR = 0x00000002, // used within relphot: skip star
     175  ID_PHOTOM_PASS_0  = 0x00000100, // average magnitude calculated in 0th pass
     176  ID_PHOTOM_PASS_1  = 0x00000200, // average magnitude calculated in 1th pass
     177  ID_PHOTOM_PASS_2  = 0x00000400, // average magnitude calculated in 2th pass
     178  ID_PHOTOM_PASS_3  = 0x00000800, // average magnitude calculated in 3th pass
     179  ID_PHOTOM_PASS_4  = 0x00001000, // average magnitude calculated in 4th pass
     180} DVOSecfiltFlags;
     181
    171182/*** general dvo structures (internal use only / not IO) ***/
    172183

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

    r33364 r33374  
    216216float         getMmos             PROTO((off_t meas, int cat));
    217217float         getMrel             PROTO((Catalog *catalog, off_t meas, int cat));
     218short         getUbercalDist      PROTO((off_t meas, int cat));
    218219Image        *getimage            PROTO((off_t N));
    219220Image        *getimages           PROTO((off_t *N, off_t **LineNumber));
    … …  
    265266int           setMrel             PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr));
    266267void          setMrelFinal        PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr));
    267 int           setMrelOutput       PROTO((Catalog *catalog, int Ncatalog, int mark, FlatCorrectionTable *flatcorr));
     268int           setMrelOutput       PROTO((Catalog *catalog, int Ncatalog, int pass, FlatCorrectionTable *flatcorr));
    268269int           setMave             PROTO((Catalog *catalog, int Ncatalog));
    269270void          set_ZP              PROTO((double ZERO));

  • branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/ImageOps.c

    r33254 r33374  
    328328}
    329329
     330// returns image.Mcal - ff(x,y)
     331short getUbercalDist (off_t meas, int cat) {
     332
     333  off_t i;
     334  short distance;
     335
     336  i = MeasureToImage[cat][meas];
     337  if (i == -1) return (1000);
     338
     339  if (image[i].flags & IMAGE_BAD) return (1000); 
     340  distance = image[i].dummy3; // XXX dummy3 == ubercalDist
     341  return (distance);
     342}
     343
    330344float getMflat (off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog) {
    331345
    … …  
    393407    }     
    394408
    395     /* we optionally do not recalibrate images with UBERCAL zero points */
    396     if (KEEP_UBERCAL && (image[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
     409    // UBERCAL image: if this is an ubercal image, set minUbercalDist to 0:
     410    // we optionally do not recalibrate images with UBERCAL zero points
     411    if (image[i].flags & ID_IMAGE_PHOTOM_UBERCAL) {
     412      image[i].dummy3 = 0; // dummy3 = ubercalDist
     413      if (KEEP_UBERCAL) continue;
     414    }
     415
     416    int minUbercalDist = 1000;
    397417
    398418    N = 0;
    … …  
    436456      }
    437457
     458      minUbercalDist = MIN (catalog[c].secfilt[n*Nsecfilt].M_80, minUbercalDist); // M_80 == ubercalDist
     459
    438460      // fprintf (stderr, "%1d, %3d : %3d, %3d : %10.6f %10.6f : %6.3f  %6.3f  %6.3f  %6.3f  %6.3f\n", (int) i, (int) j, (int) c, (int) m, catalog[c].averageT[n].R, catalog[c].averageT[n].D, Msys, Mrel, Mmos, Mgrid, Mflat);
    439461
    … …  
    473495      image[i].Mcal = 0.0;
    474496    }
     497
     498    image[i].dummy3 = minUbercalDist;
    475499  }
    476500  free (list);

  • branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/StarOps.c

    r33373 r33374  
    6161
    6262      // XXX accumulate all secfilt values in a single pass?
     63
     64      int minUbercalDist = 1000;
    6365
    6466      int Ns;
    … …  
    102104          list[N] = Msys - Mcal - Mmos - Mgrid;
    103105
     106          int myUbercalDist = getUbercalDist(m,i);
     107          minUbercalDist = MIN(minUbercalDist, myUbercalDist);
     108
    104109          // dlist gives the error, which is used as the weight in WT_MEAN.
    105110          // we can modify the resulting weight in a few ways:
    … …  
    137142
    138143        liststats (list, dlist, N, &stats);
    139 
     144       
    140145        catalog[i].secfilt[Nsecfilt*j+Nsec].M  = stats.mean;
    141146        catalog[i].secfilt[Nsecfilt*j+Nsec].dM = stats.sigma;
    142147        catalog[i].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq + 1e-4) : NAN_S_SHORT;
     148       
     149        catalog[i].secfilt[Nsecfilt*j+Nsec].M_80 = minUbercalDist; // M_80 == ubercalDist
    143150      }
    144151    }
    … …  
    159166// * setMave updates 2MASS average flags
    160167// * setMave updates average EXT flags (PS1 and 2MASS)
    161 int setMrelOutput (Catalog *catalog, int Ncatalog, int mark, FlatCorrectionTable *flatcorr) {
     168int setMrelOutput (Catalog *catalog, int Ncatalog, int pass, FlatCorrectionTable *flatcorr) {
    162169
    163170  int i, N;
    … …  
    192199        if (catalog[i].found[Nsecfilt*j+Nsec]) continue; 
    193200
     201        int Ncode = 0;
     202
    194203        N = 0;
    195204        m = catalog[i].averageT[j].measureOffset;
    … …  
    200209          if (!code) continue;
    201210          if (code->equiv != thisCode) { continue; }
     211          Ncode ++;
    202212
    203213          if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue;
    204214
    205           // XXX allow REF stars (no Image Entry) to be included in the calculation this
    206           // should be optionally set, and should allow for REF stars to be downweighted by
    207           // more than their reported errors.  how such info is carried is unclear...
    208215          if (getImageEntry (m, i) < 0) {
    209216            Mcal = Mmos = Mgrid = 0;
    … …  
    247254        if (N < 1) continue;
    248255
     256        // force WT_MEAN or MEAN here?
     257
    249258        liststats (list, dlist, N, &stats);
    250         if (mark) catalog[i].found[Nsecfilt*j+Nsec] = TRUE;
     259        catalog[i].found[Nsecfilt*j+Nsec] = TRUE;
     260        switch (pass) {
     261          case 0:
     262            catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_0;
     263            break;
     264          case 1:
     265            catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_1;
     266            break;
     267          case 2:
     268            catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_2;
     269            break;
     270          case 3:
     271            catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_3;
     272            break;
     273          case 4:
     274            catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_4;
     275            break;
     276        }
    251277
    252278        /* use sigma or error in dM for output? */
    … …  
    254280        catalog[i].secfilt[Nsecfilt*j+Nsec].dM = MAX (stats.error, stats.sigma);
    255281        catalog[i].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT;
     282        catalog[i].secfilt[Nsecfilt*j+Nsec].Ncode = Ncode;
     283        catalog[i].secfilt[Nsecfilt*j+Nsec].Nused = stats.Nmeas;
    256284      }
    257285    }
    … …  
    540568      for (Ns = 0; Ns < Nphotcodes; Ns++) {
    541569       
    542         /* on final processing, skip stars already measured */
    543         if (final && catalog[i].found[j]) continue; 
    544 
    545570        int thisCode = photcodes[Ns][0].code;
    546571        int Nsec = GetPhotcodeNsec(thisCode);
    547572       
     573        /* on final processing, skip stars already measured */
     574        if (final && catalog[i].found[Nsecfilt*j + Nsec]) continue; 
     575
    548576        /* skip bad stars to prevent them from becoming good (on inner sample) */
    549577        if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; 
    … …  
    557585          if (ecode != thisCode) { continue; }
    558586
    559           // why am I skipping this??
    560           /* if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; */
     587          // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers
    561588          Mcal  = getMcal  (m, i, flatcorr, catalog);
    562589          if (isnan(Mcal)) { Ncal ++; continue; }
    … …  
    602629          if (ecode != thisCode) { continue; }
    603630
    604           /* if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; */
     631          // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers
    605632          Mcal  = getMcal  (m, i, flatcorr, catalog);
    606633          if (isnan(Mcal)) continue;
    … …  
    725752      list[n] = pow (10.0, 0.01*Xm);
    726753      dlist[n] = 1;
    727       if (list[n] > 100000) {
    728         fprintf (stderr, "!");
    729       }
    730754      n++;
    731755    }
    … …  
    765789      dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dM;
    766790      if (isnan(dM)) continue;
    767       if (dM > 100) {
    768         fprintf (stderr, ".");
    769       }
    770791      list[n] = dM;
    771792      dlist[n] = 1;

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

    r31450 r33374  
    239239    }
    240240    if (RESET) {
    241       // XXX assignMcal (&image[nimage], (double *) NULL, -1);
    242241      // XXX this needs to be thought through a bit more
    243       image[nimage].Mcal = 0.0;
    244       image[nimage].dMcal = NAN;
     242      if (!KEEP_UBERCAL) {
     243        image[nimage].Mcal = 0.0;
     244        image[nimage].dMcal = NAN;
     245      }
    245246      image[nimage].flags &= ~ID_IMAGE_PHOTOM_POOR;
     247      image[nimage].dummy3 = 1000; // dummy3 == ubercalDist
    246248    }
    247249    line_number[nimage] = i;

  • branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/setMrelFinal.c

    r33373 r33374  
    2121      ID_MEAS_AREA;       // detetion was outside acceptable area of device
    2222   
     23    // flags used by the photometry analysis (excluding UBERCAL)
     24    unsigned int PHOTOM_PASSES =
     25      ID_PHOTOM_PASS_0 | // average measured at pass 0
     26      ID_PHOTOM_PASS_1 | // average measured at pass 1
     27      ID_PHOTOM_PASS_2 | // average measured at pass 2
     28      ID_PHOTOM_PASS_3 | // average measured at pass 3
     29      ID_PHOTOM_PASS_4;  // average measured at pass 4
     30   
    2331    int Ns;
    2432    for (Ns = 0; Ns < Nphotcodes; Ns++) {
    … …  
    3139        catalog[0].secfilt[Nsecfilt*i+Nsec].dM = NAN;
    3240        catalog[0].secfilt[Nsecfilt*i+Nsec].Xm = NAN_S_SHORT;
     41        catalog[0].secfilt[Nsecfilt*j+Nsec].Ncode = 0;
     42        catalog[0].secfilt[Nsecfilt*j+Nsec].Nused = 0;
     43
     44        catalog[0].secfilt[Nsecfilt*i+Nsec].flags &= ~PHOTOM_PASSES;
    3345
    3446        m = catalog[0].average[i].measureOffset;
    … …  
    6476  }
    6577
     78  clean_measures (catalog, 1, TRUE, flatcorr);    /* mark outliers ID_MEAS_POOR_PHOTOM */
    6679  for (i = 0; i < 5; i++) {
    6780    skip_measurements (catalog, i, flatcorr);       /* set ID_MEAS_SKIP for measures to be skipped */
    68     setMrelOutput  (catalog, 1, FALSE, flatcorr);   /* set Mrel using all measures not skipped */
    69     clean_measures (catalog, 1, TRUE, flatcorr);    /* mark outliers ID_MEAS_POOR_PHOTOM */
    70     setMrelOutput  (catalog, 1, TRUE, flatcorr);    /* set Mrel using remaining measures */
     81    setMrelOutput  (catalog, 1, i, flatcorr);    /* set Mrel using remaining measures */
    7182  }
    7283  setMcalOutput (catalog, 1, flatcorr);
    73 
    74   int Nsecfilt = GetPhotcodeNsecfilt ();
    7584
    7685  /* clear ID_STAR_POOR, ID_STAR_FEW values before writing ??? */
    7786  /* ID_MEAS_NOCAL is an internal bit, so it should be cleared */
     87  // XXX eventually: move these internal flags to an internal-only catalog->FOO element (like found)
    7888  for (i = 0; i < catalog[0].Naverage; i++) {
    7989    catalog[0].average[i].flags &= ~ID_STAR_FEW;
    … …  
    95105// acceptance on each pass. The rejections at the different passes are:
    96106
    97 // 0 : only use measurements thought to be not extended and dM < SIGMA_LIM (REMOVED)
    98 // 1 : use extended and low-sig detections (REMOVED)
    99 // 3 : accept the measurements marked as outliers; accept images marked as outliers
     107// 0 : only use measurements thought to be GOOD (photflags not POOR or BAD)
     108// 1 : accept measurements thought to be POOR based on photflags
     109// 2 : accept the measurements marked as outliers; accept images marked as outliers
     110// 3 : accept measurements thought to be BAD based on photflags
    100111// 4 : accept the measurements outside of the instrumental magnitude limits (eg, SAT)
    101112
    … …  
    106117
    107118  off_t i, k, m;
    108   int ecode, found, Ns;
    109119  off_t Ntot, Ntry, Nkeep, Nskip;
    110120  float mag;
    111121
     122  int Nsecfilt = GetPhotcodeNsecfilt ();
     123
    112124  Ntot = Ntry = Nskip = Nkeep = 0;
    113125
    114126  /* allow measures from images marked POOR and FEW */
    115   if (pass >= 3) IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL;
     127  if (pass >= 2) IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL;
    116128 
    117129  /* allow measures marked as outliers (POOR) and off image region (AREA) */
    118   if (pass >= 3) MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_SKIP_PHOTOM;
    119 
    120   /* mark measures which should be ignored on second pass */
     130  if (pass >= 2) MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_SKIP_PHOTOM;
     131
     132  /* mark measures which should be ignored on this pass */
    121133  for (i = 0; i < catalog[0].Naverage; i++) {
    122134    Ntot += catalog[0].average[i].Nmeasure;
    123135
    124     // this function does not check secfilt entries separately, but there is no harm in
    125     // twiddling the so skip a star here only we can only
    126     // have we gotten
    127 
    128     if (catalog[0].found[i]) continue;
    129   XXX : this should be done by secfilt.  think about this...
    130 
    131     m = catalog[0].average[i].measureOffset;
    132     for (k = 0; k < catalog[0].average[i].Nmeasure; k++, m++) {
    133       Ntry++;
    134 
    135       /* clear SKIP for all measures at first */
    136       catalog[0].measureT[m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
    137 
    138       /** never use these measurements (wrong photcode, bad time range) */
    139       /* skipped via NOCAL, don't mark as skipped */
    140       ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[m].photcode);
    141       found = FALSE;
    142       for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) {
    143         if (ecode == photcodes[Ns][0].code) found = TRUE;
     136    // mark measurements for each secfilt separately
     137    int Ns;
     138    for (Ns = 0; Ns < Nphotcodes; Ns++) {
     139     
     140      int thisCode = photcodes[Ns][0].code;
     141      int Nsec = GetPhotcodeNsec(thisCode);
     142     
     143      /* star/photcodes already calibrated */
     144      if (catalog[0].found[Nsecfilt*i+Nsec]) continue; 
     145
     146      m = catalog[0].average[i].measureOffset;
     147      for (k = 0; k < catalog[0].average[i].Nmeasure; k++, m++) {
     148        Ntry++;
     149
     150        // skip measurements not related to this photcode
     151        PhotCode *code = GetPhotcodebyCode (catalog[0].measure[m].photcode);
     152        if (!code) continue;
     153        if (code->equiv != thisCode) continue;
     154
     155        // skip measurements by time range (mark as skipped for this secfilt)
     156        // XXX note that this is a bit dangerous : some objects may never get calibrated
     157        if (TimeSelect) {
     158          if (catalog[0].measure[m].t < TSTART) goto skip;
     159          if (catalog[0].measure[m].t > TSTOP) goto skip;
     160        }
     161       
     162        // we now have a measurement relavant to this photcode (and time range)
     163
     164        // clear SKIP for all measures at first
     165        catalog[0].measureT[m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
     166
     167        // skip measurements from BAD images and mosaics
     168        // do NOT skip measurements without a matching image
     169        if (isnan(getMcal (m, 0, flatcorr, catalog))) goto skip;
     170        if (isnan(getMmos (m, 0))) goto skip;
     171
     172        // PASS 4 : skip measurements by inst mag limit
     173        if ((pass < 4) && ImagSelect) {
     174          mag = PhotInst (&catalog[0].measure[m]);
     175          if (mag < ImagMin) goto skip;
     176          if (mag > ImagMax) goto skip;
     177        }
     178
     179        // PASS 3 : accept bad measurements (eg, SAT, CR)
     180        if (pass < 3) {
     181          if (catalog[0].measure[m].photFlags & code->photomBadMask) goto skip;
     182          if ((catalog[i].measure[m].photcode > 10000) && (catalog[i].measure[m].photcode < 10500)) {
     183            if (catalog[i].measure[m].psfQual < 0.85) goto skip;
     184          }
     185        }
     186       
     187        // PASS 2 : internal outliers accepted
     188
     189        // PASS 1 : accept poor measurements as well (eg, POOR FIT, etc)
     190        if ((pass < 1) && (catalog[0].measure[m].photFlags & code->photomPoorMask)) goto skip;
     191
     192        // PASS 0 : skip poor, outlier, bad measurements as well (eg, POOR FIT, etc)
     193
     194        Nkeep ++;
     195       
     196        continue;
     197
     198      skip:
     199        catalog[0].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM;
     200        catalog[0].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM;
     201        Nskip ++;
    144202      }
    145       if (!found) continue;
    146 
    147       /* skip measurements by time range */
    148       if (TimeSelect) {
    149         if (catalog[0].measure[m].t < TSTART) goto skip;
    150         if (catalog[0].measure[m].t > TSTOP) goto skip;
    151       }
    152 
    153       /* skip measurements with sat. dophot values */
    154       // XXX we need to fix this: use the photFlags to ID saturation...
    155       // if ((pass < 4) && (catalog[0].measure[m].dophot == 10)) goto skip;
    156 
    157       /* skip measurements from BAD images and mosaics */
    158       /* do NOT skip measurements without a matching image */
    159       if (isnan(getMcal (m, 0, flatcorr, catalog))) goto skip;
    160       if (isnan(getMmos (m, 0))) goto skip;
    161 
    162       /* skip measurements by inst mag limit */
    163       if ((pass < 4) && ImagSelect) {
    164         mag = PhotInst (&catalog[0].measure[m]);
    165         if (mag < ImagMin) goto skip;
    166         if (mag > ImagMax) goto skip;
    167       }
    168 
    169       // d1 = (catalog[0].measure[m].dophot == 1);
    170       // d2 = (catalog[0].measure[m].dophot == 2);
    171       // fprintf (stderr, "fix the dophot / type selections!");
    172       // abort ();
    173 
    174       /* skip measurements by measurement error */
    175       // if ((pass < 2) && !(d1 || d2)) goto skip;
    176 
    177       /* if ((pass < 1) && (catalog[0].measure[m].dM > SIGMA_LIM)) goto skip; */
    178       // if ((pass < 1) && !d1) goto skip;
    179       Nkeep ++;
    180 
    181       continue;
    182 
    183     skip:
    184       catalog[0].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM;
    185       catalog[0].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM;
    186       Nskip ++;
    187203    }
    188204  }
Note: See TracChangeset for help on using the changeset viewer.