Changeset 16629
- Timestamp:
- Feb 22, 2008, 5:19:34 PM (18 years ago)
- File:
-
- 1 edited
-
trunk/psModules/src/imcombine/pmStack.c (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/imcombine/pmStack.c
r16628 r16629 8 8 * @author GLG, MHPCC 9 9 * 10 * @version $Revision: 1.2 3$ $Name: not supported by cvs2svn $11 * @date $Date: 2008-02-23 03: 01:45$10 * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ 11 * @date $Date: 2008-02-23 03:19:34 $ 12 12 * Copyright 2004-2007 Institute for Astronomy, University of Hawaii 13 13 * … … 86 86 87 87 88 // Generate a mean value for the combination88 // Determine a mean value and variance for the combination 89 89 // Not using psVectorStats because it assumes that the "weights" are errors, and weights by 1/error^2 90 static bool combinationMean(float *mean, // Mean value, to return 91 const psVector *values, // Values to combine 92 const psVector *weights, // Weights to apply 93 const psVector *masks, // Mask to apply 94 psMaskType maskVal // Value to mask 95 ) 96 { 97 assert(values && weights && masks); 90 static bool combinationMeanVariance(float *mean, // Mean value, to return 91 float *var, // Variance value, to return 92 const psVector *values, // Values to combine 93 const psVector *variances, // Pixel variances to combine 94 const psVector *weights // Weights to apply 95 ) 96 { 97 assert(mean); 98 assert(values && weights); 98 99 assert(values->n == weights->n); 99 assert(values->n == masks->n); 100 assert((var && variances) || !var); 101 assert(!variances || variances->n == values->n); 100 102 assert(values->type.type == PS_TYPE_F32); 103 assert(!values || values->type.type == PS_TYPE_F32); 101 104 assert(weights->type.type == PS_TYPE_F32); 102 assert(masks->type.type == PS_TYPE_MASK); 103 104 float sumValueWeight = 0.0; // Sum of the value multiplied by the weight 105 float sumWeight = 0.0; // Sum of the weights 105 106 // We're not using the input pixel variances to generate a weighted average for the pixel flux (because 107 // that introduces systematic biases), so the variance of the output pixel value should simply be: 108 // simga^2 = sum(weight_i^2 * sigma_i^2) / (sum(weight_i))^2 109 // This reduces, when the weights are all identically unity, to: 110 // variance_combination = sum(variance_i) / N^2 111 // and if the variances are all equal: 112 // variance_combination = variance_individual / N 113 // which makes sense --- the standard deviation of the combination is reduced by a factor of sqrt(N). 114 115 float sumValueWeight = 0.0; // Sum of the value multiplied by the weight 116 float sumVarianceWeight = 0.0; // Sum of the pixel variances multiplied by the global weights 117 float sumWeight = 0.0; // Sum of the image weights 106 118 for (int i = 0; i < values->n; i++) { 107 if (!(masks->data.PS_TYPE_MASK_DATA[i] & maskVal)) { 108 sumValueWeight += values->data.F32[i] * weights->data.F32[i]; 109 sumWeight += weights->data.F32[i]; 110 } 111 } 119 sumValueWeight += values->data.F32[i] * weights->data.F32[i]; 120 sumWeight += weights->data.F32[i]; 121 if (variances) { 122 sumVarianceWeight += variances->data.F32[i] * PS_SQR(weights->data.F32[i]); 123 } 124 } 125 112 126 if (sumWeight <= 0) { 113 127 return false; … … 115 129 116 130 *mean = sumValueWeight / sumWeight; 131 if (var) { 132 *var = sumVarianceWeight / PS_SQR(sumWeight); 133 } 117 134 return true; 118 135 } … … 124 141 const psVector *values, // Values to combine 125 142 const psVector *masks, // Mask to apply 126 psMaskType maskVal, // Value to mask127 143 psVector *sortBuffer // Buffer for sorting 128 144 ) 129 145 { 130 assert(values && masks);131 assert( values->n == masks->n);146 assert(values); 147 assert(!masks || values->n == masks->n); 132 148 assert(values->type.type == PS_TYPE_F32); 133 assert( masks->type.type == PS_TYPE_MASK);149 assert(!masks || masks->type.type == PS_TYPE_MASK); 134 150 assert(sortBuffer && sortBuffer->nalloc >= values->n && sortBuffer->type.type == PS_TYPE_F32); 135 151 … … 137 153 int num = 0; // Number of valid values 138 154 for (int i = 0; i < values->n; i++) { 139 if (! (masks->data.PS_TYPE_MASK_DATA[i])) {155 if (!masks || !masks->data.PS_TYPE_MASK_DATA[i]) { 140 156 sortBuffer->data.F32[num++] = values->data.F32[i]; 141 157 } … … 174 190 } 175 191 176 return true;177 }178 179 // Generate a variance value for the combination180 static bool combinationVariance(float *variance, // Variance value, to return181 const psVector *variances, // Pixel variances to combine182 const psVector *weights, // Image weights to apply183 const psVector *masks, // Mask to apply184 psMaskType maskVal // Value to mask185 )186 {187 assert(variances && weights && masks);188 assert(variances->n == weights->n);189 assert(variances->n == masks->n);190 assert(variances->type.type == PS_TYPE_F32);191 assert(weights->type.type == PS_TYPE_F32);192 assert(masks->type.type == PS_TYPE_MASK);193 194 // Get the variance in the combination. We're not using the input pixel variances to generate a195 // weighted average for the pixel flux (because that introduces systematic biases), so the variance196 // of the output pixel value should simply be:197 // simga^2 = sum(weight_i^2 * sigma_i^2) / (sum(weight_i))^2198 // This reduces, when the weights are all identically unity, to:199 // variance_combination = sum(variance_i) / N^2200 // and if the variances are all equal:201 // variance_combination = variance_individual / N202 // which makes sense --- the standard deviation of the combination is reduced by a factor of sqrt(N).203 204 float sumWeights = 0.0; // Sum of the global weights205 float sumVarianceWeights = 0.0; // Sum of the pixel variances multiplied by the global weights206 for (int i = 0; i < variances->n; i++) {207 if (!(masks->data.PS_TYPE_MASK_DATA[i] & maskVal)) {208 sumVarianceWeights += variances->data.F32[i] * PS_SQR(weights->data.F32[i]);209 sumWeights += weights->data.F32[i];210 }211 }212 213 if (sumWeights <= 0) {214 return false;215 }216 217 *variance = sumVarianceWeights / PS_SQR(sumWeights);218 192 return true; 219 193 } … … 263 237 psVector *pixelData = buffer->pixels; // Values for the pixel of interest 264 238 psVector *pixelMasks = buffer->masks; // Masks for the pixel of interest 265 psVector *pixelVariances = buffer->variances; // Variances for the pixel of interest239 psVector *pixelVariances = variance ? buffer->variances : NULL; // Variances for the pixel of interest 266 240 psVector *pixelWeights = buffer->weights; // Image weights for the pixel of interest 267 241 psVector *pixelSources = buffer->sources; // Sources for the pixel of interest … … 331 305 if (useVariance || !safe) { 332 306 float mean, var; // Mean and variance from combination 333 if (combinationMean(&mean, pixelData, pixelWeights, NULL, maskVal) && 334 (!variance || combinationVariance(&var, pixelVariances, pixelWeights, NULL, maskVal))) { 307 if (combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) { 335 308 imageValue = mean; 336 309 varianceValue = var; … … 354 327 // applied in the first pass might later be accepted. 355 328 float mean, var; // Mean and variance of the combination 356 if (!combinationMean(&mean, pixelData, pixelWeights, pixelMasks, maskVal) || 357 (variance && !combinationVariance(&var, pixelVariances, pixelWeights, pixelMasks, maskVal))) { 329 if (!combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) { 358 330 break; 359 331 } … … 382 354 float median, stdev; // Median and stdev of the combination, for rejection 383 355 384 if (!combinationMedianStdev(&median, useVariance ? NULL : &stdev, pixelData, pixelMasks, 385 maskVal, sort)) { 356 if (!combinationMedianStdev(&median, useVariance ? NULL : &stdev, pixelData, pixelMasks, sort)) { 386 357 psWarning("Bad median/stdev at %d,%d", x, y); 387 358 break;
Note:
See TracChangeset
for help on using the changeset viewer.
