Index: trunk/psphot/src/psphotEnsemblePSF.c
===================================================================
--- trunk/psphot/src/psphotEnsemblePSF.c	(revision 6427)
+++ trunk/psphot/src/psphotEnsemblePSF.c	(revision 6481)
@@ -1,7 +1,8 @@
 # include "psphot.h"
+bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight);
 
 // 2006.02.07 : no leaks!
 // fit all reasonable sources with the linear PSF model
-bool psphotEnsemblePSF (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf) { 
+bool psphotEnsemblePSF (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf, bool final) { 
 
     bool  status;
@@ -9,4 +10,7 @@
     float y;
     float f;
+    float r;
+
+    // psRegion allArray = psRegionSet (0, 0, 0, 0);
 
     psTimerStart ("psphot");
@@ -29,14 +33,9 @@
 
     // option to limit analysis to a specific region
-    bool UseAnalysisRegion = false;
-    psRegion AnalysisRegion = {0, 0, 0, 0};
     char *region = psMetadataLookupStr (&status, config, "ANALYSIS_REGION");
-    if (status) {
-	UseAnalysisRegion = true;
-	AnalysisRegion = psRegionFromString (region);
-	psLogMsg ("psphotEnsemblePSF", 4, "using region %f,%f - %f,%f\n", 
-		  AnalysisRegion.x0, AnalysisRegion.y0, 
-		  AnalysisRegion.x1, AnalysisRegion.y1);
-    }
+    psRegion AnalysisRegion = psRegionFromString (region);
+    AnalysisRegion = psRegionForImage (readout->image, AnalysisRegion);
+    // psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
+    if (psRegionIsNaN (AnalysisRegion)) psAbort ("psphot", "analysis region mis-defined");
 
     for (int i = 0; i < sources->n; i++) {
@@ -45,14 +44,16 @@
 	// skip non-astronomical objects (very likely defects)
 	// XXX EAM : should we try these anyway?
-	if (inSource->mode &  PM_SOURCE_BLEND) continue;
 	if (inSource->type == PM_SOURCE_DEFECT) continue; 
 	if (inSource->type == PM_SOURCE_SATURATED) continue;
-
-	if (UseAnalysisRegion) {
-	    if (inSource->moments->x < AnalysisRegion.x0) continue;
-	    if (inSource->moments->y < AnalysisRegion.y0) continue;
-	    if (inSource->moments->x > AnalysisRegion.x1) continue;
-	    if (inSource->moments->y > AnalysisRegion.y1) continue;
-	}
+	if (final) {
+	    if (inSource->mode &  PM_SOURCE_SUBTRACTED) continue;
+	} else {
+	    if (inSource->mode &  PM_SOURCE_BLEND) continue;
+	}
+
+	if (inSource->moments->x < AnalysisRegion.x0) continue;
+	if (inSource->moments->y < AnalysisRegion.y0) continue;
+	if (inSource->moments->x > AnalysisRegion.x1) continue;
+	if (inSource->moments->y > AnalysisRegion.y1) continue;
 
 	pmSource *otSource = pmSourceAlloc ();
@@ -73,6 +74,7 @@
 	    // peak-up on peak (for non-sat objects)
 
-	    int ix = inSource->peak->x;
-	    int iy = inSource->peak->y;
+	    // ix,iy must land on inSource->pixels 
+	    int ix = PS_MAX (inSource->pixels->col0 + 1, PS_MIN (inSource->pixels->col0 + inSource->pixels->numCols - 2, inSource->peak->x));
+	    int iy = PS_MAX (inSource->pixels->row0 + 1, PS_MIN (inSource->pixels->row0 + inSource->pixels->numRows - 2, inSource->peak->y));
 
 	    psPolynomial2D *bicube = psImageBicubeFit (inSource->pixels, ix, iy);
@@ -104,4 +106,5 @@
 	otSource->mask   = psImageCopy (NULL, inSource->mask,   PS_TYPE_U8);
 	otSource->pixels = psImageCopy (NULL, inSource->pixels, PS_TYPE_F32);
+	otSource->weight = psImageCopy (NULL, inSource->weight, PS_TYPE_F32);
 
 	// build the model image 
@@ -117,4 +120,8 @@
 	psImageKeepCircle (mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED);
 	pmSourceAddModel (flux, mask, model, false, false);
+
+	// calculate nDOF (nPix - 1)
+	// int Nmaskpix = psImageCountPixelMask (mask, allArray, PSPHOT_MASK_SATURATED);
+	// model->nDOF  = mask->numCols*mask->numRows - Nmaskpix - 1;
 
 	// save source in list
@@ -128,4 +135,5 @@
     // fill out the sparse matrix
     psSparse *sparse = psSparseAlloc (models->n, 100);
+    psVector *weight = psVectorAlloc (models->n, PS_TYPE_F32);
     for (int i = 0; i < models->n; i++) {
 	int N = index->data.U32[i];
@@ -133,11 +141,13 @@
 	pmSource *Mi = models->data[i];
 
-	// diagonal element (auto-cross-product)
-	f = pmSourceCrossProduct (Mi, Mi);
-	psSparseMatrixElement (sparse, i, i, f);
+	// scale by diagonal element (auto-cross-product)
+	r = pmSourceCrossProduct (Mi, Mi);
+	weight->data.F32[i] = r;
+
+	psSparseMatrixElement (sparse, i, i, 1.0);
 
 	// find the image x model value
 	f = pmSourceCrossProduct (Fi, Mi);
-	psSparseVectorElement (sparse, i, f);
+	psSparseVectorElement (sparse, i, f / r);
 
 	// loop over all other stars following this one
@@ -153,12 +163,17 @@
 	    // got an overlap; calculate cross-product and add to output array
 	    f = pmSourceCrossProduct (Mi, Mj);
-	    psSparseMatrixElement (sparse, j, i, f); 
+	    psSparseMatrixElement (sparse, j, i, f / r); 
 	}
     }
     psLogMsg ("psphot.emsemble", 4, "built matrix: %f (%d elements)\n", psTimerMark ("psphot"), sparse->Nelem);
+
+    psSparseConstraint constraint;
+    constraint.paramMin   = 0;
+    constraint.paramMax   = 1e8;
+    constraint.paramDelta = 1e8;
 
     // solve for normalization terms (need include local sky?)
     psSparseResort (sparse);
-    psVector *norm = psSparseSolve (NULL, sparse, 3);
+    psVector *norm = psSparseSolve (NULL, constraint, sparse, 3);
 
     // adjust models, set sources and subtract
@@ -168,4 +183,7 @@
 	pmSource *Mi = models->data[i];
 
+	// if we already have a PSF model, free it.
+	psFree (Fi->modelPSF);
+
 	// need to increment counter so we can free models here and sources above
 	Fi->modelPSF = psMemCopy (Mi->modelPSF);
@@ -173,16 +191,52 @@
 	// assign linearly-fitted normalization
 	Fi->modelPSF->params->data.F32[1] = norm->data.F32[i];
+	Fi->modelPSF->dparams->data.F32[1] = sqrt(sqrt(2) * norm->data.F32[i] / (sparse->Bfj->data.F32[i] * weight->data.F32[i]));
+	// XXX EAM : this factor of sqrt(2) makes the errors consistent, but I don't understand it
 
 	// subtract object
 	pmSourceSubModel (Fi->pixels, Fi->mask, Fi->modelPSF, false, false);
 	Fi->mode |= PM_SOURCE_SUBTRACTED;
-	Fi->mode |= PM_SOURCE_TEMPSUB;
-    }
+	if (!final) Fi->mode |= PM_SOURCE_TEMPSUB;
+    }
+
+    // measure chisq for each source
+    for (int i = 0; final && (i < models->n); i++) {
+	int N = index->data.U32[i];
+	pmSource *Fi = sources->data[N];
+	pmModel *model = Fi->modelPSF;
+
+	x = model->params->data.F32[2];
+	y = model->params->data.F32[3];
+
+	psImageKeepCircle (Fi->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED);
+	pmSourceChisq (model, Fi->pixels, Fi->mask, Fi->weight);
+	psImageKeepCircle (Fi->mask, x, y, model->radius, "AND", ~PSPHOT_MASK_MARKED);
+    }
+
     psFree (index);
     psFree (sparse);
     psFree (models);
     psFree (norm);
+    psFree (weight);
 
     psLogMsg ("psphot.emsemble", 3, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot"));
     return true;
 }
+
+bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight) {
+
+    double dC = 0.0;
+    int Npix = 0;
+    for (int j = 0; j < image->numRows; j++) {
+	for (int i = 0; i < image->numCols; i++) {
+	    if (mask->data.U8[j][i]) continue;
+	    if (weight->data.F32[j][i] <= 0) continue;
+	    dC += PS_SQR (image->data.F32[j][i]) / weight->data.F32[j][i];
+	    Npix ++;
+	}
+    }
+    model->nDOF = Npix - 1;
+    model->chisq = dC;
+
+    return (true);
+}