Index: trunk/psModules/src/imcombine/pmSubtractionAnalysis.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionAnalysis.c (revision 29547) +++ trunk/psModules/src/imcombine/pmSubtractionAnalysis.c (revision 29594) @@ -14,8 +14,8 @@ #include "pmSubtractionVisual.h" -#define KERNEL_MOSAIC 2 // Half-number of kernel instances in the mosaic image - //#define TESTING +// save information about the kernel in the output header. this function also generates the +// image normalization, used by ppSubMatchPSFs.c to rescale the output image. bool pmSubtractionAnalysis(psMetadata *analysis, psMetadata *header, pmSubtractionKernels *kernels, psRegion *region, @@ -36,84 +36,22 @@ } - psMetadataAddPtr(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION, - PS_DATA_REGION | PS_META_DUPLICATE_OK, - "Region over which subtraction was performed", subRegion); + psMetadataAddPtr(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION, PS_DATA_REGION | PS_META_DUPLICATE_OK, "Region over which subtraction was performed", subRegion); psString string = psRegionToString(*subRegion); psFree(subRegion); - psMetadataAddStr(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION, PS_META_DUPLICATE_OK, - "Region over which subtraction was performed", string); + psMetadataAddStr(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION, PS_META_DUPLICATE_OK, "Region over which subtraction was performed", string); psFree(string); } // 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 mode", kernels->mode); - psMetadataAddS32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE, - PS_META_DUPLICATE_OK, "Subtraction mode", kernels->mode); + 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 mode", kernels->mode); + psMetadataAddS32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE, PS_META_DUPLICATE_OK, "Subtraction mode", 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(psErrorCodeLast(), false, "Unable to generate kernel image."); - psFree(convKernels); - return false; - } - - if (psImageOverlaySection(convKernels, kernel, (i + KERNEL_MOSAIC) * fullSize, - (j + KERNEL_MOSAIC) * fullSize, "=") == 0) { - psError(psErrorCodeLast(), 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(psErrorCodeLast(), 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(psErrorCodeLast(), false, "Unable to overlay kernel image."); - psFree(kernel); - psFree(convKernels); - return false; - } - psFree(kernel); - } - } - } - - pmSubtractionVisualPlotConvKernels(convKernels); - psMetadataAddImage(analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL.IMAGE", - PS_META_DUPLICATE_OK, "Realisations of kernel", convKernels); - psFree(convKernels); - } + psImage *convKernels = pmSubtractionKernelsImageMosaic(kernels); + psMetadataAddImage(analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL.IMAGE", PS_META_DUPLICATE_OK, "Realisations of kernel", convKernels); + psFree(convKernels); // sample difference images @@ -166,5 +104,16 @@ int size = kernels->size; // Half-size of kernel int fullSize = 2 * size + 1; // Full size of kernel + + /* in the CENTRAL_DELTA case, the final convolution kernel replaces the central delta + function which was subtracted to force zero flux for each kernel. for other cases, + we have to include the zero order component of the normalization here. + */ + + //# if (CENTRAL_DELTA) +# if (1) float norm = 0.0; // Normalisation (kernel sum) +# else + float norm = 1.0; // Normalisation (kernel sum) +# endif for (int y = 0; y < fullSize; y++) { for (int x = 0; x < fullSize; x++) { @@ -172,5 +121,11 @@ } } - float max = -INFINITY; // Maximum fraction + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Kernel Integral: %f", norm); + + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_NORM, PS_META_DUPLICATE_OK, "Normalisation", norm); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_NORM, PS_META_DUPLICATE_OK, "Normalisation", norm); + + float maxCon = 0.0; // Maximum +fraction > 0.0 + float maxDec = 0.0; // Maximum -fraction < 0.0 for (int r = 1; r < size; r++) { unsigned long r2 = PS_SQR(r); // r^2 @@ -186,31 +141,44 @@ } 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); + psTrace("psModules.imcombine", 5, "(De)Convolution fraction at %d: %f\n", r, frac); + maxCon = PS_MAX(maxCon, +frac); + maxDec = PS_MAX(maxDec, -frac); } psFree(image); { - psMetadataItem *item = psMetadataLookup(analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Previous - if (item) { - max = item->data.F32 = PS_MAX(item->data.F32, max); - } else { - psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DECONV_MAX, 0, - "Maximum deconvolution fraction", max); - } - } - + psMetadataItem *item = NULL; + item = psMetadataLookup(analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Previous + if (item) { + maxDec = item->data.F32 = PS_MAX(item->data.F32, maxDec); + } else { + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DECONV_MAX, 0, "Maximum deconvolution fraction", maxDec); + } + item = psMetadataLookup(header, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Previous + if (item) { + item->data.F32 = maxDec; + } else { + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DECONV_MAX, 0, "Maximum deconvolution fraction", maxDec); + } + } { - psMetadataItem *item = psMetadataLookup(header, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Previous - if (item) { - item->data.F32 = max; - } else { - psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DECONV_MAX, 0, - "Maximum deconvolution fraction", max); - } - } - } - - // Kernel moments + psMetadataItem *item = NULL; + item = psMetadataLookup(analysis, PM_SUBTRACTION_ANALYSIS_CONVOL_MAX); // Previous + if (item) { + maxCon = item->data.F32 = PS_MAX(item->data.F32, maxCon); + } else { + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_CONVOL_MAX, 0, "Maximum convolution fraction", maxCon); + } + item = psMetadataLookup(header, PM_SUBTRACTION_ANALYSIS_CONVOL_MAX); // Previous + if (item) { + item->data.F32 = maxCon; + } else { + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_CONVOL_MAX, 0, "Maximum convolution fraction", maxCon); + } + } + } + + // Kernel moments : since the kernel can be negative, calculate the absolute flux moments, + // eg \Sum(x|f|) / Sum(|f|), \Sum(x^2|f|) / Sum(|f|) { psImage *image = pmSubtractionKernelImage(kernels, 0.5, 0.5, false); // Image of the kernel @@ -224,5 +192,5 @@ for (int y = 0, v = -size; y < fullSize; y++, v++) { for (int x = 0, u = -size; x < fullSize; x++, u++) { - float value = image->data.F32[y][x]; // Value of kernel + float value = fabs(image->data.F32[y][x]); // Value of kernel m00 += value; m10 += u * value; @@ -244,6 +212,4 @@ 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); @@ -257,6 +223,4 @@ PS_META_DUPLICATE_OK, "Moment in yy", m02); - psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_NORM, - PS_META_DUPLICATE_OK, "Normalisation", m00); psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MX, PS_META_DUPLICATE_OK, "Moment in x", m10); @@ -285,31 +249,24 @@ // Quality of fit { - psMetadataAddS32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_STAMPS, 0, "Number of stamps", - kernels->numStamps); - psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_MEAN, 0, "Mean stamp deviation", - kernels->mean); - psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_RMS, 0, "RMS stamp deviation", - kernels->rms); - - psMetadataAddS32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_STAMPS, 0, "Number of stamps", - kernels->numStamps); - psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_MEAN, 0, "Mean stamp deviation", - kernels->mean); - psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_RMS, 0, "RMS stamp deviation", - kernels->rms); - - psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_MEAN, 0, "Fractional Sigma of Residuals (Mean)", kernels->fResSigmaMean); - psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_STDEV, 0, "Fractional Sigma of Residuals (Stdev)", kernels->fResSigmaStdev); - psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_MEAN, 0, "Fractional Residual Flux (Mean, R > 2 pix)", kernels->fResOuterMean); - psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_STDEV, 0, "Fractional Residual Flux (Stdev, R > 2 pix)", kernels->fResOuterStdev); - psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_MEAN, 0, "Fractional Residual Flux (Mean, R > 0 pix)", kernels->fResTotalMean); - psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_STDEV, 0, "Fractional Residual Flux (Stdev, R > 0 pix)", kernels->fResTotalStdev); - - psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_MEAN, 0, "Fractional Sigma of Residuals (Mean)", kernels->fResSigmaMean); - psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_STDEV, 0, "Fractional Sigma of Residuals (Stdev)", kernels->fResSigmaStdev); - psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_MEAN, 0, "Fractional Residual Flux (Mean, R > 2 pix)", kernels->fResOuterMean); - psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_STDEV, 0, "Fractional Residual Flux (Stdev, R > 2 pix)", kernels->fResOuterStdev); - psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_MEAN, 0, "Fractional Residual Flux (Mean, R > 0 pix)", kernels->fResTotalMean); - psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_STDEV, 0, "Fractional Residual Flux (Stdev, R > 0 pix)", kernels->fResTotalStdev); + psMetadataAddS32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_STAMPS, PS_META_REPLACE, "Number of stamps", kernels->numStamps); + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_MEAN, PS_META_REPLACE, "Mean stamp deviation", kernels->mean); + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_RMS, PS_META_REPLACE, "RMS stamp deviation", kernels->rms); + psMetadataAddS32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_STAMPS, PS_META_REPLACE, "Number of stamps", kernels->numStamps); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_MEAN, PS_META_REPLACE, "Mean stamp deviation", kernels->mean); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_RMS, PS_META_REPLACE, "RMS stamp deviation", kernels->rms); + + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_MEAN, PS_META_REPLACE, "Fractional Sigma of Residuals (Mean)", kernels->fResSigmaMean); + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_STDEV, PS_META_REPLACE, "Fractional Sigma of Residuals (Stdev)", kernels->fResSigmaStdev); + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_MEAN, PS_META_REPLACE, "Fractional Residual Flux (Mean, R > 2 pix)", kernels->fResOuterMean); + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_STDEV, PS_META_REPLACE, "Fractional Residual Flux (Stdev, R > 2 pix)", kernels->fResOuterStdev); + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_MEAN, PS_META_REPLACE, "Fractional Residual Flux (Mean, R > 0 pix)", kernels->fResTotalMean); + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_STDEV, PS_META_REPLACE, "Fractional Residual Flux (Stdev, R > 0 pix)", kernels->fResTotalStdev); + + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_MEAN, PS_META_REPLACE, "Fractional Sigma of Residuals (Mean)", kernels->fResSigmaMean); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_STDEV, PS_META_REPLACE, "Fractional Sigma of Residuals (Stdev)", kernels->fResSigmaStdev); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_MEAN, PS_META_REPLACE, "Fractional Residual Flux (Mean, R > 2 pix)", kernels->fResOuterMean); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_STDEV, PS_META_REPLACE, "Fractional Residual Flux (Stdev, R > 2 pix)", kernels->fResOuterStdev); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_MEAN, PS_META_REPLACE, "Fractional Residual Flux (Mean, R > 0 pix)", kernels->fResTotalMean); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_STDEV, PS_META_REPLACE, "Fractional Residual Flux (Stdev, R > 0 pix)", kernels->fResTotalStdev); }