Changeset 27838 for branches/tap_branches/psphot/src/psphotSourceSize.c

Timestamp:

May 3, 2010, 8:41:49 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/tap_branches

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src (modified) (1 prop)
• psphot/src/psphotSourceSize.c (modified) (11 diffs)

branches/tap_branches/psphot/src

@@ -19 +19 @@
 psphotMomentsStudy
 psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack

@@ -19 +19 @@
 psphotMomentsStudy
 psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack

branches/tap_branches/psphot/src/psphotSourceSize.c

@@ -13 +13 @@
     float soft;
     int grow;
+    int xtest, ytest;
+    bool apply; // apply CR mask?
 } psphotSourceSizeOptions;
 
+// local functions:
 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);
+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);
 
 // we need to call this function after sources have been fitted to the PSF model and
@@ -28 +32 @@
 // deviation from the psf model at the r = FWHM/2 position
 
-// 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)
+// for now, let's store the detections on the readout->analysis for each readout
+bool psphotSourceSize (pmConfig *config, const pmFPAview *view, bool getPSFsize)
+{
+    bool status = true;
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
+    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
+
+    // skip the chisq image (optionally?)
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+        if (i == chisqNum) continue; // skip chisq image
+        if (!psphotSourceSizeReadout (config, view, "PSPHOT.INPUT", i, recipe, getPSFsize)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed on source size analysis for PSPHOT.INPUT entry %d", i);
+            return false;
+        }
+    }
+    return true;
+}
+
+// this function use an internal flag to mark sources which have already been measured
+bool psphotSourceSizeReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool getPSFsize)
 {
     bool status;
@@ -36 +67 @@
     psTimerStart ("psphot.size");
 
+    // find the currently selected readout
+    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
+    psAssert (file, "missing file?");
+
+    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
+    psAssert (readout, "missing readout?");
+
+    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+    psAssert (detections, "missing detections?");
+
+    psArray *sources = detections->allSources;
+    psAssert (sources, "missing sources?");
+
+    if (!sources->n) {
+        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size");
+        return true;
+    }
+
+    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
+    psAssert (psf, "missing psf?");
+
     // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
     options.maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
@@ -55 +107 @@
     // XXX recipe name is not great
     options.nSigmaMoments = psMetadataLookupF32 (&status, recipe, "PSPHOT.EXT.NSIGMA.MOMENTS");
+    assert (status);
+
+    // XXX recipe name is not great
+    options.xtest = psMetadataLookupS32 (&status, recipe, "PSPHOT.CRMASK.XTEST");
+    options.ytest = psMetadataLookupS32 (&status, recipe, "PSPHOT.CRMASK.YTEST");
     assert (status);
 
@@ -69 +126 @@
     }
 
+    options.apply = 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.");
+        return false;
+    }
+
     // 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
+    // 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?
-    psphotSourceSizePSF (&options, sources, psf);
+
+    // 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)
+    // 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);
 
-    psLogMsg ("psphot.size", PS_LOG_INFO, "measure source sizes for %ld sources: %f sec\n", sources->n - first, psTimerMark ("psphot.size"));
+    // XXX fix this (was source->n  - first)
+    psLogMsg ("psphot.size", PS_LOG_INFO, "measure source sizes for %ld sources: %f sec\n", sources->n, psTimerMark ("psphot.size"));
 
     psphotVisualPlotSourceSize (recipe, readout->analysis, sources);
     psphotVisualShowSourceSize (readout, sources);
     psphotVisualPlotApResid (sources, options.ApResid, options.ApSysErr);
+    psphotVisualShowSatStars (recipe, psf, sources);
 
     return true;
 }
 
-bool psphotMaskCosmicRay (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {
+// 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;
+
+    int num = 0;                        // Number of sources measured
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = sources->data[i];
+        if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
+        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); // 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, dMag);
+        psVectorAppend (ApErr, source->errMag);
+    }
+    if (num == 0) {
+        // Not raising an error, because errors aren't being checked elsewhere in this function
+        psFree(Ap);
+        psFree(ApErr);
+        return false;
+    }
+
+    // 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);
+    // XXX this is quite arbitrary...
+    if (!isfinite(options->ApSysErr)) options->ApSysErr = 0.01;
+    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, readout->analysis, "PSF.CLUMP.NREGIONS");
+    for (int i = 0; i < nRegions; i ++) {
+        snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i);
+        psMetadata *regionMD = psMetadataLookupPtr (&status, readout->analysis, 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;
+        }
+        // psphotVisualPlotSourceSize (recipe, readout->analysis, sources);
+    }
+
+    return true;
+}
+
+# define SIZE_SN_LIM 10
+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 source size since source is not subtracted\n");
+            continue;
+        }
+
+        // we are basically classifying by moments
+        psAssert (source->moments, "why is this source missing moments?");
+        if (source->mode & noMoments) {
+            Nskip ++;
+            continue;
+        }
+
+        // convert to Mmaj, Mmin:
+        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);
+        }
+
+        // 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;
+
+        // set nSigma to include both systematic and poisson error terms
+        // XXX the 'poisson error' contribution for size is probably wrong...
+        float nSigmaMAG = (dMag - options->ApResid) / hypot(source->errMag, options->ApSysErr);
+        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);
+
+        // partially-masked sources are more likely to be mis-measured PSFs
+        float sizeBias = 1.0;
+        if (source->pixWeight < 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?
+        source->extNsigma = nSigmaMAG;
+
+        // notes to clarify the source size classification rules:
+        // * a defect should be functionally equivalent to a cosmic ray
+        // * 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
+        // Defects may also be marked as SATSTAR -- XXX deactivate this flag?
+        // XXX this rule is not great
+        // XXX only accept brightish detections as CRs
+        // (nSigmaMAG < -options->nSigmaApResid) ||
+        bool isCR = isCR = (source->errMag < 1.0 / SIZE_SN_LIM) && ((Mxx < minMxx) || (Myy < minMyy));
+        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->pixWeight,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
+                    options->nSigmaApResid,sizeBias*options->nSigmaMoments);
+            source->mode |= PM_SOURCE_MODE_DEFECT;
+            source->tmpFlags |= PM_SOURCE_TMPF_SIZE_CR_CANDIDATE;
+            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);
+        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);
+
+            source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
+            source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
+            Next ++;
+            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);
+
+    return true;
+}
+
+// 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) {
+
+    psTimerStart ("psphot.cr");
+
+    int nMasked = 0;
+    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;
+        }
+
+        // only check candidates marked above
+        if (!(source->tmpFlags & PM_SOURCE_TMPF_SIZE_CR_CANDIDATE)) {
+            psTrace("psphot", 7, "Not calculating source size since it has already been measured\n");
+            continue;
+        }
+
+        // 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;
+
+        // 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;
+        }
+
+        // XXX for testing, only CRMASK a single source:
+        if (options->xtest && (fabs(source->peak->xf - options->xtest) > 5)) continue;
+        if (options->ytest && (fabs(source->peak->yf - options->ytest) > 5)) continue;
+
+        // replace object in image
+        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
+            pmSourceAdd (source, PM_MODEL_OP_FULL, options->maskVal);
+        }
+
+        // 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) {
+            psphotMaskCosmicRay(readout, source, options->crMask);
+        } else {
+            source->mode |= PM_SOURCE_MODE_CR_LIMIT;
+        }
+        nMasked ++;
+
+        // re-subtract the object, leave local sky
+        pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal);
+    }
+
+    // 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;
+    }
+
+    psLogMsg ("psphot.cr", PS_LOG_INFO, "mask CR: %d masked in %f sec\n", nMasked, psTimerMark ("psphot.cr"));
+
+    // XXX test : save the mask image
+    if (0) {
+        psphotSaveImage (NULL, readout->mask,   "mask.fits");
+    }
+
+    return true;
+}
+
+# define DUMPPICS 0
+# define LIMIT_XRANGE(X, IMAGE) { X = PS_MIN(PS_MAX(0, X), IMAGE->numCols); }
+# define LIMIT_YRANGE(Y, IMAGE) { Y = PS_MIN(PS_MAX(0, Y), IMAGE->numRows); }
+
+// Comments by CZW 20091209 : Mechanics of how to identify CR pixels taken from "Cosmic-Ray
+// Rejection by Laplacian Edge Detection" by Pieter van Dokkum, arXiv:astro-ph/0108003.  This
+// does no repair or recovery of the CR pixels, it only masks them out.  My test code can be
+// found at /data/ipp031.0/watersc1/psphot.20091209/algo_check.c
+bool psphotMaskCosmicRay (pmReadout *readout, pmSource *source, psImageMaskType maskVal) {
+
+    // Get the actual images and information about the peak.
+    psImage *mask = readout->mask;
+    pmPeak *peak = source->peak;
+    pmFootprint *footprint = peak->footprint;
+
+    // Bounding boxes are inclusive of final pixel
+    int xs = footprint->bbox.x0;
+    int xe = footprint->bbox.x1 + 1;
+    int ys = footprint->bbox.y0;
+    int ye = footprint->bbox.y1 + 1;
+
+    LIMIT_XRANGE(xs, mask);
+    LIMIT_XRANGE(xe, mask);
+    LIMIT_YRANGE(ys, mask);
+    LIMIT_YRANGE(ye, mask);
+
+    int dx = xe - xs;
+    int dy = ye - ys;
+
+    psImage *image= readout->image;
+    psImage *variance = readout->variance;
+
+    int binning = 2;
+    float sigma_thresh = 3.0;
+    int max_iter = 1; // XXX with isophot masking, we only want to do a single pass
+
+    // Temporary images.
+    psImage *mypix  = psImageAlloc(dx,dy,image->type.type);
+    psImage *myfix  = psImageAlloc(dx,dy,image->type.type);
+    psImage *myvar  = psImageAlloc(dx,dy,image->type.type);
+    psImage *binned = psImageAlloc(dx * binning,dy * binning,image->type.type);
+    psImage *conved = psImageAlloc(dx * binning,dy * binning,image->type.type);
+    psImage *edges  = psImageAlloc(dx,dy,image->type.type);
+    psImage *mymask = psImageAlloc(dx,dy,PS_TYPE_IMAGE_MASK);
+
+    // Load my copy of things.
+    for (int y = 0; y < dy; y++) {
+        for (int x = 0; x < dx; x++) {
+            mypix->data.F32[y][x] = image->data.F32[y+ys][x+xs];
+            myvar->data.F32[y][x] = variance->data.F32[y+ys][x+xs];
+            mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0x00;
+        }
+    }
+    // Mask so I can see on the output image where the footprint is.
+    for (int i = 0; i < footprint->spans->n; i++) {
+        pmSpan *sp = footprint->spans->data[i];
+        for (int j = sp->x0; j <= sp->x1; j++) {
+            int y = sp->y - ys;
+            int x = j - xs;
+            mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= 0x01;
+        }
+    }
+
+    int nCRpix = 1; // force at least one pass...
+    for (int iteration = 0; (iteration < max_iter) && (nCRpix > 0); iteration++) {
+        nCRpix = 0;
+        psImageInit (binned, 0.0);
+        psImageInit (conved, 0.0);
+        psImageInit (edges, 0.0);
+
+        // Make subsampled image. Maybe this should be called "unbinned" or something
+        for (int y = 0; y < binning * dy; y++) {
+            int yraw = y / binning;
+            for (int x = 0; x < binning * dx; x++) {
+                int xraw = x / binning;
+                binned->data.F32[y][x] = mypix->data.F32[yraw][xraw];
+            }
+        }
+
+        // Apply Laplace transform (kernel = [[0 -0.25 0][-0.25 1 -0.25][0 -0.25 0]]), clipping at zero
+        for (int y = 1; y < binning * dy - 1; y++) {
+            for (int x = 1; x < binning * dx - 1; x++) {
+                float value = binned->data.F32[y][x] - 0.25 *
+                    (binned->data.F32[y+0][x-1] + binned->data.F32[y+0][x+1] +
+                     binned->data.F32[y-1][x+0] + binned->data.F32[y+1][x+0]);
+                value = PS_MAX(0.0, value);
+
+                conved->data.F32[y][x] = value;
+            }
+        }
+
+        // Create an edge map by rebinning
+        for (int y = 0; y < binning * dy; y++) {
+            int yraw = y / binning;
+            for (int x = 0; x < binning * dx; x++) {
+                int xraw = x / binning;
+                edges->data.F32[yraw][xraw] += conved->data.F32[y][x];
+            }
+        }
+
+        // coordinate of peak in subimage pixels:
+        int xPeak = peak->x - xs;
+        int yPeak = peak->y - ys;
+
+        // Modify my mask if we're above the significance threshold, but only for connected pixels
+        nCRpix = psphotMaskCosmicRayConnected (xPeak, yPeak, mymask, myvar, edges, binning, sigma_thresh);
+
+# if DUMPPICS
+        psphotSaveImage (NULL, mypix,   "crmask.pix.fits");
+# endif
+
+// XXX do not repair the pixels in isophot version
+# if 0
+        // "Repair" Masked pixels for the next round.
+        for (int y = 1; y < dy - 1; y++) {
+            for (int x = 1; x < dx - 1; x++) {
+                if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40)) {
+                    myfix->data.F32[y][x] = mypix->data.F32[y][x];
+                    continue;
+                }
+                myfix->data.F32[y][x] = 0.25 *
+                    (mypix->data.F32[y+0][x-1] + mypix->data.F32[y+0][x+1] +
+                     mypix->data.F32[y-1][x+0] + mypix->data.F32[y+1][x+0]);
+            }
+        }
+
+        // "Repair" Masked pixels for the next round.
+        for (int y = 1; y < dy - 1; y++) {
+            for (int x = 1; x < dx - 1; x++) {
+                mypix->data.F32[y][x] = myfix->data.F32[y][x];
+            }
+        }
+# endif
+
+# if DUMPPICS
+        fprintf (stderr, "CRMASK %d %d %d %d %d\n", xs, ys, dx, dy, iteration);
+        psphotSaveImage (NULL, mypix,   "crmask.fix.fits");
+        psphotSaveImage (NULL, myvar,   "crmask.var.fits");
+        psphotSaveImage (NULL, binned,  "crmask.binn.fits");
+        psphotSaveImage (NULL, conved,  "crmask.conv.fits");
+        psphotSaveImage (NULL, edges,   "crmask.edge.fits");
+        psphotSaveImage (NULL, mymask,  "crmask.mask.fits");
+# endif
+        psTrace("psphot.czw",2,"Iter: %d Count: %d",iteration, nCRpix);
+    }
+
+# if 0
+    // A solitary masked pixel is likely a lie. Remove those
+    // XXX can't we use nCRpix == 1 to test for these?
+    for (int x = 0; x < dx; x++) {
+        for (int y = 0; y < dy; y++) {
+            if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40)) continue;
+            if ((x-1 >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x-1] & 0x40)) {
+                continue;
+            }
+            if ((y-1 >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y-1][x] & 0x40)) {
+                continue;
+            }
+            if ((x+1 < dx) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x+1] & 0x40)) {
+                continue;
+            }
+            if ((y+1 < dy) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y+1][x] & 0x40)) {
+                continue;
+            }
+            mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] ^= 0x40;
+        }
+    }
+# endif
+
+    // transfer temporary mask to real mask & count masked pixels
+    nCRpix = 0;
+    for (int x = 0; x < dx; x++) {
+        for (int y = 0; y < dy; y++) {
+            if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) {
+                mask->data.PS_TYPE_IMAGE_MASK_DATA[y+ys+mask->row0][x+xs+mask->col0] |= maskVal;
+                nCRpix ++;
+            }
+        }
+    }
+
+    // XXX if we decide this REALLY is a cosmic ray, set the CR_LIMIT bit
+    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);
+
+    psFree(mypix);
+    psFree(myfix);
+    psFree(myvar);
+    psFree(binned);
+    psFree(conved);
+    psFree(edges);
+    psFree(mymask);
+
+    return true;
+}
+
+bool psphotMaskCosmicRayFootprintCheck (psArray *sources) {
+
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = sources->data[i];
+        pmPeak *peak = source->peak;
+        pmFootprint *footprint = peak->footprint;
+        if (!footprint) continue;
+        for (int j = 0; j < footprint->spans->n; j++) {
+            pmSpan *sp = footprint->spans->data[j];
+            psAssert (sp, "missing span");
+        }
+    }
+    return true;
+}
+
+/**** ------ old versions of cosmic ray masking ----- ****/
+
+bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask);
+
+// This attempt to mask the cosmic rays used the isophotal boundary
+bool psphotMaskCosmicRay_V1 (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {
 
     // replace the source flux
@@ -103 +748 @@
     pmFootprint *footprint = peak->footprint;
     if (!footprint) {
+      psTrace("psphot.czw",2,"Using isophot CR mask code.");
+
         // if we have not footprint, use the old code to mask by isophot
         psphotMaskCosmicRayIsophot (source, maskVal, crMask);
@@ -109 +756 @@
 
     if (!footprint->spans) {
+      psTrace("psphot.czw",2,"Using isophot CR mask code.");
+
         // if we have no footprint, use the old code to mask by isophot
         psphotMaskCosmicRayIsophot (source, maskVal, crMask);
         return true;
     }
-
+    psphotMaskCosmicRayIsophot (source, maskVal, crMask);
     // mask all of the pixels covered by the spans of the footprint
     for (int j = 1; j < footprint->spans->n; j++) {
@@ -127 +776 @@
     }
     return true;
+}
+
+# define VERBOSE 0
+int psphotMaskCosmicRayConnected (int xPeak, int yPeak, psImage *mymask, psImage *myvar, psImage *edges, int binning, float sigma_thresh) {
+
+    int xLo, xRo;
+    int nCRpix = 0;
+
+    float noise_factor = 5.0 / 4.0;  // Intrinsic to the Laplacian making noise spikes spikier.
+
+    // mark the pixels in this row to the left, then the right. stay within footprint
+    int xL = xPeak; // find the range of valid pixels in this row
+    int xR = xPeak;
+    for (int ix = xPeak; (ix >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] & 0x01); ix--) {
+        float noise = binning * sqrt(noise_factor * myvar->data.F32[yPeak][ix]);
+        float value = edges->data.F32[yPeak][ix] / noise;
+        if (value < sigma_thresh ) break;
+        mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] |= 0x40;
+        xL = ix;
+        nCRpix ++;
+        if (VERBOSE) fprintf (stderr, "mark %d,%d (%d) : %d - %d\n", ix, yPeak, nCRpix, xL, xR);
+    }
+    for (int ix = xPeak; (ix < mymask->numCols) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] & 0x01); ix++) {
+        float noise = binning * sqrt(noise_factor * myvar->data.F32[yPeak][ix]);
+        float value = edges->data.F32[yPeak][ix] / noise;
+        if (value < sigma_thresh ) break;
+        mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] |= 0x40;
+        xR = ix;
+        nCRpix ++;
+        if (VERBOSE) fprintf (stderr, "mark %d,%d (%d) : %d - %d\n", ix, yPeak, nCRpix, xL, xR);
+    }
+    // xL and xR mark the first and last valid pixel in the row
+
+    // for each of the neighboring rows, mark the high pixels if they touch the range xL to xR
+    xLo = PS_MAX(xL - 1, 0);
+    xRo = PS_MIN(xR + 1, mymask->numCols);
+
+    // first go down:
+    for (int iy = yPeak - 1; iy >= 0; iy--) {
+
+        int xLn = -1;
+        int xRn = -1;
+        int newPix = 0;
+
+        // mark the pixels in the good range
+        for (int ix = xLo; ix < xRo; ix++) {
+            if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01)) continue; // only use pixels in the footprint
+            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
+            float value = edges->data.F32[iy][ix] / noise;
+            if (value < sigma_thresh ) continue;
+            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
+            if (xLn == -1) xLn = ix; // first valid pixel in this row
+            xRn = ix;                // last valid pixel in this row
+            nCRpix ++;
+            newPix ++;
+            if (VERBOSE) fprintf (stderr, "mark C %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
+        }
+
+        // mark the pixels to the left of the good range
+        for (int ix = xLo; (ix >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix--) {
+            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
+            float value = edges->data.F32[iy][ix] / noise;
+            if (value < sigma_thresh ) break;
+            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
+            if (xRn == -1) xRn = ix; // last valid pixel in this row
+            xLn = ix;
+            nCRpix ++;
+            newPix ++;
+            if (VERBOSE) fprintf (stderr, "mark L %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
+        }
+
+        // mark the pixels to the right of the good range
+        for (int ix = xRo; (ix < mymask->numCols) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix++) {
+            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
+            float value = edges->data.F32[iy][ix] / noise;
+            if (value < sigma_thresh ) break;
+            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
+            if (xLn == -1) xLn = ix; // first valid pixel in this row
+            xRn = ix;
+            nCRpix ++;
+            newPix ++;
+            if (VERBOSE) fprintf (stderr, "mark R %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
+        }
+        if (newPix == 0) break;
+        xLo = PS_MAX(xLn - 1, 0);
+        xRo = PS_MIN(xRn + 1, mymask->numCols);
+    }
+
+    xLo = PS_MAX(xL - 1, 0);
+    xRo = PS_MIN(xR + 1, mymask->numCols);
+
+    // next go up:
+    for (int iy = yPeak + 1; iy < mymask->numRows; iy++) {
+
+        int xLn = -1;
+        int xRn = -1;
+        int newPix = 0;
+
+        // mark the pixels in the good range
+        for (int ix = xLo; ix < xRo; ix++) {
+            if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01)) continue; // only use pixels in the footprint
+            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
+            float value = edges->data.F32[iy][ix] / noise;
+            if (value < sigma_thresh ) continue;
+            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
+            if (xLn == -1) xLn = ix; // first valid pixel in this row
+            xRn = ix;                // last valid pixel in this row
+            nCRpix ++;
+            newPix ++;
+            if (VERBOSE) fprintf (stderr, "mark C %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
+        }
+
+        // mark the pixels to the left of the good range
+        for (int ix = xLo; (ix >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix--) {
+            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
+            float value = edges->data.F32[iy][ix] / noise;
+            if (value < sigma_thresh ) break;
+            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
+            if (xRn == -1) xRn = ix; // last valid pixel in this row
+            xLn = ix;
+            nCRpix ++;
+            newPix ++;
+            if (VERBOSE) fprintf (stderr, "mark L %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
+        }
+
+        // mark the pixels to the right of the good range
+        for (int ix = xRo; (ix < mymask->numCols) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix++) {
+            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
+            float value = edges->data.F32[iy][ix] / noise;
+            if (value < sigma_thresh ) break;
+            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
+            if (xLn == -1) xLn = ix; // first valid pixel in this row
+            xRn = ix;
+            nCRpix ++;
+            newPix ++;
+            if (VERBOSE) fprintf (stderr, "mark R %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
+        }
+        if (newPix == 0) break;
+        xLo = PS_MAX(xLn - 1, 0);
+        xRo = PS_MIN(xRn + 1, mymask->numCols);
+    }
+
+    return nCRpix;
 }
 
@@ -192 +984 @@
         }
     }
-    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);
-        }
-
-        // 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, readout->analysis, "PSF.CLUMP.NREGIONS");
-    for (int i = 0; i < nRegions; i ++) {
-        snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i);
-        psMetadata *regionMD = psMetadataLookupPtr (&status, readout->analysis, 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 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);
-        }
-
-        // 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;
 }
@@ -486 +1061 @@
 }
 
-bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) {
+// this was an old attempt to identify cosmic rays based on the peak curvature
+bool psphotSourcePeakCurvature (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) {
 
     // classify the sources based on the CR test (place this in a function?)
@@ -618 +1194 @@
     return true;
 }
+

@@ -19 +19 @@
 psphotMomentsStudy
 psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack