Index: trunk/psModules/src/imcombine/pmStack.c =================================================================== --- trunk/psModules/src/imcombine/pmStack.c (revision 19509) +++ 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 19509) +++ 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 19509) +++ 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 19509) +++ 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 19509) +++ 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 19509) +++ 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 19509) +++ 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 19509) +++ 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 19509) +++ 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