Index: trunk/psModules/src/imcombine/pmStack.c =================================================================== --- trunk/psModules/src/imcombine/pmStack.c (revision 13457) +++ trunk/psModules/src/imcombine/pmStack.c (revision 13487) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-05-22 03:59:32 $ + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-05-23 22:25:41 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -48,5 +48,6 @@ } -static combineBuffer *combineBufferAlloc(long numImages // Number of images that will be combined +static combineBuffer *combineBufferAlloc(long numImages, // Number of images that will be combined + psStatsOptions stat // Statistic to use ) { @@ -56,5 +57,5 @@ buffer->pixels = psVectorAlloc(numImages, PS_TYPE_F32); buffer->masks = psVectorAlloc(numImages, PS_TYPE_MASK); - buffer->stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); + buffer->stats = psStatsAlloc(stat); return buffer; @@ -78,4 +79,5 @@ const psArray *inputs, // Stack data const psVector *weights, // Weights for data, or NULL + const psVector *reject, // Indices of pixels to reject, or NULL int x, int y, // Coordinates of interest psMaskType maskVal, // Value to mask @@ -96,13 +98,8 @@ int num = inputs->n; // Number of images to combine - if (buffer) { - psMemIncrRefCounter(buffer); - } else { - buffer = combineBufferAlloc(num); - } - psVector *pixelData = buffer->pixels; // Values for the pixel of interest psVector *pixelMasks = buffer->masks; // Masks for the pixel of interest psStats *stats = buffer->stats; // Statistics + // Extract the pixel and mask data @@ -114,8 +111,14 @@ pixelMasks->data.PS_TYPE_MASK_DATA[i] = mask->data.PS_TYPE_MASK_DATA[y][x]; } + if (reject) { + for (int i = 0; i < reject->n; i++) { + pixelMasks->data.PS_TYPE_MASK_DATA[reject->data.U16[i]] |= maskVal; + } + } // Record the value derived with no clipping, because pixels rejected using the harsh clipping applied in // the first pass might later be accepted. - if (!psVectorStats(stats, pixelData, pixelMasks, weights, maskVal)) { + stats->options |= PS_STAT_SAMPLE_MEAN; + if (!psVectorStats(stats, pixelData, weights, pixelMasks, maskVal)) { // Can't do anything about an error except give it a NAN and mask psErrorClear(); @@ -127,16 +130,25 @@ mask->data.PS_TYPE_MASK_DATA[y][x] = 0; + stats->options &= ~ PS_STAT_SAMPLE_MEAN; + long numClipped = LONG_MAX; // Number of pixels clipped psMaskType ignore = maskVal | bad; // Ignore values with this mask value for (int i = 0; i < numIter && numClipped > 0; i++) { +#if 1 + float mean = stats->sampleMedian; + float stdev = 0.74 * (stats->sampleUQ - stats->sampleLQ); // Rough estimate of the standard deviation +#else + float mean = stats->robustMedian; + float stdev = stats->robustStdev; +#endif + float limit = rej * stdev; // Rejection limit numClipped = 0; - float limit = rej * stats->sampleStdev; // Rejection limit for (int j = 0; j < num; j++) { if (pixelMasks->data.PS_TYPE_MASK_DATA[j] & ignore) { continue; } - float diff = fabsf(pixelData->data.F32[j] - stats->sampleMean); // Absolute difference from mean - if (diff > limit) { - pmStackData *data = inputs->data[i]; // Stack data of interest + float diff = pixelData->data.F32[j] - mean; + if (fabsf(diff) > limit) { + pmStackData *data = inputs->data[j]; // Stack data of interest data->pixels = psPixelsAdd(data->pixels, PIXEL_LIST_BUFFER, x, y); pixelMasks->data.PS_TYPE_MASK_DATA[j] |= bad; @@ -144,5 +156,5 @@ } } - if (!psVectorStats(stats, pixelData, pixelMasks, weights, maskVal)) { + if (!psVectorStats(stats, pixelData, weights, pixelMasks, maskVal)) { // Can't do anything about an error except give it a NAN and mask psErrorClear(); @@ -286,5 +298,5 @@ } - if (!psVectorStats(stats, values, valuesMask, NULL, MASK_BAD | MASK_SUSPECT)) { + if (!psVectorStats(stats, values, NULL, valuesMask, MASK_BAD | MASK_SUSPECT)) { psFree(stats); psFree(values); @@ -397,10 +409,17 @@ PS_ASSERT_INT_NONNEGATIVE(numIter, false); PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false); - - + if (havePixels) { + // 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_MASK, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(combined->image, combined->mask, false); + } if (!haveDetector && !haveSky) { psError(PS_ERR_UNEXPECTED_NULL, true, "Nothing to combine!"); return false; } + if (!haveSky) { @@ -415,5 +434,5 @@ } - // Pull out the weights + // Pull out the weightings psVector *weights = psVectorAlloc(num, PS_TYPE_F32); for (int i = 0; i < num; i++) { @@ -422,9 +441,29 @@ } - combineBuffer *buffer = combineBufferAlloc(num); // Buffer for combination + // Buffer for combination + combineBuffer *buffer = combineBufferAlloc(num, numIter == 0 ? PS_STAT_SAMPLE_MEAN : + PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN); if (havePixels) { // Only combine select pixels - + psImage *combinedImage = combined->image; // Combined image + psImage *combinedMask = combined->mask; // Combined mask + + psArray *pixelMap = pixelMapGenerate(input, numCols, numRows); // Map of pixels to source + psPixels *pixels = NULL; // Total list of pixels, with no duplicates + for (int i = 0; i < num; i++) { + pmStackData *data = input->data[i]; // Stacking data; contains the list of pixels + pixels = psPixelsConcatenate(pixels, data->pixels); + data->pixels = psPixelsRealloc(data->pixels, PIXEL_LIST_BUFFER); // Just in case more rejection + } + for (int i = 0; i < pixels->n; i++) { + int x = pixels->data[i].x, y = pixels->data[i].y; // Coordinates of interest + psVector *reject = pixelMapQuery(pixelMap, x, y); // Inspect these images closely + combinePixels(combinedImage, combinedMask, input, weights, reject, x, y, + maskVal, bad, numIter, rej, buffer); + } + psFree(pixels); + psFree(pixelMap); + psTrace("psModules.imcombine", 5, "Additional %ld pixels fixed.\n", pixels->n); } else { // Pull the products out, allocate if necessary @@ -440,11 +479,25 @@ } + // Generate the pixel lists in which to place the rejected pixels + if (numIter != 0) { + for (int i = 0; i < num; i++) { + pmStackData *data = input->data[i]; // Stack data for this input + data->pixels = psPixelsAllocEmpty(PIXEL_LIST_BUFFER); + } + } + for (int y = 0; y < numRows; y++) { for (int x = 0; x < numCols; x++) { - combinePixels(combinedImage, combinedMask, input, weights, x, y, + combinePixels(combinedImage, combinedMask, input, weights, NULL, x, y, maskVal, bad, numIter, rej, buffer); } } + if (psTraceGetLevel("psModules.imcombine") >= 5) { + for (int i = 0; i < num; i++) { + pmStackData *data = input->data[i]; // Stack data for this input + psTrace("psModules.imcombine", 5, "Image %d: %ld pixels to inspect.\n", i, data->pixels->n); + } + } } @@ -492,5 +545,5 @@ } (void)psBinaryOp(seeing, seeing, "*", seeing); - (void)psBinaryOp(seeing, seeing, "-", psScalarAlloc(PS_SQR(seeingMax), PS_TYPE_F32)); + (void)psBinaryOp(seeing, psScalarAlloc(PS_SQR(seeingMax), PS_TYPE_F32), "-", seeing); (void)psUnaryOp(seeing, seeing, "sqrt"); @@ -504,4 +557,5 @@ // Can daisy-chain multiple tests here + reject->n = 0; if (!convolveTest(reject, inspect, input, x, y, seeing, maskVal, extent, threshold)) { @@ -519,5 +573,11 @@ psPixelsAdd(data->pixels, PIXEL_LIST_BUFFER, x, y); } - + } + + if (psTraceGetLevel("psModules.imcombine") >= 5) { + for (int i = 0; i < num; i++) { + pmStackData *data = input->data[i]; // Stack data for this input + psTrace("psModules.imcombine", 5, "Image %d: %ld pixels rejected.\n", i, data->pixels->n); + } } @@ -525,388 +585,2 @@ } - - - - - - - - -#if 0 -/****************************************************************************** -XXX: Directly from Paul Price - *****************************************************************************/ -static psF32 CalcGradient( - psImage *image, - psImage *imageMask, - psS32 x, - psS32 y -) -{ - psTrace("psModules.imcombine", 4, "Calling CalcGradient(%d, %d)\n", x, y); - int num = 0; - psVector *pixels = psVectorAlloc(8, PS_TYPE_F32); // Array of pixels - psVector *mask = psVectorAlloc(8, PS_TYPE_U8); // Corresponding mask - - // Get limits - int xMin = PS_MAX(x - 1, 0); - int xMax = PS_MIN(x + 1, image->numCols - 1); - int yMin = PS_MAX(y - 1, 0); - int yMax = PS_MIN(y + 1, image->numRows - 1); - if (imageMask != NULL) { - for (int j = yMin; j <= yMax; j++) { - for (int i = xMin; i <= xMax; i++) { - if ((i != x) && (j != y) && (0 == imageMask->data.U8[j][i])) { - pixels->data.F32[num] = image->data.F32[j][i]; - mask->data.U8[num] = 0; - num++; - } else { - mask->data.U8[num] = 1; - } - } - } - } else { - // - // This code is simply the previous loop without the imageMask. - // XXX: Consider restructuring this. - // - for (int j = yMin; j <= yMax; j++) { - for (int i = xMin; i <= xMax; i++) { - if ((i != x) && (j != y)) { - pixels->data.F32[num] = image->data.F32[j][i]; - mask->data.U8[num] = 0; - num++; - } else { - mask->data.U8[num] = 1; - } - } - } - } - - pixels->n = num; - mask->n = num; - - // Get the median - psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); - psVectorStats(stats, pixels, NULL, mask, 1); - float median = stats->sampleMedian; - psFree(stats); - psFree(pixels); - psFree(mask); - - psTrace("psModules.imcombine", 4, "Exiting CalcGradient(%d, %d)\n", x, y); - return(median / image->data.F32[y][x]); -} - -/****************************************************************************** -DetermineRegion(image, myOutToIn): for a psImage and a psPlaneTransform to that -image, this routine determines the size of the input image which maps to that -image, and returns the result in a psRegion struct. - -XXX: Basically, this routine is only guaranteed to work if the transform is -linear. - -XXX: Shouldn't this functionality be part of psImageTransform()? - *****************************************************************************/ -static psRegion DetermineRegion(psImage *image, - psPlaneTransform *myOutToIn) -{ - psTrace("psModules.imcombine", 4, "Calling DetermineRegion()\n"); - psRegion myRegion; - myRegion.x0 = PS_MAX_F32; - myRegion.x1 = PS_MIN_F32; - myRegion.y0 = PS_MAX_F32; - myRegion.y1 = PS_MIN_F32; - psPlane in; - psPlane out; - - in.x = 0.0; - in.y = 0.0; - - psPlaneTransformApply(&out, myOutToIn, &in); - if (out.x < myRegion.x0) { - myRegion.x0 = out.x; - } - if (out.x > myRegion.x1) { - myRegion.x1 = out.x; - } - if (out.y < myRegion.y0) { - myRegion.y0 = out.y; - } - if (out.y > myRegion.y1) { - myRegion.y1 = out.y; - } - - in.x = (psF32) (image->numCols); - in.y = 0.0; - psPlaneTransformApply(&out, myOutToIn, &in); - if (out.x < myRegion.x0) { - myRegion.x0 = out.x; - } - if (out.x > myRegion.x1) { - myRegion.x1 = out.x; - } - if (out.y < myRegion.y0) { - myRegion.y0 = out.y; - } - if (out.y > myRegion.y1) { - myRegion.y1 = out.y; - } - - in.x = (psF32) (image->numCols); - ; - in.y = 0.0; - psPlaneTransformApply(&out, myOutToIn, &in); - if (out.x < myRegion.x0) { - myRegion.x0 = out.x; - } - if (out.x > myRegion.x1) { - myRegion.x1 = out.x; - } - if (out.y < myRegion.y0) { - myRegion.y0 = out.y; - } - if (out.y > myRegion.y1) { - myRegion.y1 = out.y; - } - - in.x = (psF32) (image->numCols); - in.y = (psF32) (image->numRows); - psPlaneTransformApply(&out, myOutToIn, &in); - if (out.x < myRegion.x0) { - myRegion.x0 = out.x; - } - if (out.x > myRegion.x1) { - myRegion.x1 = out.x; - } - if (out.y < myRegion.y0) { - myRegion.y0 = out.y; - } - if (out.y > myRegion.y1) { - myRegion.y1 = out.y; - } - - psTrace("psModules.imcombine", 4, "Exiting DetermineRegion()\n"); - return(myRegion); -} - -/****************************************************************************** -XXX: Don't we have a psLib function for this? - *****************************************************************************/ -static psImage *ImageConvertF32(psImage *image) -{ - psTrace("psModules.imcombine", 4, "Calling ImageConvertF32()\n"); - psImage *imgF32 = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32); - - for (psS32 i = 0 ; i < image->numRows ; i++) { - for (psS32 j = 0 ; j < image->numCols ; j++) { - imgF32->data.F32[i][j] = (psF32) image->data.U8[i][j]; - } - } - - psTrace("psModules.imcombine", 4, "Exiting ImageConvertF32()\n"); - return(imgF32); -} - - -// -// The following macros define how big the initial pixel list will be, and -// how much it should be incremented when realloc'ed. -// -#define PS_REJECT_PIXEL_INITIAL_PIXEL_LIST_LENGTH 100 -#define PS_REJECT_PIXEL_INITIAL_PIXEL_LIST_LENGTH_INC 100 -/****************************************************************************** -pmRejectPixels(images, errors, inToOut, outToIn, rejThreshold, -gradLimit) - -XXX: Optimization: we don't need to transform the entire mask image. -XXX: The inToOut and outToIn transforms are confusing. Verify that what - I think they mean syncs with PWP. - *****************************************************************************/ -psArray *pmRejectPixels( - const psArray *images, ///< Array of input images - const psArray *masks, ///< Array of input image masks - const psArray *errors, ///< The pixels which were rejected in the combination - const psArray *inToOut, ///< Transformation from input to output system - const psArray *outToIn, ///< Transformation from output to input system - psF32 rejThreshold, ///< Rejection threshold - psF32 gradLimit ///< Gradient limit -) -{ - psTrace("psModules.imcombine", 3, "Calling pmRejectPixels()\n"); - PS_ASSERT_PTR_NON_NULL(images, NULL); - for (psS32 im = 0 ; im < images->n ; im++) { - psImage *tmpImage = (psImage *) images->data[im]; - PS_ASSERT_IMAGE_NON_NULL(tmpImage, NULL); - PS_ASSERT_IMAGE_NON_EMPTY(tmpImage, NULL); - PS_ASSERT_IMAGE_TYPE(tmpImage, PS_TYPE_F32, NULL); - if (masks != NULL) { - PS_ASSERT_INT_EQUAL(images->n, masks->n, NULL); - psImage *tmpMask = (psImage *) masks->data[im]; - PS_ASSERT_IMAGE_NON_NULL(tmpMask, NULL); - PS_ASSERT_IMAGE_NON_EMPTY(tmpMask, NULL); - PS_ASSERT_IMAGE_TYPE(tmpMask, PS_TYPE_F32, NULL); - PS_ASSERT_IMAGES_SIZE_EQUAL(tmpImage, tmpMask, NULL); - } - PS_ASSERT_IMAGES_SIZE_EQUAL(((psImage *) images->data[0]), tmpImage, NULL); - } - PS_ASSERT_PTR_NON_NULL(errors, NULL); - PS_ASSERT_PTR_NON_NULL(inToOut, NULL); - PS_ASSERT_PTR_NON_NULL(outToIn, NULL); - // Ensure that the psArray parameters have an element for each image. - psS32 numImages = images->n; - PS_ASSERT_INT_EQUAL(numImages, errors->n, NULL); - PS_ASSERT_INT_EQUAL(numImages, inToOut->n, NULL); - PS_ASSERT_INT_EQUAL(numImages, outToIn->n, NULL); - - // - // Create the psArray of psPixelLists, one for each image, for rejected pixels. - // - psArray *rejects = psArrayAlloc(numImages); - for (psS32 im = 0 ; im < numImages ; im++) { - rejects->data[im] = (psPtr *) psPixelsAlloc(PS_REJECT_PIXEL_INITIAL_PIXEL_LIST_LENGTH); - ((psPixels *)(rejects->data[im]))->n = ((psPixels *)(rejects->data[im]))->nalloc; - psPixels *pixels = (psPixels *) rejects->data[im]; - pixels->n = 0; - } - // - // rPtr is used to maintain a count of the questionable pixels for each image. - // - psVector *rPtr = psVectorAlloc(numImages, PS_TYPE_S32); - psVectorInit(rPtr, 0); - - psS32 numCols = ((psImage *) images->data[0])->numCols; - psS32 numRows = ((psImage *) images->data[0])->numRows; - psRegion myRegion = psRegionSet(0, numCols-1, 0, numRows-1); - psU32 maskVal = 1; // XXX: Is this appropriate? - - psPlane *inCoords = psAlloc(sizeof(psPlane)); - psPlane *outCoords = psAlloc(sizeof(psPlane)); - - for (psS32 im = 0 ; im < numImages ; im++) { - // - // Extract data from psArrays. - // - psPixels *pixelList = (psPixels *) errors->data[im]; - - psImage *currImage = (psImage *) images->data[im]; - myRegion.x0 = 0; - myRegion.x1 = currImage->numCols; - myRegion.y0 = 0; - myRegion.y1 = currImage->numRows; - psPlaneTransform *myInToOut = (psPlaneTransform *) inToOut->data[im]; - psPlaneTransform *myOutToIn = (psPlaneTransform *) outToIn->data[im]; - - // - // Create a psU8 mask image from the list of cosmic pixels. - // - psImage *maskImage = NULL; - maskImage = psPixelsToMask(maskImage, pixelList, myRegion, maskVal); - psImage *maskImageF32 = ImageConvertF32(maskImage); - - // - // Transform that mask image into detector coordinate space - // - psRegion myRegionXForm = DetermineRegion(maskImageF32, myOutToIn); - psImage *transformedImage = psImageTransform(NULL, NULL, maskImageF32, NULL, - 0, myOutToIn, myRegionXForm, NULL, - PS_INTERPOLATE_BILINEAR, 0); - - // - // Loop over all cosmic pixels. Transform their coords to detector space. - // If the value of the transformed mask is larger than rejThreshold, then - // this might be a cosmic ray pixel. We then calculate the mean gradient - // in other images. - // - - psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(PS_INTERPOLATE_BILINEAR, - transformedImage, NULL, NULL, - 0, 0.0, 0.0, 0, 0, 0.0); - - for (psS32 p = 0 ; p < pixelList->n ; p++) { - inCoords->x = 0.5 + (psF32) (pixelList->data[p]).x; - inCoords->y = 0.5 + (psF32) (pixelList->data[p]).y; - psPlaneTransformApply(outCoords, myInToOut, inCoords); - double maskVal; - if (!psImageInterpolate(&maskVal, NULL, NULL, outCoords->x, outCoords->y, interp)) { - psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image."); - psFree(interp); - psFree(maskImage); - psFree(maskImageF32); - psFree(transformedImage); - psFree(inCoords); - psFree(outCoords); - psFree(rejects); - return NULL; - } - if (maskVal > rejThreshold) { - - // This is a possible cosmic array pixel. We must calculate the gradient - // at this location in all input images. - psF32 meanGrads = 0.0; - psS32 numGrads = 0; - // - // Loop through all other images, calculate their mean gradient. - // - for (psS32 otherImg = 0 ; otherImg < numImages ; otherImg++) { - if (im != otherImg) { - // Map the outCoords to inCoords that for otherImg space. - psImage *tmpMask = NULL; - if (masks != NULL) { - tmpMask = masks->data[otherImg]; - } - psPlaneTransformApply(inCoords, - (psPlaneTransform * )outToIn->data[otherImg], - outCoords); - psS32 xPix = (int)(inCoords->x + 0.5); - psS32 yPix = (int)(inCoords->y + 0.5); - if ((xPix >= 0) && (xPix <= ((psImage*)(images->data[otherImg]))->numCols - 1) && - (yPix >= 0) && (yPix <= ((psImage*)(images->data[otherImg]))->numRows - 1)) { - meanGrads += CalcGradient(images->data[otherImg], tmpMask, xPix, yPix); - numGrads++; - } - } - } - if (numGrads > 0) { - meanGrads /= (psF32) numGrads; - } else { - // XXX: my idea. Verify with PWP: - meanGrads = 1.0 + gradLimit; - } - - // XXX: The SDRS and the prototype code differ significantly here: - // if (CalcGradient(inputs->data[i], pixelList->data.x, pixelList->data.y) < (gradLimit * meanGrads)) { - if (meanGrads < gradLimit) { - // - // Add this to the list of questionable pixels. We must ensure that the - // pixelList is large enough; if not, we realloc() - // - psS32 ptr = rPtr->data.S32[im]; - psPixels *pixelListPtr = (psPixels *) rejects->data[im]; - if (ptr >= pixelListPtr->nalloc) { - rejects->data[im] = (psPtr *) psPixelsRealloc(((psPixels *) rejects->data[im]), - ((((psPixels *) rejects->data[im])->nalloc) + PS_REJECT_PIXEL_INITIAL_PIXEL_LIST_LENGTH_INC)); - // XXX: Can the realloc() fail? Must we check for NULL? - } - - ((psPixels *) rejects->data[im])->data[ptr].x = (pixelList->data[p]).x; - ((psPixels *) rejects->data[im])->data[ptr].y = (pixelList->data[p]).y; - (rPtr->data.S32[im])++; - // XXX: this pixel ->n increment is wierd - (((psPixels *) rejects->data[im])->n)++; - } - } - } - - psFree(interp); - psFree(maskImage); - psFree(maskImageF32); - psFree(transformedImage); - } - - psFree(inCoords); - psFree(outCoords); - psTrace("psModules.imcombine", 3, "Exiting pmRejectPixels()\n"); - return(rejects); -} -#endif