Changeset 25755

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

Files:

• doc/notes.20090523.txt (modified) (2 diffs)
• doc/notes.txt (modified) (1 diff)
• doc/outline.txt (modified) (1 diff)
• src (modified) (1 prop)
• src/Makefile.am (modified) (5 diffs)
• src/psphot.h (modified) (8 diffs)
• src/psphotAddNoise.c (modified) (2 diffs)
• src/psphotApResid.c (modified) (15 diffs)
• src/psphotBlendFit.c (modified) (3 diffs)
• src/psphotChoosePSF.c (modified) (8 diffs)
• src/psphotEfficiency.c (modified) (1 prop)
• src/psphotEllipticalContour.c (copied) (copied from branches/eam_branches/20090715/psphot/src/psphotEllipticalContour.c )
• src/psphotEllipticalProfile.c (copied) (copied from branches/eam_branches/20090715/psphot/src/psphotEllipticalProfile.c )
• src/psphotExtendedSourceAnalysis.c (modified) (8 diffs)
• src/psphotExtendedSourceFits.c (modified) (2 diffs)
• src/psphotFake.c (modified) (1 prop)
• src/psphotFindDetections.c (modified) (1 diff)
• src/psphotFindFootprints.c (modified) (2 diffs)
• src/psphotFindPeaks.c (modified) (1 diff)
• src/psphotFitSourcesLinear.c (modified) (5 diffs)
• src/psphotFitSourcesLinearStack.c (modified) (1 diff)
• src/psphotGuessModels.c (modified) (1 diff)
• src/psphotImageLoop.c (modified) (1 diff)
• src/psphotImageQuality.c (modified) (1 diff)
• src/psphotMagnitudes.c (modified) (4 diffs)
• src/psphotMakeFluxScale.c (modified) (1 diff)
• src/psphotMakeResiduals.c (modified) (8 diffs)
• src/psphotMaskReadout.c (modified) (1 diff)
• src/psphotOutput.c (modified) (3 diffs)
• src/psphotPSFConvModel.c (modified) (2 diffs)
• src/psphotPetrosian.c (modified) (1 diff)
• src/psphotPetrosianAnalysis.c (copied) (copied from branches/eam_branches/20090715/psphot/src/psphotPetrosianAnalysis.c )
• src/psphotPetrosianProfile.c (copied) (copied from branches/eam_branches/20090715/psphot/src/psphotPetrosianProfile.c )
• src/psphotPetrosianRadialBins.c (copied) (copied from branches/eam_branches/20090715/psphot/src/psphotPetrosianRadialBins.c )
• src/psphotPetrosianStats.c (copied) (copied from branches/eam_branches/20090715/psphot/src/psphotPetrosianStats.c )
• src/psphotPetrosianStudy.c (copied) (copied from branches/eam_branches/20090715/psphot/src/psphotPetrosianStudy.c )
• src/psphotPetrosianVisual.c (copied) (copied from branches/eam_branches/20090715/psphot/src/psphotPetrosianVisual.c )
• src/psphotRadialProfile.c (modified) (2 diffs)
• src/psphotRadialProfileByAngles.c (copied) (copied from branches/eam_branches/20090715/psphot/src/psphotRadialProfileByAngles.c )
• src/psphotRadiiFromProfiles.c (copied) (copied from branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c )
• src/psphotRadiusChecks.c (modified) (8 diffs)
• src/psphotReadout.c (modified) (9 diffs)
• src/psphotReadoutFindPSF.c (modified) (1 diff)
• src/psphotReadoutKnownSources.c (modified) (1 diff)
• src/psphotReadoutMinimal.c (modified) (1 diff)
• src/psphotReplaceUnfit.c (modified) (5 diffs)
• src/psphotRoughClass.c (modified) (3 diffs)
• src/psphotSetThreads.c (modified) (2 diffs)
• src/psphotSourceFits.c (modified) (19 diffs)
• src/psphotSourcePlots.c (modified) (3 diffs)
• src/psphotSourceSize.c (modified) (6 diffs)
• src/psphotSourceStats.c (modified) (7 diffs)
• src/psphotVisual.c (modified) (68 diffs)

trunk/psphot/doc/notes.20090523.txt

@@ +1 @@
+
+20090809 : more on extended sources:
+
+ my algorithm for getting the petrosian radii and fluxes is:
+
+  * measure radial profile by interpolation to specific locations along the radial line
+    (at low radii, this gives much more accurate results; at high radii, we may need to 
+    average a group of pixels -- say, for 15 deg separations, choose a box that is 
+    < 5 degree in width -- this is > 1 pixel at a distance of 12 pixels.
+
+    psphotRadialProfilesByAngle
+
+  * find intersection of profile with isophot
+
+    psphotRadiusFromProfile -> r50(theta)
+
+  * generate r50x, r50y and fit to ellipse
+
+    psphotEllipticalContour
+
+  * generate elliptical profile
+
+    psphotEllipticalProfile
+
+  * convert to \alpha r_i < r < \beta r_i bins
+
+    psphotPetrosianRadialBins
 
 20090725 : extended source radial profiles
@@ -54 +81 @@
     *** get rest of math from my notebook...
 
-
-
 20090525 : some clarity of issues:
 

trunk/psphot/doc/notes.txt

@@ +1 @@
+
+2009.09.26
+
+  Remaining PSPHOT Issues:
+
+  o soften errors a bit for bright stars?  (essentially allow a
+    systematic error on the flat of some percentage).  this may help
+    the very saturated stars (peaks > 10k) which seem to be
+    excessively driven by the cores.
+    ** did not seem to be a big success --> makes things worse for f =
+    0.01, about the same for f = 0.005
+
+  o clean up and remove test bits of code.
+
+  o clean up psphot.config files : drop unused, fix some conflicting
+    names, reduce camera-specific options, make S/N limits
+    consistent...
+
+  * PSF Eval & ChiSq test??
+
+  * EXT fail fit?  lack of iterations?
+
+  * extended source output
+
+  * still some concerns with the source size analysis: the parameters
+    space seems fine (Mxx, Myy, nSigma), but we need to define more
+    sophisticated boundaries (perhaps include dMag or Mag).  
+
+  * some optimization issues:         
+    * too many re-calculations of all magnitude / aperture data
+    * is FluxScale really faster?  do we need a finer mapping? (5x5
+      clearly yields too much scatter)
+    * is psphotFitSourcesLinear.c (build matrix) super-slow when optimized?
+
+  * very saturated stars are somewhat segmented.
+
+  * footprint coverage in dense simtest image 005.001?
+
+  * footprint S/N limit vs peaks S/N limit?
+
+  * warning about burntool table?
+
+  * is the residual map failing to be masked at the radius?
+
+  * PSF 2D convergence : how to choose between marginally different
+    options?  (note that the GPC1 images may go to 2x2 instead of 3x3)
+
+  * -visual needs more fine-grained control
+
+  * visualization graphs should reset the sections in general.
+
+  * some excess verbosity (my log level is 9, lower?)
+
+  ----- 
+
+  Apparently I'm getting bad psf vs ap mags because of the difference
+  between linear and non-linear fits.  I think this is due to the
+  weighting schemes:
+
+  psphotChoosePSF & psphotBlendFit depend on POISSON.ERRORS.PHOT.LMM,
+  which is set to TRUE => use Poisson errors
+
+  psphotFitSourcesLinear depends on CONSTANT_PHOTOMETRIC_WEIGHTS,
+  which is set to TRUE => do not use Poisson errors.
+
+  nothing seems to use the value POISSON.ERRORS.PHOT.LIN
+
+  ** Poisson errors seeem to be working much better than Constant
+     errors (as judged by consistency between PSF and APERTURE
+     magnitdues).  Not sure I understand this.  It could be that the
+     implementation of Constant errors in the linear fitting analysis
+     is busted.  
+
+  ** Also, I was getting some additional scatter from the FluxScale
+     (lookup-table for Io -> mags).  
+
+     * surprisingly, it is not clear that the FluxScale method is
+       faster than the more accurate integration technique...
+
+2009.09.25
+
+  Notes on pmSourceMagnitudes:
+
+  * in psphotChoosePSF, in pmPSFtry: raw PSF & raw AP
+  * in psphotChoosePSF, psphotMakeGrowthCurve: uses direct all to pmSourcePhotometryAper
+  * in psphotSourceSize, in psphotSourceSizePSF, raw PSF & raw AP (but uses moments->Sum) for PSF stars
+    NOTE: on the initial call, the PSF stars may have had multiple models tried and the psfMag
+    may not be the value for the chosen model.  
+  * in psphotSourceSize, in psphotSourceClassRegion, raw PSF & raw AP for all sources for which size is not yet measured.
+  * in psphotApResid: raw PSF & raw AP for STARS only
+  * in psphotMagnitudes; corrected PSF & AP mags
 
 2009.01.24

trunk/psphot/doc/outline.txt

@@ -2 +2 @@
 Here is a summary outline of the steps taken by psphotReadout:
 
-* setup (choose recipe, choose readout, define break-point)
-* create mask and weight images if needed, set mask for analysis region
-* psphotImageMedian : construct a background model image and subtract it
-* psphotFindPeaks   : smooth and find the peak pixels from the smoothed image
-  - no pmSource defined yet, only pmPeak
-* psphotSourceStats : create sources for each peak and measure their moments
-  - pmSource defined here, with pixels covering SKY_OUTER_RADIUS
-  - no pmModel is defined yet.
-* psphotBasicDeblend : identify blended sources by proximity and valley depth
-* psphotRoughClass : source classification guess based on moments
-  - set the source.type,mode.
-  - re-calculate moments for saturated stars
-* psphotChoosePSF : define the psf model
-  - this generates model fits to the identified psf stars
-    (one for each tested model)
-  - these models are not kept (seems like a waste, but later fits must be consistent)
-* psphotGuessModels : set the initial PSF model for each object (even EXTENDED)
-  - creates modelPSF entries
+* setup (choose recipe, readout, etc)
+
+* create mask and weight images if needed
+
+* construct a background model image and subtract it
+
+* generate the significance image (smoothed by small Gaussian)
+
+* find peaks and associated footprints
+
+* create sources for each peak and measure moments
+  - scan over several window sigma values and choose an optimal window size
+  - at this stage, the moments are aimed at identifying psf-like objects
+
+* identify blended sources by proximity and valley depth
+
+* crude source classification guess based on moments & saturated pixels
+  - major goal is to identify the psf-stars
+  - allows for 2D variations in the psf 
+
+* use the selected PSF candidate stars to generate a PSF model  
+  - fits PSF model parameters as a function of position
+  - uses psfMag - apMag to choose order of 2D  (minimize systematic error in this value)
+  - generates a PSF residual image & optional 2D variations
+
+* fit all detected sources to the psf model (linear fit for normalization only)
+
+* high-quality source classification
+  - uses Mxx, Myy, psfMag - moment->Sum (equiv to apMags) 
+  * should the extended / psf cut be a function of galactic latitude?
+  * should the cosmic ray / psf cut be a function of galactic latitude?
+
+* non-linear fit for all brighter sources to single psf, double psf, or extended source model
+  - sources identified as extended above are fitted with extended source model
+  * are group / clumps of sources being fitted in the best way?
+  * is the psf-model ap mag being measured in an appropriate-size aperture?
+  * is the criterion for choosing between 2 psfs and extended source OK?
+  * is the test for double psf OK? (uses moments to guess starting positions)
+  
+* subtract, re-detect, measure faint sources (PSF-only)
+
+* optional extended source aperture-like measurements (petrosian, etc)
+
+* optional extended source non-linear fits (Sersic, etc)
+
+* measure psfMag-apMag correcion (2D)
+  * convert analysis to use systematic error measurement 
+
+* generate magnitudes
+
+* measure detection efficiency
+
+* output

trunk/psphot/src

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

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

trunk/psphot/src/Makefile.am

@@ -25 +25 @@
 libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
 
-bin_PROGRAMS = psphot psphotTest psphotMomentsStudy
+bin_PROGRAMS = psphot psphotTest psphotMomentsStudy 
+# bin_PROGRAMS = psphotPetrosianStudy 
+# bin_PROGRAMS = psphot
 
 psphot_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
@@ -38 +40 @@
 psphotMomentsStudy_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
 psphotMomentsStudy_LDADD = libpsphot.la
+
+# psphotPetrosianStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
+# psphotPetrosianStudy_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
+# psphotPetrosianStudy_LDADD = libpsphot.la
 
 psphot_SOURCES = \
@@ -61 +67 @@
 psphotMomentsStudy_SOURCES = \
         psphotMomentsStudy.c
+
+# psphotPetrosianStudy_SOURCES = \
+#         psphotPetrosianStudy.c
 
 libpsphot_la_SOURCES = \
@@ -104 +113 @@
         psphotExtendedSourceAnalysis.c \
         psphotExtendedSourceFits.c     \
-        psphotRadialProfile.c          \
-        psphotPetrosian.c              \
-        psphotIsophotal.c              \
-        psphotAnnuli.c                 \
-        psphotKron.c                   \
         psphotKernelFromPSF.c          \
         psphotPSFConvModel.c           \
@@ -129 +133 @@
         psphotThreadTools.c            \
         psphotAddNoise.c               \
+        psphotRadialProfile.c          \
+        psphotRadialProfileByAngles.c  \
+        psphotRadiiFromProfiles.c      \
+        psphotEllipticalContour.c      \
+        psphotEllipticalProfile.c      \
+        psphotPetrosian.c              \
+        psphotPetrosianRadialBins.c    \
+        psphotPetrosianStats.c         \
         psphotEfficiency.c
 
-# dropped? psphotGrowthCurve.c
+# re-instate these
+#       psphotIsophotal.c              \
+#       psphotAnnuli.c                 \
+#       psphotKron.c                   \
+#       psphotPetrosianVisual.c        \
+#
+
+# test versions
+#       psphotPetrosianProfile.c       \
+#       psphotPetrosianAnalysis.c      \
+#
 
 include_HEADERS = \

trunk/psphot/src/psphot.h

@@ -58 +58 @@
 bool            psphotBlendFit_Threaded (psThreadJob *job);
 
-psArray        *psphotSourceStats (pmConfig *config, pmReadout *readout, pmDetections *detections);
+psArray        *psphotSourceStats (pmConfig *config, pmReadout *readout, pmDetections *detections, bool setWindow);
 bool            psphotSourceStats_Threaded (psThreadJob *job);
 
@@ -71 +71 @@
 
 bool            psphotApResid (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf);
-bool            psphotApResidMags_Threaded (psThreadJob *job);
 
 bool            psphotSkyReplace (pmConfig *config, const pmFPAview *view);
@@ -87 +86 @@
 psImage        *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, const int pass, psImageMaskType maskVal);
 psArray        *psphotFindPeaks (psImage *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax);
-bool            psphotFindFootprints (pmDetections *detections, psImage *significance, pmReadout *readout, psMetadata *recipe, const int pass, psImageMaskType maskVal);
+bool            psphotFindFootprints (pmDetections *detections, psImage *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int pass, psImageMaskType maskVal);
 psErrorCode     psphotCullPeaks(const psImage *img, const psImage *weight, const psMetadata *recipe, psArray *footprints);
 
@@ -94 +93 @@
 
 // used by ApResid
-bool            psphotMagErrorScale (float *errorScale, float *errorFloor, psVector *dMag, psVector *dap, psVector *mask, int nGroup);
-bool            psphotApResidTrend (pmReadout *readout, pmPSF *psf, int Npsf, int scale, float *errorScale, float *errorFloor, psVector *mask, psVector *xPos, psVector *yPos, psVector *apResid, psVector *dMag);
+pmTrend2D      *psphotApResidTrend (float *apResidSysErr, pmReadout *readout, int Nx, int Ny, psVector *xPos, psVector *yPos, psVector *apResid, psVector *dMag);
+bool            psphotApResidMags_Threaded (psThreadJob *job);
 
 // basic support functions
@@ -109 +108 @@
 bool            psphotInitRadiusEXT (psMetadata *recipe, pmModelType type);
 bool            psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal);
+float           psphotSetRadiusEXT (pmReadout *readout, pmSource *source, psImageMaskType markVal);
 
 // output functions
@@ -159 +159 @@
 bool            psphotSetState (pmSource *source, bool curState, psImageMaskType maskVal);
 bool            psphotDeblendSatstars (pmReadout *readout, psArray *sources, psMetadata *recipe);
-bool            psphotSourceSize (pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, long first);
+bool            psphotSourceSize (pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, long first);
 
 bool            psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf, psImageMaskType maskVal);
@@ -167 +167 @@
 psKernel       *psphotKernelFromPSF (pmSource *source, int nPix);
 
-bool            psphotRadialProfile (pmSource *source, psMetadata *recipe, psImageMaskType maskVal);
-bool            psphotPetrosian (pmSource *source, psMetadata *recipe, psImageMaskType maskVal);
-bool            psphotIsophotal (pmSource *source, psMetadata *recipe, psImageMaskType maskVal);
-bool            psphotAnnuli (pmSource *source, psMetadata *recipe, psImageMaskType maskVal);
-bool            psphotKron (pmSource *source, psMetadata *recipe, psImageMaskType maskVal);
+// functions related to extended source analysis
+bool  psphotRadialProfile (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal);
+bool  psphotRadialProfilesByAngles (pmSource *source, int Nsec, float Rmax);
+float psphotRadiusFromProfile (pmSource *source, psVector *radius, psVector *flux, float fluxMin, float fluxMax);
+bool  psphotRadiiFromProfiles (pmSource *source, float fluxMin, float fluxMax);
+bool  psphotEllipticalProfile (pmSource *source);
+bool  psphotEllipticalContour (pmSource *source);
 
 // psphotVisual functions
@@ -190 +192 @@
 bool psphotVisualShowSourceSize (pmReadout *readout, psArray *sources);
 bool psphotVisualShowResidualImage (pmReadout *readout);
-bool psphotVisualPlotApResid (psArray *sources);
-bool psphotVisualPlotSourceSize (psArray *sources);
+bool psphotVisualPlotApResid (psArray *sources, float mean, float error);
+bool psphotVisualPlotSourceSize (psMetadata *recipe, psArray *sources);
+bool psphotVisualShowPetrosians (psArray *sources);
+
+// bool psphotPetrosianAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe);
+// bool psphotPetrosianProfile (pmReadout *readout, pmSource *source, float skynoise);
+
+bool psphotPetrosian (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal);
+bool psphotPetrosianRadialBins (pmSource *source, float radiusMax, float skynoise);
+bool psphotPetrosianStats (pmSource *source);
+
+// XXX old versions, currently disabled
+// bool            psphotIsophotal (pmSource *source, psMetadata *recipe, psImageMaskType maskVal);
+// bool            psphotAnnuli (pmSource *source, psMetadata *recipe, psImageMaskType maskVal);
+// bool            psphotKron (pmSource *source, psMetadata *recipe, psImageMaskType maskVal);
+
+// XXX visualization functions related to radial profiles (disabled)
+// bool psphotPetrosianVisualProfileByAngle (psVector *radius, psVector *flux);
+// bool psphotPetrosianVisualProfileRadii (psVector *radius, psVector *flux, psVector *radiusBin, psVector *fluxBin, float peakFlux, float RadiusRef);
+// bool psphotPetrosianVisualEllipticalContour (pmPetrosian *petrosian);
+// bool psphotPetrosianVisualStats (psVector *radBin, psVector *fluxBin, 
+//                               psVector *refRadius, psVector *meanSB, 
+//                               psVector *petRatio, psVector *petRatioErr, psVector *fluxSum, 
+//                               float petRadius, float ratioForRadius,
+//                               float petFlux, float radiusForFlux);
 
 bool psphotImageQuality (psMetadata *recipe, psArray *sources);

trunk/psphot/src/psphotAddNoise.c

@@ -50 +50 @@
 
         // skip sources which were not subtracted
+        // NOTE: this bit is not modified when pmSourceOp applies to noise
         if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue;
 
@@ -79 +80 @@
         axes.minor *= SIZE;
         newshape = psEllipseAxesToShape (axes);
-        PAR[PM_PAR_I0]  = FACTOR*PS_SQR(oldI0);
+        PAR[PM_PAR_I0]  = FACTOR*oldI0;
         PAR[PM_PAR_SXX] = newshape.sx;
         PAR[PM_PAR_SYY] = newshape.sy;

trunk/psphot/src/psphotApResid.c

@@ -33 +33 @@
     if (!measureAptrend) {
         // save nan values since these were not calculated
-        psMetadataAdd (recipe, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   NAN);
-        psMetadataAdd (recipe, PS_LIST_TAIL, "SKYSAT",   PS_DATA_F32 | PS_META_REPLACE, "aperture-determined saturation",   NAN);
         psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   NAN);
         psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", NAN);
+        psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0);
         psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", NAN);
-        psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0);
         return true;
     }
@@ -53 +51 @@
     maskVal |= markVal;
 
-    // S/N limit to perform full non-linear fits
+    // clipping for extreme outliers
+    // XXX this is not currently defined in the recipe
     float MAX_AP_OFFSET = psMetadataLookupF32 (&status, recipe, "MAX_AP_OFFSET");
 
-    // this is the smallest radius allowed: need to at least extend growth curve down to this...
+    // options for how the photometry is calculated
+    // XXX are these sensible?
     bool IGNORE_GROWTH = psMetadataLookupBool (&status, recipe, "IGNORE_GROWTH");
     bool INTERPOLATE_AP = psMetadataLookupBool (&status, recipe, "INTERPOLATE_AP");
@@ -100 +100 @@
             PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32);
             PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
+            PS_ARRAY_ADD_SCALAR(job->args, markVal,  PS_TYPE_IMAGE_MASK);
 
             PS_ARRAY_ADD_SCALAR(job->args, 0,        PS_TYPE_S32); // this is used as a return value for Nskip
@@ -110 +111 @@
             }
             psFree(job);
-
-# if (0)
-                int nskip = 0;
-                int nfail = 0;
-
-                if (!psphotApResidMags_Unthreaded (&nskip, &nfail, sources, psf, photMode, maskVal)) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
-                    return false;
-                }
-                Nskip += nskip;
-                Nfail += nfail;
-# endif
-
         }
 
@@ -138 +126 @@
             } else {
                 psScalar *scalar = NULL;
-                scalar = job->args->data[4];
+                scalar = job->args->data[5];
                 Nskip += scalar->data.S32;
-                scalar = job->args->data[5];
+                scalar = job->args->data[6];
                 Nfail += scalar->data.S32;
             }
@@ -150 +138 @@
 
     // gather the stats to assess the aperture residuals
-    psVector *mask    = psVectorAllocEmpty (300, PS_TYPE_VECTOR_MASK);
     psVector *mag     = psVectorAllocEmpty (300, PS_TYPE_F32);
     psVector *xPos    = psVectorAllocEmpty (300, PS_TYPE_F32);
@@ -158 +145 @@
     Npsf = 0;
 
+# ifdef DEBUG    
+    FILE *f = fopen ("apresid.dat", "w");
+    psAssert (f, "failed open");
+# endif
+
     for (int i = 0; i < sources->n; i++) {
         source = sources->data[i];
@@ -170 +162 @@
         if (source->mode &  PM_SOURCE_MODE_EXT_LIMIT) SKIPSTAR ("EXTENDED");
         if (source->mode &  PM_SOURCE_MODE_CR_LIMIT) SKIPSTAR ("COSMIC RAY");
+        if (source->mode &  PM_SOURCE_MODE_DEFECT) SKIPSTAR ("DEFECT");
             
         if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
             continue;
         }
+
+        // XXX make this user-configurable?
+        if (source->errMag > 0.01) continue;
 
         // aperture residual for this source
@@ -185 +181 @@
         }
 
-        mag->data.F32[Npsf]     = source->psfMag;
-        apResid->data.F32[Npsf] = dap;
-        xPos->data.F32[Npsf]    = model->params->data.F32[PM_PAR_XPOS];
-        yPos->data.F32[Npsf]    = model->params->data.F32[PM_PAR_YPOS];
-
-        mask->data.PS_TYPE_VECTOR_MASK_DATA[Npsf] = 0;
-
-        dMag->data.F32[Npsf] = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
-
-        psVectorExtend (mag,     100, 1);
-        psVectorExtend (mask,    100, 1);
-        psVectorExtend (xPos,    100, 1);
-        psVectorExtend (yPos,    100, 1);
-        psVectorExtend (dMag,    100, 1);
-        psVectorExtend (apResid, 100, 1);
+# ifdef DEBUG
+        fprintf (f, "%6.1f %6.1f : %6.1f %6.1f : %8.3f %8.3f %8.3f : %f : %f %f %f : %f\n",
+                 source->peak->xf, source->peak->yf, 
+                 source->modelPSF->params->data.F32[PM_PAR_XPOS], source->modelPSF->params->data.F32[PM_PAR_YPOS], 
+                 source->psfMag, source->apMag, source->errMag,
+                 source->modelPSF->params->data.F32[PM_PAR_I0], 
+                 source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SXY], source->modelPSF->params->data.F32[PM_PAR_SYY], 
+                 source->modelPSF->params->data.F32[PM_PAR_7]);
+# endif
+        if (!isfinite(source->psfMag)) psAbort ("nan in psfMag");
+        if (!isfinite(source->errMag)) psAbort ("nan in errMag");
+        if (!isfinite(source->apMag)) psAbort ("nan in apMag");
+        if (!isfinite(model->params->data.F32[PM_PAR_XPOS])) psAbort ("nan in xPos");
+        if (!isfinite(model->params->data.F32[PM_PAR_YPOS])) psAbort ("nan in yPos");
+
+        psVectorAppend (mag, source->psfMag);
+        psVectorAppend (dMag,source->errMag);
+        psVectorAppend (apResid, dap);
+        psVectorAppend (xPos, model->params->data.F32[PM_PAR_XPOS]);
+        psVectorAppend (yPos, model->params->data.F32[PM_PAR_YPOS]);
         Npsf ++;
     }
@@ -205 +206 @@
     psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "measure aperture residuals for %d objects (%d skipped, %d failed, %ld invalid)\n",
               Npsf, Nskip, Nfail, sources->n - Npsf - Nskip - Nfail);
+
+# ifdef DEBUG
+    fclose (f);
+# endif
 
     // XXX choose a better value here?
@@ -212 +217 @@
     }
 
-    // XXX deprecating the old code which allowed the ApResid to be fitted as a function of flux and r^2/flux
-    // XXX is this asymmetric clipping still needed?  this analysis should come after neighbors are subtracted...
-    // 3hi/1lo sigma clipping on the rflux vs metric fit
-    // systematic error information
-    float errorScale = 0.0;
-    float errorFloor = 0.0;
-
+    // XXX set the min number of needed source more carefully
+    if ((Npsf < 15) && (APTREND_ORDER_MAX >= 4)) APTREND_ORDER_MAX = 3;
+    if ((Npsf < 11) && (APTREND_ORDER_MAX >= 3)) APTREND_ORDER_MAX = 2;
+    if ((Npsf <  8) && (APTREND_ORDER_MAX >= 2)) APTREND_ORDER_MAX = 1;
+
+    psFree (psf->ApTrend);
+    psf->ApTrend = NULL;
     float errorFloorMin = FLT_MAX;
-    int entryMin = -1;
-
-    // Fit out the dap vs mag trend, iterate over spatial scale until error Floor increases.
-    // Stop if Npsf / (Nx * Ny) < 3
+
+    // as we loop over orders, we need to refer to the initial selection, but we modify the
+    // option values to match the current guess: save the max values here:
+    int NX = readout->image->numCols;
+    int NY = readout->image->numRows;
     for (int i = 1; i <= APTREND_ORDER_MAX; i++) {
-
-        if (!psphotApResidTrend (readout, psf, Npsf, i, &errorScale, &errorFloor, mask, xPos, yPos, apResid, dMag)) {
-            break;
-        }
-
-        // store the resulting errorFloor values and the scales, redo the best
+        
+        int Nx, Ny;
+        if (NX > NY) {
+            Nx = i;
+            Ny = PS_MAX (1, (int)(i * (NY / (float)(NX)) + 0.5));
+        } else {
+            Ny = i;
+            Nx = PS_MAX (1, (int)(i * (NX / (float)(NY)) + 0.5));
+        }
+
+        float errorFloor;
+        pmTrend2D *apTrend = psphotApResidTrend (&errorFloor, readout, Nx, Ny, xPos, yPos, apResid, dMag);
+        if (!apTrend) {
+            continue;
+        }
+
+        // apply ApTrend results
+        // float xc = 0.5*readout->image->numCols + readout->image->col0 + 0.5;
+        // float yc = 0.5*readout->image->numRows + readout->image->row0 + 0.5;
+        // float ApResid = pmTrend2DEval (psf->ApTrend, xc, yc); // ap-fit at chip center
+        // if (!isfinite(ApResid)) psAbort("nan apresid @ center");
+
+        // store the minimum errorFloor and best ApTrend to keep
         if (errorFloor < errorFloorMin) {
             errorFloorMin = errorFloor;
-            entryMin = i;
+            psFree (psf->ApTrend);
+            psf->ApTrend = psMemIncrRefCounter(apTrend);
         }
-    }
-    if (entryMin == -1) {
+        psFree (apTrend);
+    }
+    if (psf->ApTrend == NULL) {
         psWarning("Failed to find a valid aperture residual value");
         goto escape;
     }
 
-    // XXX catch error condition
-    psphotApResidTrend (readout, psf, Npsf, entryMin, &errorScale, &errorFloor, mask, xPos, yPos, apResid, dMag);
+    // apply ApTrend results
+    float xc = 0.5*readout->image->numCols + readout->image->col0 + 0.5;
+    float yc = 0.5*readout->image->numRows + readout->image->row0 + 0.5;
+
+    psf->ApResid  = pmTrend2DEval (psf->ApTrend, xc, yc); // ap-fit at chip center
+    psf->dApResid = errorFloorMin;
+    psf->nApResid = Npsf;
+
+    // save results for later output
+    psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   psf->ApResid);
+    psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid);
+    psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", psf->nApResid);
+    psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", psf->growth->apLoss);
+
+    psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "aperture residual: %f +/- %f\n", psf->ApResid, psf->dApResid);
+    psLogMsg ("psphot.apresid", PS_LOG_INFO, "measure full-frame aperture residuals for %d sources: %f sec\n", Npsf, psTimerMark ("psphot.apresid"));
+
+    psFree (xPos);
+    psFree (yPos);
+    psFree (apResid);
+    psFree (mag);
+    psFree (dMag);
+
+    psphotVisualPlotApResid (sources, psf->ApResid, psf->dApResid);
+
+    return true;
+
+escape:
+    // save nan values since these were not calculated
+    psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   NAN);
+    psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", NAN);
+    psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0);
+    psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", NAN);
+
+    psFree (xPos);
+    psFree (yPos);
+    psFree (apResid);
+    psFree (mag);
+    psFree (dMag);
+    return false;
+}
+
+pmTrend2D *psphotApResidTrend (float *apResidSysErr, pmReadout *readout, int Nx, int Ny, psVector *xPos, psVector *yPos, psVector *apResid, psVector *dMag) {
+
+    // the mask marks the values not used to calculate the ApTrend
+    psVector *mask = psVectorAlloc(xPos->n, PS_TYPE_VECTOR_MASK);
+    psVectorInit (mask, 0);
+
+    // XXX allow user to set this optionally?
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+
+    // measure Trend2D for the current spatial scale
+    pmTrend2D *apTrend = pmTrend2DAlloc (PM_TREND_MAP, readout->image, Nx, Ny, stats);
+
+    // XXX somewhat arbitrary: soften the errors so the few bright stars do not totally dominate:
+    // XXX use this or not?  probably not, since this is the point of the systematic error analysis
+    psVector *dMagSoft = psVectorAlloc (dMag->n, PS_TYPE_F32);
+    for (int i = 0; i < dMag->n; i++) {
+        dMagSoft->data.F32[i] = hypot(dMag->data.F32[i], 0.005);
+    }
+
+    // XXX test for errors here
+    if (!pmTrend2DFit (apTrend, mask, 0xff, xPos, yPos, apResid, dMagSoft)) {
+        psWarning("Failed to fit trend for %d x %d map", Nx, Ny);
+        psFree (apTrend);
+        return NULL;
+    }
+    if (apTrend->mode == PM_TREND_MAP) {
+        // p_psImagePrint (2, apTrend->map->map, "ApTrend Before"); // XXX TEST:
+        psImageMapRepair (apTrend->map->map);
+        // p_psImagePrint (2, apTrend->map->map, "ApTrend After"); // XXX TEST:
+    }
 
     // construct the fitted values and the residuals
-    psVector *apResidFit = pmTrend2DEvalVector (psf->ApTrend, mask, 0xff, xPos, yPos);
+    psVector *apResidFit = pmTrend2DEvalVector (apTrend, mask, 0xff, xPos, yPos);
     psVector *apResidRes = (psVector *) psBinaryOp (NULL, (void *) apResid, "-", (void *) apResidFit);
-    psVector *dMagSys = (psVector *) psBinaryOp (NULL, (void *) dMag, "*", (void *) psScalarAlloc(errorScale, PS_TYPE_F32));
-
-    if (psTraceGetLevel("psphot") >= 2) {
-        FILE *dumpFile = fopen ("apresid.dat", "w");
+
+    // measure systematic error
+    *apResidSysErr = psVectorSystematicError (apResidRes, dMag, 0.10);
+    if (!isfinite(*apResidSysErr)) {
+        psWarning("Failed to find systematic error for %d x %d map", Nx, Ny);
+        psFree (apTrend);
+        return NULL;
+    }
+
+    psLogMsg ("psphot.apresid", PS_LOG_INFO, "result of %d x %d grid\n", Nx, Ny);
+    psLogMsg ("psphot.apresid", PS_LOG_INFO, "systematic scatter floor: %f\n", *apResidSysErr);
+
+    if (psTraceGetLevel("psphot") >= 4) {
+        char filename[64];
+        snprintf (filename, 64, "apresid.%dx%d.dat", Nx, Ny);
+        FILE *dumpFile = fopen (filename, "w");
         for (int i = 0; i < xPos->n; i++) {
-            fprintf (dumpFile, "%f %f  %f %f %f  %f %f  %x\n",
+            fprintf (dumpFile, "%f %f  %f %f  %f %f  %x\n",
                      xPos->data.F32[i], yPos->data.F32[i],
-                     mag->data.F32[i], dMag->data.F32[i], dMagSys->data.F32[i],
+                     dMag->data.F32[i], hypot(dMag->data.F32[i], *apResidSysErr),
                      apResid->data.F32[i], apResidRes->data.F32[i],
                      mask->data.PS_TYPE_VECTOR_MASK_DATA[i]);
@@ -261 +368 @@
     }
 
-    // apply ApTrend results
-    float xc = 0.5*readout->image->numCols + readout->image->col0 + 0.5;
-    float yc = 0.5*readout->image->numRows + readout->image->row0 + 0.5;
-
-    psf->ApResid  = pmTrend2DEval (psf->ApTrend, xc, yc); // ap-fit at chip center
-    psf->dApResid = errorFloor;
-    psf->nApResid = Npsf;
-
-    // save results for later output
-    psMetadataAdd (recipe, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   0.0);
-    psMetadataAdd (recipe, PS_LIST_TAIL, "SKYSAT",   PS_DATA_F32 | PS_META_REPLACE, "aperture-determined saturation",   0.0);
-    psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   psf->ApResid);
-    psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid);
-    psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", psf->growth->apLoss);
-    psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", psf->nApResid);
-
-    psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "aperture residual: %f +/- %f\n", psf->ApResid, psf->dApResid);
-    psLogMsg ("psphot.apresid", PS_LOG_INFO, "measure full-frame aperture residuals for %d sources: %f sec\n", Npsf, psTimerMark ("psphot.apresid"));
-
-    psFree (mag);
     psFree (mask);
-    psFree (xPos);
-    psFree (yPos);
-
-    psFree (apResid);
-    psFree (apResidFit);
-    psFree (apResidRes);
-
-    psFree (dMagSys);
-    psFree (dMag);
-
-    psphotVisualPlotApResid (sources);
-
-    return true;
-
-escape:
-    // save nan values since these were not calculated
-    psMetadataAdd (recipe, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   NAN);
-    psMetadataAdd (recipe, PS_LIST_TAIL, "SKYSAT",   PS_DATA_F32 | PS_META_REPLACE, "aperture-determined saturation",   NAN);
-    psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   NAN);
-    psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", NAN);
-    psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", NAN);
-    psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0);
-
-    psFree (mag);
-    psFree (mask);
-    psFree (xPos);
-    psFree (yPos);
-    psFree (apResid);
-    psFree (dMag);
-    return false;
-}
-
-/*
-  (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y)
-  (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y)
-  (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y)
-*/
-
-// XXX this still sucks...  need a better way to estimate the error floor...
-bool psphotMagErrorScale (float *errorScale, float *errorFloor, psVector *dMag, psVector *dap, psVector *mask, int nGroup) {
-
-    psStats *statsS = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
-    psStats *statsM = psStatsAlloc (PS_STAT_SAMPLE_MEAN);
-
-    // measure the trend in bins with 10 values each; if < 10 total, use them all
-    int nBin = PS_MAX (dMag->n / nGroup, 1);
-
-    // output vectors for ApResid trend
-    psVector *dSo = psVectorAlloc (nBin, PS_TYPE_F32);
-    psVector *dMo = psVectorAlloc (nBin, PS_TYPE_F32);
-    psVector *dRo = psVectorAlloc (nBin, PS_TYPE_F32);
-
-    // use dMag to group the dMag and dap vectors
-    psVector *index = psVectorSortIndex (NULL, dMag);
-
-    // subset vectors for dMag and dap values within the given range
-    psVector *dMSubset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
-    psVector *dASubset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
-    psVector *mkSubset = psVectorAllocEmpty (nGroup, PS_TYPE_VECTOR_MASK);
-
-    int n = 0;
-    for (int i = 0; i < dMo->n; i++) {
-        int j;
-        for (j = 0; (j < nGroup) && (n < dMag->n); j++, n++) {
-            int N = index->data.U32[n];
-            dMSubset->data.F32[j] = dMag->data.F32[N];
-            dASubset->data.F32[j] = dap->data.F32[N];
-            mkSubset->data.PS_TYPE_VECTOR_MASK_DATA[j] = mask->data.PS_TYPE_VECTOR_MASK_DATA[N];
-        }
-        dMSubset->n = j;
-        dASubset->n = j;
-        mkSubset->n = j;
-
-        psStatsInit (statsS);
-        psStatsInit (statsM);
-
-        if (j > 2) {
-            if (!psVectorStats (statsS, dASubset, NULL, mkSubset, 0xff)) {
-                psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
-                return false;
-            }
-            if (!psVectorStats (statsM, dMSubset, NULL, mkSubset, 0xff)) {
-                psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
-                return false;
-            }
-            dSo->data.F32[i] = statsS->robustStdev;
-            dMo->data.F32[i] = statsM->sampleMean;
-            dRo->data.F32[i] = statsS->robustStdev / statsM->sampleMean;
-        } else {
-            dSo->data.F32[i] = NAN;
-            dMo->data.F32[i] = NAN;
-            dRo->data.F32[i] = NAN;
-        }
-    }
-    psFree (dMSubset);
-    psFree (dASubset);
-    psFree (mkSubset);
-
-    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
-    if (!psVectorStats (stats, dRo, NULL, NULL, 0)) {
-        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
-        return false;
-    }
-
-    *errorScale = stats->sampleMedian;
-    for (int i = 0; i < dSo->n; i++) {
-        *errorFloor = dSo->data.F32[i];
-        if (fabs(*errorFloor) <= FLT_EPSILON) continue;
-        if (isfinite(*errorFloor)) break;
-    }
-
-    psFree (stats);
-    psFree (index);
-
-    psFree (dRo);
-    psFree (dMo);
-    psFree (dSo);
-
-    psFree (statsS);
-    psFree (statsM);
-
-    return true;
-}
-
-bool psphotApResidTrend (pmReadout *readout, pmPSF *psf, int Npsf, int scale, float *errorScale, float *errorFloor, psVector *mask, psVector *xPos, psVector *yPos, psVector *apResid, psVector *dMag) {
-
-    int Nx, Ny;
-
-    if (readout->image->numCols > readout->image->numRows) {
-        Nx = scale;
-        float AR = readout->image->numRows / (float) readout->image->numCols;
-        Ny = (int) (Nx * AR + 0.5);
-        Ny = PS_MAX (1, Ny);
-    } else {
-        Ny = scale;
-        float AR = readout->image->numRows / (float) readout->image->numCols;
-        Nx = (int) (Ny * AR + 0.5);
-        Nx = PS_MAX (1, Nx);
-    }
-
-    // require at least 10 stars per spatial bin
-    if (Npsf < 10*Nx*Ny) {
-        return false;
-    }
-
-    // the mask marks the values not used to calculate the ApTrend
-    psVectorInit (mask, 0);
-
-    // XXX stats structure for use by ApTrend : make parameters user setable
-    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
-    stats->min = 2.0;
-    stats->max = 3.0;
-
-    // measure Trend2D for the current spatial scale
-    psFree (psf->ApTrend);
-    psf->ApTrend = pmTrend2DAlloc (PM_TREND_MAP, readout->image, Nx, Ny, stats);
-
-    // XXX somewhat arbitrary: soften the errors so the few bright stars do not totally dominate:
-    psVector *dMagSoft = psVectorAlloc (dMag->n, PS_TYPE_F32);
-    for (int i = 0; i < dMag->n; i++) {
-        dMagSoft->data.F32[i] = hypot(dMag->data.F32[i], 0.01);
-    }
-
-    // XXX test for errors here
-    pmTrend2DFit (psf->ApTrend, mask, 0xff, xPos, yPos, apResid, dMagSoft);
-
-    // construct the fitted values and the residuals
-    psVector *apResidFit = pmTrend2DEvalVector (psf->ApTrend, mask, 0xff, xPos, yPos);
-    psVector *apResidRes = (psVector *) psBinaryOp (NULL, (void *) apResid, "-", (void *) apResidFit);
-
-    // measure systematic errorFloor & systematic / photon scale factor
-    // XXX this is a bit arbitrary, but it forces ~3 stars from the bright bin per spatial bin
-    int nGroup = PS_MAX (3*Nx*Ny, 10);
-    psphotMagErrorScale (errorScale, errorFloor, dMag, apResidRes, mask, nGroup);
-
-    psLogMsg ("psphot.apresid", PS_LOG_INFO, "result of %d x %d grid (%d stars per bin)\n", Nx, Ny, nGroup);
-    psLogMsg ("psphot.apresid", PS_LOG_INFO, "systematic error / photon error: %f\n", *errorScale);
-    psLogMsg ("psphot.apresid", PS_LOG_INFO, "systematic scatter floor: %f\n", *errorFloor);
-
     psFree (stats);
     psFree (dMagSoft);
@@ -466 +374 @@
     psFree (apResidRes);
 
-    return true;
+    return apTrend;
 }
 
@@ -480 +388 @@
     pmSourcePhotometryMode photMode = PS_SCALAR_VALUE(job->args->data[2],S32);
     psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[3],PS_TYPE_IMAGE_MASK_DATA);
+    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
 
     for (int i = 0; i < sources->n; i++) {
@@ -490 +399 @@
         if (source->mode &  PM_SOURCE_MODE_POOR) SKIPSTAR ("POOR STAR");
 
-        if (!pmSourceMagnitudes (source, psf, photMode, maskVal)) {
-            Nskip ++;
-            psTrace ("psphot", 3, "skip : bad source mag");
-            continue;
+        // replace object in image
+        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
+            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
         }
 
-        if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
-            Nfail ++;
-            psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag);
-            continue;
-        }
-        source->mode |= PM_SOURCE_MODE_AP_MAGS;
+        // clear the mask bit and set the circular mask pixels
+        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
+        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal);
+
+        bool status = pmSourceMagnitudes (source, psf, photMode, maskVal);
+
+        // clear the mask bit 
+        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
+
+        // re-subtract the object, leave local sky
+        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+
+        if (!status) {
+            Nskip ++;
+            psTrace ("psphot", 3, "skip : bad source mag");
+            continue;
+        }
+    
+        if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
+            Nfail ++;
+            psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag);
+            continue;
+        }
+        source->mode |= PM_SOURCE_MODE_AP_MAGS;
     }
 
     // change the value of a scalar on the array (wrap this and put it in psArray.h)
-    scalar = job->args->data[4];
+    scalar = job->args->data[5];
     scalar->data.S32 = Nskip;
 
-    scalar = job->args->data[5];
+    scalar = job->args->data[6];
     scalar->data.S32 = Nfail;
 
     return true;
 }
-
-# if (0)
-bool psphotApResidMags_Unthreaded (int *nskip, int *nfail, psArray *sources, pmPSF *psf, pmSourcePhotometryMode photMode, psImageMaskType maskVal) {
-
-    int Nskip = 0;
-    int Nfail = 0;
-
-    for (int i = 0; i < sources->n; i++) {
-        pmSource *source = (pmSource *) sources->data[i];
-
-        if (source->type != PM_SOURCE_TYPE_STAR) SKIPSTAR ("NOT STAR");
-        if (source->mode &  PM_SOURCE_MODE_SATSTAR) SKIPSTAR ("SATSTAR");
-        if (source->mode &  PM_SOURCE_MODE_BLEND) SKIPSTAR ("BLEND");
-        if (source->mode &  PM_SOURCE_MODE_FAIL) SKIPSTAR ("FAIL STAR");
-        if (source->mode &  PM_SOURCE_MODE_POOR) SKIPSTAR ("POOR STAR");
-
-        if (!pmSourceMagnitudes (source, psf, photMode, maskVal)) {
-            Nskip ++;
-            psTrace ("psphot", 3, "skip : bad source mag");
-            continue;
-        }
-
-        if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
-            Nfail ++;
-            psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag);
-            continue;
-        }
-    }
-
-    // change the value of a scalar on the array (wrap this and put it in psArray.h)
-    *nskip = Nskip;
-    *nfail = Nfail;
-
-    return true;
-}
-# endif

trunk/psphot/src/psphotBlendFit.c

@@ -254 +254 @@
         pmSourceCacheModel (source, maskVal);
         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-        source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
     }
 
@@ -365 +364 @@
         pmSourceCacheModel (source, maskVal);
         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-        source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
     }
 
@@ -374 +372 @@
     *nfail = Nfail;
 
+    // moments are modified by the fit; re-display
+    psphotVisualPlotMoments (recipe, sources);
+    psphotVisualShowResidualImage (readout);
+
     return true;
 }

trunk/psphot/src/psphotChoosePSF.c

@@ -1 +1 @@
 # include "psphotInternal.h"
-
-void psphotCountPSFStars (psArray *sources) {
-
-    int nPSF = 0;
-
-    // select the candidate PSF stars (pointers to original sources)
-    for (int i = 0; i < sources->n; i++) {
-        pmSource *source = sources->data[i];
-        if (source->mode & PM_SOURCE_MODE_PSFSTAR) {
-            nPSF ++;
-        }
-    }
-
-    fprintf (stderr, "N PSF: %d\n", nPSF);
-}
 
 // try PSF models and select best option
@@ -73 +58 @@
     // get the fixed PSF fit radius
     // XXX check that PSF_FIT_RADIUS < SKY_OUTER_RADIUS
-    options->radius = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS");
-    assert (status);
+    // options->radius = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS");
+    // assert (status);
+
+    // We have calculated a Gaussian window function, use that for both the PSF fit radius and
+    // the aperture radius (scaling SIGMA)
+    float gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA");
+    float fitScale = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS_SCALE");
+    float apScale = psMetadataLookupF32(&status, recipe, "PSF_APERTURE_SCALE");
+    options->fitRadius = (int)(fitScale*gaussSigma);
+    options->apRadius = (int)(apScale*gaussSigma);
+
+    // XXX use the same radii for standard analysis as for the PSF creation
+    psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_FIT_RADIUS", PS_META_REPLACE, "fit radius", options->fitRadius);
+    psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_APERTURE", PS_META_REPLACE, "psf aperture", options->apRadius);
 
     // XXX ROBUST seems to be too agressive given the small numbers.
@@ -99 +96 @@
 
     psArray *stars = psArrayAllocEmpty (sources->n);
-
-    // psphotCountPSFStars (sources);
 
     // select the candidate PSF stars (pointers to original sources)
@@ -115 +110 @@
     }
 
-    // psphotCountPSFStars (sources);
-
     // check that the identified psf stars sufficiently cover the region; if not, extend the
     // limits somewhat
     psphotCheckStarDistribution (stars, sources, options);
-
-    // psphotCountPSFStars (sources);
 
     psLogMsg ("psphot.pspsf", PS_LOG_DETAIL, "selected candidate %ld PSF objects\n", stars->n);
@@ -289 +280 @@
     // XXX test dump of psf star data and psf-subtracted image
     if (psTraceGetLevel("psphot.psfstars") > 5) {
-        psphotDumpPSFStars (readout, try, options->radius, maskVal, markVal);
+        psphotDumpPSFStars (readout, try, options->fitRadius, maskVal, markVal);
     }
 
     // save only the best model;
-    // XXX we are not saving the fitted sources
-    // XXX do we want to keep them so we may optionally write them out?
     pmPSF *psf = psMemIncrRefCounter(try->psf);
     psFree (models);
@@ -305 +294 @@
     }
 
-    // psphotCountPSFStars (sources);
-
     char *modelName = pmModelClassGetName (psf->type);
     psLogMsg ("psphot.pspsf", PS_LOG_INFO, "select psf model: %f sec\n", psTimerMark ("psphot.choose.psf"));
@@ -341 +328 @@
             // create modelPSF from this model
             pmModel *modelPSF = pmModelFromPSFforXY (psf, xc, yc, 1.0);
-            if (!modelPSF) continue;
+            if (!modelPSF) {
+                fprintf (stderr, "?");
+                continue;
+            }
 
             // get the model full-width at half-max
@@ -354 +344 @@
 
             float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
-            if (!isfinite(FWHM_MAJOR) || !isfinite(FWHM_MINOR)) continue;
+            if (!isfinite(FWHM_MAJOR) || !isfinite(FWHM_MINOR)) {
+                fprintf (stderr, "!");
+                continue;
+            }
             psVectorAppend (fwhmMajor, FWHM_MAJOR);
             psVectorAppend (fwhmMinor, FWHM_MINOR);

trunk/psphot/src/psphotExtendedSourceAnalysis.c

@@ -21 +21 @@
     }
 
-    // option to limit analysis to a specific region
-    char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
-    psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
-    if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
+    // XXX require petrosian analysis for non-linear fits? 
+
+    // XXX temporary user-supplied systematic sky noise measurement (derive from background model)
+    float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE");
+
+# if (0)
+    // if backModel or backStdev are missing, the values of sky and/or skyErr will be set to NAN
+    // XXX use this to set skynoise
+    pmReadout *backModel = psphotSelectBackground (config, view);
+    pmReadout *backStdev = psphotSelectBackgroundStdev (config, view);
+# endif
 
     // S/N limit to perform full non-linear fits
@@ -38 +45 @@
     sources = psArraySort (sources, pmSourceSortBySN);
 
-    // XXX some init functions for the extended source recipe options?
+    // option to limit analysis to a specific region
+    char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
+    psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
+    if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
 
     // choose the sources of interest
@@ -49 +59 @@
         if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
         if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
+        if (source->mode & PM_SOURCE_MODE_DEFECT) continue;
+        if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
+        if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;
 
         // limit selection to some SN limit
@@ -68 +81 @@
         // if we request any of these measurements, we require the radial profile
         if (doPetrosian || doIsophotal || doAnnuli || doKron) {
-            if (!psphotRadialProfile (source, recipe, maskVal)) {
+            if (!psphotRadialProfile (source, recipe, skynoise, maskVal)) {
                 // all measurements below require the radial profile; skip them all
                 // re-subtract the object, leave local sky
                 psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
                 pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-                source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
                 continue;
             }
@@ -79 +91 @@
         }
 
+        // Petrosian Mags
+        if (doPetrosian) {
+            if (!psphotPetrosian (source, recipe, skynoise, maskVal)) {
+                psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
+            } else {
+                psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
+                Npetro ++;
+                source->mode |= PM_SOURCE_MODE_EXTENDED_STATS;
+            }
+        }
+
+# if (0)
         // Isophotal Mags
         if (doIsophotal) {
@@ -89 +113 @@
             }
         }
-
-        // Petrosian Mags
-        if (doPetrosian) {
-            if (!psphotPetrosian (source, recipe, maskVal)) {
-                psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
-            } else {
-                psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
-                Npetro ++;
-                source->mode |= PM_SOURCE_MODE_EXTENDED_STATS;
-            }
-        }
-
         // Kron Mags
         if (doKron) {
@@ -111 +123 @@
             }
         }
-
-        // Radial Annuli
-        if (doAnnuli) {
-            if (!psphotAnnuli (source, recipe, maskVal)) {
-                psError(PSPHOT_ERR_UNKNOWN, false, "failure in Annuli analysis");
-                return false;
-            } 
-            psTrace ("psphot", 5, "measured annuli for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
-            Nannuli ++;
-            source->mode |= PM_SOURCE_MODE_EXTENDED_STATS;
-        }
+# endif
 
         // re-subtract the object, leave local sky
         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-        source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
+
+        if (source->extpars) {
+            pmSourceRadialProfileFreeVectors(source->extpars->profile);
+        }
     }
 
@@ -133 +138 @@
     psLogMsg ("psphot", PS_LOG_INFO, "  %d annuli\n", Nannuli);
     psLogMsg ("psphot", PS_LOG_INFO, "  %d kron\n", Nkron);
+
+    psphotVisualShowResidualImage (readout);
+
+    if (doPetrosian) {
+        psphotVisualShowPetrosians (sources);
+    }
+
     return true;
 }

trunk/psphot/src/psphotExtendedSourceFits.c

@@ -226 +226 @@
 
           pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-          source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
 
           psFree (modelFluxes);
@@ -253 +252 @@
         // subtract the best fit from the object, leave local sky
         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-        source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
 
         // the initial model flux is no longer needed

trunk/psphot/src/psphotFindDetections.c

@@ -52 +52 @@
     // optionally merge peaks into footprints
     if (useFootprints) {
-        psphotFindFootprints (detections, significance, readout, recipe, pass, maskVal);
+        psphotFindFootprints (detections, significance, readout, recipe, threshold, pass, maskVal);
     }
 

trunk/psphot/src/psphotFindFootprints.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
-bool psphotFindFootprints (pmDetections *detections, psImage *significance, pmReadout *readout, psMetadata *recipe, const int pass, psImageMaskType maskVal) {
+bool psphotFindFootprints (pmDetections *detections, psImage *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int pass, psImageMaskType maskVal) {
 
     bool status;
@@ -9 +9 @@
     int npixMin = psMetadataLookupS32(&status, recipe, "FOOTPRINT_NPIXMIN");
     PS_ASSERT (status, false);
-
-    float FOOTPRINT_NSIGMA_LIMIT;
-    if (pass == 1) {
-        FOOTPRINT_NSIGMA_LIMIT = psMetadataLookupS32(&status, recipe, "FOOTPRINT_NSIGMA_LIMIT");
-    } else {
-        FOOTPRINT_NSIGMA_LIMIT = psMetadataLookupS32(&status, recipe, "FOOTPRINT_NSIGMA_LIMIT_2");
-    }
-    PS_ASSERT (status, false);
-
-    // XXX do we need to use the same threshold here as for peaks?  does it make sense for
-    // these to be different?
-
-    float threshold = PS_SQR(FOOTPRINT_NSIGMA_LIMIT);
 
     int growRadius = 0;

trunk/psphot/src/psphotFindPeaks.c

@@ -31 +31 @@
         peak->SN = sqrt(peak->value);
         peak->flux = readout->image->data.F32[peak->y-row0][peak->x-col0];
+        // if (peak->flux / peak->value > 5.0/12.0) {
+        //     psWarning ("odd peak levels (1)");
+        // }
+        // if (peak->value / peak->flux > 5*12.0) {
+        //     psWarning ("odd peak levels (2)");
+        // }
         if (readout->variance && isfinite (peak->dx)) {
             peak->dx *= sqrt(readout->variance->data.F32[peak->y-row0][peak->x-col0]);

trunk/psphot/src/psphotFitSourcesLinear.c

@@ -76 +76 @@
             if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) continue;
         } else {
-            if (source->mode & PM_SOURCE_MODE_BLEND) continue;
+            // if (source->mode & PM_SOURCE_MODE_BLEND) continue;
         }
 
@@ -186 +186 @@
     if (SKY_FIT_LINEAR) {
         psSparseBorderSolve (&norm, &skyfit, constraint, border, 5);
-        fprintf (stderr, "skyfit: %f\n", skyfit->data.F32[0]);
+        psLogMsg ("psphot", PS_LOG_MINUTIA, "skyfit: %f\n", skyfit->data.F32[0]);
     } else {
         norm = psSparseSolve (NULL, constraint, sparse, 5);
@@ -215 +215 @@
         // subtract object
         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-        source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
     }
     psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "sub models: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
@@ -239 +238 @@
 
     psphotVisualShowResidualImage (readout);
-    psphotVisualShowFlags (sources);
+    // psphotVisualShowFlags (sources);
 
     return true;
@@ -264 +263 @@
         float x = model->params->data.F32[PM_PAR_XPOS];
         float y = model->params->data.F32[PM_PAR_YPOS];
-        psImageMaskCircle (source->maskView, x, y, model->radiusFit, "AND", PS_NOT_IMAGE_MASK(markVal));
+        psImageMaskCircle (source->maskView, x, y, model->fitRadius, "AND", PS_NOT_IMAGE_MASK(markVal));
     }
 

trunk/psphot/src/psphotFitSourcesLinearStack.c

@@ -168 +168 @@
         // subtract object
         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-        source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
     }
     psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "sub models: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);

trunk/psphot/src/psphotGuessModels.c

@@ -177 +177 @@
 
         // set the source PSF model
+        psAssert (source->modelPSF == NULL, "failed to free one of the models?");
         source->modelPSF = modelPSF;
         source->modelPSF->residuals = psf->residuals;
 
         pmSourceCacheModel (source, maskVal);  // ALLOC x14 (!)
-
     }
 

trunk/psphot/src/psphotImageLoop.c

@@ -67 +67 @@
 
                 // Update the header
-                {
-                    pmHDU *hdu = pmHDUGetHighest(input->fpa, chip, cell);
-                    if (hdu && hdu != lastHDU) {
-                        psphotVersionHeaderFull(hdu->header);
-                        lastHDU = hdu;
-                    }
+                pmHDU *hdu = pmHDUGetHighest(input->fpa, chip, cell);
+                if (hdu && hdu != lastHDU) {
+                    psphotVersionHeaderFull(hdu->header);
+                    lastHDU = hdu;
                 }
 
-                psImageMaskType maskSat = pmConfigMaskGet("SAT", config); // Mask value for saturated pixels
-                if (!pmReadoutMaskNonfinite(readout, maskSat)) {
-                    psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels.");
-                    psFree(view);
-                    return false;
-                }
+                // if an external mask is supplied, ensure that NAN pixels are also masked
+                if (readout->mask) {
+                    psImageMaskType maskSat = pmConfigMaskGet("SAT", config); // Mask value for saturated pixels
+                    if (!pmReadoutMaskNonfinite(readout, maskSat)) {
+                        psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels.");
+                        psFree(view);
+                        return false;
+                    }
+                }
 
                 // run the actual photometry analysis on this chip/cell/readout

trunk/psphot/src/psphotImageQuality.c

@@ -279 +279 @@
 #endif
 
-    psLogMsg ("psphot", PS_LOG_INFO, "Image Quality Stats from %ld psf stars : FWHM (major, minor) : %f, %f\n",
+    psLogMsg ("psphot", PS_LOG_INFO, "Image Quality Stats from %ld psf stars : FWHM (major, minor) [pixels]: %f, %f\n",
               M2->n, fwhm_major, fwhm_minor);
 
-    psLogMsg ("psphot", PS_LOG_INFO, "M_2 : %f +/- %f, M_3 : %f +/- %f, M_4 : %f +/- %f\n",
+    psLogMsg ("psphot", PS_LOG_INFO, "M_2 : %f +/- %f, M_3 : %f +/- %f, M_4 : %f +/- %f  [pixels^n]\n",
               vM2, dM2, vM3, dM3, vM4, dM4);
 

trunk/psphot/src/psphotMagnitudes.c

@@ -71 +71 @@
             PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32);
             PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
+            PS_ARRAY_ADD_SCALAR(job->args, markVal,  PS_TYPE_IMAGE_MASK);
             PS_ARRAY_ADD_SCALAR(job->args, 0,        PS_TYPE_S32); // this is used as a return value for nAp
 
@@ -102 +103 @@
                 fprintf (stderr, "error with job\n");
             } else {
-                psScalar *scalar = job->args->data[7];
+                psScalar *scalar = job->args->data[8];
                 Nap += scalar->data.S32;
             }
@@ -127 +128 @@
     pmSourcePhotometryMode photMode = PS_SCALAR_VALUE(job->args->data[5],S32);
     psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);
+    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[7],PS_TYPE_IMAGE_MASK_DATA);
 
     for (int i = 0; i < sources->n; i++) {
         pmSource *source = (pmSource *) sources->data[i];
-        status = pmSourceMagnitudes (source, psf, photMode, maskVal);
+
+        // replace object in image
+        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
+            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+        }
+
+        // clear the mask bit and set the circular mask pixels
+        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
+        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal);
+
+        status = pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal
         if (status && isfinite(source->apMag)) Nap ++;
+
+        // clear the mask bit 
+        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
+
+        // re-subtract the object, leave local sky
+        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
 
         if (backModel) {
@@ -155 +173 @@
 
     // change the value of a scalar on the array (wrap this and put it in psArray.h)
-    psScalar *scalar = job->args->data[7];
+    psScalar *scalar = job->args->data[8];
     scalar->data.S32 = Nap;
 

trunk/psphot/src/psphotMakeFluxScale.c

@@ -60 +60 @@
         goto DONE;
     }
+    if (trend->mode == PM_TREND_MAP) {
+        // p_psImagePrint (2, trend->map->map, "FluxScale Before"); // XXX TEST:
+        psImageMapRepair (trend->map->map);
+        // p_psImagePrint (2, trend->map->map, "FluxScale After"); // XXX TEST:
+    }
 
     // XXX do something useful to measure residual statistics

trunk/psphot/src/psphotMakeResiduals.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+
+# define RESIDUAL_SOFTENING 0.005 
 
 bool psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf, psImageMaskType maskVal) {
@@ -31 +33 @@
 
     float pixelSN = psMetadataLookupF32(&status, recipe, "PSF.RESIDUALS.PIX.SN");
+    PS_ASSERT (status, false);
+
+    float radiusMax = psMetadataLookupF32(&status, recipe, "PSF.RESIDUALS.RADIUS");
     PS_ASSERT (status, false);
 
@@ -171 +176 @@
                 bool offImage = false;
                 if (psImageInterpolate (&flux, &dflux, &mflux, ix, iy, interp) == PS_INTERPOLATE_STATUS_OFF) {
-                    // fprintf (stderr, "off image: %f %f : %f %f\n", ix, iy, flux, dflux);
                     // This pixel is off the image
                     offImage = true;
@@ -179 +183 @@
                 }
                 fluxes->data.F32[i] = flux;
-                dfluxes->data.F32[i] = dflux;
+                dfluxes->data.F32[i] = hypot(dflux, RESIDUAL_SOFTENING);
                 if (isnan(flux)) {
+                    fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i] = badMask;
+                }
+                if (isnan(dflux)) {
                     fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i] = badMask;
                 }
@@ -234 +241 @@
                 }
 
+                float radius = hypot((ox - 0.5*resid->Ro->numCols), (oy - 0.5*resid->Ro->numRows));
+                if (radius > radiusMax) {
+                  resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
+                  continue;
+                }
+
                 resid->Ro->data.F32[oy][ox] = psStatsGetValue(fluxStats, statOption);
                 resid->Rx->data.F32[oy][ox] = resid->Ry->data.F32[oy][ox] = 0.0;
@@ -248 +261 @@
                   resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
                 }
-
-                // fprintf (stderr, "res: %2d %2d : %6.4f  %6.4f  %6.4f   %3d  %1d\n", ox, oy, resid->Ro->data.F32[oy][ox], fluxStats->sampleStdev, fluxStats->sampleStdev/sqrt(nKeep), nKeep, resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox]);
-
             } else {
                 assert (SPATIAL_ORDER == 1);
+
+                float radius = hypot((ox - 0.5*resid->Ro->numCols), (oy - 0.5*resid->Ro->numRows));
+                if (radius > radiusMax) {
+                  resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
+                  continue;
+                }
+
                 psImageInit(A, 0.0);
                 psVectorInit(B, 0.0);
@@ -275 +292 @@
 
                 if (!psMatrixGJSolve(A, B)) {
-                    psError(PSPHOT_ERR_PSF, false, "Singular matrix solving for (y,x) = (%d,%d)'s residuals",
-                            oy, ox);
-                    psFree(resid); resid = NULL;
-                    break;
+                    resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
+                    psWarning("Singular matrix solving for (y,x) = (%d,%d)'s residuals, masking", oy, ox);
+                    continue;
                 }
 
@@ -286 +302 @@
 
                 float dRo = sqrt(A->data.F32[0][0]);
-                // fprintf (stderr, "res: %2d %2d : %6.4f  %6.4f  %6.4f   %3d  %1d\n",
-                // ox, oy, resid->Ro->data.F32[oy][ox], dRo, dRo/sqrt(nKeep), nKeep, resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox]);
 
                 if (fabs(resid->Ro->data.F32[oy][ox]) < pixelSN*dRo/sqrt(nKeep)) {
                   resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
                 }
-                //resid->variance->data.F32[oy][ox] = XXX;
             }
         }

trunk/psphot/src/psphotMaskReadout.c

@@ -33 +33 @@
     }
 
+    bool softenVariance = psMetadataLookupBool (&status, recipe, "SOFTEN.VARIANCE");
+    float softenFraction = psMetadataLookupF32 (&status, recipe, "SOFTEN.VARIANCE.FRACTION");
+
     // make this an option via the recipe
-    if (0) {
+    if (softenVariance) {
       psImage *im = readout->image;
       psImage *wt = readout->variance;
-      psImage *mk = readout->mask;
       for (int j = 0; j < im->numRows; j++) {
         for (int i = 0; i < im->numCols; i++) {
-          if (isfinite(im->data.F32[j][i]) && isfinite(wt->data.F32[j][i])) continue;
-          mk->data.PS_TYPE_IMAGE_MASK_DATA[j][i] |= maskBad;
+            if (!isfinite(im->data.F32[j][i])) continue;
+            if (!isfinite(wt->data.F32[j][i])) continue;
+            float sysError = softenFraction * im->data.F32[j][i];
+            wt->data.F32[j][i] += PS_SQR(sysError);
         }
       }

trunk/psphot/src/psphotOutput.c

@@ -31 +31 @@
     }
     return background;
+}
+
+// dump source stats for psf stars
+bool psphotDumpStats (psArray *sources, char *stage) {
+
+    char filename[64];
+    snprintf (filename, 64, "psf.%s.dat", stage);
+    FILE *f = fopen (filename, "w");
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = sources->data[i];
+        if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
+
+        pmModel *model = source->modelPSF;
+        if (!model) continue;
+
+        // int xc = source->peak->x - source->pixels->col0;
+        // int yc = source->peak->y - source->pixels->row0;
+        // float mcore = source->modelFlux ? source->modelFlux->data.F32[yc][xc] : NAN;
+        // float mpeak = model ? model->params->data.F32[PM_PAR_I0] : NAN;
+        // bool subtracted = source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED;
+        // fprintf (stderr, "%d %d : %d : %f %f : %f %f\n", source->peak->x, source->peak->y, subtracted, source->peak->flux, source->pixels->data.F32[yc][xc], mcore, mpeak); 
+
+        fprintf (f, "%6.1f %6.1f : %6.1f %6.1f : %8.3f %8.3f %8.3f : %f : %f %f %f : %f\n",
+                 source->peak->xf, source->peak->yf, 
+                 model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS], 
+                 source->psfMag, source->apMag, source->errMag,
+                 model->params->data.F32[PM_PAR_I0], 
+                 model->params->data.F32[PM_PAR_SXX], model->params->data.F32[PM_PAR_SXY], model->params->data.F32[PM_PAR_SYY], 
+                 model->params->data.F32[PM_PAR_7]);
+    }
+    fclose (f);
+    return true;
 }
 
@@ -160 +192 @@
     psMetadataItemSupplement (header, recipe, "DAPMIFIT");
     psMetadataItemSupplement (header, recipe, "NAPMIFIT");
-    psMetadataItemSupplement (header, recipe, "SKYBIAS");
-    psMetadataItemSupplement (header, recipe, "SKYSAT");
 
     // PSF model parameters (shape values for image center)
@@ -256 +286 @@
         psImageKeepCircle (source->maskObj, x, y, radius, "OR", markVal);
         pmModelSub (source->pixels, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL, maskVal);
-        psImageKeepCircle (source->maskObj, x, y, radius, "AND", PS_NOT_IMAGE_MASK(markVal));
+        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
     }
 

trunk/psphot/src/psphotPSFConvModel.c

@@ -37 +37 @@
     }
 
+    // adjust the pixels based on the footprint
+    float radius = psphotSetRadiusEXT (readout, source, markVal);
+    if (!pmSourceMoments (source, radius, 0.0, 0.0)) return false;
+
     // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms:
     psEllipseShape psfShape;
@@ -67 +71 @@
     psVector *params  = modelConv->params;
     psVector *dparams = modelConv->dparams;
-
-    psphotCheckRadiusEXT (readout, source, modelConv, markVal);
 
     // create the minimization constraints

trunk/psphot/src/psphotPetrosian.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
-bool psphotPetrosian (pmSource *source, psMetadata *recipe, psImageMaskType maskVal) {
+bool psphotPetrosian (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal) {
 
-  bool status;
+    // XXX these need to go into recipe values
+    float Rmax = 200;
 
-  assert (source->extpars);
-  assert (source->extpars->profile);
-  assert (source->extpars->profile->radius);
-  assert (source->extpars->profile->flux);
+    psAssert (source->extpars, "need to run psphotRadialProfile first");
+    psAssert (source->extpars->profile, "need to run psphotRadialProfile first");
 
-  psVector *radius = source->extpars->profile->radius;
-  psVector *flux = source->extpars->profile->flux;
+    // integrate the radial profile for radial bins defined for the petrosian measurement:
+    // SB_i (r_i) where \alpha r_i < r < \beta r_i
+    if (!psphotPetrosianRadialBins (source, Rmax, skynoise)) {
+        psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile");
+        return false;
+    }
+  
+    // use the SB_i from above to calculate the petrosian radius and the flux within that radius
+    if (!psphotPetrosianStats (source)) {
+        psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile");
+        return false;
+    }
+  
+    psTrace ("psphot", 3, "source at %f,%f: petrosian radius: %f, flux: %f, axis ratio: %f, angle: %f",
+             source->peak->xf, source->peak->yf, 
+             source->extpars->petrosian_80->radius, 
+             source->extpars->petrosian_80->flux, 
+             source->extpars->profile->axes.minor/source->extpars->profile->axes.major, 
+             source->extpars->profile->axes.theta*PS_DEG_RAD);
 
-  // flux at which to measure isophotal parameters
-  float PETROSIAN_R0 = psMetadataLookupF32 (&status, recipe, "PETROSIAN_R0");
-  float PETROSIAN_RF = psMetadataLookupF32 (&status, recipe, "PETROSIAN_FLUX_RATIO");
-  assert (status);
-
-  // first find flux at R0
-  int firstAbove = -1;
-  int lastBelow = -1;
-  for (int i = 0; i < radius->n; i++) {
-    if (radius->data.F32[i] < PETROSIAN_R0) lastBelow = i;
-    if ((firstAbove < 0) && (radius->data.F32[i] > PETROSIAN_R0)) firstAbove = i;
-  }
-  // if we don't go out far enough, we have a problem...
-  if (lastBelow == radius->n - 1) {
-    psTrace ("psphot", 5, "did not go out far enough to reach petrosian reference radius...");
-    // XXX skip object? raise a flag ?
-    return false;
-  }
-  if (firstAbove < 0) {
-    psTrace ("psphot", 5, "did not go out far enough to bound petrosian reference radius");
-    // XXX raise a flag ?
-    return false;
-  }
-
-  // average flux in this range
-  float fluxR0 = 0.0;
-  int fluxRn = 0;
-  for (int i = PS_MIN(firstAbove, lastBelow); i <= PS_MAX(firstAbove, lastBelow); i++) {
-    fluxR0 += flux->data.F32[i];
-    fluxRn ++;
-  }
-  fluxR0 /= (float)(fluxRn);
-
-  // target flux for petrosian radius
-  float fluxRP = fluxR0 * PETROSIAN_RF;
-
-  // find the first bin below the flux level and the last above the level
-  // XXX can this be done faster with bisection?
-  // XXX do I need to worry about crazy outliers?
-  // XXX should i be smoothing or fitting the curve?
-  int firstBelow = -1;
-  int lastAbove = -1;
-  for (int i = 0; i < flux->n; i++) {
-    if (flux->data.F32[i] > fluxRP) lastAbove = i;
-    if ((firstBelow < 0) && (flux->data.F32[i] < fluxRP)) firstBelow = i;
-  }
-  // if we don't go out far enough, we have a problem...
-  if (lastAbove == radius->n - 1) {
-    psTrace ("psphot", 5, "did not go out far enough to reach petrosian radius...");
-    // XXX skip object? raise a flag ?
-    return false;
-  }
-  if (firstBelow < 0) {
-    psTrace ("psphot", 5, "did not go out far enough to bound petrosian radius");
-    // XXX raise a flag ?
-    return false;
-  }
-
-  // need to examine pixels in this vicinity
-  float fluxFirst = 0;
-  float fluxLast = 0;
-  for (int i = 0; i <= PS_MAX(firstBelow, lastAbove); i++) {
-    if (i <= firstBelow) {
-      fluxFirst += flux->data.F32[i];
-    }
-    if (i <= lastAbove) {
-      fluxLast += flux->data.F32[i];
-    }
-  }
-  float fluxRPSum    = 0.5*(fluxLast + fluxFirst);
-  float fluxRPSumErr = 0.5*fabs(fluxLast - fluxFirst);
-  // XXX need to use the weight appropriately here...
-
-  float rad     = 0.5*(radius->data.F32[firstBelow] + radius->data.F32[lastAbove]);
-  float radErr  = 0.5*fabs(radius->data.F32[firstBelow] - radius->data.F32[lastAbove]);
-
-  if (!source->extpars->petrosian) {
-    source->extpars->petrosian = pmSourcePetrosianValuesAlloc ();
-  }
-
-  // these are uncalibrated: instrumental mags and pixel units
-  source->extpars->petrosian->mag    = -2.5*log10(fluxRPSum);
-  source->extpars->petrosian->magErr = fluxRPSumErr / fluxRPSum;
-
-  source->extpars->petrosian->rad    = rad;
-  source->extpars->petrosian->radErr = radErr;
-
-  psTrace ("psphot", 5, "Petrosian flux:%f +/- %f @ %f +/- %f for %f, %f\n",
-           source->extpars->petrosian->mag, source->extpars->petrosian->magErr,
-           source->extpars->petrosian->rad, source->extpars->petrosian->radErr,
-           source->peak->xf, source->peak->yf);
-
-  return true;
-
+    return true;
 }

trunk/psphot/src/psphotRadialProfile.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
-# define COMPARE_RADIUS(A,B) (radius->data.F32[A] < radius->data.F32[B])
-# define SWAP_RADIUS(TYPE,A,B) { \
-  float tmp; \
-  if (A != B) { \
-    tmp = radius->data.F32[A]; \
-    radius->data.F32[A] = radius->data.F32[B]; \
-    radius->data.F32[B] = tmp; \
-    tmp = flux->data.F32[A]; \
-    flux->data.F32[A] = flux->data.F32[B]; \
-    flux->data.F32[B] = tmp; \
-    tmp = variance->data.F32[A]; \
-    variance->data.F32[A] = variance->data.F32[B]; \
-    variance->data.F32[B] = tmp; \
-  } \
-}
-
-bool psphotRadialProfile (pmSource *source, psMetadata *recipe, psImageMaskType maskVal) {
+bool psphotRadialProfile (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal) {
 
     // allocate pmSourceExtendedParameters, if not already defined
@@ -28 +12 @@
     }
 
-    int nPts = source->pixels->numRows * source->pixels->numCols;
-    source->extpars->profile->radius = psVectorAllocEmpty (nPts, PS_TYPE_F32);
-    source->extpars->profile->flux   = psVectorAllocEmpty (nPts, PS_TYPE_F32);
-    source->extpars->profile->variance = psVectorAllocEmpty (nPts, PS_TYPE_F32);
+    // XXX these need to go into recipe values
+    int Nsec = 24;
+    float Rmax = 200;
+    float fluxMin = 0.0;
+    float fluxMax = source->peak->flux;
 
-    psVector *radius = source->extpars->profile->radius;
-    psVector *flux   = source->extpars->profile->flux;
-    psVector *variance = source->extpars->profile->variance;
+    // generate a series of radial profiles at Nsec evenly spaced angles.  the profile flux
+    // is measured by interpolation for small radii; for large radii, the pixels in a box
+    // are averaged to increase the S/N
+    if (!psphotRadialProfilesByAngles (source, Nsec, Rmax)) {
+        psError (PS_ERR_UNKNOWN, false, "failed to measure radial profile for petrosian");
+        return false;
+    }
 
-    // XXX use the extended source model here for Xo, Yo?
-    // XXX define a radius scaled to the elliptical contour?
+    // use the radial profiles to determine the radius of a given isophote.  this isophote
+    // is used to determine the elliptical shape of the object, so it has a relatively high
+    // value (nominally 50% of the peak)
+    if (!psphotRadiiFromProfiles (source, fluxMin, fluxMax)) {
+        psError (PS_ERR_UNKNOWN, false, "failed to measure isophotal radii from profiles");
+        return false;
+    }
 
-    int n = 0;
-
-    float Xo = 0.0;
-    float Yo = 0.0;
-
-    if (source->modelEXT) {
-      Xo = source->modelEXT->params->data.F32[PM_PAR_XPOS] - source->pixels->col0;
-      Yo = source->modelEXT->params->data.F32[PM_PAR_YPOS] - source->pixels->row0;
-    } else {
-      Xo = source->peak->xf - source->pixels->col0;
-      Yo = source->peak->yf - source->pixels->row0;
+    // convert the isophotal radius vs angle measurements to an elliptical contour
+    if (!psphotEllipticalContour (source)) {
+        psLogMsg ("psphot", 3, "failed to measure elliptical contour");
+        return false;
     }
-    for (int iy = 0; iy < source->pixels->numRows; iy++) {
-        for (int ix = 0; ix < source->pixels->numCols; ix++) {
-            if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) continue;
-            radius->data.F32[n] = hypot (ix - Xo, iy - Yo) ;
-            flux->data.F32[n]   = source->pixels->data.F32[iy][ix];
-            variance->data.F32[n] = source->variance->data.F32[iy][ix];
-            n++;
-        }
+  
+    // generate a single, normalized radial profile following the elliptical contours.
+    // the radius is normalized by the axis ratio so that on the major axis, 1 pixel = 1 pixel
+    if (!psphotEllipticalProfile (source)) {
+        psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile");
+        return false;
     }
-    radius->n = n;
-    variance->n = n;
-    flux->n = n;
-
-    // sort the vector set by the radius
-    PSSORT (radius->n, COMPARE_RADIUS, SWAP_RADIUS, NONE);
-
+  
     return true;
 }

trunk/psphot/src/psphotRadiusChecks.c

@@ -7 +7 @@
                                         // and a per-object radius is calculated)
 
+static float PSF_APERTURE = 0;  // radius to use in PSF aperture mags
+
+
 bool psphotInitRadiusPSF(const psMetadata *recipe, const pmModelType type) {
 
@@ -13 +16 @@
     PSF_FIT_NSIGMA  = psMetadataLookupF32(&status, recipe, "PSF_FIT_NSIGMA");
     PSF_FIT_PADDING = psMetadataLookupF32(&status, recipe, "PSF_FIT_PADDING");
-    PSF_FIT_RADIUS =  psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS");
+    PSF_FIT_RADIUS  =  psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS");
+    PSF_APERTURE    =  psMetadataLookupF32(&status, recipe, "PSF_APERTURE");
 
     return true;
@@ -34 +38 @@
             radiusFit = model->modelRadius(model->params, 1.0);
         }
+        model->fitRadius = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING);
+    } else {
+        model->fitRadius = radiusFit;
     }
-    model->radiusFit = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING);
-
-    if (isnan(model->radiusFit)) psAbort("error in radius");
+    if (isnan(model->fitRadius)) psAbort("error in radius");
         
     if (source->mode & PM_SOURCE_MODE_SATSTAR) {
-        model->radiusFit *= 2;
+        model->fitRadius *= 2;
     }
 
-    bool status = pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->radiusFit);
+    // radius used to measure aperture photometry
+    source->apRadius = PSF_APERTURE;
+
+    bool status = pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius);
 
     // set the mask to flag the excluded pixels
-    psImageKeepCircle (source->maskObj, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->radiusFit, "OR", markVal);
+    psImageKeepCircle (source->maskObj, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius, "OR", markVal);
     return status;
 }
@@ -58 +66 @@
 
     // set the fit radius based on the object flux limit and the model
-    model->radiusFit = (RADIUS_TYPE) (model->modelRadius (model->params, PSF_FIT_NSIGMA*moments->dSky) + dR + PSF_FIT_PADDING);
-    if (isnan(model->radiusFit)) psAbort("error in radius");
-        
+    float radiusFit = PSF_FIT_RADIUS;
+    if (radiusFit <= 0) {               // use fixed radius
+        if (moments == NULL) {
+            radiusFit = model->modelRadius(model->params, PSF_FIT_NSIGMA*moments->dSky);
+        } else {
+            radiusFit = model->modelRadius(model->params, 1.0);
+        }
+        model->fitRadius = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING);
+    } else {
+        model->fitRadius = radiusFit;
+    }
+    if (isnan(model->fitRadius)) psAbort("error in radius");
+
+    // above sets a radius for a single star, bump by blend separation
+    model->fitRadius += dR;
+
     if (source->mode &  PM_SOURCE_MODE_SATSTAR) {
-        model->radiusFit *= 2;
+        model->fitRadius *= 2;
     }
 
-    bool status = pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->radiusFit);
+    // radius used to measure aperture photometry
+    source->apRadius = PSF_APERTURE;
+
+    bool status = pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius);
 
     // set the mask to flag the excluded pixels
-    psImageKeepCircle (source->maskObj, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->radiusFit, "OR", markVal);
+    psImageKeepCircle (source->maskObj, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius, "OR", markVal);
     return status;
 }
@@ -74 +98 @@
 static float EXT_FIT_NSIGMA;
 static float EXT_FIT_PADDING;
+static float EXT_FIT_MAX_RADIUS;
 
 bool psphotInitRadiusEXT (psMetadata *recipe, pmModelType type) {
@@ -79 +104 @@
     bool status;
 
-    EXT_FIT_NSIGMA   = psMetadataLookupF32 (&status, recipe, "EXT_FIT_NSIGMA");
-    EXT_FIT_PADDING  = psMetadataLookupF32 (&status, recipe, "EXT_FIT_PADDING");
+    EXT_FIT_NSIGMA     = psMetadataLookupF32 (&status, recipe, "EXT_FIT_NSIGMA");
+    EXT_FIT_PADDING    = psMetadataLookupF32 (&status, recipe, "EXT_FIT_PADDING");
+    EXT_FIT_MAX_RADIUS = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_RADIUS");
 
     return true;
@@ -86 +112 @@
 
 // call this function whenever you (re)-define the EXT model
+float psphotSetRadiusEXT (pmReadout *readout, pmSource *source, psImageMaskType markVal) {
+
+    psAssert (source, "source not defined??");
+    psAssert (source->peak, "peak not defined??");
+
+    pmPeak *peak = source->peak;
+
+    // set the radius based on the footprint:
+    if (!peak->footprint) goto escape;
+    pmFootprint *footprint = peak->footprint;
+    if (!footprint->spans) goto escape;
+    if (footprint->spans->n < 1) goto escape;
+
+    // find the max radius
+    float radius = 0.0;
+    for (int j = 0; j < footprint->spans->n; j++) {
+        pmSpan *span = footprint->spans->data[j];
+
+        float dY  = span->y  - peak->yf;
+        float dX0 = span->x0 - peak->xf;
+        float dX1 = span->x1 - peak->xf;
+
+        radius = PS_MAX (radius, hypot(dY, dX0));
+        radius = PS_MAX (radius, hypot(dY, dX1));
+    }
+
+    radius += EXT_FIT_PADDING;
+    if (isnan(radius)) psAbort("error in radius");
+
+    radius = PS_MIN (radius, EXT_FIT_MAX_RADIUS);
+
+    // redefine the pixels if needed
+    pmSourceRedefinePixels (source, readout, peak->xf, peak->yf, radius);
+
+    // set the mask to flag the excluded pixels
+    psImageKeepCircle (source->maskObj, peak->xf, peak->yf, radius, "OR", markVal);
+    return radius;
+
+escape:
+    return NAN;
+    // bool result = psphotCheckRadiusEXT (readout, source, model, markVal);
+    // return result;
+}
+
+// alternative EXT radius based on model guess (for use without footprints)
 bool psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal) {
+
+    psAbort ("do not use this function");
 
     psF32 *PAR = model->params->data.F32;
@@ -96 +169 @@
     float rawRadius = model->modelRadius (model->params, EXT_FIT_NSIGMA*moments->dSky);
 
-    model->radiusFit = rawRadius + EXT_FIT_PADDING;
-    if (isnan(model->radiusFit)) psAbort("error in radius");
+    model->fitRadius = rawRadius + EXT_FIT_PADDING;
+    if (isnan(model->fitRadius)) psAbort("error in radius");
 
     // redefine the pixels if needed
-    bool status = pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->radiusFit);
+    bool status = pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius);
 
     // set the mask to flag the excluded pixels
-    psImageKeepCircle (source->maskObj, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->radiusFit, "OR", markVal);
+    psImageKeepCircle (source->maskObj, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius, "OR", markVal);
     return status;
 }

trunk/psphot/src/psphotReadout.c

@@ -81 +81 @@
 
     // construct sources and measure basic stats
-    psArray *sources = psphotSourceStats (config, readout, detections);
+    psArray *sources = psphotSourceStats (config, readout, detections, true);
     if (!sources) return false;
     if (!strcasecmp (breakPt, "PEAKS")) {
@@ -126 +126 @@
         return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
     }
-
     psphotVisualShowPSFModel (readout, psf);
 
@@ -132 +131 @@
     psphotLoadExtSources (config, view, sources);
 
-    // construct an initial model for each object
+    // construct an initial model for each object, set the radius to fitRadius, set circular fit mask
     psphotGuessModels (config, readout, sources, psf);
 
-    // XXX test output of models
-    // psphotTestSourceOutput (readout, sources, recipe, psf);
-
-    // linear PSF fit to source peaks
+    // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
     psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE);
 
@@ -144 +140 @@
     // psphotGuessModels or fitted until psphotFitSourcesLinear.
     psphotVisualShowPSFStars (recipe, psf, sources);
-    psphotVisualShowSatStars (recipe, psf, sources);
 
     // identify CRs and extended sources
-    psphotSourceSize (config, readout, sources, recipe, 0);
+    psphotSourceSize (config, readout, sources, recipe, psf, 0);
     if (!strcasecmp (breakPt, "ENSEMBLE")) {
         goto finish;
     }
+    psphotVisualShowSatStars (recipe, psf, sources);
 
     // non-linear PSF and EXT fit to brighter sources
+    // replace model flux, adjust mask as needed, fit, subtract the models (full stamp)
     psphotBlendFit (config, readout, sources, psf);
 
@@ -158 +155 @@
     psphotReplaceAllSources (sources, recipe);
 
-    // linear fit to include all sources
+    // linear fit to include all sources (subtract again)
     psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE);
 
@@ -165 +162 @@
         goto pass1finish;
     }
-
-    // XXX for the moment, drop the re-calc of the background (prove this works)
-    // replace background in residual image
-    // psphotSkyReplace (config, view);
-    // re-measure background model (median, smoothed image)
-    // psphotImageMedian (config, view);
+    // NOTE: possibly re-measure background model here with objects subtracted
 
     // add noise for subtracted objects
@@ -182 +174 @@
 
     // define new sources based on only the new peaks
-    psArray *newSources = psphotSourceStats (config, readout, detections);
+    psArray *newSources = psphotSourceStats (config, readout, detections, false);
 
     // set source type
@@ -190 +182 @@
     }
 
-    // create full input models
+    // create full input models, set the radius to fitRadius, set circular fit mask
     psphotGuessModels (config, readout, newSources, psf);
 
@@ -206 +198 @@
 
     // measure source size for the remaining sources
-    psphotSourceSize (config, readout, sources, recipe, 0);
+    psphotSourceSize (config, readout, sources, recipe, psf, 0);
 
     psphotExtendedSourceAnalysis (readout, sources, recipe);

trunk/psphot/src/psphotReadoutFindPSF.c

@@ -41 +41 @@
 
     // construct sources and measure basic stats (moments, local sky)
-    psArray *sources = psphotSourceStats(config, readout, detections);
+    psArray *sources = psphotSourceStats(config, readout, detections, true);
     if (!sources) return false;
 

trunk/psphot/src/psphotReadoutKnownSources.c

@@ -41 +41 @@
 
     // construct sources and measure basic stats
-    psArray *sources = psphotSourceStats (config, readout, detections);
+    psArray *sources = psphotSourceStats (config, readout, detections, true);
     if (!sources) return false;
 

trunk/psphot/src/psphotReadoutMinimal.c

@@ -56 +56 @@
 
     // construct sources and measure basic stats
-    psArray *sources = psphotSourceStats (config, readout, detections);
+    psArray *sources = psphotSourceStats (config, readout, detections, true);
     if (!sources) return false;
 

trunk/psphot/src/psphotReplaceUnfit.c

@@ -17 +17 @@
     replace:
         pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
-        source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
     }
     psLogMsg ("psphot.replace", 3, "replace unfitted models: %f sec (%ld objects)\n", psTimerMark ("psphot.replace"), sources->n);
@@ -41 +40 @@
 
       pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
-      source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
     }
     psLogMsg ("psphot.replace", PS_LOG_INFO, "replaced models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.replace"));
@@ -64 +62 @@
       if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) continue;
 
-      pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
-      source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
+      pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     }
     psLogMsg ("psphot.replace", PS_LOG_INFO, "replaced models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.replace"));
@@ -71 +68 @@
 }
 
+# if (0)
 // add source, if the source has been subtracted; do not modify state
 bool psphotAddWithTest (pmSource *source, bool useState, psImageMaskType maskVal) {
@@ -108 +106 @@
     return true;
 }
+# endif

trunk/psphot/src/psphotRoughClass.c

@@ -26 +26 @@
         for (int iy = 0; iy < NY; iy ++) {
 
-            psRegion region = psRegionSet (ix*dX, (ix + 1)*dX, iy*dY, (iy + 1)*dY);
-            if (!psphotRoughClassRegion (nRegion, &region, sources, recipe, havePSF)) {
+            psRegion *region = psRegionAlloc (ix*dX, (ix + 1)*dX, iy*dY, (iy + 1)*dY);
+            if (!psphotRoughClassRegion (nRegion, region, sources, recipe, havePSF)) {
                 psLogMsg ("psphot", 4, "Failed to determine rough classification for region %f,%f - %f,%f\n", 
-                         region.x0, region.y0, region.x1, region.y1);
+                         region->x0, region->y0, region->x1, region->y1);
+                psFree (region);
                 continue;
             }
-            
+            psFree (region);
             nRegion ++;
         }
@@ -45 +46 @@
     psphotVisualPlotMoments (recipe, sources);
     psphotVisualShowRoughClass (sources);
-    psphotVisualShowFlags (sources);
+    // XXX better visualization: psphotVisualShowFlags (sources);
 
     return true;
@@ -63 +64 @@
     psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
     if (!regionMD) {
+        // allocate the region metadata folder and add this region to it.
         regionMD = psMetadataAlloc();
         psMetadataAddMetadata (recipe, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD);
         psFree (regionMD);
     }
+    psMetadataAddPtr (regionMD, PS_LIST_TAIL, "REGION", PS_DATA_REGION | PS_META_REPLACE, "psf clump region", region);
 
     if (!havePSF) {
         // determine the PSF parameters from the source moment values
+        // XXX why not save the psfClump as a PTR?
         psfClump = pmSourcePSFClump (region, sources, recipe);
         psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X",  PS_META_REPLACE, "psf clump center", psfClump.X);

trunk/psphot/src/psphotSetThreads.c

@@ -10 +10 @@
     psFree(task);
 
-    task = psThreadTaskAlloc("PSPHOT_MAGNITUDES", 8);
+    task = psThreadTaskAlloc("PSPHOT_MAGNITUDES", 9);
     task->function = &psphotMagnitudes_Threaded;
     psThreadTaskAdd(task);
@@ -20 +20 @@
     psFree(task);
 
-    task = psThreadTaskAlloc("PSPHOT_APRESID_MAGS", 6);
+    task = psThreadTaskAlloc("PSPHOT_APRESID_MAGS", 7);
     task->function = &psphotApResidMags_Threaded;
     psThreadTaskAdd(task);

trunk/psphot/src/psphotSourceFits.c

@@ -90 +90 @@
     psphotCheckRadiusPSFBlend (readout, source, PSF, markVal, dR);
 
-    // fit PSF model (set/unset the pixel mask)
+    // fit PSF model
     pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF, maskVal);
+
+    // clear the circular mask
+    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); 
+
+    if (!isfinite(PSF->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
 
     // correct model chisq for flux trend
@@ -101 +106 @@
         pmSource *blend = sourceSet->data[i];
         pmModel *model  = modelSet->data[i];
+
+        if (!isfinite(model->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
 
         // correct model chisq for flux trend
@@ -120 +127 @@
         pmSourceCacheModel (blend, maskVal);
         pmSourceSub (blend, PM_MODEL_OP_FULL, maskVal);
-        blend->tmpFlags |=  PM_SOURCE_TMPF_SUBTRACTED;
         blend->mode |=  PM_SOURCE_MODE_BLEND_FIT;
     }
@@ -144 +150 @@
     pmSourceCacheModel (source, maskVal);
     pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-    source->tmpFlags |=  PM_SOURCE_TMPF_SUBTRACTED;
     source->mode |=  PM_SOURCE_MODE_BLEND_FIT;
     return true;
@@ -167 +172 @@
     // fit PSF model (set/unset the pixel mask)
     pmSourceFitModel (source, PSF, PM_SOURCE_FIT_PSF, maskVal);
+
+    if (!isfinite(PSF->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
+
+    // clear the circular mask
+    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); 
 
     // correct model chisq for flux trend
@@ -186 +196 @@
     pmSourceCacheModel (source, maskVal);
     pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-    source->tmpFlags |=  PM_SOURCE_TMPF_SUBTRACTED;
-    return true;
-}
-
-static float EXT_MIN_SN;
-static float EXT_MOMENTS_RAD;
+    return true;
+}
+
 static pmModelType modelTypeEXT;
 
@@ -197 +204 @@
 
     bool status;
-
-    // extended source model parameters
-    EXT_MIN_SN       = psMetadataLookupF32 (&status, recipe, "EXT_MIN_SN");
-    EXT_MOMENTS_RAD  = psMetadataLookupF32 (&status, recipe, "EXT_MOMENTS_RADIUS");
 
     // extended source model descriptions
@@ -221 +224 @@
     if (source->type == PM_SOURCE_TYPE_DEFECT) return false;
     if (source->type == PM_SOURCE_TYPE_SATURATED) return false;
-    if (source->peak->SN < EXT_MIN_SN) return false;
-
+
+    // set the radius based on the footprint (also sets the mask pixels)
+    float radius = psphotSetRadiusEXT (readout, source, markVal);
+
+    // XXX note that this changes the source moments that are published...
     // recalculate the source moments using the larger extended-source moments radius
     // at this stage, skip Gaussian windowing, and do not clip pixels by S/N
-    if (!pmSourceMoments (source, EXT_MOMENTS_RAD, 0.0, 0.0)) return false;
+    // this uses the footprint to judge both radius and aperture?
+    if (!pmSourceMoments (source, radius, 0.0, 0.0)) return false;
 
     psTrace ("psphot", 5, "trying blob...\n");
@@ -237 +244 @@
     // XXX need to handle failures better here
     pmModel *EXT = psphotFitEXT (readout, source, modelTypeEXT, maskVal, markVal);
+    if (!isfinite(EXT->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
+
     okEXT = psphotEvalEXT (tmpSrc, EXT);
     chiEXT = EXT ? EXT->chisq / EXT->nDOF : NAN;
@@ -246 +255 @@
     // XXX should I keep / save the flags set in the eval functions?
 
+    // clear the circular mask
+    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); 
+
     // correct first model chisqs for flux trend
     chiDBL = NAN;
     ONE = DBL->data[0];
     if (ONE) {
+        if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
       chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
       ONE->chisqNorm = ONE->chisq / chiTrend;
@@ -258 +271 @@
     ONE = DBL->data[1];
     if (ONE) {
+        if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
       chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
       ONE->chisqNorm = ONE->chisq / chiTrend;
@@ -277 +291 @@
 
     // both models failed; reject them both
+    // XXX -- change type flags to psf in this case and keep original moments?
     psFree (EXT);
     psFree (DBL);
@@ -287 +302 @@
     // save new model
     source->modelEXT = EXT;
+    source->modelEXT->fitRadius = radius;
     source->type = PM_SOURCE_TYPE_EXTENDED;
     source->mode |= PM_SOURCE_MODE_EXTMODEL;
@@ -293 +309 @@
     pmSourceCacheModel (source, maskVal);
     pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-    source->tmpFlags |=  PM_SOURCE_TMPF_SUBTRACTED;
+
+# if (PS_TRACE_ON)   
     psTrace ("psphot", 5, "blob as EXT: %f %f\n", EXT->params->data.F32[PM_PAR_XPOS], EXT->params->data.F32[PM_PAR_YPOS]);
+    if (psTraceGetLevel("psphot") >= 6) {
+        psLogMsg ("psphot", 1, "source 2:\n");
+        for (int i = 0; i < source->modelEXT->params->n; i++) {
+            psLogMsg ("psphot", 1, "PAR %d : %f +/- %f\n", i, source->modelEXT->params->data.F32[i], source->modelEXT->dparams->data.F32[i]);
+        }
+    }
+# endif
+
     return true;
 
@@ -304 +329 @@
     psFree (source->modelPSF);
     source->modelPSF = psMemIncrRefCounter (DBL->data[0]);
-    source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
     source->mode     |= PM_SOURCE_MODE_PAIR;
+    source->modelPSF->fitRadius = radius;
 
     // copy most data from the primary source (modelEXT, blends stay NULL)
-    // XXX use pmSourceCopy?
     pmSource *newSrc = pmSourceCopy (source);
     newSrc->modelPSF = psMemIncrRefCounter (DBL->data[1]);
+    newSrc->modelPSF->fitRadius = radius;
 
     // build cached models and subtract
@@ -317 +342 @@
     pmSourceCacheModel (newSrc, maskVal);
     pmSourceSub (newSrc, PM_MODEL_OP_FULL, maskVal);
+
+# if (PS_TRACE_ON)   
     psTrace ("psphot", 5, "blob as DBL: %f %f\n", ONE->params->data.F32[PM_PAR_XPOS], ONE->params->data.F32[PM_PAR_YPOS]);
+    if (psTraceGetLevel("psphot") >= 6) {
+        psLogMsg ("psphot", 1, "source 1:\n");
+        for (int i = 0; i < newSrc->modelPSF->params->n; i++) {
+            psLogMsg ("psphot", 1, "PAR %d : %f +/- %f\n", i, newSrc->modelPSF->params->data.F32[i], newSrc->modelPSF->dparams->data.F32[i]);
+        }
+        psLogMsg ("psphot", 1, "source 2:\n");
+        for (int i = 0; i < source->modelPSF->params->n; i++) {
+            psLogMsg ("psphot", 1, "PAR %d : %f +/- %f\n", i, source->modelPSF->params->data.F32[i], source->modelPSF->dparams->data.F32[i]);
+        }
+        psphotVisualShowResidualImage (readout);
+    }
+# endif
 
     psArrayAdd (newSources, 100, newSrc);
@@ -356 +395 @@
     // save the PSF model from the Ensemble fit
     PSF = source->modelPSF;
-    psphotCheckRadiusPSFBlend (readout, source, PSF, markVal, 8.0);
     if (isnan(PSF->params->data.F32[1])) psAbort("nan in dbl fit");
 
@@ -389 +427 @@
     PS_ASSERT (EXT, NULL);
 
-    // if (isnan(EXT->params->data.F32[1])) psAbort("nan in ext fit");
-
-    psphotCheckRadiusEXT (readout, source, EXT, markVal);
-
     if ((source->moments->Mxx < 1e-3) || (source->moments->Myy < 1e-3)) {
         psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy);
@@ -401 +435 @@
     return (EXT);
 }
-

trunk/psphot/src/psphotSourcePlots.c

@@ -111 +111 @@
             if (Xo == 0) {
                 // place source alone on this row
-                psphotAddWithTest (source, true, maskVal); // replace source if subtracted
+                bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
+                if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
                 psphotRadialPlot (&kapa, "radial.plots.ps", source);
                 psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
 
-                psphotSubWithTest (source, false, maskVal); // remove source (force)
+                pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
                 psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
 
-                psphotSetState (source, false, maskVal); // replace source (has been subtracted)
+                if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+
                 Yo += DY;
                 Xo = 0;
@@ -126 +128 @@
                 Yo += dY;
                 Xo = 0;
-                psphotAddWithTest (source, true, maskVal); // replace source if subtracted
+
+                bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
+                if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
                 psphotRadialPlot (&kapa, "radial.plots.ps", source);
                 psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
 
-                psphotSubWithTest (source, false, maskVal); // remove source (force)
+                pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
                 psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
-                psphotSetState (source, false, maskVal); // replace source (has been subtracted)
+                if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
 
                 Xo = DX;
@@ -139 +143 @@
         } else {
             // extend this row
-            psphotAddWithTest (source, true, maskVal); // replace source if subtracted
+            bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
+            if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
             psphotRadialPlot (&kapa, "radial.plots.ps", source);
             psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
 
-            psphotSubWithTest (source, false, maskVal); // remove source (force)
+            pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
             psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
-            psphotSetState (source, false, maskVal); // replace source (has been subtracted)
+            if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
 
             Xo += DX;

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;
-}

trunk/psphot/src/psphotSourceStats.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
-psArray *psphotSourceStats (pmConfig *config, pmReadout *readout, pmDetections *detections) {
+bool psphotSetMomentsWindow (psMetadata *recipe, psArray *sources);
+
+psArray *psphotSourceStats (pmConfig *config, pmReadout *readout, pmDetections *detections, bool setWindow) {
 
     bool status = false;
@@ -21 +23 @@
     float OUTER    = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
     if (!status) return NULL;
+
+    OUTER = PS_MAX(OUTER, 20.0); // XXX Guarantee that we can encompass the max moments radius
+
     char *breakPt  = psMetadataLookupStr (&status, recipe, "BREAK_POINT");
     if (!status) return NULL;
@@ -60 +65 @@
         psphotVisualShowMoments (sources);
         return sources;
+    }
+
+    if (setWindow) {
+        if (!psphotSetMomentsWindow(recipe, sources)) {
+            psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!");
+            return NULL;
+        }
     }
 
@@ -144 +156 @@
     float RADIUS       = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
     if (!status) return false;
-    float MIN_PIXEL_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_MIN_PIXEL_SN");
-    if (!status) return false;
     float SIGMA        = psMetadataLookupF32 (&status, recipe, "MOMENTS_GAUSS_SIGMA");
     if (!status) return false;
@@ -194 +204 @@
         }
 
-        // measure basic source moments
-        status = pmSourceMoments (source, RADIUS, SIGMA, MIN_PIXEL_SN);
+        // measure basic source moments (no S/N clipping on input pixels)
+        status = pmSourceMoments (source, RADIUS, SIGMA, 0.0);
         if (status) {
             Nmoments ++;
@@ -205 +215 @@
         BIG_RADIUS = PS_MIN (INNER, 3*RADIUS);
         psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y);
-        status = pmSourceMoments (source, BIG_RADIUS, 3.0*SIGMA, MIN_PIXEL_SN);
+        status = pmSourceMoments (source, BIG_RADIUS, 3.0*SIGMA, 0.0);
         if (status) {
             source->mode |= PM_SOURCE_MODE_BIG_RADIUS;
@@ -231 +241 @@
 }
 
-# if (0)
-bool psphotSourceStats_Unthreaded (int *nfail, int *nmoments, psArray *sources, psMetadata *recipe) {
-
-    bool status = false;
-    float BIG_RADIUS;
-
-    float INNER    = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
-    if (!status) return false;
-    float MIN_SN   = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
-    if (!status) return false;
-    float RADIUS   = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
-    if (!status) return false;
-
-    // bit-masks to test for good/bad pixels
-    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
-    assert (maskVal);
-
-    // bit-mask to mark pixels not used in analysis
-    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
-    assert (markVal);
-
-    // maskVal is used to test for rejected pixels, and must include markVal
-    maskVal |= markVal;
-
-    // threaded measurement of the sources moments
-    int Nfail = 0;
-    int Nmoments = 0;
-    for (int i = 0; i < sources->n; i++) {
-        pmSource *source = sources->data[i];
-
-        // skip faint sources for moments measurement
-        if (source->peak->SN < MIN_SN) {
-            continue;
-        }
-
-        // measure a local sky value
-        // the local sky is now ignored; kept here for reference only
-        status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
-        if (!status) {
-            psErrorClear(); // XXX re-consider the errors raised here
-            Nfail ++;
-            continue;
-        }
-
-        // measure the local sky variance (needed if noise is not sqrt(signal))
-        // XXX EAM : this should use ROBUST not SAMPLE median, but it is broken
-        status = pmSourceLocalSkyVariance (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
-        if (!status) {
-            Nfail ++;
-            psErrorClear();
-            continue;
-        }
-
-        // measure basic source moments
-        status = pmSourceMoments (source, RADIUS, SIGMA, MIN_PIXEL_SN);
-        if (status) {
-            Nmoments ++;
-            continue;
-        }
-
-        // if no valid pixels, or massive swing, likely saturated source,
-        // try a much larger box
-        BIG_RADIUS = PS_MIN (INNER, 3*RADIUS);
-        psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y);
-        status = pmSourceMoments (source, BIG_RADIUS, 3.0*SIGMA, MIN_PIXEL_SN);
-        if (status) {
-            Nmoments ++;
-            continue;
-        }
-
-        Nfail ++;
-        psErrorClear();
-        continue;
-    }
-
-    // change the value of a scalar on the array (wrap this and put it in psArray.h)
-    *nmoments = Nmoments;
-    *nfail = Nfail;
-
+// this function attempts to iteratively determine the best value for sigma of the moments weighting Gaussian
+bool psphotSetMomentsWindow (psMetadata *recipe, psArray *sources) {
+
+    bool status;
+
+    float MIN_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
+    if (!status) return false;
+
+    // XXX move this to a config file?
+    float sigma[4] = {1.0, 2.0, 3.0, 4.5};
+    float Sout[4];
+
+    // this sorts by peak->SN
+    sources = psArraySort (sources, pmSourceSortBySN);
+
+    // loop over radii:
+    for (int i = 0; i < 4; i++) {
+
+        // XXX move max source number to config
+        for (int j = 0; (j < sources->n) && (j < 400); j++) {
+ 
+            pmSource *source = sources->data[j];
+            psAssert (source->moments, "force moments to exist");
+            source->moments->nPixels = 0;
+
+            // skip faint sources for moments measurement
+            if (source->peak->SN < MIN_SN) {
+                continue;
+            }
+
+            // measure basic source moments (no S/N clipping on input pixels)
+            status = pmSourceMoments (source, 4*sigma[i], sigma[i], 0.0);
+        }
+
+        // choose a grid scale that is a fixed fraction of the psf sigma^2
+        psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(sigma[i]));
+        psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i]));
+        psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i]));
+
+        // determine the PSF parameters from the source moment values
+        pmPSFClump psfClump = pmSourcePSFClump (NULL, sources, recipe);
+        psLogMsg ("psphot", 3, "radius %.1f, nStars: %d, nSigma: %5.2f, X,  Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]);
+
+        psMetadataAddS32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", 1);
+        psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, "PSF.CLUMP.REGION.000");
+        if (!regionMD) {
+            regionMD = psMetadataAlloc();
+            psMetadataAddMetadata (recipe, PS_LIST_TAIL, "PSF.CLUMP.REGION.000", PS_META_REPLACE, "psf clump region", regionMD);
+            psFree (regionMD);
+        }
+        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X",  PS_META_REPLACE, "psf clump center", psfClump.X);
+        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y",  PS_META_REPLACE, "psf clump center", psfClump.Y);
+        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX);
+        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY);
+            
+        // psphotVisualPlotMoments (recipe, sources);
+
+        Sout[i] = sqrt(0.5*(psfClump.X + psfClump.Y)) / sigma[i];
+    }
+
+    // we are looking for sigma for which Sout = 0.65 (or some other value)
+
+    float Sigma = NAN;
+    float minS = Sout[0];
+    float maxS = Sout[0];
+    for (int i = 0; i < 4; i++) {
+        minS = PS_MIN(Sout[i], minS);
+        maxS = PS_MAX(Sout[i], maxS);
+    }
+    if (minS > 0.65) Sigma = sigma[3];
+    if (maxS < 0.65) Sigma = sigma[0];
+
+    for (int i = 0; i < 3; i++) {
+        if ((Sout[i] > 0.65) && (Sout[i+1] > 0.65)) continue;
+        if ((Sout[i] < 0.65) && (Sout[i+1] < 0.65)) continue;
+        Sigma = sigma[i] + (0.65 - Sout[i])*(sigma[i+1] - sigma[i])/(Sout[i+1] - Sout[i]);
+    }
+    psAssert (isfinite(Sigma), "did we miss a case?");
+        
+    // choose a grid scale that is a fixed fraction of the psf sigma^2
+    psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(Sigma));
+    psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma));
+    psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma));
+    psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_GAUSS_SIGMA", PS_META_REPLACE, "moments limit", Sigma);
+    psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_MOMENTS_RADIUS", PS_META_REPLACE, "moments limit", 4.0*Sigma);
+
+    psLogMsg ("psphot", 3, "using window function with sigma = %f\n", Sigma);
     return true;
 }
-# endif

trunk/psphot/src/psphotVisual.c

@@ -15 +15 @@
 # include <kapa.h>
 
+bool pmVisualLimitsFromVectors (Graphdata *graphdata, psVector *xVec, psVector *yVec);
+
 // functions used to visualize the analysis as it goes
 // these are invoked by the -visual options
 
-static int kapa = -1;
+static int kapa1 = -1;
 static int kapa2 = -1;
 static int kapa3 = -1;
 
+int psphotKapaChannel (int channel) {
+
+    switch (channel) {
+      case 1:
+        if (kapa1 == -1) {
+            kapa1 = KapaOpenNamedSocket ("kapa", "psphot:images");
+            if (kapa1 == -1) {
+                fprintf (stderr, "failure to open kapa; visual mode disabled\n");
+                pmVisualSetVisual(false);
+            }
+        }
+        return kapa1;
+      case 2:
+        if (kapa2 == -1) {
+            kapa2 = KapaOpenNamedSocket ("kapa", "psphot:plots");
+            if (kapa2 == -1) {
+                fprintf (stderr, "failure to open kapa; visual mode disabled\n");
+                pmVisualSetVisual(false);
+            }
+        }
+        return kapa2;
+      case 3:
+        if (kapa3 == -1) {
+            kapa3 = KapaOpenNamedSocket ("kapa", "psphot:stamps");
+            if (kapa3 == -1) {
+                fprintf (stderr, "failure to open kapa; visual mode disabled\n");
+                pmVisualSetVisual(false);
+            }
+        }
+        return kapa3;
+      default:
+        psAbort ("unknown kapa channel");
+    }
+    psAbort ("unknown kapa channel");
+}
 
 bool psphotVisualShowMask (int kapaFD, psImage *inImage, const char *name, int channel) {
@@ -131 +168 @@
     if (!pmVisualIsVisual()) return true;
 
-    if (kapa == -1) {
-        kapa = KapaOpenNamedSocket ("kapa", "psphot:images");
-        if (kapa == -1) {
-            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
-            pmVisualSetVisual(false);
-            return false;
-        }
-    }
+    int kapa = psphotKapaChannel (1);
+    if (kapa == -1) return false;
 
     // psphotVisualShowMask (kapa, readout->mask, "mask", 2);
@@ -144 +175 @@
     psphotVisualScaleImage (kapa, readout->image, "image", 0);
 
-    // pause and wait for user input:
-    // continue, save (provide name), ??
-    char key[10];
-    fprintf (stdout, "[c]ontinue? ");
-    if (!fgets(key, 8, stdin)) {
-        psWarning("Unable to read option");
-    }
+    pmVisualAskUser(NULL);
     return true;
 }
@@ -160 +185 @@
     if (!pmVisualIsVisual()) return true;
 
-    if (kapa == -1) {
-        kapa = KapaOpenNamedSocket ("kapa", "psphot:images");
-        if (kapa == -1) {
-            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
-            pmVisualSetVisual(false);
-            return false;
-        }
-    }
+    int kapa = psphotKapaChannel (1);
+    if (kapa == -1) return false;
 
     bool status = false;
@@ -181 +200 @@
     psphotVisualScaleImage (kapa, readout->image, "backsub", 0);
 
-    // pause and wait for user input:
-    // continue, save (provide name), ??
-    char key[10];
-    fprintf (stdout, "[c]ontinue? ");
-    if (!fgets(key, 8, stdin)) {
-        psWarning("Unable to read option");
-    }
+    pmVisualAskUser(NULL);
     return true;
 }
@@ -195 +208 @@
     if (!pmVisualIsVisual()) return true;
 
-    if (kapa == -1) {
-        kapa = KapaOpenNamedSocket ("kapa", "psphot:images");
-        if (kapa == -1) {
-            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
-            pmVisualSetVisual(false);
-            return false;
-        }
-    }
-
-    // XXX test: image->data.F32[10][10] = 10000;
+    int kapa = psphotKapaChannel (1);
+    if (kapa == -1) return false;
+
     psphotVisualRangeImage (kapa, image, "signif", 2, -1.0, 25.0*25.0);
 
-    // pause and wait for user input:
-    // continue, save (provide name), ??
-    char key[10];
-    fprintf (stdout, "[c]ontinue? ");
-    if (!fgets(key, 8, stdin)) {
-        psWarning("Unable to read option");
-    }
+    pmVisualAskUser(NULL);
     return true;
 }
@@ -224 +224 @@
     if (!pmVisualIsVisual()) return true;
 
-    if (kapa == -1) {
-        fprintf (stderr, "kapa not opened, skipping\n");
-        return false;
-    }
+    int kapa = psphotKapaChannel (1);
+    if (kapa == -1) return false;
 
     psArray *peaks = detections->peaks;
@@ -233 +231 @@
     // note: this uses the Ohana allocation tools:
     // ALLOCATE (overlay, KiiOverlay, 3*peaks->n + 1);
-    ALLOCATE (overlay, KiiOverlay, peaks->n);
+    ALLOCATE (overlay, KiiOverlay, peaks->n + 2);
 
     Noverlay = 0;
@@ -271 +269 @@
     }
 
-# if (0)
+# if (1)
     overlay[Noverlay].type = KII_OVERLAY_BOX;
     overlay[Noverlay].x = 10.0;
     overlay[Noverlay].y = 10.0;
-    overlay[Noverlay].dx = 0.5;
-    overlay[Noverlay].dy = 0.5;
+    overlay[Noverlay].dx = 1.0;
+    overlay[Noverlay].dy = 1.0;
     overlay[Noverlay].angle = 0.0;
     overlay[Noverlay].text = NULL;
@@ -285 +283 @@
     FREE (overlay);
 
-    // pause and wait for user input:
-    // continue, save (provide name), ??
-    char key[10];
-    fprintf (stdout, "[c]ontinue? ");
-    if (!fgets(key, 8, stdin)) {
-        psWarning("Unable to read option");
-    }
+    pmVisualAskUser(NULL);
     return true;
 }
@@ -302 +294 @@
     if (!pmVisualIsVisual()) return true;
 
-    if (kapa == -1) {
-        fprintf (stderr, "kapa not opened, skipping\n");
-        return false;
-    }
+    int kapa = psphotKapaChannel (1);
+    if (kapa == -1) return false;
 
     psArray *footprints = detections->footprints;
@@ -325 +315 @@
 
         // draw the top
+        // XXX need to allow top (and bottom) to have more than one span
         span = footprint->spans->data[0];
         overlay[Noverlay].type = KII_OVERLAY_LINE;
@@ -399 +390 @@
     FREE (overlay);
 
-    // pause and wait for user input:
-    // continue, save (provide name), ??
-    char key[10];
-    fprintf (stdout, "[c]ontinue? ");
-    if (!fgets(key, 8, stdin)) {
-        psWarning("Unable to read option");
-    }
+    pmVisualAskUser(NULL);
     return true;
 }
@@ -419 +404 @@
     if (!pmVisualIsVisual()) return true;
 
-    if (kapa == -1) {
-        fprintf (stderr, "kapa not opened, skipping\n");
-        return false;
-    }
+    int kapa = psphotKapaChannel (1);
+    if (kapa == -1) return false;
+
+    // XXX mark the different source classes with different color/shape dots
+    // XXX are moments S/N and peak S/N consistent?
 
     // note: this uses the Ohana allocation tools:
@@ -448 +434 @@
         overlay[Noverlay].dx = 2.0*axes.major;
         overlay[Noverlay].dy = 2.0*axes.minor;
-        overlay[Noverlay].angle = -axes.theta * PS_DEG_RAD;  // XXXXXXXX the axes angle is negative to display of object on kapa
+
+        overlay[Noverlay].angle = axes.theta * PS_DEG_RAD;
+
         overlay[Noverlay].text = NULL;
         Noverlay ++;
@@ -456 +444 @@
     FREE (overlay);
 
-    // pause and wait for user input:
-    // continue, save (provide name), ??
-    char key[10];
-    fprintf (stdout, "[c]ontinue? ");
-    if (!fgets(key, 8, stdin)) {
-        psWarning("Unable to read option");
-    }
-
+    pmVisualAskUser(NULL);
     return true;
 }
@@ -475 +456 @@
     if (!pmVisualIsVisual()) return true;
 
-    if (kapa3 == -1) {
-        kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots");
-        if (kapa3 == -1) {
-            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
-            pmVisualSetVisual(false);
-            return false;
-        }
-    }
-
-    KapaClearPlots (kapa3);
+    int myKapa = psphotKapaChannel (2);
+    if (myKapa == -1) return false;
+
+    KapaClearPlots (myKapa);
     KapaInitGraph (&graphdata);
-    KapaSetFont (kapa3, "courier", 14);
+    KapaSetFont (myKapa, "courier", 14);
 
     float SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM");
 
     // select the max psfX,Y values for the plot limits
-    float Xmin = 0.0, Xmax = 0.0;
-    float Ymin = 0.0, Ymax = 0.0;
+    float Xmin = 1000.0, Xmax = 0.0;
+    float Ymin = 1000.0, Ymax = 0.0;
     {
         int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS");
@@ -511 +486 @@
             float Y1 = psfY + 4.0*psfdY;
 
-            if (isfinite(X0)) { Xmin = PS_MAX(Xmin, X0); }
+            if (isfinite(X0)) { Xmin = PS_MIN(Xmin, X0); }
             if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); }
-            if (isfinite(Y0)) { Ymin = PS_MAX(Ymin, Y0); }
+            if (isfinite(Y0)) { Ymin = PS_MIN(Ymin, Y0); }
             if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); }
         }
     }
+    Xmin = PS_MAX(Xmin, -0.1); 
+    Ymin = PS_MAX(Ymin, -0.1); 
 
     // storage vectors for data to be plotted
@@ -564 +541 @@
     section.y  = 0.00;
     section.name = psStringCopy ("MxxMyy");
-    KapaSetSection (kapa3, &section);
+    KapaSetSection (myKapa, &section);
     psFree (section.name);
 
@@ -572 +549 @@
     graphdata.xmax = Xmax;
     graphdata.ymax = Ymax;
-    KapaSetLimits (kapa3, &graphdata);
-
-    KapaBox (kapa3, &graphdata);
-    KapaSendLabel (kapa3, "M_xx| (pixels)", KAPA_LABEL_XM);
-    KapaSendLabel (kapa3, "M_yy| (pixels)", KAPA_LABEL_YM);
+    KapaSetLimits (myKapa, &graphdata);
+
+    KapaBox (myKapa, &graphdata);
+    KapaSendLabel (myKapa, "M_xx| (pixels)", KAPA_LABEL_XM);
+    KapaSendLabel (myKapa, "M_yy| (pixels)", KAPA_LABEL_YM);
 
     graphdata.color = KapaColorByName ("black");
@@ -582 +559 @@
     graphdata.size = 0.3;
     graphdata.style = 2;
-    KapaPrepPlot (kapa3, nF, &graphdata);
-    KapaPlotVector (kapa3, nF, xFaint->data.F32, "x");
-    KapaPlotVector (kapa3, nF, yFaint->data.F32, "y");
+    KapaPrepPlot (myKapa, nF, &graphdata);
+    KapaPlotVector (myKapa, nF, xFaint->data.F32, "x");
+    KapaPlotVector (myKapa, nF, yFaint->data.F32, "y");
 
     graphdata.color = KapaColorByName ("red");
@@ -590 +567 @@
     graphdata.size = 0.5;
     graphdata.style = 2;
-    KapaPrepPlot (kapa3, nB, &graphdata);
-    KapaPlotVector (kapa3, nB, xBright->data.F32, "x");
-    KapaPlotVector (kapa3, nB, yBright->data.F32, "y");
+    KapaPrepPlot (myKapa, nB, &graphdata);
+    KapaPlotVector (myKapa, nB, xBright->data.F32, "x");
+    KapaPlotVector (myKapa, nB, yBright->data.F32, "y");
 
     // second section: MagMyy
@@ -600 +577 @@
     section.y  = 0.80;
     section.name = psStringCopy ("MagMyy");
-    KapaSetSection (kapa3, &section);
+    KapaSetSection (myKapa, &section);
     psFree (section.name);
 
@@ -608 +585 @@
     graphdata.ymin = Ymin;
     graphdata.ymax = Ymax;
-    KapaSetLimits (kapa3, &graphdata);
+    KapaSetLimits (myKapa, &graphdata);
 
     strcpy (graphdata.labels, "0210");
-    KapaBox (kapa3, &graphdata);
-    KapaSendLabel (kapa3, "inst mag", KAPA_LABEL_XP);
-    KapaSendLabel (kapa3, "M_yy| (pixels)", KAPA_LABEL_YM);
+    KapaBox (myKapa, &graphdata);
+    KapaSendLabel (myKapa, "inst mag", KAPA_LABEL_XP);
+    KapaSendLabel (myKapa, "M_yy| (pixels)", KAPA_LABEL_YM);
 
     graphdata.color = KapaColorByName ("black");
@@ -619 +596 @@
     graphdata.size = 0.3;
     graphdata.style = 2;
-    KapaPrepPlot (kapa3, nF, &graphdata);
-    KapaPlotVector (kapa3, nF, mFaint->data.F32, "x");
-    KapaPlotVector (kapa3, nF, yFaint->data.F32, "y");
+    KapaPrepPlot (myKapa, nF, &graphdata);
+    KapaPlotVector (myKapa, nF, mFaint->data.F32, "x");
+    KapaPlotVector (myKapa, nF, yFaint->data.F32, "y");
 
     graphdata.color = KapaColorByName ("red");
@@ -627 +604 @@
     graphdata.size = 0.5;
     graphdata.style = 2;
-    KapaPrepPlot (kapa3, nB, &graphdata);
-    KapaPlotVector (kapa3, nB, mBright->data.F32, "x");
-    KapaPlotVector (kapa3, nB, yBright->data.F32, "y");
+    KapaPrepPlot (myKapa, nB, &graphdata);
+    KapaPlotVector (myKapa, nB, mBright->data.F32, "x");
+    KapaPlotVector (myKapa, nB, yBright->data.F32, "y");
 
     // third section: MagMxx
@@ -637 +614 @@
     section.y  = 0.00;
     section.name = psStringCopy ("MagMxx");
-    KapaSetSection (kapa3, &section);
+    KapaSetSection (myKapa, &section);
     psFree (section.name);
 
@@ -645 +622 @@
     graphdata.ymin =  -7.9;
     graphdata.ymax = -17.1;
-    KapaSetLimits (kapa3, &graphdata);
+    KapaSetLimits (myKapa, &graphdata);
 
     strcpy (graphdata.labels, "2001");
-    KapaBox (kapa3, &graphdata);
-    KapaSendLabel (kapa3, "M_xx| (pixels)", KAPA_LABEL_XM);
-    KapaSendLabel (kapa3, "inst mag", KAPA_LABEL_YP);
+    KapaBox (myKapa, &graphdata);
+    KapaSendLabel (myKapa, "M_xx| (pixels)", KAPA_LABEL_XM);
+    KapaSendLabel (myKapa, "inst mag", KAPA_LABEL_YP);
 
     graphdata.color = KapaColorByName ("black");
@@ -656 +633 @@
     graphdata.size = 0.3;
     graphdata.style = 2;
-    KapaPrepPlot (kapa3, nF, &graphdata);
-    KapaPlotVector (kapa3, nF, xFaint->data.F32, "x");
-    KapaPlotVector (kapa3, nF, mFaint->data.F32, "y");
+    KapaPrepPlot (myKapa, nF, &graphdata);
+    KapaPlotVector (myKapa, nF, xFaint->data.F32, "x");
+    KapaPlotVector (myKapa, nF, mFaint->data.F32, "y");
 
     graphdata.color = KapaColorByName ("red");
@@ -664 +641 @@
     graphdata.size = 0.5;
     graphdata.style = 2;
-    KapaPrepPlot (kapa3, nB, &graphdata);
-    KapaPlotVector (kapa3, nB, xBright->data.F32, "x");
-    KapaPlotVector (kapa3, nB, mBright->data.F32, "y");
+    KapaPrepPlot (myKapa, nB, &graphdata);
+    KapaPlotVector (myKapa, nB, xBright->data.F32, "x");
+    KapaPlotVector (myKapa, nB, mBright->data.F32, "y");
 
     // draw N circles to outline the clumps
     {
-        KapaSelectSection (kapa3, "MxxMyy");
+        KapaSelectSection (myKapa, "MxxMyy");
 
         // draw a circle centered on psfX,Y with size of the psf limit
@@ -686 +663 @@
         graphdata.xmax = Xmax;
         graphdata.ymax = Ymax;
-        KapaSetLimits (kapa3, &graphdata);
+        KapaSetLimits (myKapa, &graphdata);
 
         for (int n = 0; n < nRegions; n++) {
@@ -705 +682 @@
                 yLimit->data.F32[i] = Ry*sin(i*2.0*M_PI/120.0) + psfY;
             }
-            KapaPrepPlot (kapa3, xLimit->n, &graphdata);
-            KapaPlotVector (kapa3, xLimit->n, xLimit->data.F32, "x");
-            KapaPlotVector (kapa3, yLimit->n, yLimit->data.F32, "y");
+            KapaPrepPlot (myKapa, xLimit->n, &graphdata);
+            KapaPlotVector (myKapa, xLimit->n, xLimit->data.F32, "x");
+            KapaPlotVector (myKapa, yLimit->n, yLimit->data.F32, "y");
         }
         psFree (xLimit);
         psFree (yLimit);
     }
-
-# if (0)
-    // *** make a histogram of the source counts in the x and y directions
-    psHistogram *nX = psHistogramAlloc (graphdata.xmin, graphdata.xmax, 50.0);
-    psHistogram *nY = psHistogramAlloc (graphdata.ymin, graphdata.ymax, 50.0);
-    psVectorHistogram (nX, xFaint, NULL, NULL, 0);
-    psVectorHistogram (nY, yFaint, NULL, NULL, 0);
-    psVector *dX = psVectorAlloc (nX->nums->n, PS_TYPE_F32);
-    psVector *vX = psVectorAlloc (nX->nums->n, PS_TYPE_F32);
-    psVector *dY = psVectorAlloc (nY->nums->n, PS_TYPE_F32);
-    psVector *vY = psVectorAlloc (nY->nums->n, PS_TYPE_F32);
-    for (int i = 0; i < nX->nums->n; i++) {
-        dX->data.F32[i] = nX->nums->data.S32[i];
-        vX->data.F32[i] = 0.5*(nX->bounds->data.F32[i] + nX->bounds->data.F32[i+1]);
-    }
-    for (int i = 0; i < nY->nums->n; i++) {
-        dY->data.F32[i] = nY->nums->data.S32[i];
-        vY->data.F32[i] = 0.5*(nY->bounds->data.F32[i] + nY->bounds->data.F32[i+1]);
-    }
-
-    graphdata.color = KapaColorByName ("black");
-    graphdata.ptype = 0;
-    graphdata.size = 0.0;
-    graphdata.style = 0;
-    KapaPrepPlot (kapa3, dX->n, &graphdata);
-    KapaPlotVector (kapa3, dX->n, dX->data.F32, "x");
-    KapaPlotVector (kapa3, vX->n, vX->data.F32, "y");
-
-    psFree (nX);
-    psFree (dX);
-    psFree (vX);
-
-    psFree (nY);
-    psFree (dY);
-    psFree (vY);
-# endif
 
     psFree (xBright);
@@ -756 +697 @@
     psFree (mFaint);
 
-    // pause and wait for user input:
-    // continue, save (provide name), ??
-    char key[10];
-    fprintf (stdout, "[c]ontinue? ");
-    if (!fgets(key, 8, stdin)) {
-        psWarning("Unable to read option");
-    }
+    pmVisualAskUser(NULL);
     return true;
 }
 
 // assumes 'kapa' value is checked and set
-bool psphotVisualShowRoughClass_Single (psArray *sources, pmSourceType type, pmSourceMode mode, char *color) {
+bool psphotVisualShowRoughClass_Single (int myKapa, psArray *sources, pmSourceType type, pmSourceMode mode, char *color) {
 
     int Noverlay;
@@ -802 +737 @@
         overlay[Noverlay].dx = 2.0*axes.major;
         overlay[Noverlay].dy = 2.0*axes.minor;
-        overlay[Noverlay].angle = -axes.theta * PS_DEG_RAD;
+        overlay[Noverlay].angle = axes.theta * PS_DEG_RAD;
         overlay[Noverlay].text = NULL;
         Noverlay ++;
     }
 
-    KiiLoadOverlay (kapa, overlay, Noverlay, color);
+    KiiLoadOverlay (myKapa, overlay, Noverlay, color);
     FREE (overlay);
 
@@ -817 +752 @@
     if (!pmVisualIsVisual()) return true;
 
-    if (kapa == -1) {
-        fprintf (stderr, "kapa not opened, skipping\n");
-        return false;
-    }
-
-    KiiEraseOverlay (kapa, "yellow"); // moments
-
-    psphotVisualShowRoughClass_Single (sources, PM_SOURCE_TYPE_STAR, 0, "red");
-    psphotVisualShowRoughClass_Single (sources, PM_SOURCE_TYPE_EXTENDED, 0, "blue");
-    psphotVisualShowRoughClass_Single (sources, PM_SOURCE_TYPE_DEFECT, 0, "blue");
-    psphotVisualShowRoughClass_Single (sources, PM_SOURCE_TYPE_SATURATED, 0, "red");
-    psphotVisualShowRoughClass_Single (sources, PM_SOURCE_TYPE_STAR, PM_SOURCE_MODE_PSFSTAR, "yellow");
-    psphotVisualShowRoughClass_Single (sources, PM_SOURCE_TYPE_STAR, PM_SOURCE_MODE_SATSTAR, "green");
-
-    // pause and wait for user input:
-    // continue, save (provide name), ??
-    char key[10];
+    int myKapa = psphotKapaChannel (1);
+    if (myKapa == -1) return false;
+
+    KiiEraseOverlay (myKapa, "yellow"); // moments
+
+    psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_STAR, 0, "red");
+    psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_EXTENDED, 0, "blue");
+    psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_DEFECT, 0, "blue");
+    psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_SATURATED, 0, "red");
+    psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_STAR, PM_SOURCE_MODE_PSFSTAR, "yellow");
+    psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_STAR, PM_SOURCE_MODE_SATSTAR, "green");
+
     fprintf (stdout, "red: STAR or SAT AREA; blue: EXTENDED or DEFECT; green: SATSTAR; yellow: PSFSTAR\n");
-    fprintf (stdout, "[c]ontinue? ");
-    if (!fgets(key, 8, stdin)) {
-        psWarning("Unable to read option");
-    }
+    pmVisualAskUser(NULL);
     return true;
 }
@@ -846 +773 @@
     if (!pmVisualIsVisual()) return true;
 
-    if (kapa2 == -1) {
-        kapa2 = KapaOpenNamedSocket ("kapa", "psphot:psfstars");
-        if (kapa2 == -1) {
-            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
-            pmVisualSetVisual(false);
-            return false;
-        }
-    }
+    int myKapa = psphotKapaChannel (3);
+    if (myKapa == -1) return false;
 
     int DX = 64;
@@ -898 +819 @@
 
     psImage *psfLogFlux = (psImage *) psUnaryOp (NULL, psfMosaic, "log");
-    psphotVisualRangeImage (kapa2, psfLogFlux, "psf_mosaic",    0, -2.0, 3.0);
-    psphotVisualRangeImage (kapa2, funMosaic, "psf_analytical", 1, -10.0, 100.0);
-    psphotVisualRangeImage (kapa2, resMosaic, "psf_residual",   2, -10.0, 100.0);
+    psphotVisualRangeImage (myKapa, psfLogFlux, "psf_mosaic",    0, -2.0, 3.0);
+    psphotVisualRangeImage (myKapa, funMosaic, "psf_analytical", 1, -10.0, 100.0);
+    psphotVisualRangeImage (myKapa, resMosaic, "psf_residual",   2, -10.0, 100.0);
 
     psFree (psfMosaic);
@@ -908 +829 @@
     psFree (modelRef);
 
-    // pause and wait for user input:
-    // continue, save (provide name), ??
-    char key[10];
-    fprintf (stdout, "[c]ontinue? ");
-    if (!fgets(key, 8, stdin)) {
-        psWarning("Unable to read option");
-    }
+    pmVisualAskUser(NULL);
     return true;
 }
@@ -924 +839 @@
     if (!pmVisualIsVisual()) return true;
 
-    if (kapa2 == -1) {
-        kapa2 = KapaOpenNamedSocket ("kapa", "psphot:psfstars");
-        if (kapa2 == -1) {
-            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
-            pmVisualSetVisual(false);
-            return false;
-        }
-    }
+    int myKapa = psphotKapaChannel (3);
+    if (myKapa == -1) return false;
 
     // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
@@ -1017 +926 @@
             if (Xo == 0) {
                 // place source alone on this row
-                psphotAddWithTest (source, true, maskVal); // replace source if subtracted
+                bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
+                if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
                 psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
 
-                psphotSubWithTest (source, false, maskVal); // remove source (force)
+                pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
                 psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
-                psphotSetState (source, false, maskVal); // reset source Add/Sub state to recorded
+
+                if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
 
                 Yo += DY;
@@ -1032 +943 @@
                 Xo = 0;
 
-                psphotAddWithTest (source, true, maskVal); // replace source if subtracted
+                bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
+                if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
                 psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
 
-                psphotSubWithTest (source, false, maskVal); // remove source (force)
+                pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
                 psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
-                psphotSetState (source, false, maskVal); // replace source (has been subtracted)
+
+                if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
 
                 Xo = DX;
@@ -1044 +957 @@
         } else {
             // extend this row
-            psphotAddWithTest (source, true, maskVal); // replace source if subtracted
+            bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
+            if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
             psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
 
-            psphotSubWithTest (source, false, maskVal); // remove source (force)
+            pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
             psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
-            psphotSetState (source, false, maskVal); // replace source (has been subtracted)
+            if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
 
             Xo += DX;
@@ -1056 +970 @@
     }
 
-    psphotVisualRangeImage (kapa2, outpos, "psfpos", 0, -0.05, 0.95);
-    psphotVisualRangeImage (kapa2, outsub, "psfsub", 1, -0.05, 0.95);
-
-    // pause and wait for user input:
-    // continue, save (provide name), ??
-    char key[10];
-    fprintf (stdout, "[c]ontinue? ");
-    if (!fgets(key, 8, stdin)) {
-        psWarning("Unable to read option");
-    }
-
+    psphotVisualRangeImage (myKapa, outpos, "psfpos", 0, -0.05, 0.95);
+    psphotVisualRangeImage (myKapa, outsub, "psfsub", 1, -0.05, 0.95);
+
+    pmVisualAskUser(NULL);
     psFree (outpos);
     psFree (outsub);
@@ -1084 +991 @@
     if (!pmVisualIsVisual()) return true;
 
-    if (kapa2 == -1) {
-        kapa2 = KapaOpenNamedSocket ("kapa", "psphot:images");
-        if (kapa2 == -1) {
-            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
-            pmVisualSetVisual(false);
-            return false;
-        }
-    }
+    int myKapa = psphotKapaChannel (3);
+    if (myKapa == -1) return false;
 
     // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
@@ -1119 +1020 @@
         pmSource *source = sources->data[i];
 
-        bool keep = false;
-        keep |= (source->mode & PM_SOURCE_MODE_SATSTAR);
-        if (!keep) continue;
+        // only show "real" saturated stars (not defects)
+        if (!(source->mode & PM_SOURCE_MODE_SATSTAR)) continue;;
+        if (source->mode & PM_SOURCE_MODE_DEFECT) continue;;
 
         // how does this subimage get placed into the output image?
@@ -1164 +1065 @@
         pmSource *source = sources->data[i];
 
-        bool keep = false;
-        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
-            nSAT ++;
-            keep = true;
-        }
-        if (!keep) continue;
+        // only show "real" saturated stars (not defects)
+        if (!(source->mode & PM_SOURCE_MODE_SATSTAR)) continue;;
+        if (source->mode & PM_SOURCE_MODE_DEFECT) continue;;
+        nSAT ++;
 
         if (Xo + DX > NX) {
@@ -1175 +1074 @@
             if (Xo == 0) {
                 // place source alone on this row
-                psphotAddWithTest (source, true, maskVal); // replace source if subtracted
+                bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
+                if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
                 psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
-                psphotSetState (source, true, maskVal); // reset source Add/Sub state to recorded
+                if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
 
                 Yo += DY;
@@ -1186 +1086 @@
                 Yo += dY;
                 Xo = 0;
-                psphotAddWithTest (source, true, maskVal); // replace source if subtracted
+
+                bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
+                if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
                 psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
-                psphotSetState (source, true, maskVal); // replace source (has been subtracted)
+                if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
 
                 Xo = DX;
@@ -1195 +1097 @@
         } else {
             // extend this row
-            psphotAddWithTest (source, true, maskVal); // replace source if subtracted
+            bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
+            if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
             psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
-            psphotSetState (source, true, maskVal); // replace source (has been subtracted)
+            if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
 
             Xo += DX;
@@ -1204 +1107 @@
     }
 
-    psphotVisualScaleImage (kapa2, outsat, "satstar", 2);
-
-    // pause and wait for user input:
-    // continue, save (provide name), ??
-    char key[10];
-    fprintf (stdout, "[c]ontinue? ");
-    if (!fgets(key, 8, stdin)) {
-        psWarning("Unable to read option");
-    }
-
+    psphotVisualScaleImage (myKapa, outsat, "satstar", 2);
+
+    pmVisualAskUser(NULL);
     psFree (outsat);
     return true;
@@ -1222 +1118 @@
     Graphdata graphdata;
 
-    bool state = !(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
-    psphotAddWithTest (source, true, maskVal); // replace source if subtracted
+    bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
+    if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
 
     int nPts = source->pixels->numRows * source->pixels->numCols;
@@ -1240 +1136 @@
         for (int ix = 0; ix < source->pixels->numCols; ix++) {
             if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) {
-                // rb->data.F32[nb] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ;
-                rb->data.F32[nb] = hypot (ix - Xo, iy - Yo) ;
+                rb->data.F32[nb] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ;
+                // rb->data.F32[nb] = hypot (ix - Xo, iy - Yo) ;
                 Rb->data.F32[nb] = log10(rb->data.F32[nb]);
                 fb->data.F32[nb] = log10(source->pixels->data.F32[iy][ix]);
                 nb++;
             } else {
-                // rg->data.F32[ng] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ;
-                rg->data.F32[ng] = hypot (ix - Xo, iy - Yo) ;
+                rg->data.F32[ng] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ;
+                // rg->data.F32[ng] = hypot (ix - Xo, iy - Yo) ;
                 Rg->data.F32[ng] = log10(rg->data.F32[ng]);
                 fg->data.F32[ng] = log10(source->pixels->data.F32[iy][ix]);
@@ -1256 +1152 @@
 
     // reset source Add/Sub state to recorded
-    psphotSetState (source, state, maskVal);
+    if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
 
     KapaInitGraph (&graphdata);
@@ -1337 +1233 @@
     if (!pmVisualIsVisual()) return true;
 
-    if (kapa3 == -1) {
-        kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots");
-        if (kapa3 == -1) {
-            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
-            pmVisualSetVisual(false);
-            return false;
-        }
-    }
+    int myKapa = psphotKapaChannel (2);
+    if (myKapa == -1) return false;
 
     // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
@@ -1351 +1241 @@
     assert (maskVal);
 
-    KapaClearPlots (kapa3);
+    KapaClearPlots (myKapa);
     // first section : mag vs CR nSigma
     section.dx = 1.0;
@@ -1359 +1249 @@
     section.name = NULL;
     psStringAppend (&section.name, "linlog");
-    KapaSetSection (kapa3, &section);
+    KapaSetSection (myKapa, &section);
     psFree (section.name);
 
@@ -1369 +1259 @@
     section.name = NULL;
     psStringAppend (&section.name, "loglog");
-    KapaSetSection (kapa3, &section);
+    KapaSetSection (myKapa, &section);
     psFree (section.name);
 
@@ -1378 +1268 @@
         if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
 
-        psphotVisualPlotRadialProfile (kapa3, source, maskVal);
+        psphotVisualPlotRadialProfile (myKapa, source, maskVal);
 
         // pause and wait for user input:
@@ -1388 +1278 @@
         }
         if (key[0] == 'e') {
-            KapaClearPlots (kapa3);
+            KapaClearPlots (myKapa);
         }
         if (key[0] == 's') {
@@ -1407 +1297 @@
     psEllipseAxes axes;
 
+    // XXX skip this for now: it is not very clear
+    return true;
+
     if (!pmVisualIsVisual()) return true;
 
-    if (kapa == -1) {
-        fprintf (stderr, "kapa not opened, skipping\n");
-        return false;
-    }
+    int myKapa = psphotKapaChannel (1);
+    if (myKapa == -1) return false;
 
     // note: this uses the Ohana allocation tools:
@@ -1485 +1376 @@
     }
 
-    KiiLoadOverlay (kapa, overlayE, NoverlayE, "red");
-    KiiLoadOverlay (kapa, overlayO, NoverlayO, "yellow");
+    KiiLoadOverlay (myKapa, overlayE, NoverlayE, "red");
+    KiiLoadOverlay (myKapa, overlayO, NoverlayO, "yellow");
     FREE (overlayE);
     FREE (overlayO);
 
-    // pause and wait for user input:
-    // continue, save (provide name), ??
-    char key[10];
     fprintf (stdout, "even bits (0x0001, 0x0004, ... : red\n");
     fprintf (stdout, "odd bits (0x0002, 0x0008, ... : yellow\n");
-    fprintf (stdout, "[c]ontinue? ");
-    if (!fgets(key, 8, stdin)) {
-        psWarning("Unable to read option");
-    }
-
-    return true;
-}
-
-bool psphotVisualShowSourceSize (pmReadout *readout, psArray *sources) {
-
-    int Noverlay, NOVERLAY;
+    pmVisualAskUser(NULL);
+
+    return true;
+}
+
+bool psphotVisualShowSourceSize_Single (int myKapa, psArray *sources, pmSourceMode mode, bool keep, float scale, char *color) {
+
+    int Noverlay;
     KiiOverlay *overlay;
 
-    if (!pmVisualIsVisual()) return true;
-
-    if (kapa == -1) {
-        fprintf (stderr, "kapa not opened, skipping\n");
-        return false;
-    }
+    psEllipseMoments emoments;
+    psEllipseAxes axes;
 
     // note: this uses the Ohana allocation tools:
+    ALLOCATE (overlay, KiiOverlay, sources->n);
+
     Noverlay = 0;
-    NOVERLAY = 100;
-    ALLOCATE (overlay, KiiOverlay, sources->n);
-
-    // mark CRs with red boxes
     for (int i = 0; i < sources->n; i++) {
 
@@ -1526 +1405 @@
         if (source == NULL) continue;
 
-        if (!(source->mode & PM_SOURCE_MODE_CR_LIMIT)) continue;
-
-        overlay[Noverlay].type = KII_OVERLAY_BOX;
-        overlay[Noverlay].x = source->peak->xf;
-        overlay[Noverlay].y = source->peak->yf;
-
-        overlay[Noverlay].dx = 4;
-        overlay[Noverlay].dy = 4;
-        overlay[Noverlay].angle = 0;
+        // if (source->type != type) continue;
+        if (mode) {
+            if (keep) {
+                if (!(source->mode & mode)) continue;
+            } else {
+                if (source->mode & mode) continue;
+            }
+        }
+
+        pmMoments *moments = source->moments;
+        if (moments == NULL) continue;
+
+        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
+        overlay[Noverlay].x = moments->Mx;
+        overlay[Noverlay].y = moments->My;
+
+        emoments.x2 = moments->Mxx;
+        emoments.y2 = moments->Myy;
+        emoments.xy = moments->Mxy;
+
+        axes = psEllipseMomentsToAxes (emoments, 20.0);
+
+        overlay[Noverlay].dx = scale*2.0*axes.major;
+        overlay[Noverlay].dy = scale*2.0*axes.minor;
+        overlay[Noverlay].angle = axes.theta * PS_DEG_RAD;
         overlay[Noverlay].text = NULL;
         Noverlay ++;
-        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
-    }
-    KiiLoadOverlay (kapa, overlay, Noverlay, "red");
-
-
-    Noverlay = 0;
-    for (int i = 0; i < sources->n; i++) {
-
-        pmSource *source = sources->data[i];
-        if (source == NULL) continue;
-
-        // mark EXTs with yellow circles
-        if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;
-
-        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
-        overlay[Noverlay].x = source->peak->xf;
-        overlay[Noverlay].y = source->peak->yf;
-
-        overlay[Noverlay].dx = 10;
-        overlay[Noverlay].dy = 10;
-        overlay[Noverlay].angle = 0;
-        overlay[Noverlay].text = NULL;
-        Noverlay ++;
-        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
-    }
-
-    KiiLoadOverlay (kapa, overlay, Noverlay, "blue");
+    }
+
+    KiiLoadOverlay (myKapa, overlay, Noverlay, color);
     FREE (overlay);
 
-    psphotVisualShowMask (kapa, readout->mask, "mask", 2);
-
-    // pause and wait for user input:
-    // continue, save (provide name), ??
-    char key[10];
-    fprintf (stdout, "CR: 4pix red BOX; EXT: 10pix blue circle\n");
-    fprintf (stdout, "[c]ontinue? ");
-    if (!fgets(key, 8, stdin)) {
-        psWarning("Unable to read option");
-    }
-
-    return true;
-}
-
-bool psphotVisualPlotSourceSize (psArray *sources) {
-
+    return true;
+}
+
+bool psphotVisualShowSourceSize (pmReadout *readout, psArray *sources) {
+
+    if (!pmVisualIsVisual()) return true;
+
+    int myKapa = psphotKapaChannel (1);
+    if (myKapa == -1) return false;
+
+    KiiEraseOverlay (myKapa, "red");
+    KiiEraseOverlay (myKapa, "blue");
+    KiiEraseOverlay (myKapa, "green");
+    KiiEraseOverlay (myKapa, "yellow");
+
+    psphotVisualShowSourceSize_Single (myKapa, sources, PM_SOURCE_MODE_EXT_LIMIT | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_SATSTAR, 0, 1.0, "green");
+    psphotVisualShowSourceSize_Single (myKapa, sources, PM_SOURCE_MODE_EXT_LIMIT, 1, 1.0, "blue");
+    psphotVisualShowSourceSize_Single (myKapa, sources, PM_SOURCE_MODE_CR_LIMIT, 1, 1.0, "red");
+    psphotVisualShowSourceSize_Single (myKapa, sources, PM_SOURCE_MODE_DEFECT, 1, 2.0, "red");
+    psphotVisualShowSourceSize_Single (myKapa, sources, PM_SOURCE_MODE_SATSTAR, 1, 1.0, "yellow");
+
+    fprintf (stdout, "red: CR; blue: EXTENDED; green: PSF-like; yellow: SATSTAR\n");
+    pmVisualAskUser(NULL);
+    return true;
+}
+
+bool psphotVisualPlotSourceSize (psMetadata *recipe, psArray *sources) {
+
+    bool status;
     Graphdata graphdata;
     KapaSection section;
@@ -1587 +1471 @@
     if (!pmVisualIsVisual()) return true;
 
-    if (kapa3 == -1) {
-        kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots");
-        if (kapa3 == -1) {
-            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
-            pmVisualSetVisual(false);
-            return false;
+    int myKapa = psphotKapaChannel (2);
+    if (myKapa == -1) return false;
+
+    KapaClearPlots (myKapa);
+    KapaInitGraph (&graphdata);
+    KapaSetFont (myKapa, "courier", 14);
+
+    // select the max psfX,Y values for the plot limits
+    float Xmin = 1000.0, Xmax = 0.0;
+    float Ymin = 1000.0, Ymax = 0.0;
+    {
+        int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS");
+        for (int n = 0; n < nRegions; n++) {
+
+            char regionName[64];
+            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
+            psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
+
+            float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
+            float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
+            float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");
+            float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");
+
+            float X0 = psfX - 10.0*psfdX;
+            float X1 = psfX + 10.0*psfdX;
+            float Y0 = psfY - 10.0*psfdY;
+            float Y1 = psfY + 10.0*psfdY;
+
+            if (isfinite(X0)) { Xmin = PS_MIN(Xmin, X0); }
+            if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); }
+            if (isfinite(Y0)) { Ymin = PS_MIN(Ymin, Y0); }
+            if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); }
         }
     }
-
-    KapaClearPlots (kapa3);
+    Xmin = PS_MAX(Xmin, -0.1); 
+    Ymin = PS_MAX(Ymin, -0.1); 
+
+    // storage vectors for data to be plotted
+    psVector *xSAT = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *ySAT = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *mSAT = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *sSAT = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+
+    psVector *xPSF = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *yPSF = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *mPSF = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *sPSF = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+
+    psVector *xEXT = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *yEXT = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *mEXT = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *sEXT = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+
+    psVector *xDEF = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *yDEF = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *mDEF = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *sDEF = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+
+    psVector *xCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *yCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *mCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *sCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+
+    // construct the vectors
+    int nSAT = 0;
+    int nEXT = 0;
+    int nPSF = 0;
+    int nDEF = 0;
+    int nCR  = 0;
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = sources->data[i];
+        if (source->moments == NULL) continue;
+
+        if (source->mode & PM_SOURCE_MODE_CR_LIMIT) { 
+            xCR->data.F32[nCR] = source->moments->Mxx;
+            yCR->data.F32[nCR] = source->moments->Myy;
+            mCR->data.F32[nCR] = -2.5*log10(source->moments->Sum);
+            sCR->data.F32[nCR] = source->extNsigma;
+            nCR++;
+        }
+        if (source->mode & PM_SOURCE_MODE_SATSTAR) { 
+            xSAT->data.F32[nSAT] = source->moments->Mxx;
+            ySAT->data.F32[nSAT] = source->moments->Myy;
+            mSAT->data.F32[nSAT] = -2.5*log10(source->moments->Sum);
+            sSAT->data.F32[nSAT] = source->extNsigma;
+            nSAT++;
+        }
+        if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { 
+            xEXT->data.F32[nEXT] = source->moments->Mxx;
+            yEXT->data.F32[nEXT] = source->moments->Myy;
+            mEXT->data.F32[nEXT] = -2.5*log10(source->moments->Sum);
+            sEXT->data.F32[nEXT] = source->extNsigma;
+            nEXT++;
+            continue;
+        }
+        if (source->mode & PM_SOURCE_MODE_DEFECT) { 
+            xDEF->data.F32[nDEF] = source->moments->Mxx;
+            yDEF->data.F32[nDEF] = source->moments->Myy;
+            mDEF->data.F32[nDEF] = -2.5*log10(source->moments->Sum);
+            sDEF->data.F32[nDEF] = source->extNsigma;
+            nDEF++;
+            continue;
+        }
+        if ((source->mode & PM_SOURCE_MODE_CR_LIMIT) || (source->mode & PM_SOURCE_MODE_SATSTAR)) {
+            continue;
+        }
+        xPSF->data.F32[nPSF] = source->moments->Mxx;
+        yPSF->data.F32[nPSF] = source->moments->Myy;
+        mPSF->data.F32[nPSF] = -2.5*log10(source->moments->Sum);
+        sPSF->data.F32[nPSF] = source->extNsigma;
+        nPSF++;
+    }
+    xSAT->n = nSAT;
+    ySAT->n = nSAT;
+    mSAT->n = nSAT;
+    sSAT->n = nSAT;
+
+    xPSF->n = nPSF;
+    yPSF->n = nPSF;
+    mPSF->n = nPSF;
+    sPSF->n = nPSF;
+
+    xEXT->n = nEXT;
+    yEXT->n = nEXT;
+    mEXT->n = nEXT;
+    sEXT->n = nEXT;
+
+    xCR->n = nCR;
+    yCR->n = nCR;
+    mCR->n = nCR;
+    sCR->n = nCR;
+
+    xDEF->n = nDEF;
+    yDEF->n = nDEF;
+    mDEF->n = nDEF;
+    sDEF->n = nDEF;
+
+    // four sections: MxxMyy, MagMxx, MagMyy, MagSigma
+
+    // first section: MxxMyy
+    section.dx = 0.75;
+    section.dy = 0.60;
+    section.x  = 0.00;
+    section.y  = 0.00;
+    section.name = psStringCopy ("MxxMyy");
+    KapaSetSection (myKapa, &section);
+    psFree (section.name);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.xmin = Xmin;
+    graphdata.ymin = Ymin;
+    graphdata.xmax = Xmax;
+    graphdata.ymax = Ymax;
+    KapaSetLimits (myKapa, &graphdata);
+
+    KapaBox (myKapa, &graphdata);
+    KapaSendLabel (myKapa, "M_xx| (pixels)", KAPA_LABEL_XM);
+    KapaSendLabel (myKapa, "M_yy| (pixels)", KAPA_LABEL_YM);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.ptype = 0;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+    KapaPrepPlot   (myKapa, nPSF, &graphdata);
+    KapaPlotVector (myKapa, nPSF, xPSF->data.F32, "x");
+    KapaPlotVector (myKapa, nPSF, yPSF->data.F32, "y");
+
+    graphdata.color = KapaColorByName ("blue");
+    graphdata.ptype = 0;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+    KapaPrepPlot   (myKapa, nEXT, &graphdata);
+    KapaPlotVector (myKapa, nEXT, xEXT->data.F32, "x");
+    KapaPlotVector (myKapa, nEXT, yEXT->data.F32, "y");
+
+    graphdata.color = KapaColorByName ("red");
+    graphdata.ptype = 0;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+    KapaPrepPlot   (myKapa, nDEF, &graphdata);
+    KapaPlotVector (myKapa, nDEF, xDEF->data.F32, "x");
+    KapaPlotVector (myKapa, nDEF, yDEF->data.F32, "y");
+
+    graphdata.color = KapaColorByName ("red");
+    graphdata.ptype = 7;
+    graphdata.size = 1.0;
+    graphdata.style = 2;
+    KapaPrepPlot   (myKapa, nCR, &graphdata);
+    KapaPlotVector (myKapa, nCR, xCR->data.F32, "x");
+    KapaPlotVector (myKapa, nCR, yCR->data.F32, "y");
+
+    graphdata.color = KapaColorByName ("blue");
+    graphdata.ptype = 7;
+    graphdata.size = 1.0;
+    graphdata.style = 2;
+    KapaPrepPlot   (myKapa, nSAT, &graphdata);
+    KapaPlotVector (myKapa, nSAT, xSAT->data.F32, "x");
+    KapaPlotVector (myKapa, nSAT, ySAT->data.F32, "y");
+
+    // second section: MagMyy
+    section.dx = 0.75;
+    section.dy = 0.20;
+    section.x  = 0.00;
+    section.y  = 0.80;
+    section.name = psStringCopy ("MagMyy");
+    KapaSetSection (myKapa, &section);
+    psFree (section.name);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.xmin = -17.1;
+    graphdata.xmax =  -7.9;
+    graphdata.ymin = Ymin;
+    graphdata.ymax = Ymax;
+    KapaSetLimits (myKapa, &graphdata);
+
+    strcpy (graphdata.labels, "0210");
+    KapaBox (myKapa, &graphdata);
+    KapaSendLabel (myKapa, "inst mag", KAPA_LABEL_XP);
+    KapaSendLabel (myKapa, "M_yy| (pixels)", KAPA_LABEL_YM);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.ptype = 0;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+    KapaPrepPlot   (myKapa, nPSF, &graphdata);
+    KapaPlotVector (myKapa, nPSF, mPSF->data.F32, "x");
+    KapaPlotVector (myKapa, nPSF, yPSF->data.F32, "y");
+
+    graphdata.color = KapaColorByName ("blue");
+    graphdata.ptype = 0;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+    KapaPrepPlot   (myKapa, nEXT, &graphdata);
+    KapaPlotVector (myKapa, nEXT, mEXT->data.F32, "x");
+    KapaPlotVector (myKapa, nEXT, yEXT->data.F32, "y");
+
+    graphdata.color = KapaColorByName ("red");
+    graphdata.ptype = 0;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+    KapaPrepPlot   (myKapa, nDEF, &graphdata);
+    KapaPlotVector (myKapa, nDEF, mDEF->data.F32, "x");
+    KapaPlotVector (myKapa, nDEF, yDEF->data.F32, "y");
+
+    graphdata.color = KapaColorByName ("red");
+    graphdata.ptype = 7;
+    graphdata.size = 1.0;
+    graphdata.style = 2;
+    KapaPrepPlot   (myKapa, nCR, &graphdata);
+    KapaPlotVector (myKapa, nCR, mCR->data.F32, "x");
+    KapaPlotVector (myKapa, nCR, yCR->data.F32, "y");
+
+    graphdata.color = KapaColorByName ("blue");
+    graphdata.ptype = 7;
+    graphdata.size = 1.0;
+    graphdata.style = 2;
+    KapaPrepPlot   (myKapa, nSAT, &graphdata);
+    KapaPlotVector (myKapa, nSAT, mSAT->data.F32, "x");
+    KapaPlotVector (myKapa, nSAT, ySAT->data.F32, "y");
+
+    // third section: MagMxx
+    section.dx = 0.25;
+    section.dy = 0.60;
+    section.x  = 0.80;
+    section.y  = 0.00;
+    section.name = psStringCopy ("MagMxx");
+    KapaSetSection (myKapa, &section);
+    psFree (section.name);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.xmin = Xmin;
+    graphdata.xmax = Xmax;
+    graphdata.ymin =  -7.9;
+    graphdata.ymax = -17.1;
+    KapaSetLimits (myKapa, &graphdata);
+
+    strcpy (graphdata.labels, "2001");
+    KapaBox (myKapa, &graphdata);
+    KapaSendLabel (myKapa, "M_xx| (pixels)", KAPA_LABEL_XM);
+    KapaSendLabel (myKapa, "inst mag", KAPA_LABEL_YP);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.ptype = 0;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+    KapaPrepPlot   (myKapa, nPSF, &graphdata);
+    KapaPlotVector (myKapa, nPSF, xPSF->data.F32, "x");
+    KapaPlotVector (myKapa, nPSF, mPSF->data.F32, "y");
+
+    graphdata.color = KapaColorByName ("blue");
+    graphdata.ptype = 0;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+    KapaPrepPlot   (myKapa, nEXT, &graphdata);
+    KapaPlotVector (myKapa, nEXT, xEXT->data.F32, "x");
+    KapaPlotVector (myKapa, nEXT, mEXT->data.F32, "y");
+
+    graphdata.color = KapaColorByName ("red");
+    graphdata.ptype = 0;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+    KapaPrepPlot   (myKapa, nDEF, &graphdata);
+    KapaPlotVector (myKapa, nDEF, xDEF->data.F32, "x");
+    KapaPlotVector (myKapa, nDEF, mDEF->data.F32, "y");
+
+    graphdata.color = KapaColorByName ("red");
+    graphdata.ptype = 7;
+    graphdata.size = 1.0;
+    graphdata.style = 2;
+    KapaPrepPlot   (myKapa, nCR, &graphdata);
+    KapaPlotVector (myKapa, nCR, xCR->data.F32, "x");
+    KapaPlotVector (myKapa, nCR, mCR->data.F32, "y");
+
+    graphdata.color = KapaColorByName ("blue");
+    graphdata.ptype = 7;
+    graphdata.size = 1.0;
+    graphdata.style = 2;
+    KapaPrepPlot   (myKapa, nSAT, &graphdata);
+    KapaPlotVector (myKapa, nSAT, xSAT->data.F32, "x");
+    KapaPlotVector (myKapa, nSAT, mSAT->data.F32, "y");
+
+    // fourth section: MagSigma
+    section.dx = 0.75;
+    section.dy = 0.15;
+    section.x  = 0.00;
+    section.y  = 0.65;
+    section.name = psStringCopy ("MagSigma");
+    KapaSetSection (myKapa, &section);
+    psFree (section.name);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.xmax =  -7.9;
+    graphdata.xmin = -17.1;
+    graphdata.ymin = -20.1;
+    graphdata.ymax = +20.1;
+    KapaSetLimits (myKapa, &graphdata);
+
+    strcpy (graphdata.labels, "0100");
+    KapaBox (myKapa, &graphdata);
+    // KapaSendLabel (myKapa, "inst mag", KAPA_LABEL_XM);
+    KapaSendLabel (myKapa, "EXT&ss&c", KAPA_LABEL_YM);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.ptype = 0;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+    KapaPrepPlot   (myKapa, nPSF, &graphdata);
+    KapaPlotVector (myKapa, nPSF, mPSF->data.F32, "x");
+    KapaPlotVector (myKapa, nPSF, sPSF->data.F32, "y");
+
+    graphdata.color = KapaColorByName ("blue");
+    graphdata.ptype = 0;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+    KapaPrepPlot   (myKapa, nEXT, &graphdata);
+    KapaPlotVector (myKapa, nEXT, mEXT->data.F32, "x");
+    KapaPlotVector (myKapa, nEXT, sEXT->data.F32, "y");
+
+    graphdata.color = KapaColorByName ("red");
+    graphdata.ptype = 0;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+    KapaPrepPlot   (myKapa, nDEF, &graphdata);
+    KapaPlotVector (myKapa, nDEF, mDEF->data.F32, "x");
+    KapaPlotVector (myKapa, nDEF, sDEF->data.F32, "y");
+
+    graphdata.color = KapaColorByName ("red");
+    graphdata.ptype = 7;
+    graphdata.size = 1.0;
+    graphdata.style = 2;
+    KapaPrepPlot   (myKapa, nCR, &graphdata);
+    KapaPlotVector (myKapa, nCR, mCR->data.F32, "x");
+    KapaPlotVector (myKapa, nCR, sCR->data.F32, "y");
+
+    graphdata.color = KapaColorByName ("blue");
+    graphdata.ptype = 7;
+    graphdata.size = 1.0;
+    graphdata.style = 2;
+    KapaPrepPlot   (myKapa, nSAT, &graphdata);
+    KapaPlotVector (myKapa, nSAT, mSAT->data.F32, "x");
+    KapaPlotVector (myKapa, nSAT, sSAT->data.F32, "y");
+
+    // draw N circles to outline the clumps
+    {
+        KapaSelectSection (myKapa, "MxxMyy");
+
+        // draw a circle centered on psfX,Y with size of the psf limit
+        psVector *xLimit  = psVectorAlloc (120, PS_TYPE_F32);
+        psVector *yLimit  = psVectorAlloc (120, PS_TYPE_F32);
+
+        int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS");
+        float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA");
+
+        graphdata.color = KapaColorByName ("blue");
+        graphdata.style = 0;
+
+        graphdata.xmin = Xmin;
+        graphdata.ymin = Ymin;
+        graphdata.xmax = Xmax;
+        graphdata.ymax = Ymax;
+        KapaSetLimits (myKapa, &graphdata);
+
+        for (int n = 0; n < nRegions; n++) {
+
+            char regionName[64];
+            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
+            psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
+
+            float psfX  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
+            float psfY  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
+            float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");
+            float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");
+            float Rx = psfdX * PSF_CLUMP_NSIGMA;
+            float Ry = psfdY * PSF_CLUMP_NSIGMA;
+
+            for (int i = 0; i < xLimit->n; i++) {
+                xLimit->data.F32[i] = Rx*cos(i*2.0*M_PI/120.0) + psfX;
+                yLimit->data.F32[i] = Ry*sin(i*2.0*M_PI/120.0) + psfY;
+            }
+            KapaPrepPlot (myKapa, xLimit->n, &graphdata);
+            KapaPlotVector (myKapa, xLimit->n, xLimit->data.F32, "x");
+            KapaPlotVector (myKapa, yLimit->n, yLimit->data.F32, "y");
+        }
+        psFree (xLimit);
+        psFree (yLimit);
+    }
+
+    psFree (xSAT);
+    psFree (ySAT);
+    psFree (mSAT);
+
+    psFree (xEXT);
+    psFree (yEXT);
+    psFree (mEXT);
+
+    psFree (xPSF);
+    psFree (yPSF);
+    psFree (mPSF);
+
+    psFree (xDEF);
+    psFree (yDEF);
+    psFree (mDEF);
+
+    psFree (xCR);
+    psFree (yCR);
+    psFree (mCR);
+
+    pmVisualAskUser(NULL);
+    return true;
+}
+
+bool psphotVisualShowResidualImage (pmReadout *readout) {
+
+    if (!pmVisualIsVisual()) return true;
+
+    int myKapa = psphotKapaChannel (1);
+    if (myKapa == -1) return false;
+
+    psphotVisualScaleImage (myKapa, readout->image, "resid", 1);
+
+    pmVisualAskUser(NULL);
+    return true;
+}
+
+bool psphotVisualPlotApResid (psArray *sources, float mean, float error) {
+
+    Graphdata graphdata;
+    float lineX[2], lineY[2];
+
+    if (!pmVisualIsVisual()) return true;
+
+    int myKapa = psphotKapaChannel (2);
+    if (myKapa == -1) return false;
+
+    KapaClearPlots (myKapa);
     KapaInitGraph (&graphdata);
-
-    // first section : mag vs CR nSigma
-    section.dx = 1.0;
-    section.dy = 0.5;
-    section.x = 0.0;
-    section.y = 0.0;
-    section.name = NULL;
-    psStringAppend (&section.name, "a1");
-    KapaSetSection (kapa3, &section);
-    psFree (section.name);
 
     psVector *x = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     psVector *y = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
-
-    graphdata.xmin = +32.0;
-    graphdata.xmax = -32.0;
-    graphdata.ymin = +32.0;
-    graphdata.ymax = -32.0;
-
-    // construct the plot vectors
-    int n = 0;
-    for (int i = 0; i < sources->n; i++) {
-        pmSource *source = sources->data[i];
-        if (!source) continue;
-        if (source->type != PM_SOURCE_TYPE_STAR) continue;
-        if (!isfinite (source->crNsigma)) continue;
-
-        x->data.F32[n] = -2.5*log10(source->peak->flux);
-        y->data.F32[n] = source->crNsigma;
-        graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]);
-        graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]);
-        graphdata.ymin = -0.5;
-        graphdata.ymax = 10.0;
-
-        n++;
-    }
-    x->n = y->n = n;
-
-    float range;
-    range = graphdata.xmax - graphdata.xmin;
-    graphdata.xmax += 0.05*range;
-    graphdata.xmin -= 0.05*range;
-
-    // XXX set the plot range to match the image
-    KapaSetLimits (kapa3, &graphdata);
-
-    KapaSetFont (kapa3, "helvetica", 14);
-    KapaBox (kapa3, &graphdata);
-    KapaSendLabel (kapa3, "Peak as Mag", KAPA_LABEL_XM);
-    KapaSendLabel (kapa3, "CR N Sigma", KAPA_LABEL_YM);
-
-    graphdata.color = KapaColorByName ("black");
-    graphdata.ptype = 2;
-    graphdata.size = 0.5;
-    graphdata.style = 2;
-    KapaPrepPlot (kapa3, n, &graphdata);
-    KapaPlotVector (kapa3, n, x->data.F32, "x");
-    KapaPlotVector (kapa3, n, y->data.F32, "y");
-
-    // second section : mag vs EXT nSigma
-    section.dx = 1.0;
-    section.dy = 0.5;
-    section.x = 0.0;
-    section.y = 0.5;
-    section.name = NULL;
-    psStringAppend (&section.name, "a2");
-    KapaSetSection (kapa3, &section);
-    psFree (section.name);
-
-    graphdata.xmin = +32.0;
-    graphdata.xmax = -32.0;
-    graphdata.ymin = +32.0;
-    graphdata.ymax = -32.0;
-
-    // construct the plot vectors
-    n = 0;
-    for (int i = 0; i < sources->n; i++) {
-        pmSource *source = sources->data[i];
-        if (!source) continue;
-        if (source->type != PM_SOURCE_TYPE_STAR) continue;
-        if (!isfinite (source->extNsigma)) continue;
-
-        x->data.F32[n] = -2.5*log10(source->peak->flux);
-        y->data.F32[n] = source->extNsigma;
-        graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]);
-        graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]);
-        graphdata.ymin = -0.5;
-        graphdata.ymax = 10.0;
-
-        n++;
-    }
-    x->n = y->n = n;
-
-    range = graphdata.xmax - graphdata.xmin;
-    graphdata.xmax += 0.05*range;
-    graphdata.xmin -= 0.05*range;
-
-    // XXX set the plot range to match the image
-    KapaSetLimits (kapa3, &graphdata);
-
-    KapaSetFont (kapa3, "helvetica", 14);
-    KapaBox (kapa3, &graphdata);
-    KapaSendLabel (kapa3, "EXT N Sigma", KAPA_LABEL_YM);
-
-    graphdata.color = KapaColorByName ("black");
-    graphdata.ptype = 2;
-    graphdata.size = 0.5;
-    graphdata.style = 2;
-    KapaPrepPlot (kapa3, n, &graphdata);
-    KapaPlotVector (kapa3, n, x->data.F32, "x");
-    KapaPlotVector (kapa3, n, y->data.F32, "y");
-
-    psFree (x);
-    psFree (y);
-
-    // pause and wait for user input:
-    // continue, save (provide name), ??
-    char key[10];
-    fprintf (stdout, "[c]ontinue? ");
-    if (!fgets(key, 8, stdin)) {
-        psWarning("Unable to read option");
-    }
-    return true;
-}
-
-bool psphotVisualShowResidualImage (pmReadout *readout) {
-
-    if (!pmVisualIsVisual()) return true;
-
-    if (kapa == -1) {
-        kapa = KapaOpenNamedSocket ("kapa", "psphot:images");
-        if (kapa == -1) {
-            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
-            pmVisualSetVisual(false);
-            return false;
-        }
-    }
-
-    psphotVisualScaleImage (kapa, readout->image, "resid", 1);
-
-    // pause and wait for user input:
-    // continue, save (provide name), ??
-    char key[10];
-    fprintf (stdout, "[c]ontinue? ");
-    if (!fgets(key, 8, stdin)) {
-        psWarning("Unable to read option");
-    }
-    return true;
-}
-
-bool psphotVisualPlotApResid (psArray *sources) {
-
-    Graphdata graphdata;
-
-    if (!pmVisualIsVisual()) return true;
-
-    if (kapa3 == -1) {
-        kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots");
-        if (kapa3 == -1) {
-            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
-            pmVisualSetVisual(false);
-            return false;
-        }
-    }
-
-    KapaClearPlots (kapa3);
-    KapaInitGraph (&graphdata);
-
-    psVector *x = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
-    psVector *y = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *dy = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
 
     graphdata.xmin = +32.0;
@@ -1785 +1968 @@
         x->data.F32[n] = source->psfMag;
         y->data.F32[n] = source->apMag - source->psfMag;
+        dy->data.F32[n] = source->errMag;
         graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]);
         graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]);
@@ -1792 +1976 @@
         n++;
     }
-    x->n = y->n = n;
+    x->n = y->n = dy->n = n;
 
     float range;
@@ -1802 +1986 @@
     graphdata.ymin -= 0.05*range;
 
-    // XXX set the plot range to match the image
-    KapaSetLimits (kapa3, &graphdata);
-
-    KapaSetFont (kapa3, "helvetica", 14);
-    KapaBox (kapa3, &graphdata);
-    KapaSendLabel (kapa3, "PSF Mag", KAPA_LABEL_XM);
-    KapaSendLabel (kapa3, "Ap Mag - PSF Mag", KAPA_LABEL_YM);
+    // XXX test
+    graphdata.xmin = -17.0;
+    graphdata.xmax =  -9.0;
+    graphdata.ymin = -0.31;
+    graphdata.ymax = +0.31;
+
+    KapaSetLimits (myKapa, &graphdata);
+
+    KapaSetFont (myKapa, "helvetica", 14);
+    KapaBox (myKapa, &graphdata);
+    KapaSendLabel (myKapa, "PSF Mag", KAPA_LABEL_XM);
+    KapaSendLabel (myKapa, "Ap Mag - PSF Mag", KAPA_LABEL_YM);
 
     graphdata.color = KapaColorByName ("black");
@@ -1814 +2003 @@
     graphdata.size = 0.5;
     graphdata.style = 2;
-    KapaPrepPlot (kapa3, n, &graphdata);
-    KapaPlotVector (kapa3, n, x->data.F32, "x");
-    KapaPlotVector (kapa3, n, y->data.F32, "y");
+    graphdata.etype |= 0x01;
+    KapaPrepPlot (myKapa, n, &graphdata);
+    KapaPlotVector (myKapa, n, x->data.F32, "x");
+    KapaPlotVector (myKapa, n, y->data.F32, "y");
+    KapaPlotVector (myKapa, n, dy->data.F32, "dym");
+    KapaPlotVector (myKapa, n, dy->data.F32, "dyp");
+
+    graphdata.color = KapaColorByName ("blue");
+    graphdata.ptype = 0;
+    graphdata.size = 0.5;
+    graphdata.style = 0;
+    graphdata.etype = 0;
+    lineX[0] = graphdata.xmin;
+    lineX[1] = graphdata.xmax;
+    lineY[0] = lineY[1] = mean;
+    KapaPrepPlot (myKapa, 2, &graphdata);
+    KapaPlotVector (myKapa, 2, lineX, "x");
+    KapaPlotVector (myKapa, 2, lineY, "y");
+
+    graphdata.color = KapaColorByName ("red");
+    graphdata.ptype = 0;
+    graphdata.size = 0.5;
+    graphdata.style = 0;
+    graphdata.etype = 0;
+    lineX[0] = graphdata.xmin;
+    lineX[1] = graphdata.xmax;
+    lineY[0] = lineY[1] = mean + error;
+    KapaPrepPlot (myKapa, 2, &graphdata);
+    KapaPlotVector (myKapa, 2, lineX, "x");
+    KapaPlotVector (myKapa, 2, lineY, "y");
+
+    graphdata.color = KapaColorByName ("red");
+    graphdata.ptype = 0;
+    graphdata.size = 0.5;
+    graphdata.style = 0;
+    graphdata.etype = 0;
+    lineX[0] = graphdata.xmin;
+    lineX[1] = graphdata.xmax;
+    lineY[0] = lineY[1] = mean - error;
+    KapaPrepPlot (myKapa, 2, &graphdata);
+    KapaPlotVector (myKapa, 2, lineX, "x");
+    KapaPlotVector (myKapa, 2, lineY, "y");
 
     psFree (x);
     psFree (y);
-
-    // pause and wait for user input:
-    // continue, save (provide name), ??
-    char key[10];
-    fprintf (stdout, "[c]ontinue? ");
-    if (!fgets(key, 8, stdin)) {
-        psWarning("Unable to read option");
-    }
+    psFree (dy);
+
+    pmVisualAskUser(NULL);
+    return true;
+}
+
+bool psphotVisualShowPetrosians (psArray *sources) {
+
+    int Noverlay, NOVERLAY;
+    KiiOverlay *overlay;
+
+    if (!pmVisualIsVisual()) return true;
+
+    int kapa = psphotKapaChannel (1);
+    if (kapa == -1) return false;
+
+    Noverlay = 0;
+    NOVERLAY = 100;
+    ALLOCATE (overlay, KiiOverlay, NOVERLAY);
+
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = sources->data[i];
+
+        if (!source) continue;
+        if (!source->extpars) continue;
+        if (!source->extpars->profile) continue;
+        if (!source->extpars->petrosian_80) continue;
+
+        pmSourceRadialProfile *profile = source->extpars->profile;
+        pmSourceExtendedFlux *petrosian = source->extpars->petrosian_80;
+
+        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
+        overlay[Noverlay].x = source->peak->xf;
+        overlay[Noverlay].y = source->peak->yf;
+        overlay[Noverlay].dx = 2.0*petrosian->radius;
+        overlay[Noverlay].dy = 2.0*petrosian->radius*profile->axes.minor/profile->axes.major;
+        overlay[Noverlay].angle = profile->axes.theta * PS_DEG_RAD;
+        overlay[Noverlay].text = NULL;
+        Noverlay ++;
+        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
+
+        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
+        overlay[Noverlay].x = source->peak->xf;
+        overlay[Noverlay].y = source->peak->yf;
+        overlay[Noverlay].dx = 4.0*petrosian->radius;
+        overlay[Noverlay].dy = 4.0*petrosian->radius*profile->axes.minor/profile->axes.major;
+        overlay[Noverlay].angle = profile->axes.theta * PS_DEG_RAD;
+        overlay[Noverlay].text = NULL;
+        Noverlay ++;
+        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
+    }
+
+    KiiLoadOverlay (kapa, overlay, Noverlay, "red");
+    FREE (overlay);
+
+    pmVisualAskUser(NULL);
     return true;
 }
@@ -1850 +2126 @@
 
 # endif
+
+# if (0)
+// *** make a histogram of the source counts in the x and y directions
+psHistogram *nX = psHistogramAlloc (graphdata.xmin, graphdata.xmax, 50.0);
+psHistogram *nY = psHistogramAlloc (graphdata.ymin, graphdata.ymax, 50.0);
+psVectorHistogram (nX, xFaint, NULL, NULL, 0);
+psVectorHistogram (nY, yFaint, NULL, NULL, 0);
+psVector *dX = psVectorAlloc (nX->nums->n, PS_TYPE_F32);
+psVector *vX = psVectorAlloc (nX->nums->n, PS_TYPE_F32);
+psVector *dY = psVectorAlloc (nY->nums->n, PS_TYPE_F32);
+psVector *vY = psVectorAlloc (nY->nums->n, PS_TYPE_F32);
+for (int i = 0; i < nX->nums->n; i++) {
+    dX->data.F32[i] = nX->nums->data.S32[i];
+    vX->data.F32[i] = 0.5*(nX->bounds->data.F32[i] + nX->bounds->data.F32[i+1]);
+}
+for (int i = 0; i < nY->nums->n; i++) {
+    dY->data.F32[i] = nY->nums->data.S32[i];
+    vY->data.F32[i] = 0.5*(nY->bounds->data.F32[i] + nY->bounds->data.F32[i+1]);
+}
+
+graphdata.color = KapaColorByName ("black");
+graphdata.ptype = 0;
+graphdata.size = 0.0;
+graphdata.style = 0;
+KapaPrepPlot (myKapa, dX->n, &graphdata);
+KapaPlotVector (myKapa, dX->n, dX->data.F32, "x");
+KapaPlotVector (myKapa, vX->n, vX->data.F32, "y");
+
+psFree (nX);
+psFree (dX);
+psFree (vX);
+
+psFree (nY);
+psFree (dY);
+psFree (vY);
+
+# endif
+

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