Changeset 28152 for trunk/psLib/src/imageops/psImageCovariance.c

Timestamp:

May 27, 2010, 7:34:23 PM (16 years ago)

Author:

Paul Price

Message:

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.

File:

• trunk/psLib/src/imageops/psImageCovariance.c (modified) (3 diffs)

trunk/psLib/src/imageops/psImageCovariance.c

@@ -576 +576 @@
         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)) {
@@ -586 +587 @@
     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++) {
@@ -593 +600 @@
     }
 
-    float norm = inSum / PS_SQR(scale) / outSum; // Renormalisation (to remove errors in interp.)
+    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++) {
-            out->kernel[y][x] *= norm;
+            if (x != 0 && y != 0) {
+                out->kernel[y][x] *= norm;
+            }
         }
     }