Index: trunk/psModules/src/objects/pmSource.c
===================================================================
--- trunk/psModules/src/objects/pmSource.c	(revision 24874)
+++ trunk/psModules/src/objects/pmSource.c	(revision 25754)
@@ -282,6 +282,6 @@
 
     psArray *peaks  = NULL;
-    pmPSFClump errorClump = {-1.0, -1.0, 0.0, 0.0};
-    pmPSFClump emptyClump = {+0.0, +0.0, 0.0, 0.0};
+    pmPSFClump errorClump = {-1.0, -1.0, 0.0, 0.0, 0, 0.0};
+    pmPSFClump emptyClump = {+0.0, +0.0, 0.0, 0.0, 0, 0.0};
     pmPSFClump psfClump;
 
@@ -319,5 +319,6 @@
 
         // construct a sigma-plane image
-        int numCols = SX_MAX / PSF_CLUMP_GRID_SCALE, numRows = SY_MAX / PSF_CLUMP_GRID_SCALE; // Size of sigma-plane image
+        int numCols = 1 + SX_MAX / PSF_CLUMP_GRID_SCALE; // Size of sigma-plane image
+	int numRows = 1 + SY_MAX / PSF_CLUMP_GRID_SCALE; // Size of sigma-plane image
         psTrace("psModules.objects", 10, "sigma-plane dimensions: %dx%d\n", numCols, numRows);
         psImage *splane = psImageAlloc(numCols, numRows, PS_TYPE_F32); // sigma-plane image
@@ -332,12 +333,16 @@
             }
 
-            int x = src->peak->x, y = src->peak->y; // Coordinates of peak
-            if (x < region->x0 || x > region->x1 || y < region->y0 || y > region->y1) {
-                continue;
-            }
+	    if (region) {
+		int x = src->peak->x, y = src->peak->y; // Coordinates of peak
+		if (x < region->x0 || x > region->x1 || y < region->y0 || y > region->y1) {
+		    continue;
+		}
+	    }
 
             if (src->mode & PM_SOURCE_MODE_BLEND) {
                 continue;
             }
+
+            if (!src->moments->nPixels) continue;
 
             if (src->moments->SN < PSF_SN_LIM) {
@@ -384,5 +389,5 @@
 
         // find the peak in this image
-        psStats *stats = psStatsAlloc (PS_STAT_MAX);
+        psStats *stats = psStatsAlloc (PS_STAT_MAX | PS_STAT_SAMPLE_STDEV);
         if (!psImageStats (stats, splane, NULL, 0)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to get image statistics.\n");
@@ -393,4 +398,6 @@
         peaks = pmPeaksInImage (splane, stats[0].max / 2);
         psTrace ("psModules.objects", 2, "clump threshold is %f\n", stats[0].max/2);
+
+	psfClump.nSigma = stats->sampleStdev;
 
         const bool keep_psf_clump = psMetadataLookupBool(NULL, recipe, "KEEP_PSF_CLUMP");
@@ -430,4 +437,9 @@
         psTrace ("psModules.objects", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->value);
 
+	// XXX store the mean sigma?
+	float meanSigma = psfClump.nSigma;
+	psfClump.nStars = clump->value;
+	psfClump.nSigma = clump->value / meanSigma;
+
         // define section window for clump
         minSx = clump->x * PSF_CLUMP_GRID_SCALE - 2.0*PSF_CLUMP_GRID_SCALE;
@@ -452,8 +464,10 @@
                 continue;
 
-            if (tmpSrc->peak->x < region->x0) continue;
-            if (tmpSrc->peak->x > region->x1) continue;
-            if (tmpSrc->peak->y < region->y0) continue;
-            if (tmpSrc->peak->y > region->y1) continue;
+	    if (region) {
+		if (tmpSrc->peak->x < region->x0) continue;
+		if (tmpSrc->peak->x > region->x1) continue;
+		if (tmpSrc->peak->y < region->y0) continue;
+		if (tmpSrc->peak->y > region->y1) continue;
+	    }
 
             if (tmpSrc->moments->Mxx < minSx)
@@ -919,4 +933,6 @@
     if (model == NULL) return false;  // model must be defined
 
+    bool addNoise = mode & PM_MODEL_OP_NOISE;
+
     if (source->modelFlux) {
         // add in the pixels from the modelFlux image
@@ -940,10 +956,9 @@
 
         psF32 **target = source->pixels->data.F32;
-        if (mode & PM_MODEL_OP_NOISE) {
-            // XXX need to scale by the gain...
+        if (addNoise) {
+	    // when adding noise, we assume the shape and Io have been modified
             target = source->variance->data.F32;
         }
 
-        // XXX need to respect the source and model masks?
         for (int iy = 0; iy < source->modelFlux->numRows; iy++) {
             int oy = iy + dY;
@@ -959,16 +974,27 @@
             }
         }
+	if (!addNoise) {
+	    if (add) {
+		source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
+	    } else {
+		source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
+	    }
+	}
         return true;
     }
 
     psImage *target = source->pixels;
-    if (mode & PM_MODEL_OP_NOISE) {
+    if (addNoise) {
         target = source->variance;
     }
 
-    if (add) {
-        status = pmModelAddWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
-    } else {
-        status = pmModelSubWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
+    if (!addNoise) {
+	if (add) {
+	    status = pmModelAddWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
+	    source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
+	} else {
+	    status = pmModelSubWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
+	    source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
+	}
     }