Index: /branches/eam_branches/ipp-20100823/psphot/doc/notes.20101013.txt =================================================================== --- /branches/eam_branches/ipp-20100823/psphot/doc/notes.20101013.txt (revision 29490) +++ /branches/eam_branches/ipp-20100823/psphot/doc/notes.20101013.txt (revision 29491) @@ -1,2 +1,16 @@ + +20101019 + + We can identify CRs in the (M_minor, KronMag) plane as a fairly + tight clump (M_minor ~< 1.0, KronMag ~< -8.0). The best bounds of + the clump vary somewhat from image to image. I've added a dynamic + CR limit function which looks for the peak in M_minor and in KronMag + in those regions, then finds the valley floor. It seems to work + fairly well. Outstanding issues: + + * only adjust the limits once + * record limits in the headers + * record CR counts in the headers + * mask the CRs. 20101017 Index: /branches/eam_branches/ipp-20100823/psphot/src/psphot.h =================================================================== --- /branches/eam_branches/ipp-20100823/psphot/src/psphot.h (revision 29490) +++ /branches/eam_branches/ipp-20100823/psphot/src/psphot.h (revision 29491) @@ -258,6 +258,7 @@ bool psphotVisualShowFlags (psArray *sources); bool psphotVisualShowSourceSize (pmReadout *readout, psArray *sources); +bool psphotVisualPlotSourceSizeAlt (psMetadata *recipe, psMetadata *analysis, psArray *sources); bool psphotVisualShowResidualImage (pmReadout *readout); -bool psphotVisualPlotApResid (psArray *sources, float mean, float error); +bool psphotVisualPlotApResid (psArray *sources, float mean, float error, bool useApMag); bool psphotVisualPlotChisq (psArray *sources); bool psphotVisualPlotSourceSize (psMetadata *recipe, psMetadata *analysis, psArray *sources); Index: /branches/eam_branches/ipp-20100823/psphot/src/psphotApResid.c =================================================================== --- /branches/eam_branches/ipp-20100823/psphot/src/psphotApResid.c (revision 29490) +++ /branches/eam_branches/ipp-20100823/psphot/src/psphotApResid.c (revision 29491) @@ -349,5 +349,5 @@ psFree (dMag); - psphotVisualPlotApResid (sources, psf->ApResid, psf->dApResid); + psphotVisualPlotApResid (sources, psf->ApResid, psf->dApResid, true); return true; Index: /branches/eam_branches/ipp-20100823/psphot/src/psphotImageLoop.c =================================================================== --- /branches/eam_branches/ipp-20100823/psphot/src/psphotImageLoop.c (revision 29490) +++ /branches/eam_branches/ipp-20100823/psphot/src/psphotImageLoop.c (revision 29491) @@ -30,4 +30,10 @@ // files associated with the science image if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa in psphot."); + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + psImageMaskType maskTest = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // for psphot, we force data to be read at the chip level @@ -95,5 +101,5 @@ if (readout->mask) { psImageMaskType maskSat = pmConfigMaskGet("SAT", config); // Mask value for saturated pixels - if (!pmReadoutMaskNonfinite(readout, maskSat)) { + if (!pmReadoutMaskInvalid(readout, maskTest, maskSat)) { psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels."); psFree(view); Index: /branches/eam_branches/ipp-20100823/psphot/src/psphotOutput.c =================================================================== --- /branches/eam_branches/ipp-20100823/psphot/src/psphotOutput.c (revision 29490) +++ /branches/eam_branches/ipp-20100823/psphot/src/psphotOutput.c (revision 29491) @@ -270,4 +270,9 @@ psMetadataItemSupplement (&status, header, analysis, "NDET_CR"); + psMetadataItemSupplement (&status, header, analysis, "PSPHOT.PSF.APRESID"); + psMetadataItemSupplement (&status, header, analysis, "PSPHOT.PSF.APRESID.SYSERR"); + psMetadataItemSupplement (&status, header, analysis, "PSPHOT.CR.MAX.SIZE"); + psMetadataItemSupplement (&status, header, analysis, "PSPHOT.CR.MAX.MAG"); + // sky background model statistics psMetadataItemSupplement (&status, header, analysis, "MSKY_MN"); Index: /branches/eam_branches/ipp-20100823/psphot/src/psphotSourceSize.c =================================================================== --- /branches/eam_branches/ipp-20100823/psphot/src/psphotSourceSize.c (revision 29490) +++ /branches/eam_branches/ipp-20100823/psphot/src/psphotSourceSize.c (revision 29491) @@ -1,6 +1,4 @@ # include "psphotInternal.h" # include - -# define KRON 1 typedef struct { @@ -16,15 +14,20 @@ int grow; int xtest, ytest; - bool apply; // apply CR mask? + bool applyCRmask; // apply CR mask? + bool dynamicLimitsCR; // apply CR mask? + float sizeLimitCR; + float magLimitCR; } psphotSourceSizeOptions; // local functions: -bool psphotSourceSizePSF (psphotSourceSizeOptions *options, psArray *sources, pmPSF *psf); +bool psphotSourceSizePSF (psphotSourceSizeOptions *options, pmReadout *readout, psArray *sources, pmPSF *psf, psMetadata *recipe); +bool psphotDynamicLimitsCR (psphotSourceSizeOptions *options, pmReadout *readout, psArray *sources, pmPSF *psf, psMetadata *recipe); bool psphotSourceClass (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options); bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options); -bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options); +bool psphotSourceSelectCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options); bool psphotMaskCosmicRay (pmReadout *readout, pmSource *source, psImageMaskType maskVal); bool psphotMaskCosmicRayFootprintCheck (psArray *sources); int psphotMaskCosmicRayConnected (int xPeak, int yPeak, psImage *mymask, psImage *myvar, psImage *edges, int binning, float sigma_thresh); +float psphotSourceSizeFindThreshold (psVector *value, psVector *mask, int maskValue, float minValue, float maxValue, float delta, float guess, float fraction); // we need to call this function after sources have been fitted to the PSF model and @@ -111,5 +114,5 @@ assert (status); - // XXX recipe name is not great + // location of a single test source options.xtest = psMetadataLookupS32 (&status, recipe, "PSPHOT.CRMASK.XTEST"); options.ytest = psMetadataLookupS32 (&status, recipe, "PSPHOT.CRMASK.YTEST"); @@ -128,5 +131,5 @@ } - options.apply = psMetadataLookupBool(&status, recipe, "PSPHOT.CRMASK.APPLY"); // Growth size for CRs + options.applyCRmask = psMetadataLookupBool(&status, recipe, "PSPHOT.CRMASK.APPLY"); // Growth size for CRs if (!status) { psError(PS_ERR_BAD_PARAMETER_VALUE, true, "PSPHOT.CRMASK.APPLY is not defined."); @@ -134,21 +137,16 @@ } - // We are using the value psfMag - 2.5*log10(moment.Sum) as a measure of the extendedness - // of an object. We need to model this distribution for the PSF stars before we can test - // the significance for a specific object - // XXX move this to the code that generates the PSF? - // XXX store the results on pmPSF? - - // XXX this should only be done on the first pass (ie, if we have newSources or allSources?) - if (getPSFsize) { - psphotSourceSizePSF (&options, sources, psf); - } - - // classify the sources based on ApResid and Moments (extended sources) + // determine the distribution of (PSF_mag - KRON_mag) for the PSF sources (saved on readout->analysis) + psphotSourceSizePSF (&options, readout, sources, psf, recipe); + + // adjust the user-supplied limits based on the distribution of CRs in the (Mminor, mKron) space + psphotDynamicLimitsCR(&options, readout, sources, psf, recipe); + + // classify the sources based on ApResid and Moments // NOTE: only sources not already measured !(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) psphotSourceClass(readout, sources, recipe, psf, &options); - // NOTE: only sources not already measured !(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) - psphotSourceSizeCR (readout, sources, &options); + // attempt to mask the candidate CRs; flag if CR nature is confirmed + psphotSourceSelectCR(readout, sources, &options); // XXX fix this (was source->n - first) @@ -156,6 +154,7 @@ psphotVisualPlotSourceSize (recipe, readout->analysis, sources); + psphotVisualPlotSourceSizeAlt (recipe, readout->analysis, sources); psphotVisualShowSourceSize (readout, sources); - psphotVisualPlotApResid (sources, options.ApResid, options.ApSysErr); + psphotVisualPlotApResid (sources, options.ApResid, options.ApSysErr, false); psphotVisualShowSatStars (recipe, psf, sources); @@ -163,9 +162,33 @@ } +# define SAVE_PSF_OPTIONS(RO,OPT) \ + /* save these on readout->analysis (a) for output to the header and (b) to prevent re-calculating in another pass */ \ + psMetadataAddF32(RO->analysis, PS_LIST_TAIL, "PSPHOT.PSF.APRESID", PS_META_REPLACE, "locus of PSF stars in PSF_MAG - KRON_MAG", OPT->ApResid); \ + psMetadataAddF32(RO->analysis, PS_LIST_TAIL, "PSPHOT.PSF.APRESID.SYSERR", PS_META_REPLACE, "systematic error of PSF_MAG - KRON_MAG", OPT->ApSysErr); + // model the apmifit distribution for the psf stars: -bool psphotSourceSizePSF (psphotSourceSizeOptions *options, psArray *sources, pmPSF *psf) { +bool psphotSourceSizePSF (psphotSourceSizeOptions *options, pmReadout *readout, psArray *sources, pmPSF *psf, psMetadata *recipe) { + + // We are using the value PSF_MAG - KRON_MAG as a measure of the extendedness of an object. + // We need to model this distribution for the PSF stars before we can test the significance + // for a specific object. + + // NOTE: we require that objects have had moments measured, and we also require psfMags to + // be calculated. but, we do not require aperture mags or any other photometry values that + // require pixel analysis. + + bool status1 = false; + bool status2 = false; + options->ApResid = psMetadataLookupF32 (&status1, readout->analysis, "PSPHOT.PSF.APRESID"); + options->ApSysErr = psMetadataLookupF32 (&status2, readout->analysis, "PSPHOT.PSF.APRESID.SYSERR"); + psAssert ((status1 && status2) || (!status1 && !status2), "inconsistent record of PSF ApResid values"); + + // if they are saved on readout->analysis, we have already calculated these values + if (status1 && status2) { + return true; + } // select stats from the psf stars - psVector *Ap = psVectorAllocEmpty (100, PS_TYPE_F32); + psVector *ApOff = psVectorAllocEmpty (100, PS_TYPE_F32); psVector *ApErr = psVectorAllocEmpty (100, PS_TYPE_F32); @@ -173,6 +196,6 @@ psImageMaskType maskVal = options->maskVal | options->markVal; - // XXX why PHOT_WEIGHT?? - pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT; + // with PSFONLY, we do not need modify the pixels + pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_PSFONLY; int num = 0; // Number of sources measured @@ -182,37 +205,19 @@ num++; - // replace object in image - if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { - pmSourceAdd (source, PM_MODEL_OP_FULL, options->maskVal); - } - - // clear the mask bit and set the circular mask pixels - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal)); - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal); - - // XXX can we test if psfMag is set and calculate only if needed? pmSourceMagnitudes (source, psf, photMode, maskVal, markVal); - // clear the mask bit - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal)); - - // re-subtract the object, leave local sky - pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal); - - // XXX test: switch to kron flux -# if (KRON) - float apMag = -2.5*log10(source->moments->KronFlux); -# else - float apMag = -2.5*log10(source->moments->Sum); -# endif - float dMag = source->psfMag - apMag; - - psVectorAppend (Ap, dMag); + float kMag = -2.5*log10(source->moments->KronFlux); + float dMag = source->psfMag - kMag; + + psVectorAppend (ApOff, dMag); psVectorAppend (ApErr, source->errMag); } if (num == 0) { - // Not raising an error, because errors aren't being checked elsewhere in this function - psFree(Ap); + // if we cannot determine the PSF distribution, call all objects PSFs... + options->ApResid = NAN; + options->ApSysErr = NAN; + psFree(ApOff); psFree(ApErr); + SAVE_PSF_OPTIONS(readout, options); return false; } @@ -221,23 +226,161 @@ // psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); psStats *stats = psStatsAlloc(PS_STAT_CLIPPED_MEAN); - psVectorStats (stats, Ap, NULL, NULL, 0); - - psVector *dAp = psVectorAlloc (Ap->n, PS_TYPE_F32); - for (int i = 0; i < Ap->n; i++) { - dAp->data.F32[i] = Ap->data.F32[i] - stats->clippedMean; + psVectorStats (stats, ApOff, NULL, NULL, 0); + + psVector *dAp = psVectorAlloc (ApOff->n, PS_TYPE_F32); + for (int i = 0; i < ApOff->n; i++) { + dAp->data.F32[i] = ApOff->data.F32[i] - stats->clippedMean; } options->ApResid = stats->clippedMean; options->ApSysErr = psVectorSystematicError(dAp, ApErr, 0.1); - // XXX this is quite arbitrary... - if (!isfinite(options->ApSysErr)) options->ApSysErr = 0.01; + + // this is quite arbitrary... a large value means fewer things classified as extended. + if (!isfinite(options->ApSysErr)) options->ApSysErr = 0.05; psLogMsg ("psphot", PS_LOG_DETAIL, "psf - Sum: %f +/- %f\n", options->ApResid, options->ApSysErr); - psFree (Ap); + psFree (ApOff); psFree (ApErr); psFree (stats); psFree (dAp); + SAVE_PSF_OPTIONS(readout, options); return true; +} + +# define SAVE_CR_OPTIONS(RO,OPT) \ + /* save these on readout->analysis (a) for output to the header and (b) to prevent re-calculating in another pass */ \ + psMetadataAddF32(RO->analysis, PS_LIST_TAIL, "PSPHOT.CR.MAX.SIZE", PS_META_REPLACE, "dynamically-set minor axis size limit for cosmic rays", OPT->sizeLimitCR); \ + psMetadataAddF32(RO->analysis, PS_LIST_TAIL, "PSPHOT.CR.MAX.MAG", PS_META_REPLACE, "dynamically-set kron magnitude limit for cosmic rays", OPT->magLimitCR); + +// model the size and magnitude distribution of the Cosmic Rays +// ** CRs are reliably flagged by a combination on Mminor < X && mag (or flux) > Y +bool psphotDynamicLimitsCR (psphotSourceSizeOptions *options, pmReadout *readout, psArray *sources, pmPSF *psf, psMetadata *recipe) { + + /* attempt to describe the CR sources: + - input parameters are sizeLimit, magLimit + - first try to refine the sizeLimit: + -- select objects which meet the magLimit and exceed the sizeLimit by a factor of 1.5 + -- generate the histogram + -- look for a peak in the histogram + -- look for the min between valley and upper limit + -- look for first bin within 5% of the valley floor after peak (new sizeLimit) + + - next try to refine the magLimit + -- select objects which meet the sizeLimit and go fainter than the magLimit by 1.0 mag + -- generate the histogram + -- look for a peak in the histogram + -- look for the min between valley and upper limit + -- look for first bin within 5% of the valley floor after peak (new magLimit) + */ + + bool status = false; + bool status1 = false; + bool status2 = false; + options->sizeLimitCR = psMetadataLookupF32 (&status1, readout->analysis, "PSPHOT.CR.MAX.SIZE"); + if (!status1) { + options->sizeLimitCR = psMetadataLookupF32 (&status, recipe, "PSPHOT.CR.MAX.SIZE"); + if (!status) { + options->sizeLimitCR = 1.0; + } + } + options->magLimitCR = psMetadataLookupF32 (&status2, readout->analysis, "PSPHOT.CR.MAX.MAG"); + if (!status2) { + options->magLimitCR = psMetadataLookupF32 (&status, recipe, "PSPHOT.CR.MAX.MAG"); + if (!status) { + options->magLimitCR = -8.0; + } + } + psAssert ((status1 && status2) || (!status1 && !status2), "inconsistent record of dynamic CR limits"); + + // if they are saved on readout->analysis, we have already calculated these values + if (status1 && status2) { + return true; + } + + // if we do not want to dynamically set these, save the user-supplied values and exit + options->dynamicLimitsCR = psMetadataLookupBool (&status, recipe, "PSPHOT.CR.AUTOSCALE"); + if (!status) { + options->dynamicLimitsCR = true; + } + if (!options->dynamicLimitsCR) { + SAVE_CR_OPTIONS(readout, options); // macro defined above + return true; + } + + psVector *minor = psVectorAllocEmpty (100, PS_TYPE_F32); + psVector *mKron = psVectorAllocEmpty (100, PS_TYPE_F32); + psVector *mask = NULL; + + psImageMaskType markVal = options->markVal; + psImageMaskType maskVal = options->maskVal | options->markVal; + + // with PSFONLY, we do not need modify the pixels + pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_PSFONLY; + + // generate vectors for all the objects of possible interest (so we can just access those + // vectors in the next sections) + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + + // XXX can we test if psfMag is set and calculate only if needed? + pmSourceMagnitudes (source, psf, photMode, maskVal, markVal); + + // convert to Mmaj, Mmin: + psF32 Mxx = source->moments->Mxx; + psF32 Myy = source->moments->Myy; + psF32 Mxy = source->moments->Mxy; + + float KronMag = -2.5*log10(source->moments->KronFlux); + + float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy)); + + if (Mminor > options->sizeLimitCR * 1.5) continue; + if (KronMag > options->magLimitCR + 2.5) continue; + + psVectorAppend (mKron, KronMag); + psVectorAppend (minor, Mminor); + } + + // if too few objects meet the criterion, give up.. + if (mKron->n < 50) goto escape; + + // set distinct masks for (minor > sizeLimit) or (mKron > magLimit) + mask = psVectorAlloc(mKron->n, PS_TYPE_U8); + psVectorInit(mask, 0); + for (int i = 0; i < mKron->n; i++) { + if (mKron->data.F32[i] > options->magLimitCR) { + mask->data.U8[i] |= 0x01; + } + if (minor->data.F32[i] > options->sizeLimitCR) { + mask->data.U8[i] |= 0x02; + } + } + + float delta1 = PS_MAX(0.02, PS_MIN(0.2, 5.0 * 0.5 / minor->n)); + float newSizeLimit = psphotSourceSizeFindThreshold(minor, mask, 0x01, 0.0, options->sizeLimitCR * 1.5, delta1, options->sizeLimitCR, 0.05); + if (isfinite(newSizeLimit)) { + options->sizeLimitCR = newSizeLimit; + } + + float delta2 = PS_MAX(0.02, PS_MIN(0.2, 5.0 * 2.0 / minor->n)); + float newMagLimit = psphotSourceSizeFindThreshold(mKron, mask, 0x02, -15.0, options->magLimitCR + 2.5, delta2, options->magLimitCR, 0.05); + if (isfinite(newMagLimit)) { + options->magLimitCR = newMagLimit; + } + + psLogMsg ("psphot", PS_LOG_DETAIL, "CR limits : %f mag | %f pix^2\n", options->magLimitCR, options->sizeLimitCR); + + psFree (mKron); + psFree (minor); + psFree (mask); + + // save these on readout->analysis (a) for output to the header and (b) to prevent re-calculating in another pass + SAVE_CR_OPTIONS(readout, options); // macro defined above + return true; + + escape: + SAVE_CR_OPTIONS(readout, options); // macro defined above + return false; } @@ -249,5 +392,12 @@ char regionName[64]; - psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4s %4s %4s %4s %4s %4s", "Npsf", "Next", "Nsat", "Ncr", "Nmiss", "Nskip"); + psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4s %4s %4s %4s %4s", "Npsf", "Next", "Nsat", "Ncr", "Nskip"); + + if (!psphotSourceClassRegion (NULL, &psfClump, sources, recipe, psf, options)) { + psLogMsg ("psphot", 4, "Failed to determine source classification for full image\n"); + } else { + psLogMsg ("psphot", 4, "source classification for full image\n"); + } + return true; int nRegions = psMetadataLookupS32 (&status, readout->analysis, "PSF.CLUMP.NREGIONS"); @@ -275,4 +425,5 @@ continue; } + psLogMsg ("psphot", 4, "source classification for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1); // psphotVisualPlotSourceSize (recipe, readout->analysis, sources); } @@ -281,5 +432,4 @@ } -# define SIZE_SN_LIM 10 bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options) { @@ -291,9 +441,10 @@ int Npsf = 0; int Ncr = 0; - int Nmiss = 0; int Nskip = 0; pmSourceMode noMoments = PM_SOURCE_MODE_MOMENTS_FAILURE | PM_SOURCE_MODE_SKYVAR_FAILURE | PM_SOURCE_MODE_SKY_FAILURE | PM_SOURCE_MODE_BELOW_MOMENTS_SN; - pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT; + + // request the pixWeight values as well as the magnitudes + pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT; psImageMaskType markVal = options->markVal; @@ -306,8 +457,10 @@ // psfClumps are found for image subregions: // skip sources not in this region - if (source->peak->x < region->x0) continue; - if (source->peak->x >= region->x1) continue; - if (source->peak->y < region->y0) continue; - if (source->peak->y >= region->y1) continue; + if (region) { + if (source->peak->x < region->x0) continue; + if (source->peak->x >= region->x1) continue; + if (source->peak->y < region->y0) continue; + if (source->peak->y >= region->y1) continue; + } // skip source if it was already measured @@ -321,8 +474,9 @@ source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED; psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n"); - continue; - } - - // we are basically classifying by moments + Nskip ++; + continue; + } + + // we are classifying by moments and PSF_MAG - KRON_MAG psAssert (source->moments, "why is this source missing moments?"); if (source->mode & noMoments) { @@ -335,6 +489,7 @@ psF32 Myy = source->moments->Myy; psF32 Mxy = source->moments->Mxy; - - // replace object in image + float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy)); + + // replace object in image to measure mag & psfWeights if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { pmSourceAdd (source, PM_MODEL_OP_FULL, options->maskVal); @@ -345,5 +500,4 @@ psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal); - // XXX can we test if psfMag is set and calculate only if needed? pmSourceMagnitudes (source, psf, photMode, maskVal, markVal); @@ -354,38 +508,16 @@ pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal); -# if (KRON) - float apMag = -2.5*log10(source->moments->KronFlux); -# else - float apMag = -2.5*log10(source->moments->Sum); -# endif - float dMag = source->psfMag - apMag; - - // set nSigma to include both systematic and poisson error terms - // XXX need to handle NAN psfMag and apMag (ie, skip?) - // XXX the 'poisson error' contribution for size is probably wrong... - // XXX add in a hard floor on the Ap Sys Err (to be a bit generous) - float nSigmaMAG = (dMag - options->ApResid) / hypot(source->errMag, hypot(options->ApSysErr, 0.025)); - float nSigmaMXX = (Mxx - psfClump->X) / hypot(psfClump->dX, psfClump->X*psfClump->X*source->errMag); - float nSigmaMYY = (Myy - psfClump->Y) / hypot(psfClump->dY, psfClump->Y*psfClump->Y*source->errMag); - - // XXX should I change psfClump to use minor, major? - - // 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", - // source->peak->xf, source->peak->yf, Mxx, Myy, source->peak->xf - source->moments->Mx, source->peak->yf - source->moments->My, - // source->psfMag, apMag, dMag, source->errMag, nSigmaMAG, nSigmaMXX, nSigmaMYY); - - // XXX double check on ths stuff!! partially-masked sources are more likely to be mis-measured PSFs - float sizeBias = 1.0; - if (source->pixWeightNotBad < 0.9) { - sizeBias = 3.0; - } - if (source->pixWeightNotPoor < 0.9) { - sizeBias = 3.0; - } - - float minMxx = psfClump->X - sizeBias*options->nSigmaMoments*psfClump->dX; - float minMyy = psfClump->Y - sizeBias*options->nSigmaMoments*psfClump->dY; - - // include MAG, MXX, and MYY? + float kMag = -2.5*log10(source->moments->KronFlux); + float dMag = source->psfMag - kMag; + + // sources without a valid magnitudes cannot have their size measured + if (!isfinite(kMag) || !isfinite(source->psfMag)) { + Nskip ++; + continue; + } + + // set nSigmaMAG to include both systematic and poisson error terms. we include a hard + // floor on the Ap Sys Err (to be a bit generous). XXX put the floor in the recipe... + float nSigmaMAG = (dMag - options->ApResid) / hypot(source->errMag, hypot(options->ApSysErr, 0.02)); source->extNsigma = nSigmaMAG; @@ -394,16 +526,4 @@ // * CR & defect should have a faintess limit (min S/N) // * SAT stars should not be faint, but defects may? - - // Anything within this region is a probably PSF-like object. Saturated stars may land - // in this region, but are detected elsewhere on the basis of their peak value. - bool isPSF = (fabs(nSigmaMAG) < options->nSigmaApResid) && (fabs(nSigmaMXX) < sizeBias*options->nSigmaMoments) && (fabs(nSigmaMYY) < sizeBias*options->nSigmaMoments); - if (isPSF) { - psTrace("psphotSourceClassRegion.PSF",4,"CLASS: %g %g\t%g %g %g %g %g %g\t%g %g\t%g PSF\t%g %g\n", - source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG, - options->nSigmaApResid,sizeBias*options->nSigmaMoments); - source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; - Npsf ++; - continue; - } // Defects may not always match CRs from peak curvature analysis @@ -413,42 +533,43 @@ // (nSigmaMAG < -options->nSigmaApResid) || - // ** CRs are reliably flagged by a combination on Mminor < X && mag (or flux) > Y - - float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy)); - if ((Mminor < 1.0) && (apMag < -8.0)) { - fprintf (stderr, "likely CR @ %f,%f\n", source->peak->xf, source->peak->yf); - } - - // XXX do I need to find the Mminor, Mmajor distribution? - - bool isCR = isCR = (source->errMag < 1.0 / SIZE_SN_LIM) && ((Mxx < minMxx) || (Myy < minMyy)); + // saturated star (too many saturated pixels or peak above saturation limit). These + // may also be saturated galaxies, or just large saturated regions. + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + psTrace("psphotSourceClassRegion.SAT",4,"CLASS: %g %g\t%g %g\t%g %g\t%g %g\t%g SAT\n", + source->peak->xf, source->peak->yf, Mminor, kMag, dMag, nSigmaMAG, options->sizeLimitCR, options->magLimitCR, options->nSigmaApResid); + source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; + Nsat ++; + continue; + } + + // any sources missing a large fraction should just be treated as PSFs + if ((source->pixWeightNotBad < 0.9) || (source->pixWeightNotPoor < 0.9)) { + psTrace("psphotSourceClassRegion.PSF",4,"CLASS: %g %g\t%g %g %g %g %g %g\t%g %g\t%g PSF\t%g %g\n", + source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,kMag,dMag,nSigmaMAG, + options->nSigmaApResid,options->nSigmaMoments); + source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; + Npsf ++; + continue; + } + + // CRs are flagged by a combination on Mminor < options->sizeLimitCR && kmag < options->magLimitCR + // NOTE: we only flag the CRs here; when we mask them we verify their CR nature (otherwise -> PSF) + // bool isCR = (kMag < options->magLimitCR) && (Mminor < options->sizeLimitCR); + bool isCR = (source->moments->SN > 7.0) && (Mminor < options->sizeLimitCR); if (isCR) { - 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", - source->peak->xf,source->peak->yf,source->pixWeightNotBad,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG, - options->nSigmaApResid,sizeBias*options->nSigmaMoments); + psTrace("psphotSourceClassRegion.CR",4,"CLASS: %g %g\t%g %g\t%g %g\t%g %g\t%g CR\n", + source->peak->xf, source->peak->yf, Mminor, kMag, dMag, nSigmaMAG, options->sizeLimitCR, options->magLimitCR, options->nSigmaApResid); source->mode |= PM_SOURCE_MODE_DEFECT; source->tmpFlags |= PM_SOURCE_TMPF_SIZE_CR_CANDIDATE; + source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; Ncr ++; continue; } - // saturated star (determined in PSF fit). These may also be saturated galaxies, or - // just large saturated regions. - if (source->mode & PM_SOURCE_MODE_SATSTAR) { - psTrace("psphotSourceClassRegion.SAT",4,"CLASS: %g %g\t%g %g %g %g %g %g\t%g %g\t%g SAT\t%g %g\n", - source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG, - options->nSigmaApResid,sizeBias*options->nSigmaMoments); - source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; - Nsat ++; - continue; - } - - // XXX allow the Mxx, Myy to be less than psfClump->X,Y (by some nSigma)? - bool isEXT = (nSigmaMAG > options->nSigmaApResid) || (Mxx > minMxx) || (Myy > minMyy); + // Likely extended source (PSF_MAG - KRON_MAG is larger than limit) + bool isEXT = (nSigmaMAG > options->nSigmaApResid); if (isEXT) { - psTrace("psphotSourceClassRegion.EXT",4,"CLASS: %g %g\t%g %g %g %g %g %g\t%g %g\t%g Ext\t%g %g\n", - source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG, - options->nSigmaApResid,sizeBias*options->nSigmaMoments); - + psTrace("psphotSourceClassRegion.EXT",4,"CLASS: %g %g\t%g %g\t%g %g\t%g %g\t%g EXT\n", + source->peak->xf, source->peak->yf, Mminor, kMag, dMag, nSigmaMAG, options->sizeLimitCR, options->magLimitCR, options->nSigmaApResid); source->mode |= PM_SOURCE_MODE_EXT_LIMIT; source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; @@ -456,15 +577,13 @@ continue; } - psTrace("psphotSourceClassRegion.MISS",4,"CLASS: %g %g\t%g %g %g %g %g %g\t%g %g\t%g Unk\t%g %g\n", - source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG, - options->nSigmaApResid,sizeBias*options->nSigmaMoments); - - // sources that reach here are probably too faint for a reasonable source size measurement - // psWarning ("sourse size was missed for %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigmaMAG); - source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; - Nmiss ++; - } - - psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4d %4d %4d %4d %4d %4d", Npsf, Next, Nsat, Ncr, Nmiss, Nskip); + + // Everything else should just be treated as a PSF + psTrace("psphotSourceClassRegion.PSF",4,"CLASS: %g %g\t%g %g\t%g %g\t%g %g\t%g PSF\n", + source->peak->xf, source->peak->yf, Mminor, kMag, dMag, nSigmaMAG, options->sizeLimitCR, options->magLimitCR, options->nSigmaApResid); + source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; + Npsf ++; + } + + psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4d %4d %4d %4d %4d", Npsf, Next, Nsat, Ncr, Nskip); return true; @@ -473,5 +592,5 @@ // given an object suspected to be a defect, generate a pixel mask using the Lapacian transform // if enough of the object is detected as 'sharp', consider the object a cosmic ray -bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) { +bool psphotSourceSelectCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) { psTimerStart ("psphot.cr"); @@ -481,10 +600,4 @@ pmSource *source = sources->data[i]; - // skip source if it was already measured - if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) { - psTrace("psphot", 7, "Not calculating source size since it has already been measured\n"); - continue; - } - // only check candidates marked above if (!(source->tmpFlags & PM_SOURCE_TMPF_SIZE_CR_CANDIDATE)) { @@ -493,7 +606,6 @@ } - // skip unless this source is thought to be a cosmic ray. flag the detection and mask the pixels - // XXX this may be degenerate with the above test - if (!(source->mode & PM_SOURCE_MODE_DEFECT)) continue; + // once we are here, remove the temporary flag + source->mode &= ~PM_SOURCE_TMPF_SIZE_CR_CANDIDATE; // Integer position of peak @@ -519,5 +631,5 @@ // XXX this is running slowly and is too agressive, but it more-or-less works psTrace("psphot", 6, "mask cosmic ray at %f, %f\n", source->peak->xf, source->peak->yf); - if (options->apply) { + if (options->applyCRmask) { psphotMaskCosmicRay(readout, source, options->crMask); } else { @@ -546,6 +658,6 @@ // XXX test : save the mask image - if (0) { - psphotSaveImage (NULL, readout->mask, "mask.fits"); + if (1) { + psphotSaveImage (NULL, readout->mask, "crmask.fits"); } @@ -736,5 +848,4 @@ if (nCRpix > 1) { source->mode |= PM_SOURCE_MODE_CR_LIMIT; - source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; } // fprintf (stderr, "CRMASK %d %d %d %d %d\n", peak->x, peak->y, dx, dy, nCRpix); @@ -1231,2 +1342,64 @@ } +float psphotSourceSizeFindThreshold (psVector *value, psVector *mask, int maskValue, float minValue, float maxValue, float delta, float guess, float fraction) { + + // make a histogram of the sources with mKron < magLimit + int nValue = (int)((maxValue - minValue) / delta) + 1; + + psVector *histogram = psVectorAlloc(nValue, PS_TYPE_S32); + psVectorInit (histogram, 0); + + for (int i = 0; i < value->n; i++) { + if (mask->data.U8[i] & maskValue) continue; + int bin = (value->data.F32[i] - minValue) / delta; + if (bin < 0.0) continue; + if (bin > histogram->n - 1) continue; + histogram->data.S32[bin] ++; + } + + // find the peak of this histogram, but stop search at guess + int last = (guess - minValue) / delta; + int nPeak = 0; + int vPeak = histogram->data.S32[0]; + for (int i = 1; (i < last) && (i < histogram->n); i++) { + if (histogram->data.S32[i] < vPeak) continue; + nPeak = i; + vPeak = histogram->data.S32[i]; + } + + // start at the peak and find the valley between here and the end of the histogram + int nValley = nPeak; + int vValley = histogram->data.S32[nPeak]; + for (int i = nPeak + 1; i < histogram->n; i++) { + if (histogram->data.S32[i] > vValley) continue; + nValley = i; + vValley = histogram->data.S32[i]; + } + + psLogMsg ("psphot", PS_LOG_MINUTIA, "CR limits threshold : peak %d @ %f, valley %d @ %f (%f sigma)\n", + vPeak, delta * nPeak + minValue, vValley, delta * nValley + minValue, vPeak / PS_MAX(1.0, sqrt(vValley))); + + if (nValley == nPeak) { + psFree (histogram); + return NAN; + } + + if (vPeak < 3.0*sqrt(vValley)) { + psFree (histogram); + return NAN; + } + + /// search for the first bin after the peak with f < vValley + 0.05*(vPeak - vValley) + float vLimit = vValley + fraction*(vPeak - vValley); + int nLimit = nValley; + for (int i = nPeak; i < nValley; i++) { + if (histogram->data.S32[i] > vLimit) continue; + nLimit = i; + break; + } + float result = nLimit * delta + minValue; + + psFree (histogram); + + return result; +} Index: /branches/eam_branches/ipp-20100823/psphot/src/psphotVisual.c =================================================================== --- /branches/eam_branches/ipp-20100823/psphot/src/psphotVisual.c (revision 29490) +++ /branches/eam_branches/ipp-20100823/psphot/src/psphotVisual.c (revision 29491) @@ -471,5 +471,5 @@ if (myKapa == -1) return false; - KapaClearPlots (myKapa); + KapaClearSections (myKapa); KapaInitGraph (&graphdata); KapaSetFont (myKapa, "courier", 14); @@ -1398,5 +1398,5 @@ section.bg = KapaColorByName ("none"); // XXX probably should be 'none' - KapaClearPlots (myKapa); + KapaClearSections (myKapa); // first section : mag vs CR nSigma section.dx = 1.0; @@ -1630,5 +1630,5 @@ if (myKapa == -1) return false; - KapaClearPlots (myKapa); + KapaClearSections (myKapa); KapaInitGraph (&graphdata); KapaSetFont (myKapa, "courier", 14); @@ -2134,4 +2134,271 @@ } +bool PlotSourceSizeAltSetVectors(psVector *m, psVector *k, psVector *v1, psVector *v2, psVector *v3, pmSource *source) { + + float Mxx = source->moments->Mxx; + float Myy = source->moments->Myy; + float Mxy = source->moments->Mxy; + float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy)); + float KronMag = -2.5*log10(source->moments->KronFlux); + + psVectorAppend(m, source->psfMag); + psVectorAppend(k, KronMag); + psVectorAppend(v1, Mminor); + psVectorAppend(v2, source->psfMag - KronMag); + psVectorAppend(v3, source->extNsigma); + return true; +} + +bool PlotSourceSizeAltAddPoints(Graphdata *graphdata, int myKapa, psVector *x, psVector *y, char *colorname, int ptype, float size) { + + graphdata->color = KapaColorByName (colorname); + graphdata->ptype = ptype; + graphdata->size = size; + graphdata->style = 2; + KapaPrepPlot (myKapa, x->n, graphdata); + KapaPlotVector (myKapa, x->n, x->data.F32, "x"); + KapaPlotVector (myKapa, x->n, y->data.F32, "y"); + return true; +} + +bool psphotVisualPlotSourceSizeAlt (psMetadata *recipe, psMetadata *analysis, psArray *sources) { + + Graphdata graphdata; + KapaSection section; + + if (!pmVisualTestLevel("psphot.size", 2)) return true; + + int myKapa = psphotKapaChannel (2); + if (myKapa == -1) return false; + + KapaClearSections (myKapa); + KapaInitGraph (&graphdata); + KapaSetFont (myKapa, "courier", 14); + + section.bg = KapaColorByName ("none"); // XXX probably should be 'none' + + // storage vectors for data to be plotted + psVector *SATm = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *SATk = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *SAT1 = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *SAT2 = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *SAT3 = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + + psVector *PSFm = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *PSFk = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *PSF1 = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *PSF2 = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *PSF3 = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + + psVector *EXTm = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *EXTk = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *EXT1 = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *EXT2 = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *EXT3 = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + + psVector *DEFm = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *DEFk = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *DEF1 = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *DEF2 = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *DEF3 = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + + psVector *BADm = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *BADk = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *BAD1 = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *BAD2 = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *BAD3 = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + + psVector *CRm = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *CRk = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *CR1 = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *CR2 = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *CR3 = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + + // int notPSF = PM_SOURCE_MODE_SATSTAR | PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT | PM_SOURCE_MODE_DEFECT; + int nSkip = 0; + + // construct the vectors + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + if (source->moments == NULL) { + nSkip ++; + continue; + } + + bool found = false; + + // only plot the measured sources... + if (!(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED)) { + nSkip ++; + continue; + } + + // any sources missing a large fraction should just be treated as PSFs + if ((source->pixWeightNotBad < 0.9) || (source->pixWeightNotPoor < 0.9)) { + PlotSourceSizeAltSetVectors(BADm, BADk, BAD1, BAD2, BAD3, source); + } + if ((source->mode & PM_SOURCE_MODE_CR_LIMIT) || (source->tmpFlags & PM_SOURCE_TMPF_SIZE_CR_CANDIDATE)) { + PlotSourceSizeAltSetVectors(CRm, CRk, CR1, CR2, CR3, source); + found = true; + } + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + PlotSourceSizeAltSetVectors(SATm, SATk, SAT1, SAT2, SAT3, source); + found = true; + } + if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { + PlotSourceSizeAltSetVectors(EXTm, EXTk, EXT1, EXT2, EXT3, source); + found = true; + } + if (source->mode & PM_SOURCE_MODE_DEFECT) { + PlotSourceSizeAltSetVectors(DEFm, DEFk, DEF1, DEF2, DEF3, source); + found = true; + } + if (!found) { + PlotSourceSizeAltSetVectors(PSFm, PSFk, PSF1, PSF2, PSF3, source); + } + // if (!(source->mode & notPSF) && !(source->tmpFlags & PM_SOURCE_TMPF_SIZE_CR_CANDIDATE)) { + // PlotSourceSizeAltSetVectors(PSFm, PSFk, PSF1, PSF2, PSF3, source); + // } + } + // three sections: kronMag vs Mminor, psfMag vs psfMag - KronMag, psfMag vs extNsigma + + // --- first section: kronMag vs Mminor --- + section.dx = 1.00; + section.dy = 0.33; + section.x = 0.00; + section.y = 0.00; + section.name = psStringCopy ("Mminor"); + KapaSetSection (myKapa, §ion); + psFree (section.name); + + graphdata.color = KapaColorByName ("black"); + graphdata.xmin = -17.1; + graphdata.xmax = -6.9; + graphdata.ymin = -0.5; + graphdata.ymax = +7.1; + KapaSetLimits (myKapa, &graphdata); + + graphdata.padXm = NAN; + graphdata.padYm = 5.0; + graphdata.padXp = NAN; + graphdata.padYp = NAN; + KapaBox (myKapa, &graphdata); + KapaSendLabel (myKapa, "kron mag", KAPA_LABEL_XM); + KapaSendLabel (myKapa, "M_minor| (pixels^2|)", KAPA_LABEL_YM); + + PlotSourceSizeAltAddPoints(&graphdata, myKapa, PSFk, PSF1, "black", 0, 0.5); + PlotSourceSizeAltAddPoints(&graphdata, myKapa, EXTk, EXT1, "blue", 0, 0.5); + PlotSourceSizeAltAddPoints(&graphdata, myKapa, CRk, CR1, "red", 0, 0.5); + PlotSourceSizeAltAddPoints(&graphdata, myKapa, DEFk, DEF1, "red", 7, 1.0); + PlotSourceSizeAltAddPoints(&graphdata, myKapa, SATk, SAT1, "blue", 7, 1.2); + PlotSourceSizeAltAddPoints(&graphdata, myKapa, BADk, BAD1, "green", 7, 1.4); + + // --- second section: dMag ---- + section.dx = 1.00; + section.dy = 0.33; + section.x = 0.00; + section.y = 0.33; + section.name = psStringCopy ("dMag"); + KapaSetSection (myKapa, §ion); + psFree (section.name); + + graphdata.color = KapaColorByName ("black"); + graphdata.xmin = -17.1; + graphdata.xmax = -6.9; + graphdata.ymin = -0.75; + graphdata.ymax = +1.50; + KapaSetLimits (myKapa, &graphdata); + + graphdata.padXm = NAN; + graphdata.padYm = 5.0; + graphdata.padXp = 0.0; + graphdata.padYp = NAN; + strcpy (graphdata.labels, "0200"); + KapaBox (myKapa, &graphdata); + KapaSendLabel (myKapa, "dMag", KAPA_LABEL_YM); + + PlotSourceSizeAltAddPoints(&graphdata, myKapa, PSFm, PSF2, "black", 0, 0.5); + PlotSourceSizeAltAddPoints(&graphdata, myKapa, EXTm, EXT2, "blue", 0, 0.5); + PlotSourceSizeAltAddPoints(&graphdata, myKapa, CRm, CR2, "red", 0, 0.5); + PlotSourceSizeAltAddPoints(&graphdata, myKapa, DEFm, DEF2, "red", 7, 1.0); + PlotSourceSizeAltAddPoints(&graphdata, myKapa, SATm, SAT2, "blue", 7, 1.2); + PlotSourceSizeAltAddPoints(&graphdata, myKapa, BADm, BAD2, "green", 7, 1.4); + + // --- third section: nSigma --- + section.dx = 1.00; + section.dy = 0.33; + section.x = 0.00; + section.y = 0.66; + section.name = psStringCopy ("nSigma"); + KapaSetSection (myKapa, §ion); + psFree (section.name); + + graphdata.color = KapaColorByName ("black"); + graphdata.xmin = -17.1; + graphdata.xmax = -6.9; + graphdata.ymin = -10.1; + graphdata.ymax = +10.1; + KapaSetLimits (myKapa, &graphdata); + + graphdata.padXm = 0.0; + graphdata.padYm = 5.0; + graphdata.padXp = NAN; + graphdata.padYp = NAN; + strcpy (graphdata.labels, "0210"); + KapaBox (myKapa, &graphdata); + KapaSendLabel (myKapa, "psf msg", KAPA_LABEL_XP); + KapaSendLabel (myKapa, "EXT nSigma", KAPA_LABEL_YM); + + PlotSourceSizeAltAddPoints(&graphdata, myKapa, PSFm, PSF3, "black", 0, 0.5); + PlotSourceSizeAltAddPoints(&graphdata, myKapa, EXTm, EXT3, "blue", 0, 0.5); + PlotSourceSizeAltAddPoints(&graphdata, myKapa, CRm, CR3, "red", 0, 0.5); + PlotSourceSizeAltAddPoints(&graphdata, myKapa, DEFm, DEF3, "red", 7, 1.0); + PlotSourceSizeAltAddPoints(&graphdata, myKapa, SATm, SAT3, "blue", 7, 1.0); + PlotSourceSizeAltAddPoints(&graphdata, myKapa, BADm, BAD3, "green", 7, 1.4); + + fprintf (stderr, "PSF: %ld, EXT: %ld, CR: %ld, DEF: %ld, SAT: %ld, TOTAL: %ld, nSkip: %d\n", + PSFm->n, EXTm->n, CRm->n, DEFm->n, SATm->n, + PSFm->n+ EXTm->n+ CRm->n+ DEFm->n+ SATm->n, nSkip); + + psFree (SATk); + psFree (SATm); + psFree (SAT1); + psFree (SAT2); + psFree (SAT3); + + psFree (EXTk); + psFree (EXTm); + psFree (EXT1); + psFree (EXT2); + psFree (EXT3); + + psFree (PSFk); + psFree (PSFm); + psFree (PSF1); + psFree (PSF2); + psFree (PSF3); + + psFree (DEFk); + psFree (DEFm); + psFree (DEF1); + psFree (DEF2); + psFree (DEF3); + + psFree (BADk); + psFree (BADm); + psFree (BAD1); + psFree (BAD2); + psFree (BAD3); + + psFree (CRk); + psFree (CRm); + psFree (CR1); + psFree (CR2); + psFree (CR3); + + pmVisualAskUser(NULL); + return true; +} + bool psphotVisualShowResidualImage (pmReadout *readout) { @@ -2147,5 +2414,5 @@ } -bool psphotVisualPlotApResid (psArray *sources, float mean, float error) { +bool psphotVisualPlotApResid (psArray *sources, float mean, float error, bool useApMag) { Graphdata graphdata; @@ -2157,5 +2424,5 @@ if (myKapa == -1) return false; - KapaClearPlots (myKapa); + KapaClearSections (myKapa); KapaInitGraph (&graphdata); @@ -2178,6 +2445,14 @@ if (!isfinite (source->psfMag)) continue; + float dMag; + if (useApMag) { + dMag = source->apMag - source->psfMag; + } else { + float kMag = -2.5*log10(source->moments->KronFlux); + dMag = source->psfMag - kMag; + } + x->data.F32[n] = source->psfMag; - y->data.F32[n] = source->apMag - source->psfMag; + y->data.F32[n] = dMag; dy->data.F32[n] = source->errMag; graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]); @@ -2275,5 +2550,5 @@ if (myKapa == -1) return false; - KapaClearPlots (myKapa); + KapaClearSections (myKapa); KapaInitGraph (&graphdata);