IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 29491


Ignore:
Timestamp:
Oct 20, 2010, 8:55:04 AM (16 years ago)
Author:
eugene
Message:

update the source size analysis : CR now detected using M_minor and Kron mag; EXT now detected using only psf - kron mags; update visualization of source sizes; fix apresid visualization (kron vs apmag options); for stand-alone psphot, do not mask as SAT NAN pixels which are already masked; function for dynamic calculation of CR detection parameters

Location:
branches/eam_branches/ipp-20100823/psphot
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20100823/psphot/doc/notes.20101013.txt

    r29474 r29491  
     1
     220101019
     3
     4  We can identify CRs in the (M_minor, KronMag) plane as a fairly
     5  tight clump (M_minor ~< 1.0, KronMag ~< -8.0).  The best bounds of
     6  the clump vary somewhat from image to image.  I've added a dynamic
     7  CR limit function which looks for the peak in M_minor and in KronMag
     8  in those regions, then finds the valley floor.  It seems to work
     9  fairly well.  Outstanding issues:
     10
     11  * only adjust the limits once
     12  * record limits in the headers
     13  * record CR counts in the headers
     14  * mask the CRs.
    115
    21620101017

  • branches/eam_branches/ipp-20100823/psphot/src/psphot.h

    r29464 r29491  
    258258bool            psphotVisualShowFlags (psArray *sources);
    259259bool            psphotVisualShowSourceSize (pmReadout *readout, psArray *sources);
     260bool            psphotVisualPlotSourceSizeAlt (psMetadata *recipe, psMetadata *analysis, psArray *sources);
    260261bool            psphotVisualShowResidualImage (pmReadout *readout);
    261 bool            psphotVisualPlotApResid (psArray *sources, float mean, float error);
     262bool            psphotVisualPlotApResid (psArray *sources, float mean, float error, bool useApMag);
    262263bool            psphotVisualPlotChisq (psArray *sources);
    263264bool            psphotVisualPlotSourceSize (psMetadata *recipe, psMetadata *analysis, psArray *sources);

  • branches/eam_branches/ipp-20100823/psphot/src/psphotApResid.c

    r29464 r29491  
    349349    psFree (dMag);
    350350
    351     psphotVisualPlotApResid (sources, psf->ApResid, psf->dApResid);
     351    psphotVisualPlotApResid (sources, psf->ApResid, psf->dApResid, true);
    352352
    353353    return true;

  • branches/eam_branches/ipp-20100823/psphot/src/psphotImageLoop.c

    r28421 r29491  
    3030    // files associated with the science image
    3131    if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa in psphot.");
     32
     33    // select the appropriate recipe information
     34    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     35    psAssert (recipe, "missing recipe?");
     36
     37    psImageMaskType maskTest = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
    3238
    3339    // for psphot, we force data to be read at the chip level
    … …  
    95101                if (readout->mask) {
    96102                    psImageMaskType maskSat = pmConfigMaskGet("SAT", config); // Mask value for saturated pixels
    97                     if (!pmReadoutMaskNonfinite(readout, maskSat)) {
     103                    if (!pmReadoutMaskInvalid(readout, maskTest, maskSat)) {
    98104                        psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels.");
    99105                        psFree(view);

  • branches/eam_branches/ipp-20100823/psphot/src/psphotOutput.c

    r28013 r29491  
    270270    psMetadataItemSupplement (&status, header, analysis, "NDET_CR");
    271271
     272    psMetadataItemSupplement (&status, header, analysis, "PSPHOT.PSF.APRESID");
     273    psMetadataItemSupplement (&status, header, analysis, "PSPHOT.PSF.APRESID.SYSERR");
     274    psMetadataItemSupplement (&status, header, analysis, "PSPHOT.CR.MAX.SIZE");
     275    psMetadataItemSupplement (&status, header, analysis, "PSPHOT.CR.MAX.MAG");
     276
    272277    // sky background model statistics
    273278    psMetadataItemSupplement (&status, header, analysis, "MSKY_MN");

  • branches/eam_branches/ipp-20100823/psphot/src/psphotSourceSize.c

    r29459 r29491  
    11# include "psphotInternal.h"
    22# include <gsl/gsl_sf_gamma.h>
    3 
    4 # define KRON 1
    53
    64typedef struct {
    … …  
    1614    int grow;
    1715    int xtest, ytest;
    18     bool apply; // apply CR mask?
     16    bool applyCRmask; // apply CR mask?
     17    bool dynamicLimitsCR; // apply CR mask?
     18    float sizeLimitCR;
     19    float magLimitCR;
    1920} psphotSourceSizeOptions;
    2021
    2122// local functions:
    22 bool psphotSourceSizePSF (psphotSourceSizeOptions *options, psArray *sources, pmPSF *psf);
     23bool psphotSourceSizePSF (psphotSourceSizeOptions *options, pmReadout *readout, psArray *sources, pmPSF *psf, psMetadata *recipe);
     24bool psphotDynamicLimitsCR (psphotSourceSizeOptions *options, pmReadout *readout, psArray *sources, pmPSF *psf, psMetadata *recipe);
    2325bool psphotSourceClass (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options);
    2426bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options);
    25 bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options);
     27bool psphotSourceSelectCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options);
    2628bool psphotMaskCosmicRay (pmReadout *readout, pmSource *source, psImageMaskType maskVal);
    2729bool psphotMaskCosmicRayFootprintCheck (psArray *sources);
    2830int  psphotMaskCosmicRayConnected (int xPeak, int yPeak, psImage *mymask, psImage *myvar, psImage *edges, int binning, float sigma_thresh);
     31float psphotSourceSizeFindThreshold (psVector *value, psVector *mask, int maskValue, float minValue, float maxValue, float delta, float guess, float fraction);
    2932
    3033// we need to call this function after sources have been fitted to the PSF model and
    … …  
    111114    assert (status);
    112115
    113     // XXX recipe name is not great
     116    // location of a single test source
    114117    options.xtest = psMetadataLookupS32 (&status, recipe, "PSPHOT.CRMASK.XTEST");
    115118    options.ytest = psMetadataLookupS32 (&status, recipe, "PSPHOT.CRMASK.YTEST");
    … …  
    128131    }
    129132
    130     options.apply = psMetadataLookupBool(&status, recipe, "PSPHOT.CRMASK.APPLY"); // Growth size for CRs
     133    options.applyCRmask = psMetadataLookupBool(&status, recipe, "PSPHOT.CRMASK.APPLY"); // Growth size for CRs
    131134    if (!status) {
    132135        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "PSPHOT.CRMASK.APPLY is not defined.");
    … …  
    134137    }
    135138
    136     // We are using the value psfMag - 2.5*log10(moment.Sum) as a measure of the extendedness
    137     // of an object.  We need to model this distribution for the PSF stars before we can test
    138     // the significance for a specific object
    139     // XXX move this to the code that generates the PSF?
    140     // XXX store the results on pmPSF?
    141 
    142     // XXX this should only be done on the first pass (ie, if we have newSources or allSources?)
    143     if (getPSFsize) {
    144         psphotSourceSizePSF (&options, sources, psf);
    145     }
    146 
    147     // classify the sources based on ApResid and Moments (extended sources)
     139    // determine the distribution of (PSF_mag - KRON_mag) for the PSF sources (saved on readout->analysis)
     140    psphotSourceSizePSF (&options, readout, sources, psf, recipe);
     141
     142    // adjust the user-supplied limits based on the distribution of CRs in the (Mminor, mKron) space
     143    psphotDynamicLimitsCR(&options, readout, sources, psf, recipe);
     144
     145    // classify the sources based on ApResid and Moments
    148146    // NOTE: only sources not already measured !(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED)
    149147    psphotSourceClass(readout, sources, recipe, psf, &options);
    150148
    151     // NOTE: only sources not already measured !(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED)
    152     psphotSourceSizeCR (readout, sources, &options);
     149    // attempt to mask the candidate CRs; flag if CR nature is confirmed
     150    psphotSourceSelectCR(readout, sources, &options);
    153151
    154152    // XXX fix this (was source->n  - first)
    … …  
    156154
    157155    psphotVisualPlotSourceSize (recipe, readout->analysis, sources);
     156    psphotVisualPlotSourceSizeAlt (recipe, readout->analysis, sources);
    158157    psphotVisualShowSourceSize (readout, sources);
    159     psphotVisualPlotApResid (sources, options.ApResid, options.ApSysErr);
     158    psphotVisualPlotApResid (sources, options.ApResid, options.ApSysErr, false);
    160159    psphotVisualShowSatStars (recipe, psf, sources);
    161160
    … …  
    163162}
    164163
     164# define SAVE_PSF_OPTIONS(RO,OPT)                                       \
     165    /* save these on readout->analysis (a) for output to the header and (b) to prevent re-calculating in another pass */ \
     166    psMetadataAddF32(RO->analysis, PS_LIST_TAIL, "PSPHOT.PSF.APRESID", PS_META_REPLACE, "locus of PSF stars in PSF_MAG - KRON_MAG", OPT->ApResid); \
     167    psMetadataAddF32(RO->analysis, PS_LIST_TAIL, "PSPHOT.PSF.APRESID.SYSERR",  PS_META_REPLACE, "systematic error of PSF_MAG - KRON_MAG",  OPT->ApSysErr);
     168
    165169// model the apmifit distribution for the psf stars:
    166 bool psphotSourceSizePSF (psphotSourceSizeOptions *options, psArray *sources, pmPSF *psf) {
     170bool psphotSourceSizePSF (psphotSourceSizeOptions *options, pmReadout *readout, psArray *sources, pmPSF *psf, psMetadata *recipe) {
     171
     172    // We are using the value PSF_MAG - KRON_MAG as a measure of the extendedness of an object.
     173    // We need to model this distribution for the PSF stars before we can test the significance
     174    // for a specific object.
     175
     176    // NOTE: we require that objects have had moments measured, and we also require psfMags to
     177    // be calculated.  but, we do not require aperture mags or any other photometry values that
     178    // require pixel analysis.
     179
     180    bool status1 = false;
     181    bool status2 = false;
     182    options->ApResid = psMetadataLookupF32 (&status1, readout->analysis, "PSPHOT.PSF.APRESID");
     183    options->ApSysErr = psMetadataLookupF32 (&status2, readout->analysis, "PSPHOT.PSF.APRESID.SYSERR");
     184    psAssert ((status1 && status2) || (!status1 && !status2), "inconsistent record of PSF ApResid values");
     185
     186    // if they are saved on readout->analysis, we have already calculated these values
     187    if (status1 && status2) {
     188        return true;
     189    }
    167190
    168191    // select stats from the psf stars
    169     psVector *Ap = psVectorAllocEmpty (100, PS_TYPE_F32);
     192    psVector *ApOff = psVectorAllocEmpty (100, PS_TYPE_F32);
    170193    psVector *ApErr = psVectorAllocEmpty (100, PS_TYPE_F32);
    171194
    … …  
    173196    psImageMaskType maskVal = options->maskVal | options->markVal;
    174197
    175     // XXX  why PHOT_WEIGHT??
    176     pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT;
     198    // with PSFONLY, we do not need modify the pixels
     199    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_PSFONLY;
    177200
    178201    int num = 0;                        // Number of sources measured
    … …  
    182205        num++;
    183206
    184         // replace object in image
    185         if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
    186             pmSourceAdd (source, PM_MODEL_OP_FULL, options->maskVal);
    187         }
    188 
    189         // clear the mask bit and set the circular mask pixels
    190         psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));
    191         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal);
    192 
    193         // XXX can we test if psfMag is set and calculate only if needed?
    194207        pmSourceMagnitudes (source, psf, photMode, maskVal, markVal);
    195208
    196         // clear the mask bit
    197         psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));
    198 
    199         // re-subtract the object, leave local sky
    200         pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal);
    201 
    202         // XXX test: switch to kron flux
    203 # if (KRON)
    204         float apMag = -2.5*log10(source->moments->KronFlux);
    205 # else
    206         float apMag = -2.5*log10(source->moments->Sum);
    207 # endif
    208         float dMag = source->psfMag - apMag;
    209 
    210         psVectorAppend (Ap, dMag);
     209        float kMag = -2.5*log10(source->moments->KronFlux);
     210        float dMag = source->psfMag - kMag;
     211
     212        psVectorAppend (ApOff, dMag);
    211213        psVectorAppend (ApErr, source->errMag);
    212214    }
    213215    if (num == 0) {
    214         // Not raising an error, because errors aren't being checked elsewhere in this function
    215         psFree(Ap);
     216        // if we cannot determine the PSF distribution, call all objects PSFs...
     217        options->ApResid = NAN;
     218        options->ApSysErr = NAN;
     219        psFree(ApOff);
    216220        psFree(ApErr);
     221        SAVE_PSF_OPTIONS(readout, options);
    217222        return false;
    218223    }
    … …  
    221226    // psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);   
    222227    psStats *stats = psStatsAlloc(PS_STAT_CLIPPED_MEAN);   
    223     psVectorStats (stats, Ap, NULL, NULL, 0);
    224 
    225     psVector *dAp = psVectorAlloc (Ap->n, PS_TYPE_F32);
    226     for (int i = 0; i < Ap->n; i++) {
    227         dAp->data.F32[i] = Ap->data.F32[i] - stats->clippedMean;
     228    psVectorStats (stats, ApOff, NULL, NULL, 0);
     229
     230    psVector *dAp = psVectorAlloc (ApOff->n, PS_TYPE_F32);
     231    for (int i = 0; i < ApOff->n; i++) {
     232        dAp->data.F32[i] = ApOff->data.F32[i] - stats->clippedMean;
    228233    }
    229234
    230235    options->ApResid = stats->clippedMean;
    231236    options->ApSysErr = psVectorSystematicError(dAp, ApErr, 0.1);
    232     // XXX this is quite arbitrary...
    233     if (!isfinite(options->ApSysErr)) options->ApSysErr = 0.01;
     237
     238    // this is quite arbitrary... a large value means fewer things classified as extended.
     239    if (!isfinite(options->ApSysErr)) options->ApSysErr = 0.05;
    234240    psLogMsg ("psphot", PS_LOG_DETAIL, "psf - Sum: %f +/- %f\n", options->ApResid, options->ApSysErr);
    235241
    236     psFree (Ap);
     242    psFree (ApOff);
    237243    psFree (ApErr);
    238244    psFree (stats);
    239245    psFree (dAp);
    240246
     247    SAVE_PSF_OPTIONS(readout, options);
    241248    return true;
     249}
     250
     251# define SAVE_CR_OPTIONS(RO,OPT)                                        \
     252    /* save these on readout->analysis (a) for output to the header and (b) to prevent re-calculating in another pass */  \
     253    psMetadataAddF32(RO->analysis, PS_LIST_TAIL, "PSPHOT.CR.MAX.SIZE", PS_META_REPLACE, "dynamically-set minor axis size limit for cosmic rays", OPT->sizeLimitCR); \
     254    psMetadataAddF32(RO->analysis, PS_LIST_TAIL, "PSPHOT.CR.MAX.MAG",  PS_META_REPLACE, "dynamically-set kron magnitude limit for cosmic rays",  OPT->magLimitCR);
     255
     256// model the size and magnitude distribution of the Cosmic Rays
     257// ** CRs are reliably flagged by a combination on Mminor < X && mag (or flux) > Y
     258bool psphotDynamicLimitsCR (psphotSourceSizeOptions *options, pmReadout *readout, psArray *sources, pmPSF *psf, psMetadata *recipe) {
     259
     260    /* attempt to describe the CR sources:
     261       - input parameters are sizeLimit, magLimit
     262       - first try to refine the sizeLimit:
     263       -- select objects which meet the magLimit and exceed the sizeLimit by a factor of 1.5
     264       -- generate the histogram
     265       -- look for a peak in the histogram
     266       -- look for the min between valley and upper limit
     267       -- look for first bin within 5% of the valley floor after peak (new sizeLimit)
     268       
     269       - next try to refine the magLimit
     270       -- select objects which meet the sizeLimit and go fainter than the magLimit by 1.0 mag
     271       -- generate the histogram
     272       -- look for a peak in the histogram
     273       -- look for the min between valley and upper limit
     274       -- look for first bin within 5% of the valley floor after peak (new magLimit)
     275    */
     276
     277    bool status  = false;
     278    bool status1 = false;
     279    bool status2 = false;
     280    options->sizeLimitCR = psMetadataLookupF32 (&status1, readout->analysis, "PSPHOT.CR.MAX.SIZE");
     281    if (!status1) {
     282        options->sizeLimitCR = psMetadataLookupF32 (&status, recipe, "PSPHOT.CR.MAX.SIZE");
     283        if (!status) {
     284            options->sizeLimitCR = 1.0;
     285        }
     286    }
     287    options->magLimitCR = psMetadataLookupF32 (&status2, readout->analysis, "PSPHOT.CR.MAX.MAG");
     288    if (!status2) {
     289        options->magLimitCR = psMetadataLookupF32 (&status, recipe, "PSPHOT.CR.MAX.MAG");
     290        if (!status) {
     291            options->magLimitCR = -8.0;
     292        }
     293    }
     294    psAssert ((status1 && status2) || (!status1 && !status2), "inconsistent record of dynamic CR limits");
     295
     296    // if they are saved on readout->analysis, we have already calculated these values
     297    if (status1 && status2) {
     298        return true;
     299    }
     300
     301    // if we do not want to dynamically set these, save the user-supplied values and exit
     302    options->dynamicLimitsCR = psMetadataLookupBool (&status, recipe, "PSPHOT.CR.AUTOSCALE");
     303    if (!status) {
     304        options->dynamicLimitsCR = true;
     305    }
     306    if (!options->dynamicLimitsCR) {
     307        SAVE_CR_OPTIONS(readout, options); // macro defined above
     308        return true;
     309    }
     310
     311    psVector *minor = psVectorAllocEmpty (100, PS_TYPE_F32);
     312    psVector *mKron = psVectorAllocEmpty (100, PS_TYPE_F32);
     313    psVector *mask  = NULL;
     314
     315    psImageMaskType markVal = options->markVal;
     316    psImageMaskType maskVal = options->maskVal | options->markVal;
     317
     318    // with PSFONLY, we do not need modify the pixels
     319    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_PSFONLY;
     320
     321    // generate vectors for all the objects of possible interest (so we can just access those
     322    // vectors in the next sections)
     323    for (int i = 0; i < sources->n; i++) {
     324        pmSource *source = sources->data[i];
     325
     326        // XXX can we test if psfMag is set and calculate only if needed?
     327        pmSourceMagnitudes (source, psf, photMode, maskVal, markVal);
     328
     329        // convert to Mmaj, Mmin:
     330        psF32 Mxx = source->moments->Mxx;
     331        psF32 Myy = source->moments->Myy;
     332        psF32 Mxy = source->moments->Mxy;
     333
     334        float KronMag = -2.5*log10(source->moments->KronFlux);
     335
     336        float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy));
     337
     338        if (Mminor > options->sizeLimitCR * 1.5) continue;
     339        if (KronMag > options->magLimitCR + 2.5) continue;
     340
     341        psVectorAppend (mKron, KronMag);
     342        psVectorAppend (minor, Mminor);
     343    }
     344
     345    // if too few objects meet the criterion, give up..
     346    if (mKron->n < 50) goto escape;
     347
     348    // set distinct masks for (minor > sizeLimit) or (mKron > magLimit)
     349    mask = psVectorAlloc(mKron->n, PS_TYPE_U8);
     350    psVectorInit(mask, 0);
     351    for (int i = 0; i < mKron->n; i++) {
     352        if (mKron->data.F32[i] > options->magLimitCR) {
     353            mask->data.U8[i] |= 0x01;
     354        }
     355        if (minor->data.F32[i] > options->sizeLimitCR) {
     356            mask->data.U8[i] |= 0x02;
     357        }
     358    }
     359       
     360    float delta1 = PS_MAX(0.02, PS_MIN(0.2, 5.0 * 0.5 / minor->n));
     361    float newSizeLimit = psphotSourceSizeFindThreshold(minor, mask, 0x01, 0.0, options->sizeLimitCR * 1.5, delta1, options->sizeLimitCR, 0.05);
     362    if (isfinite(newSizeLimit)) {
     363        options->sizeLimitCR = newSizeLimit;
     364    }
     365
     366    float delta2 = PS_MAX(0.02, PS_MIN(0.2, 5.0 * 2.0 / minor->n));
     367    float newMagLimit = psphotSourceSizeFindThreshold(mKron, mask, 0x02, -15.0, options->magLimitCR + 2.5, delta2, options->magLimitCR, 0.05);
     368    if (isfinite(newMagLimit)) {
     369        options->magLimitCR = newMagLimit;
     370    }
     371
     372    psLogMsg ("psphot", PS_LOG_DETAIL, "CR limits : %f mag | %f pix^2\n", options->magLimitCR, options->sizeLimitCR);
     373
     374    psFree (mKron);
     375    psFree (minor);
     376    psFree (mask);
     377
     378    // save these on readout->analysis (a) for output to the header and (b) to prevent re-calculating in another pass
     379    SAVE_CR_OPTIONS(readout, options); // macro defined above
     380    return true;
     381
     382 escape:
     383    SAVE_CR_OPTIONS(readout, options); // macro defined above
     384    return false;
    242385}
    243386
    … …  
    249392    char regionName[64];
    250393
    251     psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4s %4s %4s %4s %4s %4s", "Npsf", "Next", "Nsat", "Ncr", "Nmiss", "Nskip");
     394    psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4s %4s %4s %4s %4s", "Npsf", "Next", "Nsat", "Ncr", "Nskip");
     395
     396    if (!psphotSourceClassRegion (NULL, &psfClump, sources, recipe, psf, options)) {
     397        psLogMsg ("psphot", 4, "Failed to determine source classification for full image\n");
     398    } else {
     399        psLogMsg ("psphot", 4, "source classification for full image\n");
     400    }
     401    return true;
    252402
    253403    int nRegions = psMetadataLookupS32 (&status, readout->analysis, "PSF.CLUMP.NREGIONS");
    … …  
    275425            continue;
    276426        }
     427        psLogMsg ("psphot", 4, "source classification for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1);
    277428        // psphotVisualPlotSourceSize (recipe, readout->analysis, sources);
    278429    }
    … …  
    281432}
    282433
    283 # define SIZE_SN_LIM 10
    284434bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options) {
    285435
    … …  
    291441    int Npsf  = 0;
    292442    int Ncr   = 0;
    293     int Nmiss = 0;
    294443    int Nskip = 0;
    295444
    296445    pmSourceMode noMoments = PM_SOURCE_MODE_MOMENTS_FAILURE | PM_SOURCE_MODE_SKYVAR_FAILURE | PM_SOURCE_MODE_SKY_FAILURE | PM_SOURCE_MODE_BELOW_MOMENTS_SN;
    297     pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT;
     446
     447    // request the pixWeight values as well as the magnitudes
     448    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT;
    298449
    299450    psImageMaskType markVal = options->markVal;
    … …  
    306457        // psfClumps are found for image subregions:
    307458        // skip sources not in this region
    308         if (source->peak->x <  region->x0) continue;
    309         if (source->peak->x >= region->x1) continue;
    310         if (source->peak->y <  region->y0) continue;
    311         if (source->peak->y >= region->y1) continue;
     459        if (region) {
     460            if (source->peak->x <  region->x0) continue;
     461            if (source->peak->x >= region->x1) continue;
     462            if (source->peak->y <  region->y0) continue;
     463            if (source->peak->y >= region->y1) continue;
     464        }
    312465
    313466        // skip source if it was already measured
    … …  
    321474            source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED;
    322475            psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n");
    323             continue;
    324         }
    325 
    326         // we are basically classifying by moments
     476            Nskip ++;
     477            continue;
     478        }
     479
     480        // we are classifying by moments and PSF_MAG - KRON_MAG
    327481        psAssert (source->moments, "why is this source missing moments?");
    328482        if (source->mode & noMoments) {
    … …  
    335489        psF32 Myy = source->moments->Myy;
    336490        psF32 Mxy = source->moments->Mxy;
    337 
    338         // replace object in image
     491        float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy));
     492
     493        // replace object in image to measure mag & psfWeights
    339494        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
    340495            pmSourceAdd (source, PM_MODEL_OP_FULL, options->maskVal);
    … …  
    345500        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal);
    346501
    347         // XXX can we test if psfMag is set and calculate only if needed?
    348502        pmSourceMagnitudes (source, psf, photMode, maskVal, markVal);
    349503
    … …  
    354508        pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal);
    355509
    356 # if (KRON)
    357         float apMag = -2.5*log10(source->moments->KronFlux);
    358 # else
    359         float apMag = -2.5*log10(source->moments->Sum);
    360 # endif
    361         float dMag = source->psfMag - apMag;
    362 
    363         // set nSigma to include both systematic and poisson error terms
    364         // XXX need to handle NAN psfMag and apMag (ie, skip?)
    365         // XXX the 'poisson error' contribution for size is probably wrong...
    366         // XXX add in a hard floor on the Ap Sys Err (to be a bit generous)
    367         float nSigmaMAG = (dMag - options->ApResid) / hypot(source->errMag, hypot(options->ApSysErr, 0.025));
    368         float nSigmaMXX = (Mxx - psfClump->X) / hypot(psfClump->dX, psfClump->X*psfClump->X*source->errMag);
    369         float nSigmaMYY = (Myy - psfClump->Y) / hypot(psfClump->dY, psfClump->Y*psfClump->Y*source->errMag);
    370 
    371         // XXX should I change psfClump to use minor, major?
    372 
    373         // fprintf (stderr, "%f %f : Mxx: %f, Myy: %f, dx: %f, dy: %f, psfMag: %f, apMag: %f, dMag: %f, errMag: %f, nSigmaMag: %f, nSigmaMxx: %f, nSigmaMyy: %f\n",
    374         // source->peak->xf, source->peak->yf, Mxx, Myy, source->peak->xf - source->moments->Mx, source->peak->yf - source->moments->My,
    375         // source->psfMag, apMag, dMag, source->errMag, nSigmaMAG, nSigmaMXX, nSigmaMYY);
    376 
    377         // XXX double check on ths stuff!! partially-masked sources are more likely to be mis-measured PSFs
    378         float sizeBias = 1.0;
    379         if (source->pixWeightNotBad < 0.9) {
    380             sizeBias = 3.0;
    381         }
    382         if (source->pixWeightNotPoor < 0.9) {
    383             sizeBias = 3.0;
    384         }
    385 
    386         float minMxx = psfClump->X - sizeBias*options->nSigmaMoments*psfClump->dX;
    387         float minMyy = psfClump->Y - sizeBias*options->nSigmaMoments*psfClump->dY;
    388 
    389         // include MAG, MXX, and MYY?
     510        float kMag = -2.5*log10(source->moments->KronFlux);
     511        float dMag = source->psfMag - kMag;
     512
     513        // sources without a valid magnitudes cannot have their size measured
     514        if (!isfinite(kMag) || !isfinite(source->psfMag)) {
     515            Nskip ++;
     516            continue;
     517        }
     518
     519        // set nSigmaMAG to include both systematic and poisson error terms.  we include a hard
     520        // floor on the Ap Sys Err (to be a bit generous).  XXX put the floor in the recipe...
     521        float nSigmaMAG = (dMag - options->ApResid) / hypot(source->errMag, hypot(options->ApSysErr, 0.02));
    390522        source->extNsigma = nSigmaMAG;
    391523
    … …  
    394526        // * CR & defect should have a faintess limit (min S/N)
    395527        // * SAT stars should not be faint, but defects may?
    396 
    397         // Anything within this region is a probably PSF-like object. Saturated stars may land
    398         // in this region, but are detected elsewhere on the basis of their peak value.
    399         bool isPSF = (fabs(nSigmaMAG) < options->nSigmaApResid) && (fabs(nSigmaMXX) < sizeBias*options->nSigmaMoments) && (fabs(nSigmaMYY) < sizeBias*options->nSigmaMoments);
    400         if (isPSF) {
    401           psTrace("psphotSourceClassRegion.PSF",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g PSF\t%g %g\n",
    402                   source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
    403                   options->nSigmaApResid,sizeBias*options->nSigmaMoments);
    404           source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    405           Npsf ++;
    406           continue;
    407         }
    408528
    409529        // Defects may not always match CRs from peak curvature analysis
    … …  
    413533        // (nSigmaMAG < -options->nSigmaApResid) ||
    414534
    415         // ** CRs are reliably flagged by a combination on Mminor < X && mag (or flux) > Y
    416 
    417         float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy));
    418         if ((Mminor < 1.0) && (apMag < -8.0)) {
    419             fprintf (stderr, "likely CR @ %f,%f\n", source->peak->xf, source->peak->yf);
    420         }
    421 
    422         // XXX do I need to find the Mminor, Mmajor distribution?
    423 
    424         bool isCR = isCR = (source->errMag < 1.0 / SIZE_SN_LIM) && ((Mxx < minMxx) || (Myy < minMyy));
     535        // saturated star (too many saturated pixels or peak above saturation limit).  These
     536        // may also be saturated galaxies, or just large saturated regions.
     537        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
     538            psTrace("psphotSourceClassRegion.SAT",4,"CLASS: %g %g\t%g %g\t%g %g\t%g %g\t%g SAT\n",
     539                    source->peak->xf, source->peak->yf, Mminor, kMag, dMag, nSigmaMAG, options->sizeLimitCR, options->magLimitCR, options->nSigmaApResid);
     540            source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
     541            Nsat ++;
     542            continue;
     543        }
     544
     545        // any sources missing a large fraction should just be treated as PSFs
     546        if ((source->pixWeightNotBad < 0.9) || (source->pixWeightNotPoor < 0.9)) {
     547            psTrace("psphotSourceClassRegion.PSF",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g PSF\t%g %g\n",
     548                    source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,kMag,dMag,nSigmaMAG,
     549                    options->nSigmaApResid,options->nSigmaMoments);
     550            source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
     551            Npsf ++;
     552            continue;
     553        }
     554
     555        // CRs are flagged by a combination on Mminor < options->sizeLimitCR && kmag < options->magLimitCR
     556        // NOTE: we only flag the CRs here; when we mask them we verify their CR nature (otherwise -> PSF)
     557        // bool isCR = (kMag < options->magLimitCR) && (Mminor < options->sizeLimitCR);
     558        bool isCR = (source->moments->SN > 7.0) && (Mminor < options->sizeLimitCR);
    425559        if (isCR) {
    426             psTrace("psphotSourceClassRegion.CR",4,"CLASS: %g %g %f\t%g %g  %g %g  %g %g\t%g %g\t%g CR\t%g %g\n",
    427                     source->peak->xf,source->peak->yf,source->pixWeightNotBad,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
    428                     options->nSigmaApResid,sizeBias*options->nSigmaMoments);
     560            psTrace("psphotSourceClassRegion.CR",4,"CLASS: %g %g\t%g %g\t%g %g\t%g %g\t%g CR\n",
     561                    source->peak->xf, source->peak->yf, Mminor, kMag, dMag, nSigmaMAG, options->sizeLimitCR, options->magLimitCR, options->nSigmaApResid);
    429562            source->mode |= PM_SOURCE_MODE_DEFECT;
    430563            source->tmpFlags |= PM_SOURCE_TMPF_SIZE_CR_CANDIDATE;
     564            source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    431565            Ncr ++;
    432566            continue;
    433567        }
    434568
    435         // saturated star (determined in PSF fit).  These may also be saturated galaxies, or
    436         // just large saturated regions.
    437         if (source->mode & PM_SOURCE_MODE_SATSTAR) {
    438             psTrace("psphotSourceClassRegion.SAT",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g SAT\t%g %g\n",
    439                     source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
    440                     options->nSigmaApResid,sizeBias*options->nSigmaMoments);
    441             source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    442             Nsat ++;
    443             continue;
    444         }
    445 
    446         // XXX allow the Mxx, Myy to be less than psfClump->X,Y (by some nSigma)?
    447         bool isEXT = (nSigmaMAG > options->nSigmaApResid) || (Mxx > minMxx) || (Myy > minMyy);
     569        // Likely extended source (PSF_MAG - KRON_MAG is larger than limit)
     570        bool isEXT = (nSigmaMAG > options->nSigmaApResid);
    448571        if (isEXT) {
    449           psTrace("psphotSourceClassRegion.EXT",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Ext\t%g %g\n",
    450                   source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
    451                   options->nSigmaApResid,sizeBias*options->nSigmaMoments);
    452 
     572            psTrace("psphotSourceClassRegion.EXT",4,"CLASS: %g %g\t%g %g\t%g %g\t%g %g\t%g EXT\n",
     573                    source->peak->xf, source->peak->yf, Mminor, kMag, dMag, nSigmaMAG, options->sizeLimitCR, options->magLimitCR, options->nSigmaApResid);
    453574            source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
    454575            source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    … …  
    456577            continue;
    457578        }
    458         psTrace("psphotSourceClassRegion.MISS",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Unk\t%g %g\n",
    459                 source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
    460                 options->nSigmaApResid,sizeBias*options->nSigmaMoments);
    461 
    462         // sources that reach here are probably too faint for a reasonable source size measurement
    463         // psWarning ("sourse size was missed for %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigmaMAG);
    464         source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    465         Nmiss ++;
    466     }
    467 
    468     psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4d %4d %4d %4d %4d %4d", Npsf, Next, Nsat, Ncr, Nmiss, Nskip);
     579
     580        // Everything else should just be treated as a PSF
     581        psTrace("psphotSourceClassRegion.PSF",4,"CLASS: %g %g\t%g %g\t%g %g\t%g %g\t%g PSF\n",
     582                source->peak->xf, source->peak->yf, Mminor, kMag, dMag, nSigmaMAG, options->sizeLimitCR, options->magLimitCR, options->nSigmaApResid);
     583        source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
     584        Npsf ++;
     585    }
     586
     587    psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4d %4d %4d %4d %4d", Npsf, Next, Nsat, Ncr, Nskip);
    469588
    470589    return true;
    … …  
    473592// given an object suspected to be a defect, generate a pixel mask using the Lapacian transform
    474593// if enough of the object is detected as 'sharp', consider the object a cosmic ray
    475 bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) {
     594bool psphotSourceSelectCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) {
    476595
    477596    psTimerStart ("psphot.cr");
    … …  
    481600        pmSource *source = sources->data[i];
    482601
    483         // skip source if it was already measured
    484         if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) {
    485             psTrace("psphot", 7, "Not calculating source size since it has already been measured\n");
    486             continue;
    487         }
    488 
    489602        // only check candidates marked above
    490603        if (!(source->tmpFlags & PM_SOURCE_TMPF_SIZE_CR_CANDIDATE)) {
    … …  
    493606        }
    494607
    495         // skip unless this source is thought to be a cosmic ray.  flag the detection and mask the pixels
    496         // XXX this may be degenerate with the above test
    497         if (!(source->mode & PM_SOURCE_MODE_DEFECT)) continue;
     608        // once we are here, remove the temporary flag
     609        source->mode &= ~PM_SOURCE_TMPF_SIZE_CR_CANDIDATE;
    498610
    499611        // Integer position of peak
    … …  
    519631        // XXX this is running slowly and is too agressive, but it more-or-less works
    520632        psTrace("psphot", 6, "mask cosmic ray at %f, %f\n", source->peak->xf, source->peak->yf);
    521         if (options->apply) {
     633        if (options->applyCRmask) {
    522634            psphotMaskCosmicRay(readout, source, options->crMask);
    523635        } else {
    … …  
    546658
    547659    // XXX test : save the mask image
    548     if (0) {
    549         psphotSaveImage (NULL, readout->mask,   "mask.fits");
     660    if (1) {
     661        psphotSaveImage (NULL, readout->mask,   "crmask.fits");
    550662    }
    551663
    … …  
    736848    if (nCRpix > 1) {
    737849        source->mode |= PM_SOURCE_MODE_CR_LIMIT;
    738         source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    739850    }
    740851    // fprintf (stderr, "CRMASK %d %d %d %d %d\n", peak->x, peak->y, dx, dy, nCRpix);
    … …  
    12311342}
    12321343
     1344float psphotSourceSizeFindThreshold (psVector *value, psVector *mask, int maskValue, float minValue, float maxValue, float delta, float guess, float fraction) {
     1345
     1346    // make a histogram of the sources with mKron < magLimit
     1347    int nValue = (int)((maxValue - minValue) / delta) + 1;
     1348
     1349    psVector *histogram = psVectorAlloc(nValue, PS_TYPE_S32);
     1350    psVectorInit (histogram, 0);
     1351
     1352    for (int i = 0; i < value->n; i++) {
     1353        if (mask->data.U8[i] & maskValue) continue;
     1354        int bin = (value->data.F32[i] - minValue) / delta;
     1355        if (bin < 0.0) continue;
     1356        if (bin > histogram->n - 1) continue;
     1357        histogram->data.S32[bin] ++;
     1358    }
     1359
     1360    // find the peak of this histogram, but stop search at guess
     1361    int last = (guess - minValue) / delta;
     1362    int nPeak = 0;
     1363    int vPeak = histogram->data.S32[0];
     1364    for (int i = 1; (i < last) && (i < histogram->n); i++) {
     1365        if (histogram->data.S32[i] < vPeak) continue;
     1366        nPeak = i;
     1367        vPeak = histogram->data.S32[i];
     1368    }
     1369
     1370    // start at the peak and find the valley between here and the end of the histogram
     1371    int nValley = nPeak;
     1372    int vValley = histogram->data.S32[nPeak];
     1373    for (int i = nPeak + 1; i < histogram->n; i++) {
     1374        if (histogram->data.S32[i] > vValley) continue;
     1375        nValley = i;
     1376        vValley = histogram->data.S32[i];
     1377    }
     1378
     1379    psLogMsg ("psphot", PS_LOG_MINUTIA, "CR limits threshold : peak %d @ %f, valley %d @ %f (%f sigma)\n",
     1380              vPeak, delta * nPeak + minValue, vValley, delta * nValley + minValue, vPeak / PS_MAX(1.0, sqrt(vValley)));
     1381
     1382    if (nValley == nPeak) {
     1383        psFree (histogram);
     1384        return NAN;
     1385    }
     1386
     1387    if (vPeak < 3.0*sqrt(vValley)) {
     1388        psFree (histogram);
     1389        return NAN;
     1390    }
     1391
     1392    /// search for the first bin after the peak with f < vValley + 0.05*(vPeak - vValley)
     1393    float vLimit = vValley + fraction*(vPeak - vValley);
     1394    int nLimit = nValley;
     1395    for (int i = nPeak; i < nValley; i++) {
     1396        if (histogram->data.S32[i] > vLimit) continue;
     1397        nLimit = i;
     1398        break;
     1399    }
     1400    float result = nLimit * delta + minValue;
     1401
     1402    psFree (histogram);
     1403   
     1404    return result;
     1405}

  • branches/eam_branches/ipp-20100823/psphot/src/psphotVisual.c

    r29458 r29491  
    471471    if (myKapa == -1) return false;
    472472
    473     KapaClearPlots (myKapa);
     473    KapaClearSections (myKapa);
    474474    KapaInitGraph (&graphdata);
    475475    KapaSetFont (myKapa, "courier", 14);
    … …  
    13981398    section.bg  = KapaColorByName ("none"); // XXX probably should be 'none'
    13991399
    1400     KapaClearPlots (myKapa);
     1400    KapaClearSections (myKapa);
    14011401    // first section : mag vs CR nSigma
    14021402    section.dx = 1.0;
    … …  
    16301630    if (myKapa == -1) return false;
    16311631
    1632     KapaClearPlots (myKapa);
     1632    KapaClearSections (myKapa);
    16331633    KapaInitGraph (&graphdata);
    16341634    KapaSetFont (myKapa, "courier", 14);
    … …  
    21342134}
    21352135
     2136bool PlotSourceSizeAltSetVectors(psVector *m, psVector *k, psVector *v1, psVector *v2, psVector *v3, pmSource *source) {
     2137
     2138    float Mxx = source->moments->Mxx;
     2139    float Myy = source->moments->Myy;
     2140    float Mxy = source->moments->Mxy;
     2141    float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy));
     2142    float KronMag = -2.5*log10(source->moments->KronFlux);
     2143   
     2144    psVectorAppend(m,  source->psfMag);
     2145    psVectorAppend(k,  KronMag);
     2146    psVectorAppend(v1, Mminor);
     2147    psVectorAppend(v2, source->psfMag - KronMag);
     2148    psVectorAppend(v3, source->extNsigma);
     2149    return true;
     2150}
     2151
     2152bool PlotSourceSizeAltAddPoints(Graphdata *graphdata, int myKapa, psVector *x, psVector *y, char *colorname, int ptype, float size) {
     2153
     2154    graphdata->color = KapaColorByName (colorname);
     2155    graphdata->ptype = ptype;
     2156    graphdata->size = size;
     2157    graphdata->style = 2;
     2158    KapaPrepPlot   (myKapa, x->n, graphdata);
     2159    KapaPlotVector (myKapa, x->n, x->data.F32, "x");
     2160    KapaPlotVector (myKapa, x->n, y->data.F32, "y");
     2161    return true;
     2162}
     2163
     2164bool psphotVisualPlotSourceSizeAlt (psMetadata *recipe, psMetadata *analysis, psArray *sources) {
     2165
     2166    Graphdata graphdata;
     2167    KapaSection section;
     2168
     2169    if (!pmVisualTestLevel("psphot.size", 2)) return true;
     2170
     2171    int myKapa = psphotKapaChannel (2);
     2172    if (myKapa == -1) return false;
     2173
     2174    KapaClearSections (myKapa);
     2175    KapaInitGraph (&graphdata);
     2176    KapaSetFont (myKapa, "courier", 14);
     2177
     2178    section.bg  = KapaColorByName ("none"); // XXX probably should be 'none'
     2179
     2180    // storage vectors for data to be plotted
     2181    psVector *SATm = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2182    psVector *SATk = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2183    psVector *SAT1 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2184    psVector *SAT2 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2185    psVector *SAT3 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2186
     2187    psVector *PSFm = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2188    psVector *PSFk = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2189    psVector *PSF1 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2190    psVector *PSF2 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2191    psVector *PSF3 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2192
     2193    psVector *EXTm = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2194    psVector *EXTk = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2195    psVector *EXT1 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2196    psVector *EXT2 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2197    psVector *EXT3 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2198
     2199    psVector *DEFm = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2200    psVector *DEFk = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2201    psVector *DEF1 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2202    psVector *DEF2 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2203    psVector *DEF3 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2204
     2205    psVector *BADm = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2206    psVector *BADk = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2207    psVector *BAD1 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2208    psVector *BAD2 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2209    psVector *BAD3 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2210
     2211    psVector *CRm = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2212    psVector *CRk = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2213    psVector *CR1 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2214    psVector *CR2 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2215    psVector *CR3 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2216
     2217    // int notPSF = PM_SOURCE_MODE_SATSTAR | PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT | PM_SOURCE_MODE_DEFECT;
     2218    int nSkip = 0;
     2219
     2220    // construct the vectors
     2221    for (int i = 0; i < sources->n; i++) {
     2222        pmSource *source = sources->data[i];
     2223        if (source->moments == NULL) {
     2224            nSkip ++;
     2225            continue;
     2226        }
     2227
     2228        bool found = false;
     2229
     2230        // only plot the measured sources...
     2231        if (!(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED)) {
     2232            nSkip ++;
     2233            continue;
     2234        }
     2235
     2236        // any sources missing a large fraction should just be treated as PSFs
     2237        if ((source->pixWeightNotBad < 0.9) || (source->pixWeightNotPoor < 0.9)) {
     2238            PlotSourceSizeAltSetVectors(BADm, BADk, BAD1, BAD2, BAD3, source);
     2239        }
     2240        if ((source->mode & PM_SOURCE_MODE_CR_LIMIT) || (source->tmpFlags & PM_SOURCE_TMPF_SIZE_CR_CANDIDATE)) {
     2241            PlotSourceSizeAltSetVectors(CRm, CRk, CR1, CR2, CR3, source);
     2242            found = true;
     2243        }
     2244        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
     2245            PlotSourceSizeAltSetVectors(SATm, SATk, SAT1, SAT2, SAT3, source);
     2246            found = true;
     2247        }
     2248        if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
     2249            PlotSourceSizeAltSetVectors(EXTm, EXTk, EXT1, EXT2, EXT3, source);
     2250            found = true;
     2251        }
     2252        if (source->mode & PM_SOURCE_MODE_DEFECT) {
     2253            PlotSourceSizeAltSetVectors(DEFm, DEFk, DEF1, DEF2, DEF3, source);
     2254            found = true;
     2255        }
     2256        if (!found) {
     2257            PlotSourceSizeAltSetVectors(PSFm, PSFk, PSF1, PSF2, PSF3, source);
     2258        }
     2259        // if (!(source->mode & notPSF) && !(source->tmpFlags & PM_SOURCE_TMPF_SIZE_CR_CANDIDATE)) {
     2260        //     PlotSourceSizeAltSetVectors(PSFm, PSFk, PSF1, PSF2, PSF3, source);
     2261        // }
     2262    }
     2263    // three sections: kronMag vs Mminor, psfMag vs psfMag - KronMag, psfMag vs extNsigma
     2264
     2265    // --- first section: kronMag vs Mminor ---
     2266    section.dx = 1.00;
     2267    section.dy = 0.33;
     2268    section.x  = 0.00;
     2269    section.y  = 0.00;
     2270    section.name = psStringCopy ("Mminor");
     2271    KapaSetSection (myKapa, &section);
     2272    psFree (section.name);
     2273
     2274    graphdata.color = KapaColorByName ("black");
     2275    graphdata.xmin = -17.1;
     2276    graphdata.xmax =  -6.9;
     2277    graphdata.ymin = -0.5;
     2278    graphdata.ymax = +7.1;
     2279    KapaSetLimits (myKapa, &graphdata);
     2280
     2281    graphdata.padXm = NAN;
     2282    graphdata.padYm = 5.0;
     2283    graphdata.padXp = NAN;
     2284    graphdata.padYp = NAN;
     2285    KapaBox (myKapa, &graphdata);
     2286    KapaSendLabel (myKapa, "kron mag", KAPA_LABEL_XM);
     2287    KapaSendLabel (myKapa, "M_minor| (pixels^2|)", KAPA_LABEL_YM);
     2288
     2289    PlotSourceSizeAltAddPoints(&graphdata, myKapa, PSFk, PSF1, "black", 0, 0.5);
     2290    PlotSourceSizeAltAddPoints(&graphdata, myKapa, EXTk, EXT1, "blue",  0, 0.5);
     2291    PlotSourceSizeAltAddPoints(&graphdata, myKapa,  CRk,  CR1, "red",   0, 0.5);
     2292    PlotSourceSizeAltAddPoints(&graphdata, myKapa, DEFk, DEF1, "red",   7, 1.0);
     2293    PlotSourceSizeAltAddPoints(&graphdata, myKapa, SATk, SAT1, "blue",  7, 1.2);
     2294    PlotSourceSizeAltAddPoints(&graphdata, myKapa, BADk, BAD1, "green", 7, 1.4);
     2295
     2296    // --- second section: dMag ----
     2297    section.dx = 1.00;
     2298    section.dy = 0.33;
     2299    section.x  = 0.00;
     2300    section.y  = 0.33;
     2301    section.name = psStringCopy ("dMag");
     2302    KapaSetSection (myKapa, &section);
     2303    psFree (section.name);
     2304
     2305    graphdata.color = KapaColorByName ("black");
     2306    graphdata.xmin = -17.1;
     2307    graphdata.xmax =  -6.9;
     2308    graphdata.ymin = -0.75;
     2309    graphdata.ymax = +1.50;
     2310    KapaSetLimits (myKapa, &graphdata);
     2311
     2312    graphdata.padXm = NAN;
     2313    graphdata.padYm = 5.0;
     2314    graphdata.padXp = 0.0;
     2315    graphdata.padYp = NAN;
     2316    strcpy (graphdata.labels, "0200");
     2317    KapaBox (myKapa, &graphdata);
     2318    KapaSendLabel (myKapa, "dMag", KAPA_LABEL_YM);
     2319
     2320    PlotSourceSizeAltAddPoints(&graphdata, myKapa, PSFm, PSF2, "black", 0, 0.5);
     2321    PlotSourceSizeAltAddPoints(&graphdata, myKapa, EXTm, EXT2, "blue",  0, 0.5);
     2322    PlotSourceSizeAltAddPoints(&graphdata, myKapa,  CRm,  CR2, "red",   0, 0.5);
     2323    PlotSourceSizeAltAddPoints(&graphdata, myKapa, DEFm, DEF2, "red",   7, 1.0);
     2324    PlotSourceSizeAltAddPoints(&graphdata, myKapa, SATm, SAT2, "blue",  7, 1.2);
     2325    PlotSourceSizeAltAddPoints(&graphdata, myKapa, BADm, BAD2, "green", 7, 1.4);
     2326
     2327    // --- third section: nSigma ---
     2328    section.dx = 1.00;
     2329    section.dy = 0.33;
     2330    section.x  = 0.00;
     2331    section.y  = 0.66;
     2332    section.name = psStringCopy ("nSigma");
     2333    KapaSetSection (myKapa, &section);
     2334    psFree (section.name);
     2335
     2336    graphdata.color = KapaColorByName ("black");
     2337    graphdata.xmin = -17.1;
     2338    graphdata.xmax = -6.9;
     2339    graphdata.ymin = -10.1;
     2340    graphdata.ymax = +10.1;
     2341    KapaSetLimits (myKapa, &graphdata);
     2342
     2343    graphdata.padXm = 0.0;
     2344    graphdata.padYm = 5.0;
     2345    graphdata.padXp = NAN;
     2346    graphdata.padYp = NAN;
     2347    strcpy (graphdata.labels, "0210");
     2348    KapaBox (myKapa, &graphdata);
     2349    KapaSendLabel (myKapa, "psf msg", KAPA_LABEL_XP);
     2350    KapaSendLabel (myKapa, "EXT nSigma", KAPA_LABEL_YM);
     2351
     2352    PlotSourceSizeAltAddPoints(&graphdata, myKapa, PSFm, PSF3, "black", 0, 0.5);
     2353    PlotSourceSizeAltAddPoints(&graphdata, myKapa, EXTm, EXT3, "blue",  0, 0.5);
     2354    PlotSourceSizeAltAddPoints(&graphdata, myKapa,  CRm,  CR3, "red",   0, 0.5);
     2355    PlotSourceSizeAltAddPoints(&graphdata, myKapa, DEFm, DEF3, "red",   7, 1.0);
     2356    PlotSourceSizeAltAddPoints(&graphdata, myKapa, SATm, SAT3, "blue",  7, 1.0);
     2357    PlotSourceSizeAltAddPoints(&graphdata, myKapa, BADm, BAD3, "green", 7, 1.4);
     2358
     2359    fprintf (stderr, "PSF: %ld, EXT: %ld, CR: %ld, DEF: %ld, SAT: %ld, TOTAL: %ld, nSkip: %d\n",
     2360             PSFm->n, EXTm->n, CRm->n, DEFm->n, SATm->n,
     2361             PSFm->n+ EXTm->n+ CRm->n+ DEFm->n+ SATm->n, nSkip);
     2362
     2363    psFree (SATk);
     2364    psFree (SATm);
     2365    psFree (SAT1);
     2366    psFree (SAT2);
     2367    psFree (SAT3);
     2368
     2369    psFree (EXTk);
     2370    psFree (EXTm);
     2371    psFree (EXT1);
     2372    psFree (EXT2);
     2373    psFree (EXT3);
     2374
     2375    psFree (PSFk);
     2376    psFree (PSFm);
     2377    psFree (PSF1);
     2378    psFree (PSF2);
     2379    psFree (PSF3);
     2380
     2381    psFree (DEFk);
     2382    psFree (DEFm);
     2383    psFree (DEF1);
     2384    psFree (DEF2);
     2385    psFree (DEF3);
     2386
     2387    psFree (BADk);
     2388    psFree (BADm);
     2389    psFree (BAD1);
     2390    psFree (BAD2);
     2391    psFree (BAD3);
     2392
     2393    psFree (CRk);
     2394    psFree (CRm);
     2395    psFree (CR1);
     2396    psFree (CR2);
     2397    psFree (CR3);
     2398
     2399    pmVisualAskUser(NULL);
     2400    return true;
     2401}
     2402
    21362403bool psphotVisualShowResidualImage (pmReadout *readout) {
    21372404
    … …  
    21472414}
    21482415
    2149 bool psphotVisualPlotApResid (psArray *sources, float mean, float error) {
     2416bool psphotVisualPlotApResid (psArray *sources, float mean, float error, bool useApMag) {
    21502417
    21512418    Graphdata graphdata;
    … …  
    21572424    if (myKapa == -1) return false;
    21582425
    2159     KapaClearPlots (myKapa);
     2426    KapaClearSections (myKapa);
    21602427    KapaInitGraph (&graphdata);
    21612428
    … …  
    21782445        if (!isfinite (source->psfMag)) continue;
    21792446
     2447        float dMag;
     2448        if (useApMag) {
     2449            dMag = source->apMag - source->psfMag;
     2450        } else {
     2451            float kMag = -2.5*log10(source->moments->KronFlux);
     2452            dMag = source->psfMag - kMag;
     2453        }
     2454
    21802455        x->data.F32[n] = source->psfMag;
    2181         y->data.F32[n] = source->apMag - source->psfMag;
     2456        y->data.F32[n] = dMag;
    21822457        dy->data.F32[n] = source->errMag;
    21832458        graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]);
    … …  
    22752550    if (myKapa == -1) return false;
    22762551
    2277     KapaClearPlots (myKapa);
     2552    KapaClearSections (myKapa);
    22782553    KapaInitGraph (&graphdata);
    22792554
Note: See TracChangeset for help on using the changeset viewer.