Index: trunk/psLib/src/imageops/psImageMap.c
===================================================================
--- trunk/psLib/src/imageops/psImageMap.c	(revision 14924)
+++ trunk/psLib/src/imageops/psImageMap.c	(revision 14924)
@@ -0,0 +1,290 @@
+/** @file  psImageMap.c
+ *
+ *  @brief Functions define a 2d coarse representation of a finer 2D field
+ *
+ *  @ingroup Image
+ *
+ *  @author Eugene Magnier, IfA
+ *
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-20 23:54:25 $
+ *
+ *  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"
+
+#include "psFits.h"
+#include "psAssert.h"
+#include "psRegion.h"
+#include "psFitsImage.h"
+
+#include "psVector.h"
+#include "psImage.h"
+#include "psStats.h"
+#include "psImageBinning.h"
+#include "psImageMap.h"
+#include "psImagePixelInterpolate.h"
+#include "psImageUnbin.h"
+
+static void psImageMapFree(psImageMap *map) {
+
+    if (!map) return;
+
+    psFree (map->map);
+    psFree (map->error);
+    psFree (map->field); // XXX ??? this should be freed here, but that causes an error...
+    psFree (map->stats);
+    psFree (map->binning);
+
+    return;
+}
+
+psImageMap *psImageMapAlloc(psImage *field, psImageBinning *binning, psStats *stats) {
+
+    assert (binning);
+    assert (stats);
+
+    psImageMap *map = (psImageMap*)psAlloc(sizeof(psImageMap));
+    psMemSetDeallocator(map, (psFreeFunc)psImageMapFree);
+
+    map->binning = psMemIncrRefCounter (binning);
+    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);
+    psImageBinningSetSkip (map->binning, map->field);
+
+    return map;
+}
+
+bool psImageMapModifyScale(psImageMap *map, int nXruff, int nYruff) {
+
+    assert (map);
+
+    map->binning->nXruff = nXruff;
+    map->binning->nYruff = nYruff;
+
+    psImageRecycle (map->map, nXruff, nYruff, PS_TYPE_F32);
+    psImageRecycle (map->error, nXruff, nYruff, PS_TYPE_F32);
+
+    psImageBinningSetScale (map->binning, PS_IMAGE_BINNING_CENTER);
+    psImageBinningSetSkip (map->binning, map->field);
+
+    return true;
+}
+
+// generate a psImageMap (or NULL) with the given number of superpixels in X and Y
+// this function returns an error if the output map has impossible holes
+bool psImageMapGenerate (psImageMap *map, psVector *x, psVector *y, psVector *f, psVector *df, float badFrac) {
+
+    psImage *mask = psImageAlloc (map->map->numCols, map->map->numRows, PS_TYPE_MASK);
+    psImage *xCoord = psImageAlloc (map->map->numCols, map->map->numRows, PS_TYPE_F32);
+    psImage *yCoord = psImageAlloc (map->map->numCols, map->map->numRows, PS_TYPE_F32);
+
+    // 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);
+    }
+    // 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);
+    }
+
+    // stats structure for getting the position centers
+    psStats *meanStat = psStatsAlloc (PS_STAT_SAMPLE_MEAN);
+
+    // accumulate the x,y coords for each point to calculate the mean position.
+    int Nx = map->map->numCols;
+    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);
+	}
+    }
+    psFree (pixelSets);
+    psFree (meanStat);
+   // at this point, for each map pixel, we have (f,x,y), or the pixel is masked.
+
+    psFits *fits = NULL;
+
+    fits = psFitsOpen ("imageMap.raw.fits", "w");
+    psFitsWriteImage (fits, NULL, map->map, 0, NULL);
+    psFitsClose (fits);
+
+    // did this analysis succeed?  (enough good or OK pixels?)
+    psImage *state = psImagePixelInterpolateState (&map->nBad, &map->nPoor, mask, 0xff);
+    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;
+    }
+
+    // fit the valid pixels to (0,1,2) order polynomials, interpolate values to the pixel center
+    // XXX I need to be careful about the pixel coordinates: center is 0,0 or 0.5, 0.5?
+    psImagePixelInterpolateCenter (map->map, xCoord, yCoord, state, mask, 0xff);
+    psFree (xCoord);
+    psFree (yCoord);
+
+    fits = psFitsOpen ("imageMap.ref.fits", "w");
+    psFitsWriteImage (fits, NULL, map->map, 0, NULL);
+    psFitsClose (fits);
+
+    psImagePixelInterpolatePoor (map->map, state, mask, 0xff);
+
+    fits = psFitsOpen ("imageMap.fix.fits", "w");
+    psFitsWriteImage (fits, NULL, map->map, 0, NULL);
+    psFitsClose (fits);
+
+    psFree (state);
+    psFree (mask);
+    return true;
+}
+
+// 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);
+    }
+    return true;
+}
+
+// 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_V2(x, y, map->map, map->binning);
+
+    return result; 
+}
+
+psVector *psImageMapEvalVector (psImageMap *map, psVector *x, psVector *y) {
+
+    assert (x);
+    assert (y);
+    assert (x->n == y->n);
+    assert (map);
+
+    psVector *result = psVectorAlloc (x->n, PS_TYPE_F32);
+
+    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; 
+}
Index: trunk/psLib/src/imageops/psImageMap.h
===================================================================
--- trunk/psLib/src/imageops/psImageMap.h	(revision 14924)
+++ trunk/psLib/src/imageops/psImageMap.h	(revision 14924)
@@ -0,0 +1,58 @@
+/** @file  psImageMap.c
+ *
+ *  @brief Functions define a 2d coarse representation of a finer 2D field
+ *
+ *  @ingroup Image
+ *
+ *  @author Eugene Magnier, IfA
+ *
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-20 23:54:25 $
+ *
+ *  Copyright 2007 Institute for Astronomy, University of Hawaii
+ */
+
+#ifndef PS_IMAGE_MAP_H
+#define PS_IMAGE_MAP_H
+
+/// @addtogroup ImageOps Image Operations
+/// @{
+
+// a structure to describe the 2D variations of some quantity as a function of position the
+// variation is represented as a psImage which covers the field also represented as a psImage.
+// the map image pixels are superpixels of the field image.  values in the field are determined
+// by interpolating the map image.
+typedef struct {
+    psStats *stats;
+    psImage *map;
+    psImage *error;
+    psImage *field;
+    psImageBinning *binning;
+    int nBad;
+    int nPoor;
+    int nGood;
+} psImageMap;
+
+psImageMap *psImageMapAlloc(psImage *field, psImageBinning *binning, psStats *stats) PS_ATTR_MALLOC;
+
+bool psImageMapModifyScale(psImageMap *map, int nXruff, int nYruff);
+
+// generate a psImageMap (or NULL) with the given number of superpixels in X and Y
+bool psImageMapGenerate (psImageMap *map, psVector *x, psVector *y, psVector *f, psVector *df, float badFrac);
+
+bool psImageMapGenerateScale (psImageMap *map, psVector *x, psVector *y, psVector *f, psVector *df, float badFrac);
+
+// apply the psImageMap to the given coordinate (fine image pixels)
+double psImageMapEval (psImageMap *map, float x, float y);
+
+// apply the psImageMap to the given coordinate vectors (fine image pixels)
+psVector *psImageMapEvalVector (psImageMap *map, psVector *x, psVector *y);
+
+// fit the image map to a set of points
+bool psImageMapFit (psImageMap *map, psVector *mask, psMaskType maskValue, psVector *x, psVector *y, psVector *f, psVector *df);
+
+// fit the image map to a set of points
+bool psImageMapClipFit (psImageMap *map, psStats *stats, psVector *mask, psMaskType maskValue, psVector *x, psVector *y, psVector *f, psVector *df);
+
+/// @}
+#endif // #ifndef PS_IMAGE_MAP_H
Index: trunk/psLib/src/imageops/psImageMapFit.c
===================================================================
--- trunk/psLib/src/imageops/psImageMapFit.c	(revision 14924)
+++ trunk/psLib/src/imageops/psImageMapFit.c	(revision 14924)
@@ -0,0 +1,430 @@
+/** @file  psImageMapFit.c
+ *
+ *  @brief Functions define a 2d coarse representation of a finer 2D field
+ *
+ *  @ingroup Image
+ *
+ *  @author Eugene Magnier, IfA
+ *
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-20 23:54:25 $
+ *
+ *  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"
+#include "psTrace.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 *mask, psMaskType maskValue, psVector *x, psVector *y, psVector *f, psVector *df) {
+
+    int I, J;
+
+    // XXX Add Asserts
+
+    // dimensions of the output map image
+    int Nx = map->binning->nXruff;
+    int Ny = map->binning->nYruff;
+
+    // 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;
+    }
+
+    if (Nx == 1) {
+	psAbort ("un-implemented edge case");
+	goto insert;
+    }
+    if (Ny == 1) {
+	psAbort ("un-implemented edge case");
+	goto insert;
+    }
+
+    // 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);
+
+    psImageInit (A, 0.0);
+    psVectorInit (B, 0.0);
+
+    // 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;
+	}
+    }
+    // fprintf (stderr, "Total: %f\n", Total);
+
+    // 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;
+	}
+    }
+
+    # if (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);
+    # endif
+
+    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;
+    }
+    
+    // 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]);
+	}
+    }
+
+    psFree (A);
+    psFree (B);
+    psFree (Empty);
+
+    return true;
+}
+
+// measure residuals on each pass and clip outliers based on stats
+bool psImageMapClipFit (psImageMap *map, psStats *stats, psVector *inMask, psMaskType maskValue, psVector *x, psVector *y, psVector *f, psVector *df) {
+
+    // XXX add in full PS_ASSERTS
+    assert (map);
+    assert (stats);
+    assert (x);
+    assert (y);
+    assert (f);
+    assert (df);
+
+    // the user supplies one of various stats option pairs,
+    // determine the desired mean and stdev STATS options:
+    psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_MEAN_V3 | PS_STAT_FITTED_MEAN_V4);
+    psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV_V2 | PS_STAT_FITTED_STDEV_V3 | PS_STAT_FITTED_STDEV_V4);
+    if (!meanOption) {
+        psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected");
+        return false;
+    }
+    if (!stdevOption) {
+        psError(PS_ERR_UNKNOWN, true, "no valid stdev stats option selected");
+        return false;
+    }
+
+    // clipping range defined by min and max and/or clipSigma
+    psF32 minClipSigma;
+    psF32 maxClipSigma;
+    if (isfinite(stats->max)) {
+        maxClipSigma = fabs(stats->max);
+    } else {
+        maxClipSigma = fabs(stats->clipSigma);
+    }
+    if (isfinite(stats->min)) {
+        minClipSigma = fabs(stats->min);
+    } else {
+        minClipSigma = fabs(stats->clipSigma);
+    }
+
+    psVector *mask = inMask;
+    if (!inMask) {
+	mask = psVectorAlloc (x->n, PS_TYPE_U8);
+    }
+
+    // vector to store residuals
+    psVector *resid = psVectorAlloc(f->n, PS_TYPE_F32);
+
+    psTrace("psLib.imageops", 4, "stats->clipIter is %d\n", stats->clipIter);
+    psTrace("psLib.imageops", 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma);
+
+    for (psS32 N = 0; N < stats->clipIter; N++) {
+        psTrace("psLib.imageops", 6, "Loop iteration %d.  Calling psImageMapFit()\n", N);
+        psS32 Nkeep = 0;
+        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);
+            return false;
+        }
+
+        psVector *fit = psImageMapEvalVector(map, x, y);
+        if (fit == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Failure in psImageMapEvalVector().\n");
+            psFree(resid);
+	    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]);
+        }
+
+        if (!psVectorStats(stats, resid, NULL, mask, maskValue)) {
+            psError(PS_ERR_UNKNOWN, false, "Failure to compute statistics on the resid vector.\n");
+            psFree(resid);
+            psFree(fit);
+	    if (!inMask) psFree (mask);
+            return false;
+        }
+
+        double meanValue = psStatsGetValue (stats, meanOption);
+        double stdevValue = psStatsGetValue (stats, stdevOption);
+
+        psTrace("psLib.imageops", 5, "Mean is %f\n", meanValue);
+        psTrace("psLib.imageops", 5, "Stdev is %f\n", stdevValue);
+        psF32 minClipValue = -minClipSigma*stdevValue;
+        psF32 maxClipValue = +maxClipSigma*stdevValue;
+
+        // set mask if pts are not valid
+        // we are masking out any point which is out of range
+        // recovery is not allowed with this scheme
+        for (psS32 i = 0; i < resid->n; i++) {
+	    // XXX this prevents recovery of previously masked values
+            if (mask->data.U8[i] & maskValue) {
+                continue;
+            }
+
+            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;
+                continue;
+            }
+            Nkeep++;
+        }
+
+        // We should probably exit this loop if no new elements were masked since the fit won't
+        // change.
+        psTrace("psLib.imageops", 6, "keeping %d of %ld pts for fit\n", Nkeep, x->n);
+        stats->clippedNvalues = Nkeep;
+        psFree(fit);
+    }
+
+    // Free local temporary variables
+    psFree(resid);
+    if (!inMask) psFree (mask);
+    return true;
+}
Index: trunk/psLib/src/imageops/psImagePixelInterpolate.c
===================================================================
--- trunk/psLib/src/imageops/psImagePixelInterpolate.c	(revision 14924)
+++ trunk/psLib/src/imageops/psImagePixelInterpolate.c	(revision 14924)
@@ -0,0 +1,379 @@
+/** @file  psImagePixelInterpolate.c
+ *
+ *  @brief Functions for interpolating bad pixels in images
+ *
+ *  these functions test and set masked pixels in an image.  These functions are complementary
+ *  to the psImageInterpolate functions, which perform sub-pixel interpolation.  Those
+ *  functions require all pixels surrounding the sub-pixel interpolation to have values which
+ *  are valid.  These functions enable interpolation of complete missing pixels, potentially
+ *  across large spans.
+ *
+ *  @author Eugene Magnier, IfA
+ *
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-20 23:54:25 $
+ *
+ *  Copyright 2007 Institute for Astronomy, University of Hawaii
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+# include "config.h"
+#endif
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <strings.h>
+#include <string.h>
+
+#include "psAbort.h"
+#include "psMemory.h"
+#include "psError.h"
+#include "psAssert.h"
+#include "psString.h"
+#include "psPolynomial.h"
+#include "psPolynomialUtils.h"
+#include "psMinimizePolyFit.h"
+#include "psImage.h"
+#include "psImageInterpolate.h"
+#include "psImagePixelInterpolate.h"
+
+#include "psFits.h"
+#include "psFitsImage.h"
+
+# define PS_IMAGE_ITER_STATE(XS,XE,YS,YE,N_MIN,TYPE) \
+	    nGood = 0; \
+	    for (int jy = YS; jy <= YE; jy++) { \
+		/* stick to pixels in image grid */ \
+		if (jy + iy < 0) { continue; } \
+		if (jy + iy >= mask->numRows) { continue; } \
+		for (int jx = XS; jx <= XE; jx++) { \
+		    /* stick to pixels in image grid */ \
+		    if (jx + ix < 0) { continue; } \
+		    if (jx + ix >= mask->numCols) { continue; } \
+		    /* do not test self */ \
+		    if (!jx && !jy) { continue; } \
+		    if (mask->data.PS_TYPE_MASK_DATA[iy+jy][ix+jx] & maskVal) { continue; } \
+		    nGood ++; \
+		} \
+	    } \
+	    if (nGood >= N_MIN) {  \
+		nPoor ++; \
+		result->data.S32[iy][ix] = TYPE; \
+		continue; \
+	    }
+
+// count and mark pixels based on their potential for being interpolated.  the input image is
+// just the mask, the output image contains enum values which define the type of interpolation which
+// can be performed
+psImage *psImagePixelInterpolateState (int *nBad, int *nPoor, psImage *mask, psMaskType maskVal) {
+
+    psImage *result = psImageAlloc (mask->numCols, mask->numRows, PS_TYPE_S32);
+    psImageInit (result, 0);
+    
+    *nPoor = 0;
+    *nBad = 0;
+
+    for (int iy = 0; iy < mask->numRows; iy++) {
+	for (int ix = 0; ix < mask->numCols; ix++) {
+
+	    // state of the good pixels (unmasked)
+	    if (!(mask->data.PS_TYPE_MASK_DATA[iy][ix] & maskVal)) { 
+		// count good neighbor pixels (+ self)
+		int nGood = 0;
+		int minX = +1;
+		int maxX = -1;
+		int minY = +1;
+		int maxY = -1;
+		for (int jy = -1; jy <= +1; jy++) {
+		    /* stick to pixels in image grid */
+		    if (jy + iy < 0) { continue; }
+		    if (jy + iy >= mask->numRows) { continue; }
+		    for (int jx = -1; jx <= +1; jx++) {
+			/* stick to pixels in image grid */
+			if (jx + ix < 0) { continue; }
+			if (jx + ix >= mask->numCols) { continue; }
+			if (mask->data.PS_TYPE_MASK_DATA[iy+jy][ix+jx] & maskVal) { continue; }
+			nGood ++;
+			minX = PS_MIN (minX, jx);
+			maxX = PS_MAX (maxX, jx);
+			minY = PS_MIN (minY, jy);
+			maxY = PS_MAX (maxY, jy);
+		    }
+		}
+		int dX = maxX - minX;
+		int dY = maxY - minY;
+		// what type of local interpolation can we use?
+		if ((nGood >= 6) && (dX == 2) && (dY == 2)) {
+		    result->data.S32[iy][ix] = PS_IMAGE_INTERPOLATE_GOOD2;
+		    continue; 
+		}
+		if (nGood >= 3) {
+		    result->data.S32[iy][ix] = PS_IMAGE_INTERPOLATE_GOOD1;
+		    continue; 
+		}
+		result->data.S32[iy][ix] = PS_IMAGE_INTERPOLATE_GOOD0;
+		continue; 
+	    }
+
+	    // examine the neighbors.  If at least 6 of the 8 surrounding, or 3 of the 4 corner
+	    // neighbors are valid, this is a pixel which can be interpolated.
+
+	    int nGood;
+
+	    // check for poor pixels
+	    PS_IMAGE_ITER_STATE (-1,+1,-1,+1,6, PS_IMAGE_INTERPOLATE_CENTER);
+	    PS_IMAGE_ITER_STATE (-1,+0,-1,+0,3, PS_IMAGE_INTERPOLATE_UR);
+	    PS_IMAGE_ITER_STATE (-1,+0,+0,+1,3, PS_IMAGE_INTERPOLATE_LR);
+	    PS_IMAGE_ITER_STATE (+0,+1,-1,+0,3, PS_IMAGE_INTERPOLATE_UL);
+	    PS_IMAGE_ITER_STATE (+0,+1,+0,+1,3, PS_IMAGE_INTERPOLATE_LL);
+
+	    nBad ++;
+	    result->data.S32[iy][ix] = PS_IMAGE_INTERPOLATE_BAD;
+	}	    
+    }
+    return result;
+}
+
+// interpolate the poor pixels using the available options
+bool psImagePixelInterpolatePoor (psImage *image, psImage *state, psImage *mask, psMaskType maskVal) {
+
+    assert (image->numCols == state->numCols);
+    assert (image->numRows == state->numRows);
+    assert (image->numCols == mask->numCols);
+    assert (image->numRows == mask->numRows);
+
+    // allocate the vectors for the 2nd order fit below
+    psVector *f  = psVectorAlloc (9, PS_TYPE_F32);
+    // XXX if we add the weight above, include df
+    // psVector *df = psVectorAlloc (9, PS_TYPE_F32); 
+    psVector *x  = psVectorAlloc (9, PS_TYPE_F32);
+    psVector *y  = psVectorAlloc (9, PS_TYPE_F32);
+
+    // allocate a 2D polynomial to fit a quadratic to the valid neighbor pixels.
+    psPolynomial2D *poly = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 2, 2);
+    poly->mask[2][2] = 1;
+    poly->mask[2][1] = 1;
+    poly->mask[1][2] = 1;
+
+    for (int iy = 0; iy < state->numRows; iy++) {
+	for (int ix = 0; ix < state->numCols; ix++) {
+
+	    switch (state->data.S32[iy][ix]) {
+	      case PS_IMAGE_INTERPOLATE_GOOD0: 
+	      case PS_IMAGE_INTERPOLATE_GOOD1: 
+	      case PS_IMAGE_INTERPOLATE_GOOD2: 
+		// skip the good pixels
+		break; 
+
+	      case PS_IMAGE_INTERPOLATE_BAD:
+		// skip the bad pixels
+		break; 
+
+	      case PS_IMAGE_INTERPOLATE_CENTER: {
+		  // XXX is there a fit-image-region function?
+		  int n = 0;
+		  for (int jy = -1; jy <= +1; jy++) {
+		      // skip invalid pixels 
+		      if (jy + iy < 0) { continue; }
+		      if (jy + iy >= image->numRows) { continue; }
+		      for (int jx = -1; jx <= +1; jx++) {
+			  // skip invalid pixels 
+			  if (jx + ix < 0) { continue; } 
+			  if (jx + ix >= image->numCols) { continue; } 
+			  // skip self 
+			  if (!jx && !jy) { continue; } 
+			  // skip masked pixels
+			  if (mask->data.PS_TYPE_MASK_DATA[iy+jy][ix+jx] & maskVal) { continue; }
+			  x->data.F32[n] = jx;
+			  y->data.F32[n] = jy;
+			  f->data.F32[n] = image->data.F32[iy+jy][ix+jx];
+			  // df->data.F32[n] = weight->data.F32[iy+jy][ix+jx];
+			  n++;
+		      }
+		  }
+		  // set vector lengths here
+		  x->n = n;
+		  y->n = n;
+		  f->n = n;
+		  // df->n = n;
+		  // psVectorFitPolynomial2D (poly, NULL, 0xff, f, df, x, y);
+		  psVectorFitPolynomial2D (poly, NULL, 0xff, f, NULL, x, y);
+		  // apply the fitted quadratic to get the poor pixel value
+		  image->data.F32[iy][ix] = poly->coeff[0][0];
+		  break; }
+
+		// XXX should I use 1 1D polynomial fitting all unmasked pixels in the 3x3 grid?
+		// XXX that would automatically extend to regions where only 2 pixels are valid...
+	      case PS_IMAGE_INTERPOLATE_LL: {
+		  // fit a plane to the 3 pixels at (0,1),(1,0),(1,1), extend to pixel at (0,0)
+		  image->data.F32[iy][ix] = image->data.F32[iy+1][ix+1] - image->data.F32[iy+0][ix+1] - image->data.F32[iy+1][ix+0];
+		  break; }
+
+	      case PS_IMAGE_INTERPOLATE_LR: {
+		  // fit a plane to the 3 pixels at (0,1),(-1,0),(-1,1), extend to pixel at (0,0)
+		  image->data.F32[iy][ix] = image->data.F32[iy+1][ix-1] - image->data.F32[iy+0][ix-1] - image->data.F32[iy+1][ix+0];
+		  break; }
+
+	      case PS_IMAGE_INTERPOLATE_UL: {
+		  // fit a plane to the 3 pixels at (0,-1),(1,0),(1,-1), extend to pixel at (0,0)
+		  image->data.F32[iy][ix] = image->data.F32[iy-1][ix+1] - image->data.F32[iy+0][ix+1] - image->data.F32[iy-1][ix+0];
+		  break; }
+
+	      case PS_IMAGE_INTERPOLATE_UR: {
+		  // fit a plane to the 3 pixels at (0,-1),(-1,0),(-1,-1), extend to pixel at (0,0)
+		  image->data.F32[iy][ix] = image->data.F32[iy-1][ix-1] - image->data.F32[iy+0][ix-1] - image->data.F32[iy-1][ix+0];
+		  break; }
+
+	      default:
+		psAbort("impossible case in __func__");
+	    }
+	}	    
+    }
+
+    psFree (x);
+    psFree (y);
+    psFree (f);
+
+    psFree (poly);
+    return true;
+}
+    
+// interpolate the good pixels to their true centers
+bool psImagePixelInterpolateCenter (psImage *value, psImage *xCoord, psImage *yCoord, psImage *state, psImage *mask, psMaskType maskVal) {
+
+    assert (value->numCols == state->numCols);
+    assert (value->numRows == state->numRows);
+    assert (value->numCols == mask->numCols);
+    assert (value->numRows == mask->numRows);
+
+# if (0)
+    psFits *fits = NULL;
+
+    fits = psFitsOpen ("xcoords.fits", "w");
+    psFitsWriteImage (fits, NULL, xCoord, 0, NULL);
+    psFitsClose (fits);
+
+    fits = psFitsOpen ("ycoords.fits", "w");
+    psFitsWriteImage (fits, NULL, yCoord, 0, NULL);
+    psFitsClose (fits);
+
+    fits = psFitsOpen ("value.fits", "w");
+    psFitsWriteImage (fits, NULL, value, 0, NULL);
+    psFitsClose (fits);
+# endif
+
+    psImage *output = psImageAlloc (value->numCols, value->numRows, PS_TYPE_F32);
+
+    // allocate the vectors for the 2nd order fit below
+    psVector *f  = psVectorAlloc (9, PS_TYPE_F32);
+    // XXX if we add the weight above, include df
+    // psVector *df = psVectorAlloc (9, PS_TYPE_F32); 
+    psVector *x  = psVectorAlloc (9, PS_TYPE_F32);
+    psVector *y  = psVectorAlloc (9, PS_TYPE_F32);
+
+    // allocate a 2D polynomial to fit a quadratic to the valid neighbor pixels.
+    psPolynomial2D *poly2o = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 2, 2);
+    poly2o->mask[2][2] = 1;
+    poly2o->mask[2][1] = 1;
+    poly2o->mask[1][2] = 1;
+
+    // allocate a 2D polynomial to fit a plane to the valid neighbor pixels.
+    psPolynomial2D *poly1o = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 1, 1);
+    poly2o->mask[1][1] = 1;
+
+    for (int iy = 0; iy < state->numRows; iy++) {
+	for (int ix = 0; ix < state->numCols; ix++) {
+
+	    switch (state->data.S32[iy][ix]) {
+	      case PS_IMAGE_INTERPOLATE_GOOD2: {
+		  // XXX is there a fit-image-region function?
+		  int n = 0;
+		  for (int jy = -1; jy <= +1; jy++) {
+		      // skip invalid pixels 
+		      if (jy + iy < 0) { continue; }
+		      if (jy + iy >= value->numRows) { continue; }
+		      for (int jx = -1; jx <= +1; jx++) {
+			  // skip invalid pixels 
+			  if (jx + ix < 0) { continue; } 
+			  if (jx + ix >= value->numCols) { continue; } 
+			  // skip masked pixels
+			  if (mask->data.PS_TYPE_MASK_DATA[iy+jy][ix+jx] & maskVal) { continue; }
+			  x->data.F32[n] = xCoord->data.F32[iy+jy][ix+jx];
+			  y->data.F32[n] = yCoord->data.F32[iy+jy][ix+jx];
+			  f->data.F32[n] = value->data.F32[iy+jy][ix+jx];
+			  // df->data.F32[n] = weight->data.F32[iy+jy][ix+jx];
+			  n++;
+		      }
+		  }
+		  // set vector lengths here
+		  x->n = n;
+		  y->n = n;
+		  f->n = n;
+		  // df->n = n;
+		  // psVectorFitPolynomial2D (poly, NULL, 0xff, f, df, x, y);
+		  psVectorFitPolynomial2D (poly2o, NULL, 0xff, f, NULL, x, y);
+		  // apply the fitted quadratic to get the poor pixel value
+		  // center of pixel is 0.5,0.5
+		  output->data.F32[iy][ix] = psPolynomial2DEval (poly2o, ix + 0.5, iy + 0.5);
+		  break; }
+
+	      case PS_IMAGE_INTERPOLATE_GOOD1: {
+		  // XXX is there a fit-image-region function?
+		  int n = 0;
+		  for (int jy = -1; jy <= +1; jy++) {
+		      // skip invalid pixels 
+		      if (jy + iy < 0) { continue; }
+		      if (jy + iy >= value->numRows) { continue; }
+		      for (int jx = -1; jx <= +1; jx++) {
+			  // skip invalid pixels 
+			  if (jx + ix < 0) { continue; } 
+			  if (jx + ix >= value->numCols) { continue; } 
+			  // skip masked pixels
+			  if (mask->data.PS_TYPE_MASK_DATA[iy+jy][ix+jx] & maskVal) { continue; }
+			  x->data.F32[n] = xCoord->data.F32[iy+jy][ix+jx];
+			  y->data.F32[n] = yCoord->data.F32[iy+jy][ix+jx];
+			  f->data.F32[n] = value->data.F32[iy+jy][ix+jx];
+			  // df->data.F32[n] = weight->data.F32[iy+jy][ix+jx];
+			  n++;
+		      }
+		  }
+		  // set vector lengths here
+		  x->n = n;
+		  y->n = n;
+		  f->n = n;
+		  // df->n = n;
+		  // psVectorFitPolynomial2D (poly1o, NULL, 0xff, f, df, x, y);
+		  psVectorFitPolynomial2D (poly1o, NULL, 0xff, f, NULL, x, y);
+		  // apply the fitted quadratic to get the poor pixel value
+		  output->data.F32[iy][ix] = psPolynomial2DEval (poly1o, ix + 0.5, iy + 0.5);
+		  break; }
+
+	      case PS_IMAGE_INTERPOLATE_GOOD0: {
+		  output->data.F32[iy][ix] = value->data.F32[iy][ix];
+		  break; }
+
+	      default:
+		// skip poor or bad pixels (interpolate later)
+		break;
+	    }
+	}	    
+    }
+
+    for (int iy = 0; iy < value->numRows; iy++) {
+	for (int ix = 0; ix < value->numCols; ix++) {
+	  if (mask->data.PS_TYPE_MASK_DATA[iy][ix] & maskVal) { continue; }
+	  value->data.F32[iy][ix] = output->data.F32[iy][ix];
+	}
+    }
+
+    psFree (x);
+    psFree (y);
+    psFree (f);
+
+    psFree (poly2o);
+    psFree (poly1o);
+
+    psFree (output);
+    return true;
+}
Index: trunk/psLib/src/imageops/psImagePixelInterpolate.h
===================================================================
--- trunk/psLib/src/imageops/psImagePixelInterpolate.h	(revision 14924)
+++ trunk/psLib/src/imageops/psImagePixelInterpolate.h	(revision 14924)
@@ -0,0 +1,41 @@
+/** @file  psImagePixelInterpolate.c
+ *
+ *  @brief Functions for interpolating bad pixels in images
+ *
+ *  @ingroup Image
+ *
+ *  @author Eugene Magnier, IfA
+ *
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-20 23:54:25 $
+ *
+ *  Copyright 2007 Institute for Astronomy, University of Hawaii
+ */
+
+#ifndef PS_IMAGE_PIXEL_INTERPOLATE_H
+#define PS_IMAGE_PIXEL_INTERPOLATE_H
+
+/// @addtogroup ImageOps Image Operations
+/// @{
+
+// XXX make these all bit values?
+typedef enum {
+    PS_IMAGE_INTERPOLATE_GOOD   = 0x10,
+    PS_IMAGE_INTERPOLATE_GOOD0  = 0x11,
+    PS_IMAGE_INTERPOLATE_GOOD1  = 0x12,
+    PS_IMAGE_INTERPOLATE_GOOD2  = 0x13,
+    PS_IMAGE_INTERPOLATE_BAD    = 0x01,
+    PS_IMAGE_INTERPOLATE_CENTER = 0x02,
+    PS_IMAGE_INTERPOLATE_CORNER = 0x04,
+    PS_IMAGE_INTERPOLATE_UR     = 0x04,
+    PS_IMAGE_INTERPOLATE_UL     = 0x05,
+    PS_IMAGE_INTERPOLATE_LR     = 0x06,
+    PS_IMAGE_INTERPOLATE_LL     = 0x07,
+} psImagePixelInterpolateType;
+
+psImage *psImagePixelInterpolateState (int *nBad, int *nPoor, psImage *mask, psMaskType maskVal);
+bool psImagePixelInterpolatePoor (psImage *image, psImage *state, psImage *mask, psMaskType maskVal);
+bool psImagePixelInterpolateCenter (psImage *value, psImage *xCoord, psImage *yCoord, psImage *state, psImage *mask, psMaskType maskVal);
+
+/// @}
+#endif // #ifndef PS_IMAGE_MAP_H