Changeset 26762

Timestamp:

Feb 2, 2010, 11:45:40 AM (16 years ago)

Author:

Paul Price

Message:

Adding function, psKernelTruncate, to truncate a kernel so as to lose less than some nominated fraction of the total absolute flux.

Location:

branches/eam_branches/20091201/psLib/src/imageops

Files:

• psImageConvolve.c (modified) (1 diff)
• psImageConvolve.h (modified) (1 diff)

branches/eam_branches/20091201/psLib/src/imageops/psImageConvolve.c

@@ -246 +246 @@
     return kernel;
 }
+
+bool psKernelTruncate(psKernel *kernel, float frac)
+{
+    PS_ASSERT_KERNEL_NON_NULL(kernel, false);
+    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(frac, 0.0, false);
+    PS_ASSERT_FLOAT_LESS_THAN(frac, 1.0, false);
+
+    if (frac == 0.0) {
+        // Nothing to do
+        return true;
+    }
+
+    int xMin = kernel->xMin, xMax = kernel->xMax, yMin = kernel->yMin, yMax = kernel->yMax; // Bounds
+    int maxSize = PS_MAX(PS_MAX(PS_MAX(-xMin, xMax), -yMin), yMax); // Maximum size
+
+    // Determine the threshold
+    // Summing absolute values because large negative deviations have power as well
+    double sumKernel = 0.0;             // Sum of the kernel
+    for (int y = yMin; y <= yMax; y++) {
+        for (int x = xMin; x <= xMax; x++) {
+            sumKernel += fabsf(kernel->kernel[y][x]);
+        }
+    }
+
+    float threshold = sumKernel * (1.0 - frac); // Threshold for truncation
+
+    // Find truncation size
+    int truncate = 0;                     // Truncation radius
+    for (int radius = 1; truncate == 0 && radius < maxSize; radius++) {
+        int uMin = PS_MAX(-radius, xMin);
+        int uMax = PS_MIN(radius, xMax);
+        int vMin = PS_MAX(-radius, yMin);
+        int vMax = PS_MIN(radius, yMax);
+        int r2 = PS_SQR(radius);
+        double sum = 0.0;
+        for (int v = vMin; v <= vMax; v++) {
+            int v2 = PS_SQR(v);
+            for (int u = uMin; u <= uMax; u++) {
+                int u2 = PS_SQR(u);
+                if (u2 + v2 <= r2) {
+                    sum += fabsf(kernel->kernel[v][u]);
+                }
+            }
+        }
+        if (sum > threshold) {
+            // This is the truncation radius
+            truncate = radius;
+        }
+    }
+    if (truncate == maxSize) {
+        // No truncation possible
+        return true;
+    }
+
+    // Truncate the kernel
+    {
+        int uMin = PS_MAX(-truncate, xMin);
+        int uMax = PS_MIN(truncate, xMax);
+        int vMin = PS_MAX(-truncate, yMin);
+        int vMax = PS_MIN(truncate, yMax);
+        int r2 = PS_SQR(truncate);
+        for (int v = vMin; v <= vMax; v++) {
+            int v2 = PS_SQR(v);
+            for (int u = uMin; u <= uMax; u++) {
+                int u2 = PS_SQR(u);
+                if (u2 + v2 > r2) {
+                    kernel->kernel[v][u] = 0.0;
+                }
+            }
+        }
+    }
+    kernel->xMin = PS_MAX(-truncate, kernel->xMin);
+    kernel->xMax = PS_MIN(truncate, kernel->xMax);
+    kernel->yMin = PS_MAX(-truncate, kernel->yMin);
+    kernel->yMax = PS_MIN(truncate, kernel->yMax);
+
+    return true;
+    }
+
+
 
 psImage *psImageConvolveDirect(psImage *out, const psImage *in, const psKernel *kernel)

branches/eam_branches/20091201/psLib/src/imageops/psImageConvolve.h

@@ -138 +138 @@
 );
 
+/// Truncate a kernel
+///
+/// Truncates the outer parts of the kernel where the contribution is below the nominated fraction of the
+/// total kernel.
+bool psKernelTruncate(
+    psKernel *in,                       ///< Kernel to be truncated
+    float frac                          ///< Fraction for truncation threshold
+    );
+
+
 /// Convolve an image with a kernel, using a direct convolution
 ///