Index: trunk/psModules/src/imcombine/pmStack.c =================================================================== --- trunk/psModules/src/imcombine/pmStack.c (revision 14200) +++ trunk/psModules/src/imcombine/pmStack.c (revision 14272) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-07-11 01:32:41 $ + * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-07-17 23:33:09 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -41,5 +41,6 @@ psVector *pixels; // Pixel values psVector *masks; // Pixel masks - psStats *stats; // Statistics + psVector *weights; // Pixel weights + psVector *sort; // Buffer for sorting (to get a robust estimator of the standard dev) } combineBuffer; @@ -48,5 +49,6 @@ psFree(buffer->pixels); psFree(buffer->masks); - psFree(buffer->stats); + psFree(buffer->weights); + psFree(buffer->sort); return; } @@ -61,5 +63,6 @@ buffer->pixels = psVectorAlloc(numImages, PS_TYPE_F32); buffer->masks = psVectorAlloc(numImages, PS_TYPE_MASK); - buffer->stats = psStatsAlloc(stat); + buffer->weights = psVectorAlloc(numImages, PS_TYPE_F32); + buffer->sort = psVectorAlloc(numImages, PS_TYPE_F32); return buffer; @@ -77,10 +80,115 @@ +// Generate a mean value 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); + assert(values->n == weights->n); + assert(values->n == masks->n); + assert(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 + 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]; + } + } + if (sumWeight <= 0) { + return false; + } + + *mean = sumValueWeight / sumWeight; + return true; +} + +// Generate a standard deviation for the combination +// Not using psVectorStats because it has additional allocations which slow things down +static bool combinationStdev(float *stdev, // Mean value, to return + 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->type.type == PS_TYPE_F32); + assert(masks->type.type == PS_TYPE_MASK); + assert(sortBuffer && sortBuffer->nalloc >= values->n && sortBuffer->type.type == PS_TYPE_F32); + + int num = 0; // Number of valid values + for (int i = 0; i < values->n; i++) { + if (!(masks->data.PS_TYPE_MASK_DATA[i] & maskVal)) { + sortBuffer->data.F32[num++] = values->data.F32[i]; + } + } + sortBuffer->n = num; + if (!psVectorSortInPlace(sortBuffer)) { + return false; + } + *stdev = 0.74 * (sortBuffer->data.F32[(int)(0.75 * num)] - sortBuffer->data.F32[(int)(0.25 * num)]); + + 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; +} + // Given a stack of images, combine with optional rejection. // Pixels in the stack that are rejected are marked for subsequent static bool combinePixels(psImage *image, // Combined image, for output psImage *mask, // Combined mask, for output + psImage *weight, // Combined variance map, for output const psArray *inputs, // Stack data - const psVector *weights, // Weights for data, or NULL + const psVector *weights, // Global (single value) weights for data, or NULL const psVector *reject, // Indices of pixels to reject, or NULL int x, int y, // Coordinates of interest @@ -104,6 +212,6 @@ psVector *pixelData = buffer->pixels; // Values for the pixel of interest psVector *pixelMasks = buffer->masks; // Masks for the pixel of interest - psStats *stats = buffer->stats; // Statistics - + psVector *pixelWeights = buffer->weights; // Weights for the pixel of interest + psVector *sort = buffer->sort; // Sort buffer // Extract the pixel and mask data @@ -112,6 +220,8 @@ pmStackData *data = inputs->data[i]; // Stack data of interest psImage *image = data->sky->image; // Image of interest + psImage *weight = data->sky->weight; // Weight map of interest psImage *mask = data->sky->mask; // Mask of interest pixelData->data.F32[i] = image->data.F32[y][x]; + pixelWeights->data.F32[i] = weight->data.F32[y][x]; pixelMasks->data.PS_TYPE_MASK_DATA[i] = mask->data.PS_TYPE_MASK_DATA[y][x]; if (pixelMasks->data.PS_TYPE_MASK_DATA[i] & maskVal) { @@ -122,4 +232,7 @@ // Catch the case where everything is masked image->data.F32[y][x] = NAN; + if (weight) { + weight->data.F32[y][x] = NAN; + } mask->data.PS_TYPE_MASK_DATA[y][x] = bad; return false; @@ -134,27 +247,36 @@ // Record the value derived with no clipping, because pixels rejected using the harsh clipping applied in // the first pass might later be accepted. - stats->options |= PS_STAT_SAMPLE_MEAN; - if (!psVectorStats(stats, pixelData, weights, pixelMasks, maskVal)) { - // Can't do anything about an error except give it a NAN and mask - psErrorClear(); + float mean, variance; // Mean and variance of the combination + if (!combinationMean(&mean, pixelData, weights, pixelMasks, maskVal) || + (weight && !combinationVariance(&variance, weights, pixelWeights, pixelMasks, maskVal))) { image->data.F32[y][x] = NAN; mask->data.PS_TYPE_MASK_DATA[y][x] = bad; - return false; - } - image->data.F32[y][x] = stats->sampleMean; + if (weight) { + weight->data.F32[y][x] = NAN; + } + return false; + } + + image->data.F32[y][x] = mean; + if (weight) { + weight->data.F32[y][x] = variance; + } mask->data.PS_TYPE_MASK_DATA[y][x] = 0; - stats->options &= ~ PS_STAT_SAMPLE_MEAN; - + // The clipping that follows is solely to identify suspect pixels. + // These suspect pixels will be inspected in more detail by other functions. long numClipped = LONG_MAX; // Number of pixels clipped psMaskType ignore = maskVal | bad; // Ignore values with this mask value for (int i = 0; i < numIter && numClipped > 0; i++) { -#if 1 - float mean = stats->sampleMedian; - float stdev = 0.74 * (stats->sampleUQ - stats->sampleLQ); // Rough estimate of the standard deviation -#else - float mean = stats->robustMedian; - float stdev = stats->robustStdev; -#endif + float stdev; // Standard deviation of the combination, for rejection + // Only get the mean if it's not the first iteration (because we got the mean earlier) + if ((i != 0 && !combinationMean(&mean, pixelData, weights, pixelMasks, maskVal)) || + !combinationStdev(&stdev, pixelData, pixelMasks, maskVal, sort)) { + image->data.F32[y][x] = NAN; + mask->data.PS_TYPE_MASK_DATA[y][x] = bad; + weight->data.F32[y][x] = NAN; + return false; + } + float limit = rej * stdev; // Rejection limit numClipped = 0; @@ -165,17 +287,12 @@ float diff = pixelData->data.F32[j] - mean; if (fabsf(diff) > limit) { + // Add the pixel as one to inspect pmStackData *data = inputs->data[j]; // Stack data of interest data->pixels = psPixelsAdd(data->pixels, PIXEL_LIST_BUFFER, x, y); + // Mask it so it's not considered in other iterations within this function pixelMasks->data.PS_TYPE_MASK_DATA[j] |= bad; numClipped++; } } - if (!psVectorStats(stats, pixelData, weights, pixelMasks, maskVal)) { - // Can't do anything about an error except give it a NAN and mask - psErrorClear(); - image->data.F32[y][x] = NAN; - mask->data.PS_TYPE_MASK_DATA[y][x] = bad; - return false; - } } @@ -186,9 +303,10 @@ // Ensure the input array of pmStackData is valid, and get some details out of it static bool validateInputData(bool *haveDetector, // Do we have the detector images? - bool *haveSky, // Do we have the sky images? - bool *havePixels, // Do we have lists of pixels? - int *num, // Number of inputs - int *numCols, int *numRows, // Size of (sky) images - psArray *input // Input array of pmStackData to validate + bool *haveSky, // Do we have the sky images? + bool *haveSkyWeights, // Do we have weights in the sky images? + bool *havePixels, // Do we have lists of pixels? + int *num, // Number of inputs + int *numCols, int *numRows, // Size of (sky) images + psArray *input // Input array of pmStackData to validate ) { @@ -203,11 +321,24 @@ PS_ASSERT_IMAGE_NON_NULL(data->detector->image, false); PS_ASSERT_IMAGE_NON_NULL(data->detector->mask, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(data->detector->image, data->detector->mask, false); + PS_ASSERT_IMAGE_TYPE(data->detector->image, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_TYPE(data->detector->mask, PS_TYPE_MASK, false); } *haveSky = (data->sky != NULL); + *haveSkyWeights = false; if (*haveSky) { PS_ASSERT_IMAGE_NON_NULL(data->sky->image, false); + PS_ASSERT_IMAGE_TYPE(data->sky->image, PS_TYPE_F32, false); PS_ASSERT_IMAGE_NON_NULL(data->sky->mask, false); + PS_ASSERT_IMAGE_TYPE(data->sky->mask, PS_TYPE_MASK, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(data->sky->image, data->sky->mask, false); *numCols = data->sky->image->numCols; *numRows = data->sky->image->numRows; + if (data->sky->weight) { + *haveSkyWeights = true; + PS_ASSERT_IMAGE_NON_NULL(data->sky->weight, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(data->sky->image, data->sky->weight, false); + PS_ASSERT_IMAGE_TYPE(data->sky->weight, PS_TYPE_F32, false); + } } *havePixels = (data->pixels != NULL); @@ -243,4 +374,9 @@ PS_ASSERT_IMAGE_SIZE(data->sky->image, *numCols, *numRows, false); PS_ASSERT_IMAGES_SIZE_EQUAL(data->sky->image, data->sky->mask, false); + if (*haveSkyWeights) { + PS_ASSERT_IMAGE_NON_NULL(data->sky->weight, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(data->sky->image, data->sky->weight, false); + PS_ASSERT_IMAGE_TYPE(data->sky->weight, PS_TYPE_F32, false); + } } } @@ -412,8 +548,10 @@ bool haveDetector; // Do we have the detector images? bool haveSky; // Do we have the sky images? + bool haveSkyWeights; // Do we have the sky weight maps? bool havePixels; // Do we have lists of pixels? int num; // Number of inputs int numCols, numRows; // Size of (sky) images - if (!validateInputData(&haveDetector, &haveSky, &havePixels, &num, &numCols, &numRows, input)) { + if (!validateInputData(&haveDetector, &haveSky, &haveSkyWeights, &havePixels, &num, + &numCols, &numRows, input)) { return false; } @@ -461,4 +599,5 @@ psImage *combinedImage = combined->image; // Combined image psImage *combinedMask = combined->mask; // Combined mask + psImage *combinedWeight = combined->weight; // Combined mask psArray *pixelMap = pixelMapGenerate(input, numCols, numRows); // Map of pixels to source @@ -472,5 +611,5 @@ int x = pixels->data[i].x, y = pixels->data[i].y; // Coordinates of interest psVector *reject = pixelMapQuery(pixelMap, x, y); // Inspect these images closely - combinePixels(combinedImage, combinedMask, input, weights, reject, x, y, + combinePixels(combinedImage, combinedMask, combinedWeight, input, weights, reject, x, y, maskVal, bad, numIter, rej, buffer); } @@ -491,4 +630,10 @@ } + psImage *combinedWeights = combined->weight; // Combined weight map + if (haveSkyWeights && !combinedWeights) { + combined->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); + combinedWeights = combined->weight; + } + // Generate the pixel lists in which to place the rejected pixels if (numIter != 0) { @@ -501,5 +646,5 @@ for (int y = 0; y < numRows; y++) { for (int x = 0; x < numCols; x++) { - combinePixels(combinedImage, combinedMask, input, weights, NULL, x, y, + combinePixels(combinedImage, combinedMask, combinedWeights, input, weights, NULL, x, y, maskVal, bad, numIter, rej, buffer); } @@ -538,8 +683,10 @@ bool haveDetector; // Do we have the detector images? bool haveSky; // Do we have the sky images? + bool haveSkyWeights; // Do we have the sky weight maps? bool havePixels; // Do we have lists of pixels? int num; // Number of inputs int numCols, numRows; // Size of (sky) images - if (!validateInputData(&haveDetector, &haveSky, &havePixels, &num, &numCols, &numRows, input)) { + if (!validateInputData(&haveDetector, &haveSky, &haveSkyWeights, &havePixels, + &num, &numCols, &numRows, input)) { return false; }