Changeset 28304 for branches/czw_branch/20100519/psLib/src/imageops/psImageCovariance.c

Timestamp:

Jun 10, 2010, 6:28:51 PM (16 years ago)

Author:

watersc1

Message:

Skycell Summary and stack Association stuff should be finished. I'll merge and do final tests on monday.

Location:

branches/czw_branch/20100519

Files:

• . (modified) (1 prop)
• psLib/src/imageops/psImageCovariance.c (modified) (2 diffs)

branches/czw_branch/20100519/psLib/src/imageops/psImageCovariance.c

@@ -571 +571 @@
 
     psKernel *out = psKernelAlloc(xMinOut, xMaxOut, yMinOut, yMaxOut); // Output covariance
+    double outSum = 0.0;                                               // Sum of covariance
     for (int y = yMinOut; y <= yMaxOut; y++) {
         float yIn = y * scale + 0.5 - yMinIn + 1; // Position on input image (not the kernel)
         for (int x = xMinOut; x <= xMaxOut; x++) {
             float xIn = x * scale + 0.5 - xMinIn + 1; // Position on input (not the kernel)
+
             double value;                                     // Value on output
             if (!psImageInterpolate(&value, NULL, NULL, xIn, yIn, interp)) {
@@ -580 +582 @@
                 return false;
             }
-            out->kernel[y][x] = value;
-        }
-    }
-
+            outSum += out->kernel[y][x] = value;
+        }
+    }
     psFree(interp);
+
+    // Problem: the interpolation has introduced power into the covariance matrix that shouldn't be there.  We
+    // can't scale the sum of the matrix because that would throw off the central value and affect the noise
+    // calculation for this image.  But we can't scale just by the central value because that would throw off
+    // the noise calculation for convolutions of this image.  We choose to scale by the sum of the non-central
+    // elements.  This is almost having the best of both worlds.
+
+    double inSum = 0.0;                 // Sum of covariance
+    for (int y = yMinIn; y <= yMaxIn; y++) {
+        for (int x = xMinIn; x <= xMaxIn; x++) {
+            inSum += in->kernel[y][x];
+        }
+    }
+
+    float norm = (inSum - in->kernel[0][0]) / (outSum - out->kernel[0][0]) / PS_SQR(scale); // Renormalisation
+    for (int y = yMinOut; y <= yMaxOut; y++) {
+        for (int x = xMinOut; x <= xMaxOut; x++) {
+            if (x != 0 && y != 0) {
+                out->kernel[y][x] *= norm;
+            }
+        }
+    }
 
     return out;