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_V3.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     psF32 *PAR, *dPAR;
    65     psEllipseAxes axes;
    66     psF32 xPos, yPos;
    67     psF32 xErr, yErr;
    68     psF32 chisq, apRadius;
    69     psS32 nPix, nDOF;
    7065
    7166    pmChip *chip = readout->parent->parent;
    72     pmFPA  *fpa  = chip->parent;
    73 
    74     bool status1 = false;
    75     bool status2 = false;
    76     float magOffset = NAN;
    77     float exptime   = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE");
    78     float zeropt    = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP");
    79     if (!isfinite(zeropt)) {
    80         zeropt    = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS");
    81     }
    82     if (status1 && status2 && (exptime > 0.0)) {
    83         magOffset = zeropt + 2.5*log10(exptime);
    84     }
    85     float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR");
    86 
    87     // we need a measure of the image quality (FWHM) for this image, in order to get the positional errors
    88     float fwhmMajor = psMetadataLookupF32(&status1, readout->analysis, "FWHM_MAJ");
    89     if (!status1) {
    90         fwhmMajor = psMetadataLookupF32(&status1, readout->analysis, "IQ_FW1");
    91         if (!status1) {
    92             fwhmMajor = 5.0; // XXX just a guess!
    93         }
    94     }
    95     float fwhmMinor = psMetadataLookupF32(&status1, readout->analysis, "FWHM_MIN");
    96     if (!status1) {
    97         fwhmMinor = psMetadataLookupF32(&status1, readout->analysis, "IQ_FW2");
    98         if (!status1) {
    99             fwhmMinor = 5.0; // XXX just a guess!
    100         }
    101     }
    10267
    10368    // if the sequence is defined, write these in seq order; otherwise
    … …  
    10772        if (source->seq == -1) {
    10873            // let's write these out in S/N order
    109             sources = psArraySort (sources, pmSourceSortBySN);
     74            sources = psArraySort (sources, pmSourceSortByFlux);
    11075        } else {
    11176            sources = psArraySort (sources, pmSourceSortBySeq);
    … …  
    11479
    11580    table = psArrayAllocEmpty (sources->n);
     81
     82    short nImageOverlap;
     83    float magOffset;
     84    float zeroptErr;
     85    float fwhmMajor;
     86    float fwhmMinor;
     87    pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
    11688
    11789    // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
    … …  
    128100        }
    129101
    130         // no difference between PSF and non-PSF model
    131         pmModel *model = source->modelPSF;
    132 
    133         if (model != NULL) {
    134             PAR = model->params->data.F32;
    135             dPAR = model->dparams->data.F32;
    136             xPos = PAR[PM_PAR_XPOS];
    137             yPos = PAR[PM_PAR_YPOS];
    138             if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) {
    139                 xErr = dPAR[PM_PAR_XPOS];
    140                 yErr = dPAR[PM_PAR_YPOS];
    141             } else {
    142                 xErr = fwhmMajor * source->errMag / 2.35;
    143                 yErr = fwhmMinor * source->errMag / 2.35;
    144             }
    145             if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
    146                 axes = pmPSF_ModelToAxes (PAR, 20.0);
    147             } else {
    148                 axes.major = NAN;
    149                 axes.minor = NAN;
    150                 axes.theta = NAN;
    151             }
    152             chisq = model->chisq;
    153             nDOF = model->nDOF;
    154             nPix = model->nPix;
    155             apRadius = source->apRadius;
    156         } else {
    157             bool useMoments = true;
    158             useMoments = (useMoments && source->moments);          // can't if there are no moments
    159             useMoments = (useMoments && source->moments->nPixels); // can't if the moments were not measured
    160             useMoments = (useMoments && !(source->mode && PM_SOURCE_MODE_MOMENTS_FAILURE)); // can't if the moments failed...
    161 
    162             if (useMoments) {
    163                 xPos = source->moments->Mx;
    164                 yPos = source->moments->My;
    165                 xErr = fwhmMajor * source->errMag / 2.35;
    166                 yErr = fwhmMinor * source->errMag / 2.35;
    167             } else {
    168                 xPos = source->peak->xf;
    169                 yPos = source->peak->yf;
    170                 xErr = source->peak->dx;
    171                 yErr = source->peak->dy;
    172             }
    173             axes.major = NAN;
    174             axes.minor = NAN;
    175             axes.theta = NAN;
    176             chisq = NAN;
    177             nDOF = 0;
    178             nPix = 0;
    179             apRadius = NAN;
    180         }
    181 
    182         float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN;
    183         float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN;
    184         psS16 nImageOverlap = 1;
    185 
    186         psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
    187         float posAngle = 0.0;
    188         float pltScale = 0.0;
    189         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);
    190108
    191109        row = psMetadataAlloc ();
    192110        psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
    193         psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           xPos);
    194         psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           yPos);
    195         psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  xErr);
    196         psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  yErr);
    197         psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         posAngle*PS_DEG_RAD);
    198         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);
     114        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  outputs.yErr);
     115        psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         outputs.posAngle);
     116        psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       outputs.pltScale);
    199117        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfMag);
    200118        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->errMag);
    … …  
    203121        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
    204122        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in reported aperture",             source->apMagRaw);
    205         psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              apRadius);
    206 
    207         psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   calMag);
     123        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              outputs.apRadius);
     124
     125        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   outputs.calMag);
    208126        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr);
    209127       
    210128        // NOTE: RA & DEC (both double) need to be on an 8-byte boundary...
    211         psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                ptSky.r*PS_DEG_RAD);
    212         psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               ptSky.d*PS_DEG_RAD);
    213 
    214         psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           peakMag);
     129        psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                outputs.ra);
     130        psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               outputs.dec);
     131
     132        psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
    215133        psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
    216134        psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",                         source->skyErr);
    217135
    218         psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           chisq);
     136        psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           outputs.chisq);
    219137        psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA",        PS_DATA_F32, "Nsigma deviations from PSF to CF",           source->crNsigma);
    220138        psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA",       PS_DATA_F32, "Nsigma deviations from PSF to EXT",          source->extNsigma);
    221139
    222         psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     axes.major);
    223         psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     axes.minor);
    224         psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      axes.theta);
     140        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     outputs.psfMajor);
     141        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     outputs.psfMinor);
     142        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      outputs.psfTheta);
    225143        psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor (bad)",          source->pixWeightNotBad);
    226144        psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",         source->pixWeightNotPoor);
    227         psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         nDOF);
    228         psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    nPix);
    229 
    230         // distinguish moments measure from window vs S/N > XX ??
    231         float Mxx = source->moments ? source->moments->Mxx : NAN;
    232         float Mxy = source->moments ? source->moments->Mxy : NAN;
    233         float Myy = source->moments ? source->moments->Myy : NAN;
    234 
    235         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                      Mxx);
    236         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                      Mxy);
    237         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                      Myy);
    238 
    239         float M_c3 = source->moments ? 1.0*source->moments->Mxxx - 3.0*source->moments->Mxyy : NAN;
    240         float M_s3 = source->moments ? 3.0*source->moments->Mxxy - 1.0*source->moments->Myyy : NAN;
    241         float M_c4 = source->moments ? 1.0*source->moments->Mxxxx - 6.0*source->moments->Mxxyy + 1.0*source->moments->Myyyy : NAN;
    242         float M_s4 = source->moments ? 4.0*source->moments->Mxxxy - 4.0*source->moments->Mxyyy : NAN;
    243 
    244         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C",      PS_DATA_F32, "third momemt cos theta",                    M_c3);
    245         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S",      PS_DATA_F32, "third momemt sin theta",                    M_s3);
    246         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                   M_c4);
    247         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                   M_s4);
    248 
    249         float Mrf  = source->moments ? source->moments->Mrf : NAN;
    250         float Mrh  = source->moments ? source->moments->Mrh : NAN;
    251         float Krf  = source->moments ? source->moments->KronFlux : NAN;
    252         float dKrf = source->moments ? source->moments->KronFluxErr : NAN;
    253 
    254         float Kinner = source->moments ? source->moments->KronFinner : NAN;
    255         float Kouter = source->moments ? source->moments->KronFouter : NAN;
    256 
    257         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                       Mrf);
    258         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",                        Mrh);
    259         psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",                     Krf);
    260         psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",                          dKrf);
    261 
    262         psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                     Kinner);
    263         psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                     Kouter);
    264 
    265         psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                     source->mode);
    266         psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                     source->mode2);
     145        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         outputs.nDOF);
     146        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    outputs.nPix);
     147
     148        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                       moments.Mxx);
     149        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                       moments.Mxy);
     150        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                       moments.Myy);
     151
     152        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C",      PS_DATA_F32, "third momemt cos theta",                     moments.M_c3);
     153        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S",      PS_DATA_F32, "third momemt sin theta",                     moments.M_s3);
     154        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                    moments.M_c4);
     155        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                    moments.M_s4);
     156
     157        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                        moments.Mrf);
     158        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",                         moments.Mrh);
     159        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Krf);
     160        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",                            moments.dKrf);
     161        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kinner);
     162        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kouter);
     163        // psMetadataAdd (row, PS_LIST_TAIL, "KRON_CORE_FLUX",   PS_DATA_F32, "Kron Flux (in 1.0 R1)",                      moments.KronCore);
     164        // psMetadataAdd (row, PS_LIST_TAIL, "KRON_CORE_ERROR",  PS_DATA_F32, "Kron Error (in 1.0 R1)",                     moments.KronCoreErr);
     165        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                      source->mode);
     166        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                      source->mode2);
    267167        psMetadataAdd (row, PS_LIST_TAIL, "PADDING2",         PS_DATA_S32, "more padding", 0);
    268168
    … …  
    384284        // recreate the peak to match (xPos, yPos) +/- (xErr, yErr)
    385285        source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE);
    386         source->peak->flux = peakFlux;
     286        source->peak->rawFlux = peakFlux;
     287        source->peak->smoothFlux = peakFlux;
    387288        source->peak->xf   = PAR[PM_PAR_XPOS]; // more accurate position
    388289        source->peak->yf   = PAR[PM_PAR_YPOS]; // more accurate position
    389290        source->peak->dx   = dPAR[PM_PAR_XPOS];
    390291        source->peak->dy   = dPAR[PM_PAR_YPOS];
    391         if (isfinite (source->errMag) && (source->errMag > 0.0)) {
    392           source->peak->SN = 1.0 / source->errMag;
    393         } else {
    394           source->peak->SN = sqrt(source->peak->flux); // an alternate proxy: various functions sort by peak S/N
    395         }
     292
     293        // we no longer sort by S/N, only flux
     294        // if (isfinite (source->errMag) && (source->errMag > 0.0)) {
     295        //   source->peak->SN = 1.0 / source->errMag;
     296        // } else {
     297        //   source->peak->SN = sqrt(source->peak->flux); // an alternate proxy: various functions sort by peak S/N
     298        // }
    396299
    397300        source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF");
    … …  
    407310
    408311        source->moments = pmMomentsAlloc ();
     312        source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf
     313        source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf
     314
    409315        source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX");
    410316        source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY");
    … …  
    477383
    478384    // let's write these out in S/N order
    479     sources = psArraySort (sources, pmSourceSortBySN);
     385    sources = psArraySort (sources, pmSourceSortByFlux);
    480386
    481387    table = psArrayAllocEmpty (sources->n);
    … …  
    649555
    650556    // let's write these out in S/N order
    651     sources = psArraySort (sources, pmSourceSortBySN);
     557    sources = psArraySort (sources, pmSourceSortByFlux);
    652558
    653559    // 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.