Index: trunk/psModules/src/imcombine/pmStack.c =================================================================== --- trunk/psModules/src/imcombine/pmStack.c (revision 16628) +++ trunk/psModules/src/imcombine/pmStack.c (revision 16629) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-02-23 03:01:45 $ + * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-02-23 03:19:34 $ * Copyright 2004-2007 Institute for Astronomy, University of Hawaii * @@ -86,28 +86,42 @@ -// Generate a mean value for the combination +// Determine a mean value and variance for the combination // Not using psVectorStats because it assumes that the "weights" are errors, and weights by 1/error^2 -static bool combinationMean(float *mean, // Mean value, to return - const psVector *values, // Values to combine - const psVector *weights, // Weights to apply - const psVector *masks, // Mask to apply - psMaskType maskVal // Value to mask - ) -{ - assert(values && weights && masks); +static bool combinationMeanVariance(float *mean, // Mean value, to return + float *var, // Variance value, to return + const psVector *values, // Values to combine + const psVector *variances, // Pixel variances to combine + const psVector *weights // Weights to apply + ) +{ + assert(mean); + assert(values && weights); assert(values->n == weights->n); - assert(values->n == masks->n); + assert((var && variances) || !var); + assert(!variances || variances->n == values->n); assert(values->type.type == PS_TYPE_F32); + assert(!values || values->type.type == PS_TYPE_F32); assert(weights->type.type == PS_TYPE_F32); - assert(masks->type.type == PS_TYPE_MASK); - - float sumValueWeight = 0.0; // Sum of the value multiplied by the weight - float sumWeight = 0.0; // Sum of the weights + + // We're not using the input pixel variances to generate a weighted average for the pixel flux (because + // that introduces systematic biases), so the variance of the output pixel value should simply be: + // simga^2 = sum(weight_i^2 * sigma_i^2) / (sum(weight_i))^2 + // This reduces, when the weights are all identically unity, to: + // variance_combination = sum(variance_i) / N^2 + // and if the variances are all equal: + // variance_combination = variance_individual / N + // which makes sense --- the standard deviation of the combination is reduced by a factor of sqrt(N). + + float sumValueWeight = 0.0; // Sum of the value multiplied by the weight + float sumVarianceWeight = 0.0; // Sum of the pixel variances multiplied by the global weights + float sumWeight = 0.0; // Sum of the image weights for (int i = 0; i < values->n; i++) { - if (!(masks->data.PS_TYPE_MASK_DATA[i] & maskVal)) { - sumValueWeight += values->data.F32[i] * weights->data.F32[i]; - sumWeight += weights->data.F32[i]; - } - } + sumValueWeight += values->data.F32[i] * weights->data.F32[i]; + sumWeight += weights->data.F32[i]; + if (variances) { + sumVarianceWeight += variances->data.F32[i] * PS_SQR(weights->data.F32[i]); + } + } + if (sumWeight <= 0) { return false; @@ -115,4 +129,7 @@ *mean = sumValueWeight / sumWeight; + if (var) { + *var = sumVarianceWeight / PS_SQR(sumWeight); + } return true; } @@ -124,12 +141,11 @@ const psVector *values, // Values to combine const psVector *masks, // Mask to apply - psMaskType maskVal, // Value to mask psVector *sortBuffer // Buffer for sorting ) { - assert(values && masks); - assert(values->n == masks->n); + assert(values); + assert(!masks || values->n == masks->n); assert(values->type.type == PS_TYPE_F32); - assert(masks->type.type == PS_TYPE_MASK); + assert(!masks || masks->type.type == PS_TYPE_MASK); assert(sortBuffer && sortBuffer->nalloc >= values->n && sortBuffer->type.type == PS_TYPE_F32); @@ -137,5 +153,5 @@ int num = 0; // Number of valid values for (int i = 0; i < values->n; i++) { - if (!(masks->data.PS_TYPE_MASK_DATA[i])) { + if (!masks || !masks->data.PS_TYPE_MASK_DATA[i]) { sortBuffer->data.F32[num++] = values->data.F32[i]; } @@ -174,46 +190,4 @@ } - return true; -} - -// Generate a variance value for the combination -static bool combinationVariance(float *variance, // Variance value, to return - const psVector *variances, // Pixel variances to combine - const psVector *weights, // Image weights to apply - const psVector *masks, // Mask to apply - psMaskType maskVal // Value to mask - ) -{ - assert(variances && weights && masks); - assert(variances->n == weights->n); - assert(variances->n == masks->n); - assert(variances->type.type == PS_TYPE_F32); - assert(weights->type.type == PS_TYPE_F32); - assert(masks->type.type == PS_TYPE_MASK); - - // Get the variance in the combination. We're not using the input pixel variances to generate a - // weighted average for the pixel flux (because that introduces systematic biases), so the variance - // of the output pixel value should simply be: - // simga^2 = sum(weight_i^2 * sigma_i^2) / (sum(weight_i))^2 - // This reduces, when the weights are all identically unity, to: - // variance_combination = sum(variance_i) / N^2 - // and if the variances are all equal: - // variance_combination = variance_individual / N - // which makes sense --- the standard deviation of the combination is reduced by a factor of sqrt(N). - - float sumWeights = 0.0; // Sum of the global weights - float sumVarianceWeights = 0.0; // Sum of the pixel variances multiplied by the global weights - for (int i = 0; i < variances->n; i++) { - if (!(masks->data.PS_TYPE_MASK_DATA[i] & maskVal)) { - sumVarianceWeights += variances->data.F32[i] * PS_SQR(weights->data.F32[i]); - sumWeights += weights->data.F32[i]; - } - } - - if (sumWeights <= 0) { - return false; - } - - *variance = sumVarianceWeights / PS_SQR(sumWeights); return true; } @@ -263,5 +237,5 @@ psVector *pixelData = buffer->pixels; // Values for the pixel of interest psVector *pixelMasks = buffer->masks; // Masks for the pixel of interest - psVector *pixelVariances = buffer->variances; // Variances for the pixel of interest + psVector *pixelVariances = variance ? buffer->variances : NULL; // Variances for the pixel of interest psVector *pixelWeights = buffer->weights; // Image weights for the pixel of interest psVector *pixelSources = buffer->sources; // Sources for the pixel of interest @@ -331,6 +305,5 @@ if (useVariance || !safe) { float mean, var; // Mean and variance from combination - if (combinationMean(&mean, pixelData, pixelWeights, NULL, maskVal) && - (!variance || combinationVariance(&var, pixelVariances, pixelWeights, NULL, maskVal))) { + if (combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) { imageValue = mean; varianceValue = var; @@ -354,6 +327,5 @@ // applied in the first pass might later be accepted. float mean, var; // Mean and variance of the combination - if (!combinationMean(&mean, pixelData, pixelWeights, pixelMasks, maskVal) || - (variance && !combinationVariance(&var, pixelVariances, pixelWeights, pixelMasks, maskVal))) { + if (!combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) { break; } @@ -382,6 +354,5 @@ float median, stdev; // Median and stdev of the combination, for rejection - if (!combinationMedianStdev(&median, useVariance ? NULL : &stdev, pixelData, pixelMasks, - maskVal, sort)) { + if (!combinationMedianStdev(&median, useVariance ? NULL : &stdev, pixelData, pixelMasks, sort)) { psWarning("Bad median/stdev at %d,%d", x, y); break;