Changeset 19765 for trunk/psModules/src/imcombine/pmSubtractionKernels.c

Timestamp:

Sep 26, 2008, 4:03:24 PM (18 years ago)

Author:

Paul Price

Message:

Two major changes. 1. Optimising ISIS kernels: instead of doing an FFT convolution, they're separable. 2. Trying to make the convolved weight map account for systematic errors in the convolution kernel by assuming a relative error in each kernel pixel (currently hardwired with a #define, but we can turn it into a recipe value later, or get an estimate from the data) and propagating that into the variance map. The effect is to increase the variance near the bright stars, which is hopefully sufficient to decrease the number of false sources detected in bad subtractions.

File:

• trunk/psModules/src/imcombine/pmSubtractionKernels.c (modified) (3 diffs)

trunk/psModules/src/imcombine/pmSubtractionKernels.c

@@ -42 +42 @@
     }
     return result;
+}
+
+// Generate 1D convolution kernel for ISIS
+static psVector *subtractionKernelISIS(float sigma, // Gaussian width
+                                       int order, // Polynomial order
+                                       int size // Kernel half-size
+    )
+{
+    int fullSize = 2 * size + 1;        // Full size of kernel
+    psVector *kernel = psVectorAlloc(fullSize, PS_TYPE_F32); // Kernel to return
+
+    float expNorm = -0.5 / PS_SQR(sigma); // Normalisation for exponential
+    float norm = 1.0 / (M_2_PI * sqrtf(sigma)); // Normalisation for Gaussian
+    for (int i = 0, x = -size; x <= size; i++, x++) {
+        kernel->data.F32[i] = norm * power(x, order) * expf(expNorm * PS_SQR(x));
+    }
+
+    return kernel;
 }
 
@@ -116 +134 @@
 
     // Set the kernel parameters
+    int fullSize = 2 * size + 1;        // Full size of kernels
     for (int i = 0, index = 0; i < numGaussians; i++) {
         float sigma = fwhms->data.F32[i] / (2.0 * sqrtf(2.0 * logf(2.0))); // Gaussian sigma
-        float norm = 1.0 / (M_2_PI * sqrtf(sigma)); // Normalisation for Gaussian
         // Iterate over (u,v) order
         for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
             for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
-                // Set the pre-calculated kernel
-                psKernel *preCalc = psKernelAlloc(-size, size, -size, size);
-                double sum = 0.0;       // Normalisation
+                psArray *preCalc = psArrayAlloc(2); // Array to hold precalculated values
+                psVector *xKernel = preCalc->data[0] = subtractionKernelISIS(sigma, uOrder, size); // x Kernel
+                psVector *yKernel = preCalc->data[1] = subtractionKernelISIS(sigma, vOrder, size); // y Kernel
+
+                // Calculate moments
+                double sum = 0.0;       // Sum of kernel component, for normalisation
                 double moment = 0.0;    // Moment, for penalty
-                for (int v = -size; v <= size; v++) {
-                    for (int u = -size; u <= size; u++) {
-                        sum += preCalc->kernel[v][u] = norm * power(u, uOrder) * power(v, vOrder) *
-                            expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigma));
-                        moment += preCalc->kernel[v][u] * (PS_SQR(u) + PS_SQR(v));
+                for (int v = -size, y = 0; v <= size; v++, y++) {
+                    for (int u = -size, x = 0; u <= size; u++, x++) {
+                        double value = xKernel->data.F32[x] * yKernel->data.F32[y]; // Value of kernel
+                        sum += value;
+                        moment += value * (PS_SQR(u) + PS_SQR(v));
                     }
                 }
+                moment = 0.0;
 
                 // Normalise sum of kernel component to unity for even functions
                 if (uOrder % 2 == 0 && vOrder % 2 == 0) {
-                    for (int v = -size; v <= size; v++) {
-                        for (int u = -size; u <= size; u++) {
-                            preCalc->kernel[v][u] = preCalc->kernel[v][u] / sum;
+                    double sum = 0.0;   // Sum of kernel component
+                    for (int v = 0; v < fullSize; v++) {
+                        for (int u = 0; u < fullSize; u++) {
+                            sum += xKernel->data.F32[u] * yKernel->data.F32[v];
                         }
                     }
-                    preCalc->kernel[0][0] -= 1.0;
+                    sum = 1.0 / sum;
+                    psBinaryOp(xKernel, xKernel, "*", psScalarAlloc(sum, PS_TYPE_F32));
+                    psBinaryOp(yKernel, yKernel, "*", psScalarAlloc(sum, PS_TYPE_F32));
+                    moment *= sum;
                 }
 
@@ -702 +728 @@
             // Count the number of Gaussians
             int numGauss = 0;
-            for (char *string = ptr; string; string = strchr(string, '(')) {
+            for (char *string = ptr; string; string = strchr(string + 1, '(')) {
                 numGauss++;
             }