Index: /branches/eam_02_branch/psphot/src/psphotMakeResiduals.c
===================================================================
--- /branches/eam_02_branch/psphot/src/psphotMakeResiduals.c	(revision 12947)
+++ /branches/eam_02_branch/psphot/src/psphotMakeResiduals.c	(revision 12948)
@@ -1,3 +1,4 @@
 # include "psphotInternal.h"
+# define ZERO_ORDER 0
 
 bool psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf) {
@@ -9,4 +10,6 @@
     psTimerStart ("residuals");
 
+    if (!psMetadataLookupBool(&status, recipe, "PSF.RESIDUALS")) return true;
+
     int xBin = psMetadataLookupS32(&status, recipe, "PSF.RESIDUALS.XBIN");
     PS_ASSERT (status, false);
@@ -14,4 +17,25 @@
     int yBin = psMetadataLookupS32(&status, recipe, "PSF.RESIDUALS.YBIN");
     PS_ASSERT (status, false);
+
+    float nSigma = psMetadataLookupF32(&status, recipe, "PSF.RESIDUALS.NSIGMA");
+    PS_ASSERT (status, false);
+
+    char *modeString = psMetadataLookupStr(&status, recipe, "PSF.RESIDUALS.INTERPOLATION");
+    PS_ASSERT (status, false);
+
+    psImageInterpolateMode mode = psImageInterpolateModeFromString (modeString);
+    if (mode == PS_INTERPOLATE_NONE) {
+	psError(PSPHOT_ERR_CONFIG, false, "invalid interpolation in psphot.config");
+	return false;
+    }
+
+    char *statString = psMetadataLookupStr(&status, recipe, "PSF.RESIDUALS.STATISTIC");
+    PS_ASSERT (status, false);
+
+    psStatsOptions statOption = psStatsOptionFromString (statString);
+    if (!statOption) {
+	psError(PSPHOT_ERR_CONFIG, false, "invalid residual statistic in psphot.config");
+	return false;
+    }
 
     // user parameters:
@@ -40,4 +64,5 @@
     psVector *yC = psVectorAllocEmpty (100, PS_TYPE_F32);
 
+    // build (DATA - MODEL) [an image] for each psf star
     psArray *input = psArrayAllocEmpty (100);
     for (int i = 0; i < sources->n; i++) {
@@ -55,5 +80,4 @@
         psImage *mask   = psImageCopy (NULL, source->mask,   PS_TYPE_U8);
         pmModelSub (image, mask, model, false, false);
-	psFree (mask);
 	
 	// re-normalize image and weight
@@ -62,7 +86,7 @@
 	psBinaryOp (weight, weight, "/", psScalarAlloc(Io*Io, PS_TYPE_F32));
 
-	psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(
-	    PS_INTERPOLATE_BILINEAR,
-	    image, weight, NULL, 0, 0.0, 0.0, 1, 0, 0.0);
+	// we will interpolate the image and weight - include the mask or not?
+	// XXX consider better values for the mask bits
+	psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(mode, image, weight, NULL, 0xff, 0.0, 0.0, 1, 2, 0.0);
 	psArrayAdd (input,  100, interp);
 
@@ -77,4 +101,5 @@
 
 	// free up the excess references 
+	psFree (mask);
 	psFree (image);
 	psFree (weight);
@@ -89,14 +114,20 @@
     psVector *fmasks  = psVectorAlloc (input->n, PS_TYPE_U8);
 
+    // statistic to use to determine baseline for clipping
     psStats *fluxClip     = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     psStats *fluxClipDef  = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
-    psStats *fluxStats    = psStatsAlloc (PS_STAT_SAMPLE_MEAN   | PS_STAT_SAMPLE_STDEV);
-    psStats *fluxStatsDef = psStatsAlloc (PS_STAT_SAMPLE_MEAN   | PS_STAT_SAMPLE_STDEV);
-
-    // build the residual image pixel-by-pixel
-    for (int oy = 0; oy < resid->image->numRows; oy++) {
+    // statistic to use to determine output flux
+    // XXX make API to convert statOption for MEAN/MEDIAN in to corresponding STDEV?
+    psStats *fluxStats    = psStatsAlloc (statOption | PS_STAT_SAMPLE_STDEV);
+    psStats *fluxStatsDef = psStatsAlloc (statOption | PS_STAT_SAMPLE_STDEV);
+
+// this section builds just the 0th order term Ro
+# if (ZERO_ORDER)
+    // build Ro = DATA - MODEL (rebinned image) pixel-by-pixel
+    for (int oy = 0; oy < resid->Ro->numRows; oy++) {
 	fprintf (stderr, ".");
-	for (int ox = 0; ox < resid->image->numCols; ox++) {
+	for (int ox = 0; ox < resid->Ro->numCols; ox++) {
 	    
+	    // build the vector of data values for this output pixel
 	    for (int i = 0; i < input->n; i++) {
 
@@ -104,6 +135,6 @@
 		
 		// fractional image position
-		float ix = (ox - resid->xCenter) / (float) xBin + xC->data.F32[i] - interp->image->col0;
-		float iy = (oy - resid->yCenter) / (float) yBin + yC->data.F32[i] - interp->image->row0;		
+		float ix = (ox + 0.5 - resid->xCenter) / (float) xBin + xC->data.F32[i] - interp->image->col0;
+		float iy = (oy + 0.5 - resid->yCenter) / (float) yBin + yC->data.F32[i] - interp->image->row0;
 
 		mflux = 0;
@@ -115,7 +146,7 @@
 	    }
 
+	    // measure the robust median to determine a baseline reference value
 	    *fluxClip = *fluxClipDef;
 	    psVectorStats (fluxClip, fluxes, NULL, fmasks, 0xff);
-	    psErrorClear(); 
 
 	    // mark input pixels which are more than N sigma from the median
@@ -123,27 +154,102 @@
 		float delta = fluxes->data.F32[i] - fluxClip->robustMedian;
 		float sigma = sqrt (dfluxes->data.F32[i]);
-		
 		float swing = fabs(delta) / sigma;
 
 		// make this a user option
-		if (swing > 3.0) {
+		if (swing > nSigma) {
 		    fmasks->data.U8[i] = 1;
 		}
 	    }		    
 
+	    // measure the desired statistic on the unclipped pixels
 	    *fluxStats = *fluxStatsDef;
 	    psVectorStats (fluxStats, fluxes, NULL, fmasks, 0xff);
+
+	    resid->Ro->data.F32[oy][ox] = psStatsGetValue (fluxStats, statOption);
+	    resid->weight->data.F32[oy][ox] = fluxStats->sampleStdev;
+
+	    // clear (ignore) any outstanding errors 
 	    psErrorClear(); 
-
-	    resid->image->data.F32[oy][ox] = fluxStats->sampleMean;
-	    resid->weight->data.F32[oy][ox] = fluxStats->sampleStdev;
 	}
     }
-
-    psphotSaveImage (NULL, resid->image, "resid.im.fits");
-    psphotSaveImage (NULL, resid->weight, "resid.wt.fits");
-    psphotSaveImage (NULL, resid->mask, "resid.mk.fits");
-
-    psLogMsg ("psphot.pspsf", PS_LOG_INFO, "generate residuals for %ld objects: %f sec\n", input->n, psTimerMark ("residuals"));
+    psLogMsg ("psphot.pspsf", PS_LOG_INFO, "generated 0-th order residuals for %ld objects: %f sec\n", input->n, psTimerMark ("residuals"));
+
+# else
+
+    psImage *A = psImageAlloc(3, 3, PS_TYPE_F64); // Least-squares matrix
+    psVector *B = psVectorAlloc(3, PS_TYPE_F64); // Least-squares vector
+
+    // build (x,y)*(DATA - MODEL - Ro) pixel-by-pixel
+    for (int oy = 0; oy < resid->Ro->numRows; oy++) {
+	fprintf (stderr, ".");
+	for (int ox = 0; ox < resid->Ro->numCols; ox++) {
+	    
+	    // build the vector of data values for this output pixel
+	    // XXX this is identical to the pass above: we could cache the results for speed
+	    for (int i = 0; i < input->n; i++) {
+
+		psImageInterpolateOptions *interp = input->data[i];
+		
+		// fractional image position
+		float ix = (ox + 0.5 - resid->xCenter) / (float) xBin + xC->data.F32[i] - interp->image->col0;
+		float iy = (oy + 0.5 - resid->yCenter) / (float) yBin + yC->data.F32[i] - interp->image->row0;
+
+		mflux = 0;
+		psImageInterpolate (&flux, &dflux, &mflux, ix, iy, interp);
+		fluxes->data.F32[i] = flux;
+		dfluxes->data.F32[i] = dflux;
+		fmasks->data.U8[i] = mflux;
+		// fprintf (stderr, "%f %f : %f %f (%d)\n", ix, iy, flux, dflux, fmasks->data.U8[i]);
+	    }
+
+	    // measure the robust median to determine a baseline reference value
+	    *fluxClip = *fluxClipDef;
+	    psVectorStats (fluxClip, fluxes, NULL, fmasks, 0xff);
+	    psErrorClear();		// clear (ignore) any outstanding errors 
+
+	    // mark input pixels which are more than N sigma from the median
+	    for (int i = 0; i < fluxes->n; i++) {
+		float delta = fluxes->data.F32[i] - fluxClip->robustMedian;
+		float sigma = sqrt (dfluxes->data.F32[i]);
+		float swing = fabs(delta) / sigma;
+
+		// make this a user option
+		if (swing > nSigma) {
+		    fmasks->data.U8[i] = 1;
+		}
+	    }		    
+
+	    psImageInit(A, 0.0);
+	    psVectorInit(B, 0.0);
+	    for (int i = 0; i < fluxes->n; i++) {
+		if (fmasks->data.U8[i]) continue;
+		B->data.F64[0] += fluxes->data.F32[i]/dfluxes->data.F32[i];
+		B->data.F64[1] += fluxes->data.F32[i]*xC->data.F32[i]/dfluxes->data.F32[i];
+		B->data.F64[2] += fluxes->data.F32[i]*yC->data.F32[i]/dfluxes->data.F32[i];
+
+		A->data.F64[0][0] += 1.0/dfluxes->data.F32[i];
+		A->data.F64[1][0] += xC->data.F32[i]/dfluxes->data.F32[i];
+		A->data.F64[2][0] += yC->data.F32[i]/dfluxes->data.F32[i];
+
+		A->data.F64[1][1] += PS_SQR(xC->data.F32[i])/dfluxes->data.F32[i];
+		A->data.F64[2][2] += PS_SQR(yC->data.F32[i])/dfluxes->data.F32[i];
+		A->data.F64[1][2] += xC->data.F32[i]*yC->data.F32[i]/dfluxes->data.F32[i];
+	    }
+
+	    A->data.F64[0][1] = A->data.F64[1][0];
+	    A->data.F64[0][2] = A->data.F64[2][0];
+	    A->data.F64[2][1] = A->data.F64[1][2];
+	    psMatrixGJSolve(A, B);
+
+	    resid->Ro->data.F32[oy][ox] = B->data.F64[0];
+	    resid->Rx->data.F32[oy][ox] = B->data.F64[1];
+	    resid->Ry->data.F32[oy][ox] = B->data.F64[2];
+	}
+    }
+
+    psFree (A);
+    psFree (B);
+
+# endif
 
     psFree (xC);
@@ -160,4 +266,12 @@
     psFree (fluxClipDef);
 
+    psphotSaveImage (NULL, resid->Ro,     "resid.ro.fits");
+    psphotSaveImage (NULL, resid->Rx,     "resid.rx.fits");
+    psphotSaveImage (NULL, resid->Ry,     "resid.ry.fits");
+    psphotSaveImage (NULL, resid->weight, "resid.wt.fits");
+    psphotSaveImage (NULL, resid->mask,   "resid.mk.fits");
+
+    psLogMsg ("psphot.pspsf", PS_LOG_INFO, "generate residuals for %ld objects: %f sec\n", input->n, psTimerMark ("residuals"));
+
     psf->residuals = resid;
     return true;