Changeset 29491 for branches/eam_branches/ipp-20100823/psphot/src/psphotSourceSize.c

Timestamp:

Oct 20, 2010, 8:55:04 AM (16 years ago)

Author:

eugene

Message:

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

File:

• branches/eam_branches/ipp-20100823/psphot/src/psphotSourceSize.c (modified) (29 diffs)

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

@@ -1 +1 @@
 # include "psphotInternal.h"
 # include <gsl/gsl_sf_gamma.h>
-
-# define KRON 1
 
 typedef struct {
@@ -16 +14 @@
     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 +114 @@
     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 +131 @@
     }
 
-    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 +137 @@
     }
 
-    // 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 +154 @@
 
     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 +162 @@
 }
 
+# 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 +196 @@
     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 +205 @@
         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 +226 @@
     // 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 +392 @@
     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 +425 @@
             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 +432 @@
 }
 
-# define SIZE_SN_LIM 10
 bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options) {
 
@@ -291 +441 @@
     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 +457 @@
         // 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 +474 @@
             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 +489 @@
         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 +500 @@
         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 +508 @@
         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 +526 @@
         // * 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 +533 @@
         // (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 +577 @@
             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 +592 @@
 // 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 +600 @@
         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 +606 @@
         }
 
-        // 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 +631 @@
         // 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 +658 @@
 
     // XXX test : save the mask image
-    if (0) {
-        psphotSaveImage (NULL, readout->mask,   "mask.fits");
+    if (1) {
+        psphotSaveImage (NULL, readout->mask,   "crmask.fits");
     }
 
@@ -736 +848 @@
     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 +1342 @@
 }
 
+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;
+}