Index: trunk/psLib/src/imageops/psImageMapFit.c
===================================================================
--- trunk/psLib/src/imageops/psImageMapFit.c	(revision 15001)
+++ trunk/psLib/src/imageops/psImageMapFit.c	(revision 15041)
@@ -7,6 +7,6 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-24 21:29:10 $
+ *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-27 04:27:03 $
  *
  *  Copyright 2007 Institute for Astronomy, University of Hawaii
@@ -68,10 +68,12 @@
 
     if (Nx == 1) {
-        psAbort ("un-implemented edge case");
-        goto insert;
+	bool status;
+	status = psImageMapFit1DinY (map, mask, maskValue, x, y, f, df);
+	return status;
     }
     if (Ny == 1) {
-        psAbort ("un-implemented edge case");
-        goto insert;
+	bool status;
+	status = psImageMapFit1DinX (map, mask, maskValue, x, y, f, df);
+	return status;
     }
 
@@ -231,5 +233,4 @@
             sA[+1][+1] = dx_rx_dy_ry;
 
-        insert:
             // I[ 0][ 0] = index for this n,m element:
             I = n + Nx * m;
@@ -271,5 +272,5 @@
     }
 
-    # if (0)
+# if (0)
     psFits *fits = psFitsOpen ("Agj.fits", "w");
     psFitsWriteImage (fits, NULL, A, 0, NULL);
@@ -285,5 +286,5 @@
     psFitsClose (fits);
     psFree (vector);
-    # endif
+# endif
 
     if (!psMatrixGJSolveF32(A, B)) {
@@ -433,2 +434,299 @@
     return true;
 }
+
+// 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;
+
+    // XXX Add Asserts
+    assert (map->binning->nXruff == 1);
+
+    // dimensions of the output map image
+    int Ny = map->binning->nYruff;
+
+    // set up the redirection table so we can use sA[-1][-1], etc
+    float SAv[3], *sA;
+
+    sA = SAv + 1;
+
+    // elements of the matrix equation Ax = B; we are solving for the vector x
+    psImage *A = psImageAlloc (Ny, Ny, PS_TYPE_F32);
+    psVector *B = psVectorAlloc (Ny, PS_TYPE_F32);
+
+    psImageInit (A, 0.0);
+    psVectorInit (B, 0.0);
+
+    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];
+	}
+    }
+
+    // test for empty diagonal elements (unconstained cells), mark, and set pivots to 1.0
+    psVector *Empty = psVectorAlloc (Ny, PS_TYPE_S8);
+    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 (!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;
+    }
+
+    // set bad values to NaN
+    for (int i = 0; i < Ny; i++) {
+	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]);
+    }
+
+    psFree (A);
+    psFree (B);
+    psFree (Empty);
+
+    return true;
+}
+
+// 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;
+
+    // XXX Add Asserts
+    assert (map->binning->nXruff == 1);
+
+    // dimensions of the output map image
+    int Nx = map->binning->nXruff;
+
+    // set up the redirection table so we can use sA[-1][-1], etc
+    float SAv[3], *sA;
+
+    sA = SAv + 1;
+
+    // elements of the matrix equation Ax = B; we are solving for the vector x
+    psImage *A = psImageAlloc (Nx, Nx, PS_TYPE_F32);
+    psVector *B = psVectorAlloc (Nx, PS_TYPE_F32);
+
+    psImageInit (A, 0.0);
+    psVectorInit (B, 0.0);
+
+    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];
+	}
+    }
+
+    // test for empty diagonal elements (unconstained cells), mark, and set pivots to 1.0
+    psVector *Empty = psVectorAlloc (Nx, PS_TYPE_S8);
+    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 (!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;
+    }
+
+    // set bad values to NaN
+    for (int i = 0; i < Nx; i++) {
+	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]);
+    }
+
+    psFree (A);
+    psFree (B);
+    psFree (Empty);
+
+    return true;
+}
+