Changeset 33304 for trunk/psModules/src/detrend/pmPattern.c

Timestamp:

Feb 17, 2012, 4:44:49 PM (14 years ago)

Author:

watersc1

Message:

Update to continuity code that sets a limit on how large a jump will be allowed into the continuity matrix. This means we will leave some cells uncorrected, which may cause issues around them, but will keep the rest of the cells in a single (hopefully valid) solution.

File:

• trunk/psModules/src/detrend/pmPattern.c (modified) (10 diffs)

trunk/psModules/src/detrend/pmPattern.c

@@ -565 +565 @@
     psImage *C = psImageAlloc(8,8,PS_TYPE_F64); // Right edge
     psImage *D = psImageAlloc(8,8,PS_TYPE_F64); // Left edge
+    psImageInit(A,0.0);
+    psImageInit(B,0.0);
+    psImageInit(C,0.0);
+    psImageInit(D,0.0);
     
     for (int i = 0; i < numCells; i++) {
@@ -611 +615 @@
             continue; // Move on to next edge, as only part of this cell may be a problem
           }
-
+ 
           // If the returned value is zero, assume something is wrong.  Do I still need this?
           if (psStatsGetValue(bgStats,bgStat) < 1e-6) {
@@ -632 +636 @@
             else if (j == 3) { D->data.F64[y][x] = psStatsGetValue(bgStats,bgStat); }
           }
+
+          for (int u = 0; u < subset->numCols; u++) {
+            for (int v = 0; v < subset->numRows; v++) {
+              psTrace("psModules.detrend.cont",10,"BOX: %d %d (%d %d) (%d %d) %f %d",
+                      i,j,x,y,u,v,subset->data.F32[v][u],submask->data.PS_TYPE_IMAGE_MASK_DATA[v][u]);
+            }
+          }       
+          
           psFree(subset);
           psFree(submask);
@@ -638 +650 @@
         psTrace("psModules.detrend.cont",5, "OTA: %d (%d %d) A: %f B: %f C: %f D: %f",
                 i,x,y,
-                A->data.F32[y][x],B->data.F32[y][x],C->data.F32[y][x],D->data.F32[y][x]);               
+                A->data.F64[y][x],B->data.F64[y][x],C->data.F64[y][x],D->data.F64[y][x]);               
     }
     psFree(bgStats);
@@ -661 +673 @@
       int x = (i - y) / 8;
       int j;
+      double critical_value = 0.0;
       if (meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PM_PATTERN_IGNORE) {
         continue;
@@ -666 +679 @@
       if (x + 1 < 8) {  // We have a neighbor adjacent in the +x direction
         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]);
+        }
+        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
                 i,x,y,j,
@@ -672 +691 @@
                 );
         if (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_PATTERN_IGNORE)&&  // If there are no errors with the neighbor,
-            (isfinite(C->data.F64[y][x]))&&(isfinite(D->data.F64[y][x+1]))) {    // and all edges have valid values:
+            (isfinite(C->data.F64[y][x]))&&(isfinite(D->data.F64[y][x+1]))&&     // and all edges have valid values,
+            (fabs(C->data.F64[y][x] - D->data.F64[y][x+1]) < critical_value)     // and there are no large discontinuities,
+            ) {    
           solution->data.F64[i] += C->data.F64[y][x] - D->data.F64[y][x+1];     // Take the difference
           XX->data.F64[i][i] += 1;                                              // increment our relation with ourself
@@ -685 +706 @@
                 C->data.F64[y][x-1]
                 );
-
+        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]))) {
+            (isfinite(D->data.F64[y][x]))&&(isfinite(C->data.F64[y][x-1]))&&
+            (fabs(D->data.F64[y][x] - C->data.F64[y][x-1]) < critical_value)
+            ) {
           solution->data.F64[i] += D->data.F64[y][x] - C->data.F64[y][x-1];
           XX->data.F64[i][i] += 1;
@@ -700 +728 @@
                 B->data.F64[y+1][x]
                 );
-
+        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]);
+        }
+        else {
+          critical_value = 2.0 * fabs(A->data.F64[y][x] / B->data.F64[y+1][x]);
+        }
         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]))) {
+            (isfinite(A->data.F64[y][x]))&&(isfinite(B->data.F64[y+1][x]))&&
+            (fabs(A->data.F64[y][x] - B->data.F64[y+1][x]) < critical_value)
+            ) {
           solution->data.F64[i] += A->data.F64[y][x] - B->data.F64[y+1][x];
           XX->data.F64[i][i] += 1;
@@ -715 +750 @@
                 A->data.F64[y-1][x]
                 );
-
+        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]);
+        }
+        else {
+          critical_value = 2.0 * fabs(A->data.F64[y-1][x] / B->data.F64[y][x]);
+        }
         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]))) {
+            (isfinite(B->data.F64[y][x]))&&(isfinite(A->data.F64[y-1][x]))&&
+            (fabs(B->data.F64[y][x] - A->data.F64[y-1][x]) < critical_value)
+            ) {
           solution->data.F64[i] += B->data.F64[y][x] - A->data.F64[y-1][x];
           XX->data.F64[i][i] += 1;