Index: /branches/eam_branch_20070830/psLib/src/imageops/psImageMap.c
===================================================================
--- /branches/eam_branch_20070830/psLib/src/imageops/psImageMap.c	(revision 14863)
+++ /branches/eam_branch_20070830/psLib/src/imageops/psImageMap.c	(revision 14864)
@@ -7,6 +7,6 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.1.2.5 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-07 20:19:35 $
+ *  @version $Revision: 1.1.2.6 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-17 20:17:21 $
  *
  *  Copyright 2007 Institute for Astronomy, University of Hawaii
@@ -39,5 +39,6 @@
 
     psFree (map->map);
-    // psFree (map->field); XXX ??? this should be freed here, but that causes an error...
+    psFree (map->error);
+    psFree (map->field); // XXX ??? this should be freed here, but that causes an error...
     psFree (map->stats);
     psFree (map->binning);
@@ -57,5 +58,10 @@
     map->field   = psMemIncrRefCounter (field);
     map->stats   = psMemIncrRefCounter (stats);
+
     map->map     = psImageAlloc (binning->nXruff, binning->nYruff, PS_TYPE_F32);
+    psImageInit (map->map, 0.0);
+
+    map->error   = psImageAlloc (binning->nXruff, binning->nYruff, PS_TYPE_F32);
+    psImageInit (map->error, 0.0);
 
     psImageBinningSetScale (map->binning, PS_IMAGE_BINNING_CENTER);
@@ -73,4 +79,5 @@
 
     psImageRecycle (map->map, nXruff, nYruff, PS_TYPE_F32);
+    psImageRecycle (map->error, nXruff, nYruff, PS_TYPE_F32);
 
     psImageBinningSetScale (map->binning, PS_IMAGE_BINNING_CENTER);
@@ -257,9 +264,10 @@
 }
 
+// x,y are in fractional pixel coords of the fine image (pixel center: 0.5)
 double psImageMapEval (psImageMap *map, float x, float y) {
 
     double result;
 
-    result = psImageUnbinPixel(x, y, map->map, map->binning);
+    result = psImageUnbinPixel_V2(x, y, map->map, map->binning);
 
     return result; 
@@ -276,5 +284,5 @@
 
     for (int i = 0; i < x->n; i++) {
-	result->data.F32[i] = psImageUnbinPixel(x->data.F32[i], y->data.F32[i], map->map, map->binning);
+	result->data.F32[i] = psImageUnbinPixel_V2(x->data.F32[i], y->data.F32[i], map->map, map->binning);
     }
 
Index: /branches/eam_branch_20070830/psLib/src/imageops/psImageMap.h
===================================================================
--- /branches/eam_branch_20070830/psLib/src/imageops/psImageMap.h	(revision 14863)
+++ /branches/eam_branch_20070830/psLib/src/imageops/psImageMap.h	(revision 14864)
@@ -7,6 +7,6 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.1.2.4 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-15 19:43:43 $
+ *  @version $Revision: 1.1.2.5 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-17 20:17:21 $
  *
  *  Copyright 2007 Institute for Astronomy, University of Hawaii
@@ -26,4 +26,5 @@
     psStats *stats;
     psImage *map;
+    psImage *error;
     psImage *field;
     psImageBinning *binning;
Index: /branches/eam_branch_20070830/psLib/src/imageops/psImageMapFit.c
===================================================================
--- /branches/eam_branch_20070830/psLib/src/imageops/psImageMapFit.c	(revision 14863)
+++ /branches/eam_branch_20070830/psLib/src/imageops/psImageMapFit.c	(revision 14864)
@@ -7,6 +7,6 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-17 01:14:52 $
+ *  @version $Revision: 1.1.2.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-17 20:17:21 $
  *
  *  Copyright 2007 Institute for Astronomy, University of Hawaii
@@ -22,6 +22,6 @@
 
 // XXX for testing
-// #include "psFits.h"
-// #include "psFitsImage.h"
+#include "psFits.h"
+#include "psFitsImage.h"
 
 #include "psAssert.h"
@@ -238,24 +238,32 @@
     // fprintf (stderr, "Total: %f\n", Total);
 
-    // write 1s in unused diagonal elements
-    // XXX for test only
-    # if (0)
-    for (int n = 0; n < Nx; n+=Nx-1) {
-	for (int m = 0; m < Ny; m++) {
-	    I = n + Nx * m;
-	    A->data.F32[I][I] = 1.0;
-	}
-    }
-    for (int n = 0; n < Nx; n++) {
-	for (int m = 0; m < Ny; m+=Ny-1) {
-	    I = n + Nx * m;
-	    A->data.F32[I][I] = 1.0;
-	}
-    }
-    # endif
-
-    // psFits *fits = psFitsOpen ("Agj.fits", "w");
-    // psFitsWriteImage (fits, NULL, A, 0, NULL);
-    // psFitsClose (fits);
+    // test for empty diagonal elements (unconstained cells), mark, and set pivots to 1.0
+    psVector *Empty = psVectorAlloc (Nx*Ny, PS_TYPE_S8);
+    psVectorInit (Empty, 0);
+    for (int i = 0; i < Nx*Ny; i++) {
+	if (A->data.F32[i][i] == 0.0) {
+	    Empty->data.S8[i] = 1;
+	    for (int j = 0; j < Nx*Ny; j++) {
+		A->data.F32[i][j] = 0.0;
+		A->data.F32[j][i] = 0.0;
+	    }
+	    A->data.F32[i][i] = 1.0;
+	    B->data.F32[i] = 0.0;
+	}
+    }
+
+    psFits *fits = psFitsOpen ("Agj.fits", "w");
+    psFitsWriteImage (fits, NULL, A, 0, NULL);
+    psFitsClose (fits);
+
+    psImage *vector = psImageAlloc (1, B->n, PS_TYPE_F32);
+    for (int n = 0; n < B->n; n++) {
+	vector->data.F32[0][n] = B->data.F32[n];
+    }
+
+    fits = psFitsOpen ("Bgj.fits", "w");
+    psFitsWriteImage (fits, NULL, vector, 0, NULL);
+    psFitsClose (fits);
+    psFree (vector);
 
     if (!psMatrixGJSolveF32(A, B)) {
@@ -267,8 +275,17 @@
     }
     
+    // set bad values to NaN
+    for (int i = 0; i < Nx*Ny; i++) {
+	if (Empty->data.S8[i]) {
+	    B->data.F32[i] = NAN;
+	}
+    }
+
+
     for (int n = 0; n < Nx; n++) {
 	for (int m = 0; m < Ny; m++) {
     	    I = n + Nx * m;
 	    map->map->data.F32[m][n] = B->data.F32[I];
+	    map->error->data.F32[m][n] = sqrt(A->data.F32[I][I]);
 	}
     }
@@ -276,4 +293,5 @@
     psFree (A);
     psFree (B);
+    psFree (Empty);
 
     return true;
Index: /branches/eam_branch_20070830/psLib/test/imageops/tap_psImageMapFit.c
===================================================================
--- /branches/eam_branch_20070830/psLib/test/imageops/tap_psImageMapFit.c	(revision 14863)
+++ /branches/eam_branch_20070830/psLib/test/imageops/tap_psImageMapFit.c	(revision 14864)
@@ -25,11 +25,87 @@
 # define TEST_9PT_3 1
 # define TEST_9PT_4 1
+# define TEST_9PT_5 1
 
 int main (void)
 {
 
-    plan_tests(1);
+    plan_tests(219);
     
     // *** tests to demonstrate the validity of the algorithm or concept ***
+
+    // test with unconstrained cell: 3x3 grid fitted to 8 points with simple slope and scale difference
+    # if (TEST_9PT_5)
+    {
+        psMemId id = psMemGetId();
+
+	psImageBinning *binning = psImageBinningAlloc();
+	binning->nXfine = 6;
+	binning->nYfine = 6;
+	binning->nXruff = 3;
+	binning->nYruff = 3;
+
+	// generate a grid of test data points
+	psVector *x = psVectorAlloc (50, PS_TYPE_F32);
+	psVector *y = psVectorAlloc (50, PS_TYPE_F32);
+	psVector *f = psVectorAlloc (50, PS_TYPE_F32);
+
+	// the underlying field is f = ix + iy, where ix,iy are fine pixel coordinates
+	// place the measurement points exactly on the ruff reference pixel centers
+	int n = 0;
+	for (float ix = 1.0; ix < 6.0; ix += 2.0) {
+	    for (float iy = 1.0; iy < 6.0; iy += 2.0) {
+		if ((ix == 1.0) && (iy == 1.0)) continue;
+		x->data.F32[n] = ix;
+		y->data.F32[n] = iy;
+		f->data.F32[n] = x->data.F32[n] + y->data.F32[n];
+		n++;
+	    }
+	}
+	x->n = n;
+	y->n = n;
+	f->n = n;
+
+	psImage *field = psImageAlloc(6, 6, PS_TYPE_F32);
+	for (int ix = 0; ix < 6; ix++) {
+	    for (int iy = 0; iy < 6; iy++) {
+		field->data.F32[iy][ix] = (ix + 0.5) + (iy + 0.5);
+	    }
+	}
+
+	psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+
+	// scale defines both field and map image sizes (nXfine, nXruff)
+	psImageMap *map = psImageMapAlloc (NULL, binning, stats);
+
+	// fit the data to the map
+	psImageMapFit (map, x, y, f, NULL);
+
+	psImage *model = psImageAlloc(field->numCols, field->numRows, PS_TYPE_F32);
+	for (int ix = 0; ix < model->numCols; ix++) {
+	    for (int iy = 0; iy < model->numRows; iy++) {
+		model->data.F32[iy][ix] = psImageUnbinPixel_V2 (ix + 0.5, iy + 0.5, map->map, map->binning);
+		if ((ix < 3) && (iy < 3)) {
+		    is_float (model->data.F32[iy][ix], NAN, "model matches expected NaN");
+		} else {
+		    is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], 1e-5, "model matches inputs");
+		}
+	    }
+	}
+
+	SaveImage (NULL, map->map, "map.fits");
+	SaveImage (NULL, field, "field.fits");
+	SaveImage (NULL, model, "model.fits");
+
+	psFree (model);
+	psFree (binning);
+	psFree (map);
+	psFree (stats);
+	psFree (field);
+	psFree (x);
+	psFree (y);
+	psFree (f);
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    # endif
 
     // test for more points: 3x3 grid fitted to 9 points with simple slope and scale difference