Index: trunk/ppStack/src/ppStack.h =================================================================== --- trunk/ppStack/src/ppStack.h (revision 19505) +++ trunk/ppStack/src/ppStack.h (revision 19532) @@ -89,5 +89,6 @@ const psArray *readouts, // Input readouts const psArray *regions, // Array with array of regions used in each PSF match - const psArray *kernels // Array with array of kernels used in each PSF match + const psArray *kernels, // Array with array of kernels used in each PSF match + const psVector *addVariance // Additional variance for rejection ); @@ -128,7 +129,8 @@ /// Convolve image to match specified seeing -bool ppStackMatch(pmReadout *readout, // Readout to be convolved; replaced with output - psArray **regions, // Array of regions used in each PSF matching, returned - psArray **kernels, // Array of kernels used in each PSF matching, returned +bool ppStackMatch(pmReadout *readout, // Readout to be convolved; replaced with output + psArray **regions, // Array of regions used in each PSF matching, returned + psArray **kernels, // Array of kernels used in each PSF matching, returned + float *chi2, // Chi^2 from the stamps const psArray *sources, // Array of sources const pmPSF *psf, // Target PSF Index: trunk/ppStack/src/ppStackLoop.c =================================================================== --- trunk/ppStack/src/ppStackLoop.c (revision 19505) +++ trunk/ppStack/src/ppStackLoop.c (revision 19532) @@ -390,4 +390,6 @@ psVector *inputMask = psVectorAlloc(num, PS_TYPE_U8); // Mask for inputs psVectorInit(inputMask, 0); + psVector *matchChi2 = psVectorAlloc(num, PS_TYPE_F32); // chi^2 for stamps when matching + psVectorInit(matchChi2, NAN); for (int i = 0; i < num; i++) { psTrace("ppStack", 2, "Convolving input %d of %d to target PSF....\n", i, num); @@ -402,4 +404,5 @@ psFree(rng); psFree(inputMask); + psFree(matchChi2); return false; } @@ -419,4 +422,5 @@ psFree(rng); psFree(inputMask); + psFree(matchChi2); return false; } @@ -425,5 +429,6 @@ psArray *regions = NULL, *kernels = NULL; // Regions and kernels used in subtraction psArray *sources = haveSources ? globalSources : indSources->data[i]; // Sources for matching - if (!ppStackMatch(readout, ®ions, &kernels, sources, targetPSF, rng, config)) { + if (!ppStackMatch(readout, ®ions, &kernels, &matchChi2->data.F32[i], + sources, targetPSF, rng, config)) { psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i); inputMask->data.U8[i] = PPSTACK_MASK_MATCH; @@ -459,4 +464,6 @@ psFree(cells); psFree(inputMask); + psFree(inputMask); + psFree(matchChi2); return false; } @@ -473,5 +480,5 @@ memDump("preinitial"); - // Stack the convolvd files + // Stack the convolved files psTrace("ppStack", 1, "Initial stack of convolved images....\n"); pmReadout *outRO = NULL; // Output readout @@ -489,4 +496,5 @@ psFree(stack); psFree(inputMask); + psFree(matchChi2); psFree(exptimes); psFree(cells); @@ -502,4 +510,5 @@ psFree(stack); psFree(inputMask); + psFree(matchChi2); psFree(exptimes); psFree(cells); @@ -518,4 +527,5 @@ psFree(stack); psFree(inputMask); + psFree(matchChi2); psFree(view); psFree(outRO); @@ -546,5 +556,5 @@ psFree(cellList); - psFree(cells); + psFree(cells); bool status; // Status of read @@ -559,4 +569,5 @@ psFree(stack); psFree(inputMask); + psFree(matchChi2); psFree(view); psFree(outRO); @@ -576,4 +587,5 @@ psArrayAdd(job->args, 1, subRegions); psArrayAdd(job->args, 1, subKernels); + psArrayAdd(job->args, 1, matchChi2); if (!psThreadJobAddPending(job)) { psFree(job); @@ -582,4 +594,5 @@ psFree(stack); psFree(inputMask); + psFree(matchChi2); psFree(view); psFree(outRO); @@ -596,4 +609,5 @@ psFree(stack); psFree(inputMask); + psFree(matchChi2); psFree(view); psFree(outRO); @@ -601,4 +615,5 @@ return false; } + psFree(matchChi2); // Harvest the jobs, gathering the inspection lists Index: trunk/ppStack/src/ppStackMatch.c =================================================================== --- trunk/ppStack/src/ppStackMatch.c (revision 19505) +++ trunk/ppStack/src/ppStackMatch.c (revision 19532) @@ -15,5 +15,5 @@ PM_SOURCE_MODE_CR_LIMIT) // Mask to apply to input sources -//#define TESTING // Enable debugging output +#define TESTING // Enable debugging output @@ -48,5 +48,5 @@ #endif -bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels, +bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels, float *chi2, const psArray *sources, const pmPSF *psf, psRandom *rng, const pmConfig *config) { @@ -111,5 +111,5 @@ PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels float vf = pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); // Variance factor - psMetadataItem *vfItem = psMetadataLookup(output->parent->concepts, "CELL.VARFACTOR"); + psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR"); if (!isfinite(vf)) { vf = 1.0; @@ -341,4 +341,25 @@ assert((*regions)->n == (*kernels)->n); + // Correct the variance for the chi^2 + { + *chi2 = 0.0; + int num = 0; // Number of measurements of chi^2 + psString regex = NULL; // Regular expression + psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL); + psMetadataIterator *iter = psMetadataIteratorAlloc(output->analysis, PS_LIST_HEAD, regex); // Iterator + psFree(regex); + psMetadataItem *item = NULL;// Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + assert(item->type == PS_DATA_UNKNOWN); + pmSubtractionKernels *kernels = item->data.V; // Convolution kernels + *chi2 += kernels->mean; + num++; + } + psFree(iter); + + float vf = psMetadataLookupF32(NULL, readout->parent->concepts, "CELL.VARFACTOR"); // Variance factor + *chi2 /= vf * num; + } + // Renormalise the variances if desired if (renorm) { Index: trunk/ppStack/src/ppStackReadout.c =================================================================== --- trunk/ppStack/src/ppStackReadout.c (revision 19505) +++ trunk/ppStack/src/ppStackReadout.c (revision 19532) @@ -22,6 +22,8 @@ psArray *subRegions = args->data[3]; // Regions for PSF-matching psArray *subKernels = args->data[4]; // Kernels for PSF-matching - - psArray *inspect = ppStackReadoutInitial(config, outRO, thread->readouts, subRegions, subKernels); + psVector *addVariance = args->data[5]; // Additional variance when rejecting + + psArray *inspect = ppStackReadoutInitial(config, outRO, thread->readouts, + subRegions, subKernels, addVariance); job->results = inspect; @@ -75,5 +77,5 @@ psArray *ppStackReadoutInitial(const pmConfig *config, pmReadout *outRO, const psArray *readouts, - const psArray *regions, const psArray *kernels) + const psArray *regions, const psArray *kernels, const psVector *addVariance) { assert(config); @@ -84,4 +86,5 @@ assert(readouts->n == regions->n); assert(regions->n == kernels->n); + assert(addVariance && addVariance->n == readouts->n && addVariance->type.type == PS_TYPE_F32); static int sectionNum = 0; // Section number; for debugging outputs @@ -127,5 +130,5 @@ } - stack->data[i] = pmStackDataAlloc(ro, weighting); + stack->data[i] = pmStackDataAlloc(ro, weighting, addVariance->data.F32[i]); } @@ -221,5 +224,5 @@ } - pmStackData *data = pmStackDataAlloc(ro, weighting); + pmStackData *data = pmStackDataAlloc(ro, weighting, NAN); data->reject = psMemIncrRefCounter(rejected->data[i]); stack->data[i] = data; Index: trunk/ppStack/src/ppStackThread.c =================================================================== --- trunk/ppStack/src/ppStackThread.c (revision 19505) +++ trunk/ppStack/src/ppStackThread.c (revision 19532) @@ -235,5 +235,5 @@ { - psThreadTask *task = psThreadTaskAlloc("PPSTACK_INITIAL_COMBINE", 5); + psThreadTask *task = psThreadTaskAlloc("PPSTACK_INITIAL_COMBINE", 6); task->function = &ppStackReadoutInitialThread; psThreadTaskAdd(task);