Index: /branches/eam_branch_20070830/psLib/src/imageops/psImagePixelInterpolate.c
===================================================================
--- /branches/eam_branch_20070830/psLib/src/imageops/psImagePixelInterpolate.c	(revision 14775)
+++ /branches/eam_branch_20070830/psLib/src/imageops/psImagePixelInterpolate.c	(revision 14776)
@@ -11,6 +11,6 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.1.2.5 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-03 20:29:07 $
+ *  @version $Revision: 1.1.2.6 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-07 20:21:51 $
  *
  *  Copyright 2007 Institute for Astronomy, University of Hawaii
@@ -38,4 +38,7 @@
 #include "psImageInterpolate.h"
 #include "psImagePixelInterpolate.h"
+
+#include "psFits.h"
+#include "psFitsImage.h"
 
 # define PS_IMAGE_ITER_STATE(XS,XE,YS,YE,N_MIN,TYPE) \
@@ -75,6 +78,42 @@
 	for (int ix = 0; ix < mask->numCols; ix++) {
 
-	    // skip the good pixels (unmasked)
-	    if (!(mask->data.PS_TYPE_MASK_DATA[iy][ix] & maskVal)) { continue; }
+	    // 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
@@ -85,8 +124,8 @@
 	    // 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_LL);
-	    PS_IMAGE_ITER_STATE (-1,+0,+0,+1,3, PS_IMAGE_INTERPOLATE_UL);
-	    PS_IMAGE_ITER_STATE (+0,+1,-1,+0,3, PS_IMAGE_INTERPOLATE_LR);
-	    PS_IMAGE_ITER_STATE (+0,+1,+0,+1,3, PS_IMAGE_INTERPOLATE_UR);
+	    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 ++;
@@ -122,5 +161,7 @@
 
 	    switch (state->data.S32[iy][ix]) {
-	      case PS_IMAGE_INTERPOLATE_GOOD: 
+	      case PS_IMAGE_INTERPOLATE_GOOD0: 
+	      case PS_IMAGE_INTERPOLATE_GOOD1: 
+	      case PS_IMAGE_INTERPOLATE_GOOD2: 
 		// skip the good pixels
 		break; 
@@ -167,20 +208,20 @@
 	      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+0] - image->data.F32[iy+0][ix+1] - image->data.F32[iy+1][ix+1];
+		  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+0] - image->data.F32[iy+0][ix-1] - image->data.F32[iy+1][ix-1];
+		  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+0] - image->data.F32[iy+0][ix+1] - image->data.F32[iy-1][ix+1];
+		  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+0] - image->data.F32[iy+0][ix-1] - image->data.F32[iy-1][ix-1];
+		  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; }
 
@@ -207,4 +248,20 @@
     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);
 
@@ -217,12 +274,12 @@
 
     // allocate a 2D polynomial to fit a quadratic to the valid neighbor pixels.
-    psPolynomial2D *poly2D = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 2, 2);
-    poly2D->mask[2][2] = 1;
-    poly2D->mask[2][1] = 1;
-    poly2D->mask[1][2] = 1;
+    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 *poly1D = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 1, 1);
-    poly2D->mask[1][1] = 1;
+    psPolynomial2D *poly1o = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 1, 1);
+    poly2o->mask[1][1] = 1;
 
     for (int iy = 0; iy < state->numRows; iy++) {
@@ -256,7 +313,8 @@
 		  // df->n = n;
 		  // psVectorFitPolynomial2D (poly, NULL, 0xff, f, df, x, y);
-		  psVectorFitPolynomial2D (poly2D, NULL, 0xff, f, NULL, x, y);
+		  psVectorFitPolynomial2D (poly2o, NULL, 0xff, f, NULL, x, y);
 		  // apply the fitted quadratic to get the poor pixel value
-		  output->data.F32[iy][ix] = psPolynomial2DEval (poly2D, ix, iy);
+		  // center of pixel is 0.5,0.5
+		  output->data.F32[iy][ix] = psPolynomial2DEval (poly2o, ix + 0.5, iy + 0.5);
 		  break; }
 
@@ -286,16 +344,17 @@
 		  f->n = n;
 		  // df->n = n;
-		  // psVectorFitPolynomial2D (poly1D, NULL, 0xff, f, df, x, y);
-		  psVectorFitPolynomial2D (poly1D, NULL, 0xff, f, NULL, x, y);
+		  // 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 (poly1D, ix, iy);
+		  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+jy][ix+jx];
+		  output->data.F32[iy][ix] = value->data.F32[iy][ix];
 		  break; }
 
 	      default:
-		psAbort("impossible case in __func__");
+		// skip poor or bad pixels (interpolate later)
+		break;
 	    }
 	}	    
@@ -313,6 +372,6 @@
     psFree (f);
 
-    psFree (poly2D);
-    psFree (poly1D);
+    psFree (poly2o);
+    psFree (poly1o);
 
     psFree (output);
Index: /branches/eam_branch_20070830/psLib/src/imageops/psImagePixelInterpolate.h
===================================================================
--- /branches/eam_branch_20070830/psLib/src/imageops/psImagePixelInterpolate.h	(revision 14775)
+++ /branches/eam_branch_20070830/psLib/src/imageops/psImagePixelInterpolate.h	(revision 14776)
@@ -7,6 +7,6 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-02 20:29:50 $
+ *  @version $Revision: 1.1.2.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-07 20:21:53 $
  *
  *  Copyright 2007 Institute for Astronomy, University of Hawaii
@@ -21,5 +21,8 @@
 // XXX make these all bit values?
 typedef enum {
-    PS_IMAGE_INTERPOLATE_GOOD   = 0x00,
+    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,
@@ -33,4 +36,5 @@
 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);
 
 /// @}