Index: /trunk/psphot/doc/notes.txt
===================================================================
--- /trunk/psphot/doc/notes.txt	(revision 6426)
+++ /trunk/psphot/doc/notes.txt	(revision 6427)
@@ -1,2 +1,17 @@
+
+2006.02.12
+
+  Last week, I made a new tag (psphot_dev_08) after finishing the work
+  to clean memory leaks and to modify psphot to work with Paul's
+  ppImage infrastructure.  
+
+  I have added a new background model function which generates a
+  median map based on Robert's suggestion of a) half-step windows and
+  b) linear interpolation.
+
+  Now, I am working to remove the constant sky level concepts and make
+  every relevant test use the signal/noise ratio as the trigger.  This
+  has impacts on the peak detections, the threshold for blends, and
+  the threshold for choosing pixels in the fitting routines.
 
 2005.12.23
Index: /trunk/psphot/src/Makefile.am
===================================================================
--- /trunk/psphot/src/Makefile.am	(revision 6426)
+++ /trunk/psphot/src/Makefile.am	(revision 6427)
@@ -21,8 +21,5 @@
 	psphotChoosePSF.c	\
 	psphotEnsemblePSF.c	\
-	psphotFullFit.c		\
 	psphotBlendFit.c	\
-	psphotApplyPSF.c	\
-	psphotFitGalaxies.c	\
 	psphotReplaceUnfit.c	\
 	psphotApResid.c		\
@@ -65,2 +62,7 @@
 tags:
 	etags `find . -name \*.[ch] -print`
+
+
+#	psphotFullFit.c		
+#	psphotApplyPSF.c	
+#	psphotFitGalaxies.c
Index: /trunk/psphot/src/psphot.h
===================================================================
--- /trunk/psphot/src/psphot.h	(revision 6426)
+++ /trunk/psphot/src/psphot.h	(revision 6427)
@@ -33,20 +33,22 @@
 psStats        *psphotImageStats (pmReadout *readout, psMetadata *config);
 psPolynomial2D *psphotImageBackground (pmReadout *readout, psMetadata *config, psStats *sky);
-psArray        *psphotFindPeaks (pmReadout *readout, psMetadata *config, psStats *sky);
+psArray        *psphotFindPeaks (pmReadout *readout, psMetadata *config);
 psArray        *psphotSourceStats (pmReadout *readout, psMetadata *config, psArray *allpeaks);
 bool            psphotRoughClass (psArray *sources, psMetadata *config);
-bool            psphotBasicDeblend (psArray *sources, psMetadata *config, psStats *sky);
-pmPSF          *psphotChoosePSF (psMetadata *config, psArray *sources, psStats *sky);
+bool            psphotBasicDeblend (psArray *sources, psMetadata *config);
+pmPSF          *psphotChoosePSF (psMetadata *config, psArray *sources);
 void            psphotOutput (pmReadout *readout, psMetadata *arguments);
-psImage        *psphotImageMedian (pmReadout *readout, psMetadata *config, psStats *sky);
+psImage        *psphotImageMedian (pmReadout *readout, psMetadata *config);
 
 // optional object analysis steps
-bool 	        psphotEnsemblePSF (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
-bool            psphotFullFit (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
-bool            psphotBlendFit (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
+bool 	        psphotEnsemblePSF (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf);
+bool            psphotBlendFit (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf);
 bool            psphotReplaceUnfit (psArray *sources);
-bool            psphotApplyPSF (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
-bool            psphotFitExtended (pmReadout *readout, psMetadata *config, psArray *sources, psStats *skyStats);
 bool            psphotApResid (pmReadout *readout, psArray *sources, psMetadata *config, pmPSF *psf);
+
+// XXX deprecate these?
+// bool            psphotFullFit (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
+// bool            psphotApplyPSF (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
+// bool            psphotFitExtended (pmReadout *readout, psMetadata *config, psArray *sources, psStats *skyStats);
 
 // basic support functions
@@ -63,7 +65,7 @@
 
 // functions to set the correct source pixels
-bool 		psphotInitRadiusPSF (psMetadata *config, psStats *sky, pmModelType type);
+bool 		psphotInitRadiusPSF (psMetadata *config, pmModelType type);
 bool 		psphotCheckRadiusPSF (pmReadout *readout, pmSource *source, pmModel *model);
-bool 		psphotInitRadiusEXT (psMetadata *config, psStats *sky, pmModelType type);
+bool 		psphotInitRadiusEXT (psMetadata *config, pmModelType type);
 bool 		psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model);
 bool 	        psphotDefinePixels (pmSource *mySource, const pmReadout *readout, psF32 x, psF32 y, psF32 Radius);
@@ -111,5 +113,5 @@
 //  functions to support the source fitting process
 bool 		psphotInitLimitsPSF (psMetadata *config);
-bool		psphotInitLimitsEXT (psMetadata *config, psStats *sky);
+bool		psphotInitLimitsEXT (psMetadata *config);
 bool		psphotFitBlend (pmReadout *readout, pmSource *source);
 bool		psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources);
Index: /trunk/psphot/src/psphotBasicDeblend.c
===================================================================
--- /trunk/psphot/src/psphotBasicDeblend.c	(revision 6426)
+++ /trunk/psphot/src/psphotBasicDeblend.c	(revision 6427)
@@ -2,5 +2,5 @@
 
 // 2006.02.07 : no leaks
-bool psphotBasicDeblend (psArray *sources, psMetadata *config, psStats *sky) { 
+bool psphotBasicDeblend (psArray *sources, psMetadata *config) { 
 
     int N;
@@ -18,6 +18,4 @@
     float NSIGMA   = psMetadataLookupF32 (&status, config, "DEBLEND_SKY_NSIGMA");
     if (!status) NSIGMA = 5.0;
-
-    float minThreshold = NSIGMA*sky->sampleStdev;
 
     // we need sources spatially-sorted to find overlaps
@@ -83,6 +81,10 @@
 	// generate source contour (1/4 peak counts)
 	// set the threshold based on user inputs
-	threshold = FRACTION * (source->peak->counts - source->moments->Sky) + source->moments->Sky; 
-	threshold = PS_MAX (threshold, minThreshold);
+	
+	// threshold is fraction of the source peak flux
+	// image is background subtracted; source->moments->Sky should always be 0.0
+	threshold = FRACTION * source->moments->Peak;
+	// threshold is no less than NSIGMA above the local median sigma?
+	threshold = PS_MAX (threshold, NSIGMA*source->moments->dSky);
 
 	// generate a basic contour (set of x,y coordinates at-or-below flux level)
Index: /trunk/psphot/src/psphotBlendFit.c
===================================================================
--- /trunk/psphot/src/psphotBlendFit.c	(revision 6426)
+++ /trunk/psphot/src/psphotBlendFit.c	(revision 6427)
@@ -2,5 +2,5 @@
 
 // XXX I don't like this name
-bool psphotBlendFit (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) { 
+bool psphotBlendFit (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf) { 
 
     psTimerStart ("psphot");
@@ -10,6 +10,6 @@
     
     psphotInitLimitsPSF (config);
-    psphotInitLimitsEXT (config, sky);
-    psphotInitRadiusPSF (config, sky, psf->type);
+    psphotInitLimitsEXT (config);
+    psphotInitRadiusPSF (config, psf->type);
 
     for (int i = 0; i < sources->n; i++) {
Index: /trunk/psphot/src/psphotChoosePSF.c
===================================================================
--- /trunk/psphot/src/psphotChoosePSF.c	(revision 6426)
+++ /trunk/psphot/src/psphotChoosePSF.c	(revision 6427)
@@ -3,5 +3,5 @@
 // 2006.02.07 : no leaks!
 // try PSF models and select best option
-pmPSF *psphotChoosePSF (psMetadata *config, psArray *sources, psStats *skystats) { 
+pmPSF *psphotChoosePSF (psMetadata *config, psArray *sources) { 
 
     bool       	    status;
@@ -52,10 +52,7 @@
     psListIterator *iter = psListIteratorAlloc (list, PS_LIST_HEAD, FALSE); 
     for (int i = 0; i < models->n; i++) { 
-	// XXX psListGetAndIncrement does not increment ref counter!?
-	// XXX if it did, i should free 'item' below
 	psMetadataItem *item = psListGetAndIncrement (iter);
 	modelName = item->data.V;
 	models->data[i] = pmPSFtryModel (stars, modelName, RADIUS);
-	// psFree (item);
     }
     psFree (iter);
Index: /trunk/psphot/src/psphotEnsemblePSF.c
===================================================================
--- /trunk/psphot/src/psphotEnsemblePSF.c	(revision 6426)
+++ /trunk/psphot/src/psphotEnsemblePSF.c	(revision 6427)
@@ -3,5 +3,5 @@
 // 2006.02.07 : no leaks!
 // fit all reasonable sources with the linear PSF model
-bool psphotEnsemblePSF (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) { 
+bool psphotEnsemblePSF (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf) { 
 
     bool  status;
@@ -15,11 +15,10 @@
     sources = psArraySort (sources, psphotSortByY);
 
+    // radius for source moments if star is saturated
+    // XXX we already do this in psphotSourceStats, is it still needed?
     float INNER_RADIUS     = psMetadataLookupF32 (&status, config, "SKY_INNER_RADIUS");
-    float PSF_FIT_NSIGMA   = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA");
-    float PSF_FIT_PADDING  = psMetadataLookupF32 (&status, config, "PSF_FIT_PADDING");
 
-    // set the object surface-brightness limit for fitted pixels
-    float FLUX_LIMIT  = PSF_FIT_NSIGMA * sky->sampleStdev;
-    psLogMsg ("psphot.ensemble", 4, "fitting pixels with at least %f object counts\n", FLUX_LIMIT);
+    // setup the PSF fit details
+    psphotInitRadiusPSF (config, psf->type);
 
     // pre-calculate all model pixels
@@ -29,4 +28,5 @@
     index->n = 0;
 
+    // option to limit analysis to a specific region
     bool UseAnalysisRegion = false;
     psRegion AnalysisRegion = {0, 0, 0, 0};
@@ -98,13 +98,6 @@
 	y = model->params->data.F32[3];
 
-	// get function which specifies the radius at which the model hits the given flux
-	pmModelRadius modelRadius = pmModelRadius_GetFunction (psf->type);
-
 	// set the fit radius based on the object flux limit and the model
-	model->radius = modelRadius (model->params, FLUX_LIMIT) + PSF_FIT_PADDING;
-	if (isnan(model->radius)) psAbort ("psphotEnsemblePSF", "error in radius");
-
-	// if needed, ask for more object pixels
-	psphotRedefinePixels (inSource, readout, x, y, model->radius);
+	psphotCheckRadiusPSF (readout, inSource, model);
 
 	// make temporary copies of the image pixels and mask
@@ -116,5 +109,5 @@
 	psImage *mask = otSource->mask;
 
-	// set model to unit peak, zero sky (we assume sky is constant)
+	// set model to unit peak, zero sky (we assume sky is subtracted)
 	model->params->data.F32[0] = 0.0;
 	model->params->data.F32[1] = 1.0;
@@ -165,5 +158,5 @@
     psLogMsg ("psphot.emsemble", 4, "built matrix: %f (%d elements)\n", psTimerMark ("psphot"), sparse->Nelem);
 
-    // solve for normalization terms (need include local sky)
+    // solve for normalization terms (need include local sky?)
     psSparseResort (sparse);
     psVector *norm = psSparseSolve (NULL, sparse, 3);
Index: /trunk/psphot/src/psphotFindPeaks.c
===================================================================
--- /trunk/psphot/src/psphotFindPeaks.c	(revision 6426)
+++ /trunk/psphot/src/psphotFindPeaks.c	(revision 6427)
@@ -2,5 +2,5 @@
 
 // 2006.02.02 : no leaks
-psArray *psphotFindPeaks (pmReadout *readout, psMetadata *config, psStats *sky) {
+psArray *psphotFindPeaks (pmReadout *readout, psMetadata *config) {
 
     bool  status = false;
@@ -9,5 +9,7 @@
     float threshold;
 
-    // smooth the image 
+    // smooth the image and weight map
+
+    psphotSaveImage (NULL, readout->image, "image.fits");
 
     psTimerStart ("psphot");
@@ -16,25 +18,38 @@
     NSIGMA = psMetadataLookupF32 (&status, config, "PEAKS_SMOOTH_NSIGMA");
 
-    psImage *smooth = psImageCopy (NULL, readout->image, PS_TYPE_F32);
-    psphotSaveImage (NULL, smooth, "presmooth.fits");
+    psImage *smooth_im = psImageCopy (NULL, readout->image, PS_TYPE_F32);
+    psImageSmooth (smooth_im, SIGMA, NSIGMA);
+    psLogMsg ("psphot", 4, "smooth image: %f sec\n", psTimerMark ("psphot"));
 
-    psImageSmooth (smooth, SIGMA, NSIGMA);
-    psLogMsg ("psphot", 4, "smooth: %f sec\n", psTimerMark ("psphot"));
+    psphotSaveImage (NULL, smooth_im, "smooth_im.fits");
 
-    psphotSaveImage (NULL, smooth, "smooth.fits");
+    psImage *smooth_wt = psImageCopy (NULL, readout->weight, PS_TYPE_F32);
+    psImageSmooth (smooth_wt, SIGMA, NSIGMA);
+    psLogMsg ("psphot", 4, "smooth weight: %f sec\n", psTimerMark ("psphot"));
+
+    psphotSaveImage (NULL, smooth_wt, "smooth_wt.fits");
+
+    psTimerStart ("psphot");
+    for (int j = 0; j < smooth_im->numRows; j++) {
+	for (int i = 0; i < smooth_im->numCols; i++) {
+	    if (smooth_wt->data.F32[j][i] <= 0) continue;
+	    smooth_im->data.F32[j][i] = PS_SQR(smooth_im->data.F32[j][i]) / smooth_wt->data.F32[j][i];
+	}
+    }
+    psLogMsg ("psphot", 4, "built S/N map: %f sec\n", psTimerMark ("psphot"));
 
     psTimerStart ("psphot");
 
     // set peak threshold
-    NSIGMA    = psMetadataLookupF32 (&status, config, "PEAKS_NSIGMA_LIMIT");
-    
-    threshold = NSIGMA*sky->sampleStdev + sky->sampleMean; 
-    // threshold = NSIGMA*sky->sampleStdev;
+    NSIGMA = psMetadataLookupF32 (&status, config, "PEAKS_NSIGMA_LIMIT");
+    // threshold = NSIGMA*sky->sampleStdev + sky->sampleMean; 
+    threshold = PS_SQR(NSIGMA);
     psLogMsg ("psphot", 4, "threshold: %f DN\n", threshold);
 
     // find the peaks in the smoothed image 
-    psArray *peaks = pmFindImagePeaks (smooth, threshold);
+    psArray *peaks = pmFindImagePeaks (smooth_im, threshold);
     if (peaks == NULL) psAbort ("find peaks", "no peaks found");
-    psFree (smooth);
+    psFree (smooth_im);
+    psFree (smooth_wt);
 
     // optional dump of all peak data
Index: /trunk/psphot/src/psphotImageMedian.c
===================================================================
--- /trunk/psphot/src/psphotImageMedian.c	(revision 6426)
+++ /trunk/psphot/src/psphotImageMedian.c	(revision 6427)
@@ -11,5 +11,5 @@
 // generate the median in NxN boxes, clipping heavily
 // linear interpolation to generate full-scale model
-psImage *psphotImageMedian (pmReadout *readout, psMetadata *config, psStats *sky) 
+psImage *psphotImageMedian (pmReadout *readout, psMetadata *config) 
 { 
     bool status;
@@ -20,5 +20,4 @@
     if (!status) MAX_SAMPLE_PIXELS = 1000;
 
-    // fit background to subset of image points within 3 sigma of sky
     psTimerStart ("psphot");
 
@@ -262,4 +261,6 @@
     value = value / npts;
     
+    // XXX correct for selection bias??
+
     psFree (values);
     return value;
Index: /trunk/psphot/src/psphotRadiusChecks.c
===================================================================
--- /trunk/psphot/src/psphotRadiusChecks.c	(revision 6426)
+++ /trunk/psphot/src/psphotRadiusChecks.c	(revision 6427)
@@ -1,17 +1,14 @@
 # include "psphot.h"
 
-static float PSF_OUTER_RADIUS;
+static float PSF_FIT_NSIGMA;
 static float PSF_FIT_PADDING;
-static float PSF_FLUX_LIMIT;
 static pmModelRadius modelRadiusPSF;
 
-bool psphotInitRadiusPSF (psMetadata *config, psStats *sky, pmModelType type) {
+bool psphotInitRadiusPSF (psMetadata *config, pmModelType type) {
 
     bool status;
 
-    float PSF_FIT_NSIGMA = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA");
-    PSF_OUTER_RADIUS     = psMetadataLookupF32 (&status, config, "PSF_OUTER_RADIUS");
-    PSF_FIT_PADDING      = psMetadataLookupF32 (&status, config, "PSF_FIT_PADDING");
-    PSF_FLUX_LIMIT       = PSF_FIT_NSIGMA * sky->sampleStdev;
+    PSF_FIT_NSIGMA   = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA");
+    PSF_FIT_PADDING  = psMetadataLookupF32 (&status, config, "PSF_FIT_PADDING");
 
     // this function specifies the radius at this the model hits the given flux
@@ -22,7 +19,9 @@
 bool psphotCheckRadiusPSF (pmReadout *readout, pmSource *source, pmModel *model) {
 
+    pmMoments *moments = source->moments;
+    if (moments == NULL) return false;
+
     // set the fit radius based on the object flux limit and the model
-    // XXX EAM : FLUX_LIMIT should be set based on local sky model (not global median)
-    model->radius = modelRadiusPSF (model->params, PSF_FLUX_LIMIT) + PSF_FIT_PADDING;
+    model->radius = modelRadiusPSF (model->params, PSF_FIT_NSIGMA*moments->dSky) + PSF_FIT_PADDING;
     if (isnan(model->radius)) psAbort ("apply_psf_model", "error in radius");
 	
@@ -35,17 +34,14 @@
 }
 
-static float EXT_OUTER_RADIUS;
+static float EXT_FIT_NSIGMA;
 static float EXT_FIT_PADDING;
-static float EXT_FLUX_LIMIT;
 static pmModelRadius modelRadiusEXT;
 
-bool psphotInitRadiusEXT (psMetadata *config, psStats *sky, pmModelType type) {
+bool psphotInitRadiusEXT (psMetadata *config, pmModelType type) {
 
     bool status;
 
-    float EXT_FIT_NSIGMA = psMetadataLookupF32 (&status, config, "EXT_FIT_NSIGMA");
-    EXT_OUTER_RADIUS	 = psMetadataLookupF32 (&status, config, "EXT_OUTER_RADIUS");
-    EXT_FIT_PADDING      = psMetadataLookupF32 (&status, config, "EXT_FIT_PADDING");
-    EXT_FLUX_LIMIT       = EXT_FIT_NSIGMA * sky->sampleStdev;
+    EXT_FIT_NSIGMA   = psMetadataLookupF32 (&status, config, "EXT_FIT_NSIGMA");
+    EXT_FIT_PADDING  = psMetadataLookupF32 (&status, config, "EXT_FIT_PADDING");
 
     // this function specifies the radius at this the model hits the given flux
@@ -56,6 +52,9 @@
 bool psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model) {
 
+    pmMoments *moments = source->moments;
+    if (moments == NULL) return false;
+
     // set the fit radius based on the object flux limit and the model
-    model->radius = modelRadiusEXT (model->params, EXT_FLUX_LIMIT) + EXT_FIT_PADDING;
+    model->radius = modelRadiusEXT (model->params, EXT_FIT_NSIGMA*moments->dSky) + EXT_FIT_PADDING;
     if (isnan(model->radius)) psAbort (__func__, "error in radius");
 
Index: /trunk/psphot/src/psphotReadout.c
===================================================================
--- /trunk/psphot/src/psphotReadout.c	(revision 6426)
+++ /trunk/psphot/src/psphotReadout.c	(revision 6427)
@@ -4,22 +4,15 @@
 bool psphotReadout (pmReadout *readout, psMetadata *config) {
 
-    psArray     *sources = NULL;
-    psArray     *peaks   = NULL;
-    pmPSF       *psf     = NULL;
-    psStats     *sky     = NULL;
+    psArray     *sources  = NULL;
+    psArray     *peaks    = NULL;
+    psImage     *skymodel = NULL;
+    pmPSF       *psf      = NULL;
     bool         status;
-    psPolynomial2D *model = NULL;
 
-    psphotImageMedian (readout, config, sky);
-
-    // measure image stats for initial guess 
-    sky = psphotImageStats (readout, config);
-
-    // generate a background model (currently, 2D polynomial)
-    // XXX this should be available to be re-added to the original image
-    model = psphotImageBackground (readout, config, sky);
+    // generate a background model (median, smoothed image)
+    skymodel = psphotImageMedian (readout, config);
 
     // find the peaks in the image
-    peaks = psphotFindPeaks (readout, config, sky);
+    peaks = psphotFindPeaks (readout, config);
 
     // construct sources and measure basic stats
@@ -30,32 +23,20 @@
 
     // mark blended peaks PS_SOURCE_BLEND
-    psphotBasicDeblend (sources, config, sky);
+    psphotBasicDeblend (sources, config);
 
     // use bright stellar objects to measure PSF
-    psf = psphotChoosePSF (config, sources, sky);
+    psf = psphotChoosePSF (config, sources);
 
-    // XXX change FITMODE to a string
+    // select analysis method
     char *FITMODE = psMetadataLookupStr (&status, config, "FITMODE");
     if (!strcasecmp(FITMODE, "ENSEMBLE")) {
-	psphotEnsemblePSF (readout, config, sources, psf, sky);
-    }
-
-    if (!strcasecmp(FITMODE, "FULL")) {
-	psphotEnsemblePSF (readout, config, sources, psf, sky);
-	psphotFullFit (readout, config, sources, psf, sky);
-	psphotReplaceUnfit (sources);
-	psphotApResid (readout, sources, config, psf);
+	psphotEnsemblePSF (readout, config, sources, psf);
     }
 
     if (!strcasecmp(FITMODE, "BLEND")) {
-	psphotEnsemblePSF (readout, config, sources, psf, sky);
-	psphotBlendFit (readout, config, sources, psf, sky);
+	psphotEnsemblePSF (readout, config, sources, psf);
+	psphotBlendFit (readout, config, sources, psf);
 	psphotReplaceUnfit (sources);
 	psphotApResid (readout, sources, config, psf);
-    }
-
-    if (!strcasecmp(FITMODE, "BASIC")) {
-	psphotApplyPSF (readout, config, sources, psf, sky);
-	psphotFitExtended (readout, config, sources, sky);
     }
 
@@ -65,17 +46,35 @@
     psphotUpdateHeader (header, config);
 
+    // psphotSkyReplace (readout, skymodel);
+
     // need to do something with the sources, psf, and sky
     status = psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES",   PS_DATA_ARRAY,   "psphot sources", sources);
     status = psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.PSF",       PS_DATA_UNKNOWN, "psphot psf", psf);
-    status = psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SKY.MEAN",  PS_DATA_F32,     "psphot sky mean", sky->sampleMean);
-    status = psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SKY.SIGMA", PS_DATA_F32,     "psphot sky stdev", sky->sampleStdev);
+    // XXX : replace? status = psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SKY.MEAN",  PS_DATA_F32,     "psphot sky mean", sky->sampleMean);
+    // XXX : replace? status = psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SKY.SIGMA", PS_DATA_F32,     "psphot sky stdev", sky->sampleStdev);
 
     // free up the local copies of the data
     psFree (psf);
-    psFree (sky);
-    psFree (model);
     psFree (peaks);
     psFree (sources);
+    psFree (skymodel);
 
     return true;
 }
+
+
+// XXX Deprecate or allow these versions?
+# if (0)
+    if (!strcasecmp(FITMODE, "FULL")) {
+	psphotEnsemblePSF (readout, config, sources, psf);
+	psphotFullFit (readout, config, sources, psf);
+	psphotReplaceUnfit (sources);
+	psphotApResid (readout, sources, config, psf);
+    }
+
+    if (!strcasecmp(FITMODE, "BASIC")) {
+	psphotApplyPSF (readout, config, sources, psf);
+	psphotFitExtended (readout, config, sources);
+    }
+
+# endif
Index: /trunk/psphot/src/psphotSourceFits.c
===================================================================
--- /trunk/psphot/src/psphotSourceFits.c	(revision 6426)
+++ /trunk/psphot/src/psphotSourceFits.c	(revision 6427)
@@ -45,5 +45,5 @@
 static pmModelType modelTypeEXT;
 
-bool psphotInitLimitsEXT (psMetadata *config, psStats *sky) {
+bool psphotInitLimitsEXT (psMetadata *config) {
 
     bool status;
@@ -56,5 +56,5 @@
     char *modelNameEXT = psMetadataLookupStr (&status, config, "EXT_MODEL");
     modelTypeEXT = pmModelSetType (modelNameEXT);
-    psphotInitRadiusEXT (config, sky, modelTypeEXT);
+    psphotInitRadiusEXT (config, modelTypeEXT);
 
     return true;
@@ -250,5 +250,6 @@
 	}
 
-	model->params->data.F32[1] = blend->peak->counts;
+	// XXX assume local sky is 0.0?
+	model->params->data.F32[1] = blend->moments->Peak - blend->moments->Sky;
 	model->params->data.F32[2] = blend->peak->x;
 	model->params->data.F32[3] = blend->peak->y;
Index: /trunk/psphot/src/psphotSourceStats.c
===================================================================
--- /trunk/psphot/src/psphotSourceStats.c	(revision 6426)
+++ /trunk/psphot/src/psphotSourceStats.c	(revision 6427)
@@ -25,5 +25,4 @@
 
 	// allocate image, weight, mask arrays for each peak (square of radius OUTER)
-	// XXX EAM : this has 12 leaks (v.5)
 	psphotDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER);
 
@@ -31,4 +30,12 @@
 	// XXX EAM : this should use ROBUST not SAMPLE median, but it is broken
 	status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER);
+	if (!status) {
+	  psFree (source);
+	  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);
 	if (!status) {
 	  psFree (source);