Index: branches/pap/psModules/src/imcombine/pmStack.c =================================================================== --- branches/pap/psModules/src/imcombine/pmStack.c (revision 25975) +++ branches/pap/psModules/src/imcombine/pmStack.c (revision 26007) @@ -34,7 +34,7 @@ -#define TESTING // Enable test output -#define TEST_X 4177-1 // x coordinate to examine -#define TEST_Y 2964-1 // y coordinate to examine +//#define TESTING // Enable test output +//#define TEST_X 3122-1 // x coordinate to examine +//#define TEST_Y 1028-1 // y coordinate to examine @@ -43,9 +43,9 @@ typedef struct { psVector *pixels; // Pixel values - psVector *masks; // Pixel masks psVector *variances; // Pixel variances psVector *weights; // Pixel weightings psVector *sources; // Pixel sources (which image did they come from?) psVector *limits; // Rejection limits + psVector *suspects; // Pixel is suspect? psVector *sort; // Buffer for sorting (to get a robust estimator of the standard dev) } combineBuffer; @@ -54,9 +54,9 @@ { psFree(buffer->pixels); - psFree(buffer->masks); psFree(buffer->variances); psFree(buffer->weights); psFree(buffer->sources); psFree(buffer->limits); + psFree(buffer->suspects); psFree(buffer->sort); return; @@ -70,9 +70,9 @@ buffer->pixels = psVectorAlloc(numImages, PS_TYPE_F32); - buffer->masks = psVectorAlloc(numImages, PS_TYPE_VECTOR_MASK); buffer->variances = psVectorAlloc(numImages, PS_TYPE_F32); buffer->weights = psVectorAlloc(numImages, PS_TYPE_F32); buffer->sources = psVectorAlloc(numImages, PS_TYPE_U16); buffer->limits = psVectorAlloc(numImages, PS_TYPE_F32); + buffer->suspects = psVectorAlloc(numImages, PS_TYPE_U8); buffer->sort = psVectorAlloc(numImages, PS_TYPE_F32); return buffer; @@ -144,23 +144,16 @@ float *stdev, // Standard deviation value, to return const psVector *values, // Values to combine - const psVector *masks, // Mask to apply psVector *sortBuffer // Buffer for sorting ) { assert(values); - assert(!masks || values->n == masks->n); assert(values->type.type == PS_TYPE_F32); - assert(!masks || masks->type.type == PS_TYPE_VECTOR_MASK); assert(sortBuffer && sortBuffer->nalloc >= values->n && sortBuffer->type.type == PS_TYPE_F32); - // Need to filter out clipped values - int num = 0; // Number of valid values - for (int i = 0; i < values->n; i++) { - if (!masks || !masks->data.PS_TYPE_VECTOR_MASK_DATA[i]) { - sortBuffer->data.F32[num++] = values->data.F32[i]; - } - } - sortBuffer->n = num; - if (!psVectorSortInPlace(sortBuffer)) { + int num = values->n; // Number of values + sortBuffer = psVectorSortIndex(sortBuffer, values); + if (!sortBuffer) { + *median = NAN; + *stdev = NAN; return false; } @@ -168,14 +161,15 @@ if (num == 3) { // Attempt to measure standard deviation with only three values (and one of those possibly corrupted) - *median = sortBuffer->data.F32[1]; + *median = values->data.F32[sortBuffer->data.S32[1]]; if (stdev) { - float diff1 = sortBuffer->data.F32[0] - *median; - float diff2 = sortBuffer->data.F32[2] - *median; + float diff1 = values->data.F32[sortBuffer->data.S32[0]] - *median; + float diff2 = values->data.F32[sortBuffer->data.S32[2]] - *median; // This factor of sqrt(2) might not be exact, but it's about right *stdev = M_SQRT2 * PS_MIN(fabsf(diff1), fabsf(diff2)); } } else { - *median = num % 2 ? sortBuffer->data.F32[num / 2] : - (sortBuffer->data.F32[num / 2 - 1] + sortBuffer->data.F32[num / 2]) / 2.0; + *median = num % 2 ? values->data.F32[sortBuffer->data.S32[num / 2]] : + (values->data.F32[sortBuffer->data.S32[num / 2 - 1]] + + values->data.F32[sortBuffer->data.S32[num / 2]]) / 2.0; if (stdev) { if (num <= NUM_DIRECT_STDEV) { @@ -183,11 +177,11 @@ double sum = 0.0; for (int i = 0; i < num; i++) { - sum += PS_SQR(sortBuffer->data.F32[i] - *median); + sum += PS_SQR(values->data.F32[sortBuffer->data.S32[i]] - *median); } *stdev = sqrt(sum / (double)(num - 1)); } else { // Standard deviation from the interquartile range - *stdev = 0.74 * (sortBuffer->data.F32[(int)(0.75 * num)] - - sortBuffer->data.F32[(int)(0.25 * num)]); + *stdev = 0.74 * (values->data.F32[sortBuffer->data.S32[(int)(0.75 * num)]] - + values->data.F32[sortBuffer->data.S32[(int)(0.25 * num)]]); } } @@ -196,72 +190,13 @@ return true; } - -#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 numGood = values->n; // Number of good values int numBad = frac * numGood + 0.5; // Number of bad values @@ -273,7 +208,4 @@ 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) { @@ -297,4 +229,9 @@ ) { +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Marking image %d for inspection\n", source); + } +#endif pmStackData *data = inputs->data[source]; // Stack data of interest if (!data) { @@ -314,4 +251,9 @@ ) { +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Marking pixel image %d for rejection\n", source); + } +#endif pmStackData *data = inputs->data[source]; // Stack data of interest if (!data) { @@ -321,145 +263,4 @@ data->reject = psPixelsAdd(data->reject, data->reject->nalloc, x, y); return; -} - -// General test for multiple pixels -// Returns false if we need to re-run without suspect pixels -static bool combineTestGeneral(int num, // Number of good pixels - bool suspect, // Are there suspect pixels in the list? - psArray *inputs, // Original inputs (for flagging) - combineBuffer *buffer, // Buffer with vectors - int x, int y, // Coordinates of interest; frame of output image - int numIter, // Number of rejection iterations - float rej, // Number of standard deviations at which to reject - float sys, // Relative systematic error - float olympic,// Fraction of values to discard (Olympic weighted mean) - bool useVariance, // Use variance for rejection when combining? - bool safe, // Combine safely? - bool allowSuspect // Allow suspect values? - ) -{ - 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 *pixelWeights = buffer->weights; // Image weights for the pixel of interest - psVector *pixelSources = buffer->sources; // Sources for the pixel of interest - psVector *pixelLimits = buffer->limits; // Limits for the pixel of interest - psVector *sort = buffer->sort; // Sort buffer - - - pixelMasks->n = num; - psVectorInit(pixelMasks, 0); - - // Set up rejection limits - if (useVariance) { - // Convert to rejection limits --- saves doing it later. - // 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 - 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 - pixelLimits->data.F32[i] = rej2 * (pixelVariances->data.F32[i] + sysVar); - } - } - - int numClipped = INT_MAX; // Number of pixels clipped per iteration - int totalClipped = 0; // Total number of pixels clipped - for (int i = 0; i < numIter && numClipped > 0 && num - totalClipped > 2; i++) { - numClipped = 0; - 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)) { - if (i == 0 && suspect) { - // Something's not right --- repeat - return false; - } - for (int j = 0; j < num; j++) { - combineMarkReject(inputs, x, y, pixelSources->data.U16[j]); - } - return true; - } - limit = rej * stdev; -#ifdef TESTING - if (x == TEST_X && y == TEST_Y) { - fprintf(stderr, "Flag without variance; limit: %f\n", limit); - } -#endif - } else { -#ifdef TESTING - if (x == TEST_X && y == TEST_Y) { - fprintf(stderr, "Flag with variance...\n"); - } -#endif - median = combinationWeightedOlympic(pixelData, pixelWeights, pixelMasks, olympic, sort); - } - -#ifdef TESTING - if (x == TEST_X && y == TEST_Y) { - fprintf(stderr, "Median: %f\n", median); - } -#endif - - // Mask a pixel for inspection -#define MASK_PIXEL_FOR_INSPECTION() \ - if (i == 0 && suspect) { \ - /* Something's inconsistent, so want repeat, throwing out all suspect pixels */ \ - return false; \ - } \ - pixelMasks->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xFF; \ - combineMarkInspect(inputs, x, y, pixelSources->data.U16[j]); \ - numClipped++; \ - totalClipped++; - // End - - for (int j = 0; j < num; j++) { - if (pixelMasks->data.PS_TYPE_VECTOR_MASK_DATA[j]) { - continue; - } - float diff = pixelData->data.F32[j] - median; // Difference from expected -#ifdef TESTING - if (x == TEST_X && y == TEST_Y) { - fprintf(stderr, "Testing input %d: %f\n", j, diff); - } -#endif - if (useVariance) { - // Comparing squares --- cheaper than lots of sqrts - // pixelVariances includes the rejection limit, from above - if (PS_SQR(diff) > pixelLimits->data.F32[j]) { - MASK_PIXEL_FOR_INSPECTION(); -#ifdef TESTING - if (x == TEST_X && y == TEST_Y) { - fprintf(stderr, "Flagging input %d based on variance: %f > %f\n", - j, diff, sqrtf(pixelLimits->data.F32[j])); - } -#endif - } - } else if (fabsf(diff) > limit) { - MASK_PIXEL_FOR_INSPECTION(); -#ifdef TESTING - if (x == TEST_X && y == TEST_Y) { - fprintf(stderr, "Flagging input %d based on distribution: %f > %f\n", - j, diff, limit); - } -#endif - } - } - } - - return true; } @@ -478,6 +279,5 @@ int x, int y, // Coordinates of interest; frame of output image psImageMaskType maskVal, // Value to mask - psImageMaskType maskSuspect, // Value to suspect - bool allowSuspect // Allow suspect values? + psImageMaskType maskSuspect // Value to suspect ) { @@ -493,4 +293,9 @@ psVector *pixelSources = buffer->sources; // Sources for the pixel of interest psVector *pixelLimits = buffer->limits; // Limits for the pixel of interest + psVector *pixelSuspects = buffer->suspects; // Is the pixel suspect? + + if (suspect) { + *suspect = false; + } // Extract the pixel and mask data @@ -514,13 +319,7 @@ continue; } - if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & maskSuspect) { - if (!allowSuspect) { - combineMarkReject(inputs, x, y, i); - continue; - } - if (suspect) { - *suspect = true; - } - } + + pixelSuspects->data.U8[numGood] = mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & maskSuspect ? + true : false; psImage *image = data->readout->image; // Image of interest @@ -541,4 +340,5 @@ pixelSources->n = numGood; pixelLimits->n = numGood; + pixelSuspects->n = numGood; *num = numGood; @@ -546,7 +346,7 @@ if (x == TEST_X && y == TEST_Y) { for (int i = 0; i < numGood; i++) { - fprintf(stderr, "Input %d (%" PRIu16 "): %f %f (%f) %f\n", + fprintf(stderr, "Input %d (%" PRIu16 "): %f %f (%f) %f %d\n", i, pixelSources->data.U16[i], pixelData->data.F32[i], pixelVariances->data.F32[i], - addVariance->data.F32[i], pixelWeights->data.F32[i]); + addVariance->data.F32[i], pixelWeights->data.F32[i], pixelSuspects->data.U8[i]); } } @@ -654,4 +454,9 @@ } +#ifdef TESTING + mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = num; +#endif + + return; } @@ -665,155 +470,322 @@ combineBuffer *buffer, // Buffer with vectors int x, int y, // Coordinates of interest; frame of output image - int numIter, // Number of rejection iterations + float iter, // Number of rejection iterations per input float rej, // Number of standard deviations at which to reject float sys, // Relative systematic error float olympic,// Fraction of values to discard (Olympic weighted mean) bool useVariance, // Use variance for rejection when combining? - bool safe, // Combine safely? - bool allowSuspect // Allow suspect pixels in the combination? + bool safe // Combine safely? ) { - if (numIter <= 0) { + if (iter <= 0) { return true; } + int numIter = PS_MAX(iter * num, 1); // Number of iterations + +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Testing pixel %d,%d: %d %f %f %f %d %d\n", + x, y, numIter, rej, sys, olympic, useVariance, safe); + } +#endif + psVector *pixelData = buffer->pixels; // Values for the pixel of interest + psVector *pixelWeights = buffer->weights; // Is the pixel suspect? psVector *pixelVariances = buffer->variances; // Variances for the pixel of interest psVector *pixelSources = buffer->sources; // Sources for the pixel of interest - - switch (num) { - case 0: - break; - case 1: - if (safe) { - combineMarkReject(inputs, x, y, pixelSources->data.U16[0]); + psVector *pixelSuspects = buffer->suspects; // Is the pixel suspect? + psVector *pixelLimits = buffer->limits; // Is the pixel suspect? + + // Set up rejection limits + float rej2 = PS_SQR(rej); // Rejection level squared + if (num > 2 && useVariance) { + // Convert rejection limits --- saves doing it later multiple times + // Using squared rejection limit because it's cheaper than sqrts + 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 + 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 + pixelLimits->data.F32[i] = rej2 * (pixelVariances->data.F32[i] + sysVar); + } + } + + int maskIndex = 0; // Index of pixel to mask + int totalClipped = 0; // Total number of pixels clipped + for (int i = 0; i < numIter && maskIndex >= 0; i++) { + maskIndex = -1; // Nothing to reject + + switch (num) { + case 0: + break; + case 1: + if (i == 0 && safe) { + combineMarkReject(inputs, x, y, pixelSources->data.U16[0]); + } + break; + case 2: { + if (useVariance) { + // Use variance to check that the two are consistent + float diff = 0.5 * (pixelData->data.F32[0] - pixelData->data.F32[1]); // Mean flux difference + float var1 = pixelVariances->data.F32[0]; // Variance of first + float var2 = pixelVariances->data.F32[1]; // Variance of second + // 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; // Combined variance + if (PS_SQR(diff) > rej2 * sigma2) { + // Not consistent: don't believe either! + if (i == 0 && suspect) { + combineMarkReject(inputs, x, y, pixelSources->data.U16[0]); + combineMarkReject(inputs, x, y, pixelSources->data.U16[1]); + } else { + combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]); + combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]); + } +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Flagged both pixels (%f > %f x %f\n)", + diff, rej, sqrtf(sigma2)); + } +#endif + } + } else if (i == 0 && safe) { + // Can't test them, and we want to be safe, so reject + combineMarkReject(inputs, x, y, pixelSources->data.U16[0]); + combineMarkReject(inputs, x, y, pixelSources->data.U16[1]); + } + break; } - break; - case 2: { - if (useVariance) { - // Use variance to check that the two are consistent - float diff = 0.5 * (pixelData->data.F32[0] - pixelData->data.F32[1]); // Mean flux difference - float var1 = pixelVariances->data.F32[0]; // Variance of first - float var2 = pixelVariances->data.F32[1]; // Variance of second - // 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; // Combined variance - if (PS_SQR(diff) > PS_SQR(rej) * sigma2) { - // Not consistent: don't believe either! - if (allowSuspect && suspect) { - combineMarkReject(inputs, x, y, pixelSources->data.U16[0]); - combineMarkReject(inputs, x, y, pixelSources->data.U16[1]); - } else { +#if 0 + case 3: { + // Want to be a bit careful on the rejection than for a larger number of inputs + if (!useVariance) { + return combineTestGeneral(num, suspect, inputs, buffer, x, y, numIter, rej, sys, + olympic, useVariance, safe, allowSuspect); + } + + // Differences between pixel values + float diff01 = pixelData->data.F32[0] - pixelData->data.F32[1]; + float diff12 = pixelData->data.F32[1] - pixelData->data.F32[2]; + float diff20 = pixelData->data.F32[2] - pixelData->data.F32[0]; + // Variance for each pixel + float var0 = pixelVariances->data.F32[0] + PS_SQR(sys * pixelData->data.F32[0]); + float var1 = pixelVariances->data.F32[1] + PS_SQR(sys * pixelData->data.F32[1]); + float var2 = pixelVariances->data.F32[2] + PS_SQR(sys * pixelData->data.F32[2]); + // Errors in pixel differences + float err01 = var0 + var1; + float err12 = var1 + var2; + float err20 = var2 + var0; + +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Diff 0-1: %f %f\n", diff01, err01); + fprintf(stderr, "Diff 1-2: %f %f\n", diff12, err12); + fprintf(stderr, "Diff 2-0: %f %f\n", diff20, err20); + } +#endif + + int badPairs = 0; // Number of bad pairs + bool bad01 = false, bad12 = false, bad20 = false; // Pair is bad? + if (PS_SQR(diff01) > rej2 * err01) { + bad01 = true; + badPairs++; + } + if (PS_SQR(diff12) > rej2 * err12) { + bad12 = true; + badPairs++; + } + if (PS_SQR(diff20) > rej2 * err20) { + bad20 = true; + badPairs++; + } + + if (badPairs > 0 && allowSuspect && suspect) { + return false; + } + + switch (badPairs) { + case 0: + // Nothing to worry about! + break; + case 1: + // Can't tell which image is bad, so be sure to get it + if (bad01) { combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]); combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]); + break; } -#ifdef TESTING - if (x == TEST_X && y == TEST_Y) { - fprintf(stderr, "Flagged both pixels (%f > %f x %f\n)", - diff, rej, sqrtf(sigma2)); + if (bad12) { + combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]); + combineMarkInspect(inputs, x, y, pixelSources->data.U16[2]); + break; } -#endif - } - } else if (safe) { - // Can't test them, and we want to be safe, so reject - combineMarkReject(inputs, x, y, pixelSources->data.U16[0]); - combineMarkReject(inputs, x, y, pixelSources->data.U16[1]); - } - break; - } - case 3: { - // Want to be a bit careful on the rejection than for a larger number of inputs - if (!useVariance) { - return combineTestGeneral(num, suspect, inputs, buffer, x, y, numIter, rej, sys, - olympic, useVariance, safe, allowSuspect); - } - - // Differences between pixel values - float diff01 = pixelData->data.F32[0] - pixelData->data.F32[1]; - float diff12 = pixelData->data.F32[1] - pixelData->data.F32[2]; - float diff20 = pixelData->data.F32[2] - pixelData->data.F32[0]; - // Variance for each pixel - float var0 = pixelVariances->data.F32[0] + PS_SQR(sys * pixelData->data.F32[0]); - float var1 = pixelVariances->data.F32[1] + PS_SQR(sys * pixelData->data.F32[1]); - float var2 = pixelVariances->data.F32[2] + PS_SQR(sys * pixelData->data.F32[2]); - float rej2 = PS_SQR(rej); // Rejection level squared - // Errors in pixel differences - float err01 = var0 + var1; - float err12 = var1 + var2; - float err20 = var2 + var0; - - int badPairs = 0; // Number of bad pairs - bool bad01 = false, bad12 = false, bad20 = false; // Pair is bad? - if (PS_SQR(diff01) > rej2 * err01) { - bad01 = true; - badPairs++; - } - if (PS_SQR(diff12) > rej2 * err12) { - bad12 = true; - badPairs++; - } - if (PS_SQR(diff20) > rej2 * err20) { - bad20 = true; - badPairs++; - } - - if (badPairs > 0 && allowSuspect && suspect) { - return false; - } - - switch (badPairs) { - case 0: - // Nothing to worry about! - break; - case 1: - // Can't tell which image is bad, so be sure to get it - if (bad01) { + if (bad20) { + combineMarkInspect(inputs, x, y, pixelSources->data.U16[2]); + combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]); + break; + } + psAbort("Should never get here"); + case 2: + if (bad01 && bad12) { + // 2 and 0 are good + combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]); + break; + } + if (bad12 && bad20) { + // 0 and 1 are good + combineMarkInspect(inputs, x, y, pixelSources->data.U16[2]); + break; + } + if (bad20 && bad01) { + // 1 and 2 are good + combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]); + break; + } + psAbort("Should never get here"); + case 3: + // Everything's bad combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]); - combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]); - break; - } - if (bad12) { combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]); combineMarkInspect(inputs, x, y, pixelSources->data.U16[2]); break; } - if (bad20) { - combineMarkInspect(inputs, x, y, pixelSources->data.U16[2]); - combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]); - break; - } - psAbort("Should never get here"); - case 2: - if (bad01 && bad12) { - // 2 and 0 are good - combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]); - break; - } - if (bad12 && bad20) { - // 0 and 1 are good - combineMarkInspect(inputs, x, y, pixelSources->data.U16[2]); - break; - } - if (bad20 && bad01) { - // 1 and 2 are good - combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]); - break; - } - psAbort("Should never get here"); - case 3: - // Everything's bad - combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]); - combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]); - combineMarkInspect(inputs, x, y, pixelSources->data.U16[2]); break; } - break; - } - default: { - return combineTestGeneral(num, suspect, inputs, buffer, x, y, numIter, rej, sys, - olympic, useVariance, safe, allowSuspect); - } +#endif + default: { + if (useVariance) { + float median = combinationWeightedOlympic(pixelData, pixelWeights, + olympic, buffer->sort); // Median for stack +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Flag with variance, median = %f\n", median); + } +#endif + float worst = -INFINITY; // Largest deviation + for (int j = 0; j < num; j++) { + float diff = pixelData->data.F32[j] - median; // Difference from expected +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Testing input %d: %f\n", j, diff); + } +#endif + + // Comparing squares --- cheaper than lots of sqrts + // pixelVariances includes the rejection limit, from above + float diff2 = PS_SQR(diff); // Square difference + if (diff2 > pixelLimits->data.F32[j]) { + float dev = diff2 / pixelLimits->data.F32[j]; // Deviation + if (dev > worst) { + worst = dev; + maskIndex = j; + } + } + } + } else { + float median = NAN, stdev = NAN; // Median and stdev of the combination, for rejection + combinationMedianStdev(&median, &stdev, pixelData, buffer->sort); + float limit = rej * stdev; // Rejection limit +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Flag without variance; median = %f, stdev = %f, limit = %f\n", + median, stdev, limit); + } +#endif + float worst = -INFINITY; // Largest deviation + for (int j = 0; j < num; j++) { + float diff = fabsf(pixelData->data.F32[j] - median); // Difference from expected + + if (diff > limit) { + float dev = diff / limit; // Deviation + if (dev > worst) { + worst = dev; + maskIndex = j; + } + } + } + } + } + } + + // Do the actual rejection of the pixel + if (maskIndex >= 0) { + if (suspect) { +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Throwing out all suspect pixels\n"); + } +#endif + // Throw out all suspect pixels + int numGood = 0; // Number of good pixels + for (int j = 0; j < num; j++) { + if (pixelSuspects->data.U8[j]) { + combineMarkReject(inputs, x, y, pixelSources->data.U16[j]); + continue; + } + if (numGood == j) { + numGood++; + continue; + } + pixelData->data.F32[numGood] = pixelData->data.F32[j]; + pixelWeights->data.F32[numGood] = pixelWeights->data.F32[j]; + pixelSources->data.U16[numGood] = pixelSources->data.U16[j]; + pixelLimits->data.F32[numGood] = pixelLimits->data.F32[j]; + pixelVariances->data.F32[numGood] = pixelVariances->data.F32[j]; + numGood++; + } + pixelData->n = numGood; + pixelWeights->n = numGood; + pixelSources->n = numGood; + pixelLimits->n = numGood; + pixelVariances->n = numGood; + totalClipped += num - numGood; + num = numGood; + suspect = false; + } else { + // Throw out masked pixel +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Throwing out input %d\n", maskIndex); + } +#endif + combineMarkInspect(inputs, x, y, pixelSources->data.U16[maskIndex]); + int numGood = 0; // Number of good pixels + for (int j = 0; j < num; j++) { + if (j == maskIndex) { + continue; + } + if (numGood == j) { + numGood++; + continue; + } + pixelData->data.F32[numGood] = pixelData->data.F32[j]; + pixelWeights->data.F32[numGood] = pixelWeights->data.F32[j]; + pixelSources->data.U16[numGood] = pixelSources->data.U16[j]; + pixelLimits->data.F32[numGood] = pixelLimits->data.F32[j]; + pixelVariances->data.F32[numGood] = pixelVariances->data.F32[j]; + numGood++; + } + pixelData->n = numGood; + pixelWeights->n = numGood; + pixelSources->n = numGood; + pixelLimits->n = numGood; + pixelVariances->n = numGood; + totalClipped++; + num--; + } + } } @@ -821,262 +793,4 @@ } - - -#if 0 -// Given a stack of images, combine with optional rejection. -// Pixels in the stack that are rejected are marked for subsequent inspection -static void combinePixels(psImage *image, // Combined image, for output - psImage *mask, // Combined mask, 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 - 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 - psImageMaskType maskVal, // Value to mask - psImageMaskType suspect, // Value to suspect - psImageMaskType bad, // Value to give bad pixels - int numIter, // Number of rejection iterations - float rej, // Number of standard deviations at which to reject - float sys, // Relative systematic error - float olympic,// Fraction of values to discard (Olympic weighted mean) - bool useVariance, // Use variance for rejection when combining? - bool safe, // Combine safely? - bool rejection, // Reject values marked for inspection from combination? - combineBuffer *buffer, // Buffer for combination; to avoid multiple allocations - bool allowSuspect // Allow suspect pixels in the combination? - ) -{ - // Rudimentary error checking - assert(image); - assert(mask); - assert(inputs); - assert(numIter >= 0); - assert(buffer); - assert(addVariance); - 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 *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 - psVector *pixelLimits = buffer->limits; // Limits for the pixel of interest - psVector *sort = buffer->sort; // Sort buffer - - - // The sensible thing to do varies according to how many good pixels there are. - // Default option is that the pixel is bad - float imageValue = NAN, varianceValue = NAN; // Value for combined image and variance map - psImageMaskType maskValue = bad; // Value for combined mask - switch (num) { - 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) { - varianceValue = pixelVariances->data.F32[0]; - } - maskValue = 0; - } else { - if (!rejection) { - combineMarkReject(inputs, x, y, pixelSources->data.U16[0]); - } -#ifdef TESTING - numRejected = 1; - if (x == TEST_X && y == TEST_Y) { - fprintf(stderr, "Single input to combine, safety on, pixel is bad.\n"); - } -#endif - } - break; - } - case 2: { - // Accept the mean of the pixels only if we're going to reject based on the variance, or we're not - // playing safe - if (useVariance || !safe) { - float mean, var; // Mean and variance from combination - if (combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) { - imageValue = mean; - 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 - } - } - if (useVariance && numIter > 0) { - // Use variance to check that the two are consistent - float diff = 0.5 * (pixelData->data.F32[0] - pixelData->data.F32[1]); // Mean flux difference - float var1 = pixelVariances->data.F32[0]; // Variance of first - float var2 = pixelVariances->data.F32[1]; // Variance of second - // 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; // Combined variance - if (PS_SQR(diff) > PS_SQR(rej) * sigma2) { - // Not consistent: don't believe either! - if (rejection) { - imageValue = NAN; - varianceValue = NAN; - maskValue = bad; - } else if (allowSuspect && numSuspect > 0) { - combineMarkReject(inputs, x, y, pixelSources->data.U16[0]); - combineMarkReject(inputs, x, y, pixelSources->data.U16[1]); - } else { - combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]); - combineMarkInspect(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)); - } - numRejected = 2; -#endif - } - } - break; - } - case 3: { - // Can combine without too much worrying, but want to be a bit careful on the rejection - float mean, var; // Mean and variance of the combination - if (!combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) { - break; - } - imageValue = mean; - varianceValue = var; - maskValue = 0; -#ifdef TESTING - if (x == TEST_X && y == TEST_Y) { - fprintf(stderr, "Combined inputs: %f %f\n", mean, var); - } -#endif - - if (numIter > 0) { - if (!useVariance) { - combineRejectionGeneral(); - } else { - // Differences between pixel values - float diff01 = pixelData->data.F32[0] - pixelData->data.F32[1]; - float diff12 = pixelData->data.F32[1] - pixelData->data.F32[2]; - float diff20 = pixelData->data.F32[2] - pixelData->data.F32[0]; - // Variance for each pixel - float var0 = pixelVariances->data.F32[0] + PS_SQR(sys * pixelData->data.F32[0]); - float var1 = pixelVariances->data.F32[1] + PS_SQR(sys * pixelData->data.F32[1]); - float var2 = pixelVariances->data.F32[2] + PS_SQR(sys * pixelData->data.F32[2]); - float rej2 = PS_SQR(rej); // Rejection level squared - // Errors in pixel differences - float err01 = var0 + var1; - float err12 = var1 + var2; - float err20 = var2 + var0; - - int badPairs = 0; // Number of bad pairs - bool bad01 = false, bad12 = false, bad20 = false; // Pair is bad? - if (PS_SQR(diff01) > rej2 * err01) { - bad01 = true; - badPairs++; - } - if (PS_SQR(diff12) > rej2 * err12) { - bad12 = true; - badPairs++; - } - if (PS_SQR(diff20) > rej2 * err20) { - bad20 = true; - badPairs++; - } - - switch (badPairs) { - case 0: - // Nothing to worry about! - break; - case 1: - // Can't tell which image is bad, so be sure to get it - if (bad01) { - combineInspect(inputs, x, y, pixelSources->data.U16[0]); - combineInspect(inputs, x, y, pixelSources->data.U16[1]); - break; - } - if (bad12) { - combineInspect(inputs, x, y, pixelSources->data.U16[1]); - combineInspect(inputs, x, y, pixelSources->data.U16[2]); - break; - } - if (bad20) { - combineInspect(inputs, x, y, pixelSources->data.U16[2]); - combineInspect(inputs, x, y, pixelSources->data.U16[0]); - break; - } - psAbort("Should never get here"); - case 2: - if (bad01 && bad12) { - // 2 and 0 are good - combineInspect(inputs, x, y, pixelSources->data.U16[1]); - break; - } - if (bad12 && bad20) { - // 0 and 1 are good - combineInspect(inputs, x, y, pixelSources->data.U16[2]); - break; - } - if (bad20 && bad01) { - // 1 and 2 are good - combineInspect(inputs, x, y, pixelSources->data.U16[0]); - break; - } - psAbort("Should never get here"); - case 3: - // Everything's bad - combineInspect(inputs, x, y, pixelSources->data.U16[0]); - combineInspect(inputs, x, y, pixelSources->data.U16[1]); - combineInspect(inputs, x, y, pixelSources->data.U16[2]); - break; - } - } - } - break; - } - default: { - // 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, var; // Mean and variance of the combination - if (!combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) { - break; - } - imageValue = mean; - 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 - - break; - } - } - - return; -} -#endif // Ensure the input array of pmStackData is valid, and get some details out of it @@ -1235,5 +949,5 @@ /// Stack input images bool pmStackCombine(pmReadout *combined, psArray *input, psImageMaskType maskVal, psImageMaskType maskSuspect, - psImageMaskType bad, int kernelSize, int numIter, float rej, float sys, float olympic, + psImageMaskType bad, int kernelSize, float iter, float rej, float sys, float olympic, bool useVariance, bool safe, bool rejection) { @@ -1246,8 +960,8 @@ PS_ASSERT_INT_NONNEGATIVE(kernelSize, false); PS_ASSERT_INT_POSITIVE(bad, false); - PS_ASSERT_INT_NONNEGATIVE(numIter, false); if (isnan(rej)) { - PS_ASSERT_INT_EQUAL(numIter, 0, false); + PS_ASSERT_FLOAT_EQUAL(iter, 0, false); } else { + PS_ASSERT_FLOAT_LARGER_THAN(iter, 0, false); PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false); } @@ -1340,6 +1054,6 @@ #ifdef TESTING if (x == TEST_X && y == TEST_Y) { - fprintf(stderr, "Combining pixel %d,%d: %x %x %d %f %f %f %d %d %d\n", - x, y, maskVal, bad, numIter, rej, sys, olympic, useVariance, safe, rejection); + fprintf(stderr, "Combining pixel %d,%d: %x %x %f %f %f %f %d %d %d\n", + x, y, maskVal, bad, iter, rej, sys, olympic, useVariance, safe, rejection); } #endif @@ -1365,17 +1079,10 @@ bool suspect; // Suspect pixels in stack? combineExtract(&num, &suspect, buffer, combinedImage, combinedMask, combinedVariance, - input, weights, addVariance, reject, x, y, maskVal, maskSuspect, true); - if (numIter == 0) { - combinePixels(combinedImage, combinedMask, combinedVariance, - num, buffer, x, y, bad, safe); - } else { - if (combineTest(num, suspect, input, buffer, x, y, numIter, rej, sys, olympic, - useVariance, safe, true)) { - // Need to repeat without suspect pixels - combineExtract(&num, &suspect, buffer, combinedImage, combinedMask, combinedVariance, - input, weights, addVariance, reject, x, y, maskVal, maskSuspect, false); - combineTest(num, suspect, input, buffer, x, y, numIter, rej, sys, olympic, - useVariance, safe, false); - } + input, weights, addVariance, reject, x, y, maskVal, maskSuspect); + combinePixels(combinedImage, combinedMask, combinedVariance, num, buffer, x, y, bad, safe); + + if (iter > 0) { + combineTest(num, suspect, input, buffer, x, y, iter, rej, sys, olympic, + useVariance, safe); } } @@ -1385,4 +1092,6 @@ psFree(weights); psFree(buffer); + + #ifndef PS_NO_TRACE