IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Apr 4, 2011, 1:04:41 PM (15 years ago)
Author:
eugene
Message:

updates to pmPeak to better distinguish peak flux versions; updates to visualization; add bits for substantial suspect masking; consolidate assignment of source position and flux based on peak, moments, etc; improve footprint culling process; fix PSF_QF and PSF_QF_PERFECT calculations; fix source model chisq values

Location:
trunk/psModules
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules

    • Property svn:ignore

      •  

        old new  
        2828ChangeLog
        2929psmodules-*.tar.*
         30a.out.dSYM

  • trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c

    r30621 r31153  
    4949
    5050#include "pmSourceIO.h"
     51#include "pmSourceOutputs.h"
    5152
    5253// panstarrs-style FITS table output (header + table in 1st extension)
    … …  
    6263    psArray *table;
    6364    psMetadata *row;
    64     int i;
    65     psF32 *PAR, *dPAR;
    66     psEllipseAxes axes;
    67     psF32 xPos, yPos;
    68     psF32 xErr, yErr;
    69     psF32 errMag, chisq, apRadius;
    70     psS32 nPix, nDOF;
    7165
    7266    pmChip *chip = readout->parent->parent;
    73     pmFPA  *fpa  = chip->parent;
    74 
    75     bool status1 = false;
    76     bool status2 = false;
    77     float magOffset = NAN;
    78     float exptime   = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE");
    79     float zeropt    = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP");
    80     if (!isfinite(zeropt)) {
    81         zeropt    = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS");
    82     }
    83     if (status1 && status2 && (exptime > 0.0)) {
    84         magOffset = zeropt + 2.5*log10(exptime);
    85     }
    86     float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR");
    8767
    8868    // if the sequence is defined, write these in seq order; otherwise
    … …  
    9272        if (source->seq == -1) {
    9373          // let's write these out in S/N order
    94           sources = psArraySort (sources, pmSourceSortBySN);
     74          sources = psArraySort (sources, pmSourceSortByFlux);
    9575        } else {
    9676          sources = psArraySort (sources, pmSourceSortBySeq);
    … …  
    10080    table = psArrayAllocEmpty (sources->n);
    10181
     82    short nImageOverlap;
     83    float magOffset;
     84    float zeroptErr;
     85    float fwhmMajor;
     86    float fwhmMinor;
     87    pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
     88
    10289    // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
    10390    // by the time we call this function, all values should be assigned.  let's use asserts to be sure in some cases.
    104     for (i = 0; i < sources->n; i++) {
     91    for (int i = 0; i < sources->n; i++) {
    10592        pmSource *source = (pmSource *) sources->data[i];
    10693
    … …  
    113100        }
    114101
    115         // no difference between PSF and non-PSF model
    116         pmModel *model = source->modelPSF;
    117 
    118         if (model != NULL) {
    119             PAR = model->params->data.F32;
    120             dPAR = model->dparams->data.F32;
    121             xPos = PAR[PM_PAR_XPOS];
    122             yPos = PAR[PM_PAR_YPOS];
    123             if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) {
    124               xErr = dPAR[PM_PAR_XPOS];
    125               yErr = dPAR[PM_PAR_YPOS];
    126             } else {
    127               // in linear-fit mode, there is no error on the centroid
    128               xErr = source->peak->dx;
    129               yErr = source->peak->dy;
    130             }
    131             if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
    132                 axes = pmPSF_ModelToAxes (PAR, 20.0);
    133             } else {
    134                 axes.major = NAN;
    135                 axes.minor = NAN;
    136                 axes.theta = NAN;
    137             }
    138             chisq = model->chisq;
    139             nDOF = model->nDOF;
    140             nPix = model->nPix;
    141             apRadius = source->apRadius;
    142             errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
    143         } else {
    144             xPos = source->peak->xf;
    145             yPos = source->peak->yf;
    146             xErr = source->peak->dx;
    147             yErr = source->peak->dy;
    148             axes.major = NAN;
    149             axes.minor = NAN;
    150             axes.theta = NAN;
    151             chisq = NAN;
    152             nDOF = 0;
    153             nPix = 0;
    154             apRadius = NAN;
    155             errMag = NAN;
    156         }
    157 
    158         float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN;
    159         float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN;
    160         psS16 nImageOverlap = 1;
    161 
    162         psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
    163         float posAngle = 0.0;
    164         float pltScale = 0.0;
    165         pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos);
     102        // set the 'best' values for various output fields:
     103        pmSourceOutputs outputs;
     104        pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset);
     105
     106        pmSourceOutputsMoments moments;
     107        pmSourceOutputsSetMoments (&moments, source);
    166108
    167109        row = psMetadataAlloc ();
    168110        psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
    169         psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           xPos);
    170         psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           yPos);
    171         psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  xErr); // XXX this is only measured for non-linear fits
    172         psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  yErr); // XXX this is only measured for non-linear fits
    173         psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F32, "PSF RA coordinate (degrees)",                ptSky.r*PS_DEG_RAD);
    174         psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F32, "PSF DEC coordinate (degrees)",               ptSky.d*PS_DEG_RAD);
    175         psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         posAngle*PS_DEG_RAD);
    176         psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       pltScale*PS_DEG_RAD*3600.0);
     111        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           outputs.xPos);
     112        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           outputs.yPos);
     113        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  outputs.xErr); // XXX this is only measured for non-linear fits
     114        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  outputs.yErr); // XXX this is only measured for non-linear fits
     115        psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F32, "PSF RA coordinate (degrees)",                outputs.ra);
     116        psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F32, "PSF DEC coordinate (degrees)",               outputs.dec);
     117        psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         outputs.posAngle);
     118        psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       outputs.pltScale);
    177119        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfMag);
    178         psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        errMag);
     120        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->errMag);
    179121        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
    180         psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              apRadius);
    181         psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           peakMag);
    182         psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   calMag);
     122        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              outputs.apRadius);
     123        psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
     124        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   outputs.calMag);
    183125        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr);
    184126        psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
    185127        psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",                         source->skyErr);
    186128
    187         psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           chisq);
     129        psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           outputs.chisq);
    188130        psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA",        PS_DATA_F32, "Nsigma deviations from PSF to CF",           source->crNsigma);
    189131        psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA",       PS_DATA_F32, "Nsigma deviations from PSF to EXT",          source->extNsigma);
    190132
    191         psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     axes.major);
    192         psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     axes.minor);
    193         psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      axes.theta);
     133        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     outputs.psfMajor);
     134        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     outputs.psfMinor);
     135        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      outputs.psfTheta);
    194136        psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor",                source->pixWeightNotBad);
    195         psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         nDOF);
    196         psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    nPix);
    197 
    198         // distinguish moments measure from window vs S/N > XX ??
    199         float mxx = source->moments ? source->moments->Mxx : NAN;
    200         float mxy = source->moments ? source->moments->Mxy : NAN;
    201         float myy = source->moments ? source->moments->Myy : NAN;
    202         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                      mxx);
    203         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                      mxy);
    204         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                      myy);
     137        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         outputs.nDOF);
     138        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    outputs.nPix);
     139
     140        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                       moments.Mxx);
     141        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                       moments.Mxy);
     142        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                       moments.Myy);
    205143
    206144        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                      source->mode);
    … …  
    317255        // recreate the peak to match (xPos, yPos) +/- (xErr, yErr)
    318256        source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE);
    319         source->peak->flux = peakFlux;
     257        source->peak->rawFlux = peakFlux;
     258        source->peak->smoothFlux = peakFlux;
    320259        source->peak->xf   = PAR[PM_PAR_XPOS]; // more accurate position
    321260        source->peak->yf   = PAR[PM_PAR_YPOS]; // more accurate position
    322261        source->peak->dx   = dPAR[PM_PAR_XPOS];
    323262        source->peak->dy   = dPAR[PM_PAR_YPOS];
    324         if (isfinite (source->errMag) && (source->errMag > 0.0)) {
    325           source->peak->SN = 1.0 / source->errMag;
    326         } else {
    327           source->peak->SN = sqrt(source->peak->flux); // an alternate proxy: various functions sort by peak S/N
    328         }
    329263
    330264        source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF");
    … …  
    339273
    340274        source->moments = pmMomentsAlloc ();
     275        source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf
     276        source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf
     277
    341278        source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX");
    342279        source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY");
    … …  
    371308
    372309    // let's write these out in S/N order
    373     sources = psArraySort (sources, pmSourceSortBySN);
     310    sources = psArraySort (sources, pmSourceSortByFlux);
    374311
    375312    table = psArrayAllocEmpty (sources->n);
    … …  
    549486
    550487    // let's write these out in S/N order
    551     sources = psArraySort (sources, pmSourceSortBySN);
     488    sources = psArraySort (sources, pmSourceSortByFlux);
    552489
    553490    // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams
Note: See TracChangeset for help on using the changeset viewer.