Index: trunk/psModules/src/imcombine/pmStack.c =================================================================== --- trunk/psModules/src/imcombine/pmStack.c (revision 16624) +++ trunk/psModules/src/imcombine/pmStack.c (revision 16626) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.21 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-02-23 02:33:59 $ + * @version $Revision: 1.22 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-02-23 03:00:30 $ * Copyright 2004-2007 Institute for Astronomy, University of Hawaii * @@ -43,5 +43,6 @@ psVector *pixels; // Pixel values psVector *masks; // Pixel masks - psVector *weights; // Pixel weights + psVector *variances; // Pixel variances + psVector *weights; // Pixel weightings psVector *sources; // Pixel sources (which image did they come from?) psVector *sort; // Buffer for sorting (to get a robust estimator of the standard dev) @@ -52,4 +53,5 @@ psFree(buffer->pixels); psFree(buffer->masks); + psFree(buffer->variances); psFree(buffer->weights); psFree(buffer->sources); @@ -67,4 +69,5 @@ buffer->pixels = psVectorAlloc(numImages, PS_TYPE_F32); buffer->masks = psVectorAlloc(numImages, PS_TYPE_MASK); + buffer->variances = psVectorAlloc(numImages, PS_TYPE_F32); buffer->weights = psVectorAlloc(numImages, PS_TYPE_F32); buffer->sources = psVectorAlloc(numImages, PS_TYPE_U16); @@ -84,5 +87,5 @@ // Generate a mean value for the combination -// Not using psVectorStats because it assumes that the weights are errors, and weights by 1/error^2 +// 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 @@ -216,19 +219,4 @@ } -// Common path for a bad combined pixel -static inline bool combineBad(psImage *image, // Combined image - psImage *mask, // Combined mask - psImage *weight, // Combined variance map - int x, int y, // Coordinates of interest - psMaskType bad // Value for bad pixels - ) -{ - 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; -} // Mark a pixel for inspection @@ -250,5 +238,5 @@ static void combinePixels(psImage *image, // Combined image, for output psImage *mask, // Combined mask, for output - psImage *weight, // Combined variance map, for output + psImage *variance, // Combined variance map, for output const psArray *inputs, // Stack data const psVector *weights, // Global (single value) weights for data, or NULL @@ -271,9 +259,10 @@ assert(rej > 0); assert(buffer); - assert((useVariance && weight) || !useVariance); + assert((useVariance && variance) || !useVariance); psVector *pixelData = buffer->pixels; // Values for the pixel of interest psVector *pixelMasks = buffer->masks; // Masks for the pixel of interest - psVector *pixelWeights = buffer->weights; // Weights for the pixel of interest + psVector *pixelVariances = buffer->variances; // 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 psVector *sort = buffer->sort; // Sort buffer @@ -301,17 +290,20 @@ psImage *image = data->readout->image; // Image of interest - psImage *weight = data->readout->weight; // Weight map of interest + psImage *variance = data->readout->weight; // Variance ("weight") map of interest pixelData->data.F32[num] = image->data.F32[yIn][xIn]; - if (weight) { - pixelWeights->data.F32[num] = weight->data.F32[yIn][xIn]; - } + if (variance) { + pixelVariances->data.F32[num] = variance->data.F32[yIn][xIn]; + } + pixelWeights->data.F32[num] = data->weight; pixelSources->data.U16[num] = i; num++; } pixelData->n = num; - if (weight) { - pixelWeights->n = num; - } + if (variance) { + pixelVariances->n = num; + } + pixelWeights->n = num; pixelSources->n = num; + // The sensible thing to do varies according to how many good pixels there are. @@ -327,6 +319,6 @@ if (!safe) { imageValue = pixelData->data.F32[0]; - if (weight) { - varianceValue = pixelWeights->data.F32[0]; + if (variance) { + varianceValue = pixelVariances->data.F32[0]; } maskValue = 0; @@ -338,9 +330,9 @@ // playing safe if (useVariance || !safe) { - float mean, variance; // Mean and variance from combination - if (combinationMean(&mean, pixelData, weights, pixelMasks, maskVal) && - (!weight || combinationVariance(&variance, pixelWeights, weights, pixelMasks, maskVal))) { + float mean, var; // Mean and variance from combination + if (combinationMean(&mean, pixelData, pixelWeights, NULL, maskVal) && + (!variance || combinationVariance(&var, pixelVariances, pixelWeights, NULL, maskVal))) { imageValue = mean; - varianceValue = variance; + varianceValue = var; maskValue = 0; } @@ -349,5 +341,5 @@ // Use variance to check that the two are consistent float diff = pixelData->data.F32[0] - pixelData->data.F32[1]; - float sigma = sqrtf(pixelWeights->data.F32[0] + pixelWeights->data.F32[1]); + float sigma = sqrtf(pixelVariances->data.F32[0] + pixelVariances->data.F32[1]); if (fabs(diff) > rej * sigma) { // Not consistent: mark both for inspection @@ -361,11 +353,11 @@ // Record the value derived with no clipping, because pixels rejected using the harsh clipping // applied in the first pass might later be accepted. - float mean, variance; // Mean and variance of the combination - if (!combinationMean(&mean, pixelData, weights, pixelMasks, maskVal) || - (weight && !combinationVariance(&variance, pixelWeights, weights, pixelMasks, maskVal))) { + float mean, var; // Mean and variance of the combination + if (!combinationMean(&mean, pixelData, pixelWeights, pixelMasks, maskVal) || + (variance && !combinationVariance(&var, pixelVariances, pixelWeights, pixelMasks, maskVal))) { break; } imageValue = mean; - varianceValue = variance; + varianceValue = var; maskValue = 0; @@ -377,5 +369,5 @@ // Convert to rejection limits for (int i = 0; i < num; i++) { - pixelWeights->data.F32[i] = rej * sqrtf(pixelWeights->data.F32[i]); + pixelVariances->data.F32[i] = rej * sqrtf(pixelVariances->data.F32[i]); } } @@ -406,5 +398,5 @@ } float diff = fabsf(pixelData->data.F32[j] - median); // Difference from expected - if (diff > (useVariance ? pixelWeights->data.F32[i] : limit)) { + if (diff > (useVariance ? pixelVariances->data.F32[i] : limit)) { pixelMasks->data.PS_TYPE_MASK_DATA[j] = 0xff; combineInspect(inputs, x, y, pixelSources->data.U16[j]); @@ -423,5 +415,5 @@ // Ensure the input array of pmStackData is valid, and get some details out of it -static bool validateInputData(bool *haveWeights, // Do we have weights in the sky images? +static bool validateInputData(bool *haveVariances, // Do we have variance maps in the sky images? bool *havePixels, // Do we have lists of pixels? int *num, // Number of inputs @@ -439,5 +431,5 @@ PS_ASSERT_PTR_NON_NULL(data, false); assert(psMemGetDeallocator(data) == (psFreeFunc)stackDataFree); // Ensure it's the right type - *haveWeights = false; + *haveVariances = false; PS_ASSERT_IMAGE_NON_NULL(data->readout->image, false); PS_ASSERT_IMAGE_TYPE(data->readout->image, PS_TYPE_F32, false); @@ -448,5 +440,5 @@ *numRows = data->readout->image->numRows; if (data->readout->weight) { - *haveWeights = true; + *haveVariances = true; PS_ASSERT_IMAGE_NON_NULL(data->readout->weight, false); PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout->weight, false); @@ -476,5 +468,5 @@ PS_ASSERT_IMAGE_SIZE(data->readout->image, *numCols, *numRows, false); PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout->mask, false); - if (*haveWeights) { + if (*haveVariances) { PS_ASSERT_IMAGE_NON_NULL(data->readout->weight, false); PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout->weight, false); @@ -651,9 +643,9 @@ { PS_ASSERT_PTR_NON_NULL(combined, false); - bool haveWeights; // Do we have the weight maps? + bool haveVariances; // Do we have the variance maps? bool havePixels; // Do we have lists of pixels? int num; // Number of inputs int numCols, numRows; // Size of (sky) images - if (!validateInputData(&haveWeights, &havePixels, &num, &numCols, &numRows, input)) { + if (!validateInputData(&haveVariances, &havePixels, &num, &numCols, &numRows, input)) { return false; } @@ -669,10 +661,10 @@ PS_ASSERT_IMAGES_SIZE_EQUAL(combined->image, combined->mask, false); } - if (useVariance && !haveWeights) { + if (useVariance && !haveVariances) { psWarning("Unable to use variance in rejection if no variance maps supplied --- option turned off"); useVariance = false; } - // Pull out the weightings + // Pull out the image weightings psVector *weights = psVectorAlloc(num, PS_TYPE_F32); for (int i = 0; i < num; i++) { @@ -712,5 +704,5 @@ psImage *combinedImage = combined->image; // Combined image psImage *combinedMask = combined->mask; // Combined mask - psImage *combinedWeight = combined->weight; // Combined mask + psImage *combinedVariance = combined->weight; // Combined variance map psArray *pixelMap = pixelMapGenerate(input, maxInputCols, maxInputRows); // Map of pixels to source @@ -731,5 +723,5 @@ } psVector *reject = pixelMapQuery(pixelMap, x, y); // Inspect these images closely - combinePixels(combinedImage, combinedMask, combinedWeight, input, weights, reject, x, y, + combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, reject, x, y, maskVal, bad, numIter, rej, useVariance, safe, buffer); } @@ -750,8 +742,8 @@ } - psImage *combinedWeights = combined->weight; // Combined weight map - if (haveWeights && !combinedWeights) { + psImage *combinedVariance = combined->weight; // Combined variance map + if (haveVariances && !combinedVariance) { combined->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); - combinedWeights = combined->weight; + combinedVariance = combined->weight; } @@ -769,5 +761,5 @@ for (int y = minInputRows; y < maxInputRows; y++) { for (int x = minInputCols; x < maxInputCols; x++) { - combinePixels(combinedImage, combinedMask, combinedWeights, input, weights, NULL, x, y, + combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, NULL, x, y, maskVal, bad, numIter, rej, useVariance, safe, buffer); }