Index: trunk/psModules/src/imcombine/Makefile.am =================================================================== --- trunk/psModules/src/imcombine/Makefile.am (revision 20048) +++ trunk/psModules/src/imcombine/Makefile.am (revision 20049) @@ -9,4 +9,5 @@ pmStackReject.c \ pmSubtraction.c \ + pmSubtractionAnalysis.c \ pmSubtractionEquation.c \ pmSubtractionIO.c \ @@ -25,4 +26,5 @@ pmStackReject.h \ pmSubtraction.h \ + pmSubtractionAnalysis.h \ pmSubtractionEquation.h \ pmSubtractionIO.h \ Index: trunk/psModules/src/imcombine/pmSubtractionAnalysis.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionAnalysis.c (revision 20049) +++ trunk/psModules/src/imcombine/pmSubtractionAnalysis.c (revision 20049) @@ -0,0 +1,224 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "pmFPA.h" +#include "pmSubtraction.h" +#include "pmSubtractionKernels.h" +#include "pmSubtractionEquation.h" + +#include "pmSubtractionAnalysis.h" + +#define KERNEL_MOSAIC 2 // Half-number of kernel instances in the mosaic image + + +bool pmSubtractionAnalysis(psMetadata *analysis, pmSubtractionKernels *kernels, psRegion *region, + int numCols, int numRows) +{ + if (analysis) { + PS_ASSERT_METADATA_NON_NULL(analysis, false); + } + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); + PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, false); + + // Record region for subtraction + { + psRegion *subRegion; // Region over which subtraction was performed + if (region) { + subRegion = psMemIncrRefCounter(region); + } else { + subRegion = psRegionAlloc(0, numCols, 0, numRows); + } + + psMetadataAddPtr(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION, + PS_DATA_REGION | PS_META_DUPLICATE_OK, + "Region over which subtraction was performed", subRegion); + psFree(subRegion); + } + + // Record kernel + psMetadataAddPtr(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_KERNEL, + PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels); + psMetadataAddS32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE, + PS_META_DUPLICATE_OK, "Subtraction kernels", kernels->mode); + + // Realisations of kernel + { + psTrace("psModules.imcombine", 2, "Generating diagnostics...\n"); + // Generate image with convolution kernels + int size = kernels->size; // Half-size of kernel + int fullSize = 2 * size + 1 + 1; // Full size of kernel + int imageSize = (2 * KERNEL_MOSAIC + 1) * fullSize; + psImage *convKernels = psImageAlloc((kernels->mode == PM_SUBTRACTION_MODE_DUAL ? 2 : 1) * + imageSize - 1 + + (kernels->mode == 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); + return false; + } + + 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); + return false; + } + psFree(kernel); + + if (kernels->mode == 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); + return false; + } + + 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); + return false; + } + psFree(kernel); + } + } + } + + psMetadataAddImage(analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL.IMAGE", + PS_META_DUPLICATE_OK, "Realisations of kernel", convKernels); + 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 + + + // Set the variance factors + float vf1 = 1.0, vf2 = 1.0; // Variance factors for each image + switch (kernels->mode) { + case PM_SUBTRACTION_MODE_1: + vf1 = pmSubtractionVarianceFactor(kernels, 0.5, 0.5, false); + break; + case PM_SUBTRACTION_MODE_2: + vf2 = pmSubtractionVarianceFactor(kernels, 0.5, 0.5, false); + break; + case PM_SUBTRACTION_MODE_DUAL: + vf1 = pmSubtractionVarianceFactor(kernels, 0.5, 0.5, false); + vf2 = pmSubtractionVarianceFactor(kernels, 0.5, 0.5, true); + break; + default: + psAbort("Invalid subtraction mode: %x", kernels->mode); + } + + // Weight by the area + if (region) { + float norm = (region->x1 - region->x0 + 1) * (region->y1 - region->y0 + 1) / (numCols * numRows); + vf1 *= norm; + vf2 *= norm; + } + + // Update the variance factor +#define UPDATE_VARFACTOR(VF, ANALYSIS) { \ + psMetadataItem *vfItem = psMetadataLookup(analysis, ANALYSIS); \ + if (vfItem) { \ + psAssert(vfItem->type == PS_TYPE_F32, "Should be the type we said."); \ + vfItem->data.F32 += VF; \ + } else { \ + psMetadataAddF32(analysis, PS_LIST_TAIL, ANALYSIS, 0, "Variance factor weighted by the area", VF); \ + } \ +} + + UPDATE_VARFACTOR(vf1, PM_SUBTRACTION_ANALYSIS_VARFACTOR_1); + UPDATE_VARFACTOR(vf2, PM_SUBTRACTION_ANALYSIS_VARFACTOR_2); + + // Kernel moments + { + 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; + } + double m00 = 0, m10 = 0, m01 = 0, m20 = 0, m11 = 0, m02 = 0; // Moments to calculate + int size = kernels->size; // Half-size of kernel + for (int y = -size; y <= size; y++) { + for (int x = -size; x <= size; x++) { + float value = image->data.F32[y][x]; // Value of kernel + m00 += value; + m10 += x * value; + m01 += y * value; + m20 += x * x * value; + m11 += x * y * value; + m02 += y * y * value; + } + } + psFree(image); + + // Convert first moments to centroids + m10 /= m00; + m01 /= m00; + + // Convert second moments to covariance + m20 = m20 / m00 - PS_SQR(m10); + m02 = m02 / m00 - PS_SQR(m01); + m11 = m11 / m00 - m10 * m01; + + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_NORM, + PS_META_DUPLICATE_OK, "Normalisation", m00); + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MX, + PS_META_DUPLICATE_OK, "Moment in x", m10); + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MY, + PS_META_DUPLICATE_OK, "Moment in y", m01); + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MXX, + PS_META_DUPLICATE_OK, "Moment in xx", m20); + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MXY, + PS_META_DUPLICATE_OK, "Moment in xy", m11); + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MYY, + PS_META_DUPLICATE_OK, "Moment in yy", m02); + } + + // Difference in background + { + psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, numCols, numRows, + numCols / 2.0, numRows / 2.0); // Polynomial + float bg = p_pmSubtractionSolutionBackground(kernels, polyValues); // Background difference + + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_BGDIFF, + PS_META_DUPLICATE_OK, "Background difference", bg); + psFree(polyValues); + } + + return true; +} Index: trunk/psModules/src/imcombine/pmSubtractionAnalysis.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionAnalysis.h (revision 20049) +++ trunk/psModules/src/imcombine/pmSubtractionAnalysis.h (revision 20049) @@ -0,0 +1,31 @@ +#ifndef PM_SUBTRACTION_ANALYSIS_H +#define PM_SUBTRACTION_ANALYSIS_H + +#include +#include + +// Names for things put on the readout analysis metadata +#define PM_SUBTRACTION_ANALYSIS_KERNEL "SUBTRACTION.KERNEL" // Kernel used for convolving +#define PM_SUBTRACTION_ANALYSIS_KERNEL_IMAGE "SUBTRACTION.KERNEL.IMAGE" // Image with kernel realisations +#define PM_SUBTRACTION_ANALYSIS_MODE "SUBTRACTION.MODE" // Subtraction mode +#define PM_SUBTRACTION_ANALYSIS_REGION "SUBTRACTION.REGION" // Subtraction region +#define PM_SUBTRACTION_ANALYSIS_VARFACTOR_1 "SUBTRACTION.VARFACTOR.1" // Variance factor for image 1 +#define PM_SUBTRACTION_ANALYSIS_VARFACTOR_2 "SUBTRACTION.VARFACTOR.2" // Variance factor for image 2 +#define PM_SUBTRACTION_ANALYSIS_NORM "SUBTRACTION.NORM" // Normalisation +#define PM_SUBTRACTION_ANALYSIS_BGDIFF "SUBTRACTION.BGDIFF" // Background difference +#define PM_SUBTRACTION_ANALYSIS_MX "SUBTRACTION.MX" // Kernel moment in x +#define PM_SUBTRACTION_ANALYSIS_MY "SUBTRACTION.MY" // Kernel moment in y +#define PM_SUBTRACTION_ANALYSIS_MXX "SUBTRACTION.MXX" // Kernel moment in xx +#define PM_SUBTRACTION_ANALYSIS_MXY "SUBTRACTION.MXY" // Kernel moment in xy +#define PM_SUBTRACTION_ANALYSIS_MYY "SUBTRACTION.MYY" // Kernel moment in yy + +// Derive QA information about the subtraction +bool pmSubtractionAnalysis( + psMetadata *analysis, ///< Metadata container for QA information + pmSubtractionKernels *kernels, ///< Kernels + psRegion *region, ///< Region for subtraction + int numCols, int numRows ///< Size of image + ); + + +#endif Index: trunk/psModules/src/imcombine/pmSubtractionIO.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionIO.c (revision 20048) +++ trunk/psModules/src/imcombine/pmSubtractionIO.c (revision 20049) @@ -13,4 +13,7 @@ #include "pmSubtractionKernels.h" #include "pmSubtractionMatch.h" +#include "pmSubtractionAnalysis.h" + +#include "pmSubtractionIO.h" #define ARRAY_BUFFER 16 // Number to add to array at a time @@ -140,6 +143,6 @@ // CVS tags, used to identify the version of this file (in case incompatibilities are introduced) psString cvsFile = psStringCopy("$RCSfile: pmSubtractionIO.c,v $"); - psString cvsRev = psStringCopy("$Revision: 1.7 $"); - psString cvsDate = psStringCopy("$Date: 2008-10-01 19:55:09 $"); + psString cvsRev = psStringCopy("$Revision: 1.8 $"); + psString cvsDate = psStringCopy("$Date: 2008-10-10 23:54:33 $"); psStringSubstitute(&cvsFile, NULL, "RCSfile: "); psStringSubstitute(&cvsRev, NULL, "Revision: "); Index: trunk/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 20048) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 20049) @@ -2,4 +2,5 @@ #define PM_SUBTRACTION_KERNELS_H +#include #include Index: trunk/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 20048) +++ 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); Index: trunk/psModules/src/imcombine/pmSubtractionMatch.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMatch.h (revision 20048) +++ trunk/psModules/src/imcombine/pmSubtractionMatch.h (revision 20049) @@ -8,11 +8,4 @@ #include #include - -// Names for things put on the readout analysis metadata -#define PM_SUBTRACTION_ANALYSIS_KERNEL "SUBTRACTION.KERNEL" // Kernel used for convolving -#define PM_SUBTRACTION_ANALYSIS_MODE "SUBTRACTION.MODE" // Subtraction mode -#define PM_SUBTRACTION_ANALYSIS_REGION "SUBTRACTION.REGION" // Subtraction region -#define PM_SUBTRACTION_ANALYSIS_VARFACTOR "SUBTRACTION.VARFACTOR"// Variance factor - /// Match two images Index: trunk/psModules/src/psmodules.h =================================================================== --- trunk/psModules/src/psmodules.h (revision 20048) +++ trunk/psModules/src/psmodules.h (revision 20049) @@ -89,4 +89,5 @@ #include #include +#include #include #include