Index: /branches/eam_branch_20070830/psLib/src/imageops/psImagePixelInterpolate.c
===================================================================
--- /branches/eam_branch_20070830/psLib/src/imageops/psImagePixelInterpolate.c	(revision 14723)
+++ /branches/eam_branch_20070830/psLib/src/imageops/psImagePixelInterpolate.c	(revision 14724)
@@ -11,6 +11,6 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-02 02:03:58 $
+ *  @version $Revision: 1.1.2.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-02 02:42:55 $
  *
  *  Copyright 2007 Institute for Astronomy, University of Hawaii
@@ -34,4 +34,5 @@
 #include "psPolynomial.h"
 #include "psPolynomialUtils.h"
+#include "psMinimizePolyFit.h"
 #include "psImage.h"
 #include "psImageInterpolate.h"
@@ -49,6 +50,6 @@
 		    if (jx + ix >= mask->numCols) { continue; } \
 		    /* do not test self */ \
-		    if (!jx && jy) { continue; } \
-		    if (mask->data.PS_TYPE_MASK_DATA[iy][ix] & maskVal) { continue; } \
+		    if (!jx && !jy) { continue; } \
+		    if (mask->data.PS_TYPE_MASK_DATA[iy+jy][ix+jx] & maskVal) { continue; } \
 		    nGood ++; \
 		} \
@@ -115,34 +116,60 @@
 		  // fit a quadratic to the valid neighbor pixels
 		  // apply the quadratic to get the poor pixel value
-		  psPolynomial2D *poly = psImageBicubeFit (image, mask, maskVal, ix, iy);
+		  psPolynomial2D *poly = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 2, 2);
+		  poly->mask[2][2] = 1;
+		  poly->mask[2][1] = 1;
+		  poly->mask[1][2] = 1;
+		  psVector *f = psVectorAlloc (9, PS_TYPE_F32);
+		  psVector *df = psVectorAlloc (9, PS_TYPE_F32);
+		  psVector *x = psVectorAlloc (9, PS_TYPE_F32);
+		  psVector *y = psVectorAlloc (9, PS_TYPE_F32);
+
+		  int n = 0;
+		  for (int jy = -1; jy <= +1; jy++) {
+		      // skip invalid pixels (jx < 0, >= Nx), etc
+		      if (jy + iy < 0) { continue; }
+		      if (jy + iy >= image->numRows) { continue; }
+		      for (int jx = -1; jx <= +1; jx++) {
+			  // skip invalid pixels (jx < 0, >= Nx), etc
+			  if (jx + ix < 0) { continue; } 
+			  if (jx + ix >= mask->numCols) { continue; } 
+			  /* skip self */
+			  if (!jx && !jy) { continue; } 
+			  // skip masked pixels
+			  if (mask->data.PS_TYPE_MASK_DATA[iy+jy][ix+jx] & maskVal) { continue; } \
+			  f->data.F32[n] = image->data.F32[iy+jy][ix+jx];
+			  // df->data.F32[n] = weight->data.F32[jy][jx];
+			  x->data.F32[n] = jx;
+			  y->data.F32[n] = jy;
+			  n++;
+		      }
+		  }
+		  // XXX set vector lengths here
+		  x->n = n;
+		  y->n = n;
+		  f->n = n;
+		  // df->n = n;
+		  psVectorFitPolynomial2D (poly, mask, maskValue, f, df, x, y);
 		  image->data.F32[iy][ix] = poly->coeff[0][0];
 		  break; }
 
 	      case PS_IMAGE_INTERPOLATE_LL: {
-		  // fit a plane to the valid neighbor pixels
-		  // apply the plane to get the poor pixel value
-		  psPolynomial2D *poly = psImageBilinearFit (image, mask, maskVal, ix - 1, iy - 1);
-		  image->data.F32[iy][ix] = poly->coeff[0][0] + poly->coeff[1][0] + poly->coeff[0][1];
+		  // 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];
 		  break; }
 
 	      case PS_IMAGE_INTERPOLATE_LR: {
-		  // fit a plane to the valid neighbor pixels
-		  // apply the plane to get the poor pixel value
-		  psPolynomial2D *poly = psImageBilinearFit (image, mask, maskVal, ix - 0, iy - 1);
-		  image->data.F32[iy][ix] = poly->coeff[0][0] + poly->coeff[0][1];
+		  // 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];
 		  break; }
 
 	      case PS_IMAGE_INTERPOLATE_UL: {
-		  // fit a plane to the valid neighbor pixels
-		  // apply the plane to get the poor pixel value
-		  psPolynomial2D *poly = psImageBilinearFit (image, mask, maskVal, ix - 1, iy - 0);
-		  image->data.F32[iy][ix] = poly->coeff[0][0] + poly->coeff[1][0];
+		  // 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];
 		  break; }
 
 	      case PS_IMAGE_INTERPOLATE_UR: {
-		  // fit a plane to the valid neighbor pixels
-		  // apply the plane to get the poor pixel value
-		  psPolynomial2D *poly = psImageBilinearFit (image, mask, maskVal, ix - 0, iy - 0);
-		  image->data.F32[iy][ix] = poly->coeff[0][0];
+		  // 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];
 		  break; }