Index: /trunk/ppStack/src/ppStack.h =================================================================== --- /trunk/ppStack/src/ppStack.h (revision 19531) +++ /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 19531) +++ /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 19531) +++ /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 19531) +++ /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 19531) +++ /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); Index: /trunk/psModules/src/imcombine/pmStack.c =================================================================== --- /trunk/psModules/src/imcombine/pmStack.c (revision 19531) +++ /trunk/psModules/src/imcombine/pmStack.c (revision 19532) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.40 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-09-11 20:52:16 $ + * @version $Revision: 1.41 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-09-12 04:12:42 $ * Copyright 2004-2007 Institute for Astronomy, University of Hawaii * @@ -30,5 +30,5 @@ #define PIXEL_LIST_BUFFER 100 // Number of entries to add to pixel list at a time #define PIXEL_MAP_BUFFER 2 // Number of entries to add to pixel map at a time -//#define VARIANCE_FACTORS // Use variance factors when calculating the variances? +#define VARIANCE_FACTORS // Use variance factors when calculating the variances? #define NUM_DIRECT_STDEV 5 // For less than this number of values, measure stdev directly @@ -312,5 +312,5 @@ // Use variance to check that the two are consistent float diff = pixelData->data.F32[0] - pixelData->data.F32[1]; -#if VARIANCE_FACTORS +#ifdef VARIANCE_FACTORS float sigma2 = pixelVariances->data.F32[0] * varFactors->data.F32[pixelSources->data.U16[0]] + pixelVariances->data.F32[1] * varFactors->data.F32[pixelSources->data.U16[1]]; @@ -346,5 +346,5 @@ float rej2 = PS_SQR(rej); // Rejection level squared for (int i = 0; i < num; i++) { -#if VARIANCE_FACTORS +#ifdef VARIANCE_FACTORS pixelVariances->data.F32[i] *= rej2 * varFactors->data.F32[pixelSources->data.U16[i]]; #else @@ -538,5 +538,5 @@ /// Constructor -pmStackData *pmStackDataAlloc(pmReadout *readout, float weight) +pmStackData *pmStackDataAlloc(pmReadout *readout, float weight, float addVariance) { pmStackData *data = psAlloc(sizeof(pmStackData)); // Stack data, to return @@ -547,4 +547,5 @@ data->inspect = NULL; data->weight = weight; + data->addVariance = addVariance; return data; @@ -602,4 +603,7 @@ } varFactors->data.F32[i] = vf; + if (isfinite(data->addVariance)) { + varFactors->data.F32[i] *= data->addVariance; + } if (!haveRejects && !data->inspect) { data->inspect = psPixelsAllocEmpty(PIXEL_LIST_BUFFER); Index: /trunk/psModules/src/imcombine/pmStack.h =================================================================== --- /trunk/psModules/src/imcombine/pmStack.h (revision 19531) +++ /trunk/psModules/src/imcombine/pmStack.h (revision 19532) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-03-17 21:38:43 $ + * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-09-12 04:12:42 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -31,9 +31,11 @@ psPixels *inspect; ///< Pixels to inspect float weight; ///< Relative weighting for image + float addVariance; ///< Additional variance when rejecting } pmStackData; /// Constructor pmStackData *pmStackDataAlloc(pmReadout *readout, ///< Warped readout (sky cell) - float weight ///< Weight to apply + float weight, ///< Weight to apply + float addVariance ///< Additional variance when rejecting ); Index: /trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtraction.c (revision 19531) +++ /trunk/psModules/src/imcombine/pmSubtraction.c (revision 19532) @@ -665,8 +665,8 @@ -int pmSubtractionRejectStamps(float *rmsPtr, int *numPtr, pmSubtractionStampList *stamps, - const psVector *deviations, psImage *subMask, float sigmaRej, - int footprint) -{ +int pmSubtractionRejectStamps(pmSubtractionKernels *kernels, pmSubtractionStampList *stamps, + const psVector *deviations, psImage *subMask, float sigmaRej, int footprint) +{ + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, -1); PS_ASSERT_VECTOR_NON_NULL(deviations, -1); @@ -678,4 +678,8 @@ // the distribution is actually something like a chi^2 or Student's t, both of which become Gaussian-like // with large N. Therefore, let's just treat this as a Gaussian distribution. + + kernels->mean = NAN; + kernels->rms = NAN; + kernels->numStamps = -1; int numStamps = 0; // Number of used stamps @@ -721,15 +725,13 @@ psTrace("psModules.imcombine", 1, "RMS deviation: %f\n", rms); - if (rmsPtr) { - *rmsPtr = rms; - } - if (numPtr) { - *numPtr = numStamps; - } - - psLogMsg("psModules.imcombine", PS_LOG_INFO, "RMS deviation from %d stamps: %lf", numStamps, rms); + kernels->mean = mean; + kernels->rms = rms; + kernels->numStamps = numStamps; + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Mean deviation from %d stamps: %lf +/- %lf", + numStamps, mean, rms); if (!isfinite(sigmaRej) || sigmaRej <= 0.0) { - // User just wanted to calculate and record the RMS for posterity + // User just wanted to calculate and record the deviation for posterity return 0; } Index: /trunk/psModules/src/imcombine/pmSubtraction.h =================================================================== --- /trunk/psModules/src/imcombine/pmSubtraction.h (revision 19531) +++ /trunk/psModules/src/imcombine/pmSubtraction.h (revision 19532) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.31 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-08-30 02:28:07 $ + * @version $Revision: 1.32 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-09-12 04:12:42 $ * Copyright 2004-207 Institute for Astronomy, University of Hawaii */ @@ -64,6 +64,5 @@ /// Reject stamps -int pmSubtractionRejectStamps(float *rms, ///< RMS deviation, to return - int *num, ///< Number of good stamps, to return +int pmSubtractionRejectStamps(pmSubtractionKernels *kernels, ///< Kernel parameters to update pmSubtractionStampList *stamps, ///< Stamps const psVector *deviations, ///< Deviations for each stamp Index: /trunk/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 19531) +++ /trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 19532) @@ -1085,5 +1085,5 @@ } } - deviations->data.F32[i] = sqrtf(devNorm * deviation); + deviations->data.F32[i] = devNorm * deviation; psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n", i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]); Index: /trunk/psModules/src/imcombine/pmSubtractionIO.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionIO.c (revision 19531) +++ /trunk/psModules/src/imcombine/pmSubtractionIO.c (revision 19532) @@ -32,4 +32,7 @@ #define NAME_SOL1 "SOLUTION_1" // Solution for convolving image 1 #define NAME_SOL2 "SOLUTION_2" // Solution for convolving image 2 +#define NAME_MEAN "MEAN" // Mean of chi^2 from stamps +#define NAME_RMS "RMS" // RMS of chi^2 from stamps +#define NAME_NUMSTAMPS "NUMSTAMPS" // Number of good stamps #define FALLBACK_EXTNAME "SUBTRACTION_KERNEL" // Fallback extension name @@ -159,4 +162,7 @@ psMetadataAddVector(row, PS_LIST_TAIL, NAME_SOL2, 0, "Solution vector 2", kernel->solution2); } + psMetadataAddF32(row, PS_LIST_TAIL, NAME_MEAN, 0, "Mean of chi^2 from stamps", kernel->mean); + psMetadataAddF32(row, PS_LIST_TAIL, NAME_RMS, 0, "RMS of chi^2 from stamps", kernel->rms); + psMetadataAddF32(row, PS_LIST_TAIL, NAME_NUMSTAMPS, 0, "Number of good stamps", kernel->numStamps); } psFree(regions); @@ -167,6 +173,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.3 $"); - psString cvsDate = psStringCopy("$Date: 2008-08-27 02:54:44 $"); + psString cvsRev = psStringCopy("$Revision: 1.4 $"); + psString cvsDate = psStringCopy("$Date: 2008-09-12 04:12:42 $"); psStringSubstitute(&cvsFile, NULL, "RCSfile: "); psStringSubstitute(&cvsRev, NULL, "Revision: "); @@ -338,7 +344,14 @@ TABLE_LOOKUP(int, S32, numRows, NAME_ROWS); + TABLE_LOOKUP(float, F32, mean, NAME_MEAN); + TABLE_LOOKUP(float, F32, rms, NAME_RMS); + TABLE_LOOKUP(int, S32, numStamps, NAME_NUMSTAMPS); + pmSubtractionKernels *kernels = pmSubtractionKernelsFromDescription(description, bg, mode); kernels->numCols = numCols; kernels->numRows = numRows; + kernels->mean = mean; + kernels->rms = rms; + kernels->numStamps = numStamps; bool mdok; // Status of MD lookup Index: /trunk/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 19531) +++ /trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 19532) @@ -42,4 +42,7 @@ int numCols, numRows; ///< Size of image (for normalisation), or zero to use image provided psVector *solution1, *solution2; ///< Solution for the PSF matching + // Quality information + float mean, rms; ///< Mean and RMS of chi^2 from stamps + int numStamps; ///< Number of good stamps } pmSubtractionKernels; Index: /trunk/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 19531) +++ /trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 19532) @@ -370,6 +370,4 @@ memCheck("kernels"); - float rmsStamps = NAN; // RMS for stamps - int numStamps = 0; // Number of good stamps int numRejected = -1; // Number of rejected stamps in each iteration for (int k = 0; k < iter && numRejected != 0; k++) { @@ -407,6 +405,5 @@ psTrace("psModules.imcombine", 3, "Rejecting stamps...\n"); - numRejected = pmSubtractionRejectStamps(&rmsStamps, &numStamps, stamps, deviations, - subMask, rej, footprint); + numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej, footprint); if (numRejected < 0) { psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps."); @@ -430,23 +427,9 @@ goto MATCH_ERROR; } - pmSubtractionRejectStamps(&rmsStamps, &numStamps, stamps, deviations, - subMask, NAN, footprint); + pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN, footprint); psFree(deviations); } psFree(stamps); stamps = NULL; - - if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_DUAL) { - psMetadataAddS32(conv1->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_STAMPS_NUM, - PS_META_DUPLICATE_OK, "Number of good stamps", numStamps); - psMetadataAddF32(conv1->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_STAMPS_RMS, - PS_META_DUPLICATE_OK, "RMS deviation of stamps", rmsStamps); - } - if (subMode == PM_SUBTRACTION_MODE_2 || subMode == PM_SUBTRACTION_MODE_DUAL) { - psMetadataAddS32(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_STAMPS_NUM, - PS_META_DUPLICATE_OK, "Number of good stamps", numStamps); - psMetadataAddF32(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_STAMPS_RMS, - PS_META_DUPLICATE_OK, "RMS deviation of stamps", rmsStamps); - } memCheck("solution"); Index: /trunk/psModules/src/imcombine/pmSubtractionMatch.h =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionMatch.h (revision 19531) +++ /trunk/psModules/src/imcombine/pmSubtractionMatch.h (revision 19532) @@ -9,10 +9,12 @@ #include -#define PM_SUBTRACTION_ANALYSIS_KERNEL "SUBTRACTION.KERNEL" // Name of kernel in the analysis metadata -#define PM_SUBTRACTION_ANALYSIS_MODE "SUBTRACTION.MODE" // Name of subtraction mode in the analysis metadata -#define PM_SUBTRACTION_ANALYSIS_REGION "SUBTRACTION.REGION" // Name of subtraction region in the analysis MD -#define PM_SUBTRACTION_ANALYSIS_STAMPS_RMS "SUBTRACTION.RMS" // Name of stamp rms in the analysis metadata -#define PM_SUBTRACTION_ANALYSIS_STAMPS_NUM "SUBTRACTION.NUM"// Name of the number of stamps in the analysis MD -#define PM_SUBTRACTION_ANALYSIS_VARFACTOR "SUBTRACTION.VARFACTOR"// Name of variance factor in the analysis MD +// 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_STAMPS_MEAN "SUBTRACTION.MEAN" // Stamp mean deviation +#define PM_SUBTRACTION_ANALYSIS_STAMPS_RMS "SUBTRACTION.RMS" // Stamp RMS deviation +#define PM_SUBTRACTION_ANALYSIS_STAMPS_NUM "SUBTRACTION.NUM" // Number of stamps +#define PM_SUBTRACTION_ANALYSIS_VARFACTOR "SUBTRACTION.VARFACTOR"// Variance factor