Index: branches/pap/psModules/src/imcombine/pmStack.c =================================================================== --- branches/pap/psModules/src/imcombine/pmStack.c (revision 25967) +++ branches/pap/psModules/src/imcombine/pmStack.c (revision 25975) @@ -30,10 +30,11 @@ #define PIXEL_LIST_BUFFER 100 // Number of entries to add to pixel list at a time #define PIXEL_MAP_BUFFER 2 // Number of entries to add to pixel map at a time -#define ADD_VARIANCE // Allow additional variance (besides variance factor)? +//#define ADD_VARIANCE // Allow additional variance (besides variance factor)? #define NUM_DIRECT_STDEV 5 // For less than this number of values, measure stdev directly + #define TESTING // Enable test output -#define TEST_X 1500-1 // x coordinate to examine -#define TEST_Y 4000-1 // y coordinate to examine +#define TEST_X 4177-1 // x coordinate to examine +#define TEST_Y 2964-1 // y coordinate to examine @@ -291,8 +292,8 @@ // Mark a pixel for inspection // Value in pixel doesn't seem to agree with the stack, so need to look closer -static inline void combineInspect(const psArray *inputs, // Stack data - int x, int y, // Pixel - int source // Source image index - ) +static inline void combineMarkInspect(const psArray *inputs, // Stack data + int x, int y, // Pixel + int source // Source image index + ) { pmStackData *data = inputs->data[source]; // Stack data of interest @@ -308,8 +309,8 @@ // Cannot possibly inspect this pixel and confirm that it's good. // e.g., Only a single input -static inline void combineReject(const psArray *inputs, // Stack data - int x, int y, // Pixel - int source // Source image index - ) +static inline void combineMarkReject(const psArray *inputs, // Stack data + int x, int y, // Pixel + int source // Source image index + ) { pmStackData *data = inputs->data[source]; // Stack data of interest @@ -322,4 +323,505 @@ } +// 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; +} + + +// Extract vectors for simple combination/rejection operations +static void combineExtract(int *num, // Number of good pixels + bool *suspect, // Any suspect pixels? + combineBuffer *buffer, // Buffer with vectors + 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 maskSuspect, // Value to suspect + bool allowSuspect // Allow suspect values? + ) +{ + // Rudimentary error checking + assert(buffer); + assert(image); + assert(mask); + assert(inputs); + + psVector *pixelData = buffer->pixels; // Values 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 + + // Extract the pixel and mask data + int numGood = 0; // Number of good pixels + for (int i = 0, j = 0; i < inputs->n; i++) { + // Check if this pixel has been rejected. Assumes that the rejection pixel list is sorted --- it + // should be because of how pixelMapGenerate works + if (reject && reject->data.U16[j] == i) { + j++; + continue; + } + + pmStackData *data = inputs->data[i]; // Stack data of interest + if (!data) { + continue; + } + + int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout + psImage *mask = data->readout->mask; // Mask of interest + if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & maskVal) { + continue; + } + if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & maskSuspect) { + if (!allowSuspect) { + combineMarkReject(inputs, x, y, i); + continue; + } + if (suspect) { + *suspect = true; + } + } + + psImage *image = data->readout->image; // Image of interest + psImage *variance = data->readout->variance; // Variance map of interest + pixelData->data.F32[numGood] = image->data.F32[yIn][xIn]; + if (variance) { + pixelVariances->data.F32[numGood] = variance->data.F32[yIn][xIn] * addVariance->data.F32[i]; + } + pixelWeights->data.F32[numGood] = data->weight; + pixelSources->data.U16[numGood] = i; + numGood++; + } + pixelData->n = numGood; + if (variance) { + pixelVariances->n = numGood; + } + pixelWeights->n = numGood; + pixelSources->n = numGood; + pixelLimits->n = numGood; + *num = numGood; + +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + for (int i = 0; i < numGood; i++) { + fprintf(stderr, "Input %d (%" PRIu16 "): %f %f (%f) %f\n", + i, pixelSources->data.U16[i], pixelData->data.F32[i], pixelVariances->data.F32[i], + addVariance->data.F32[i], pixelWeights->data.F32[i]); + } + } +#endif + + return; +} + + +// Combine pixels +static void combinePixels(psImage *image, // Combined image, for output + psImage *mask, // Combined mask, for output + psImage *variance, // Combined variance map, for output + int num, // Number of good pixels + combineBuffer *buffer, // Buffer with vectors + int x, int y, // Coordinates of interest; frame of output image + psImageMaskType bad, // Value for bad pixels + bool safe // Safe combination? + ) +{ + psVector *pixelData = buffer->pixels; // Values 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 + + // 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; + } +#ifdef TESTING + else if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Single input to combine, safety on, pixel is bad.\n"); + } +#endif + break; + } + case 2: { + // Automatically accept the mean of the pixels only if we're not playing safe + if (!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 unsafe --> %f %f\n", mean, var); + } +#endif + } + } +#ifdef TESTING + else { + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Two inputs to combine, safety on, pixel is bad\n"); + } + } +#endif + break; + } + default: { + // Can combine without too much worrying + 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 + break; + } + } + + image->data.F32[y][x] = imageValue; + mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = maskValue; + if (variance) { + variance->data.F32[y][x] = varianceValue; + } + + return; +} + + +// Test pixels to be combined +// Returns false to repeat without suspect pixels +static bool combineTest(int num, // Number of good pixels + bool suspect, // Does the stack contain suspect pixels? + 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 pixels in the combination? + ) +{ + if (numIter <= 0) { + return true; + } + + psVector *pixelData = buffer->pixels; // Values for the pixel of interest + 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]); + } + 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 { + 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 (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) { + 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); + } + } + + return true; +} + + + +#if 0 // Given a stack of images, combine with optional rejection. // Pixels in the stack that are rejected are marked for subsequent inspection @@ -333,4 +835,5 @@ 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 @@ -341,5 +844,6 @@ bool safe, // Combine safely? bool rejection, // Reject values marked for inspection from combination? - combineBuffer *buffer // Buffer for combination; to avoid multiple allocations + combineBuffer *buffer, // Buffer for combination; to avoid multiple allocations + bool allowSuspect // Allow suspect pixels in the combination? ) { @@ -361,53 +865,4 @@ psVector *sort = buffer->sort; // Sort buffer - // Extract the pixel and mask data - int num = 0; // Number of good images - for (int i = 0, j = 0; i < inputs->n; i++) { - // Check if this pixel has been rejected. Assumes that the rejection pixel list is sorted --- it - // should be because of how pixelMapGenerate works - if (reject && reject->data.U16[j] == i) { - j++; - continue; - } - - pmStackData *data = inputs->data[i]; // Stack data of interest - if (!data) { - continue; - } - - int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout - psImage *mask = data->readout->mask; // Mask of interest - if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & maskVal) { - continue; - } - - psImage *image = data->readout->image; // Image of interest - psImage *variance = data->readout->variance; // Variance map of interest - pixelData->data.F32[num] = image->data.F32[yIn][xIn]; - if (variance) { - pixelVariances->data.F32[num] = variance->data.F32[yIn][xIn] * addVariance->data.F32[i]; - } - pixelWeights->data.F32[num] = data->weight; - pixelSources->data.U16[num] = i; - num++; - } - pixelData->n = num; - if (variance) { - pixelVariances->n = num; - } - pixelWeights->n = num; - pixelSources->n = num; - pixelLimits->n = num; - -#ifdef TESTING - if (x == TEST_X && y == TEST_Y) { - for (int i = 0; i < num; i++) { - fprintf(stderr, "Input %d (%" PRIu16 "): %f %f %f\n", - i, pixelSources->data.U16[i], pixelData->data.F32[i], pixelVariances->data.F32[i], - pixelWeights->data.F32[i]); - } - } - int numRejected = 0; // Number of rejected inputs -#endif // The sensible thing to do varies according to how many good pixels there are. @@ -439,5 +894,5 @@ } else { if (!rejection) { - combineReject(inputs, x, y, pixelSources->data.U16[0]); + combineMarkReject(inputs, x, y, pixelSources->data.U16[0]); } #ifdef TESTING @@ -478,12 +933,15 @@ float sigma2 = var1 + var2; // Combined variance if (PS_SQR(diff) > PS_SQR(rej) * sigma2) { - // Not consistent: mark both for inspection + // 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 { - combineInspect(inputs, x, y, pixelSources->data.U16[0]); - combineInspect(inputs, x, y, pixelSources->data.U16[1]); + combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]); + combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]); } #ifdef TESTING @@ -494,4 +952,102 @@ 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; + } } } @@ -515,143 +1071,4 @@ // Prepare for clipping iteration - if (numIter > 0) { - pixelMasks->n = num; - psVectorInit(pixelMasks, 0); - 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); - } - } - } - - // The clipping that follows is solely to identify suspect pixels. - // These suspect pixels will be inspected in more detail by other functions. - 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)) { - psWarning("Bad median/stdev at %d,%d", x, y); - break; - } - limit = rej * stdev; -#ifdef TESTING - if (x == TEST_X && y == TEST_Y) { - fprintf(stderr, "Rejecting without variance; rejection limit: %f\n", limit); - } -#endif - } else { -#ifdef TESTING - if (x == TEST_X && y == TEST_Y) { - fprintf(stderr, "Rejecting 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() \ - pixelMasks->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xff; \ - if (!rejection) { \ - combineInspect(inputs, x, y, pixelSources->data.U16[j]); \ - } \ - numClipped++; \ - totalClipped++; - - 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 for rejection: %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, "Rejecting input %d based on variance: %f > %f\n", - j, diff, sqrtf(pixelLimits->data.F32[j])); - } - numRejected++; -#endif - } - } else if (fabsf(diff) > limit) { - MASK_PIXEL_FOR_INSPECTION(); -#ifdef TESTING - if (x == TEST_X && y == TEST_Y) { - fprintf(stderr, "Rejecting input %d based on distribution: %f > %f\n", - j, diff, limit); - } - numRejected++; -#endif - } - } - } - - if (rejection && totalClipped > 0) { - // Get rid of the masked values - // The alternative to this is to make combinationMeanVariance() accept a mask - int good = 0; // Index of good value - for (int i = 0; i < num; i++) { - if (pixelMasks->data.PS_TYPE_VECTOR_MASK_DATA[i]) { - continue; - } - if (i != good) { - pixelData->data.F32[good] = pixelData->data.F32[i]; - pixelVariances->data.F32[good] = pixelVariances->data.F32[i]; - pixelWeights->data.F32[good] = pixelWeights->data.F32[i]; - pixelData->data.F32[good] = pixelData->data.F32[i]; - } - good++; - } - pixelData->n = good; - pixelVariances->n = good; - pixelWeights->n = good; - if (combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) { - imageValue = mean; - varianceValue = var; - maskValue = 0; - } else { - imageValue = NAN; - varianceValue = NAN; - maskValue = bad; - } - } break; @@ -659,20 +1076,7 @@ } - image->data.F32[y][x] = imageValue; - mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = maskValue; - if (variance) { - variance->data.F32[y][x] = varianceValue; - } - -#ifdef TESTING - mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = num; - if (variance) { - variance->data.F32[y][x] = num > 0 ? 1.0 - (float)numRejected / (float)num : 0.0; - } -#endif - return; } - +#endif // Ensure the input array of pmStackData is valid, and get some details out of it @@ -830,6 +1234,6 @@ /// Stack input images -bool pmStackCombine(pmReadout *combined, psArray *input, psImageMaskType maskVal, psImageMaskType bad, - int kernelSize, int numIter, float rej, float sys, float olympic, +bool pmStackCombine(pmReadout *combined, psArray *input, psImageMaskType maskVal, psImageMaskType maskSuspect, + psImageMaskType bad, int kernelSize, int numIter, float rej, float sys, float olympic, bool useVariance, bool safe, bool rejection) { @@ -946,14 +1350,33 @@ if (x == TEST_X && y == TEST_Y) { fprintf(stderr, "Rejected inputs: "); - for (int i = 0; i < reject->n; i++) { - fprintf(stderr, "%d ", reject->data.U16[i]); + if (!reject) { + fprintf(stderr, "\n"); + } else { + for (int i = 0; i < reject->n; i++) { + fprintf(stderr, "%d ", reject->data.U16[i]); + } + fprintf(stderr, "\n"); } - fprintf(stderr, "\n"); } #endif } - combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, - addVariance, reject, x, y, maskVal, bad, numIter, rej, sys, olympic, - useVariance, safe, rejection, buffer); + + int num; // Number of good pixels + 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); + } + } } }