Changeset 21382

Timestamp:

Feb 6, 2009, 12:39:15 PM (17 years ago)

Author:

Paul Price

Message:

More efficient determination of truncation radius.
Sums need to be of absolute values, since covariance can be negative.

File:

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

trunk/psLib/src/imageops/psImageCovariance.c

@@ -4 +4 @@
 
 #include <stdio.h>
+#include <stdlib.h>
 #include <limits.h>
 #include <string.h>
@@ -181 +182 @@
 
     int xMin = covar->xMin, xMax = covar->xMax, yMin = covar->yMin, yMax = covar->yMax; // Range
+    int maxRadius = PS_MAX(PS_MAX(PS_MAX(xMax, yMax), -xMin), yMin); // Maximum radius of covariance matrix
 
     double sum = 0.0;                   // Sum of covariance
+    psVector *radiusSum = psVectorAlloc(maxRadius, PS_TYPE_F64); // Totals within (square) radius
+    psVectorInit(radiusSum, 0.0);
     for (int y = yMin; y <= yMax; y++) {
         for (int x = xMin; x <= xMax; x++) {
-            sum += covar->kernel[y][x];
+            int radius = PS_MAX(abs(x), abs(y)); // Squarish radius
+            radiusSum->data.F64[radius] += fabsf(covar->kernel[y][x]);
+            sum += fabsf(covar->kernel[y][x]);
         }
     }
@@ -191 +197 @@
     // 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 = 0.0;              // 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];
-            }
-        }
-    }
+        enclosed += radiusSum->data.F64[radius];
+    }
+    psFree(radiusSum);
 
     // Generate truncated version