Index: trunk/psLib/src/imageops/psImageMapFit.c
===================================================================
--- trunk/psLib/src/imageops/psImageMapFit.c	(revision 15041)
+++ trunk/psLib/src/imageops/psImageMapFit.c	(revision 15841)
@@ -7,6 +7,6 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-27 04:27:03 $
+ *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-12-15 01:20:03 $
  *
  *  Copyright 2007 Institute for Astronomy, University of Hawaii
@@ -37,4 +37,7 @@
 // #include "psImageUnbin.h"
 
+#include "psImageMapFit.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
@@ -44,8 +47,7 @@
 
 // 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;
-
+bool psImageMapFit(psImageMap *map, const psVector *mask, psMaskType maskValue,
+                   const psVector *x, const psVector *y, const psVector *f, const psVector *df)
+{
     // XXX Add Asserts
 
@@ -56,27 +58,39 @@
     // no spatial information, just calculate mean & stdev
     if ((Nx == 1) && (Ny == 1)) {
-        psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+        psStatsInit(map->stats);
+
+        // the user has supplied one of various stats option pairs,
+        psStatsOptions mean = psStatsMeanOption(map->stats->options);
+        psStatsOptions stdev = psStatsStdevOption(map->stats->options);
+        if (!psStatsSingleOption(mean)) {
+            psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected");
+            return false;
+        }
+        if (!psStatsSingleOption(stdev)) {
+            psError(PS_ERR_UNKNOWN, true, "no valid stdev stats option selected");
+            return false;
+        }
 
         // 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);
+        psVectorStats(map->stats, f, NULL, mask, maskValue);
+
+        map->map->data.F32[0][0]   = psStatsGetValue(map->stats, mean);
+        map->error->data.F32[0][0] = psStatsGetValue(map->stats, stdev);
         return true;
     }
 
     if (Nx == 1) {
-	bool status;
-	status = psImageMapFit1DinY (map, mask, maskValue, x, y, f, df);
-	return status;
+        bool status;
+        status = psImageMapFit1DinY (map, mask, maskValue, x, y, f, df);
+        return status;
     }
     if (Ny == 1) {
-	bool status;
-	status = psImageMapFit1DinX (map, mask, maskValue, x, y, f, df);
-	return status;
+        bool status;
+        status = psImageMapFit1DinX (map, mask, maskValue, x, y, f, df);
+        return status;
     }
 
     // set up the redirection table so we can use sA[-1][-1], etc
+    // XXX psKernel does this for you --- PAP.
     float SAm[3][3], *SAv[3], **sA;
     // float TAm[3][3], *TAv[3], **tA;
@@ -182,8 +196,8 @@
                 if (m == Ny - 1) Qy = 0;
 
-		assert (isfinite(fi));
-		assert (isfinite(wt));
-		assert (isfinite(rx));
-		assert (isfinite(ry));
+                assert (isfinite(fi));
+                assert (isfinite(wt));
+                assert (isfinite(rx));
+                assert (isfinite(ry));
 
                 // points at offset 1,1
@@ -234,5 +248,5 @@
 
             // I[ 0][ 0] = index for this n,m element:
-            I = n + Nx * m;
+            int I = n + Nx * m;
             B->data.F32[I] = fi_rx_ry + fi_rx_py + fi_px_ry + fi_px_py;
 
@@ -245,5 +259,5 @@
                     if (m + jm <   0) continue;
                     if (m + jm >= Ny) continue;
-                    J = (n + jn) + Nx * (m + jm);
+                    int 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]);
@@ -306,5 +320,5 @@
     for (int n = 0; n < Nx; n++) {
         for (int m = 0; m < Ny; m++) {
-            I = n + Nx * m;
+            int 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]);
@@ -320,22 +334,23 @@
 
 // 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) {
-
+bool psImageMapClipFit(psImageMap *map, psStats *stats, psVector *inMask, psMaskType maskValue,
+                       const psVector *x, const psVector *y, const psVector *f, const psVector *df)
+{
     // XXX add in full PS_ASSERTS
-    assert (map);
-    assert (stats);
-    assert (x);
-    assert (y);
-    assert (f);
+    assert(map);
+    assert(stats);
+    assert(x);
+    assert(y);
+    assert(f);
 
     // 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) {
+    psStatsOptions meanOption = psStatsMeanOption(stats->options);
+    psStatsOptions stdevOption = psStatsStdevOption(stats->options);
+    if (!psStatsSingleOption(meanOption)) {
         psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected");
         return false;
     }
-    if (!stdevOption) {
+    if (!psStatsSingleOption(stdevOption)) {
         psError(PS_ERR_UNKNOWN, true, "no valid stdev stats option selected");
         return false;
@@ -436,8 +451,7 @@
 
 // map defines the output image dimensions and scaling.
-bool psImageMapFit1DinY (psImageMap *map, psVector *mask, psMaskType maskValue, psVector *x, psVector *y, psVector *f, psVector *df) {
-
-    int I, J;
-
+bool psImageMapFit1DinY(psImageMap *map, const psVector *mask, psMaskType maskValue,
+                        const psVector *x, const psVector *y, const psVector *f, const psVector *df)
+{
     // XXX Add Asserts
     assert (map->binning->nXruff == 1);
@@ -459,84 +473,84 @@
 
     for (int m = 0; m < Ny; m++) {
-	// define & init summing variables
-	float ry_ry = 0;
-	float dy_ry = 0;
-	float fi_ry = 0;
-	float py_py = 0;
-	float qy_py = 0;
-	float fi_py = 0;
-
-	// generate the sums for the fitting matrix element I,J
-	// I = m
-	// J = m + jm
-	for (int i = 0; i < y->n; i++) {
-
-	    if (mask && (mask->data.U8[i] & maskValue)) continue;
-
-	    float dy = psImageBinningGetRuffY (map->binning, y->data.F32[i]) - (m + 0.5);
-
-	    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 (!edgeY && (fabs(dy) > 1.0)) continue;
-
-	    // related offset values
-	    float ry = 1.0 - dy;
-	    float py = 1.0 + dy;
-	    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 Qy = (dy >= 0) ? 1 : 0;
-	    if (m ==      0) Qy = 1;
-	    if (m == Ny - 1) Qy = 0;
-
-	    assert (isfinite(fi));
-	    assert (isfinite(wt));
-	    assert (isfinite(ry));
-
-	    // points at offset 1,1
-	    if (Qy == 1) {
-		ry_ry += ry*ry*wt;
-		dy_ry += dy*ry*wt;
-		fi_ry += fi*ry*wt;
-	    }
-	    // points at offset 1,0
-	    if (Qy == 0) {
-		py_py += py*py*wt;
-		qy_py += qy*py*wt;
-		fi_py += fi*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] = qy_py;
-	sA[ 0] = ry_ry + py_py;
-	sA[+1] = dy_ry;
-
-	// I[ 0][ 0] = index for this n,m element:
-	I = m;
-	B->data.F32[I] = fi_ry + fi_py;
-
-	// insert these values into their corresponding locations in A, B
-	for (int jm = -1; jm <= +1; jm++) {
-	    if (m + jm <   0) continue;
-	    if (m + jm >= Ny) continue;
-	    J = (m + jm);
-	    A->data.F32[J][I] = sA[jm];
-	}
+        // define & init summing variables
+        float ry_ry = 0;
+        float dy_ry = 0;
+        float fi_ry = 0;
+        float py_py = 0;
+        float qy_py = 0;
+        float fi_py = 0;
+
+        // generate the sums for the fitting matrix element I,J
+        // I = m
+        // J = m + jm
+        for (int i = 0; i < y->n; i++) {
+
+            if (mask && (mask->data.U8[i] & maskValue)) continue;
+
+            float dy = psImageBinningGetRuffY (map->binning, y->data.F32[i]) - (m + 0.5);
+
+            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 (!edgeY && (fabs(dy) > 1.0)) continue;
+
+            // related offset values
+            float ry = 1.0 - dy;
+            float py = 1.0 + dy;
+            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 Qy = (dy >= 0) ? 1 : 0;
+            if (m ==      0) Qy = 1;
+            if (m == Ny - 1) Qy = 0;
+
+            assert (isfinite(fi));
+            assert (isfinite(wt));
+            assert (isfinite(ry));
+
+            // points at offset 1,1
+            if (Qy == 1) {
+                ry_ry += ry*ry*wt;
+                dy_ry += dy*ry*wt;
+                fi_ry += fi*ry*wt;
+            }
+            // points at offset 1,0
+            if (Qy == 0) {
+                py_py += py*py*wt;
+                qy_py += qy*py*wt;
+                fi_py += fi*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] = qy_py;
+        sA[ 0] = ry_ry + py_py;
+        sA[+1] = dy_ry;
+
+        // I[ 0][ 0] = index for this n,m element:
+        int I = m;
+        B->data.F32[I] = fi_ry + fi_py;
+
+        // insert these values into their corresponding locations in A, B
+        for (int jm = -1; jm <= +1; jm++) {
+            if (m + jm <   0) continue;
+            if (m + jm >= Ny) continue;
+            int J = (m + jm);
+            A->data.F32[J][I] = sA[jm];
+        }
     }
 
@@ -545,33 +559,33 @@
     psVectorInit (Empty, 0);
     for (int i = 0; i < Ny; i++) {
-	if (A->data.F32[i][i] == 0.0) {
-	    Empty->data.S8[i] = 1;
-	    for (int j = 0; j < 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 < 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 (!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 < Ny; i++) {
-	if (Empty->data.S8[i]) {
-	    B->data.F32[i] = NAN;
-	}
+        if (Empty->data.S8[i]) {
+            B->data.F32[i] = NAN;
+        }
     }
 
     for (int m = 0; m < Ny; m++) {
-	map->map->data.F32[m][0] = B->data.F32[m];
-	map->error->data.F32[m][0] = sqrt(A->data.F32[m][m]);
+        map->map->data.F32[m][0] = B->data.F32[m];
+        map->error->data.F32[m][0] = sqrt(A->data.F32[m][m]);
     }
 
@@ -584,8 +598,7 @@
 
 // map defines the output image dimensions and scaling.
-bool psImageMapFit1DinX (psImageMap *map, psVector *mask, psMaskType maskValue, psVector *x, psVector *y, psVector *f, psVector *df) {
-
-    int I, J;
-
+bool psImageMapFit1DinX(psImageMap *map, const psVector *mask, psMaskType maskValue,
+                        const psVector *x, const psVector *y, const psVector *f, const psVector *df)
+{
     // XXX Add Asserts
     assert (map->binning->nXruff == 1);
@@ -607,84 +620,84 @@
 
     for (int m = 0; m < Nx; m++) {
-	// define & init summing variables
-	float rx_rx = 0;
-	float dx_rx = 0;
-	float fi_rx = 0;
-	float px_px = 0;
-	float qx_px = 0;
-	float fi_px = 0;
-
-	// generate the sums for the fitting matrix element I,J
-	// I = m
-	// J = m + jm
-	for (int i = 0; i < x->n; i++) {
-
-	    if (mask && (mask->data.U8[i] & maskValue)) continue;
-
-	    float dx = psImageBinningGetRuffX (map->binning, x->data.F32[i]) - (m + 0.5);
-
-	    bool edgeX = false;
-	    edgeX |= ((m == 1) && (dx < -1.0));
-	    edgeX |= ((m == Nx - 2) && (dx > +1.0));
-
-	    // skip points outside of 2x2 grid centered on n,m:
-	    if (!edgeX && (fabs(dx) > 1.0)) continue;
-
-	    // related offset values
-	    float rx = 1.0 - dx;
-	    float px = 1.0 + dx;
-	    float qx = -dx;
-
-	    // 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 (m ==      0) Qx = 1;
-	    if (m == Nx - 1) Qx = 0;
-
-	    assert (isfinite(fi));
-	    assert (isfinite(wt));
-	    assert (isfinite(rx));
-
-	    // points at offset 1,1
-	    if (Qx == 1) {
-		rx_rx += rx*rx*wt;
-		dx_rx += dx*rx*wt;
-		fi_rx += fi*rx*wt;
-	    }
-	    // points at offset 1,0
-	    if (Qx == 0) {
-		px_px += px*px*wt;
-		qx_px += qx*px*wt;
-		fi_px += fi*px*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] = qx_px;
-	sA[ 0] = rx_rx + px_px;
-	sA[+1] = dx_rx;
-
-	// I[ 0][ 0] = index for this n,m element:
-	I = m;
-	B->data.F32[I] = fi_rx + fi_px;
-
-	// insert these values into their corresponding locations in A, B
-	for (int jm = -1; jm <= +1; jm++) {
-	    if (m + jm <   0) continue;
-	    if (m + jm >= Nx) continue;
-	    J = (m + jm);
-	    A->data.F32[J][I] = sA[jm];
-	}
+        // define & init summing variables
+        float rx_rx = 0;
+        float dx_rx = 0;
+        float fi_rx = 0;
+        float px_px = 0;
+        float qx_px = 0;
+        float fi_px = 0;
+
+        // generate the sums for the fitting matrix element I,J
+        // I = m
+        // J = m + jm
+        for (int i = 0; i < x->n; i++) {
+
+            if (mask && (mask->data.U8[i] & maskValue)) continue;
+
+            float dx = psImageBinningGetRuffX (map->binning, x->data.F32[i]) - (m + 0.5);
+
+            bool edgeX = false;
+            edgeX |= ((m == 1) && (dx < -1.0));
+            edgeX |= ((m == Nx - 2) && (dx > +1.0));
+
+            // skip points outside of 2x2 grid centered on n,m:
+            if (!edgeX && (fabs(dx) > 1.0)) continue;
+
+            // related offset values
+            float rx = 1.0 - dx;
+            float px = 1.0 + dx;
+            float qx = -dx;
+
+            // 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 (m ==      0) Qx = 1;
+            if (m == Nx - 1) Qx = 0;
+
+            assert (isfinite(fi));
+            assert (isfinite(wt));
+            assert (isfinite(rx));
+
+            // points at offset 1,1
+            if (Qx == 1) {
+                rx_rx += rx*rx*wt;
+                dx_rx += dx*rx*wt;
+                fi_rx += fi*rx*wt;
+            }
+            // points at offset 1,0
+            if (Qx == 0) {
+                px_px += px*px*wt;
+                qx_px += qx*px*wt;
+                fi_px += fi*px*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] = qx_px;
+        sA[ 0] = rx_rx + px_px;
+        sA[+1] = dx_rx;
+
+        // I[ 0][ 0] = index for this n,m element:
+        int I = m;
+        B->data.F32[I] = fi_rx + fi_px;
+
+        // insert these values into their corresponding locations in A, B
+        for (int jm = -1; jm <= +1; jm++) {
+            if (m + jm <   0) continue;
+            if (m + jm >= Nx) continue;
+            int J = (m + jm);
+            A->data.F32[J][I] = sA[jm];
+        }
     }
 
@@ -693,33 +706,34 @@
     psVectorInit (Empty, 0);
     for (int i = 0; i < Nx; i++) {
-	if (A->data.F32[i][i] == 0.0) {
-	    Empty->data.S8[i] = 1;
-	    for (int j = 0; j < Nx; 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; 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 (!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; i++) {
-	if (Empty->data.S8[i]) {
-	    B->data.F32[i] = NAN;
-	}
+        if (Empty->data.S8[i]) {
+            B->data.F32[i] = NAN;
+        }
     }
 
     for (int m = 0; m < Nx; m++) {
-	map->map->data.F32[0][m] = B->data.F32[I];
-	map->error->data.F32[0][m] = sqrt(A->data.F32[I][I]);
+        int I = m; // XXX I'm not entirely sure about this; it wasn't set for this scope --- PAP.
+        map->map->data.F32[0][m] = B->data.F32[I];
+        map->error->data.F32[0][m] = sqrt(A->data.F32[I][I]);
     }