Changeset 33316

Timestamp:

Feb 21, 2012, 11:15:01 AM (14 years ago)

Author:

watersc1

Message:

Fix pattern cotinuity correction when one cell has impossible values.

Location:

tags/ipp-20120216/psModules/src/detrend

Files:

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

tags/ipp-20120216/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;