Changeset 21149 for trunk/psModules/src/imcombine/pmSubtractionAnalysis.c

Timestamp:

Jan 22, 2009, 9:03:46 AM (17 years ago)

Author:

Paul Price

Message:

Adding test for deconvolution: measures the fraction of the kernel flux within a range of apertures about the centre and takes the maximum; if significantly in excess of 1 (e.g., 1.1), the convolved images are generally bad.

File:

• trunk/psModules/src/imcombine/pmSubtractionAnalysis.c (modified) (1 diff)

trunk/psModules/src/imcombine/pmSubtractionAnalysis.c

@@ -165 +165 @@
     UPDATE_VARFACTOR(vf2, PM_SUBTRACTION_ANALYSIS_VARFACTOR_2);
 
+    // Kernel shape
+    {
+        psImage *image = pmSubtractionKernelImage(kernels, 0.5, 0.5, false); // Image of the kernel
+        if (!image) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to generate image of kernel.");
+            return false;
+        }
+        int size = kernels->size;       // Half-size of kernel
+        int fullSize = 2 * size + 1;    // Full size of kernel
+        float norm = 0.0;               // Normalisation (kernel sum)
+        for (int y = 0; y < fullSize; y++) {
+            for (int x = 0; x < fullSize; x++) {
+                norm += image->data.F32[y][x];
+            }
+        }
+        float max = -INFINITY;          // Maximum fraction
+        for (int r = 1; r < size; r++) {
+            unsigned long r2 = PS_SQR(r); // r^2
+            float sum = 0.0;            // Sum within circle
+            for (int y = 0, v = -size; y < fullSize; y++, v++) {
+                unsigned long v2 = PS_SQR(v); // y^2
+                for (int x = 0, u = -size; x < fullSize; x++, u++) {
+                    unsigned long u2 = PS_SQR(u); // u^2
+                    if (u2 + v2 <= r2) {
+                        sum += image->data.F32[y][x];
+                    }
+                }
+            }
+            float frac = sum / norm;    // Fraction of flux moving towards centre
+            psTrace("psModules.imcombine", 5, "Deconvolution fraction at %d: %f\n", r, frac);
+            max = PS_MAX(max, frac);
+        }
+        psFree(image);
+
+        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DECONV_MAX,
+                         PS_META_DUPLICATE_OK, "Maximum deconvolution fraction", max);
+    }
+
     // Kernel moments
     {