Index: trunk/psModules/src/pmImageCombine.c =================================================================== --- trunk/psModules/src/pmImageCombine.c (revision 4185) +++ trunk/psModules/src/pmImageCombine.c (revision 4217) @@ -1,11 +1,15 @@ /** @file pmImageCombine.c * - * This file will ... + * This file will perform image combination of several images of the + * same field, produce a list of questionable pixels, then tag some + * of those pixels as cosmic rays. * - * @author Paul Price, IfA + * @author Paul Price, IfA (original prototype) * @author GLG, MHPCC * - * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-06-09 06:16:43 $ + * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-06-13 19:36:12 $ + * + * XXX: pmRejectPixels() has a known bug with the pmImageTransform() call. * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -17,7 +21,4 @@ #include #include "pslib.h" -#include "psConstants.h" -// XXX: Remove this, once psPixels.h is handled properly. -#include "psPixels.h" // @@ -27,7 +28,10 @@ #define PS_COMBINE_IMAGE_INITIAL_PIXEL_LIST_LENGTH 100 #define PS_COMBINE_IMAGE_INITIAL_PIXEL_LIST_LENGTH_INC 100 +#define PS_COMBINE_IMAGE_MAX_QUESTIONABLE_PIXELS 1000 /****************************************************************************** pmCombineImages(combine, questionablePixels, images, errors, masks, maskVal, pixels, numIter, sigmaClip, stats) + +XXX: Allocate a dummy psStats structure so that we don't destroy away its data. *****************************************************************************/ psImage *pmCombineImages(psImage *combine, ///< Combined image (output) @@ -70,24 +74,25 @@ for (psS32 i = 0 ; i < numImages ; i++) { psImage *tmpDataImg = (psImage *) images->data[i]; - psImage *tmpErrorImg = (psImage *) errors->data[i]; - psImage *tmpMaskImg = (psImage *) masks->data[i]; - PS_ASSERT_IMAGE_NON_NULL(tmpDataImg, combine); - PS_ASSERT_IMAGE_NON_NULL(tmpErrorImg, combine); - PS_ASSERT_IMAGE_NON_NULL(tmpMaskImg, combine); - PS_ASSERT_IMAGE_TYPE(tmpDataImg, PS_TYPE_F32, combine); - PS_ASSERT_IMAGE_TYPE(tmpErrorImg, PS_TYPE_F32, combine); - PS_ASSERT_IMAGE_TYPE(tmpMaskImg, PS_TYPE_U8, combine); - if ((tmpDataImg->numRows != numRows) || (tmpDataImg->numCols != numCols)) { - psError(PS_ERR_UNKNOWN, true, "error image %d has size (%d, %d); should be (%d, %d).\n", + psError(PS_ERR_UNKNOWN, true, "image %d has size (%d, %d); should be (%d, %d).\n", i, tmpDataImg->numRows, tmpDataImg->numCols, numRows, numCols); } - PS_ASSERT_IMAGES_SIZE_EQUAL(tmpDataImg, tmpErrorImg, NULL); - PS_ASSERT_IMAGES_SIZE_EQUAL(tmpDataImg, tmpMaskImg, NULL); - } - PS_ASSERT_PTR_NON_NULL(pixels, combine); + if (errors != NULL) { + psImage *tmpErrorImg = (psImage *) errors->data[i]; + PS_ASSERT_IMAGE_NON_NULL(tmpErrorImg, combine); + PS_ASSERT_IMAGE_TYPE(tmpErrorImg, PS_TYPE_F32, combine); + PS_ASSERT_IMAGES_SIZE_EQUAL(tmpDataImg, tmpErrorImg, NULL); + } + + if (masks != NULL) { + psImage *tmpMaskImg = (psImage *) masks->data[i]; + PS_ASSERT_IMAGE_NON_NULL(tmpMaskImg, combine); + PS_ASSERT_IMAGE_TYPE(tmpMaskImg, PS_TYPE_U8, combine); + PS_ASSERT_IMAGES_SIZE_EQUAL(tmpDataImg, tmpMaskImg, NULL); + } + } PS_ASSERT_PTR_NON_NULL(stats, combine); @@ -112,4 +117,5 @@ psFree((*questionablePixels)->data[im]); ((*questionablePixels)->data[im]) = (psPtr *) psPixelsAlloc(PS_COMBINE_IMAGE_INITIAL_PIXEL_LIST_LENGTH); + ((psPixels *) ((*questionablePixels)->data[im]))->n = 0; } // @@ -128,5 +134,5 @@ psVector *pixelMask = NULL; if (masks != NULL) { - pixelMask = psVectorAlloc(numImages, PS_TYPE_F32); + pixelMask = psVectorAlloc(numImages, PS_TYPE_U8); PS_VECTOR_SET_U8(pixelMask, 0); } @@ -140,5 +146,7 @@ if (pixels != NULL) { // Only those specified pixels should be combined. + psStats *stdevStats = psStatsAlloc(PS_STAT_SAMPLE_STDEV); + for (psS32 p = 0 ; p < pixels->n ; p++) { // Must initialize the mask to 0 for every pixel. @@ -154,8 +162,7 @@ // Set the pixel data pixelData->data.F32[im] = ((psImage *) images->data[im])->data.F32[row][col]; - // Set the pixel mask data, if necessary if (masks != NULL) { - pixelMask->data.F32[im] = ((psImage *) masks->data[im])->data.F32[row][col]; + pixelMask->data.U8[im] = ((psImage *) masks->data[im])->data.F32[row][col]; } @@ -171,9 +178,10 @@ for (psS32 iter = 0 ; iter < numIter ; iter++) { // Combine all the pixels, using the specified stat. + stats = psVectorStats((psStats *) stats, pixelData, pixelErrors, pixelMask, maskVal); psF64 combinedPixel; psBool rc = p_psGetStatValue(stats, &combinedPixel); if (rc != true) { - // XXX: Warning. + psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not combine pixels (%d, %d) with the specified stat.\n", row, col); } if (iter == 0) { @@ -185,27 +193,61 @@ // the combined pixel value. // + psS32 numRejects = 0; for (psS32 im = 0 ; im < numImages ; im++) { - stats = psVectorStats(stdevStats, pixelData, pixelErrors, pixelMask, maskVal); + stdevStats = psVectorStats(stdevStats, pixelData, pixelErrors, pixelMask, maskVal); psF64 stdev; psBool rc = p_psGetStatValue(stdevStats, &stdev); if (rc != true) { - // XXX: Warning. + psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not compute the standard deviation of pixel (%d, %d).\n", row, col); + psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not combine pixels (%d, %d) with the specified stat.\n", row, col); } - if (fabs(pixelData->data.F32[im] - combinedPixel) > fabs(sigmaClip * stdev)) { - pixelMask->data.F32[im] = maskVal; - // - // XXX: These data structures indirections are getting complicated. - // - psS32 ptr = qpPtr->data.S32[im]; - psPixels *pixelListPtr = ((psPixels *) ((*questionablePixels)->data[im])); - if (ptr >= pixelListPtr->n) { - (*questionablePixels)->data[im] = (psPtr *) psPixelsRealloc(((psPixels *) ((*questionablePixels)->data[im])), ((((psPixels *) ((*questionablePixels)->data[im]))->n) + PS_COMBINE_IMAGE_INITIAL_PIXEL_LIST_LENGTH_INC)); + if (!(pixelMask->data.U8[im] & maskVal)) { + if (fabs(pixelData->data.F32[im]) > + (fabs(sigmaClip * stdev) + fabs(combinedPixel))) { + numRejects++; + //printf("Rejecting pixel. Image %d. (row, col) (%d, %d)\n", im, row, col); + pixelMask->data.U8[im] = maskVal; + // + // XXX: These data structures indirections are getting complicated. + // + psS32 ptr = qpPtr->data.S32[im]; + psPixels *pixelListPtr = ((psPixels *) ((*questionablePixels)->data[im])); + + if (ptr >= pixelListPtr->nalloc) { + (*questionablePixels)->data[im] = + (psPtr *) psPixelsRealloc(((psPixels *) ((*questionablePixels)->data[im])), + ((((psPixels *) ((*questionablePixels)->data[im]))->nalloc) + + PS_COMBINE_IMAGE_INITIAL_PIXEL_LIST_LENGTH_INC)); + // XXX: Can the realloc() fail? Must we check for NULL? + } + ((psPixels *) ((*questionablePixels)->data[im]))->data[ptr].x = col; + ((psPixels *) ((*questionablePixels)->data[im]))->data[ptr].y = row; + (qpPtr->data.S32[im])++; + ((psPixels *) ((*questionablePixels)->data[im]))->n = qpPtr->data.S32[im]; } - // XXX: Can the realloc() fail? Must we check for NULL? - ((psPixels *) ((*questionablePixels)->data[im]))->data[ptr].x = col; - ((psPixels *) ((*questionablePixels)->data[im]))->data[ptr].y = row; - (qpPtr->data.S32[im])++; } + } + + // + // If the number of rejected pixels is zero, then there's no point + // continuing the loop. + // + if (numRejects == 0) { + break; + } + psS32 totalRejects = 0; + for (psS32 im = 0 ; im < numImages ; im++) { + if (pixelMask->data.U8[im] & maskVal) { + totalRejects++; + } + } + + // + // XXX: Is it possible to have all pixels rejected? If so, we should + // exit the loop. + // + if (totalRejects == numImages) { + break; } } @@ -231,5 +273,5 @@ // Set the pixel mask data, if necessary if (masks != NULL) { - pixelMask->data.F32[im] = ((psImage *) masks->data[im])->data.F32[row][col]; + pixelMask->data.U8[im] = ((psImage *) masks->data[im])->data.F32[row][col]; } @@ -246,5 +288,5 @@ combine->data.F32[row][col] = (psF32) tmpF64; if (rc != true) { - // XXX: Warning. + psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not combine pixels (%d, %d) with the specified stat.\n", row, col); } } @@ -299,4 +341,113 @@ 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) +{ + 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; + } + + return(myRegion); +} + +/****************************************************************************** +XXX: Don't we have a psLib function for this? + *****************************************************************************/ +static psImage *ImageConvertF32(psImage *image) +{ + 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]; + } + } + + return(imgF32); +} + // @@ -314,4 +465,5 @@ I think they mean syncs with PWP. XXX: Must add mask parameter, use it in gradient calculation. +XXX: This function does not work. *****************************************************************************/ psArray *pmRejectPixels(const psArray *images, ///< Array of input images @@ -323,4 +475,5 @@ ) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: pmRejectPixels() has known bugs. Specifically, in the psImageTransform() call.\n"); PS_ASSERT_PTR_NON_NULL(images, NULL); PS_ASSERT_PTR_NON_NULL(errors, NULL); @@ -332,5 +485,5 @@ PS_ASSERT_INT_EQUAL(numImages, inToOut->n, NULL); PS_ASSERT_INT_EQUAL(numImages, outToIn->n, NULL); - // XXX: Loop through all images, ensure their sizes are equal? + // XXX: Loop through all images, ensure their sizes are equal // @@ -340,4 +493,6 @@ for (psS32 im = 0 ; im < numImages ; im++) { rejects->data[im] = (psPtr *) psPixelsAlloc(PS_REJECT_PIXEL_INITIAL_PIXEL_LIST_LENGTH); + psPixels *pixels = (psPixels *) rejects->data[im]; + pixels->n = 0; } // @@ -346,10 +501,9 @@ psVector *rPtr = psVectorAlloc(numImages, PS_TYPE_S32); PS_VECTOR_SET_S32(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; + psU32 maskVal = 1; // XXX: Is this appropriate? psPlane *inCoords = psAlloc(sizeof(psPlane)); @@ -361,9 +515,10 @@ // psPixels *pixelList = (psPixels *) errors->data[im]; + psImage *currImage = (psImage *) images->data[im]; myRegion.x0 = 0; - myRegion.x1 = currImage->numCols - 1; + myRegion.x1 = currImage->numCols; myRegion.y0 = 0; - myRegion.y1 = currImage->numRows - 1; + myRegion.y1 = currImage->numRows; psPlaneTransform *myInToOut = (psPlaneTransform *) inToOut->data[im]; psPlaneTransform *myOutToIn = (psPlaneTransform *) outToIn->data[im]; @@ -374,11 +529,21 @@ psImage *maskImage = NULL; maskImage = psPixelsToMask(maskImage, pixelList, myRegion, maskVal); + psImage *maskImageF32 = ImageConvertF32(maskImage); // // Transform that mask image into detector coordinate space // - // XXX: conform to new psImageTransform() - // psImage *transformedImage = psImageTransform(NULL, maskImage, NULL, 0, myOutToIn, NULL, 0); - psImage *transformedImage = NULL; + psRegion myRegion = DetermineRegion(maskImageF32, myOutToIn); + psImage *transformedImage = psImageTransform(NULL, NULL, maskImageF32, NULL, + 0, myOutToIn, myRegion, NULL, + PS_INTERPOLATE_BILINEAR, 0); + // + // XXX: Currently, a possibly buggy psImageTransform() corrupts the data in + // in outToIn, and possibly other places. The following printf() demonstrates + // this in conjunction with the test code. + // + // psPlaneTransform *tmpOutToIn2 = (psPlaneTransform *) outToIn->data[0]; + // printf("HMMM3: (%d, %d) (%d %d)\n", tmpOutToIn2->x->nX, tmpOutToIn2->x->nY, tmpOutToIn2->y->nX, tmpOutToIn2->y->nY); + // // @@ -392,5 +557,4 @@ inCoords->y = 0.5 + (psF32) (pixelList->data[p]).y; psPlaneTransformApply(outCoords, myInToOut, inCoords); - psF32 maskVal = (psF32) psImagePixelInterpolate(transformedImage, outCoords->x, outCoords->y, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); @@ -408,6 +572,7 @@ if (im != otherImg) { // Map the outCoords to inCoords that for otherImg space. - psPlaneTransformApply(inCoords, (psPlaneTransform * )outToIn->data[otherImg], outCoords); - + psPlaneTransformApply(inCoords, + (psPlaneTransform * )outToIn->data[otherImg], + outCoords); psS32 xPix = (int)(inCoords->x + 0.5); psS32 yPix = (int)(inCoords->y + 0.5); @@ -417,5 +582,4 @@ numGrads++; } - } } @@ -434,11 +598,13 @@ // pixelList is large enough; if not, we realloc() // + psS32 ptr = rPtr->data.S32[im]; psPixels *pixelListPtr = (psPixels *) rejects->data[im]; - if (ptr >= pixelListPtr->n) { + if (ptr >= pixelListPtr->nalloc) { rejects->data[im] = (psPtr *) psPixelsRealloc(((psPixels *) rejects->data[im]), - ((((psPixels *) rejects->data[im])->n) + PS_REJECT_PIXEL_INITIAL_PIXEL_LIST_LENGTH_INC)); + ((((psPixels *) rejects->data[im])->nalloc) + PS_REJECT_PIXEL_INITIAL_PIXEL_LIST_LENGTH_INC)); + // XXX: Can the realloc() fail? Must we check for NULL? } - // 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; @@ -450,11 +616,11 @@ psFree(myOutToIn); psFree(maskImage); + psFree(maskImageF32); psFree(transformedImage); } - // psFree(myRegion); psFree(inCoords); psFree(outCoords); - return(NULL); + return(rejects); }