Changeset 25755 for trunk/psphot/src/psphotSourceSize.c

Timestamp:

Oct 2, 2009, 3:12:47 PM (17 years ago)

Author:

eugene

Message:

check in changes from gene@development branch : extensive changes to moments calculation, psf model generation, aperture residuals; improvements to extended source analysis

Location:

trunk/psphot/src

Files:

• . (modified) (1 prop)
• psphotSourceSize.c (modified) (6 diffs)

trunk/psphot/src

@@ -18 +18 @@
 psphotVersionDefinitions.h
 psphotMomentsStudy
+psphotPetrosianStudy

@@ -18 +18 @@
 psphotVersionDefinitions.h
 psphotMomentsStudy
+psphotPetrosianStudy

trunk/psphot/src/psphotSourceSize.c

@@ -2 +2 @@
 # include <gsl/gsl_sf_gamma.h>
 
-static float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask,
-                                psImageMaskType maskVal, const pmModel *model, float Ro);
-
-bool psphotMaskCosmicRay_Old (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask);
-bool psphotMaskCosmicRay_New (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask);
+typedef struct {
+    psImageMaskType maskVal;
+    psImageMaskType markVal;
+    psImageMaskType crMask;
+    float ApResid;
+    float ApSysErr;
+    float nSigmaApResid;
+    float nSigmaMoments;
+    float nSigmaCR;
+    float soft;
+    int grow;
+} psphotSourceSizeOptions;
+
+bool psphotSourceSizePSF (psphotSourceSizeOptions *options, psArray *sources, pmPSF *psf);
+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 psphotMaskCosmicRay (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask);
+bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask);
 
 // we need to call this function after sources have been fitted to the PSF model and
@@ -14 +28 @@
 // deviation from the psf model at the r = FWHM/2 position
 
-bool psphotSourceSize(pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, long first)
+// XXX use an internal flag to mark sources which have already been measured
+bool psphotSourceSize(pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, long first)
 {
     bool status;
+    psphotSourceSizeOptions options;
 
     psTimerStart ("psphot.size");
 
     // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
-    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
-    assert (maskVal);
+    options.maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
+    assert (options.maskVal);
+
+    options.markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); // Mask value for bad pixels
+    assert (options.markVal);
 
     // bit to mask the cosmic-ray pixels
-    psImageMaskType crMask  = pmConfigMaskGet("CR", config); // Mask value for cosmic rays
-
-    float CR_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.CR.NSIGMA.LIMIT");
+    options.crMask  = pmConfigMaskGet("CR", config); // Mask value for cosmic rays
+
+    options.nSigmaCR = psMetadataLookupF32 (&status, recipe, "PSPHOT.CR.NSIGMA.LIMIT");
     assert (status);
 
-    float EXT_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.EXT.NSIGMA.LIMIT");
+    // XXX recipe name is not great
+    options.nSigmaApResid = psMetadataLookupF32 (&status, recipe, "PSPHOT.EXT.NSIGMA.LIMIT");
     assert (status);
 
-    int grow = psMetadataLookupS32(&status, recipe, "PSPHOT.CR.GROW"); // Growth size for CRs
-    if (!status || grow < 0) {
+    // XXX recipe name is not great
+    options.nSigmaMoments = psMetadataLookupF32 (&status, recipe, "PSPHOT.EXT.NSIGMA.MOMENTS");
+    assert (status);
+
+    options.grow = psMetadataLookupS32(&status, recipe, "PSPHOT.CR.GROW"); // Growth size for CRs
+    if (!status || options.grow < 0) {
         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "PSPHOT.CR.GROW is not positive.");
         return false;
     }
 
-    float soft = psMetadataLookupF32(&status, recipe, "PSPHOT.CR.NSIGMA.SOFTEN"); // Softening parameter
-    if (!status || !isfinite(soft) || soft < 0.0) {
+    options.soft = psMetadataLookupF32(&status, recipe, "PSPHOT.CR.NSIGMA.SOFTEN"); // Softening parameter
+    if (!status || !isfinite(options.soft) || options.soft < 0.0) {
         psWarning("PSPHOT.CR.NSIGMA.SOFTEN not set; defaulting to zero.");
-        soft = 0.0;
-    }
-
-    // loop over all source
-    for (int i = first; i < sources->n; i++) {
-        pmSource *source = sources->data[i];
-
-        // skip source if it was already measured
-        if (isfinite(source->crNsigma)) {
-            psTrace("psphot", 7, "Not calculating extNsigma,crNsigma since already measured\n");
-            continue;
+        options.soft = 0.0;
+    }
+
+    // We are using the value psfMag - 2.5*log10(moment.Sum) as a measure of the extendedness
+    // of and 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?
+    psphotSourceSizePSF (&options, sources, psf);
+
+    // classify the sources based on ApResid and Moments (extended sources)
+    psphotSourceClass(readout, sources, recipe, psf, &options);
+
+    psphotSourceSizeCR (readout, sources, &options);
+
+    psLogMsg ("psphot.size", PS_LOG_INFO, "measure source sizes for %ld sources: %f sec\n", sources->n - first, psTimerMark ("psphot.size"));
+
+    psphotVisualPlotSourceSize (recipe, sources);
+    psphotVisualShowSourceSize (readout, sources);
+    psphotVisualPlotApResid (sources, options.ApResid, options.ApSysErr);
+
+    return true;
+}
+
+bool psphotMaskCosmicRay (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {
+
+    // replace the source flux
+    pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+
+    // flag this as a CR
+    source->mode |= PM_SOURCE_MODE_CR_LIMIT;
+    pmPeak *peak = source->peak;
+    psAssert (peak, "NULL peak");
+
+    // grab the matching footprint
+    pmFootprint *footprint = peak->footprint;
+    if (!footprint) {
+        // if we have not footprint, use the old code to mask by isophot
+        psphotMaskCosmicRayIsophot (source, maskVal, crMask);
+        return true;
+    }
+
+    if (!footprint->spans) {
+        // if we have no footprint, use the old code to mask by isophot
+        psphotMaskCosmicRayIsophot (source, maskVal, crMask);
+        return true;
+    }
+
+    // mask all of the pixels covered by the spans of the footprint
+    for (int j = 1; j < footprint->spans->n; j++) {
+        pmSpan *span1 = footprint->spans->data[j];
+
+        int iy = span1->y;
+        int xs = span1->x0;
+        int xe = span1->x1;
+
+        for (int ix = xs; ix < xe; ix++) {
+            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
+        }
+    }
+    return true;
+}
+
+bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {
+
+    source->mode |= PM_SOURCE_MODE_CR_LIMIT;
+    pmPeak *peak = source->peak;
+    psAssert (peak, "NULL peak");
+
+    psImage *mask   = source->maskView;
+    psImage *pixels = source->pixels;
+    psImage *variance = source->variance;
+
+    // XXX This should be a recipe variable
+# define SN_LIMIT 5.0
+
+    int xo = peak->x - pixels->col0;
+    int yo = peak->y - pixels->row0;
+
+    // mark the pixels in this row to the left, then the right
+    for (int ix = xo; ix >= 0; ix--) {
+        float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
+        if (SN > SN_LIMIT) {
+            mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
+        }
+    }
+    for (int ix = xo + 1; ix < pixels->numCols; ix++) {
+        float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
+        if (SN > SN_LIMIT) {
+            mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
+        }
+    }
+
+    // for each of the neighboring rows, mark the high pixels if they have a marked neighbor
+    // first go up:
+    for (int iy = PS_MIN(yo, mask->numRows-2); iy >= 0; iy--) {
+        // mark the pixels in this row to the left, then the right
+        for (int ix = 0; ix < pixels->numCols; ix++) {
+            float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
+            if (SN < SN_LIMIT) continue;
+
+            bool valid = false;
+            valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask);
+            valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0;
+            valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0;
+
+            if (!valid) continue;
+            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
+        }
+    }
+    // next go down:
+    for (int iy = PS_MIN(yo+1, mask->numRows-1); iy < pixels->numRows; iy++) {
+        // mark the pixels in this row to the left, then the right
+        for (int ix = 0; ix < pixels->numCols; ix++) {
+            float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
+            if (SN < SN_LIMIT) continue;
+
+            bool valid = false;
+            valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask);
+            valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0;
+            valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0;
+
+            if (!valid) continue;
+            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
+        }
+    }
+    return true;
+}
+
+// model the apmifit distribution for the psf stars:
+bool psphotSourceSizePSF (psphotSourceSizeOptions *options, psArray *sources, pmPSF *psf) {
+
+    // select stats from the psf stars
+    psVector *Ap = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *ApErr = psVectorAllocEmpty (100, PS_TYPE_F32);
+    
+    psImageMaskType maskVal = options->maskVal | options->markVal;
+
+    // XXX  why PHOT_WEIGHT??
+    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT;
+
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = sources->data[i];
+        if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
+
+        // replace object in image
+        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
+            pmSourceAdd (source, PM_MODEL_OP_FULL, options->maskVal);
         }
 
-        // source must have been subtracted
-        if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) {
+        // 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); // maskVal includes 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);
+
+        float apMag = -2.5*log10(source->moments->Sum);
+        float dMag = source->psfMag - apMag;
+        
+        psVectorAppend (Ap, 100, dMag);
+        psVectorAppend (ApErr, 100, source->errMag);
+    }
+
+    // model the distribution as a mean or median value and a systematic error from that value:
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
+    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->robustMedian;
+    }
+
+    options->ApResid = stats->robustMedian;
+    options->ApSysErr = psVectorSystematicError(dAp, ApErr, 0.05);
+    psLogMsg ("psphot", PS_LOG_DETAIL, "psf - Sum: %f +/- %f\n", options->ApResid, options->ApSysErr);
+
+    psFree (Ap);
+    psFree (ApErr);
+    psFree (stats);
+    psFree (dAp);
+
+    return true;
+}
+
+// classify sources based on the combination of psf-mag, Mxx, Myy
+bool psphotSourceClass (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options) {
+
+    bool status;
+    pmPSFClump psfClump;
+    char regionName[64];
+
+    psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4s %4s %4s %4s %4s %4s", "Npsf", "Next", "Nsat", "Ncr", "Nmiss", "Nskip");
+
+    int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS");
+    for (int i = 0; i < nRegions; i ++) {
+        snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i);
+        psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
+        psAssert (regionMD, "regions must be defined by earlier call to psphotRoughClassRegion");
+
+        psRegion *region = psMetadataLookupPtr (&status, regionMD, "REGION");
+        psAssert (region, "regions must be defined by earlier call to psphotRoughClassRegion");
+
+        // pull FWHM_X,Y from the recipe, use to define psfClump.X,Y
+        psfClump.X  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");   psAssert (status, "missing PSF.CLUMP.X");
+        psfClump.Y  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");   psAssert (status, "missing PSF.CLUMP.Y");
+        psfClump.dX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");  psAssert (status, "missing PSF.CLUMP.DX");
+        psfClump.dY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");  psAssert (status, "missing PSF.CLUMP.DY");
+
+        if ((psfClump.X < 0) || (psfClump.Y < 0) || !psfClump.X || !psfClump.Y || isnan(psfClump.X) || isnan(psfClump.Y)) {
+            psLogMsg ("psphot", 4, "Failed to find a valid PSF clump for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1);
+            continue;
+        }
+        
+        if (!psphotSourceClassRegion (region, &psfClump, sources, recipe, psf, options)) {
+            psLogMsg ("psphot", 4, "Failed to determine source classification for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1);
+            continue;
+        }
+    }   
+
+    return true;
+}
+
+bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options) {
+
+    PS_ASSERT_PTR_NON_NULL(sources, false);
+    PS_ASSERT_PTR_NON_NULL(recipe, false);
+
+    int Nsat  = 0;
+    int Next  = 0;
+    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;
+
+    psImageMaskType maskVal = options->maskVal | options->markVal;
+
+    for (psS32 i = 0 ; i < sources->n ; i++) {
+
+        pmSource *source = (pmSource *) sources->data[i];
+
+        // 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;
+
+        // 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;
+        }
+
+        // source must have been subtracted
+        if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) {
             source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED;
-            psTrace("psphot", 7, "Not calculating extNsigma,crNsigma since source is not subtracted\n");
-            continue;
+            psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n");
+            continue;
+        }
+
+        // we are basically classifying by moments; use the default if not found
+        psAssert (source->moments, "why is this source missing moments?");
+        if (source->mode & noMoments) { 
+            Nskip ++;
+            continue;
+        }
+
+        psF32 Mxx = source->moments->Mxx;
+        psF32 Myy = source->moments->Myy;
+
+        // replace object in image
+        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
+            pmSourceAdd (source, PM_MODEL_OP_FULL, options->maskVal);
         }
 
-        psF32 **resid  = source->pixels->data.F32;
-        psF32 **variance = source->variance->data.F32;
-        psImageMaskType **mask = source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
-
-        // check for extendedness: measure the delta flux significance at the 1 sigma contour
-        source->extNsigma = psphotModelContour(source->pixels, source->variance, source->maskObj, maskVal,
-                                               source->modelPSF, 1.0);
-
-        // XXX prevent a source from being both CR and EXT?
-        if (source->extNsigma > EXT_NSIGMA_LIMIT) {
-            source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
-        }
-
-        // Integer position of peak
-        int xPeak = source->peak->xf - source->pixels->col0 + 0.5;
-        int yPeak = source->peak->yf - source->pixels->row0 + 0.5;
-
-        // XXX for now, skip sources which are too close to a boundary
-        // XXX raise a flag?
-        if (xPeak < 1 || xPeak > source->pixels->numCols - 2 ||
-            yPeak < 1 || yPeak > source->pixels->numRows - 2) {
-            source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED;
-            psTrace("psphot", 7, "Not calculating crNsigma due to edge\n");
-            continue;
-        }
-
-        // XXX for now, just skip any sources with masked pixels
-        // XXX raise a flag?
-        bool keep = true;
-        for (int iy = -1; (iy <= +1) && keep; iy++) {
-            for (int ix = -1; (ix <= +1) && keep; ix++) {
-                if (mask[yPeak+iy][xPeak+ix] & maskVal) {
-                    keep = false;
-                }
-            }
-        }
-        if (!keep) {
-            psTrace("psphot", 7, "Not calculating crNsigma due to masked pixels\n");
-            source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED;
-            continue;
-        }
-
-        // Compare the central pixel with those on either side, for the four possible lines through it.
-
-        // Soften variances (add systematic error)
-        float softening = soft * PS_SQR(source->peak->flux); // Softening for variances
-
-        // Across the middle: y = 0
-        float cX = 2*resid[yPeak][xPeak]   - resid[yPeak+0][xPeak-1]  - resid[yPeak+0][xPeak+1];
-        float dcX = 4*variance[yPeak][xPeak] + variance[yPeak+0][xPeak-1] + variance[yPeak+0][xPeak+1];
-        float nX = cX / sqrtf(dcX + softening);
-
-        // Up the centre: x = 0
-        float cY = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak+0]  - resid[yPeak+1][xPeak+0];
-        float dcY = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak+0] + variance[yPeak+1][xPeak+0];
-        float nY = cY / sqrtf(dcY + softening);
-
-        // Diagonal: x = y
-        float cL = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak-1]  - resid[yPeak+1][xPeak+1];
-        float dcL = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak-1] + variance[yPeak+1][xPeak+1];
-        float nL = cL / sqrtf(dcL + softening);
-
-        // Diagonal: x = - y
-        float cR = 2*resid[yPeak][xPeak]   - resid[yPeak+1][xPeak-1]  - resid[yPeak-1][xPeak+1];
-        float dcR = 4*variance[yPeak][xPeak] + variance[yPeak+1][xPeak-1] + variance[yPeak-1][xPeak+1];
-        float nR = cR / sqrtf(dcR + softening);
-
-        // P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2)
-        // Ndof = 4 ? (four measurements, no free parameters)
-        // XXX this value is going to be biased low because of systematic errors.
-        // we need to calibrate it somehow
-        // source->psfProb = gsl_sf_gamma_inc_Q (2, 0.5*chisq);
-
-        // not strictly accurate: overcounts the chisq contribution from the center pixel (by
-        // factor of 4); also biases a bit low if any pixels are masked
-        // XXX I am not sure I want to keep this value...
-        source->psfChisq = PS_SQR(nX) + PS_SQR(nY) + PS_SQR(nL) + PS_SQR(nR);
-
-        float fCR = 0.0;
-        int nCR = 0;
-        if (nX > 0.0) {
-            fCR += nX;
-            nCR ++;
-        }
-        if (nY > 0.0) {
-            fCR += nY;
-            nCR ++;
-        }
-        if (nL > 0.0) {
-            fCR += nL;
-            nCR ++;
-        }
-        if (nR > 0.0) {
-            fCR += nR;
-            nCR ++;
-        }
-        source->crNsigma  = (nCR > 0)  ? fCR / nCR : 0.0;
-        if (!isfinite(source->crNsigma)) {
-            continue;
-        }
-
-        // this source is thought to be a cosmic ray.  flag the detection and mask the pixels
-        if (source->crNsigma > CR_NSIGMA_LIMIT) {
-            // XXX still testing... : psphotMaskCosmicRay_New (readout->mask, source, maskVal, crMask);
-            psphotMaskCosmicRay_Old (source, maskVal, crMask);
-        }
-    }
-
-    // now that we have masked pixels associated with CRs, we can grow the mask
-    if (grow > 0) {
-        bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading for psImageConvolveMask
-        psImage *newMask = psImageConvolveMask(NULL, readout->mask, crMask, crMask, -grow, grow, -grow, grow);
-        psImageConvolveSetThreads(oldThreads);
-        if (!newMask) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to grow CR mask");
-            return false;
-        }
-        psFree(readout->mask);
-        readout->mask = newMask;
-    }
-
-    psLogMsg ("psphot.size", PS_LOG_INFO, "measure source sizes for %ld sources: %f sec\n",
-              sources->n - first, psTimerMark ("psphot.size"));
-
-    psphotVisualPlotSourceSize (sources);
-    psphotVisualShowSourceSize (readout, sources);
+        // 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); // maskVal includes 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);
+
+        float apMag = -2.5*log10(source->moments->Sum);
+        float dMag = source->psfMag - apMag;
+        float nSigma = (dMag - options->ApResid) / hypot(source->errMag, options->ApSysErr);
+
+        source->extNsigma = nSigma;
+        source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
+
+        // 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(nSigma) < options->nSigmaApResid) && (fabs(Mxx - psfClump->X) < options->nSigmaMoments*psfClump->dX) && (fabs(Myy - psfClump->Y) < options->nSigmaMoments*psfClump->dY);
+        if (isPSF) {
+            Npsf ++;
+            continue;
+        }
+
+        // Defects may not always match CRs from peak curvature analysis
+        // Defects may also be marked as SATSTAR -- XXX deactivate this flag?
+        // XXX this rule is not great
+        if ((Mxx < psfClump->X) || (Myy < psfClump->Y)) {
+            source->mode |= PM_SOURCE_MODE_DEFECT;
+            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) { 
+            Nsat ++;
+            continue;
+        }
+
+        // XXX allow the Mxx, Myy to be less than psfClump->X,Y (by some nSigma)?
+        bool isEXT = (nSigma > options->nSigmaApResid) || ((Mxx > psfClump->X) && (Myy > psfClump->Y));
+        if (isEXT) {
+            source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
+            Next ++;
+            continue;
+        }
+
+        psWarning ("sourse size was missed for %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigma);
+        Nmiss ++;
+    }
+
+    psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4d %4d %4d %4d %4d %4d", Npsf, Next, Nsat, Ncr, Nmiss, Nskip);
 
     return true;
@@ -194 +414 @@
 // given the PSF ellipse parameters, navigate around the 1sigma contour, return the total
 // deviation in sigmas.  This is measured on the residual image - should we ignore negative
-// deviations?
-static float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask,
-                                psImageMaskType maskVal, const pmModel *model, float Ro)
+// deviations?  NOTE: This function was an early attempt to classify extended objects, and is
+// no longer used by psphot.
+float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask,
+                         psImageMaskType maskVal, const pmModel *model, float Ro)
 {
     psF32 *PAR = model->params->data.F32; // Model parameters
@@ -211 +432 @@
     float Q = Ro * PS_SQR(sxx) / (1.0 - PS_SQR(sxx * syy * sxy) / 4.0);
     if (Q < 0.0) {
-        // ellipse is imaginary
-        return NAN;
+        // ellipse is imaginary
+        return NAN;
     }
 
@@ -220 +441 @@
 
     for (int x = -radius; x <= radius; x++) {
-        // Polynomial coefficients
-        // XXX Should we be using the centre of the pixel as x or x+0.5?
-        float A = PS_SQR (1.0 / syy);
-        float B = x * sxy;
-        float C = PS_SQR (x / sxx) - Ro;
-        float T = PS_SQR(B) - 4*A*C;
-        if (T < 0.0) {
-            continue;
-        }
-
-        // y position in source frame
-        float yP = (-B + sqrt (T)) / (2.0 * A);
-        float yM = (-B - sqrt (T)) / (2.0 * A);
-
-        // Get the closest pixel positions (image frame)
-        int xPix  = x  + PAR[PM_PAR_XPOS] - image->col0 + 0.5;
-        int yPixM = yM + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
-        int yPixP = yP + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
-
-        if (xPix < 0 || xPix >= image->numCols) {
-            continue;
-        }
-
-        if (yPixM >= 0 && yPixM < image->numRows &&
-            !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixM][xPix] & maskVal))) {
-            float dSigma = image->data.F32[yPixM][xPix] / sqrtf(variance->data.F32[yPixM][xPix]);
-            nSigma += dSigma;
-            nPts++;
-        }
-
-        if (yPixM == yPixP) {
-            continue;
-        }
-
-        if (yPixP >= 0 && yPixP < image->numRows &&
-            !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixP][xPix] & maskVal))) {
-            float dSigma = image->data.F32[yPixP][xPix] / sqrtf(variance->data.F32[yPixP][xPix]);
-            nSigma += dSigma;
-            nPts++;
-        }
+        // Polynomial coefficients
+        // XXX Should we be using the centre of the pixel as x or x+0.5?
+        float A = PS_SQR (1.0 / syy);
+        float B = x * sxy;
+        float C = PS_SQR (x / sxx) - Ro;
+        float T = PS_SQR(B) - 4*A*C;
+        if (T < 0.0) {
+            continue;
+        }
+
+        // y position in source frame
+        float yP = (-B + sqrt (T)) / (2.0 * A);
+        float yM = (-B - sqrt (T)) / (2.0 * A);
+
+        // Get the closest pixel positions (image frame)
+        int xPix  = x  + PAR[PM_PAR_XPOS] - image->col0 + 0.5;
+        int yPixM = yM + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
+        int yPixP = yP + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
+
+        if (xPix < 0 || xPix >= image->numCols) {
+            continue;
+        }
+
+        if (yPixM >= 0 && yPixM < image->numRows &&
+            !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixM][xPix] & maskVal))) {
+            float dSigma = image->data.F32[yPixM][xPix] / sqrtf(variance->data.F32[yPixM][xPix]);
+            nSigma += dSigma;
+            nPts++;
+        }
+
+        if (yPixM == yPixP) {
+            continue;
+        }
+
+        if (yPixP >= 0 && yPixP < image->numRows &&
+            !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixP][xPix] & maskVal))) {
+            float dSigma = image->data.F32[yPixP][xPix] / sqrtf(variance->data.F32[yPixP][xPix]);
+            nSigma += dSigma;
+            nPts++;
+        }
     }
     nSigma /= nPts;
@@ -265 +486 @@
 }
 
-bool psphotMaskCosmicRay_New (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {
-
-    // replace the source flux
-    pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
-    source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
-
-    // flag this as a CR
-    source->mode |= PM_SOURCE_MODE_CR_LIMIT;
-    pmPeak *peak = source->peak;
-    psAssert (peak, "NULL peak");
-
-    // grab the matching footprint
-    pmFootprint *footprint = peak->footprint;
-    if (!footprint) {
-        // if we have not footprint, use the old code to mask by isophot
-        psphotMaskCosmicRay_Old (source, maskVal, crMask);
-        return true;
-    }
-
-    if (!footprint->spans) {
-        // if we have not footprint, use the old code to mask by isophot
-        psphotMaskCosmicRay_Old (source, maskVal, crMask);
-        return true;
-    }
-
-    // mask all of the pixels covered by the spans of the footprint
-    for (int j = 1; j < footprint->spans->n; j++) {
-        pmSpan *span1 = footprint->spans->data[j];
-
-        int iy = span1->y;
-        int xs = span1->x0;
-        int xe = span1->x1;
-
-        for (int ix = xs; ix < xe; ix++) {
-            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
-        }
+bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) {
+
+    // classify the sources based on the CR test (place this in a function?)
+    // XXX use an internal flag to mark sources which have already been measured
+    for (int i = 0; i < sources->n; i++) {
+        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;
+        }
+
+        // source must have been subtracted
+        if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) {
+            source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED;
+            psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n");
+            continue;
+        }
+
+        psF32 **resid  = source->pixels->data.F32;
+        psF32 **variance = source->variance->data.F32;
+        psImageMaskType **mask = source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
+
+        // Integer position of peak
+        int xPeak = source->peak->xf - source->pixels->col0 + 0.5;
+        int yPeak = source->peak->yf - source->pixels->row0 + 0.5;
+
+        // Skip sources which are too close to a boundary.  These are mostly caught as DEFECT
+        if (xPeak < 1 || xPeak > source->pixels->numCols - 2 ||
+            yPeak < 1 || yPeak > source->pixels->numRows - 2) {
+            psTrace("psphot", 7, "Not calculating crNsigma due to edge\n");
+            continue;
+        }
+
+        // Skip sources with masked pixels.  These are mostly caught as DEFECT
+        bool keep = true;
+        for (int iy = -1; (iy <= +1) && keep; iy++) {
+            for (int ix = -1; (ix <= +1) && keep; ix++) {
+                if (mask[yPeak+iy][xPeak+ix] & options->maskVal) {
+                    keep = false;
+                }
+            }
+        }
+        if (!keep) {
+            psTrace("psphot", 7, "Not calculating crNsigma due to masked pixels\n");
+            continue;
+        }
+
+        // Compare the central pixel with those on either side, for the four possible lines through it.
+
+        // Soften variances (add systematic error)
+        float softening = options->soft * PS_SQR(source->peak->flux); // Softening for variances
+
+        // Across the middle: y = 0
+        float cX = 2*resid[yPeak][xPeak]   - resid[yPeak+0][xPeak-1]  - resid[yPeak+0][xPeak+1];
+        float dcX = 4*variance[yPeak][xPeak] + variance[yPeak+0][xPeak-1] + variance[yPeak+0][xPeak+1];
+        float nX = cX / sqrtf(dcX + softening);
+
+        // Up the centre: x = 0
+        float cY = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak+0]  - resid[yPeak+1][xPeak+0];
+        float dcY = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak+0] + variance[yPeak+1][xPeak+0];
+        float nY = cY / sqrtf(dcY + softening);
+
+        // Diagonal: x = y
+        float cL = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak-1]  - resid[yPeak+1][xPeak+1];
+        float dcL = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak-1] + variance[yPeak+1][xPeak+1];
+        float nL = cL / sqrtf(dcL + softening);
+
+        // Diagonal: x = - y
+        float cR = 2*resid[yPeak][xPeak]   - resid[yPeak+1][xPeak-1]  - resid[yPeak-1][xPeak+1];
+        float dcR = 4*variance[yPeak][xPeak] + variance[yPeak+1][xPeak-1] + variance[yPeak-1][xPeak+1];
+        float nR = cR / sqrtf(dcR + softening);
+
+        // P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2)
+        // Ndof = 4 ? (four measurements, no free parameters)
+        // XXX this value is going to be biased low because of systematic errors.
+        // we need to calibrate it somehow
+        // source->psfProb = gsl_sf_gamma_inc_Q (2, 0.5*chisq);
+
+        // not strictly accurate: overcounts the chisq contribution from the center pixel (by
+        // factor of 4); also biases a bit low if any pixels are masked
+        // XXX I am not sure I want to keep this value...
+        source->psfChisq = PS_SQR(nX) + PS_SQR(nY) + PS_SQR(nL) + PS_SQR(nR);
+
+        float fCR = 0.0;
+        int nCR = 0;
+        if (nX > 0.0) {
+            fCR += nX;
+            nCR ++;
+        }
+        if (nY > 0.0) {
+            fCR += nY;
+            nCR ++;
+        }
+        if (nL > 0.0) {
+            fCR += nL;
+            nCR ++;
+        }
+        if (nR > 0.0) {
+            fCR += nR;
+            nCR ++;
+        }
+        source->crNsigma  = (nCR > 0)  ? fCR / nCR : 0.0;
+        source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
+
+        if (!isfinite(source->crNsigma)) {
+            continue;
+        }
+
+        // this source is thought to be a cosmic ray.  flag the detection and mask the pixels
+        if (source->crNsigma > options->nSigmaCR) {
+            source->mode |= PM_SOURCE_MODE_CR_LIMIT;
+            // XXX still testing... : psphotMaskCosmicRay (readout->mask, source, maskVal, crMask);
+            // XXX acting strange... psphotMaskCosmicRay_Old (source, maskVal, crMask);
+        }
+    }
+
+    // now that we have masked pixels associated with CRs, we can grow the mask
+    if (options->grow > 0) {
+        bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading for psImageConvolveMask
+        psImage *newMask = psImageConvolveMask(NULL, readout->mask, options->crMask, options->crMask, -options->grow, options->grow, -options->grow, options->grow);
+        psImageConvolveSetThreads(oldThreads);
+        if (!newMask) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to grow CR mask");
+            return false;
+        }
+        psFree(readout->mask);
+        readout->mask = newMask;
     }
     return true;
 }
-
-bool psphotMaskCosmicRay_Old (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {
-
-    source->mode |= PM_SOURCE_MODE_CR_LIMIT;
-    pmPeak *peak = source->peak;
-    psAssert (peak, "NULL peak");
-
-    psImage *mask   = source->maskView;
-    psImage *pixels = source->pixels;
-    psImage *variance = source->variance;
-
-    // XXX This should be a recipe variable
-# define SN_LIMIT 5.0
-
-    int xo = peak->x - pixels->col0;
-    int yo = peak->y - pixels->row0;
-
-    // mark the pixels in this row to the left, then the right
-    for (int ix = xo; ix >= 0; ix--) {
-        float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
-        if (SN > SN_LIMIT) {
-            mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
-        }
-    }
-    for (int ix = xo + 1; ix < pixels->numCols; ix++) {
-        float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
-        if (SN > SN_LIMIT) {
-            mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
-        }
-    }
-
-    // for each of the neighboring rows, mark the high pixels if they have a marked neighbor
-    // first go up:
-    for (int iy = PS_MIN(yo, mask->numRows-2); iy >= 0; iy--) {
-        // mark the pixels in this row to the left, then the right
-        for (int ix = 0; ix < pixels->numCols; ix++) {
-            float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
-            if (SN < SN_LIMIT) continue;
-
-            bool valid = false;
-            valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask);
-            valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0;
-            valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0;
-
-            if (!valid) continue;
-            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
-        }
-    }
-    // next go down:
-    for (int iy = PS_MIN(yo+1, mask->numRows-1); iy < pixels->numRows; iy++) {
-        // mark the pixels in this row to the left, then the right
-        for (int ix = 0; ix < pixels->numCols; ix++) {
-            float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
-            if (SN < SN_LIMIT) continue;
-
-            bool valid = false;
-            valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask);
-            valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0;
-            valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0;
-
-            if (!valid) continue;
-            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
-        }
-    }
-    return true;
-}

@@ -18 +18 @@
 psphotVersionDefinitions.h
 psphotMomentsStudy
+psphotPetrosianStudy