Index: trunk/psLib/src/imageops/psImageMapFit.c
===================================================================
--- trunk/psLib/src/imageops/psImageMapFit.c	(revision 14960)
+++ 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;
             }