Index: /branches/eam_branch_20070830/psLib/src/imageops/psImageMapFit.c
===================================================================
--- /branches/eam_branch_20070830/psLib/src/imageops/psImageMapFit.c	(revision 14861)
+++ /branches/eam_branch_20070830/psLib/src/imageops/psImageMapFit.c	(revision 14862)
@@ -7,6 +7,6 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-15 19:43:43 $
+ *  @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-17 01:14:52 $
  *
  *  Copyright 2007 Institute for Astronomy, University of Hawaii
@@ -22,6 +22,6 @@
 
 // XXX for testing
-#include "psFits.h"
-#include "psFitsImage.h"
+// #include "psFits.h"
+// #include "psFitsImage.h"
 
 #include "psAssert.h"
@@ -54,16 +54,19 @@
     // 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;
+    // float TAm[3][3], *TAv[3], **tA;
 
     for (int i = 0; i < 3; i++) {
 	SAv[i] = SAm[i] + 1;
-	TAv[i] = TAm[i] + 1;
+	// TAv[i] = TAm[i] + 1;
     }
     sA = SAv + 1;
-    tA = TAv + 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; 
@@ -72,5 +75,5 @@
     // for (int m = 1; m < Ny - 1; m++) {
     
-    float Total = 0.0;
+    // float Total = 0.0;
     for (int n = 0; n < Nx; n++) {
 	for (int m = 0; m < Ny; m++) {
@@ -108,7 +111,16 @@
 		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 (fabs(dx) > 1.0) continue;
-		if (fabs(dy) > 1.0) continue;
+		if (!edgeX && (fabs(dx) > 1.0)) continue;
+		if (!edgeY && (fabs(dy) > 1.0)) continue;
 
 		// related offset values
@@ -133,6 +145,14 @@
 		// 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 ((dx >= 0) && (dy >= 0)) {
+		if ((Qx == 1) && (Qy == 1)) {
 		    rx_rx_ry_ry += rx*rx*ry*ry*wt;
 		    rx_rx_dy_ry += rx*rx*dy*ry*wt;
@@ -142,5 +162,5 @@
 		}
 		// points at offset 1,0
-		if ((dx >= 0) && (dy <  0)) {
+		if ((Qx == 1) && (Qy == 0)) {
 		    rx_rx_py_py += rx*rx*py*py*wt;
 		    rx_rx_qy_py += rx*rx*qy*py*wt;
@@ -150,5 +170,5 @@
 		}
 		// points at offset 0,1
-		if ((dx <  0) && (dy >= 0)) {
+		if ((Qx == 0) && (Qy == 1)) {
 		    px_px_ry_ry += px*px*ry*ry*wt;
 		    px_px_dy_ry += px*px*dy*ry*wt;
@@ -158,5 +178,5 @@
 		}
 		// points at offset 0,0
-		if ((dx <  0) && (dy <  0)) {
+		if ((Qx == 0) && (Qy == 0)) {
 		    px_px_py_py += px*px*py*py*wt;
 		    px_px_qy_py += px*px*qy*py*wt;
@@ -178,18 +198,4 @@
 	    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:
@@ -207,85 +213,4 @@
 	    }
 
-	    // 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:
@@ -294,5 +219,5 @@
 
 	    // insert these values into their corresponding locations in A, B
-	    float Sum = 0.0;
+	    // float Sum = 0.0;
 	    for (int jn = -1; jn <= +1; jn++) {
 		if (n + jn <   0) continue;
@@ -302,15 +227,14 @@
 		    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);
+		    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);
 
     // write 1s in unused diagonal elements
@@ -331,7 +255,7 @@
     # endif
 
-    psFits *fits = psFitsOpen ("Agj.fits", "w");
-    psFitsWriteImage (fits, NULL, A, 0, NULL);
-    psFitsClose (fits);
+    // psFits *fits = psFitsOpen ("Agj.fits", "w");
+    // psFitsWriteImage (fits, NULL, A, 0, NULL);
+    // psFitsClose (fits);
 
     if (!psMatrixGJSolveF32(A, B)) {