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

Timestamp:

Feb 5, 2009, 10:15:32 AM (17 years ago)

Author:

Paul Price

Message:

Adding function to truncate a covariance matrix. This should help to keep them from growing too large (we're using kernels of half-size ~ 15 pixels for stack, which makes the covariance matrix grow by 60 pixels, which puts a big stress on psImageCovarianceCalculate).

File:

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

trunk/psLib/src/imageops/psImageCovariance.c

@@ -5 +5 @@
 #include <stdio.h>
 #include <limits.h>
+#include <string.h>
 
 #include "psMemory.h"
@@ -172 +173 @@
     return sum;
 }
+
+
+psKernel *psImageCovarianceTruncate(const psKernel *covar, float frac)
+{
+    PS_ASSERT_KERNEL_NON_NULL(covar, NULL);
+    PS_ASSERT_FLOAT_WITHIN_RANGE(frac, 0, 1, NULL);
+
+    int xMin = covar->xMin, xMax = covar->xMax, yMin = covar->yMin, yMax = covar->yMax; // Range
+
+    double sum = 0.0;                   // Sum of covariance
+    for (int y = yMin; y <= yMax; y++) {
+        for (int x = xMin; x <= xMax; x++) {
+            sum += covar->kernel[y][x];
+        }
+    }
+
+    // Determine radius for truncation
+    double threshold = (1.0 - frac) * sum; // Threshold for truncation
+    int maxRadius = PS_MAX(PS_MAX(PS_MAX(xMax, yMax), -xMin), yMin); // Maximum radius of covariance matrix
+    double enclosed = NAN;              // Enclosed value
+    int radius;                         // Radius at which to truncate
+    for (radius = 0; radius <= maxRadius && enclosed < threshold; radius++) {
+        // Not so efficient to execute, but easy to code; it should be fast anyway
+        enclosed = 0.0;
+        for (int y = PS_MAX(yMin, -radius); y <= PS_MIN(yMax, radius); y++) {
+            for (int x = PS_MAX(xMin, -radius); x <= PS_MIN(xMax, radius); x++) {
+                enclosed += covar->kernel[y][x];
+            }
+        }
+    }
+
+    // Generate truncated version
+    psKernel *trunc = psKernelAlloc(-radius, radius, -radius, radius); // Truncated covariance matrix
+    int numBytes = (2 * radius + 1) * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy
+    for (int y = -radius; y <= radius; y++) {
+        memcpy(trunc->kernel[y], covar->kernel[y], numBytes);
+    }
+
+    return trunc;
+}