Index: /trunk/psModules/src/pmSubtractSky.c =================================================================== --- /trunk/psModules/src/pmSubtractSky.c (revision 2754) +++ /trunk/psModules/src/pmSubtractSky.c (revision 2755) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-16 00:47:01 $ + * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-18 02:27:42 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -17,12 +17,8 @@ #include "pslib.h" #include "psConstants.h" - -// XXX: this is pmFit in pmSubtractBias.c, named psFit here. -typedef enum { - PM_FIT_NONE, ///< No fit - PM_FIT_POLYNOMIAL, ///< Fit polynomial - PM_FIT_SPLINE ///< Fit cubic splines -} psFit; - +#include "pmSubtractSky.h" + +/****************************************************************************** + *****************************************************************************/ int p_psDetermineNumBits(unsigned int data) { @@ -40,4 +36,6 @@ } +/****************************************************************************** + *****************************************************************************/ psU64 getHighestPriorityStatOption(psU64 statOptions) { @@ -60,5 +58,4 @@ } - /****************************************************************************** psImage *binImage(origImage, binFactor, statOptions): This routine takes an @@ -67,4 +64,7 @@ XXX: use static vectors for myStats, binVector and binMask. +XXX: I coded this before I was aware of a psLib reBin function. I don't +use this function in this module. I'm keeping it here in the event that +requirements change and we might need a custom reBin function. *****************************************************************************/ psImage *binImage(psImage *origImage, @@ -72,16 +72,16 @@ psStatsOptions statOptions) { + if (binFactor <= 0) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: binImage(): binFactor is %d\n", binFactor); + return(origImage); + } + if (binFactor == 1) { + return(origImage); + } + psVector *binVector = psVectorAlloc(binFactor * binFactor, PS_TYPE_F32); psVector *binMask = psVectorAlloc(binFactor * binFactor, PS_TYPE_U8); psStats *myStats = psStatsAlloc(statOptions); - - if (binFactor <= 0) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: binImage(): binFactor is %d\n", binFactor); - return(origImage); - } - if (binFactor == 1) { - return(origImage); - } for (int row = 0; row < origImage->numRows ; row+=binFactor) { @@ -115,7 +115,7 @@ } } - psFree(myStats); psFree(binVector); psFree(binMask); + psFree(myStats); return(origImage); @@ -123,8 +123,27 @@ /****************************************************************************** -buildPolyTerms(): this routine computes a 2-D array polyTerms[] that holds + *****************************************************************************/ +psS32 CalculatePolyTerms(psS32 xOrder, psS32 yOrder) +{ + psS32 maxOrder = PS_MAX(xOrder, yOrder); + psS32 localPolyTerms = 0; + psS32 order = 0; + psS32 num=0; + + for (order=0;order<=maxOrder;order++) { + for (num=0;num<=order;num++) { + if (((order-num) <= xOrder) && (num <= yOrder)) { + localPolyTerms++; + } + } + } + return(localPolyTerms); +} + +/****************************************************************************** +buildPolyTerms(): this routine computes a 2-D array polyTerms[][] that holds terms for the polynomial that is used to model the sky background. We use -this array primarily for convenience in many computations involving that sky -model polynomials. It is defined as: +this array primarily for convenience in computations involving sky model +polynomials. It is defined as: polyTerms[i][0] = the power to which X is raised in the i-th term of in an poly-order sky background polynomial. @@ -133,12 +152,13 @@ poly-order sky background polynomial. *****************************************************************************/ -psS32 **buildPolyTerms(psS32 polyOrder) +psS32 **buildPolyTerms(psS32 xOrder, psS32 yOrder) { psS32 i=0; psS32 order = 0; psS32 num=0; - psS32 localPolyTerms= (((polyOrder+1) * (polyOrder + 2)) / 2); - - // We create the data structure which we hold the xy order of each coeff. + psS32 localPolyTerms = CalculatePolyTerms(xOrder, yOrder); + psS32 maxOrder = PS_MAX(xOrder, yOrder); + + // Create the data structure which we hold the xy order of each coeff. psS32 **polyTerms = (psS32 **) psAlloc(localPolyTerms * sizeof(psS32 *)); for (i=0; i < localPolyTerms ; i++) { @@ -150,10 +170,18 @@ // calculates the power to which x/y are raised in that i-th term. // We first do the 0-order terms, then the 1-order terms, etc. - for (order=0;order<=polyOrder;order++) { + for (order=0;order<=maxOrder;order++) { for (num=0;num<=order;num++) { - polyTerms[i][0] = order-num; - polyTerms[i][1] = num; - i++; - } + if (((order-num) <= xOrder) && (num <= yOrder)) { + polyTerms[i][0] = order-num; + polyTerms[i][1] = num; + i++; + } + } + } + + // XXX: Get rid of this. + for (i=0; i < localPolyTerms ; i++) { + psTrace(".psModule.pmSubtractSky.buildPolyTerms", 6, + "x^%d * y^%d\n", polyTerms[i][0], polyTerms[i][1]); } @@ -161,15 +189,17 @@ } -#define PS_TWENTY 20 - +/****************************************************************************** +This procedure calculates various combinations of powers of x and y and stores +them in the data structure sums[][]. After it completes: + + sums[i][j] == x^i * y^j + *****************************************************************************/ // XXX: Use variable size arrays for polynomial sums. -/** @brief This procedure calculates various combinations of powers of x and y - * and stores them in the data structure sums[][]. After it completes: - * sums[i][j] == x^i * y^j - */ +#define PS_MAX_POLYNOMIAL_ORDER 20 void buildSums(psF64 x, psF64 y, - psF64 sums[PS_TWENTY][PS_TWENTY], - psS32 polyOrder) + psF64 sums[PS_MAX_POLYNOMIAL_ORDER][PS_MAX_POLYNOMIAL_ORDER], + psS32 xOrder, + psS32 yOrder) { psS32 i = 0; @@ -180,7 +210,7 @@ xSum = 1.0; ySum = 1.0; - for(i=0;i<=polyOrder;i++) { + for(i=0;i<=xOrder;i++) { ySum = xSum; - for(j=0;j<=polyOrder;j++) { + for(j=0;j<=yOrder;j++) { sums[i][j] = ySum; ySum*= y; @@ -191,5 +221,4 @@ /****************************************************************************** - *****************************************************************************/ psPolynomial2D *ImageFitPolynomial(psPolynomial2D *myPoly, @@ -197,7 +226,8 @@ psImage *maskImage) { + psTrace(".psModule.pmSubtractSky.ImageFitPolynomial", 4, + "---- ImageFitPolynomial() begin ----\n"); PS_POLY_CHECK_NULL(myPoly, NULL); PS_POLY_CHECK_TYPE(myPoly, PS_POLYNOMIAL_ORD, NULL); - PS_INT_CHECK_NON_EQUALS(myPoly->nX, myPoly->nY, NULL); PS_IMAGE_CHECK_NULL(binnedImage, NULL); PS_IMAGE_CHECK_EMPTY(binnedImage, NULL); @@ -207,4 +237,24 @@ PS_IMAGE_CHECK_TYPE(maskImage, PS_TYPE_U8, NULL); PS_IMAGE_CHECK_SIZE_EQUAL(binnedImage, maskImage, NULL); + psS32 oldPolyX = -1; + psS32 oldPolyY = -1; + + // The matrix equations become singular if there are more powers of X + // in myPoly then there are rows of the image. I think. Similarly for + // powers of Y and columns. So. Here we reduce the complexity of the + // polynomial if there are not enough rows/columns in the input image. + + if (myPoly->nX > binnedImage->numRows) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: ImageFitPolynomial(): Reducing polynomial complexity in x-dimension.\n"); + oldPolyX = myPoly->nX; + myPoly->nX = binnedImage->numRows; + } + if (myPoly->nY > binnedImage->numCols) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: ImageFitPolynomial(): Reducing polynomial complexity in y-dimension.\n"); + oldPolyY = myPoly->nY; + myPoly->nY = binnedImage->numCols; + } psS32 i; psS32 j; @@ -213,11 +263,13 @@ psS32 aRow; psS32 aCol; - psS32 polyOrder = myPoly->nX; - psS32 localPolyTerms= (((polyOrder+1) * (polyOrder + 2)) / 2); - psS32 **polyTerms = buildPolyTerms(polyOrder); - psF64 sums[PS_TWENTY][PS_TWENTY]; - psImage *A = psImageAlloc(localPolyTerms+1, localPolyTerms+1, PS_TYPE_F64); - psVector *B = psVectorAlloc(localPolyTerms+1, PS_TYPE_F64); - psVector *outPerm = NULL; + // XXX: Document this. The myPoly->nX and ->nY terms are actual 1 larger + // than the order of the polynomial. + psS32 **polyTerms = buildPolyTerms(myPoly->nX-1, myPoly->nY-1); + psS32 localPolyTerms = CalculatePolyTerms(myPoly->nX-1, myPoly->nY-1); + psF64 sums[PS_MAX_POLYNOMIAL_ORDER][PS_MAX_POLYNOMIAL_ORDER]; + psImage *A = psImageAlloc(localPolyTerms, localPolyTerms, PS_TYPE_F64); + psImage *Aout = psImageAlloc(localPolyTerms, localPolyTerms, PS_TYPE_F64); + psVector *B = psVectorAlloc(localPolyTerms, PS_TYPE_F64); + psVector *outPerm = psVectorAlloc(localPolyTerms, PS_TYPE_F64); for(i=0;inumRows;i++) { @@ -232,9 +284,9 @@ for (x=0;xnumRows;x++) { for (y=0;ynumCols;y++) { - if (maskImage->data.U8[x][y] != 0) { - buildSums((psF64) x, (psF64) y, sums, polyOrder); + if (maskImage->data.U8[x][y] == 0) { + buildSums((psF64) x, (psF64) y, sums, myPoly->nX-1, myPoly->nY-1); /************************************************************ - Equation (7) from the ADD describes 16 linear equations. + Equation (7) from the pilot ADD describes 16 linear equations. The i-th equation is simply the partial derivative of the sky background polynomial (1) w.r.t. to the i-th term in @@ -251,7 +303,6 @@ } } - // Build the B[] vector, which is the right-hand side of (7). - for (i=0;i<=localPolyTerms;i++) { + for (i=0;idata.F64[i]+= binnedImage->data.F32[x][y] * sums[ polyTerms[i][0] ][ polyTerms[i][1] ]; @@ -260,11 +311,32 @@ } } - psImage *ALUD = psMatrixLUD(NULL, outPerm, A); - psVector *C = psMatrixLUSolve(C, ALUD, B, outPerm); - + + // XXX: Put this loop inside a psTrace conditional, somehow. + for (aRow=0;aRowdata.F64[aRow][aCol]); + } + } + // XXX: Put this loop inside a psTrace conditional, somehow. + for (i=0;i<=localPolyTerms;i++) { + psTrace(".psModule.pmSubtractSky.ImageFitPolynomial", 6, + "B[%d] is %f\n", i, B->data.F64[i]); + } + + + // 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]; - } - + psTrace(".psModule.pmSubtractSky.ImageFitPolynomial", 6, + "myPoly->coeff[%d][%d] is %f\n", polyTerms[i][0], polyTerms[i][1], myPoly->coeff[ polyTerms[i][0] ][ polyTerms[i][1] ]); + } + + // Free data structures that were allocated in this module. for (i=0;inX = oldPolyX; + for (i=oldPolyX ; i < myPoly->nX ; i++) { + for (j=0;jnY ; j++) { + myPoly->coeff[i][j] = 0.0; + } + } + } + if (oldPolyY != -1) { + myPoly->nY = oldPolyY; + for (i=0 ; i < myPoly->nX ; i++) { + for (j=oldPolyY;jnY ; j++) { + myPoly->coeff[i][j] = 0.0; + } + } + } + + psTrace(".psModule.pmSubtractSky.ImageFitPolynomial", 4, + "---- ImageFitPolynomial() end successfully ----\n"); return(myPoly); } @@ -293,4 +386,7 @@ float clipSD) { + psTrace(".psModule.pmSubtractSky", 4, + "---- pmSubtractSky() begin ----\n"); + // Return the original input readout if the fit specs are poorly defined. if ((fitSpec == NULL) || @@ -299,12 +395,14 @@ } psImage *origImage = in->image; - psImage *maskImage = in->mask; psImage *binnedImage = NULL; psPolynomial2D *myPoly; + psImage *binnedMaskImage = NULL; 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. + stats->options = statOptions; if (statOptions <= 0) { psLogMsg(__func__, PS_LOG_WARN, @@ -314,22 +412,38 @@ } + // + // Generate required warning messages. + // + if (0 == p_psDetermineNumBits(statOptions)) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmSubtractSky(): no stats->options was requested\n"); + } + if (binFactor <= 0) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmSubtractSky(): binFactor is %d\n", binFactor); + } + if (stats == NULL) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmSubtractSky(): input parameter stats is NULL\n"); + } + + // // Bin the input image according to input parameters. - if ((binFactor <= 0) || (stats == NULL)) { - // Simply use the original image: no binning. - if (binFactor <= 0) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractSky(): binFactor is %d\n", binFactor); - } - if (stats == NULL) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractSky(): input parameter stats is NULL\n"); - } - binnedImage = origImage; + // + if ((binFactor <= 1) || + (stats == NULL) || + (0 == p_psDetermineNumBits(statOptions))) { + binnedImage = psImageCopy(binnedImage, origImage, PS_TYPE_F32); } else { - // Call a private function to do the binning. - binnedImage = binImage(origImage, binFactor, statOptions); - } - + // Add the original image mask in here. + binnedImage = psImageCopy(binnedImage, origImage, PS_TYPE_F32); + binnedImage = psImageRebin(NULL, binnedImage, NULL, 0, binFactor, stats); + } + psTrace(".psModule.pmSubtractSky", 4, + "binnedImage size is (%d, %d)\n", binnedImage->numRows, binnedImage->numCols); + + // // Clip pixels that are outside the acceptable range. + // if (clipSD <= 0.0) { psLogMsg(__func__, PS_LOG_WARN, @@ -342,4 +456,6 @@ myStats = psImageStats(myStats, binnedImage, NULL, 0); p_psGetStatValue(myStats, &binnedMean); + psTrace(".psModule.pmSubtractSky", 8, + "binned Mean is %f\n", binnedMean); myStats->options = PS_STAT_SAMPLE_STDEV; @@ -347,4 +463,6 @@ p_psGetStatValue(myStats, &binnedStdev); psFree(myStats); + psTrace(".psModule.pmSubtractSky", 8, + "binned StDev is %f\n", binnedStdev); // Clip all pixels which are more than clipSD sigmas from the mean. @@ -352,9 +470,17 @@ // 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, + "clipSD is %f\n", clipSD); + for (int row = 0; row < binnedImage->numRows ; row++) { for (int col = 0; col < binnedImage->numCols ; col++) { if (fabs(binnedImage->data.F32[row][col] - binnedMean) > (clipSD * binnedStdev)) { - maskImage->data.U8[row][col] = 1; + binnedMaskImage->data.U8[row][col] = 1; } } @@ -362,17 +488,42 @@ } - // XXX: fit the polynomial to the binned image + // + // Fit the polynomial to the binned image + // if (fit == PM_FIT_POLYNOMIAL) { myPoly = (psPolynomial2D *) fitSpec; - - // XXX Ensure that the polynomial is of type Chebyshev. - myPoly = ImageFitPolynomial(myPoly, binnedImage, maskImage); - - // XXX Do we need to do something with ordinate scaling if Chebyshev? - binnedImage = psImageEvalPolynomial(binnedImage, myPoly); - } - + PS_POLY_CHECK_NULL(myPoly, NULL); + PS_POLY_CHECK_TYPE(myPoly, PS_POLYNOMIAL_ORD, NULL); + + myPoly = ImageFitPolynomial(myPoly, binnedImage, binnedMaskImage); + + 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]); + } + } + } else { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmSubtractSky(): could not model sky with a polynomial.\n"); + psFree(binnedMaskImage); + if (!((binFactor <= 1) || (stats == NULL))) { + psFree(binnedImage); + } + return(in); + } + } else { + //XXX: error + } + + // //Subtract the polynomially fitted image from the original image - if (binFactor <= 0) { + // + if (binFactor <= 1) { // The binned image is the same size as the original image. for (int row = 0; row < origImage->numRows ; row++) { @@ -384,16 +535,37 @@ for (int row = 0; row < origImage->numRows ; row++) { for (int col = 0; col < origImage->numCols ; col++) { - int binRow = row / binFactor; - int binCol = col / binFactor; - origImage->data.F32[row][col]-= binnedImage->data.F32[binRow][binCol]; - } - } - - } - - if (!((binFactor <= 0) || (stats == NULL))) { + // We calculate the F32 value of the pixel coordinates in the + // binned image and then use a pixel interpolation routine to + // determine the value of the pixel at that location. + psF32 binRowF64 = ((psF32) row) / ((psF32) binFactor); + psF32 binColF64 = ((psF32) col) / ((psF32) binFactor); + + // We add 0.5 to the pixel locations since the pixel + // interpolation routine defines the location of pixel + // (i, j) as (i+0.5, j+0.5). + binRowF64+= 0.5; + binColF64+= 0.5; + + psF32 binPixel = (psF32) psImagePixelInterpolate( + binnedImage, binRowF64, binColF64, + 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]); + } + } + + } + + psFree(binnedMaskImage); + if (!((binFactor <= 1) || (stats == NULL))) { psFree(binnedImage); } + psTrace(".psModule.pmSubtractSky", 4, + "---- pmSubtractSky() exit successfully ----\n"); return(in); } Index: /trunk/psModules/src/pmSubtractSky.h =================================================================== --- /trunk/psModules/src/pmSubtractSky.h (revision 2754) +++ /trunk/psModules/src/pmSubtractSky.h (revision 2755) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-01 21:20:24 $ + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-18 02:27:42 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -26,9 +26,16 @@ #include "../../psLib/src/dataManip/psConstants.h" -psImage *pmSubtractSky(psReadout *in, - void *fitSpec, - psFit fit, - int binFactor, - psStats *stats, - float clipSD); +// XXX: this is pmFit in pmSubtractBias.c, named psFit here. +typedef enum { + PM_FIT_NONE, ///< No fit + PM_FIT_POLYNOMIAL, ///< Fit polynomial + PM_FIT_SPLINE ///< Fit cubic splines +} psFit; + +psReadout *pmSubtractSky(psReadout *in, + void *fitSpec, + psFit fit, + int binFactor, + psStats *stats, + float clipSD); #endif