Index: /trunk/psModules/src/pmSubtractSky.c =================================================================== --- /trunk/psModules/src/pmSubtractSky.c (revision 2759) +++ /trunk/psModules/src/pmSubtractSky.c (revision 2760) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-18 02:27:42 $ + * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-20 21:38:46 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -20,12 +20,14 @@ /****************************************************************************** - *****************************************************************************/ -int p_psDetermineNumBits(unsigned int data) -{ - int i; - unsigned int tmpData = data; - int numBits = 0; - - for (i=0;ioptions with multiple options set and returns one with a single +option set according to the precedence set in the SDRS. *****************************************************************************/ psU64 getHighestPriorityStatOption(psU64 statOptions) @@ -85,9 +90,9 @@ psStats *myStats = psStatsAlloc(statOptions); - for (int row = 0; row < origImage->numRows ; row+=binFactor) { - for (int col = 0; col < origImage->numCols ; col+=binFactor) { - int count = 0; - for (int binRow = 0; binRow <= binFactor ; binRow++) { - for (int binCol = 0; binCol <= binFactor ; binCol++) { + for (psS32 row = 0; row < origImage->numRows ; row+=binFactor) { + for (psS32 col = 0; col < origImage->numCols ; col+=binFactor) { + psS32 count = 0; + for (psS32 binRow = 0; binRow <= binFactor ; binRow++) { + for (psS32 binCol = 0; binCol <= binFactor ; binCol++) { if (((row + binRow) < origImage->numRows) && ((col + binCol) < origImage->numCols)) { @@ -123,4 +128,6 @@ /****************************************************************************** +CalculatePolyTerms(xOrder, yOrder): this routine will calculate the number of +coefficients (or terms) in a polynomial of order (xOrder, yOrder). *****************************************************************************/ psS32 CalculatePolyTerms(psS32 xOrder, psS32 yOrder) @@ -194,6 +201,10 @@ sums[i][j] == x^i * y^j - *****************************************************************************/ -// XXX: Use variable size arrays for polynomial sums. + +XXX: Use a psImage for the sums data structure? +XXX: Check for non-NULL sums argument? +XXX: Check for positive x- and yOrder. +XXX: Use variable size arrays for polynomial sums. + *****************************************************************************/ #define PS_MAX_POLYNOMIAL_ORDER 20 void buildSums(psF64 x, @@ -221,4 +232,12 @@ /****************************************************************************** +ImageFitPolynomial(myPoly, binnedImage, maskImage): this private routine takes +an input image along with a mask and fits a polynomial to it. The degree of +the polynomial is specified by input parameter myPoly, and need not be +symmetrical in orders of X and Y. The polynomial must be type +PS_POLYNOMIAL_ORD. If there are not enough rows or columns in the input image +for the order of the polynomial, then that order is reduced. The algorithm +used in this routine is based on that of the pilot project ADD, but is not +documented anywhere. *****************************************************************************/ psPolynomial2D *ImageFitPolynomial(psPolynomial2D *myPoly, @@ -263,5 +282,5 @@ psS32 aRow; psS32 aCol; - // XXX: Document this. The myPoly->nX and ->nY terms are actual 1 larger + // XXX: Document this. The myPoly->nX and ->nY terms are actually 1 larger // than the order of the polynomial. psS32 **polyTerms = buildPolyTerms(myPoly->nX-1, myPoly->nY-1); @@ -273,4 +292,7 @@ psVector *outPerm = psVectorAlloc(localPolyTerms, PS_TYPE_F64); + // + // Initialize A matrix and B vector. + // for(i=0;inumRows;i++) { for(j=0;jnumCols;j++) { @@ -282,4 +304,8 @@ } + + // + // We build the A matrix and B vector. + // for (x=0;xnumRows;x++) { for (y=0;ynumCols;y++) { @@ -315,5 +341,5 @@ for (aRow=0;aRowdata.F64[aRow][aCol]); } @@ -325,11 +351,14 @@ } - + // // Solve the matrix equations for the polynomial coefficients C. + // Aout = psMatrixLUD(Aout, outPerm, A); psVector *C = psVectorAlloc(localPolyTerms, PS_TYPE_F64); psMatrixLUSolve(C, Aout, B, outPerm); + // // Set the appropriate coefficients in the myPoly structure. + // for (i=0;icoeff[ polyTerms[i][0] ][ polyTerms[i][1] ] = C->data.F64[i]; @@ -338,5 +367,7 @@ } + // // Free data structures that were allocated in this module. + // for (i=0;inX = oldPolyX; @@ -382,8 +415,11 @@ void *fitSpec, psFit fit, - int binFactor, + psS32 binFactor, psStats *stats, - float clipSD) -{ + psF32 clipSD) +{ + PS_READOUT_CHECK_NULL(in, NULL); + PS_READOUT_CHECK_EMPTY(in, NULL); + PS_READOUT_CHECK_TYPE(in, PS_TYPE_F32, NULL); psTrace(".psModule.pmSubtractSky", 4, "---- pmSubtractSky() begin ----\n"); @@ -396,18 +432,23 @@ psImage *origImage = in->image; psImage *binnedImage = NULL; - psPolynomial2D *myPoly; + psPolynomial2D *myPoly = NULL; psImage *binnedMaskImage = NULL; - + psU32 oldStatOptions = 0; + + // + // Determine which statistic to use when binning pixels, if any. + // psStatsOptions statOptions = stats->options; - if (1 > p_psDetermineNumBits(statOptions)) { + if (1 < p_psDetermineNumBits(statOptions)) { //XXX psWarning(PS_ERR_UNKNOWN,true, "Multiple statistical options have been requested.\n"); statOptions = getHighestPriorityStatOption(statOptions); // XXX: Don't modify input parameter. + oldStatOptions = stats->options; stats->options = statOptions; - if (statOptions <= 0) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractSky(): unallowable stats->options was requested\n"); - return(in); - } + } + if (0 >= p_psDetermineNumBits(statOptions)) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmSubtractSky(): no stats->options was requested\n"); + return(in); } @@ -430,13 +471,24 @@ // // Bin the input image according to input parameters. - // - if ((binFactor <= 1) || - (stats == NULL) || - (0 == p_psDetermineNumBits(statOptions))) { + // Create a new binned image mask. + // + if ((binFactor <= 1) || (stats == NULL) || (0 == p_psDetermineNumBits(statOptions))) { + // No binning is required here. Simply create a copy of the image + // and a mask. binnedImage = psImageCopy(binnedImage, origImage, PS_TYPE_F32); + if (in->mask != NULL) { + binnedMaskImage = psImageCopy(binnedMaskImage, in->mask, PS_TYPE_U8); + } else { + binnedMaskImage = psImageAlloc(binnedImage->numCols, + binnedImage->numRows, + PS_TYPE_U8); + PS_IMAGE_SET_U8(binnedMaskImage, 0); + } } else { - // Add the original image mask in here. - binnedImage = psImageCopy(binnedImage, origImage, PS_TYPE_F32); - binnedImage = psImageRebin(NULL, binnedImage, NULL, 0, binFactor, stats); + binnedImage = psImageRebin(NULL, origImage, in->mask, 0, binFactor, stats); + binnedMaskImage = psImageAlloc(binnedImage->numCols, + binnedImage->numRows, + PS_TYPE_U8); + PS_IMAGE_SET_U8(binnedMaskImage, 0); } psTrace(".psModule.pmSubtractSky", 4, @@ -456,5 +508,5 @@ myStats = psImageStats(myStats, binnedImage, NULL, 0); p_psGetStatValue(myStats, &binnedMean); - psTrace(".psModule.pmSubtractSky", 8, + psTrace(".psModule.pmSubtractSky", 6, "binned Mean is %f\n", binnedMean); @@ -463,21 +515,13 @@ p_psGetStatValue(myStats, &binnedStdev); psFree(myStats); - psTrace(".psModule.pmSubtractSky", 8, + psTrace(".psModule.pmSubtractSky", 6, "binned StDev is %f\n", binnedStdev); // Clip all pixels which are more than clipSD sigmas from the mean. - // XXX: Is this correct? We simply set the mask. - // XXX: we must unset this later since we modify the image mask. - // XXX: Determine which pixels, mask or image, should be clipped. - - binnedMaskImage = psImageAlloc(binnedImage->numCols, - binnedImage->numRows, - PS_TYPE_U8); - - psTrace(".psModule.pmSubtractSky", 8, + psTrace(".psModule.pmSubtractSky", 6, "clipSD is %f\n", clipSD); - for (int row = 0; row < binnedImage->numRows ; row++) { - for (int col = 0; col < binnedImage->numCols ; col++) { + for (psS32 row = 0; row < binnedImage->numRows ; row++) { + for (psS32 col = 0; col < binnedImage->numCols ; col++) { if (fabs(binnedImage->data.F32[row][col] - binnedMean) > (clipSD * binnedStdev)) { @@ -499,14 +543,6 @@ if (myPoly != NULL) { - // XXX:Add ordinary polynomials to psImageEvalPolynomial() - for (int row = 0; row < binnedImage->numRows ; row++) { - for (int col = 0; col < binnedImage->numCols ; col++) { - binnedImage->data.F32[row][col] = - psPolynomial2DEval(myPoly, (psF32) row, (psF32) col); - psTrace(".psModule.pmSubtractSky", 8, - "binned Image[%d][%d] is %f\n", - row, col, binnedImage->data.F32[row][col]); - } - } + // Set the pixels in the binned image to that of the polynomial. + binnedImage = psImageEvalPolynomial(binnedImage, myPoly); } else { psLogMsg(__func__, PS_LOG_WARN, @@ -516,4 +552,7 @@ psFree(binnedImage); } + if (oldStatOptions != 0) { + stats->options = statOptions; + } return(in); } @@ -527,12 +566,12 @@ if (binFactor <= 1) { // The binned image is the same size as the original image. - for (int row = 0; row < origImage->numRows ; row++) { - for (int col = 0; col < origImage->numCols ; col++) { + for (psS32 row = 0; row < origImage->numRows ; row++) { + for (psS32 col = 0; col < origImage->numCols ; col++) { origImage->data.F32[row][col]-= binnedImage->data.F32[row][col]; } } } else { - for (int row = 0; row < origImage->numRows ; row++) { - for (int col = 0; col < origImage->numCols ; col++) { + for (psS32 row = 0; row < origImage->numRows ; row++) { + for (psS32 col = 0; col < origImage->numCols ; col++) { // We calculate the F32 value of the pixel coordinates in the // binned image and then use a pixel interpolation routine to @@ -548,20 +587,19 @@ psF32 binPixel = (psF32) psImagePixelInterpolate( - binnedImage, binRowF64, binColF64, + binnedImage, binColF64, binRowF64, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); origImage->data.F32[row][col]-= binPixel; psTrace(".psModule.pmSubtractSky", 8, - "image[%d][%d] <--> binnedImage[%.2f][%.2f]: %f (%f)\n", - row, col, binRowF64-0.5, binColF64-0.5, binPixel, - binnedImage->data.F32[(psS32)binRowF64][(psS32)binColF64]); - } - } - - } - + "image[%d][%d] <--> binnedImage[%.2f][%.2f]: %f\n", + row, col, binRowF64-0.5, binColF64-0.5, binPixel); + } + } + + } psFree(binnedMaskImage); - if (!((binFactor <= 1) || (stats == NULL))) { - psFree(binnedImage); + psFree(binnedImage); + if (oldStatOptions != 0) { + stats->options = statOptions; }