Index: trunk/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 19532) +++ trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 20049) @@ -15,4 +15,5 @@ #include "pmSubtractionEquation.h" #include "pmSubtraction.h" +#include "pmSubtractionAnalysis.h" #include "pmSubtractionMask.h" #include "pmSubtractionThreads.h" @@ -20,5 +21,4 @@ -#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,27 +91,4 @@ -// 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,4 +231,6 @@ } + psMetadata *analysis = psMetadataAlloc(); // QA data + // Iterate over iso-kernel regions for (int j = 0; j < yRegions; j++) { @@ -338,34 +317,4 @@ } - // 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,85 +384,8 @@ 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,25 +396,4 @@ 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,4 +434,11 @@ memCheck("convolution"); + if (conv1) { + psMetadataCopy(conv1->analysis, analysis); + } + if (conv2) { + psMetadataCopy(conv2->analysis, analysis); + } + psFree(analysis); #ifdef TESTING @@ -603,4 +461,5 @@ MATCH_ERROR: + psFree(analysis); psFree(region); psFree(regionString);