Changeset 14864 for branches/eam_branch_20070830/psLib/src/imageops/psImageMapFit.c

Timestamp:

Sep 17, 2007, 10:17:21 AM (19 years ago)

Author:

magnier

Message:

added errors to psImageMap, NaN for unconstrained pixels in output map

File:

• branches/eam_branch_20070830/psLib/src/imageops/psImageMapFit.c (modified) (5 diffs)

branches/eam_branch_20070830/psLib/src/imageops/psImageMapFit.c

@@ -7 +7 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-17 01:14:52 $
+ *  @version $Revision: 1.1.2.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-17 20:17:21 $
  *
  *  Copyright 2007 Institute for Astronomy, University of Hawaii
@@ -22 +22 @@
 
 // XXX for testing
-// #include "psFits.h"
-// #include "psFitsImage.h"
+#include "psFits.h"
+#include "psFitsImage.h"
 
 #include "psAssert.h"
@@ -238 +238 @@
     // fprintf (stderr, "Total: %f\n", Total);
 
-    // write 1s in unused diagonal elements
-    // XXX for test only
-    # if (0)
-    for (int n = 0; n < Nx; n+=Nx-1) {
-        for (int m = 0; m < Ny; m++) {
-            I = n + Nx * m;
-            A->data.F32[I][I] = 1.0;
-        }
-    }
-    for (int n = 0; n < Nx; n++) {
-        for (int m = 0; m < Ny; m+=Ny-1) {
-            I = n + Nx * m;
-            A->data.F32[I][I] = 1.0;
-        }
-    }
-    # endif
-
-    // psFits *fits = psFitsOpen ("Agj.fits", "w");
-    // psFitsWriteImage (fits, NULL, A, 0, NULL);
-    // psFitsClose (fits);
+    // test for empty diagonal elements (unconstained cells), mark, and set pivots to 1.0
+    psVector *Empty = psVectorAlloc (Nx*Ny, PS_TYPE_S8);
+    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;
+        }
+    }
+
+    psFits *fits = psFitsOpen ("Agj.fits", "w");
+    psFitsWriteImage (fits, NULL, A, 0, NULL);
+    psFitsClose (fits);
+
+    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];
+    }
+
+    fits = psFitsOpen ("Bgj.fits", "w");
+    psFitsWriteImage (fits, NULL, vector, 0, NULL);
+    psFitsClose (fits);
+    psFree (vector);
 
     if (!psMatrixGJSolveF32(A, B)) {
@@ -267 +275 @@
     }
     
+    // set bad values to NaN
+    for (int i = 0; i < Nx*Ny; i++) {
+        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]);
         }
     }
@@ -276 +293 @@
     psFree (A);
     psFree (B);
+    psFree (Empty);
 
     return true;