Index: trunk/psLib/src/imageops/psImageBinning.c
===================================================================
--- trunk/psLib/src/imageops/psImageBinning.c	(revision 14982)
+++ trunk/psLib/src/imageops/psImageBinning.c	(revision 14983)
@@ -8,6 +8,6 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-20 23:53:46 $
+ *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-22 03:05:50 $
  *
  *  Copyright 2007 Institute for Astronomy, University of Hawaii
@@ -19,4 +19,5 @@
 
 #include <stdio.h>
+#include "psMemory.h"
 #include "psError.h"
 #include "psAbort.h"
@@ -38,5 +39,5 @@
 
 void psImageBinningSetRuffSize(psImageBinning *binning, psImageBinningAlign align) {
-    
+
     assert (binning->nXfine > 0);
     assert (binning->nYfine > 0);
@@ -54,17 +55,17 @@
     switch (align) {
       case PS_IMAGE_BINNING_LEFT:
-	binning->nXoff = 0;
-	binning->nYoff = 0;
-	break;
+        binning->nXoff = 0;
+        binning->nYoff = 0;
+        break;
       case PS_IMAGE_BINNING_CENTER:
-	binning->nXoff = (binning->nXruff * binning->nXbin - binning->nXfine) / 2;
-	binning->nYoff = (binning->nYruff * binning->nYbin - binning->nYfine) / 2;
-	break;
+        binning->nXoff = (binning->nXruff * binning->nXbin - binning->nXfine) / 2;
+        binning->nYoff = (binning->nYruff * binning->nYbin - binning->nYfine) / 2;
+        break;
       case PS_IMAGE_BINNING_RIGHT:
-	binning->nXoff = (binning->nXruff * binning->nXbin - binning->nXfine);
-	binning->nYoff = (binning->nYruff * binning->nYbin - binning->nYfine);
-	break;
+        binning->nXoff = (binning->nXruff * binning->nXbin - binning->nXfine);
+        binning->nYoff = (binning->nYruff * binning->nYbin - binning->nYfine);
+        break;
       default:
-	psAbort ("programming error in %s: impossible case\n", __func__);
+        psAbort ("programming error in %s: impossible case\n", __func__);
     }
     return;
@@ -81,9 +82,9 @@
 
     if (image != NULL) {
-	binning->nXskip = image->col0 - binning->nXoff;
-	binning->nYskip = image->row0 - binning->nYoff;
+        binning->nXskip = image->col0 - binning->nXoff;
+        binning->nYskip = image->row0 - binning->nYoff;
     } else {
-	binning->nXskip = 0 - binning->nXoff;
-	binning->nYskip = 0 - binning->nYoff;
+        binning->nXskip = 0 - binning->nXoff;
+        binning->nYskip = 0 - binning->nYoff;
     }
     return;
@@ -91,5 +92,5 @@
 
 void psImageBinningSetScale(psImageBinning *binning, psImageBinningAlign align) {
-    
+
     assert (binning->nXfine > 0);
     assert (binning->nYfine > 0);
@@ -107,17 +108,17 @@
     switch (align) {
       case PS_IMAGE_BINNING_LEFT:
-	binning->nXoff = 0;
-	binning->nYoff = 0;
-	break;
+        binning->nXoff = 0;
+        binning->nYoff = 0;
+        break;
       case PS_IMAGE_BINNING_CENTER:
-	binning->nXoff = (binning->nXruff * binning->nXbin - binning->nXfine) / 2;
-	binning->nYoff = (binning->nYruff * binning->nYbin - binning->nYfine) / 2;
-	break;
+        binning->nXoff = (binning->nXruff * binning->nXbin - binning->nXfine) / 2;
+        binning->nYoff = (binning->nYruff * binning->nYbin - binning->nYfine) / 2;
+        break;
       case PS_IMAGE_BINNING_RIGHT:
-	binning->nXoff = (binning->nXruff * binning->nXbin - binning->nXfine);
-	binning->nYoff = (binning->nYruff * binning->nYbin - binning->nYfine);
-	break;
+        binning->nXoff = (binning->nXruff * binning->nXbin - binning->nXfine);
+        binning->nYoff = (binning->nYruff * binning->nYbin - binning->nYfine);
+        break;
       default:
-	psAbort ("programming error in %s: impossible case\n", __func__);
+        psAbort ("programming error in %s: impossible case\n", __func__);
     }
     return;
Index: trunk/psLib/src/imageops/psImageMap.c
===================================================================
--- trunk/psLib/src/imageops/psImageMap.c	(revision 14982)
+++ trunk/psLib/src/imageops/psImageMap.c	(revision 14983)
@@ -7,6 +7,6 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-20 23:54:25 $
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-22 03:05:50 $
  *
  *  Copyright 2007 Institute for Astronomy, University of Hawaii
@@ -18,4 +18,5 @@
 
 #include <stdio.h>
+
 #include "psError.h"
 #include "psAbort.h"
@@ -26,4 +27,5 @@
 #include "psFitsImage.h"
 
+#include "psMemory.h"
 #include "psVector.h"
 #include "psImage.h"
@@ -96,22 +98,22 @@
 
     // accumulate the values for each map pixel
-    
+
     // we can do this by accumulating a vector of pixel indexes for each cell
     psArray *pixelSets = psArrayAlloc (map->map->numCols*map->map->numRows);
     for (int i = 0; i < pixelSets->n; i++) {
-	pixelSets->data[i] = psVectorAllocEmpty (4, PS_TYPE_S32);
+        pixelSets->data[i] = psVectorAllocEmpty (4, PS_TYPE_S32);
     }
     // associate each value with a cell
     for (int i = 0; i < x->n; i++) {
-	int xRuff = psImageBinningGetRuffX (map->binning, x->data.F32[i]);
-	int yRuff = psImageBinningGetRuffY (map->binning, y->data.F32[i]);
-
-	int bin = xRuff + yRuff*map->map->numCols;
-	assert (bin >= 0);
-	assert (bin < pixelSets->n);
-	
-	psVector *pixels = pixelSets->data[bin];
-	pixels->data.S32[pixels->n] = i;
-	psVectorExtend (pixels, 4, 1);
+        int xRuff = psImageBinningGetRuffX (map->binning, x->data.F32[i]);
+        int yRuff = psImageBinningGetRuffY (map->binning, y->data.F32[i]);
+
+        int bin = xRuff + yRuff*map->map->numCols;
+        assert (bin >= 0);
+        assert (bin < pixelSets->n);
+
+        psVector *pixels = pixelSets->data[bin];
+        pixels->data.S32[pixels->n] = i;
+        psVectorExtend (pixels, 4, 1);
     }
 
@@ -123,59 +125,59 @@
     int Ny = map->map->numRows;
     for (int iy = 0; iy < Ny; iy++) {
-	for (int ix = 0; ix < Nx; ix++) {
-    
-	    // pixel index for this cell
-	    psVector *pixels = pixelSets->data[ix + iy*Nx];
-
-	    // storage vectors
-	    psVector *xCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
-	    psVector *yCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
-	    psVector *fCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
-
-	    // error vector, if needed
-	    psVector *dfCell = NULL;
-	    if (df) {
-		dfCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
-	    }		
-
-	    // collect data for this cell
-	    for (int i = 0; i < pixels->n; i++) {
-		int bin = pixels->data.S32[i];
-		// convert x,y in the fine image to the ruff image
-		xCell->data.F32[i]  = psImageBinningGetRuffX (map->binning, x->data.F32[bin]);
-		yCell->data.F32[i]  = psImageBinningGetRuffY (map->binning, y->data.F32[bin]);
-		fCell->data.F32[i]  = f->data.F32[bin];
-		if (df) {
-		    dfCell->data.F32[i] = df->data.F32[bin];
-		}
-	    }
-
-	    // reset the stats to avoid contamination from the previous loop
-	    psStatsInit (map->stats);
-
-	    // get the value
-	    // XXX need to supply a mask and skip the masked pixels when calculating the centroid
-	    // this will not in general be properly weighted...
-	    if (psVectorStats (map->stats, fCell, dfCell, NULL, 0)) {
-		mask->data.U8[iy][ix] = 0;
-		// XXX ensure only one option is selected, or save both position and width
-		map->map->data.F32[iy][ix] = psStatsGetValue (map->stats, map->stats->options); 
-
-		// calculate the mean position and save:
-		psStatsInit (meanStat);
-		psVectorStats (meanStat, xCell, NULL, NULL, 0);
-		xCoord->data.F32[iy][ix] = psStatsGetValue (meanStat, meanStat->options); 
-		psStatsInit (meanStat);
-		psVectorStats (meanStat, yCell, NULL, NULL, 0);
-		yCoord->data.F32[iy][ix] = psStatsGetValue (meanStat, meanStat->options); 
-	    } else {
-		mask->data.U8[iy][ix] = 1;
-	    }
-	    
-	    psFree (xCell);
-	    psFree (yCell);
-	    psFree (fCell);
-	    psFree (dfCell);
-	}
+        for (int ix = 0; ix < Nx; ix++) {
+
+            // pixel index for this cell
+            psVector *pixels = pixelSets->data[ix + iy*Nx];
+
+            // storage vectors
+            psVector *xCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
+            psVector *yCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
+            psVector *fCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
+
+            // error vector, if needed
+            psVector *dfCell = NULL;
+            if (df) {
+                dfCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
+            }
+
+            // collect data for this cell
+            for (int i = 0; i < pixels->n; i++) {
+                int bin = pixels->data.S32[i];
+                // convert x,y in the fine image to the ruff image
+                xCell->data.F32[i]  = psImageBinningGetRuffX (map->binning, x->data.F32[bin]);
+                yCell->data.F32[i]  = psImageBinningGetRuffY (map->binning, y->data.F32[bin]);
+                fCell->data.F32[i]  = f->data.F32[bin];
+                if (df) {
+                    dfCell->data.F32[i] = df->data.F32[bin];
+                }
+            }
+
+            // reset the stats to avoid contamination from the previous loop
+            psStatsInit (map->stats);
+
+            // get the value
+            // XXX need to supply a mask and skip the masked pixels when calculating the centroid
+            // this will not in general be properly weighted...
+            if (psVectorStats (map->stats, fCell, dfCell, NULL, 0)) {
+                mask->data.U8[iy][ix] = 0;
+                // XXX ensure only one option is selected, or save both position and width
+                map->map->data.F32[iy][ix] = psStatsGetValue (map->stats, map->stats->options);
+
+                // calculate the mean position and save:
+                psStatsInit (meanStat);
+                psVectorStats (meanStat, xCell, NULL, NULL, 0);
+                xCoord->data.F32[iy][ix] = psStatsGetValue (meanStat, meanStat->options);
+                psStatsInit (meanStat);
+                psVectorStats (meanStat, yCell, NULL, NULL, 0);
+                yCoord->data.F32[iy][ix] = psStatsGetValue (meanStat, meanStat->options);
+            } else {
+                mask->data.U8[iy][ix] = 1;
+            }
+
+            psFree (xCell);
+            psFree (yCell);
+            psFree (fCell);
+            psFree (dfCell);
+        }
     }
     psFree (pixelSets);
@@ -193,8 +195,8 @@
     map->nGood = mask->numCols * mask->numRows - map->nBad - map->nPoor;
     if (map->nBad > badFrac * mask->numCols * mask->numRows) {
-	psFree (xCoord);
-	psFree (yCoord);
-	psFree (mask);
-	return false;
+        psFree (xCoord);
+        psFree (yCoord);
+        psFree (mask);
+        return false;
     }
 
@@ -222,42 +224,42 @@
 // using the points given, generate a map with maximum resolution that yields only good and ok pixels
 bool psImageMapGenerateScale (psImageMap *map, psVector *x, psVector *y, psVector *f, psVector *df, float badFrac) {
-    
+
     int nXruff, nYruff;
-	
+
     while (!psImageMapGenerate (map, x, y, f, df, badFrac)) {
-	// if we failed to build an acceptable map, decrease nXruff, nYruff as appropriate, and
-	// try again...  try to keep the aspect ratio.
-
-	// if both axes are at 1, give up
-	if ((map->binning->nXruff == 1) && (map->binning->nYruff == 1)) {
-	    return false;
-	}
-
-	// if one axis is at 1, decrement the other
-	if (map->binning->nXruff == 1) {
-	    nXruff = map->binning->nXruff;
-	    nYruff = map->binning->nYruff - 1;
-	    psImageMapModifyScale (map, nXruff, nYruff);
-	    continue;
-	}
-	if (map->binning->nYruff == 1) {
-	    nYruff = map->binning->nYruff;
-	    nXruff = map->binning->nXruff - 1;
-	    psImageMapModifyScale (map, nXruff, nYruff);
-	    continue;
-	}
-
-	// otherwise, decrement the larger axis, and set the smaller based
-	// on the aspect ratio
-	float aRatio = map->binning->nXruff / map->binning->nYruff;
-	if (map->binning->nXruff > map->binning->nYruff) {
-	    nXruff = map->binning->nXruff - 1;
-	    nYruff = (int)(0.5 + (nXruff / aRatio));
-	} else {
-	    nYruff = map->binning->nYruff - 1;
-	    nXruff = (int)(0.5 + (nYruff * aRatio));
-	}
-
-	psImageMapModifyScale (map, nXruff, nYruff);
+        // if we failed to build an acceptable map, decrease nXruff, nYruff as appropriate, and
+        // try again...  try to keep the aspect ratio.
+
+        // if both axes are at 1, give up
+        if ((map->binning->nXruff == 1) && (map->binning->nYruff == 1)) {
+            return false;
+        }
+
+        // if one axis is at 1, decrement the other
+        if (map->binning->nXruff == 1) {
+            nXruff = map->binning->nXruff;
+            nYruff = map->binning->nYruff - 1;
+            psImageMapModifyScale (map, nXruff, nYruff);
+            continue;
+        }
+        if (map->binning->nYruff == 1) {
+            nYruff = map->binning->nYruff;
+            nXruff = map->binning->nXruff - 1;
+            psImageMapModifyScale (map, nXruff, nYruff);
+            continue;
+        }
+
+        // otherwise, decrement the larger axis, and set the smaller based
+        // on the aspect ratio
+        float aRatio = map->binning->nXruff / map->binning->nYruff;
+        if (map->binning->nXruff > map->binning->nYruff) {
+            nXruff = map->binning->nXruff - 1;
+            nYruff = (int)(0.5 + (nXruff / aRatio));
+        } else {
+            nYruff = map->binning->nYruff - 1;
+            nXruff = (int)(0.5 + (nYruff * aRatio));
+        }
+
+        psImageMapModifyScale (map, nXruff, nYruff);
     }
     return true;
@@ -271,5 +273,5 @@
     result = psImageUnbinPixel_V2(x, y, map->map, map->binning);
 
-    return result; 
+    return result;
 }
 
@@ -284,7 +286,7 @@
 
     for (int i = 0; i < x->n; i++) {
-	result->data.F32[i] = psImageUnbinPixel_V2(x->data.F32[i], y->data.F32[i], map->map, map->binning);
-    }
-
-    return result; 
-}
+        result->data.F32[i] = psImageUnbinPixel_V2(x->data.F32[i], y->data.F32[i], map->map, map->binning);
+    }
+
+    return result;
+}
Index: trunk/psLib/src/imageops/psImageMapFit.c
===================================================================
--- trunk/psLib/src/imageops/psImageMapFit.c	(revision 14982)
+++ trunk/psLib/src/imageops/psImageMapFit.c	(revision 14983)
@@ -7,6 +7,6 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-21 02:45:33 $
+ *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-22 03:05:50 $
  *
  *  Copyright 2007 Institute for Astronomy, University of Hawaii
@@ -18,4 +18,5 @@
 
 #include <stdio.h>
+#include "psMemory.h"
 #include "psError.h"
 #include "psAbort.h"
@@ -57,22 +58,22 @@
     // no spatial information, just calculate mean & stdev
     if ((Nx == 1) && (Ny == 1)) {
-	psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
-	
-	// XXX does ROBUST_MEDIAN work with weight?
-	psVectorStats (stats, f, NULL, mask, maskValue);
-
-	map->map->data.F32[0][0]   = stats->robustMedian;
-	map->error->data.F32[0][0] = stats->robustStdev;
-	psFree (stats);
-	return true;
+        psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+
+        // XXX does ROBUST_MEDIAN work with weight?
+        psVectorStats (stats, f, NULL, mask, maskValue);
+
+        map->map->data.F32[0][0]   = stats->robustMedian;
+        map->error->data.F32[0][0] = stats->robustStdev;
+        psFree (stats);
+        return true;
     }
 
     if (Nx == 1) {
-	psAbort ("un-implemented edge case");
-	goto insert;
+        psAbort ("un-implemented edge case");
+        goto insert;
     }
     if (Ny == 1) {
-	psAbort ("un-implemented edge case");
-	goto insert;
+        psAbort ("un-implemented edge case");
+        goto insert;
     }
 
@@ -82,11 +83,11 @@
 
     for (int i = 0; i < 3; i++) {
-	SAv[i] = SAm[i] + 1;
-	// TAv[i] = TAm[i] + 1;
+        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 
+    // 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);
@@ -95,158 +96,158 @@
     psVectorInit (B, 0.0);
 
-    // we are looping over the Nx,Ny image map elements; 
-    // the matrix equation contains Nx*Ny rows and columns 
+    // 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++) {
-
-		if (mask && (mask->data.U8[i] & maskValue)) continue;
-
-		// 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);
-
-		// edge cases to include:
-		bool edgeX = false;
-		edgeX |= ((n == 1) && (dx < -1.0));
-		edgeX |= ((n == Nx - 2) && (dx > +1.0));
-
-		bool edgeY = false;
-		edgeY |= ((m == 1) && (dy < -1.0));
-		edgeY |= ((m == Ny - 2) && (dy > +1.0));
-		
-		// skip points outside of 2x2 grid centered on n,m:
-		if (!edgeX && (fabs(dx) > 1.0)) continue;
-		if (!edgeY && (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
-
-		int Qx = (dx >= 0) ? 1 : 0;
-		if (n ==      0) Qx = 1;
-		if (n == Nx - 1) Qx = 0;
-
-		int Qy = (dy >= 0) ? 1 : 0;
-		if (m ==      0) Qy = 1;
-		if (m == Ny - 1) Qy = 0;
-
-		// points at offset 1,1
-		if ((Qx == 1) && (Qy == 1)) {
-		    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 ((Qx == 1) && (Qy == 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 ((Qx == 0) && (Qy == 1)) {
-		    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 ((Qx == 0) && (Qy == 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;
-
-	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] = sA[jn][jm];
-		    // fprintf (stderr, "A %d %d (%d %d : %d %d): %f\n", I, J, n, m, n + jn, m + jm, sA[jn][jm]);
-		    // Sum += sA[jn][jm];
-		}
-	    }
-	    // fprintf (stderr, "B %d (%d %d) : %f  :  %f\n", I, n, m, B->data.F32[I], Sum);
-	    // Total += Sum;
-	}
+        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++) {
+
+                if (mask && (mask->data.U8[i] & maskValue)) continue;
+
+                // 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);
+
+                // edge cases to include:
+                bool edgeX = false;
+                edgeX |= ((n == 1) && (dx < -1.0));
+                edgeX |= ((n == Nx - 2) && (dx > +1.0));
+
+                bool edgeY = false;
+                edgeY |= ((m == 1) && (dy < -1.0));
+                edgeY |= ((m == Ny - 2) && (dy > +1.0));
+
+                // skip points outside of 2x2 grid centered on n,m:
+                if (!edgeX && (fabs(dx) > 1.0)) continue;
+                if (!edgeY && (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
+
+                int Qx = (dx >= 0) ? 1 : 0;
+                if (n ==      0) Qx = 1;
+                if (n == Nx - 1) Qx = 0;
+
+                int Qy = (dy >= 0) ? 1 : 0;
+                if (m ==      0) Qy = 1;
+                if (m == Ny - 1) Qy = 0;
+
+                // points at offset 1,1
+                if ((Qx == 1) && (Qy == 1)) {
+                    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 ((Qx == 1) && (Qy == 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 ((Qx == 0) && (Qy == 1)) {
+                    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 ((Qx == 0) && (Qy == 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;
+
+        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] = sA[jn][jm];
+                    // fprintf (stderr, "A %d %d (%d %d : %d %d): %f\n", I, J, n, m, n + jn, m + jm, sA[jn][jm]);
+                    // Sum += sA[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);
@@ -256,13 +257,13 @@
     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;
-	}
+        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;
+        }
     }
 
@@ -274,5 +275,5 @@
     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];
+        vector->data.F32[0][n] = B->data.F32[n];
     }
 
@@ -284,25 +285,25 @@
 
     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;
-    }
-    
+        psAbort ("failed on linear equations");
+        psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
+        psFree (A);
+        psFree (B);
+        return false;
+    }
+
     // set bad values to NaN
     for (int i = 0; i < Nx*Ny; i++) {
-	if (Empty->data.S8[i]) {
-	    B->data.F32[i] = NAN;
-	}
+        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]);
-	}
+        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]);
+        }
     }
 
@@ -353,6 +354,6 @@
     psVector *mask = inMask;
     if (!inMask) {
-	mask = psVectorAlloc (x->n, PS_TYPE_U8);
-	psVectorInit (mask, 0);
+        mask = psVectorAlloc (x->n, PS_TYPE_U8);
+        psVectorInit (mask, 0);
     }
 
@@ -368,6 +369,6 @@
         if (!psImageMapFit(map, mask, maskValue, x, y, f, df)) {
             psError(PS_ERR_UNKNOWN, false, "Could not fit image map.\n");
-	    psFree(resid);
-	    if (!inMask) psFree (mask);
+            psFree(resid);
+            if (!inMask) psFree (mask);
             return false;
         }
@@ -377,9 +378,9 @@
             psError(PS_ERR_UNKNOWN, false, "Failure in psImageMapEvalVector().\n");
             psFree(resid);
-	    if (!inMask) psFree (mask);
+            if (!inMask) psFree (mask);
             return false;
         }
         for (int i = 0 ; i < f->n ; i++) {
-	    resid->data.F32[i] = (f->data.F32[i] - fit->data.F32[i]);
+            resid->data.F32[i] = (f->data.F32[i] - fit->data.F32[i]);
         }
 
@@ -388,5 +389,5 @@
             psFree(resid);
             psFree(fit);
-	    if (!inMask) psFree (mask);
+            if (!inMask) psFree (mask);
             return false;
         }
@@ -404,5 +405,5 @@
         // recovery is not allowed with this scheme
         for (psS32 i = 0; i < resid->n; i++) {
-	    // XXX this prevents recovery of previously masked values
+            // XXX this prevents recovery of previously masked values
             if (mask->data.U8[i] & maskValue) {
                 continue;
@@ -410,6 +411,6 @@
 
             if ((resid->data.F32[i] - meanValue > maxClipValue) || (resid->data.F32[i] - meanValue < minClipValue)) {
-		psTrace("psLib.imageops", 6, "Masking element %d  : %f vs %f : resid is %f\n", i, f->data.F32[i], fit->data.F32[i], resid->data.F32[i]);
-		mask->data.U8[i] |= 0x01;
+                psTrace("psLib.imageops", 6, "Masking element %d  : %f vs %f : resid is %f\n", i, f->data.F32[i], fit->data.F32[i], resid->data.F32[i]);
+                mask->data.U8[i] |= 0x01;
                 continue;
             }
Index: trunk/psLib/src/types/psMetadataConfig.c
===================================================================
--- trunk/psLib/src/types/psMetadataConfig.c	(revision 14982)
+++ trunk/psLib/src/types/psMetadataConfig.c	(revision 14983)
@@ -11,6 +11,6 @@
 *  @author Joshua Hoblitt, University of Hawaii 2006-2007
 *
-*  @version $Revision: 1.139 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2007-08-29 00:27:13 $
+*  @version $Revision: 1.140 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2007-09-22 03:05:50 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -31,4 +31,5 @@
 #include <unistd.h>
 
+#include "psMemory.h"
 #include "psMetadataConfig.h"
 #include "psAssert.h"