Changeset 20049 for trunk/psModules/src/imcombine/pmSubtractionMatch.c

Timestamp:

Oct 10, 2008, 1:54:33 PM (18 years ago)

Author:

Paul Price

Message:

Adding new file with function to add stuff to a common 'analysis' metadata container which is applied to the convolved outputs when done.

File:

• trunk/psModules/src/imcombine/pmSubtractionMatch.c (modified) (9 diffs)

trunk/psModules/src/imcombine/pmSubtractionMatch.c

@@ -15 +15 @@
 #include "pmSubtractionEquation.h"
 #include "pmSubtraction.h"
+#include "pmSubtractionAnalysis.h"
 #include "pmSubtractionMask.h"
 #include "pmSubtractionThreads.h"
@@ -20 +21 @@
 
 
-#define KERNEL_MOSAIC 2                 // Half-number of kernel instances in the mosaic image
 #define BG_STAT PS_STAT_ROBUST_MEDIAN   // Statistic to use for background
 
@@ -91 +91 @@
 
 
-// Record the variance factor in a readout
-static bool recordVarianceFactor(pmReadout *ro, // Readout of interest
-                                 float varFactor, // Variance factor
-                                 const psRegion *region // Region of interest
-                                 )
-{
-    if (region) {
-        varFactor *= (region->x1 - region->x0 + 1) * (region->y1 - region->y0 + 1) /
-            (ro->image->numCols * ro->image->numRows);
-    }
-
-    psMetadataItem *vfItem = psMetadataLookup(ro->analysis, PM_SUBTRACTION_ANALYSIS_VARFACTOR);
-    if (vfItem) {
-        psAssert(vfItem->type == PS_TYPE_F32, "Should be the type we said.");
-        vfItem->data.F32 += varFactor;
-    } else {
-        psMetadataAddF32(ro->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_VARFACTOR, 0,
-                         "Variance factor weighted by the area", varFactor);
-    }
-
-    return true;
-}
-
 bool pmSubtractionMatch(pmReadout *conv1, pmReadout *conv2, const pmReadout *ro1, const pmReadout *ro2,
                         int footprint, float regionSize, float stampSpacing, float threshold,
@@ -254 +231 @@
     }
 
+    psMetadata *analysis = psMetadataAlloc(); // QA data
+
     // Iterate over iso-kernel regions
     for (int j = 0; j < yRegions; j++) {
@@ -338 +317 @@
             }
 
-            // Add analysis metadata
-            {
-                psRegion *subRegion;    // Region over which subtraction was performed
-                if (region) {
-                    subRegion = psMemIncrRefCounter(region);
-                } else {
-                    subRegion = psRegionAlloc(0, numCols, 0, numRows);
-                }
-
-                if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_DUAL) {
-                    psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_KERNEL,
-                                     PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels);
-                    psMetadataAddS32(conv1->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE,
-                                     PS_META_DUPLICATE_OK, "Subtraction kernels", subMode);
-                    psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION,
-                                     PS_DATA_REGION | PS_META_DUPLICATE_OK,
-                                     "Region over which subtraction was performed", subRegion);
-                }
-                if (subMode == PM_SUBTRACTION_MODE_2 || subMode == PM_SUBTRACTION_MODE_DUAL) {
-                    psMetadataAddPtr(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_KERNEL,
-                                     PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels);
-                    psMetadataAddS32(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE,
-                                     PS_META_DUPLICATE_OK, "Subtraction kernels", subMode);
-                    psMetadataAddPtr(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION,
-                                     PS_DATA_REGION | PS_META_DUPLICATE_OK,
-                                     "Region over which subtraction was performed", subRegion);
-                }
-                psFree(subRegion);
-            }
-
             memCheck("kernels");
 
@@ -435 +384 @@
             memCheck("solution");
 
-            {
-                psTrace("psModules.imcombine", 2, "Generating diagnostics...\n");
-                // Generate image with convolution kernels
-                int fullSize = 2 * size + 1 + 1;    // Full size of kernel
-                int imageSize = (2 * KERNEL_MOSAIC + 1) * fullSize;
-                psImage *convKernels = psImageAlloc((subMode == PM_SUBTRACTION_MODE_DUAL ? 2 : 1) *
-                                                    imageSize - 1 +
-                                                    (subMode == PM_SUBTRACTION_MODE_DUAL ? 4 : 0),
-                                                    imageSize - 1, PS_TYPE_F32);
-                psImageInit(convKernels, NAN);
-                for (int j = -KERNEL_MOSAIC; j <= KERNEL_MOSAIC; j++) {
-                    for (int i = -KERNEL_MOSAIC; i <= KERNEL_MOSAIC; i++) {
-                        psImage *kernel = pmSubtractionKernelImage(kernels, (float)i / (float)KERNEL_MOSAIC,
-                                                                   (float)j / (float)KERNEL_MOSAIC,
-                                                                   false); // Image of the kernel
-                        if (!kernel) {
-                            psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");
-                            psFree(convKernels);
-                            goto MATCH_ERROR;
-                        }
-
-                        if (psImageOverlaySection(convKernels, kernel, (i + KERNEL_MOSAIC) * fullSize,
-                                                  (j + KERNEL_MOSAIC) * fullSize, "=") == 0) {
-                            psError(PS_ERR_UNKNOWN, false, "Unable to overlay kernel image.");
-                            psFree(kernel);
-                            psFree(convKernels);
-                            goto MATCH_ERROR;
-                        }
-                        psFree(kernel);
-
-                        if (subMode == PM_SUBTRACTION_MODE_DUAL) {
-                            kernel = pmSubtractionKernelImage(kernels, (float)i / (float)KERNEL_MOSAIC,
-                                                              (float)j / (float)KERNEL_MOSAIC,
-                                                              true); // Image of the kernel
-                            if (!kernel) {
-                                psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");
-                                psFree(convKernels);
-                                goto MATCH_ERROR;
-                            }
-
-                            if (psImageOverlaySection(convKernels, kernel,
-                                                      (2 * KERNEL_MOSAIC + 1 + i + KERNEL_MOSAIC) * fullSize +
-                                                      4,
-                                                      (j + KERNEL_MOSAIC) * fullSize, "=") == 0) {
-                                psError(PS_ERR_UNKNOWN, false, "Unable to overlay kernel image.");
-                                psFree(kernel);
-                                psFree(convKernels);
-                                goto MATCH_ERROR;
-                            }
-                            psFree(kernel);
-                        }
-
-                    }
-                }
-
-                psString comment = NULL; // Comment for metadata
-                psStringAppend(&comment, "Subtraction kernel for region %s", regionString);
-                psMetadataAddImage(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL.IMAGE",
-                                   PS_META_DUPLICATE_OK, comment, convKernels);
-                psFree(comment);
-                psFree(convKernels);
-            }
-
-#if 0
-            {
-                // Generate images of the kernel components
-                psMetadata *header = psMetadataAlloc(); // Header
-                for (int i = 0; i < solution->n; i++) {
-                    psString name = NULL;       // Header keyword
-                    psStringAppend(&name, "SOLN%04d", i);
-                    psMetadataAddF64(header, PS_LIST_TAIL, name, 0, NULL, solution->data.F64[i]);
-                    psFree(name);
-                }
-                psArray *kernelImages = pmSubtractionKernelSolutions(solution, kernels, 0.0, 0.0);
-                psFits *kernelFile = psFitsOpen("kernels.fits", "w");
-                (void)psFitsWriteImageCube(kernelFile, header, kernelImages, NULL);
-                psFitsClose(kernelFile);
-                psFree(kernelImages);
-                psFree(header);
-            }
-#endif
+            if (!pmSubtractionAnalysis(analysis, kernels, region, numCols, numRows)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to generate QA data");
+                goto MATCH_ERROR;
+            }
 
             memCheck("diag outputs");
@@ -524 +396 @@
                 psError(PS_ERR_UNKNOWN, false, "Unable to convolve image.");
                 goto MATCH_ERROR;
-            }
-
-            // Set the variance factors
-            switch (subMode) {
-              case PM_SUBTRACTION_MODE_1: {
-                  recordVarianceFactor(conv1, pmSubtractionVarianceFactor(kernels, 0.5, 0.5, false),
-                                       region);
-                  break;
-              }
-              case PM_SUBTRACTION_MODE_2:
-                recordVarianceFactor(conv1, pmSubtractionVarianceFactor(kernels, 0.5, 0.5, false),
-                                     region);
-                break;
-              case PM_SUBTRACTION_MODE_DUAL:
-                recordVarianceFactor(conv1, pmSubtractionVarianceFactor(kernels, 0.5, 0.5, false),
-                                     region);
-                recordVarianceFactor(conv2, pmSubtractionVarianceFactor(kernels, 0.5, 0.5, true),
-                                     region);
-                break;
-              default:
-                psAbort("Should never reach here.");
             }
 
@@ -583 +434 @@
     memCheck("convolution");
 
+    if (conv1) {
+        psMetadataCopy(conv1->analysis, analysis);
+    }
+    if (conv2) {
+        psMetadataCopy(conv2->analysis, analysis);
+    }
+    psFree(analysis);
 
 #ifdef TESTING
@@ -603 +461 @@
 
 MATCH_ERROR:
+    psFree(analysis);
     psFree(region);
     psFree(regionString);