Index: /trunk/psphot/src/LocalSky.c
===================================================================
--- /trunk/psphot/src/LocalSky.c	(revision 4215)
+++ /trunk/psphot/src/LocalSky.c	(revision 4216)
@@ -1,35 +1,18 @@
 # include "psphot.h"
 
-psSource *pmSourceDefinePixels(const psImageData *imdata,
-			       const psPeak *peak,
-			       psF32 Radius)
+bool pmSourceDefinePixels(psSource *mySource, 
+			  const psImageData *imdata,
+			  psF32 x, 
+			  psF32 y,
+			  psF32 Radius)
 {
-    // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1)
-    psS32 SubImageStartRow  = PS_MAX (0, peak->y - Radius);
-    psS32 SubImageEndRow    = PS_MIN (image->numRows, peak->y + Radius + 1);
-    psS32 SubImageStartCol  = PS_MAX (0, peak->x - Radius);
-    psS32 SubImageEndCol    = PS_MIN (image->numCols, peak->x + Radius + 1);
+    // Grab a subimage of the original image of size (2 * outerRadius).
+    srcRegion = psRegionSquare (x, y, Radius);
+    srcRegion = psRegionForImage (srcRegion, imdata->image);
 
-    // Grab a subimage of the original image of size (2 * outerRadius).
-    psImage *subImage = psImageSubset(imdata->image,
-                                      SubImageStartCol,
-                                      SubImageStartRow,
-                                      SubImageEndCol,
-                                      SubImageEndRow);
+    psImage *subImage = psImageSubset(imdata->image, srcRegion);
+    psImage *subImageMask = psImageSubset(imdata->mask, srcRegion);
+    psImage *subImageNoise = psImageSubset(imdata->noise, srcRegion);
 
-    psImage *subImageMask = psImageSubset(imdata->mask,
-					  SubImageStartCol,
-					  SubImageStartRow,
-					  SubImageEndCol,
-					  SubImageEndRow);
-
-    psImage *subImageNoise = psImageSubset(imdata->noise,
-					   SubImageStartCol,
-					   SubImageStartRow,
-					   SubImageEndCol,
-					   SubImageEndRow);
-
-    psSource *mySource = pmSourceAlloc();
-    mySource->peak = peak;
     mySource->pixels = subImage;
     mySource->mask = subImageMask;
@@ -49,25 +32,11 @@
     psImage *peak  = source->peak;
 
-    // Loop through the subimage, mask off pixels in the inner square.
-    // this uses a static mask value of 0x08
-    psS32 StartCol  = PS_MAX (0, peak->x - mask->col0 - Radius);
-    psS32 StartRow  = PS_MAX (0, peak->y - mask->row0 - Radius);
-    psS32 EndCol    = PS_MIN (mask->numCols, peak->x - mask->col0 + Radius + 1);
-    psS32 EndRow    = PS_MIN (mask->numRows, peak->y - mask->row0 + Radius + 1);
-    for (psS32 row = StartRow; row < EndRow; row++) {
-        for (psS32 col = StartCol; col < EndCol; col++) {
-	    mask->data.U8[row][col] |= 0x80;
-        }
-    }
+    srcRegion = psRegionSquare (peak->x, peak->y, Radius);
+    srcRegion = psRegionForImage (srcRegion, imdata->mask);
+    psImageMaskRegion (imdata->mask, region, OR, 0x80);
 
     psStats *myStats = psStatsAlloc(statsOptions);
     myStats = psImageStats(myStats, image, mask, 0xff);
-    
-    // clear the mask on the inner pixels
-    for (psS32 row = StartRow; row < EndRow; row++) {
-        for (psS32 col = StartCol; col < EndCol; col++) {
-	    mask->data.U8[row][col] &= 0x7f;
-        }
-    }
+    psImageMaskRegion (imdata->mask, region, AND, 0x7f);
 
     mySource->moments = pmMomentsAlloc();
@@ -79,34 +48,7 @@
 }
 
-// set mask pixels within the circle
-bool psImageMaskCircle(psImage *mask,
-		       psF32 x, 
-		       psF32 y,
-		       psF32 radius,
-		       int maskValue)
-{
-    // XXX EAM : for the circle to stay on the image
-    // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1)
-    psS32 SubImageStartRow  = PS_MAX (0, x + mask->row0 - radiusS32);
-    psS32 SubImageEndRow    = PS_MIN (mask->numRows, SubImageCenterRow + radiusS32 + 1);
-    psS32 SubImageStartCol  = PS_MAX (0, SubImageCenterCol - radiusS32);
-    psS32 SubImageEndCol    = PS_MIN (mask->numCols, SubImageCenterCol + radiusS32 + 1);
-
-    // Loop through the subimage mask, set maskValue bits
-    for (psS32 row = 0 ; row < source->mask->numRows; row++) {
-        for (psS32 col = 0 ; col < source->mask->numCols; col++) {
-
-            if (CheckRadius2((psF32) radiusS32,
-                             (psF32) radiusS32,
-                             radius,
-                             (psF32) col,
-                             (psF32) row)) {
-                source->mask->data.U8[row][col] = 0;
-            } else {
-                source->mask->data.U8[row][col] = 1;
-            }
-        }
-    }
-    return(true);
+// define a square region centered on the given coordinate
+psRegion *psRegionSquare (psF32 x, psF32 y, psF32 radius) {
+    region = psRegionAlloc (x - Radius, y - Radius, x + Radius + 1, y + Radius + 1);
+    return (region);
 }
-
Index: /trunk/psphot/src/apply_psf_model.c
===================================================================
--- /trunk/psphot/src/apply_psf_model.c	(revision 4215)
+++ /trunk/psphot/src/apply_psf_model.c	(revision 4216)
@@ -4,5 +4,4 @@
 // PSFSTAR objects will be refitted
 // run this function to a specific flux limit?
-
 
 bool apply_psf_model (psImage *image, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) 
@@ -18,12 +17,7 @@
     float FIT_NSIGMA  = psMetadataLookupF32 (&status, config, "FIT_NSIGMA");
     float FIT_PADDING = psMetadataLookupF32 (&status, config, "FIT_PADDING");
-
-    float XBORDER  = psMetadataLookupF32 (&status, config, "XBORDER");
-    float YBORDER  = psMetadataLookupF32 (&status, config, "YBORDER");
-    psRegion *keep = psRegionAlloc (XBORDER, -XBORDER, YBORDER, -YBORDER);
-    keep           = psRegionForImage (keep, image, keep);
-
     float shapeNsigma = psMetadataLookupF32 (&status, config, "PSF_SHAPE_NSIGMA");
-    // float snFaint     = psMetadataLookupF32 (&status, config, "FAINT_SN_LIM");
+    float snFaint     = psMetadataLookupF32 (&status, config, "FAINT_SN_LIM");
+    float OUTER       = psMetadataLookupF32 (&status, config, "OUTER_RADIUS");
 
     // set the object surface-brightness limit for fitted pixels
@@ -46,4 +40,6 @@
 	// set PSF parameters for this model
 	model  = psModelFromPSF (model, psf);
+	x = model->params->data.F32[2];
+	y = model->params->data.F32[3];
 
 	// set the fit radius based on the object flux limit and the model
@@ -54,12 +50,20 @@
 	  continue;
 	}
-	    
-	// set the valid/invalid pixel mask
-	pmSourceSetPixelCircle (source, image, model->radius);
-	pmSourceMaskRegion (source, keep);
+	
+	// mask off saturated pixels (move this to setup?)
 	pmSourceMaskSaturated (source, SATURATE);
 
-	// fit as PSF, not FLT (skip poor fits)
-	if (!pmSourceFitModel (source, model, true)) continue;
+	// check if we need to redefine the pixels
+	if (model->radius > OUTER) {
+	  // allocate image, noise, mask arrays for each peak (square of radius OUTER)
+	  pmSourceDefinePixels (source, imdata, x, y, OUTER);
+	}
+
+	// set the valid/invalid pixel mask
+	psImageKeepCircle (source->mask, x, y, model->radius, OR, 0x80);
+	status = pmSourceFitModel (source, model, true);
+	psImageKeepCircle (source->mask, x, y, model->radius, AND, 0x7f);
+	if (!status) continue;
+
 	source->modelPSF = model;
 	Niter += model[0].nIter;
Index: /trunk/psphot/src/choose_psf_model.c
===================================================================
--- /trunk/psphot/src/choose_psf_model.c	(revision 4215)
+++ /trunk/psphot/src/choose_psf_model.c	(revision 4216)
@@ -23,11 +23,6 @@
 
     // define model fit pixels
-    // XXX use a different radius - PSF_FIT_RADIUS?
+    // check that FIT_RADIUS < OUTER_RADIUS
     float RADIUS   = psMetadataLookupF32 (&status, config, "FIT_RADIUS");
-
-    for (int i = 0; i < stars->n; i++) {
-	psSource *source = stars->data[i];
-	pmSourceSetPixelCircle (source, image, RADIUS);
-    }
 
     // the user specifies 1 or more PSF model options
@@ -38,5 +33,5 @@
 	sprintf (key, "PSF_MODEL_%d", i);
 	modelName = psMetadataLookupPtr (&status, config, key);
-	tests->data[i] = pmPSF_TestModel (stars, modelName);
+	tests->data[i] = pmPSF_TestModel (stars, modelName, RADIUS);
     }
 
Index: /trunk/psphot/src/find_peaks.c
===================================================================
--- /trunk/psphot/src/find_peaks.c	(revision 4215)
+++ /trunk/psphot/src/find_peaks.c	(revision 4216)
@@ -16,10 +16,8 @@
     NSIGMA = psMetadataLookupF32 (&status, config, "SMOOTH_NSIGMA");
 
-    psImage *image  = imdata->image;
-    psImage *smooth = psImageCopy (NULL, image, PS_TYPE_F32);
+    psImage *smooth = psImageCopy (NULL, imdata->image, PS_TYPE_F32);
     psImageSmooth (smooth, SIGMA, NSIGMA);
 
     psLogMsg ("psphot", 4, "smooth: %f sec\n", psTimerMark ("psphot"));
-
 
     // find the peaks in the smoothed image 
@@ -33,4 +31,5 @@
     psArray *peaks = pmFindImagePeaks (smooth, threshold);
     if (peaks == NULL) psAbort ("find peaks", "no peaks found");
+    psFree (smooth);
     // XXX Need to gracefully handle no source detections
 
Index: /trunk/psphot/src/image_stats.c
===================================================================
--- /trunk/psphot/src/image_stats.c	(revision 4215)
+++ /trunk/psphot/src/image_stats.c	(revision 4216)
@@ -42,4 +42,5 @@
     float GAIN  = psMetadataLookupF32 (&status, config, "GAIN");
 
+    // convert instrumental background to poisson stats
     sky = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
     sky->sampleMean   = stats->sampleMedian;
Index: /trunk/psphot/src/psphot-utils.c
===================================================================
--- /trunk/psphot/src/psphot-utils.c	(revision 4215)
+++ /trunk/psphot/src/psphot-utils.c	(revision 4216)
@@ -129,4 +129,13 @@
 }
 
+bool DumpImage (psImage *image, char *filename) {
+
+    unlink (filename);
+    psFits *fits = psFitsAlloc (filename);
+    psFitsWriteImage (fits, NULL, image, 0, NULL);
+    psFree (fits);
+    return;
+}
+
 void pmSourceMaskSaturated (psSource *source, float saturate) {
 
@@ -142,8 +151,26 @@
 }
 
-void pmSourceMaskRegion (psSource *source, psRegion *region) {
-
-    psImage *image = source->pixels;
-    psImage *mask  = source->mask;
+// mask the area contained by the region
+// the region is defined wrt the parent image
+void psImageMaskRegion (psImage *image, psRegion *region, bool logical_and, int maskValue) {
+
+    for (int iy = 0; iy < image->numRows; iy++) {
+	for (int ix = 0; ix < image->numCols; ix++) {
+	    if (ix + image->col0 <  region->x0) continue;
+	    if (ix + image->col0 >= region->x1) continue;
+	    if (iy + image->row0 <  region->y0) continue;
+	    if (iy + image->row0 >= region->y1) continue;
+	    if (logical_and) {
+		image->data.U8[iy][ix] &= maskValue;
+	    } else {
+		image->data.U8[iy][ix] |= maskValue;
+	    }
+	}
+    }
+}
+
+// mask the area not contained by the region
+// the region is defined wrt the parent image
+void psImageKeepRegion (psImage *image, psRegion *region, bool logical_and, int maskValue) {
 
     for (int iy = 0; iy < image->numRows; iy++) {
@@ -155,5 +182,9 @@
 	    continue;
 	maskit:
-	    mask->data.U8[iy][ix] |= 0x02;
+	    if (logical_and) {
+		image->data.U8[iy][ix] &= maskValue;
+	    } else {
+		image->data.U8[iy][ix] |= maskValue;
+	    }
 	}
     }
@@ -161,29 +192,41 @@
 
 // mask the area contained by the region
-void psImageMaskRegion (psImage *image, psRegion *region, int maskValue) {
-
-    for (int iy = 0; iy < image->numRows; iy++) {
-	for (int ix = 0; ix < image->numCols; ix++) {
-	    if (ix + image->col0 <  region->x0) continue;
-	    if (ix + image->col0 >= region->x1) continue;
-	    if (iy + image->row0 <  region->y0) continue;
-	    if (iy + image->row0 >= region->y1) continue;
-	    image->data.U8[iy][ix] |= maskValue;
-	}
-    }
-}
-
-// mask the area not contained by the region
-void psImageKeepRegion (psImage *image, psRegion *region, int maskValue) {
-
-    for (int iy = 0; iy < image->numRows; iy++) {
-	for (int ix = 0; ix < image->numCols; ix++) {
-	    if (ix + image->col0 <  region->x0) goto maskit;
-	    if (ix + image->col0 >= region->x1) goto maskit;
-	    if (iy + image->row0 <  region->y0) goto maskit;
-	    if (iy + image->row0 >= region->y1) goto maskit;
-	    continue;
-	maskit:
-	    image->data.U8[iy][ix] |= maskValue;
+// the region is defined wrt the parent image
+void psImageMaskCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue) {
+
+    double R2 = PS_SQR(radius);
+
+    for (int iy = 0; iy < image->numRows; iy++) {
+	for (int ix = 0; ix < image->numCols; ix++) {
+	    dx = ix + image->col0 - x;
+	    dy = iy + image->row0 - y;
+	    r2 = PS_SQR(dx) + PS_SQR(dy);
+	    if (r2 > R2) continue;
+	    if (logical_and) {
+		image->data.U8[iy][ix] &= maskValue;
+	    } else {
+		image->data.U8[iy][ix] |= maskValue;
+	    }
+	}
+    }
+}
+
+// mask the area contained by the region
+// the region is defined wrt the parent image
+void psImageKeepCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue) {
+
+    double R2 = PS_SQR(radius);
+
+    for (int iy = 0; iy < image->numRows; iy++) {
+	for (int ix = 0; ix < image->numCols; ix++) {
+	    dx = ix + image->col0 - x;
+	    dy = iy + image->row0 - y;
+	    r2 = PS_SQR(dx) + PS_SQR(dy);
+	    if (r2 < R2) continue;
+	    if (logical_and) {
+		image->data.U8[iy][ix] &= maskValue;
+	    } else {
+		image->data.U8[iy][ix] |= maskValue;
+	    }
 	}
     }
Index: /trunk/psphot/src/psphot.c
===================================================================
--- /trunk/psphot/src/psphot.c	(revision 4215)
+++ /trunk/psphot/src/psphot.c	(revision 4216)
@@ -1,4 +1,3 @@
 # include "psphot.h"
-void dump_image (psImage *image, char *filename);
 
 int main (int argc, char **argv) {
@@ -32,16 +31,16 @@
 
     // use stellar objects SN > PSF_SN_LIM
-    psf = choose_psf_model (image, config, sources);
+    psf = choose_psf_model (imdata, config, sources);
 
     if (0) { 
-      psImage *zap = psImageAlloc (image->numCols, image->numRows, PS_TYPE_F32);
+      psImage *zap = psImageAlloc (imdata->image->numCols, imdata->image->numRows, PS_TYPE_F32);
       psImageInit (zap, 0);
       find_defects (zap, sources, config, psf);
-      dump_image (zap, argv[2]);
+      DumpImage (zap, argv[2]);
       exit (0);
     }
 
     // test PSF on all except SATURATE and DEFECT (artifacts)
-    apply_psf_model (image, config, sources, psf, sky);
+    apply_psf_model (imdata, config, sources, psf, sky);
 
     // XXX this test compares dSx,dSy distribution with model.
@@ -49,20 +48,10 @@
     // test_psf_scatter (sources);
 
-    // fit_galaxies (image, config, sources);
+    // fit_galaxies (imdata, config, sources);
 
     // subtract_galaxies (sources, config);
  
-    dump_image (image, argv[2]);
+    DumpImage (imdata->image, argv[2]);
     DumpModelPSF (sources, "sources.dat");
     exit (0);
 }
-
-void dump_image (psImage *image, char *filename) {
-
-    unlink (filename);
-    psFits *fits = psFitsAlloc (filename);
-    psFitsWriteImage (fits, NULL, image, 0, NULL);
-    psFree (fits);
-    return;
-}
-
Index: /trunk/psphot/src/psphot.h
===================================================================
--- /trunk/psphot/src/psphot.h	(revision 4215)
+++ /trunk/psphot/src/psphot.h	(revision 4216)
@@ -36,4 +36,5 @@
 bool DumpModelFits (psArray *models, char *filename);
 bool DumpModelPSF (psArray *sources, char *filename);
+bool DumpImage (psImage *image, char *filename);
 
 int usage ();
Index: /trunk/psphot/src/pspsf.c
===================================================================
--- /trunk/psphot/src/pspsf.c	(revision 4215)
+++ /trunk/psphot/src/pspsf.c	(revision 4216)
@@ -46,5 +46,5 @@
 }
 
-pmPSF_Test *pmPSF_TestModel (psArray *sources, char *modelName) {
+pmPSF_Test *pmPSF_TestModel (psArray *sources, char *modelName, float RADIUS) {
   
     int Nflt = 0;
@@ -58,7 +58,14 @@
 	psSource *source = test->sources->data[i];
 	psModel  *model  = pmSourceModelGuess (source, test->modelType); 
-
+	x = source->peak->x;
+	y = source->peak->y;
+
+	// set temporary object mask and fit object
 	// fit model as FLT, not PSF (mask & skip poor fits)
-	if (!pmSourceFitModel (source, model, false)) {
+	psImageKeepCircle (source->mask, x, y, RADIUS, OR, 0x80);
+	status = pmSourceFitModel (source, model, false);
+	psImageKeepCircle (source->mask, x, y, RADIUS, AND, 0x7f);
+
+	if (!status) {
 	  test->mask->data.U8[i] = 1;
 	  continue;
@@ -96,6 +103,10 @@
 	psModel  *modelPSF = psModelFromPSF (modelFLT, test->psf); // set shape for this model
 
-	// fit model as PSF, not FLT (skip poor fits)
-	if (!pmSourceFitModel (source, modelPSF, true)) {
+	psImageKeepCircle (source->mask, RADIUS, OR, 0x80);
+	status = pmSourceFitModel (source, modelPSF, true);
+	psImageKeepCircle (source->mask, RADIUS, AND, 0x7f);
+
+	// skip poor fits
+	if (!status) {
 	  test->mask->data.U8[i] = 1;
 	  continue;
@@ -224,20 +235,2 @@
 }
 
-// this was an attempt to measure ap-fit using image stats.  it was giving
-// strangely wrong answers.  
-# if (0)
-	psImage  *pixels = psImageCopy (NULL, source->pixels, source->pixels->type.type);
-	pixels->col0 = source->pixels->col0;
-	pixels->row0 = source->pixels->row0;
-	// psImageCopy does not retain the parent offset
-	// this copy is relative to the parent array.
-
-	pmSourceSubModel (pixels, source->mask, model, false);
-
-	// use a robust statistic or not?
-	psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN);
-	psImageStats (stats, pixels, NULL, 0);
-	// psImageStats (stats, pixels, source->mask, 0);
-	// XXX this should be the SUM, but we don't have such a statistic
-	// XXX given the mask, we need stat functions to return Npixels used
-# endif
Index: /trunk/psphot/src/setup.c
===================================================================
--- /trunk/psphot/src/setup.c	(revision 4215)
+++ /trunk/psphot/src/setup.c	(revision 4216)
@@ -1,19 +1,31 @@
 # include "psphot.h"
+
+// load the image
+// load or construct the mask
+//   apply X/Y border info
+//   apply the SATURATE mask?
+// load or construct the noise 
 
 psImageData *setup (psMetadata *config) 
 {
+    psMetadata *header = NULL;
+    psImage *image = NULL;
+    psImage *noise = NULL;
+    psImage *mask = NULL;
+
     psRegion  region = {0,0,0,0};	// a region representing the entire array
 
     psTimerStart ("psphot");
 
+    // setup header template, specify COMMENT and HISTORY as multis
+    header = psMetadataAlloc ();
+    psMetadataAdd (header, PS_LIST_HEAD, "COMMENT", PS_META_MULTI, "folder for comment", NULL);
+    psMetadataAdd (header, PS_LIST_HEAD, "HISTORY", PS_META_MULTI, "folder for history", NULL);
+
+    // read header and image data from INPUT 
     char *input = psMetadataLookupSTR (&status, config, "INPUT");
     psFits *file = psFitsAlloc (input);
-
-    psImageData *imdata = psAlloc(sizeof(psImageData));
-    psMetadata *header = psMetadataAlloc ();
-    psMetadataAdd (header, PS_LIST_HEAD, "COMMENT", PS_META_MULTI, "folder for comment", NULL);
-    psMetadataAdd (header, PS_LIST_HEAD, "HISTORY", PS_META_MULTI, "folder for history", NULL);
     header = psFitsReadHeader (header, file);
-    psImage *image = psFitsReadImage (NULL, file, region, 0);
+    image = psFitsReadImage (NULL, file, region, 0);
     psFree (file);
 
@@ -25,20 +37,17 @@
     psMetadataAdd (config, PS_LIST_TAIL, "GAIN",    PS_META_F32 | PS_META_REPLACE, "gain value used", GAIN);
 
-    imdata->image = image;
-    imdata->header = header;
-
     // load the noise image if it is supplied, otherwise build from input
     char *noiseName = psMetadataLookupSTR (&status, config, "NOISE");
     if (status == true) {
       file = psFitsAlloc (noiseName);
-      imdata->noise = psFitsReadImage  (NULL, file, region, 0);
+      noise = psFitsReadImage  (NULL, file, region, 0);
       psFree (file);
     } else {
       psScalar *value = psScalarAlloc (1.0 / GAIN, PS_TYPE_F64);
-      imdata->noise = BinaryOp (NULL, imdata->image, '/', value);
+      noise = BinaryOp (NULL, image, '/', value);
       psFree (value);
 
       psScalar *value = psScalarAlloc (PS_SQR(RDNOISE/GAIN), PS_TYPE_F64);
-      imdata->noise = BinaryOp (imdata->noise, imdata->noise, '+', value);
+      noise = BinaryOp (noise, noise, '+', value);
       psFree (value);
     }
@@ -47,9 +56,9 @@
     if (status == true) {
       file = psFitsAlloc (maskName);
-      imdata->mask  = psFitsReadImage  (NULL, file, region, 0);
-      // check if U8 & force?
+      mask  = psFitsReadImage  (NULL, file, region, 0);
+      // require U8
       psFree (file);
     } else {
-      imdata->mask = psImageAlloc (image->numCols, image->numRows, PS_TYPE_U8);
+      mask = psImageAlloc (image->numCols, image->numRows, PS_TYPE_U8);
     }
 
@@ -58,8 +67,24 @@
     psRegion *keep = psRegionAlloc (XBORDER, -XBORDER, YBORDER, -YBORDER);
     keep           = psRegionForImage (keep, image, keep);
-    psImageKeepRegion (imdata->mask, keep, 0x01);
+    psImageKeepRegion (mask, keep, 0x01);
+
+    // mask the saturated pixels
+    float SATURATE = psMetadataLookupF32 (&status, config, "SATURATE");
+    for (int i = 0; i < image->numRows; i++) {
+	for (int j = 0; j < image->numCols; j++) {
+	    if (image->data.F32[i][j] >= SATURATE) {
+		mask->data.U8[i][j] |= 0x02;
+	    }
+	}
+    }
 
     psLogMsg ("psphot", 3, "load data: %f sec\n", psTimerMark ("psphot"));
 
+    psImageData *imdata = psAlloc(sizeof(psImageData));
+    imdata->image = image;
+    imdata->header = header;
+    imdata->noise = noise;
+    imdata->mask = mask;
+
     return (true);
 }
Index: /trunk/psphot/src/source_moments.c
===================================================================
--- /trunk/psphot/src/source_moments.c	(revision 4215)
+++ /trunk/psphot/src/source_moments.c	(revision 4216)
@@ -17,10 +17,13 @@
     sources->n = 0;
 
-    // this should use ROBUST not SAMPLE median, but it is broken
     for (int i = 0; i < peaks->n; i++) {
-	// provide the user arguments as metadata?
-	psSource *source  = pmSourceDefinePixels (imdata, (psPeak *)peaks->data[i], OUTER);
-	if (source == NULL) continue;
+	// create a new source, add peak
+	psSource *source = pmSourceAlloc();
+	source->peak = (psPeak *)peaks->data[i];
 
+	// allocate image, noise, mask arrays for each peak (square of radius OUTER)
+	pmSourceDefinePixels (source, imdata, source->peak->x, source->peak->y, OUTER);
+
+	// this should use ROBUST not SAMPLE median, but it is broken
 	pmSourceLocalSky_EAM (source, PS_STAT_SAMPLE_MEDIAN, INNER);
 	pmSourceMoments (source, RADIUS);