Changeset 33341

Timestamp:

Feb 22, 2012, 3:05:56 PM (14 years ago)

Author:

watersc1

Message:

Updated PATTERN.CONTINUITY correction, with better maximum diff limits and constrained fitting to force one cell to have zero offset.

Location:

tags/ipp-20120216/psModules/src/detrend

Files:

• . (modified) (1 prop)
• pmPattern.c (modified) (6 diffs)

tags/ipp-20120216/psModules/src/detrend/pmPattern.c

@@ -680 +680 @@
         j = 8 * (x + 1) + y; // Determine that neighbor's index
         if (fabs(C->data.F64[y][x]) > fabs(D->data.F64[y][x+1])) {
-          critical_value = 2.0 * fabs(D->data.F64[y][x+1] / C->data.F64[y][x]);
+          critical_value = 2.0 * fabs(D->data.F64[y][x+1]);
         }
         else {
-          critical_value = 2.0 * fabs(C->data.F64[y][x] / D->data.F64[y][x+1]);
-        }
-        psTrace("psModules.detrend.cont",5,"CmD %d %d %d %d %g %g", // diagnostic
+          critical_value = 2.0 * fabs(C->data.F64[y][x]);
+        }
+        if (critical_value < 25) { critical_value = 25; }
+        psTrace("psModules.detrend.cont",5,"CmD %d %d %d %d %g %g %g", // diagnostic
                 i,x,y,j,
                 C->data.F64[y][x],
-                D->data.F64[y][x+1]
+                D->data.F64[y][x+1],
+                critical_value
                 );
         if (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_PATTERN_IGNORE)&&  // If there are no errors with the neighbor,
@@ -701 +703 @@
       if (x - 1 > -1) { // etc.
         j = 8 * (x - 1) + y;
-        psTrace("psModules.detrend.cont",5,"DmC %d %d %d %d %g %g",
+        if (fabs(C->data.F64[y][x-1]) > fabs(D->data.F64[y][x])) {
+          critical_value = 2.0 * fabs(D->data.F64[y][x]);
+        }
+        else {
+          critical_value = 2.0 * fabs(C->data.F64[y][x-1]);
+        }
+        if (critical_value < 25) { critical_value = 25; }
+        psTrace("psModules.detrend.cont",5,"DmC %d %d %d %d %g %g %g",
                 i,x,y,j,
                 D->data.F64[y][x],
-                C->data.F64[y][x-1]
+                C->data.F64[y][x-1],
+                critical_value
                 );
-        if (fabs(C->data.F64[y][x-1]) > fabs(D->data.F64[y][x])) {
-          critical_value = 2.0 * fabs(D->data.F64[y][x] / C->data.F64[y][x-1]);
-        }
-        else {
-          critical_value = 2.0 * fabs(C->data.F64[y][x-1] / D->data.F64[y][x]);
-        }
+
         if (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_PATTERN_IGNORE)&&
             (isfinite(D->data.F64[y][x]))&&(isfinite(C->data.F64[y][x-1]))&&
@@ -729 +734 @@
                 );
         if (fabs(A->data.F64[y][x]) > fabs(B->data.F64[y+1][x])) {
-          critical_value = 2.0 * fabs(B->data.F64[y+1][x] / A->data.F64[y][x]);
+          critical_value = 2.0 * fabs(B->data.F64[y+1][x]);
         }
         else {
-          critical_value = 2.0 * fabs(A->data.F64[y][x] / B->data.F64[y+1][x]);
-        }
+          critical_value = 2.0 * fabs(A->data.F64[y][x]);
+        }
+        if (critical_value < 25) { critical_value = 25; }
         if (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_PATTERN_IGNORE)&&
             (isfinite(A->data.F64[y][x]))&&(isfinite(B->data.F64[y+1][x]))&&
@@ -751 +757 @@
                 );
         if (fabs(A->data.F64[y-1][x]) > fabs(B->data.F64[y][x])) {
-          critical_value = 2.0 * fabs(B->data.F64[y][x] / A->data.F64[y-1][x]);
+          critical_value = 2.0 * fabs(B->data.F64[y][x]);
         }
         else {
-          critical_value = 2.0 * fabs(A->data.F64[y-1][x] / B->data.F64[y][x]);
-        }
+          critical_value = 2.0 * fabs(A->data.F64[y-1][x]);
+        }
+        if (critical_value < 25) { critical_value = 25; }
         if (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_PATTERN_IGNORE)&&
             (isfinite(B->data.F64[y][x]))&&(isfinite(A->data.F64[y-1][x]))&&
@@ -766 +773 @@
       }
     }
+    double max_XX = 0;
+    double solution_V = 0;
+    int i_peak = -1;
     for (int i = 0; i < numCells; i++) { // If any cells have no value of themself, set the matrix to 1.0.
       if (XX->data.F64[i][i] == 0.0) {
         XX->data.F64[i][i] = 1.0;
       }
-    }
-#if (0)
+      if (XX->data.F64[i][i] > max_XX) {
+        max_XX = XX->data.F64[i][i];
+        solution_V = solution->data.F64[i];
+        i_peak = i;
+      }
+    }
+    psTrace("psModules.detrend.cont",5,"fixed point: %d %g\n",
+            i_peak,solution_V);
+
+    for (int i = 0; i < numCells; i++) {
+/*        if (!((XX->data.F64[i][i] == 1.0)&& */
+/*          (solution->data.F64[i] == 0.0))) { */
+        solution->data.F64[i] -= solution_V;
+        if (i != i_peak) {
+          for (int j = 0; j < numCells; j++) {
+            XX->data.F64[i][j] -= XX->data.F64[i_peak][j];
+          }
+        }
+/*        } */
+    }
+    for (int i = 0; i < numCells; i++) {
+      XX->data.F64[i_peak][i] = 0.0;
+    }
+    XX->data.F64[i_peak][i_peak] = 1.0;
+    
+    
+#if (1)
     for (int i = 0; i < numCells; i++) { // print matrix A
       psTrace("psModules.detrend.cont",5,"A: %3d % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f",
@@ -805 +840 @@
     //    psMatrixLUSolve(XX,solution);
     psMatrixGJSolve(XX,solution);
-#if (0)
+#if (1)
     for (int i = 0; i < numCells; i++) { // print vector b
-      psTrace("psModules.detrend.cont",5,"b: %d %f",
+      psTrace("psModules.detrend.cont",5,"x: %d %f",
               i,
               solution->data.F64[i]