Index: trunk/psModules/src/imcombine/pmStack.c =================================================================== --- trunk/psModules/src/imcombine/pmStack.c (revision 25380) +++ trunk/psModules/src/imcombine/pmStack.c (revision 26076) @@ -33,7 +33,8 @@ #define NUM_DIRECT_STDEV 5 // For less than this number of values, measure stdev directly + //#define TESTING // Enable test output -//#define TEST_X 1085 // x coordinate to examine -//#define TEST_Y 3371 // y coordinate to examine +//#define TEST_X 3122-1 // x coordinate to examine +//#define TEST_Y 1028-1 // y coordinate to examine @@ -42,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; @@ -53,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; @@ -69,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; @@ -143,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; } @@ -167,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) { @@ -182,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)]]); } } @@ -195,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 @@ -272,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) { @@ -290,9 +223,15 @@ // Mark a pixel for inspection -static inline void combineInspect(const psArray *inputs, // Stack data - int x, int y, // Pixel - int source // Source image index - ) -{ +// Value in pixel doesn't seem to agree with the stack, so need to look closer +static inline void combineMarkInspect(const psArray *inputs, // Stack data + int x, int y, // Pixel + int source // Source image index + ) +{ +#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) { @@ -304,45 +243,62 @@ } -// 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 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 discard,// Fraction of values to discard (Olympic weighted mean) - bool useVariance, // Use variance for rejection when combining? - bool safe, // Combine safely? - bool rejectInspect, // Reject values marked for inspection from combination? - combineBuffer *buffer // Buffer for combination; to avoid multiple allocations - ) +// Mark a pixel for rejection +// Cannot possibly inspect this pixel and confirm that it's good. +// e.g., Only a single input +static inline void combineMarkReject(const psArray *inputs, // Stack data + int x, int y, // Pixel + int source // Source image index + ) +{ +#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) { + psWarning("Can't find input data for source %d", source); + return; + } + data->reject = psPixelsAdd(data->reject, data->reject->nalloc, x, y); + return; +} + + +// 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 + ) { // Rudimentary error checking + assert(buffer); 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 + psVector *pixelSuspects = buffer->suspects; // Is the pixel suspect? + + if (suspect) { + *suspect = false; + } // Extract the pixel and mask data - int num = 0; // Number of good images + 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 @@ -364,45 +320,69 @@ } + pixelSuspects->data.U8[numGood] = mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & maskSuspect ? + true : false; + 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]; + pixelData->data.F32[numGood] = 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; + 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 = num; - } - pixelWeights->n = num; - pixelSources->n = num; - pixelLimits->n = num; + pixelVariances->n = numGood; + } + pixelWeights->n = numGood; + pixelSources->n = numGood; + pixelLimits->n = numGood; + pixelSuspects->n = numGood; + *num = numGood; #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", + for (int i = 0; i < numGood; i++) { + fprintf(stderr, "Input %d (%" PRIu16 "): %f %f (%f) %f %d\n", i, pixelSources->data.U16[i], pixelData->data.F32[i], pixelVariances->data.F32[i], - pixelWeights->data.F32[i]); - } - } -#endif - - // The sensible thing to do varies according to how many good pixels there are. + addVariance->data.F32[i], pixelWeights->data.F32[i], pixelSuspects->data.U8[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 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 @@ -427,7 +407,6 @@ } 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) { + // 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)) { @@ -437,43 +416,20 @@ #ifdef TESTING if (x == TEST_X && y == TEST_Y) { - fprintf(stderr, "Two inputs to combine using variance/unsafe --> %f %f\n", - mean, var); + fprintf(stderr, "Two inputs to combine using 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: mark both for inspection - if (rejectInspect) { - imageValue = NAN; - varianceValue = NAN; - maskValue = bad; - } else { - combineInspect(inputs, x, y, pixelSources->data.U16[0]); - combineInspect(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)); - } -#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: { - // Record the value derived with no clipping, because pixels rejected using the harsh clipping - // applied in the first pass might later be accepted. + // Can combine without too much worrying float mean, var; // Mean and variance of the combination if (!combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) { @@ -488,139 +444,4 @@ } #endif - - // 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, discard, 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 (!rejectInspect) { \ - 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 - 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])); - } -#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); - } -#endif - } - } - } - - if (rejectInspect && 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; } @@ -633,5 +454,341 @@ } +#ifdef TESTING + mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = num; +#endif + + 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 + 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? + ) +{ + 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 + 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; + } +#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; + } + 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; + } +#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--; + } + } + } + + return true; } @@ -639,8 +796,8 @@ // Ensure the input array of pmStackData is valid, and get some details out of it static bool validateInputData(bool *haveVariances, // Do we have variance maps in the sky images? - bool *haveRejects, // Do we have lists of rejected pixels? int *num, // Number of inputs int *numCols, int *numRows, // Size of (sky) images - psArray *input // Input array of pmStackData to validate + const psArray *input, // Input array of pmStackData to validate + const pmReadout *output // Output readout ) { @@ -668,5 +825,5 @@ PS_ASSERT_IMAGE_TYPE(data->readout->variance, PS_TYPE_F32, false); } - *haveRejects = (data->reject != NULL); + bool haveRejects = (data->reject != NULL); // Do we have rejected pixels? // Make sure the rest correspond with the first @@ -681,5 +838,5 @@ return false; } - if ((*haveRejects && !data->reject) || (data->reject && !*haveRejects)) { + if ((haveRejects && !data->reject) || (data->reject && !haveRejects)) { psError(PS_ERR_UNEXPECTED_NULL, true, "The rejected pixels are specified in some but not all inputs."); @@ -697,4 +854,13 @@ PS_ASSERT_IMAGE_TYPE(data->readout->variance, PS_TYPE_F32, false); } + } + + PM_ASSERT_READOUT_NON_NULL(output, false); + if (output->image) { + PS_ASSERT_IMAGE_NON_NULL(output->image, false); + PS_ASSERT_IMAGE_TYPE(output->image, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_NON_NULL(output->mask, false); + PS_ASSERT_IMAGE_TYPE(output->mask, PS_TYPE_IMAGE_MASK, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(output->image, output->mask, false); } @@ -782,31 +948,21 @@ /// Stack input images -bool pmStackCombine(pmReadout *combined, psArray *input, psImageMaskType maskVal, psImageMaskType bad, - int kernelSize, int numIter, float rej, float sys, float discard, - bool entire, bool useVariance, bool safe, bool rejectInspect) -{ - PS_ASSERT_PTR_NON_NULL(combined, false); +bool pmStackCombine(pmReadout *combined, psArray *input, psImageMaskType maskVal, psImageMaskType maskSuspect, + psImageMaskType bad, int kernelSize, float iter, float rej, float sys, float olympic, + bool useVariance, bool safe, bool rejection) +{ bool haveVariances; // Do we have the variance maps? - bool haveRejects; // Do we have lists of rejected pixels? int num; // Number of inputs int numCols, numRows; // Size of (sky) images - if (!validateInputData(&haveVariances, &haveRejects, &num, &numCols, &numRows, input)) { + if (!validateInputData(&haveVariances, &num, &numCols, &numRows, input, combined)) { return false; } 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); - } - if (haveRejects) { - // This is a subsequent combination, so expect that the image and mask already exist - PS_ASSERT_IMAGE_NON_NULL(combined->image, false); - PS_ASSERT_IMAGE_TYPE(combined->image, PS_TYPE_F32, false); - PS_ASSERT_IMAGE_NON_NULL(combined->mask, false); - PS_ASSERT_IMAGE_TYPE(combined->mask, PS_TYPE_IMAGE_MASK, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(combined->image, combined->mask, false); } if (useVariance && !haveVariances) { @@ -833,6 +989,12 @@ } #endif - if (!haveRejects && !data->inspect) { - data->inspect = psPixelsAllocEmpty(PIXEL_LIST_BUFFER); + if (!rejection) { + // Ensure pixels can be put on the appropriate list + if (!data->inspect) { + data->inspect = psPixelsAllocEmpty(PIXEL_LIST_BUFFER); + } + if (!data->reject) { + data->reject = psPixelsAllocEmpty(PIXEL_LIST_BUFFER); + } } } @@ -863,91 +1025,85 @@ combineBuffer *buffer = combineBufferAlloc(num); - if (haveRejects) { - psImage *combinedImage = combined->image; // Combined image - psImage *combinedMask = combined->mask; // Combined mask - psImage *combinedVariance = combined->variance; // Combined variance map - - psArray *pixelMap = pixelMapGenerate(input, minInputCols, maxInputCols, - minInputRows, maxInputRows); // Map of pixels to source - psPixels *pixels = NULL; // Total list of pixels, with no duplicates + // Pull the products out, allocate if necessary + psImage *combinedImage = combined->image; // Combined image + if (!combinedImage) { + combined->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); + combinedImage = combined->image; + } + psImage *combinedMask = combined->mask; // Combined mask + if (!combinedMask) { + combined->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); + combinedMask = combined->mask; + } + + psImage *combinedVariance = combined->variance; // Combined variance map + if (haveVariances && !combinedVariance) { + combined->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); + combinedVariance = combined->variance; + } + + // Set up rejection list + psArray *pixelMap = NULL; // Map of pixels to source + if (rejection) { + pixelMap = pixelMapGenerate(input, minInputCols, maxInputCols, minInputRows, maxInputRows); + } + + // Combine each pixel + for (int y = minInputRows; y < maxInputRows; y++) { + for (int x = minInputCols; x < maxInputCols; x++) { +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + 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 + psVector *reject = NULL; // Images to reject for this pixel + if (rejection) { + reject = pixelMapQuery(pixelMap, minInputCols, minInputRows, x, y); +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Rejected inputs: "); + if (!reject) { + fprintf(stderr, "\n"); + } else { + for (int i = 0; i < reject->n; i++) { + fprintf(stderr, "%d ", reject->data.U16[i]); + } + fprintf(stderr, "\n"); + } + } +#endif + } + + 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); + 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); + } + } + } + + psFree(pixelMap); + psFree(weights); + psFree(buffer); + + + +#ifndef PS_NO_TRACE + if (!rejection && psTraceGetLevel("psModules.imcombine") >= 5) { for (int i = 0; i < num; i++) { - pmStackData *data = input->data[i]; // Stacking data; contains the list of pixels - if (!data) { + pmStackData *data = input->data[i]; // Stack data for this input + if (!data || !data->inspect) { continue; } - pixels = psPixelsConcatenate(pixels, data->reject); - } - pixels = psPixelsDuplicates(pixels, pixels); - - if (entire) { - // Combine entire image - for (int y = minInputRows; y < maxInputRows; y++) { - for (int x = minInputCols; x < maxInputCols; x++) { - psVector *reject = pixelMapQuery(pixelMap, minInputCols, minInputRows, - x, y); // Reject these images - combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, - addVariance, reject, x, y, maskVal, bad, numIter, rej, sys, discard, - useVariance, safe, rejectInspect, buffer); - } - } - } else { - // Only combine previously rejected pixels - for (int i = 0; i < pixels->n; i++) { - // Pixel coordinates are in the frame of the original image - int x = pixels->data[i].x, y = pixels->data[i].y; // Coordinates of interest - if (x < minInputCols || x >= maxInputCols || y < minInputRows || y >= maxInputRows) { - continue; - } - psVector *reject = pixelMapQuery(pixelMap, minInputCols, minInputRows, - x, y); // Reject these images - combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, - addVariance, reject, x, y, maskVal, bad, numIter, rej, sys, discard, - useVariance, safe, rejectInspect, buffer); - } - } - psFree(pixels); - psFree(pixelMap); - } else { - // Pull the products out, allocate if necessary - psImage *combinedImage = combined->image; // Combined image - if (!combinedImage) { - combined->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); - combinedImage = combined->image; - } - psImage *combinedMask = combined->mask; // Combined mask - if (!combinedMask) { - combined->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); - combinedMask = combined->mask; - } - - psImage *combinedVariance = combined->variance; // Combined variance map - if (haveVariances && !combinedVariance) { - combined->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); - combinedVariance = combined->variance; - } - - for (int y = minInputRows; y < maxInputRows; y++) { - for (int x = minInputCols; x < maxInputCols; x++) { - combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, - addVariance, NULL, x, y, maskVal, bad, numIter, rej, sys, discard, - useVariance, safe, rejectInspect, buffer); - } - } - -#ifndef PS_NO_TRACE - if (psTraceGetLevel("psModules.imcombine") >= 5) { - for (int i = 0; i < num; i++) { - pmStackData *data = input->data[i]; // Stack data for this input - if (!data || !data->inspect) { - continue; - } - psTrace("psModules.imcombine", 5, "Image %d: %ld pixels to inspect.\n", i, data->inspect->n); - } - } -#endif - } - - psFree(weights); - psFree(buffer); + psTrace("psModules.imcombine", 5, "Image %d: %ld pixels to inspect.\n", i, data->inspect->n); + } + } +#endif return true;