Index: trunk/psModules/src/imcombine/pmStack.c =================================================================== --- trunk/psModules/src/imcombine/pmStack.c (revision 21474) +++ trunk/psModules/src/imcombine/pmStack.c (revision 21476) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.47 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-06 02:31:25 $ + * @version $Revision: 1.48 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-13 23:52:14 $ * Copyright 2004-2007 Institute for Astronomy, University of Hawaii * @@ -30,7 +30,10 @@ #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 ADD_VARIANCE // Allow additional variance (besides variance factor)? +#define ADD_VARIANCE // Allow additional variance (besides variance factor)? #define NUM_DIRECT_STDEV 5 // For less than this number of values, measure stdev directly + +//#define TESTING // Enable test output +//#define TEST_X 2318 // x coordinate to examine +//#define TEST_Y 2306 // y coordinate to examine @@ -190,4 +193,96 @@ } +#if 0 +// Return the weighted median for the pixels +// This does not appear to produce as clean images as the weighted Olympic mean +static float combinationWeightedMedian(const psVector *values, // Values to combine + const psVector *weights, // Weights to combine + const psVector *masks, // Mask to apply + psVector *sortBuffer // Buffer for sorting + ) +{ + double sumWeight = 0.0; // Sum of weights + for (int j = 0; j < values->n; j++) { + if (masks->data.PS_TYPE_VECTOR_MASK_DATA[j]) { + continue; + } + sumWeight += weights->data.F32[j]; + } + + sortBuffer = psVectorSortIndex(sortBuffer, values); + double target = sumWeight / 2.0; // Target weight + + int dominant = -1; // Index of dominant value, if any + double cumulativeWeight = 0.0; // Sum of weights + for (int j = 0; j < values->n; j++) { + if (masks->data.PS_TYPE_VECTOR_MASK_DATA[j]) { + continue; + } + int index = sortBuffer->data.S32[j]; // Index of value of interest + float weight = weights->data.F32[index]; // Weight for value of interest + if (weight >= target) { + // Get the weighted median of the rest + dominant = index; + sumWeight -= weight; + target = sumWeight / 2.0; + continue; + } + cumulativeWeight += weight; + if (cumulativeWeight >= target) { + float median = values->data.F32[index]; // Weighted median median + if (dominant != -1) { + // In the case that a single value contains a disproportionate weight compared to the rest, + // we use a weighted mean between that dominant value and the weighted median of the rest. + return (values->data.F32[dominant] * weights->data.F32[dominant] + median * sumWeight) / + (weights->data.F32[dominant] + sumWeight); + } + return median; + } + } + + return NAN; +} +#endif + +// Return the weighted Olympic mean for the pixels +static float combinationWeightedOlympic(const psVector *values, // Values to combine + const psVector *weights, // Weights to combine + const psVector *masks, // Mask to apply + float frac, // Fraction to discard + psVector *sortBuffer // Buffer for sorting + ) +{ + int numGood = 0; // Number of good values + for (int i = 0; i < values->n; i++) { + if (masks->data.PS_TYPE_VECTOR_MASK_DATA[i]) { + continue; + } + numGood++; + } + + int numBad = frac * numGood + 0.5; // Number of bad values + int low = numBad / 2, high = low + numGood; // Indices (modulo masked pixels) delimiting range of interest + + sortBuffer = psVectorSortIndex(sortBuffer, values); + + double sumValues = 0.0, sumWeight = 0.0; // Sums for weighted mean + for (int i = 0, j = 0; i < values->n; i++) { + int index = sortBuffer->data.S32[i]; // Index of interest + if (masks->data.PS_TYPE_VECTOR_MASK_DATA[index]) { + continue; + } + j++; + if (j > high) { + break; + } + if (j <= low) { + continue; + } + sumValues += values->data.F32[index] * weights->data.F32[index]; + sumWeight += weights->data.F32[index]; + } + + return sumValues / sumWeight; +} // Mark a pixel for inspection @@ -213,5 +308,5 @@ const psArray *inputs, // Stack data const psVector *weights, // Global (single value) weights for data, or NULL - const psVector *varFactors, // Variance factors for data + const psVector *addVariance, // Additional variance for data const psVector *reject, // Indices of pixels to reject, or NULL int x, int y, // Coordinates of interest; frame of output image @@ -221,4 +316,5 @@ float rej, // Number of standard deviations at which to reject float sys, // Relative systematic error + float discard,// Fraction of values to discard (Olympic weighted mean) bool useVariance, // Use variance for rejection when combining? bool safe, // Combine safely? @@ -232,5 +328,5 @@ assert(numIter >= 0); assert(buffer); - assert(varFactors); + assert(addVariance); assert((useVariance && variance) || !useVariance); @@ -267,5 +363,5 @@ pixelData->data.F32[num] = image->data.F32[yIn][xIn]; if (variance) { - pixelVariances->data.F32[num] = variance->data.F32[yIn][xIn]; + pixelVariances->data.F32[num] = variance->data.F32[yIn][xIn] * addVariance->data.F32[i]; } pixelWeights->data.F32[num] = data->weight; @@ -280,4 +376,13 @@ pixelSources->n = num; +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + for (int i = 0; i < num; i++) { + fprintf(stderr, "Input %d (%" PRIu16 "): %f %f %f\n", + i, pixelSources->data.U16[i], pixelData->data.F32[i], pixelVariances->data.F32[i], + pixelWeights->data.F32[i]); + } + } +#endif // The sensible thing to do varies according to how many good pixels there are. @@ -288,8 +393,18 @@ case 0: // Nothing to combine: it's bad +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "No inputs to combine, pixel is bad.\n"); + } +#endif break; case 1: { // Accept the single pixel unless we have to be safe if (!safe) { +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Single input to combine, safety off.\n"); + } +#endif imageValue = pixelData->data.F32[0]; if (variance) { @@ -298,4 +413,9 @@ maskValue = 0; } +#ifdef TESTING + else if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Single input to combine, safety on, pixel is bad.\n"); + } +#endif break; } @@ -309,4 +429,10 @@ varianceValue = var; maskValue = 0; +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Two inputs to combine using variance/unsafe --> %f %f\n", + mean, var); + } +#endif } } @@ -316,18 +442,19 @@ float var1 = pixelVariances->data.F32[0]; // Variance of first float var2 = pixelVariances->data.F32[1]; // Variance of second -#ifdef VARIANCE_FACTORS - // Variance factor contributes to the rejection level - var1 *= varFactors->data.F32[pixelSources->data.U16[0]]; - var2 *= varFactors->data.F32[pixelSources->data.U16[1]]; -#endif // Systematic error contributes to the rejection level var1 += PS_SQR(sys * pixelData->data.F32[0]); var2 += PS_SQR(sys * pixelData->data.F32[1]); - float sigma2 = var1 + var2; // + float sigma2 = var1 + var2; // Combined variance if (PS_SQR(diff) > PS_SQR(rej) * sigma2) { // Not consistent: mark both for inspection combineInspect(inputs, x, y, pixelSources->data.U16[0]); combineInspect(inputs, x, y, pixelSources->data.U16[1]); +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Both pixels marked for inspection (%f > %f x %f\n)", + diff, rej, sqrtf(sigma2)); + } +#endif } } @@ -344,4 +471,9 @@ varianceValue = var; maskValue = 0; +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Combined inputs: %f %f\n", mean, var); + } +#endif // Prepare for clipping iteration @@ -353,12 +485,21 @@ // Using squared rejection limit because it's cheaper than sqrts float rej2 = PS_SQR(rej); // Rejection level squared + double sumWeights = 0.0; + for (int i = 0; i < num; i++) { + sumWeights += pixelWeights->data.F32[i]; + } for (int i = 0; i < num; i++) { // Systematic error contributes to the rejection level - pixelVariances->data.F32[i] += PS_SQR(sys * pixelData->data.F32[i]); -#ifdef VARIANCE_FACTORS - // Variance factor contributes to the rejection level - pixelVariances->data.F32[i] *= varFactors->data.F32[pixelSources->data.U16[i]]; -#endif - pixelVariances->data.F32[i] *= rej2; + float sysVar = PS_SQR(sys * pixelData->data.F32[i]); +#ifdef TESTING + // Correct variance for comparison against weighted mean including itself + float compare = 1.0 - pixelWeights->data.F32[i] / sumWeights; + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Variance %d (%d): %f %f %f\n", i, pixelSources->data.U16[i], + pixelVariances->data.F32[i], sysVar, compare); + } +#endif + + pixelVariances->data.F32[i] = rej2 * (pixelVariances->data.F32[i] + sysVar); } } @@ -371,13 +512,29 @@ for (int i = 0; i < numIter && numClipped > 0 && num - totalClipped > 2; i++) { numClipped = 0; - float median, stdev; // Median and stdev of the combination, for rejection - - if (!combinationMedianStdev(&median, useVariance ? NULL : &stdev, - pixelData, pixelMasks, sort)) { - psWarning("Bad median/stdev at %d,%d", x, y); - break; + float median = NAN; // Middle of distribution + float limit = NAN; // Rejection limit + if (!useVariance) { + float stdev; // Median and stdev of the combination, for rejection + if (!combinationMedianStdev(&median, useVariance ? NULL : &stdev, + pixelData, pixelMasks, sort)) { + psWarning("Bad median/stdev at %d,%d", x, y); + break; + } + limit = rej * stdev; +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Rejection limit: %f\n", limit); + } +#endif + } else { + median = combinationWeightedOlympic(pixelData, pixelWeights, pixelMasks, discard, sort); } - float limit = rej * stdev; // Rejection limit, when rejecting based on data properties +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Median: %f\n", median); + } +#endif + // Mask a pixel for inspection @@ -395,10 +552,22 @@ if (useVariance) { // Comparing squares --- cheaper than lots of sqrts - // pixelVariances includes the variance factor and the rejection limit, from above + // pixelVariances includes the rejection limit, from above if (PS_SQR(diff) > pixelVariances->data.F32[j]) { MASK_PIXEL_FOR_INSPECTION(); +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Rejecting input %d based on variance: %f > %f\n", + j, diff, sqrtf(pixelVariances->data.F32[j])); + } +#endif } } else if (fabsf(diff) > limit) { MASK_PIXEL_FOR_INSPECTION(); +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Rejecting input %d based on distribution: %f > %f\n", + j, diff, limit); + } +#endif } } @@ -564,5 +733,5 @@ /// Stack input images bool pmStackCombine(pmReadout *combined, psArray *input, psImageMaskType maskVal, psImageMaskType bad, - int kernelSize, int numIter, float rej, float sys, + int kernelSize, int numIter, float rej, float sys, float discard, bool entire, bool useVariance, bool safe) { @@ -596,5 +765,5 @@ } - psVector *varFactors = psVectorAlloc(num, PS_TYPE_F32); // Variance factors for each image + psVector *addVariance = psVectorAlloc(num, PS_TYPE_F32); // Additional variance for each image psVector *weights = psVectorAlloc(num, PS_TYPE_F32); // Relative weighting for each image psArray *stack = psArrayAlloc(num); // Stack of readouts @@ -606,10 +775,10 @@ } weights->data.F32[i] = data->weight; - stack->data[i] = psMemIncrRefCounter(data->readout); - varFactors->data.F32[i] = data->readout->covariance ? - psImageCovarianceFactor(data->readout->covariance) : 1.0; -#if 0 + pmReadout *ro = data->readout; // Readout of interest + stack->data[i] = psMemIncrRefCounter(ro); + addVariance->data.F32[i] = ro->covariance ? psImageCovarianceFactor(ro->covariance) : 1.0; +#ifdef ADD_VARIANCE if (isfinite(data->addVariance)) { - varFactors->data.F32[i] *= data->addVariance; + addVariance->data.F32[i] *= data->addVariance; } #endif @@ -667,6 +836,7 @@ psVector *reject = pixelMapQuery(pixelMap, minInputCols, minInputRows, x, y); // Reject these images - combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, varFactors, - reject, x, y, maskVal, bad, numIter, rej, sys, useVariance, safe, buffer); + combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, + addVariance, reject, x, y, maskVal, bad, numIter, rej, sys, discard, + useVariance, safe, buffer); } } @@ -681,6 +851,7 @@ psVector *reject = pixelMapQuery(pixelMap, minInputCols, minInputRows, x, y); // Reject these images - combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, varFactors, - reject, x, y, maskVal, bad, numIter, rej, sys, useVariance, safe, buffer); + combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, + addVariance, reject, x, y, maskVal, bad, numIter, rej, sys, discard, + useVariance, safe, buffer); } } @@ -708,6 +879,7 @@ for (int y = minInputRows; y < maxInputRows; y++) { for (int x = minInputCols; x < maxInputCols; x++) { - combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, varFactors, - NULL, x, y, maskVal, bad, numIter, rej, sys, useVariance, safe, buffer); + combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, + addVariance, NULL, x, y, maskVal, bad, numIter, rej, sys, discard, + useVariance, safe, buffer); } } Index: trunk/psModules/src/imcombine/pmStack.h =================================================================== --- trunk/psModules/src/imcombine/pmStack.h (revision 21474) +++ trunk/psModules/src/imcombine/pmStack.h (revision 21476) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-01-27 06:39:38 $ + * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-13 23:52:14 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -49,4 +49,5 @@ float rej, ///< Rejection limit (standard deviations) float sys, ///< Relative systematic error + float discard, ///< Fraction of values to discard for Olympic weighted mean bool entire, ///< Combine entire image even if rejection lists provided? bool useVariance, ///< Use variance values for rejection?