Changeset 15041 for trunk/psLib/src/imageops/psImageMapFit.c

Timestamp:

Sep 26, 2007, 6:27:03 PM (19 years ago)

Author:

eugene

Message:

added case of 1xN or Nx1 maps

File:

• trunk/psLib/src/imageops/psImageMapFit.c (modified) (6 diffs)

trunk/psLib/src/imageops/psImageMapFit.c

@@ -7 +7 @@
  *  @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 +68 @@
 
     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 +233 @@
             sA[+1][+1] = dx_rx_dy_ry;
 
-        insert:
             // I[ 0][ 0] = index for this n,m element:
             I = n + Nx * m;
@@ -271 +272 @@
     }
 
-    # if (0)
+# if (0)
     psFits *fits = psFitsOpen ("Agj.fits", "w");
     psFitsWriteImage (fits, NULL, A, 0, NULL);
@@ -285 +286 @@
     psFitsClose (fits);
     psFree (vector);
-    # endif
+# endif
 
     if (!psMatrixGJSolveF32(A, B)) {
@@ -433 +434 @@
     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;
+}
+