Index: /branches/eam_branch_20070830/psLib/src/imageops/Makefile.am
===================================================================
--- /branches/eam_branch_20070830/psLib/src/imageops/Makefile.am	(revision 14858)
+++ /branches/eam_branch_20070830/psLib/src/imageops/Makefile.am	(revision 14859)
@@ -16,4 +16,5 @@
 	psImageBinning.c \
 	psImageMap.c \
+	psImageMapFit.c \
 	psImagePixelInterpolate.c \
 	psImageUnbin.c
Index: /branches/eam_branch_20070830/psLib/src/imageops/psImageBinning.c
===================================================================
--- /branches/eam_branch_20070830/psLib/src/imageops/psImageBinning.c	(revision 14858)
+++ /branches/eam_branch_20070830/psLib/src/imageops/psImageBinning.c	(revision 14859)
@@ -8,6 +8,6 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.2.8.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-07 20:19:02 $
+ *  @version $Revision: 1.2.8.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-15 19:43:43 $
  *
  *  Copyright 2007 Institute for Astronomy, University of Hawaii
@@ -148,20 +148,44 @@
 
 // convert the fine coordinate to the ruff coordinate
-float psImageBinningGetRuffX (psImageBinning *binning, float xFine) {
-    float xRuff = (xFine - binning->nXskip)/binning->nXbin;
+double psImageBinningGetRuffX (const psImageBinning *binning, const double xFine) {
+
+    PS_ASSERT_INT_POSITIVE(binning->nXbin, NAN);
+    PS_ASSERT_INT_POSITIVE(binning->nYbin, NAN);
+    PS_ASSERT_INT_LESS_THAN_OR_EQUAL(binning->nXskip, binning->nXbin, NAN);
+    PS_ASSERT_INT_LESS_THAN_OR_EQUAL(binning->nYskip, binning->nYbin, NAN);
+
+    double xRuff = (xFine - binning->nXskip)/binning->nXbin;
     return xRuff;
 }
-float psImageBinningGetRuffY (psImageBinning *binning, float yFine) {
-    float yRuff = (yFine - binning->nYskip)/binning->nYbin;
+double psImageBinningGetRuffY (const psImageBinning *binning, const double yFine) {
+
+    PS_ASSERT_INT_POSITIVE(binning->nXbin, NAN);
+    PS_ASSERT_INT_POSITIVE(binning->nYbin, NAN);
+    PS_ASSERT_INT_LESS_THAN_OR_EQUAL(binning->nXskip, binning->nXbin, NAN);
+    PS_ASSERT_INT_LESS_THAN_OR_EQUAL(binning->nYskip, binning->nYbin, NAN);
+
+    double yRuff = (yFine - binning->nYskip)/binning->nYbin;
     return yRuff;
 }
 
 // convert the ruff coordinate to the fine coordinate
-float psImageBinningGetFineX (psImageBinning *binning, float xRuff) {
-    float xFine = xRuff * binning->nXbin + binning->nXskip;
+double psImageBinningGetFineX (const psImageBinning *binning, const double xRuff) {
+
+    PS_ASSERT_INT_POSITIVE(binning->nXbin, NAN);
+    PS_ASSERT_INT_POSITIVE(binning->nYbin, NAN);
+    PS_ASSERT_INT_LESS_THAN_OR_EQUAL(binning->nXskip, binning->nXbin, NAN);
+    PS_ASSERT_INT_LESS_THAN_OR_EQUAL(binning->nYskip, binning->nYbin, NAN);
+
+    double xFine = xRuff * binning->nXbin + binning->nXskip;
     return xFine;
 }
-float psImageBinningGetFineY (psImageBinning *binning, float yRuff) {
-    float yFine = yRuff * binning->nYbin + binning->nYskip;
+double psImageBinningGetFineY (const psImageBinning *binning, const double yRuff) {
+
+    PS_ASSERT_INT_POSITIVE(binning->nXbin, NAN);
+    PS_ASSERT_INT_POSITIVE(binning->nYbin, NAN);
+    PS_ASSERT_INT_LESS_THAN_OR_EQUAL(binning->nXskip, binning->nXbin, NAN);
+    PS_ASSERT_INT_LESS_THAN_OR_EQUAL(binning->nYskip, binning->nYbin, NAN);
+
+    double yFine = yRuff * binning->nYbin + binning->nYskip;
     return yFine;
 }
Index: /branches/eam_branch_20070830/psLib/src/imageops/psImageBinning.h
===================================================================
--- /branches/eam_branch_20070830/psLib/src/imageops/psImageBinning.h	(revision 14858)
+++ /branches/eam_branch_20070830/psLib/src/imageops/psImageBinning.h	(revision 14859)
@@ -8,6 +8,6 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.2.2.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-07 20:19:05 $
+ *  @version $Revision: 1.2.2.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-15 19:43:43 $
  *
  *  Copyright 2007 Institute for Astronomy, University of Hawaii
@@ -51,8 +51,8 @@
 psRegion psImageBinningSetRuffRegion (psImageBinning *binning, psRegion fineRegion);
 
-float psImageBinningGetRuffX (psImageBinning *binning, float xFine);
-float psImageBinningGetRuffY (psImageBinning *binning, float yFine);
-float psImageBinningGetFineX (psImageBinning *binning, float xRuff);
-float psImageBinningGetFineY (psImageBinning *binning, float yRuff);
+double psImageBinningGetRuffX (const psImageBinning *binning, const double xFine);
+double psImageBinningGetRuffY (const psImageBinning *binning, const double yFine);
+double psImageBinningGetFineX (const psImageBinning *binning, const double xRuff);
+double psImageBinningGetFineY (const psImageBinning *binning, const double yRuff);
 
 /// @}
Index: /branches/eam_branch_20070830/psLib/src/imageops/psImageMap.h
===================================================================
--- /branches/eam_branch_20070830/psLib/src/imageops/psImageMap.h	(revision 14858)
+++ /branches/eam_branch_20070830/psLib/src/imageops/psImageMap.h	(revision 14859)
@@ -7,6 +7,6 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.1.2.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-07 20:19:35 $
+ *  @version $Revision: 1.1.2.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-15 19:43:43 $
  *
  *  Copyright 2007 Institute for Astronomy, University of Hawaii
@@ -48,4 +48,7 @@
 psVector *psImageMapEvalVector (psImageMap *map, psVector *x, psVector *y);
 
+// fit the image map to a set of points
+bool psImageMapFit (psImageMap *map, psVector *x, psVector *y, psVector *f, psVector *df);
+
 /// @}
 #endif // #ifndef PS_IMAGE_MAP_H
Index: /branches/eam_branch_20070830/psLib/src/imageops/psImageMapFit.c
===================================================================
--- /branches/eam_branch_20070830/psLib/src/imageops/psImageMapFit.c	(revision 14859)
+++ /branches/eam_branch_20070830/psLib/src/imageops/psImageMapFit.c	(revision 14859)
@@ -0,0 +1,356 @@
+/** @file  psImageMapFit.c
+ *
+ *  @brief Functions define a 2d coarse representation of a finer 2D field
+ *
+ *  @ingroup Image
+ *
+ *  @author Eugene Magnier, IfA
+ *
+ *  @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-15 19:43:43 $
+ *
+ *  Copyright 2007 Institute for Astronomy, University of Hawaii
+ */
+
+#ifdef HAVE_CONFIG_H
+# include "config.h"
+#endif
+
+#include <stdio.h>
+#include "psError.h"
+#include "psAbort.h"
+
+// XXX for testing
+#include "psFits.h"
+#include "psFitsImage.h"
+
+#include "psAssert.h"
+#include "psRegion.h"
+
+#include "psVector.h"
+#include "psImage.h"
+#include "psMatrix.h"
+#include "psStats.h"
+#include "psImageBinning.h"
+#include "psImageMap.h"
+// #include "psImagePixelInterpolate.h"
+// #include "psImageUnbin.h"
+
+// given a randomly-sampled field of values & weights at points: (f, df) @ (x, y), find the
+// best fit image from which Bilinear interpolation yields the input field.  The fitted image
+// consists of a grid of values g(n,m) at coordinates (n,m).
+
+// relationship between x,y and n,m coordinates:
+
+// map defines the output image dimensions and scaling.
+bool psImageMapFit (psImageMap *map, psVector *x, psVector *y, psVector *f, psVector *df) {
+
+    int I, J;
+
+    // dimensions of the output map image
+    int Nx = map->binning->nXruff;
+    int Ny = map->binning->nYruff;
+
+    // set up the redirection table so we can use sA[-1][-1], etc
+    float SAm[3][3], *SAv[3], **sA;
+    float TAm[3][3], *TAv[3], **tA;
+
+    for (int i = 0; i < 3; i++) {
+	SAv[i] = SAm[i] + 1;
+	TAv[i] = TAm[i] + 1;
+    }
+    sA = SAv + 1;
+    tA = TAv + 1;
+
+    // elements of the matrix equation Ax = B; we are solving for the vector x 
+    psImage *A = psImageAlloc (Nx*Ny, Nx*Ny, PS_TYPE_F32);
+    psVector *B = psVectorAlloc (Nx*Ny, PS_TYPE_F32);
+
+    // we are looping over the Nx,Ny image map elements; 
+    // the matrix equation contains Nx*Ny rows and columns 
+    // for (int n = 1; n < Nx - 1; n++) {
+    // for (int m = 1; m < Ny - 1; m++) {
+    
+    float Total = 0.0;
+    for (int n = 0; n < Nx; n++) {
+	for (int m = 0; m < Ny; m++) {
+	    // define & init summing variables
+	    float rx_rx_ry_ry = 0;
+	    float rx_rx_dy_ry = 0;
+	    float dx_rx_ry_ry = 0;
+	    float dx_rx_dy_ry = 0;
+	    float fi_rx_ry    = 0;
+	    float rx_rx_py_py = 0;
+	    float rx_rx_qy_py = 0;
+	    float dx_rx_py_py = 0;
+	    float dx_rx_qy_py = 0;
+	    float fi_rx_py    = 0;
+	    float px_px_ry_ry = 0;
+	    float px_px_dy_ry = 0;
+	    float qx_px_ry_ry = 0;
+	    float qx_px_dy_ry = 0;
+	    float fi_px_ry    = 0;
+	    float px_px_py_py = 0;
+	    float px_px_qy_py = 0;
+	    float qx_px_py_py = 0;
+	    float qx_px_qy_py = 0;
+	    float fi_px_py    = 0;
+
+	    // generate the sums for the fitting matrix element I,J
+	    // I = n + nX*m
+	    // J = (n + jn) + nX*(m + jm)
+	    for (int i = 0; i < x->n; i++) {
+		// base coordinate offset for this point (x,y) relative to this map element (n,m)
+		// float dx = x->data.F32[i] - psImageBinningGetFineX (map->binning, n + 0.5);
+		// float dy = y->data.F32[i] - psImageBinningGetFineY (map->binning, m + 0.5);
+
+		float dx = psImageBinningGetRuffX (map->binning, x->data.F32[i]) - (n + 0.5);
+		float dy = psImageBinningGetRuffY (map->binning, y->data.F32[i]) - (m + 0.5);
+
+		// skip points outside of 2x2 grid centered on n,m:
+		if (fabs(dx) > 1.0) continue;
+		if (fabs(dy) > 1.0) continue;
+
+		// related offset values
+		float rx = 1.0 - dx;
+		float ry = 1.0 - dy;
+		float px = 1.0 + dx;
+		float py = 1.0 + dy;
+		float qx = -dx;
+		float qy = -dy;
+
+		// data value & weight for this point
+		float fi = f->data.F32[i];
+		float wt = 1.0;
+		if (df != NULL) {
+		    if (df->data.F32[i] == 0.0) {
+			wt = 0.0;
+		    } else {
+			wt = 1.0 / PS_SQR(df->data.F32[i]); // XXX test for dz == NULL or dz_i = 0
+		    }
+		}
+
+		// sum the appropriate elements for the different quadrants
+
+		// points at offset 1,1
+		if ((dx >= 0) && (dy >= 0)) {
+		    rx_rx_ry_ry += rx*rx*ry*ry*wt;
+		    rx_rx_dy_ry += rx*rx*dy*ry*wt;
+		    dx_rx_ry_ry += dx*rx*ry*ry*wt;
+		    dx_rx_dy_ry += dx*rx*dy*ry*wt;
+		    fi_rx_ry    += fi*rx*ry*wt;
+		}
+		// points at offset 1,0
+		if ((dx >= 0) && (dy <  0)) {
+		    rx_rx_py_py += rx*rx*py*py*wt;
+		    rx_rx_qy_py += rx*rx*qy*py*wt;
+		    dx_rx_py_py += dx*rx*py*py*wt;
+		    dx_rx_qy_py += dx*rx*qy*py*wt;
+		    fi_rx_py    += fi*rx*py*wt;
+		}
+		// points at offset 0,1
+		if ((dx <  0) && (dy >= 0)) {
+		    px_px_ry_ry += px*px*ry*ry*wt;
+		    px_px_dy_ry += px*px*dy*ry*wt;
+		    qx_px_ry_ry += qx*px*ry*ry*wt;
+		    qx_px_dy_ry += qx*px*dy*ry*wt;
+		    fi_px_ry    += fi*px*ry*wt;
+		}
+		// points at offset 0,0
+		if ((dx <  0) && (dy <  0)) {
+		    px_px_py_py += px*px*py*py*wt;
+		    px_px_qy_py += px*px*qy*py*wt;
+		    qx_px_py_py += qx*px*py*py*wt;
+		    qx_px_qy_py += qx*px*qy*py*wt;
+		    fi_px_py    += fi*px*py*wt;
+		}
+	    }		
+
+	    // the chi-square derivatives have elements of the form g(n+jn,m+jm)*A(jn,jm),
+	    // jn,jm = -1 to +1. Convert the sums above into the correct coefficients
+	    sA[-1][-1] = qx_px_qy_py;
+	    sA[-1][ 0] = qx_px_ry_ry + qx_px_py_py;
+	    sA[-1][+1] = qx_px_dy_ry;
+	    sA[ 0][-1] = rx_rx_qy_py + px_px_qy_py;
+	    sA[ 0][ 0] = rx_rx_ry_ry + px_px_ry_ry + rx_rx_py_py + px_px_py_py;
+	    sA[ 0][+1] = rx_rx_dy_ry + px_px_dy_ry;
+	    sA[+1][-1] = dx_rx_qy_py;
+	    sA[+1][ 0] = dx_rx_ry_ry + dx_rx_py_py;
+	    sA[+1][+1] = dx_rx_dy_ry;
+	    
+	    // in for edge cases, we only have certain grid points available.  Depending on location
+	    // in the plane, we need to re-arrange the terms to account for this.  we
+	    // implicitly extrapolate across the g[n][m] grid point, assuming, eg g[n+1][m] =
+	    // 2*g[n][m] - g[n-1][m].  generate the true output coefficients for the different
+	    // cases.  start by zeroing out this array.
+	    for (int jn = -1; jn <= +1; jn++) {
+		for (int jm = -1; jm <= +1; jm++) {
+		    tA[jn][jm] = 0;
+		}
+	    }
+
+	    // XXX test:
+	    // goto skip_ahead;
+
+	    // XXX need to add three additional edge cases:
+	    if ((Nx == 1) && (Ny == 1)) {
+		psAbort ("un-implemented edge case");
+		goto insert;
+	    }
+	    if (Nx == 1) {
+		psAbort ("un-implemented edge case");
+		goto insert;
+	    }
+	    if (Ny == 1) {
+		psAbort ("un-implemented edge case");
+		goto insert;
+	    }
+
+	    // UPPER-RIGHT Corner
+	    if ((n == Nx - 1) && (m == Ny - 1)) {
+		tA[-1][-1] = sA[-1][-1] -   sA[+1][+1] -   sA[+1][-1] -   sA[-1][+1];
+		tA[ 0][-1] = sA[ 0][-1] + 2*sA[+1][-1] -   sA[ 0][+1];
+		tA[-1][ 0] = sA[-1][ 0] + 2*sA[-1][+1] -   sA[+1][ 0];
+		tA[ 0][ 0] = sA[ 0][ 0] + 2*sA[+1][ 0] + 2*sA[ 0][+1] + 2*sA[+1][+1];
+		goto insert;
+	    }
+	    // LOWER-RIGHT Corner 
+	    if ((n == Nx - 1) && (m == 0)) {
+		tA[-1][+1] = sA[-1][+1] -   sA[+1][-1] -   sA[+1][+1] -   sA[-1][-1];
+		tA[ 0][+1] = sA[ 0][+1] + 2*sA[+1][+1] -   sA[ 0][-1];
+		tA[-1][ 0] = sA[-1][ 0] + 2*sA[-1][-1] -   sA[+1][ 0];
+		tA[ 0][ 0] = sA[ 0][ 0] + 2*sA[+1][ 0] + 2*sA[ 0][-1] + 2*sA[+1][-1];
+		goto insert;
+	    }
+	    // LOWER-LEFT Corner 
+	    if ((n == 0) && (m == 0)) {
+		tA[+1][+1] = sA[+1][+1] -   sA[-1][-1] -   sA[-1][+1] -   sA[+1][-1];
+		tA[ 0][+1] = sA[ 0][+1] + 2*sA[-1][+1] -   sA[ 0][-1];
+		tA[+1][ 0] = sA[+1][ 0] + 2*sA[+1][-1] -   sA[-1][ 0];
+		tA[ 0][ 0] = sA[ 0][ 0] + 2*sA[-1][ 0] + 2*sA[ 0][-1] + 2*sA[-1][-1];
+		goto insert;
+	    }
+	    // UPPER-LEFT Corner 
+	    if ((n == 0) && (m == Ny - 1)) {
+		tA[+1][-1] = sA[+1][-1] -   sA[-1][+1] -   sA[-1][-1] -   sA[+1][+1];
+		tA[ 0][-1] = sA[ 0][-1] + 2*sA[-1][-1] -   sA[ 0][+1];
+		tA[+1][ 0] = sA[+1][ 0] + 2*sA[+1][+1] -   sA[-1][ 0];
+		tA[ 0][ 0] = sA[ 0][ 0] + 2*sA[-1][ 0] + 2*sA[ 0][+1] + 2*sA[-1][+1];
+		goto insert;
+	    }
+	    // LEFT Edge 
+	    if (n == 0) {
+		tA[+1][-1] = sA[+1][-1] - sA[-1][-1];
+		tA[+1][ 0] = sA[+1][ 0] - sA[-1][ 0];
+		tA[+1][+1] = sA[+1][+1] - sA[-1][+1];
+		tA[ 0][-1] = sA[ 0][-1] + 2.0*sA[-1][-1];
+		tA[ 0][ 0] = sA[ 0][ 0] + 2.0*sA[-1][ 0];
+		tA[ 0][+1] = sA[ 0][+1] + 2.0*sA[-1][+1];
+		goto insert;
+	    }
+	    // RIGHT Edge 
+	    if (n == Nx - 1) {
+		tA[-1][-1] = sA[-1][-1] - sA[+1][-1];
+		tA[-1][ 0] = sA[-1][ 0] - sA[+1][ 0];
+		tA[-1][+1] = sA[-1][+1] - sA[+1][+1];
+		tA[ 0][-1] = sA[ 0][-1] + 2.0*sA[+1][-1];
+		tA[ 0][ 0] = sA[ 0][ 0] + 2.0*sA[+1][ 0];
+		tA[ 0][+1] = sA[ 0][+1] + 2.0*sA[+1][+1];
+		goto insert;
+	    }
+	    // BOTTOM Edge 
+	    if (m == 0) {
+		tA[-1][+1] = sA[-1][+1] - sA[-1][-1];
+		tA[ 0][+1] = sA[ 0][+1] - sA[ 0][-1];
+		tA[+1][+1] = sA[+1][+1] - sA[+1][-1];
+		tA[-1][ 0] = sA[-1][ 0] + 2.0*sA[-1][-1];
+		tA[ 0][ 0] = sA[ 0][ 0] + 2.0*sA[ 0][-1];
+		tA[+1][ 0] = sA[+1][ 0] + 2.0*sA[+1][-1];
+		goto insert;
+	    }
+	    // TOP Edge 
+	    if (m == Ny - 1) {
+		tA[-1][-1] = sA[-1][-1] - sA[-1][+1];
+		tA[ 0][-1] = sA[ 0][-1] - sA[ 0][+1];
+		tA[+1][-1] = sA[+1][-1] - sA[+1][+1];
+		tA[-1][ 0] = sA[-1][ 0] + 2.0*sA[-1][+1];
+		tA[ 0][ 0] = sA[ 0][ 0] + 2.0*sA[ 0][+1];
+		tA[+1][ 0] = sA[+1][ 0] + 2.0*sA[+1][+1];
+		goto insert;
+	    }
+
+	    // skip_ahead:
+	    // center pixels:
+	    for (int jn = -1; jn <= +1; jn++) {
+		for (int jm = -1; jm <= +1; jm++) {
+		    tA[jn][jm] = sA[jn][jm];
+		}
+	    }
+
+	insert:
+	    // I[ 0][ 0] = index for this n,m element:
+	    I = n + Nx * m;
+	    B->data.F32[I] = fi_rx_ry + fi_rx_py + fi_px_ry + fi_px_py;
+
+	    // insert these values into their corresponding locations in A, B
+	    float Sum = 0.0;
+	    for (int jn = -1; jn <= +1; jn++) {
+		if (n + jn <   0) continue;
+		if (n + jn >= Nx) continue;
+		for (int jm = -1; jm <= +1; jm++) {
+		    if (m + jm <   0) continue;
+		    if (m + jm >= Ny) continue;
+		    J = (n + jn) + Nx * (m + jm);
+		    A->data.F32[J][I] = tA[jn][jm];
+		    fprintf (stderr, "A %d %d (%d %d : %d %d): %f\n", I, J, n, m, n + jn, m + jm, tA[jn][jm]);
+		    Sum += tA[jn][jm];
+		}
+	    }
+	    fprintf (stderr, "B %d (%d %d) : %f  :  %f\n", I, n, m, B->data.F32[I], Sum);
+	    Total += Sum;
+	}
+    }
+
+    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);
+
+    if (!psMatrixGJSolveF32(A, B)) {
+	psAbort ("failed on linear equations");
+	psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
+	psFree (A);
+	psFree (B);
+	return false;
+    }
+    
+    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];
+	}
+    }
+
+    psFree (A);
+    psFree (B);
+
+    return true;
+}
Index: /branches/eam_branch_20070830/psLib/src/imageops/psImageUnbin.c
===================================================================
--- /branches/eam_branch_20070830/psLib/src/imageops/psImageUnbin.c	(revision 14858)
+++ /branches/eam_branch_20070830/psLib/src/imageops/psImageUnbin.c	(revision 14859)
@@ -7,6 +7,6 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-03-27 02:43:22 $
+ *  @version $Revision: 1.7.8.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-15 19:43:43 $
  *
  *  Copyright 2007 Institute for Astronomy, University of Hawaii
@@ -205,4 +205,7 @@
  * N.b. This code only works for the central part of the image; the edge
  * cases should be added
+
+ * XXXX this is bilinear interpolation, but written sub-optimally
+ * XXXX this function should be taking float input coordinates!!!
  */
 double psImageUnbinPixel(const int ix, const int iy, // desired Unbinned point (parent coords)
@@ -249,2 +252,65 @@
     return Vxs + (Vxe - Vxs)*(ix - xs)/DX; // value at [iy][ix]
 }
+
+double psImageUnbinPixel_V2(const double xFine, const double yFine, // desired Unbinned point (parent coords)
+			    const psImage *in, // binned image
+			    const psImageBinning *binning)   //!< Overhang
+{
+    PS_ASSERT_IMAGE_NON_NULL(in, NAN);
+    assert (in->type.type == PS_TYPE_F32);
+
+    const float xRuff = psImageBinningGetRuffX (binning, xFine);
+    const float yRuff = psImageBinningGetRuffY (binning, yFine);
+
+    const double value = psImageInterpolatePixelBilinear (xRuff, yRuff, in);
+
+    return value;
+}
+
+// fast & simple API to interpolate to a subpixel position using bilinear interpolation
+// x,y in parent image coordinates (pixel centers at 0.5, 0.5)
+double psImageInterpolatePixelBilinear (const double xIn, const double yIn, const psImage *in) {
+
+    PS_ASSERT_PTR_NON_NULL(in, PS_ERR_BAD_PARAMETER_VALUE);
+    assert (in->type.type == PS_TYPE_F32);
+
+    const double x = xIn - in->col0;
+    const double y = yIn - in->row0;
+
+    // allow extrapolation to edge of valid pixels, but not beyond
+    if ((x < 0) || (x >= in->numCols) || (y < 0) || (y >= in->numRows)) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Point (%f,%f) lies outside binned image", x, y);
+        return NAN;
+    }
+
+    // handle edge cases with extrapolation
+
+    const int ix = x - 0.5; // index of reference pixel
+    const int iy = y - 0.5; // index of reference pixel
+
+    const int Xs = PS_MAX (PS_MIN (ix, in->numCols - 2), 0);
+    const int Ys = PS_MAX (PS_MIN (iy, in->numRows - 2), 0);
+
+    const int Xe = Xs + 1;
+    const int Ye = Ys + 1;
+
+    // dx,dy range from 0.0 to 1.0 for interpolated pixels, and -0.5 to 1.5 for extrapolation
+    const double dx = x - 0.5 - Xs;
+    const double dy = y - 0.5 - Ys;
+
+    const double rx = 1.0 - dx;
+    const double ry = 1.0 - dy;
+
+    // Vxy 
+    double V00 = in->data.F32[Ys][Xs];
+    double V10 = in->data.F32[Ys][Xe];
+    double V01 = in->data.F32[Ye][Xs];
+    double V11 = in->data.F32[Ye][Xe];
+
+    // bilinear interpolation
+    const double value = V00*rx*ry + V10*rx*dy + V01*dx*ry + V11*dx*dy;
+
+    return value;
+}
+
+    
Index: /branches/eam_branch_20070830/psLib/src/imageops/psImageUnbin.h
===================================================================
--- /branches/eam_branch_20070830/psLib/src/imageops/psImageUnbin.h	(revision 14858)
+++ /branches/eam_branch_20070830/psLib/src/imageops/psImageUnbin.h	(revision 14859)
@@ -6,6 +6,6 @@
  *  @author Robert DeSonia, MHPCC
  *
- *  $Revision: 1.4 $ $Name: not supported by cvs2svn $
- *  $Date: 2007-03-27 02:43:22 $
+ *  $Revision: 1.4.8.1 $ $Name: not supported by cvs2svn $
+ *  $Date: 2007-09-15 19:43:43 $
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
  */
@@ -28,4 +28,11 @@
                         );
 
+double psImageUnbinPixel_V2(const double xFine, const double yFine, // desired Unbinned point (parent coords)
+			    const psImage *in, // binned image
+			    const psImageBinning *binning   //!< Overhang
+    );
+
+double psImageInterpolatePixelBilinear (const double xIn, const double yIn, const psImage *in);
+
 /// @}
 #endif
Index: /branches/eam_branch_20070830/psLib/test/imageops/Makefile.am
===================================================================
--- /branches/eam_branch_20070830/psLib/test/imageops/Makefile.am	(revision 14858)
+++ /branches/eam_branch_20070830/psLib/test/imageops/Makefile.am	(revision 14859)
@@ -23,4 +23,5 @@
 	tap_psImageInterpolate2 \
 	tap_psImageMap \
+	tap_psImageMapFit \
 	tap_psImageMaskOps
 
Index: /branches/eam_branch_20070830/psLib/test/imageops/tap_psImageMapFit.c
===================================================================
--- /branches/eam_branch_20070830/psLib/test/imageops/tap_psImageMapFit.c	(revision 14859)
+++ /branches/eam_branch_20070830/psLib/test/imageops/tap_psImageMapFit.c	(revision 14859)
@@ -0,0 +1,531 @@
+#include <stdio.h>
+#include <string.h>
+#include <pslib.h>
+
+#include "tap.h"
+#include "pstap.h"
+
+int SaveImage (psMetadata *header, psImage *image, char *filename) {
+
+    psFits *fits = psFitsOpen (filename, "w");
+    psFitsWriteImage (fits, NULL, image, 0, NULL);
+    psFitsClose (fits);
+    return (TRUE);
+}
+
+# define TEST_4PT_0 0
+# define TEST_4PT_1 0
+# define TEST_4PT_2 0
+# define TEST_4PT_3 0
+
+# define TEST_9PT_0 0
+# define TEST_9PT_1 0
+# define TEST_9PT_2 0
+
+# define TEST_6PT_0 1
+
+# define TEST1 0
+# define TEST2 0
+# define TEST3 0
+
+# define C00 +0.0
+# define C01 +0.0
+# define C10 +1.0
+
+int main (void)
+{
+
+    plan_tests(1);
+    
+    // *** tests to demonstrate the validity of the algorithm or concept ***
+
+    // test for more points: 3x3 grid fitted to 9 points with simple slope and scale difference
+    # if (TEST_6PT_0)
+    {
+	psImageBinning *binning = psImageBinningAlloc();
+	binning->nXfine = 3;
+	binning->nYfine = 2;
+	binning->nXruff = 3;
+	binning->nYruff = 2;
+
+	// 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 = 0.5; ix < 3.0; ix += 1.0) {
+	    for (int iy = 0.5; iy < 2.0; iy += 1.0) {
+		x->data.F32[n] = ix;
+		y->data.F32[n] = iy;
+		# if (0)
+		if ((ix == 3.0) && (iy == 1.0)) {
+		    x->data.F32[n] = ix + 0.0;
+		    y->data.F32[n] = iy + 0.1; // add in both points.  
+		    f->data.F32[n] = ix + iy;
+		    n++;
+		    x->data.F32[n] = ix + 0.0;
+		    y->data.F32[n] = iy - 0.1;
+		} else {
+		    x->data.F32[n] = ix;
+		    y->data.F32[n] = iy;
+		}
+		# endif
+		f->data.F32[n] = ix + iy;
+		n++;
+	    }
+	}
+	x->n = n;
+	y->n = n;
+	f->n = n;
+
+	psImage *field = psImageAlloc(3, 2, PS_TYPE_F32);
+	for (int ix = 0; ix < 3; ix++) {
+	    for (int iy = 0; iy < 2; 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);
+	psFree (binning);
+
+	SaveImage (NULL, map->map, "map.fits");
+	
+	SaveImage (NULL, field, "field.fits");
+
+	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);
+		is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], FLT_EPSILON, "model matches inputs");
+	    }
+	}
+	SaveImage (NULL, model, "model.fits");
+        // 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
+    # if (TEST_9PT_2)
+    {
+	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; ix += 2.0) {
+	    for (int iy = 1.0; iy < 6; iy += 2.0) {
+		if ((ix == 3.0) && (iy == 1.0)) {
+		    x->data.F32[n] = ix + 0.0;
+		    y->data.F32[n] = iy + 0.1; // add in both points.  
+		    f->data.F32[n] = ix + iy;
+		    n++;
+		    x->data.F32[n] = ix + 0.0;
+		    y->data.F32[n] = iy - 0.1;
+		} else {
+		    x->data.F32[n] = ix;
+		    y->data.F32[n] = iy;
+		}
+		f->data.F32[n] = ix + iy;
+		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);
+	psFree (binning);
+
+	SaveImage (NULL, map->map, "map.fits");
+	
+	SaveImage (NULL, field, "field.fits");
+
+	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);
+		is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], FLT_EPSILON, "model matches inputs");
+	    }
+	}
+	SaveImage (NULL, model, "model.fits");
+        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    # endif
+
+    // still a simple test: 3x3 grid fitted to 9 points with simple slope and scale difference
+    # if (TEST_9PT_1)
+    {
+	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 (9, PS_TYPE_F32);
+	psVector *y = psVectorAlloc (9, PS_TYPE_F32);
+	psVector *f = psVectorAlloc (9, 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 (int ix = 1; ix < 6; ix += 2) {
+	    for (int iy = 1; iy < 6; iy += 2) {
+		x->data.F32[n] = ix;
+		y->data.F32[n] = iy;
+		f->data.F32[n] = ix + iy;
+		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);
+	psFree (binning);
+
+	SaveImage (NULL, map->map, "map.fits");
+	
+	SaveImage (NULL, field, "field.fits");
+
+	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);
+		is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], FLT_EPSILON, "model matches inputs");
+	    }
+	}
+	SaveImage (NULL, model, "model.fits");
+        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    # endif
+
+    // very simple test: 3x3 grid fitted to 9 points with simple slope
+    # if (TEST_9PT_0)
+    {
+	// function is defined over the range 0-1000, 0-1000
+	psImageBinning *binning = psImageBinningAlloc();
+	binning->nXfine = 3;
+	binning->nYfine = 3;
+	binning->nXruff = 3;
+	binning->nYruff = 3;
+
+	// generate a grid of test data points
+	psVector *x = psVectorAlloc (9, PS_TYPE_F32);
+	psVector *y = psVectorAlloc (9, PS_TYPE_F32);
+	psVector *f = psVectorAlloc (9, PS_TYPE_F32);
+
+	int n = 0;
+	psImage *field = psImageAlloc(3, 3, PS_TYPE_F32);
+	for (int ix = 0; ix < 3; ix++) {
+	    for (int iy = 0; iy < 3; iy++) {
+		x->data.F32[n] = ix + 0.5;
+		y->data.F32[n] = iy + 0.5;
+		f->data.F32[n] = ix + iy;
+		field->data.F32[iy][ix] = f->data.F32[n];
+		n++;
+	    }
+	}
+
+	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);
+	psFree (binning);
+
+	SaveImage (NULL, map->map, "map.fits");
+	
+	SaveImage (NULL, field, "field.fits");
+
+	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);
+		is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], FLT_EPSILON, "model matches inputs");
+	    }
+	}
+	SaveImage (NULL, model, "model.fits");
+        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    # endif
+
+    // more complex test: 2x2 grid fitted to 16 points with simple slope offset from grid centers
+    // this one uses points inset relative to the reference points so they are all visible
+    // to the reference
+    # if (TEST_4PT_3)
+    {
+	// function is defined over the range 0-1000, 0-1000
+	psImageBinning *binning = psImageBinningAlloc();
+	binning->nXfine = 4;
+	binning->nYfine = 4;
+	binning->nXruff = 2;
+	binning->nYruff = 2;
+
+	// generate a grid of test data points
+	psVector *x = psVectorAlloc (16, PS_TYPE_F32);
+	psVector *y = psVectorAlloc (16, PS_TYPE_F32);
+	psVector *f = psVectorAlloc (16, PS_TYPE_F32);
+
+	int n = 0;
+	// actual field is f = x + y, where x & y are the subpixel positions
+	for (float ix = 0.5; ix < 4; ix += 1.0) {
+	    for (float iy = 0.5; iy < 4; iy += 1.0) {
+		x->data.F32[n] = ix;
+		y->data.F32[n] = iy;
+		f->data.F32[n] = x->data.F32[n] + y->data.F32[n];
+		n++;
+	    }
+	}
+
+	psImage *field = psImageAlloc(4, 4, PS_TYPE_F32);
+	for (int ix = 0; ix < 4; ix++) {
+	    for (int iy = 0; iy < 4; 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);
+	psFree (binning);
+
+	SaveImage (NULL, map->map, "map.fits");
+	
+	SaveImage (NULL, field, "field.fits");
+
+	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);
+		is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], 10*FLT_EPSILON, "model matches inputs");
+	    }
+	}
+	SaveImage (NULL, model, "model.fits");
+        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    # endif
+
+    // more complex test: 2x2 grid fitted to 4 points with simple slope offset from grid centers
+    // this one uses points inset relative to the reference points so they are all visible
+    // to the reference
+    # if (TEST_4PT_2)
+    {
+	// function is defined over the range 0-1000, 0-1000
+	psImageBinning *binning = psImageBinningAlloc();
+	binning->nXfine = 2;
+	binning->nYfine = 2;
+	binning->nXruff = 2;
+	binning->nYruff = 2;
+
+	// generate a grid of test data points
+	psVector *x = psVectorAlloc (4, PS_TYPE_F32);
+	psVector *y = psVectorAlloc (4, PS_TYPE_F32);
+	psVector *f = psVectorAlloc (4, PS_TYPE_F32);
+
+	int n = 0;
+	psImage *field = psImageAlloc(2, 2, PS_TYPE_F32);
+	// actual field is f = x + y, where x & y are the subpixel positions
+	for (int ix = 0; ix < 2; ix++) {
+	    for (int iy = 0; iy < 2; iy++) {
+		// add in one extra point on the surface
+		if (ix) {
+		    x->data.F32[n] = ix + 0.5 - 0.1;
+		} else {
+		    x->data.F32[n] = ix + 0.5 + 0.1;
+		}
+		if (iy) {
+		    y->data.F32[n] = iy + 0.5 - 0.1;
+		} else {
+		    y->data.F32[n] = iy + 0.5 + 0.1;
+		}
+		f->data.F32[n] = x->data.F32[n] + y->data.F32[n];
+		field->data.F32[iy][ix] = ix + 0.5 + iy + 0.5;
+		n++;
+	    }
+	}
+
+	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);
+	psFree (binning);
+
+	SaveImage (NULL, map->map, "map.fits");
+	
+	SaveImage (NULL, field, "field.fits");
+
+	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);
+		is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], 5*FLT_EPSILON, "model matches inputs");
+	    }
+	}
+	SaveImage (NULL, model, "model.fits");
+        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    # endif
+
+    // still a simple test: 2x2 grid fitted to 4 points with simple slope and scale difference
+    # if (TEST_4PT_1)
+    {
+	// function is defined over the range 0-1000, 0-1000
+	psImageBinning *binning = psImageBinningAlloc();
+	binning->nXfine = 4;
+	binning->nYfine = 4;
+	binning->nXruff = 2;
+	binning->nYruff = 2;
+
+	// generate a grid of test data points
+	psVector *x = psVectorAlloc (4, PS_TYPE_F32);
+	psVector *y = psVectorAlloc (4, PS_TYPE_F32);
+	psVector *f = psVectorAlloc (4, PS_TYPE_F32);
+
+	// the underlying field is f = ix + iy, where ix,iy are fine pixel coordinates
+	int n = 0;
+	for (int ix = 1; ix < 4; ix += 2) {
+	    for (int iy = 1; iy < 4; iy += 2) {
+		x->data.F32[n] = ix;
+		y->data.F32[n] = iy;
+		f->data.F32[n] = ix + iy;
+		n++;
+	    }
+	}
+
+	psImage *field = psImageAlloc(4, 4, PS_TYPE_F32);
+	for (int ix = 0; ix < 4; ix++) {
+	    for (int iy = 0; iy < 4; 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);
+	psFree (binning);
+
+	SaveImage (NULL, map->map, "map.fits");
+	
+	SaveImage (NULL, field, "field.fits");
+
+	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);
+		is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], FLT_EPSILON, "model matches inputs");
+	    }
+	}
+	SaveImage (NULL, model, "model.fits");
+        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    # endif
+
+    // simplest possible test: 2x2 grid fitted to 4 points with simple slope
+    # if (TEST_4PT_0)
+    {
+	// function is defined over the range 0-1000, 0-1000
+	psImageBinning *binning = psImageBinningAlloc();
+	binning->nXfine = 2;
+	binning->nYfine = 2;
+	binning->nXruff = 2;
+	binning->nYruff = 2;
+
+	// generate a grid of test data points
+	psVector *x = psVectorAlloc (4, PS_TYPE_F32);
+	psVector *y = psVectorAlloc (4, PS_TYPE_F32);
+	psVector *f = psVectorAlloc (4, PS_TYPE_F32);
+
+	int n = 0;
+	psImage *field = psImageAlloc(2, 2, PS_TYPE_F32);
+	for (int ix = 0; ix < 2; ix++) {
+	    for (int iy = 0; iy < 2; iy++) {
+		x->data.F32[n] = ix + 0.5;
+		y->data.F32[n] = iy + 0.5;
+		f->data.F32[n] = x->data.F32[n] + y->data.F32[n];
+		field->data.F32[iy][ix] = f->data.F32[n];
+		n++;
+	    }
+	}
+
+	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);
+	psFree (binning);
+
+	SaveImage (NULL, map->map, "map.fits");
+	
+	SaveImage (NULL, field, "field.fits");
+
+	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);
+		is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], FLT_EPSILON, "model matches inputs");
+	    }
+	}
+	SaveImage (NULL, model, "model.fits");
+        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    # endif
+
+    return exit_status();
+}   
Index: /branches/eam_branch_20070830/psLib/test/imageops/tap_psImageMapFit2.c
===================================================================
--- /branches/eam_branch_20070830/psLib/test/imageops/tap_psImageMapFit2.c	(revision 14859)
+++ /branches/eam_branch_20070830/psLib/test/imageops/tap_psImageMapFit2.c	(revision 14859)
@@ -0,0 +1,216 @@
+
+    // make a model for a well-sampled field of a simple function (f = ax + by + c)
+    # if (TEST1)
+    {
+	// function is defined over the range 0-1000, 0-1000
+	psImageBinning *binning = psImageBinningAlloc();
+	binning->nXfine = 1000;
+	binning->nYfine = 1000;
+	binning->nXruff = 10;
+	binning->nYruff = 10;
+
+	// generate a grid of test data points
+	psVector *x = psVectorAllocEmpty (100, PS_TYPE_F32);
+	psVector *y = psVectorAllocEmpty (100, PS_TYPE_F32);
+	psVector *f = psVectorAllocEmpty (100, PS_TYPE_F32);
+
+	for (int ix = 0; ix < 1000; ix += 20) {
+	    for (int iy = 0; iy < 1000; iy += 20) {
+		x->data.F32[x->n] = ix;
+		y->data.F32[y->n] = iy;
+		f->data.F32[f->n] = C00 + C10*ix + C01*iy;
+		psVectorExtend (x, 100, 1);
+		psVectorExtend (y, 100, 1);
+		psVectorExtend (f, 100, 1);
+	    }
+	}
+
+	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);
+	psFree (binning);
+
+	// fprintf (stderr, "nGood: %d\n", map->nGood);
+	// fprintf (stderr, "nPoor: %d\n", map->nPoor);
+	// fprintf (stderr, "nBad:  %d\n", map->nBad);
+	
+	SaveImage (NULL, map->map, "map.fits");
+	
+	psImage *field = psImageAlloc(1000, 1000, PS_TYPE_F32);
+	for (int ix = 0; ix < 1000; ix++) {
+	    for (int iy = 0; iy < 1000; iy++) {
+		field->data.F32[iy][ix] = C00 + C10*ix + C01*iy;
+	    }
+	}
+	SaveImage (NULL, field, "field.fits");
+
+	// measure difference between model (map) and data (field)
+	for (int ix = 0; ix < map->map->numCols; ix++) {
+	    for (int iy = 0; iy < map->map->numRows; iy++) {
+		int xo = psImageBinningGetFineX(map->binning, ix + 0.5);
+		int yo = psImageBinningGetFineY(map->binning, iy + 0.5);
+		float df = field->data.F32[yo][xo] - map->map->data.F32[iy][ix];
+		fprintf (stderr, "%d %d -> %d %d  :  %f = %f - %f\n", ix, iy, xo, yo, df, field->data.F32[yo][xo], map->map->data.F32[iy][ix]);
+	    }
+	}
+
+	psImage *model = psImageAlloc(1000, 1000, PS_TYPE_F32);
+	for (int ix = 0; ix < 1000; ix++) {
+	    for (int iy = 0; iy < 1000; iy++) {
+		model->data.F32[iy][ix] = psImageUnbinPixel_V2 (ix + 0.5, iy + 0.5, map->map, map->binning);
+	    }
+	}
+	SaveImage (NULL, model, "model.fits");
+    }
+    # endif
+
+    // make a model for a poorly-sampled field of a simple function (f = ax + by + c)
+
+    // choose a model scale for a poorly-sampled field of a simple function (f = ax + by + c)
+
+    // make a model for a well-sampled field of a high-order polynomial function (f = (a(x-xo)^2 + b(y-yo)^2) + c)
+    # if (TEST2)
+    {
+	// function is defined over the range 0-1000, 0-1000
+	psImageBinning *binning = psImageBinningAlloc();
+	binning->nXfine = 1000;
+	binning->nYfine = 1000;
+	binning->nXruff = 5;
+	binning->nYruff = 5;
+
+	// generate a grid of test data points
+	psVector *x = psVectorAllocEmpty (100, PS_TYPE_F32);
+	psVector *y = psVectorAllocEmpty (100, PS_TYPE_F32);
+	psVector *f = psVectorAllocEmpty (100, PS_TYPE_F32);
+
+	for (int ix = 0; ix < 1000; ix += 20) {
+	    for (int iy = 0; iy < 1000; iy += 20) {
+		x->data.F32[x->n] = ix;
+		y->data.F32[y->n] = iy;
+		f->data.F32[f->n] = C10*PS_SQR(ix-500) + C01*PS_SQR(iy-500);
+		psVectorExtend (x, 100, 1);
+		psVectorExtend (y, 100, 1);
+		psVectorExtend (f, 100, 1);
+	    }
+	}
+
+	psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+
+	// scale defines both field and map image sizes (nXfine, nXruff)
+	psImageMap *map = psImageMapAlloc (NULL, binning, stats);
+
+	// allocate a map, but we have no field image to supply
+	// XXX this function needs to correct for the mean superpixel position 
+	// which is sampled...
+	psImageMapGenerate (map, x, y, f, 0.1);
+	psFree (binning);
+
+	fprintf (stderr, "nGood: %d\n", map->nGood);
+	fprintf (stderr, "nPoor: %d\n", map->nPoor);
+	fprintf (stderr, "nBad:  %d\n", map->nBad);
+	
+	SaveImage (NULL, map->map, "map.fits");
+	
+	psImage *field = psImageAlloc(1000, 1000, PS_TYPE_F32);
+	for (int ix = 0; ix < 1000; ix++) {
+	    for (int iy = 0; iy < 1000; iy++) {
+		field->data.F32[iy][ix] = C10*PS_SQR(ix-500) + C01*PS_SQR(iy-500);
+	    }
+	}
+	SaveImage (NULL, field, "field.fits");
+
+	// measure difference between model (map) and data (field)
+	for (int ix = 0; ix < map->map->numCols; ix++) {
+	    for (int iy = 0; iy < map->map->numRows; iy++) {
+		int xo = psImageBinningGetFineX(map->binning, ix + 0.5);
+		int yo = psImageBinningGetFineY(map->binning, iy + 0.5);
+		float df = field->data.F32[yo][xo] - map->map->data.F32[iy][ix];
+		fprintf (stderr, "%d %d -> %d %d  :  %f = %f - %f\n", ix, iy, xo, yo, df, field->data.F32[yo][xo], map->map->data.F32[iy][ix]);
+	    }
+	}
+    }
+    # endif
+
+    // make a model for a well-sampled field of a high-order polynomial function (f = (a(x-xo)^2 + b(y-yo)^2) + c)
+    # if (TEST3)
+    {
+	// function is defined over the range 0-1000, 0-1000
+	psImageBinning *binning = psImageBinningAlloc();
+	binning->nXfine = 1000;
+	binning->nYfine = 1000;
+	binning->nXruff = 5;
+	binning->nYruff = 5;
+
+	// generate a grid of test data points
+	psVector *x = psVectorAllocEmpty (100, PS_TYPE_F32);
+	psVector *y = psVectorAllocEmpty (100, PS_TYPE_F32);
+	psVector *f = psVectorAllocEmpty (100, PS_TYPE_F32);
+
+	for (int ix = 0; ix < 1000; ix += 20) {
+	    for (int iy = 0; iy < 1000; iy += 20) {
+		x->data.F32[x->n] = ix;
+		y->data.F32[y->n] = iy;
+		f->data.F32[f->n] = C10*(ix-500) + PS_SQR(C01*(iy-500));
+		psVectorExtend (x, 100, 1);
+		psVectorExtend (y, 100, 1);
+		psVectorExtend (f, 100, 1);
+	    }
+	}
+
+	psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+
+	// scale defines both field and map image sizes (nXfine, nXruff)
+	psImageMap *map = psImageMapAlloc (NULL, binning, stats);
+
+	// allocate a map, but we have no field image to supply
+	// XXX this function needs to correct for the mean superpixel position 
+	// which is sampled...
+	psImageMapGenerate (map, x, y, f, 0.1);
+	psFree (binning);
+
+	fprintf (stderr, "nGood: %d\n", map->nGood);
+	fprintf (stderr, "nPoor: %d\n", map->nPoor);
+	fprintf (stderr, "nBad:  %d\n", map->nBad);
+	
+	SaveImage (NULL, map->map, "map.fits");
+	
+	psImage *field = psImageAlloc(1000, 1000, PS_TYPE_F32);
+	for (int ix = 0; ix < 1000; ix++) {
+	    for (int iy = 0; iy < 1000; iy++) {
+		field->data.F32[iy][ix] = C10*(ix-500) + PS_SQR(C01*(iy-500));
+	    }
+	}
+	SaveImage (NULL, field, "field.fits");
+
+	// measure difference between model (map) and data (field)
+	for (int ix = 0; ix < map->map->numCols; ix++) {
+	    for (int iy = 0; iy < map->map->numRows; iy++) {
+		int xo = psImageBinningGetFineX(map->binning, ix + 0.5);
+		int yo = psImageBinningGetFineY(map->binning, iy + 0.5);
+		float df = field->data.F32[yo][xo] - map->map->data.F32[iy][ix];
+		fprintf (stderr, "%d %d -> %d %d  :  %f = %f - %f\n", ix, iy, xo, yo, df, field->data.F32[yo][xo], map->map->data.F32[iy][ix]);
+	    }
+	}
+
+
+	psImage *model = psImageAlloc(1000, 1000, PS_TYPE_F32);
+	for (int ix = 0; ix < 1000; ix++) {
+	    for (int iy = 0; iy < 1000; iy++) {
+		// XXX the binned coordinates are probably off by 0.5 pix, or we are interpolating 
+		// to the wrong binned coordinate.  
+		// XXX fix edge cases
+		model->data.F32[iy][ix] = psImageUnbinPixel_V2 (ix + 0.5, iy + 0.5, map->map, map->binning);
+	    }
+	}
+	SaveImage (NULL, model, "model.fits");
+    }
+    # endif
+
+    // make a model for a poorly-sampled field of a non-polynomial function (f = (a(x-xo)^2 + b(y-yo)^2)^-1)
+
+    // choose a model scale for a poorly-sampled field of a non-polynomial function (f = (a(x-xo)^2 + b(y-yo)^2)^-1)
+