Index: trunk/psModules/src/imcombine/pmStack.c =================================================================== --- trunk/psModules/src/imcombine/pmStack.c (revision 36710) +++ trunk/psModules/src/imcombine/pmStack.c (revision 41892) @@ -756,4 +756,7 @@ // *goodMask &= mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn]; // save the mask bits still used // check for set bits and increment counter as appropriate + // count the number of times a given mask bit is set in the input pixels. + // NOTE: since we have explicitly skipped the pixels with any bad bits, these are only + // the suspect bits (nGoodBits is a bit of a misnomer: it is more like 'nSuspectBitsForGoodInputs' { psImageMaskType value = 0x0001; @@ -961,5 +964,5 @@ // KMM values; - float Punimodal = 1.0,KMMmean = NAN,KMMsigma = NAN,KMMpi = NAN; + float Punimodal = 1.0, KMMmean = NAN, KMMsigma = NAN, KMMpi = NAN; int KMM_MINIMUM_INPUTS = 6; bool useKMM = false; @@ -1510,4 +1513,130 @@ } +# define MIN_GOOD_PERCENTILE 2 +bool pmStackCombineByPercentile( + pmReadout *combined, + psArray *stackData, + psF64 minRange, + psF64 maxRange, + psImageMaskType badMaskBits, // treat these bits as 'bad' + psImageMaskType suspectMaskBits, // treat these bits as 'suspect' + psImageMaskType blankMaskBits // use this mask value for pixels missing input data (distinguish between Ninput = 0 and Ngood = 0?) +) { + int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine + int xSize, ySize; // Size of the output image + + // we need to copy the readouts to their own array for validation + psArray *stackReadouts = psArrayAlloc(stackData->n); + for (int i = 0; i < stackData->n; i++) { + pmStackData *data = stackData->data[i]; // Stack data for this input + pmReadout *ro = data->readout; // Readout of interest + stackReadouts->data[i] = psMemIncrRefCounter(ro); + } + + if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize, stackReadouts)) { + psError(psErrorCodeLast(), false, "Input stack is not valid."); + psFree(stackReadouts); + return false; + } + psFree(stackReadouts); + + psVector *pixelData = psVectorAlloc(stackData->n, PS_TYPE_F32); + psVector *pixelMask = psVectorAlloc(stackData->n, PS_TYPE_VECTOR_MASK); + + for (int y = minInputRows; y < maxInputRows; y++) { + for (int x = minInputCols; x < maxInputCols; x++) { + + int nGood = 0; + for (int i = 0; i < stackData->n; i++) { + + pmStackData *data = stackData->data[i]; // Stack data for this input + pmReadout *ro = data->readout; // Readout of interest + psImage *image = ro->image; + psImage *mask = ro->mask; + + int xIn = x - data->readout->col0; + int yIn = y - data->readout->row0; // Coordinates on input readout + + if (!isfinite(image->data.F32[yIn][xIn])) continue; + if (fabs(image->data.F32[yIn][xIn]) > 1e5) continue; + // XXX need to test the input mask as well here + if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & badMaskBits); + + // XXX save the count of suspect bits + +# if (0) + // count the number of times a given mask bit is set in the input pixels. + // NOTE: since we have explicitly skipped the pixels with any bad bits, these are only + // the suspect bits (nGoodBits is a bit of a misnomer: it is more like 'nSuspectBitsForGoodInputs' + if (1) { + psImageMaskType value = 0x0001; + for (int nbit = 0; nbit < 16; nbit ++) { + if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & value) { + // XXX nGoodBits[nbit] ++; + } + value <<= 1; + } + } +# endif + + pixelData->data.F32[nGood] = image->data.F32[yIn][xIn]; + nGood ++; + } + pixelData->n = nGood; + psVectorSortInPlace (pixelData); + + if (nGood < MIN_GOOD_PERCENTILE) { + combined->image->data.F32[y][x] = NAN; + combined->variance->data.F32[y][x] = NAN; + combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 1; + continue; + } + +# if (0) +# define SUSPECT_FRACTION 0.65 + // set the mask bits if nGoodBits[i] > f*numGood + if (0) { + psImageMaskType value = 0x0001; + for (int nbit = 0; nbit < 16; nbit ++) { + if (nGoodBits[nbit] > SUSPECT_FRACTION*numGood) { + *goodMask |= value; + } + value <<= 1; + } + } +# endif + + int Ns = MIN(MAX(0, minRange * nGood),nGood); // e.g., 0.1 * 50 = 5 + int Ne = MIN(MAX(0, maxRange * nGood),nGood); // e.g., 0.9 * 50 = 45 + int Npt = Ne - Ns; + + float sum = 0.0; + for (int n = Ns; n < Ne; n++) { + sum += pixelData->data.F32[n]; + } + float mean = sum / (float) Npt; + + float varSum = 0.0; + for (int n = Ns; n < Ne; n++) { + varSum += SQ(pixelData->data.F32[n] - mean); + } + // variance on the mean (stdev / sqrt(N))^2 + float varValue = varSum / (Npt - 1) / Npt; + + // XXX this is probably not the correct variance to save + // the predicted variance should be the 1 / sum (1/var) + + combined->image->data.F32[y][x] = mean; + combined->variance->data.F32[y][x] = varValue; + combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0; + } + } + + psFree(pixelData); + psFree(pixelMask); + + return (true); +} + /// Stack input images bool pmStackCombine( @@ -1644,10 +1773,4 @@ expweight = expmaps->variance; } - -/* // Pre-reject inputs using KMM bimodality test. */ -/* if (1) { */ -/* /\* KMMRejectUnpopular(input,x,y); *\/ */ -/* rejection = true; */ -/* } */ // Set up rejection list @@ -1660,20 +1783,6 @@ for (int y = minInputRows; y < maxInputRows; y++) { for (int x = minInputCols; x < maxInputCols; x++) { + // this is only for TESTING: CHECKPIX(x, y, "Combining pixel %d,%d: %x %x %f %f %f %f %d %d %d\n", x, y, badMaskBits, blankMaskBits, iter, rej, sys, olympic, useVariance, safe, rejection); - -/* // Pre-reject inputs using KMM bimodality test. */ -/* if (1) { */ -/* KMMRejectUnpopular(input,x,y); */ -/* /\* rejection = true; *\/ */ -/* } */ -/* else { */ -/* KMMRejectBright(input,x,y); */ -/* } */ - -#ifdef TESTING - if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) { - fprintf (stderr, "combine\n"); - } -# endif psVector *reject = NULL; // Images to reject for this pixel