IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 29548


Ignore:
Timestamp:
Oct 25, 2010, 3:16:58 PM (16 years ago)
Author:
eugene
Message:

merge changes from eam_branches/ipp-20100823

Location:
trunk/psphot
Files:
23 edited
5 copied

Legend:

Unmodified
Added
Removed

  • trunk/psphot/doc/notes.20100715.txt

    r29004 r29548  
     1
     22010.08.24
     3
     4  Remaining work to be done:
     5
     6  * test and turn on CR masking
     7  * double-check source size analysis:
     8    * are we doing the right thing with sources flagged as bad in 'rough'?
     9    * are we doing the right thing for SAT, CR, EXT, PSF after SourceSize?
     10  * convert CR_SIGMA and EXT_SIGMA to probabilities
     11  * example results of model fits vs fake inputs
     12  * what is the right solution for the stack PSF photometry?
     13    * Nigel sees 4-6%
     14    * we see 1.9% at the worst
     15    * non-poisson errors are clearly better -- check on the chisq values
     16    * why is QGAUSS so good for the one example (better for poisson, if trendy)
    117
    2182010.08.12

  • trunk/psphot/src/psphot.h

    r29004 r29548  
    167167
    168168// used by psphotFindDetections
    169 psImage        *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, const int pass, psImageMaskType maskVal);
     169psImage        *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, psImageMaskType maskVal);
    170170psArray        *psphotFindPeaks (psImage *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax);
    171171bool            psphotFindFootprints (pmDetections *detections, psImage *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int pass, psImageMaskType maskVal);
    … …  
    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);
    264265bool            psphotVisualShowPetrosians (psArray *sources);
    265266bool            psphotVisualEraseOverlays (int channel, char *overlay);
     267bool            psphotVisualClose(void);
    266268
    267269bool psphotPetrosian (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal);

  • trunk/psphot/src/psphotApResid.c

    r29004 r29548  
    339339    psFree (dMag);
    340340
    341     psphotVisualPlotApResid (sources, psf->ApResid, psf->dApResid);
     341    psphotVisualPlotApResid (sources, psf->ApResid, psf->dApResid, true);
    342342
    343343    return true;

  • trunk/psphot/src/psphotCleanup.c

    r28013 r29548  
    1414    psFree (config);
    1515
     16    psphotVisualClose();
    1617    psMemCheckCorruption (stderr, true);
    1718    pmModelClassCleanup ();
    … …  
    1920    pmConceptsDone ();
    2021    pmConfigDone ();
     22    pmVisualCleanup ();
    2123    psLibFinalize();
    2224    // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psphot");

  • trunk/psphot/src/psphotExtendedSourceAnalysisByObject.c

    r28013 r29548  
    3838    // which extended source analyses should we perform?
    3939    bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
    40     bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
    4140    bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
    42     bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
    4341
    4442    // number of images used to define sources
    … …  
    117115
    118116            // if we request any of these measurements, we require the radial profile
    119             if (doPetrosian || doIsophotal || doAnnuli || doKron) {
     117            if (doPetrosian || doAnnuli) {
    120118                if (!psphotRadialProfile (source, recipe, skynoise, maskVal)) {
    121119                    // all measurements below require the radial profile; skip them all
    122120                    // re-subtract the object, leave local sky
    123                     psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
     121                    psTrace ("psphot", 5, "FAILED radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
    124122                    pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    125123                    continue;
     124                } else {
     125                    psTrace ("psphot", 5, "measured radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
     126                    Nannuli ++;
     127                    source->mode |= PM_SOURCE_MODE_RADIAL_FLUX;
    126128                }
    127                 source->mode |= PM_SOURCE_MODE_RADIAL_FLUX;
    128129            }
    129130
    … …  
    138139                }
    139140            }
    140 
    141141
    142142            // re-subtract the object, leave local sky

  • trunk/psphot/src/psphotExtendedSourceFits.c

    r29139 r29548  
    8585        return true;
    8686    }
     87
     88    psphotInitRadiusEXT (recipe, readout);
    8789
    8890    // validate the model entries
    … …  
    167169            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for NplainPass
    168170
    169             if (!psThreadJobAddPending(job)) {
     171            if (false && !psThreadJobAddPending(job)) {
    170172                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    171173                psFree(AnalysisRegion);
    172174                return false;
    173             }
     175            } else {
     176                if (!psphotExtendedSourceFits_Threaded(job)) {
     177                    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     178                    psFree(AnalysisRegion);
     179                    return false;
     180                }
     181                psScalar *scalar = NULL;
     182                scalar = job->args->data[7];
     183                Next += scalar->data.S32;
     184                scalar = job->args->data[8];
     185                Nconvolve += scalar->data.S32;
     186                scalar = job->args->data[9];
     187                NconvolvePass += scalar->data.S32;
     188                scalar = job->args->data[10];
     189                Nplain += scalar->data.S32;
     190                scalar = job->args->data[11];
     191                NplainPass += scalar->data.S32;
     192                psFree(job);
     193            }
    174194        }
    175195
    … …  
    272292            fprintf (stderr, "skipping (1) %f, %f\n", source->peak->xf, source->peak->yf);
    273293            pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    274             // XXX raise an error flag of some kind
     294            // XXX raise an error of some kind
    275295            continue;
    276296        }
    … …  
    296316            if (!modelFluxStart) {
    297317                fprintf (stderr, "skipping (3) %f, %f\n", source->peak->xf, source->peak->yf);
    298                 // XXX raise an error flag of some kind?
     318                pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     319                // XXX raise an error of some kind?
    299320                continue;
    300321            }

  • trunk/psphot/src/psphotFindDetections.c

    r29004 r29548  
    8585
    8686    // generate the smoothed significance image
    87     psImage *significance = psphotSignificanceImage (readout, recipe, pass, maskVal);
     87    psImage *significance = psphotSignificanceImage (readout, recipe, maskVal);
    8888
    8989    // display the significance image

  • trunk/psphot/src/psphotFitSourcesLinear.c

    r29004 r29548  
    104104    float MIN_VALID_FLUX = psMetadataLookupF32(&status, recipe, "PSF_FIT_MIN_VALID_FLUX");
    105105    if (!status) {
    106         MIN_VALID_FLUX = 1e-8;
     106        MIN_VALID_FLUX = 0.0;
    107107    }
    108108    float MAX_VALID_FLUX = psMetadataLookupF32(&status, recipe, "PSF_FIT_MAX_VALID_FLUX");

  • trunk/psphot/src/psphotFitSourcesLinearStack.c

    r29004 r29548  
    3535    bool CONSTANT_PHOTOMETRIC_WEIGHTS = psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
    3636    psAssert (status, "You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS");
     37
     38    float MIN_VALID_FLUX = psMetadataLookupF32(&status, recipe, "PSF_FIT_MIN_VALID_FLUX");
     39    if (!status) {
     40        MIN_VALID_FLUX = 0.0;
     41    }
     42    float MAX_VALID_FLUX = psMetadataLookupF32(&status, recipe, "PSF_FIT_MAX_VALID_FLUX");
     43    if (!status) {
     44        MAX_VALID_FLUX = 1e+8;
     45    }
    3746
    3847    // XXX store a local static array of covar factors for each of the images (by image ID)
    … …  
    123132
    124133    psSparseConstraint constraint;
    125     constraint.paramMin   = 0.0;
    126     constraint.paramMax   = 1e8;
    127     constraint.paramDelta = 1e8;
     134    constraint.paramMin   = MIN_VALID_FLUX;
     135    constraint.paramMax   = MAX_VALID_FLUX;
     136    constraint.paramDelta = 1e7;
    128137
    129138    // solve for normalization terms (need include local sky?)

  • trunk/psphot/src/psphotGuessModels.c

    r29004 r29548  
    176176
    177177        // We have two options to get a guess for the object position: the position from the
    178         // peak and the position from the moments.  Use the peak position if (a) there are no
    179         // moments and (b) the sources is not saturated
    180 
     178        // peak and the position from the moments.  Use the peak position if there are no
     179        // moments
     180
     181# if (0)
     182        bool useMoments = true;
     183# else
    181184        bool useMoments = false;
    182185        useMoments = (source->mode & PM_SOURCE_MODE_SATSTAR);  // we only want to try if SATSTAR is set, but..
     186# endif
     187
    183188        useMoments = (useMoments && source->moments);          // can't if there are no moments
    184189        useMoments = (useMoments && source->moments->nPixels); // can't if the moments were not measured

  • trunk/psphot/src/psphotImageLoop.c

    r28421 r29548  
    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);

  • trunk/psphot/src/psphotModelTest.c

    r28013 r29548  
    226226    // measure the source mags
    227227    pmSourcePhotometryModel (&fitMag, model);
    228     pmSourcePhotometryAper  (&obsMag, model, source->pixels, source->maskObj, maskVal);
     228    pmSourcePhotometryAper  (&obsMag, NULL, NULL, model, source->pixels, source->variance, source->maskObj, maskVal);
    229229    fprintf (stderr, "ap: %f, fit: %f, apmifit: %f, nIter: %d\n", obsMag, fitMag, obsMag - fitMag, model->nIter);
    230230

  • trunk/psphot/src/psphotOutput.c

    r28013 r29548  
    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");

  • trunk/psphot/src/psphotSignificanceImage.c

    r29004 r29548  
    44// (S/N)^2.  If FWMH_X,Y have been recorded, use them, otherwise use PEAKS_SMOOTH_SIGMA for the
    55// smoothing kernel.
    6 psImage *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, const int pass, psImageMaskType maskVal) {
     6psImage *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, psImageMaskType maskVal) {
    77
    88    float SIGMA_SMTH, NSIGMA_SMTH;
    … …  
    6767    // XXX change these to recipe value checks
    6868    if (psTraceGetLevel("psphot") > 5) {
     69        static int pass = 0;
    6970        char name[64];
    7071        sprintf (name, "imsmooth.v%d.fits", pass);
    … …  
    7273        sprintf (name, "wtsmooth.v%d.fits", pass);
    7374        psphotSaveImage(NULL, smooth_wt, name);
     75        pass ++;
    7476    }
    7577
    … …  
    108110    if (psTraceGetLevel("psphot") > 5) {
    109111        char name[64];
     112        static int pass = 0;
    110113        sprintf (name, "snsmooth.v%d.fits", pass);
    111114        psphotSaveImage (NULL, smooth_im, name);
     115        pass ++;
    112116    }
    113117

  • trunk/psphot/src/psphotSourceFits.c

    r29017 r29548  
    551551
    552552        if (!psphotFitSersicIndexPCM (pcm, readout, source, fitOptions, maskVal, markVal, psfSize)) {
     553            psFree (pcm);
    553554            model->flags |= PM_MODEL_STATUS_BADARGS;
    554555            return model;

  • trunk/psphot/src/psphotSourceMatch.c

    r28013 r29548  
    228228            pmSource *source = pmSourceAlloc();
    229229            source->imageID = index;
     230            source->mode2 |= PM_SOURCE_MODE2_MATCHED; // source is generated based on another image
    230231
    231232            // add the peak

  • trunk/psphot/src/psphotSourceSize.c

    r29017 r29548  
    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);
     32
     33// save test output images?
     34# define DUMPPICS 0
    2935
    3036// we need to call this function after sources have been fitted to the PSF model and
    … …  
    111117    assert (status);
    112118
    113     // XXX recipe name is not great
     119    // location of a single test source
    114120    options.xtest = psMetadataLookupS32 (&status, recipe, "PSPHOT.CRMASK.XTEST");
    115121    options.ytest = psMetadataLookupS32 (&status, recipe, "PSPHOT.CRMASK.YTEST");
    … …  
    128134    }
    129135
    130     options.apply = psMetadataLookupBool(&status, recipe, "PSPHOT.CRMASK.APPLY"); // Growth size for CRs
     136    options.applyCRmask = psMetadataLookupBool(&status, recipe, "PSPHOT.CRMASK.APPLY"); // Growth size for CRs
    131137    if (!status) {
    132138        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "PSPHOT.CRMASK.APPLY is not defined.");
    … …  
    134140    }
    135141
    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)
     142    // determine the distribution of (PSF_mag - KRON_mag) for the PSF sources (saved on readout->analysis)
     143    psphotSourceSizePSF (&options, readout, sources, psf, recipe);
     144
     145    // adjust the user-supplied limits based on the distribution of CRs in the (Mminor, mKron) space
     146    psphotDynamicLimitsCR(&options, readout, sources, psf, recipe);
     147
     148    // classify the sources based on ApResid and Moments
    148149    // NOTE: only sources not already measured !(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED)
    149150    psphotSourceClass(readout, sources, recipe, psf, &options);
    150151
    151     // NOTE: only sources not already measured !(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED)
    152     psphotSourceSizeCR (readout, sources, &options);
     152    // attempt to mask the candidate CRs; flag if CR nature is confirmed
     153    psphotSourceSelectCR(readout, sources, &options);
    153154
    154155    // XXX fix this (was source->n  - first)
    … …  
    156157
    157158    psphotVisualPlotSourceSize (recipe, readout->analysis, sources);
     159    psphotVisualPlotSourceSizeAlt (recipe, readout->analysis, sources);
    158160    psphotVisualShowSourceSize (readout, sources);
    159     psphotVisualPlotApResid (sources, options.ApResid, options.ApSysErr);
     161    psphotVisualPlotApResid (sources, options.ApResid, options.ApSysErr, false);
    160162    psphotVisualShowSatStars (recipe, psf, sources);
    161163
    … …  
    163165}
    164166
     167# define SAVE_PSF_OPTIONS(RO,OPT)                                       \
     168    /* save these on readout->analysis (a) for output to the header and (b) to prevent re-calculating in another pass */ \
     169    psMetadataAddF32(RO->analysis, PS_LIST_TAIL, "PSPHOT.PSF.APRESID", PS_META_REPLACE, "locus of PSF stars in PSF_MAG - KRON_MAG", OPT->ApResid); \
     170    psMetadataAddF32(RO->analysis, PS_LIST_TAIL, "PSPHOT.PSF.APRESID.SYSERR",  PS_META_REPLACE, "systematic error of PSF_MAG - KRON_MAG",  OPT->ApSysErr);
     171
    165172// model the apmifit distribution for the psf stars:
    166 bool psphotSourceSizePSF (psphotSourceSizeOptions *options, psArray *sources, pmPSF *psf) {
     173bool psphotSourceSizePSF (psphotSourceSizeOptions *options, pmReadout *readout, psArray *sources, pmPSF *psf, psMetadata *recipe) {
     174
     175    // We are using the value PSF_MAG - KRON_MAG as a measure of the extendedness of an object.
     176    // We need to model this distribution for the PSF stars before we can test the significance
     177    // for a specific object.
     178
     179    // NOTE: we require that objects have had moments measured, and we also require psfMags to
     180    // be calculated.  but, we do not require aperture mags or any other photometry values that
     181    // require pixel analysis.
     182
     183    bool status1 = false;
     184    bool status2 = false;
     185    options->ApResid = psMetadataLookupF32 (&status1, readout->analysis, "PSPHOT.PSF.APRESID");
     186    options->ApSysErr = psMetadataLookupF32 (&status2, readout->analysis, "PSPHOT.PSF.APRESID.SYSERR");
     187    psAssert ((status1 && status2) || (!status1 && !status2), "inconsistent record of PSF ApResid values");
     188
     189    // if they are saved on readout->analysis, we have already calculated these values
     190    if (status1 && status2) {
     191        return true;
     192    }
    167193
    168194    // select stats from the psf stars
    169     psVector *Ap = psVectorAllocEmpty (100, PS_TYPE_F32);
     195    psVector *ApOff = psVectorAllocEmpty (100, PS_TYPE_F32);
    170196    psVector *ApErr = psVectorAllocEmpty (100, PS_TYPE_F32);
    171197
    … …  
    173199    psImageMaskType maskVal = options->maskVal | options->markVal;
    174200
    175     // XXX  why PHOT_WEIGHT??
    176     pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT;
     201    // with PSFONLY, we do not need modify the pixels
     202    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_PSFONLY;
    177203
    178204    int num = 0;                        // Number of sources measured
    … …  
    182208        num++;
    183209
    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?
    194210        pmSourceMagnitudes (source, psf, photMode, maskVal, markVal);
    195211
    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);
     212        float kMag = -2.5*log10(source->moments->KronFlux);
     213        float dMag = source->psfMag - kMag;
     214
     215        psVectorAppend (ApOff, dMag);
    211216        psVectorAppend (ApErr, source->errMag);
    212217    }
    213218    if (num == 0) {
    214         // Not raising an error, because errors aren't being checked elsewhere in this function
    215         psFree(Ap);
     219        // if we cannot determine the PSF distribution, call all objects PSFs...
     220        options->ApResid = NAN;
     221        options->ApSysErr = NAN;
     222        psFree(ApOff);
    216223        psFree(ApErr);
     224        SAVE_PSF_OPTIONS(readout, options);
    217225        return false;
    218226    }
    219227
    220228    // model the distribution as a mean or median value and a systematic error from that value:
    221     psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
    222     psVectorStats (stats, Ap, NULL, NULL, 0);
    223 
    224     psVector *dAp = psVectorAlloc (Ap->n, PS_TYPE_F32);
    225     for (int i = 0; i < Ap->n; i++) {
    226         dAp->data.F32[i] = Ap->data.F32[i] - stats->robustMedian;
    227     }
    228 
    229     options->ApResid = stats->robustMedian;
    230     options->ApSysErr = psVectorSystematicError(dAp, ApErr, 0.05);
    231     // XXX this is quite arbitrary...
    232     if (!isfinite(options->ApSysErr)) options->ApSysErr = 0.01;
     229    // psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);   
     230    psStats *stats = psStatsAlloc(PS_STAT_CLIPPED_MEAN);   
     231    psVectorStats (stats, ApOff, NULL, NULL, 0);
     232
     233    psVector *dAp = psVectorAlloc (ApOff->n, PS_TYPE_F32);
     234    for (int i = 0; i < ApOff->n; i++) {
     235        dAp->data.F32[i] = ApOff->data.F32[i] - stats->clippedMean;
     236    }
     237
     238    options->ApResid = stats->clippedMean;
     239    options->ApSysErr = psVectorSystematicError(dAp, ApErr, 0.1);
     240
     241    // this is quite arbitrary... a large value means fewer things classified as extended.
     242    if (!isfinite(options->ApSysErr)) options->ApSysErr = 0.05;
    233243    psLogMsg ("psphot", PS_LOG_DETAIL, "psf - Sum: %f +/- %f\n", options->ApResid, options->ApSysErr);
    234244
    235     psFree (Ap);
     245    psFree (ApOff);
    236246    psFree (ApErr);
    237247    psFree (stats);
    238248    psFree (dAp);
    239249
     250    SAVE_PSF_OPTIONS(readout, options);
    240251    return true;
     252}
     253
     254# define SAVE_CR_OPTIONS(RO,OPT)                                        \
     255    /* save these on readout->analysis (a) for output to the header and (b) to prevent re-calculating in another pass */  \
     256    psMetadataAddF32(RO->analysis, PS_LIST_TAIL, "PSPHOT.CR.MAX.SIZE", PS_META_REPLACE, "dynamically-set minor axis size limit for cosmic rays", OPT->sizeLimitCR); \
     257    psMetadataAddF32(RO->analysis, PS_LIST_TAIL, "PSPHOT.CR.MAX.MAG",  PS_META_REPLACE, "dynamically-set kron magnitude limit for cosmic rays",  OPT->magLimitCR);
     258
     259// model the size and magnitude distribution of the Cosmic Rays
     260// ** CRs are reliably flagged by a combination on Mminor < X && mag (or flux) > Y
     261bool psphotDynamicLimitsCR (psphotSourceSizeOptions *options, pmReadout *readout, psArray *sources, pmPSF *psf, psMetadata *recipe) {
     262
     263    /* attempt to describe the CR sources:
     264       - input parameters are sizeLimit, magLimit
     265       - first try to refine the sizeLimit:
     266       -- select objects which meet the magLimit and exceed the sizeLimit by a factor of 1.5
     267       -- generate the histogram
     268       -- look for a peak in the histogram
     269       -- look for the min between valley and upper limit
     270       -- look for first bin within 5% of the valley floor after peak (new sizeLimit)
     271       
     272       - next try to refine the magLimit
     273       -- select objects which meet the sizeLimit and go fainter than the magLimit by 1.0 mag
     274       -- generate the histogram
     275       -- look for a peak in the histogram
     276       -- look for the min between valley and upper limit
     277       -- look for first bin within 5% of the valley floor after peak (new magLimit)
     278    */
     279
     280    bool status  = false;
     281    bool status1 = false;
     282    bool status2 = false;
     283    options->sizeLimitCR = psMetadataLookupF32 (&status1, readout->analysis, "PSPHOT.CR.MAX.SIZE");
     284    if (!status1) {
     285        options->sizeLimitCR = psMetadataLookupF32 (&status, recipe, "PSPHOT.CR.MAX.SIZE");
     286        if (!status) {
     287            options->sizeLimitCR = 1.0;
     288        }
     289    }
     290    options->magLimitCR = psMetadataLookupF32 (&status2, readout->analysis, "PSPHOT.CR.MAX.MAG");
     291    if (!status2) {
     292        options->magLimitCR = psMetadataLookupF32 (&status, recipe, "PSPHOT.CR.MAX.MAG");
     293        if (!status) {
     294            options->magLimitCR = -8.0;
     295        }
     296    }
     297    psAssert ((status1 && status2) || (!status1 && !status2), "inconsistent record of dynamic CR limits");
     298
     299    // if they are saved on readout->analysis, we have already calculated these values
     300    if (status1 && status2) {
     301        return true;
     302    }
     303
     304    // if we do not want to dynamically set these, save the user-supplied values and exit
     305    options->dynamicLimitsCR = psMetadataLookupBool (&status, recipe, "PSPHOT.CR.AUTOSCALE");
     306    if (!status) {
     307        options->dynamicLimitsCR = true;
     308    }
     309    if (!options->dynamicLimitsCR) {
     310        SAVE_CR_OPTIONS(readout, options); // macro defined above
     311        return true;
     312    }
     313
     314    psVector *minor = psVectorAllocEmpty (100, PS_TYPE_F32);
     315    psVector *mKron = psVectorAllocEmpty (100, PS_TYPE_F32);
     316    psVector *mask  = NULL;
     317
     318    psImageMaskType markVal = options->markVal;
     319    psImageMaskType maskVal = options->maskVal | options->markVal;
     320
     321    // with PSFONLY, we do not need modify the pixels
     322    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_PSFONLY;
     323
     324    // generate vectors for all the objects of possible interest (so we can just access those
     325    // vectors in the next sections)
     326    for (int i = 0; i < sources->n; i++) {
     327        pmSource *source = sources->data[i];
     328
     329        // XXX can we test if psfMag is set and calculate only if needed?
     330        pmSourceMagnitudes (source, psf, photMode, maskVal, markVal);
     331
     332        // convert to Mmaj, Mmin:
     333        psF32 Mxx = source->moments->Mxx;
     334        psF32 Myy = source->moments->Myy;
     335        psF32 Mxy = source->moments->Mxy;
     336
     337        float KronMag = -2.5*log10(source->moments->KronFlux);
     338
     339        float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy));
     340
     341        if (Mminor > options->sizeLimitCR * 1.5) continue;
     342        if (KronMag > options->magLimitCR + 2.5) continue;
     343
     344        psVectorAppend (mKron, KronMag);
     345        psVectorAppend (minor, Mminor);
     346    }
     347
     348    // if too few objects meet the criterion, give up..
     349    if (mKron->n < 50) goto escape;
     350
     351    // set distinct masks for (minor > sizeLimit) or (mKron > magLimit)
     352    mask = psVectorAlloc(mKron->n, PS_TYPE_U8);
     353    psVectorInit(mask, 0);
     354    for (int i = 0; i < mKron->n; i++) {
     355        if (mKron->data.F32[i] > options->magLimitCR) {
     356            mask->data.U8[i] |= 0x01;
     357        }
     358        if (minor->data.F32[i] > options->sizeLimitCR) {
     359            mask->data.U8[i] |= 0x02;
     360        }
     361    }
     362       
     363    float delta1 = PS_MAX(0.02, PS_MIN(0.2, 5.0 * 0.5 / minor->n));
     364    float newSizeLimit = psphotSourceSizeFindThreshold(minor, mask, 0x01, 0.0, options->sizeLimitCR * 1.5, delta1, options->sizeLimitCR, 0.05);
     365    if (isfinite(newSizeLimit)) {
     366        options->sizeLimitCR = newSizeLimit;
     367    }
     368
     369    float delta2 = PS_MAX(0.02, PS_MIN(0.2, 5.0 * 2.0 / minor->n));
     370    float newMagLimit = psphotSourceSizeFindThreshold(mKron, mask, 0x02, -15.0, options->magLimitCR + 2.5, delta2, options->magLimitCR, 0.05);
     371    if (isfinite(newMagLimit)) {
     372        options->magLimitCR = newMagLimit;
     373    }
     374
     375    psLogMsg ("psphot", PS_LOG_DETAIL, "CR limits : %f mag | %f pix^2\n", options->magLimitCR, options->sizeLimitCR);
     376
     377    psFree (mKron);
     378    psFree (minor);
     379    psFree (mask);
     380
     381    // save these on readout->analysis (a) for output to the header and (b) to prevent re-calculating in another pass
     382    SAVE_CR_OPTIONS(readout, options); // macro defined above
     383    return true;
     384
     385 escape:
     386    SAVE_CR_OPTIONS(readout, options); // macro defined above
     387    return false;
    241388}
    242389
    … …  
    248395    char regionName[64];
    249396
    250     psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4s %4s %4s %4s %4s %4s", "Npsf", "Next", "Nsat", "Ncr", "Nmiss", "Nskip");
     397    psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4s %4s %4s %4s %4s", "Npsf", "Next", "Nsat", "Ncr", "Nskip");
     398
     399    if (!psphotSourceClassRegion (NULL, &psfClump, sources, recipe, psf, options)) {
     400        psLogMsg ("psphot", 4, "Failed to determine source classification for full image\n");
     401    } else {
     402        psLogMsg ("psphot", 4, "source classification for full image\n");
     403    }
     404    return true;
    251405
    252406    int nRegions = psMetadataLookupS32 (&status, readout->analysis, "PSF.CLUMP.NREGIONS");
    … …  
    274428            continue;
    275429        }
     430        psLogMsg ("psphot", 4, "source classification for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1);
    276431        // psphotVisualPlotSourceSize (recipe, readout->analysis, sources);
    277432    }
    … …  
    280435}
    281436
    282 # define SIZE_SN_LIM 10
    283437bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options) {
    284438
    … …  
    290444    int Npsf  = 0;
    291445    int Ncr   = 0;
    292     int Nmiss = 0;
    293446    int Nskip = 0;
    294447
    295448    pmSourceMode noMoments = PM_SOURCE_MODE_MOMENTS_FAILURE | PM_SOURCE_MODE_SKYVAR_FAILURE | PM_SOURCE_MODE_SKY_FAILURE | PM_SOURCE_MODE_BELOW_MOMENTS_SN;
    296     pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT;
     449
     450    // request the pixWeight values as well as the magnitudes
     451    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT;
    297452
    298453    psImageMaskType markVal = options->markVal;
    … …  
    305460        // psfClumps are found for image subregions:
    306461        // skip sources not in this region
    307         if (source->peak->x <  region->x0) continue;
    308         if (source->peak->x >= region->x1) continue;
    309         if (source->peak->y <  region->y0) continue;
    310         if (source->peak->y >= region->y1) continue;
     462        if (region) {
     463            if (source->peak->x <  region->x0) continue;
     464            if (source->peak->x >= region->x1) continue;
     465            if (source->peak->y <  region->y0) continue;
     466            if (source->peak->y >= region->y1) continue;
     467        }
    311468
    312469        // skip source if it was already measured
    … …  
    320477            source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED;
    321478            psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n");
    322             continue;
    323         }
    324 
    325         // we are basically classifying by moments
     479            Nskip ++;
     480            continue;
     481        }
     482
     483        // we are classifying by moments and PSF_MAG - KRON_MAG
    326484        psAssert (source->moments, "why is this source missing moments?");
    327485        if (source->mode & noMoments) {
    … …  
    333491        psF32 Mxx = source->moments->Mxx;
    334492        psF32 Myy = source->moments->Myy;
    335 
    336         // replace object in image
     493        psF32 Mxy = source->moments->Mxy;
     494        float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy));
     495
     496        // replace object in image to measure mag & psfWeights
    337497        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
    338498            pmSourceAdd (source, PM_MODEL_OP_FULL, options->maskVal);
    … …  
    343503        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal);
    344504
    345         // XXX can we test if psfMag is set and calculate only if needed?
    346505        pmSourceMagnitudes (source, psf, photMode, maskVal, markVal);
    347506
    … …  
    352511        pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal);
    353512
    354 # if (KRON)
    355         float apMag = -2.5*log10(source->moments->KronFlux);
    356 # else
    357         float apMag = -2.5*log10(source->moments->Sum);
    358 # endif
    359         float dMag = source->psfMag - apMag;
    360 
    361         // set nSigma to include both systematic and poisson error terms
    362         // XXX the 'poisson error' contribution for size is probably wrong...
    363         // XXX add in a hard floor on the Ap Sys Err (to be a bit generous)
    364         float nSigmaMAG = (dMag - options->ApResid) / hypot(source->errMag, hypot(options->ApSysErr, 0.025));
    365         float nSigmaMXX = (Mxx - psfClump->X) / hypot(psfClump->dX, psfClump->X*psfClump->X*source->errMag);
    366         float nSigmaMYY = (Myy - psfClump->Y) / hypot(psfClump->dY, psfClump->Y*psfClump->Y*source->errMag);
    367 
    368         // 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",
    369         // source->peak->xf, source->peak->yf, Mxx, Myy, source->peak->xf - source->moments->Mx, source->peak->yf - source->moments->My,
    370         // source->psfMag, apMag, dMag, source->errMag, nSigmaMAG, nSigmaMXX, nSigmaMYY);
    371 
    372         // XXX double check on ths stuff!! partially-masked sources are more likely to be mis-measured PSFs
    373         float sizeBias = 1.0;
    374         if (source->pixWeightNotBad < 0.9) {
    375             sizeBias = 3.0;
    376         }
    377         if (source->pixWeightNotPoor < 0.9) {
    378             sizeBias = 3.0;
    379         }
    380 
    381         float minMxx = psfClump->X - sizeBias*options->nSigmaMoments*psfClump->dX;
    382         float minMyy = psfClump->Y - sizeBias*options->nSigmaMoments*psfClump->dY;
    383 
    384         // include MAG, MXX, and MYY?
     513        float kMag = -2.5*log10(source->moments->KronFlux);
     514        float dMag = source->psfMag - kMag;
     515
     516        // sources without a valid magnitudes cannot have their size measured
     517        if (!isfinite(kMag) || !isfinite(source->psfMag)) {
     518            Nskip ++;
     519            continue;
     520        }
     521
     522        // set nSigmaMAG to include both systematic and poisson error terms.  we include a hard
     523        // floor on the Ap Sys Err (to be a bit generous).  XXX put the floor in the recipe...
     524        float nSigmaMAG = (dMag - options->ApResid) / hypot(source->errMag, hypot(options->ApSysErr, 0.02));
    385525        source->extNsigma = nSigmaMAG;
    386526
    … …  
    389529        // * CR & defect should have a faintess limit (min S/N)
    390530        // * SAT stars should not be faint, but defects may?
    391 
    392         // Anything within this region is a probably PSF-like object. Saturated stars may land
    393         // in this region, but are detected elsewhere on the basis of their peak value.
    394         bool isPSF = (fabs(nSigmaMAG) < options->nSigmaApResid) && (fabs(nSigmaMXX) < sizeBias*options->nSigmaMoments) && (fabs(nSigmaMYY) < sizeBias*options->nSigmaMoments);
    395         if (isPSF) {
    396           psTrace("psphotSourceClassRegion.PSF",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g PSF\t%g %g\n",
    397                   source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
    398                   options->nSigmaApResid,sizeBias*options->nSigmaMoments);
    399           source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    400           Npsf ++;
    401           continue;
    402         }
    403531
    404532        // Defects may not always match CRs from peak curvature analysis
    … …  
    407535        // XXX only accept brightish detections as CRs
    408536        // (nSigmaMAG < -options->nSigmaApResid) ||
    409         bool isCR = isCR = (source->errMag < 1.0 / SIZE_SN_LIM) && ((Mxx < minMxx) || (Myy < minMyy));
     537
     538        // saturated star (too many saturated pixels or peak above saturation limit).  These
     539        // may also be saturated galaxies, or just large saturated regions.
     540        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
     541            psTrace("psphotSourceClassRegion.SAT",4,"CLASS: %g %g\t%g %g\t%g %g\t%g %g\t%g SAT\n",
     542                    source->peak->xf, source->peak->yf, Mminor, kMag, dMag, nSigmaMAG, options->sizeLimitCR, options->magLimitCR, options->nSigmaApResid);
     543            source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
     544            Nsat ++;
     545            continue;
     546        }
     547
     548        // any sources missing a large fraction should just be treated as PSFs
     549        if ((source->pixWeightNotBad < 0.9) || (source->pixWeightNotPoor < 0.9)) {
     550            psTrace("psphotSourceClassRegion.PSF",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g PSF\t%g %g\n",
     551                    source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,kMag,dMag,nSigmaMAG,
     552                    options->nSigmaApResid,options->nSigmaMoments);
     553            source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
     554            Npsf ++;
     555            continue;
     556        }
     557
     558        // CRs are flagged by a combination on Mminor < options->sizeLimitCR && kmag < options->magLimitCR
     559        // NOTE: we only flag the CRs here; when we mask them we verify their CR nature (otherwise -> PSF)
     560        // bool isCR = (kMag < options->magLimitCR) && (Mminor < options->sizeLimitCR);
     561        // XXX skip if we have already marked it??
     562        bool isCR = (source->moments->SN > 7.0) && (Mminor < options->sizeLimitCR);
    410563        if (isCR) {
    411             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",
    412                     source->peak->xf,source->peak->yf,source->pixWeightNotBad,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
    413                     options->nSigmaApResid,sizeBias*options->nSigmaMoments);
     564            psTrace("psphotSourceClassRegion.CR",4,"CLASS: %g %g\t%g %g\t%g %g\t%g %g\t%g CR\n",
     565                    source->peak->xf, source->peak->yf, Mminor, kMag, dMag, nSigmaMAG, options->sizeLimitCR, options->magLimitCR, options->nSigmaApResid);
    414566            source->mode |= PM_SOURCE_MODE_DEFECT;
    415567            source->tmpFlags |= PM_SOURCE_TMPF_SIZE_CR_CANDIDATE;
     568            source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    416569            Ncr ++;
    417570            continue;
    418571        }
    419572
    420         // saturated star (determined in PSF fit).  These may also be saturated galaxies, or
    421         // just large saturated regions.
    422         if (source->mode & PM_SOURCE_MODE_SATSTAR) {
    423             psTrace("psphotSourceClassRegion.SAT",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g SAT\t%g %g\n",
    424                     source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
    425                     options->nSigmaApResid,sizeBias*options->nSigmaMoments);
    426             source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    427             Nsat ++;
    428             continue;
    429         }
    430 
    431         // XXX allow the Mxx, Myy to be less than psfClump->X,Y (by some nSigma)?
    432         bool isEXT = (nSigmaMAG > options->nSigmaApResid) || (Mxx > minMxx) || (Myy > minMyy);
     573        // Likely extended source (PSF_MAG - KRON_MAG is larger than limit)
     574        bool isEXT = (nSigmaMAG > options->nSigmaApResid);
    433575        if (isEXT) {
    434           psTrace("psphotSourceClassRegion.EXT",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Ext\t%g %g\n",
    435                   source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
    436                   options->nSigmaApResid,sizeBias*options->nSigmaMoments);
    437 
     576            psTrace("psphotSourceClassRegion.EXT",4,"CLASS: %g %g\t%g %g\t%g %g\t%g %g\t%g EXT\n",
     577                    source->peak->xf, source->peak->yf, Mminor, kMag, dMag, nSigmaMAG, options->sizeLimitCR, options->magLimitCR, options->nSigmaApResid);
    438578            source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
    439579            source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    … …  
    441581            continue;
    442582        }
    443         psTrace("psphotSourceClassRegion.MISS",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Unk\t%g %g\n",
    444                 source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
    445                 options->nSigmaApResid,sizeBias*options->nSigmaMoments);
    446 
    447         // sources that reach here are probably too faint for a reasonable source size measurement
    448         // psWarning ("sourse size was missed for %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigmaMAG);
    449         source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    450         Nmiss ++;
    451     }
    452 
    453     psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4d %4d %4d %4d %4d %4d", Npsf, Next, Nsat, Ncr, Nmiss, Nskip);
     583
     584        // Everything else should just be treated as a PSF
     585        psTrace("psphotSourceClassRegion.PSF",4,"CLASS: %g %g\t%g %g\t%g %g\t%g %g\t%g PSF\n",
     586                source->peak->xf, source->peak->yf, Mminor, kMag, dMag, nSigmaMAG, options->sizeLimitCR, options->magLimitCR, options->nSigmaApResid);
     587        source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
     588        Npsf ++;
     589    }
     590
     591    psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4d %4d %4d %4d %4d", Npsf, Next, Nsat, Ncr, Nskip);
    454592
    455593    return true;
    … …  
    458596// given an object suspected to be a defect, generate a pixel mask using the Lapacian transform
    459597// if enough of the object is detected as 'sharp', consider the object a cosmic ray
    460 bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) {
     598bool psphotSourceSelectCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) {
    461599
    462600    psTimerStart ("psphot.cr");
    … …  
    466604        pmSource *source = sources->data[i];
    467605
    468         // skip source if it was already measured
    469         if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) {
    470             psTrace("psphot", 7, "Not calculating source size since it has already been measured\n");
    471             continue;
    472         }
    473 
    474606        // only check candidates marked above
    475607        if (!(source->tmpFlags & PM_SOURCE_TMPF_SIZE_CR_CANDIDATE)) {
    … …  
    478610        }
    479611
    480         // skip unless this source is thought to be a cosmic ray.  flag the detection and mask the pixels
    481         // XXX this may be degenerate with the above test
    482         if (!(source->mode & PM_SOURCE_MODE_DEFECT)) continue;
     612        // once we are here, remove the temporary flag (this allows us to try again if we do not mark it now)
     613        source->tmpFlags &= ~PM_SOURCE_TMPF_SIZE_CR_CANDIDATE;
    483614
    484615        // Integer position of peak
    … …  
    504635        // XXX this is running slowly and is too agressive, but it more-or-less works
    505636        psTrace("psphot", 6, "mask cosmic ray at %f, %f\n", source->peak->xf, source->peak->yf);
    506         if (options->apply) {
     637        if (options->applyCRmask) {
    507638            psphotMaskCosmicRay(readout, source, options->crMask);
    508639        } else {
    … …  
    531662
    532663    // XXX test : save the mask image
    533     if (0) {
    534         psphotSaveImage (NULL, readout->mask,   "mask.fits");
    535     }
     664# if (DUMPPICS)
     665    psphotSaveImage (NULL, readout->mask,   "crmask.fits");
     666# endif
    536667
    537668    return true;
    538669}
    539670
    540 # define DUMPPICS 0
    541671# define LIMIT_XRANGE(X, IMAGE) { X = PS_MIN(PS_MAX(0, X), IMAGE->numCols); }
    542672# define LIMIT_YRANGE(Y, IMAGE) { Y = PS_MIN(PS_MAX(0, Y), IMAGE->numRows); }
    … …  
    721851    if (nCRpix > 1) {
    722852        source->mode |= PM_SOURCE_MODE_CR_LIMIT;
    723         source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    724853    }
    725854    // fprintf (stderr, "CRMASK %d %d %d %d %d\n", peak->x, peak->y, dx, dy, nCRpix);
    … …  
    12161345}
    12171346
     1347float psphotSourceSizeFindThreshold (psVector *value, psVector *mask, int maskValue, float minValue, float maxValue, float delta, float guess, float fraction) {
     1348
     1349    // make a histogram of the sources with mKron < magLimit
     1350    int nValue = (int)((maxValue - minValue) / delta) + 1;
     1351
     1352    psVector *histogram = psVectorAlloc(nValue, PS_TYPE_S32);
     1353    psVectorInit (histogram, 0);
     1354
     1355    for (int i = 0; i < value->n; i++) {
     1356        if (mask->data.U8[i] & maskValue) continue;
     1357        int bin = (value->data.F32[i] - minValue) / delta;
     1358        if (bin < 0.0) continue;
     1359        if (bin > histogram->n - 1) continue;
     1360        histogram->data.S32[bin] ++;
     1361    }
     1362
     1363    // find the peak of this histogram, but stop search at guess
     1364    int last = (guess - minValue) / delta;
     1365    int nPeak = 0;
     1366    int vPeak = histogram->data.S32[0];
     1367    for (int i = 1; (i < last) && (i < histogram->n); i++) {
     1368        if (histogram->data.S32[i] < vPeak) continue;
     1369        nPeak = i;
     1370        vPeak = histogram->data.S32[i];
     1371    }
     1372
     1373    // start at the peak and find the valley between here and the end of the histogram
     1374    int nValley = nPeak;
     1375    int vValley = histogram->data.S32[nPeak];
     1376    for (int i = nPeak + 1; i < histogram->n; i++) {
     1377        if (histogram->data.S32[i] > vValley) continue;
     1378        nValley = i;
     1379        vValley = histogram->data.S32[i];
     1380    }
     1381
     1382    psLogMsg ("psphot", PS_LOG_MINUTIA, "CR limits threshold : peak %d @ %f, valley %d @ %f (%f sigma)\n",
     1383              vPeak, delta * nPeak + minValue, vValley, delta * nValley + minValue, vPeak / PS_MAX(1.0, sqrt(vValley)));
     1384
     1385    if (nValley == nPeak) {
     1386        psFree (histogram);
     1387        return NAN;
     1388    }
     1389
     1390    if (vPeak < 3.0*sqrt(vValley)) {
     1391        psFree (histogram);
     1392        return NAN;
     1393    }
     1394
     1395    /// search for the first bin after the peak with f < vValley + 0.05*(vPeak - vValley)
     1396    float vLimit = vValley + fraction*(vPeak - vValley);
     1397    int nLimit = nValley;
     1398    for (int i = nPeak; i < nValley; i++) {
     1399        if (histogram->data.S32[i] > vLimit) continue;
     1400        nLimit = i;
     1401        break;
     1402    }
     1403    float result = nLimit * delta + minValue;
     1404
     1405    psFree (histogram);
     1406   
     1407    return result;
     1408}

  • trunk/psphot/src/psphotSourceStats.c

    r29004 r29548  
    4545    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
    4646    psAssert (readout, "missing readout?");
     47
     48    if (psTraceGetLevel("psphot") > 5) {
     49        static int pass = 0;
     50        char name[64];
     51        sprintf (name, "srstats.v%d.fits", pass);
     52        psphotSaveImage(NULL, readout->image, name);
     53        pass ++;
     54    }
    4755
    4856    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
    … …  
    418426           
    419427            float apMag = NAN;
    420             pmSourcePhotometryAper (&apMag, NULL, source->pixels, source->maskObj, maskVal);
     428            pmSourcePhotometryAper (&apMag, NULL, NULL, NULL, source->pixels, source->variance, source->maskObj, maskVal);
    421429            fprintf (stderr, "apMag: %f, kronMag: %f\n", apMag, -2.5*log10(source->moments->KronFlux));
    422430

  • trunk/psphot/src/psphotStackMatchPSFs.c

    r28013 r29548  
    8585
    8686    // set NAN pixels to 'SAT'
     87    // XXX replace this is pmReadoutMaskInvalid?
    8788    psImageMaskType maskVal = pmConfigMaskGet("SAT", config);
    8889    if (!pmReadoutMaskNonfinite(readoutSrc, maskVal)) {
    … …  
    106107    rescaleData(readoutOut, config, options, index);
    107108
    108     dumpImage(readoutOut, readoutSrc, index, "convolved");
     109    // dumpImage(readoutOut, readoutSrc, index, "convolved");
    109110
    110111    return true;

  • trunk/psphot/src/psphotStackMatchPSFsUtils.c

    r29027 r29548  
    361361        // Scale the input parameters
    362362        widthsCopy = psVectorCopy(NULL, widths, PS_TYPE_F32); // Copy of kernel widths
    363         if (scale && !pmSubtractionParamsScale(&size, &footprint, widthsCopy, options->inputSeeing->data.F32[index], options->targetSeeing, scaleRef, scaleMin, scaleMax)) {
     363
     364        // we need to register the FWHM values for use downstream
     365        pmSubtractionSetFWHMs(options->inputSeeing->data.F32[index], options->targetSeeing);
     366
     367        if (scale && !pmSubtractionParamsScale(&size, &footprint, widthsCopy, scaleRef, scaleMin, scaleMax)) {
    364368            psError(psErrorCodeLast(), false, "Unable to scale kernel parameters");
    365369            goto escape;

  • trunk/psphot/src/psphotStackParseCamera.c

    r28013 r29548  
    9191
    9292        if (!rawInputFile && !cnvInputFile) {
    93             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Component %s (%d) lacks IMAGE of type STR", item->name, i);
     93            psError(PSPHOT_ERR_CONFIG, true, "Component %s (%d) lacks both RAW:IMAGE and CNV:IMAGE of type STR", item->name, i);
    9494            return false;
    9595        }

  • trunk/psphot/src/psphotStackReadout.c

    r28013 r29548  
    154154
    155155    psphotStackObjectsUnifyPosition (objects);
     156
     157    // measure circular, radial apertures (objects sorted by S/N)
    156158    psphotRadialAperturesByObject (config, objects, view, STACK_OUT);
    157159
     160    // measure elliptical apertures, petrosians (objects sorted by S/N)
    158161    psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_OUT); // pass 1 (detections->allSources)
     162
     163    // measure non-linear extended source models (exponential, deVaucouleur, Sersic) (sources sorted by S/N)
    159164    psphotExtendedSourceFits (config, view, STACK_OUT); // pass 1 (detections->allSources)
    160165
    … …  
    162167    psphotMagnitudes(config, view, STACK_OUT);
    163168
    164     if (!psphotEfficiency(config, view, STACK_OUT)) {
     169    if (0 && !psphotEfficiency(config, view, STACK_OUT)) {
    165170        psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
    166171        psErrorClear();

  • trunk/psphot/src/psphotVisual.c

    r29004 r29548  
    2525static int kapa2 = -1;
    2626static int kapa3 = -1;
     27
     28/** destroy windows at the end of a run*/
     29bool psphotVisualClose(void)
     30{
     31    if(kapa1 != -1) KiiClose(kapa1);
     32    if(kapa2 != -1) KiiClose(kapa2);
     33    if(kapa3 != -1) KiiClose(kapa3);
     34    return true;
     35}
    2736
    2837int psphotKapaChannel (int channel) {
    … …  
    462471    if (myKapa == -1) return false;
    463472
    464     KapaClearPlots (myKapa);
     473    KapaClearSections (myKapa);
    465474    KapaInitGraph (&graphdata);
    466475    KapaSetFont (myKapa, "courier", 14);
    … …  
    13891398    section.bg  = KapaColorByName ("none"); // XXX probably should be 'none'
    13901399
    1391     KapaClearPlots (myKapa);
     1400    KapaClearSections (myKapa);
    13921401    // first section : mag vs CR nSigma
    13931402    section.dx = 1.0;
    … …  
    16211630    if (myKapa == -1) return false;
    16221631
    1623     KapaClearPlots (myKapa);
     1632    KapaClearSections (myKapa);
    16241633    KapaInitGraph (&graphdata);
    16251634    KapaSetFont (myKapa, "courier", 14);
    … …  
    21252134}
    21262135
     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
    21272403bool psphotVisualShowResidualImage (pmReadout *readout) {
    21282404
    … …  
    21382414}
    21392415
    2140 bool psphotVisualPlotApResid (psArray *sources, float mean, float error) {
     2416bool psphotVisualPlotApResid (psArray *sources, float mean, float error, bool useApMag) {
    21412417
    21422418    Graphdata graphdata;
    … …  
    21482424    if (myKapa == -1) return false;
    21492425
    2150     KapaClearPlots (myKapa);
     2426    KapaClearSections (myKapa);
    21512427    KapaInitGraph (&graphdata);
    21522428
    … …  
    21692445        if (!isfinite (source->psfMag)) continue;
    21702446
     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
    21712455        x->data.F32[n] = source->psfMag;
    2172         y->data.F32[n] = source->apMag - source->psfMag;
     2456        y->data.F32[n] = dMag;
    21732457        dy->data.F32[n] = source->errMag;
    21742458        graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]);
    … …  
    22662550    if (myKapa == -1) return false;
    22672551
    2268     KapaClearPlots (myKapa);
     2552    KapaClearSections (myKapa);
    22692553    KapaInitGraph (&graphdata);
    22702554
Note: See TracChangeset for help on using the changeset viewer.