Index: /branches/eam_branches/ipp-20120805/psphot/src/psphotBlendFit.c
===================================================================
--- /branches/eam_branches/ipp-20120805/psphot/src/psphotBlendFit.c	(revision 34379)
+++ /branches/eam_branches/ipp-20120805/psphot/src/psphotBlendFit.c	(revision 34380)
@@ -1,3 +1,4 @@
 # include "psphotInternal.h"
+bool psphotBlendFitSetSource(pmSource *source, psImage *IDimage, float sigSigma);
 
 // for now, let's store the detections on the readout->analysis for each readout
@@ -92,4 +93,15 @@
     float skySig = psMetadataLookupF32(&status, recipe, "SKY_SIG");
     assert (status && isfinite(skySig) && skySig > 0);
+
+    // **** test block : generate an ID image where pixels are set based the source models (flux > 0.1 peak && flux > 2.0 skySig)
+    psImage *IDimage = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_S32);
+    psImageInit (IDimage, 0);
+
+    // start with the currently known moments (Mxx, Myy, Mxy) and generate a window image
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = sources->data[i];
+	psphotBlendFitSetSource (source, IDimage, skySig);
+    }
+    psphotSaveImage (NULL, IDimage, "idimage.fits");
 
     // Define source fitting parameters for extended source fits
@@ -359,2 +371,82 @@
 }
 
+bool psphotBlendFitSetSource(pmSource *source, psImage *IDimage, float skySigma) {
+
+    if (!source) return false;
+    if (!source->peak) return false; // XXX how can we have a peak-less source?
+    if (!source->moments) return false;
+    if (source->type == PM_SOURCE_TYPE_DEFECT) return false;
+    if (source->type == PM_SOURCE_TYPE_SATURATED) return false;
+    if (source->mode2 == PM_SOURCE_MODE2_MATCHED) return false;
+    psAssert(IDimage, "need a window");
+
+    if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) {
+	// find the radius at which we hit something like 0.1*SATURATION
+	return false;
+    }
+
+    if (!isfinite(source->moments->Mrf) || source->moments->Mrf < 0 ) return false;
+
+    // we have a source with moments Mx, My, Mxx, Mxy, Myy.  we just need to define a Gaussian that has 
+    // these values (and peak of 1.0)
+
+    int Nx = IDimage->numCols;
+    int Ny = IDimage->numRows;
+
+    float Xo = source->moments->Mx;
+    float Yo = source->moments->My;
+
+    psEllipseMoments moments;
+    moments.x2 = source->moments->Mxx;
+    moments.y2 = source->moments->Myy;
+    moments.xy = source->moments->Mxy;
+
+    psEllipseAxes axes = psEllipseMomentsToAxes(moments, 20.0);
+    if (! (isfinite(axes.major) && isfinite(axes.minor) && isfinite(axes.theta)) ) {
+        fprintf (stderr, "invalid axes found id: %4d major: %f minor: %f theta: %f\n", source->id, axes.major, axes.minor, axes.theta);
+        return false;
+    }
+
+    // Use 1st radial moment as size, not sigma.  Why this factor of 0.5 ?
+    float scale = 0.5 * source->moments->Mrf / axes.major;
+    axes.major *= scale;
+    axes.minor *= scale;
+
+    psEllipseShape shape = psEllipseAxesToShape(axes);
+    if (! (isfinite(shape.sx) && isfinite(shape.sy) && isfinite(shape.sxy)) ) {
+        fprintf (stderr, "invalid shape found id: %d sx: %f sy %f sxy: %f\n", source->id, shape.sx, shape.sy, shape.sxy);
+        return false;
+    }
+
+    float Smajor = axes.major;
+
+    // we want to go to 2.15 sigma (0.1 Io)
+    int minX = PS_MIN(PS_MAX(Xo - 2.15*Smajor, 0), Nx - 1);
+    int maxX = PS_MIN(PS_MAX(Xo + 2.15*Smajor, 0), Nx - 1);
+    int minY = PS_MIN(PS_MAX(Yo - 2.15*Smajor, 0), Ny - 1);
+    int maxY = PS_MIN(PS_MAX(Yo + 2.15*Smajor, 0), Ny - 1);
+
+    float rMxx = 0.5 / PS_SQR(shape.sx);
+    float rMyy = 0.5 / PS_SQR(shape.sy);
+    float Sxy = -1. * shape.sxy;    // factor of -1 is included to match the previous window function
+                                    // implementation. XXX: Is this correct?
+
+    int ID = source->id;
+    float Io = source->peak->rawFlux;
+    for (int iy = minY; iy < maxY; iy++) {
+	for (int ix = minX; ix < maxX; ix++) {
+
+	    float dX = (ix + 0.5 - Xo);
+	    float dY = (iy + 0.5 - Yo);
+
+	    float z = rMxx * PS_SQR(dX) + rMyy * PS_SQR(dY) + Sxy*dX*dY;
+
+	    float f = Io*exp(-z);
+	    
+	    if (f < 2.0*skySigma) continue;
+            IDimage->data.S32[iy][ix] = ID;
+	}
+    }
+
+    return true;
+}