Index: /trunk/psModules/src/pmSubtractSky.c =================================================================== --- /trunk/psModules/src/pmSubtractSky.c (revision 2718) +++ /trunk/psModules/src/pmSubtractSky.c (revision 2719) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-01 21:46:06 $ + * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-16 00:46:22 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -17,5 +17,5 @@ #include "psConstants.h" -// XXX: this is pmFit in pmSubtractBias.c, and psFit here. +// XXX: this is pmFit in pmSubtractBias.c, named psFit here. typedef enum { PM_FIT_NONE, ///< No fit @@ -38,6 +38,4 @@ return(numBits); } - - psU64 getHighestPriorityStatOption(psU64 statOptions) @@ -123,9 +121,167 @@ } +/****************************************************************************** +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: + polyTerms[i][0] = the power to which X is raised in the i-th term of in an + poly-order sky background polynomial. + + polyTerms[i][1] = the power to which Y is raised in the i-th term of in an + poly-order sky background polynomial. + *****************************************************************************/ +psS32 **buildPolyTerms(psS32 polyOrder) +{ + 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 **polyTerms = (psS32 **) psAlloc(localPolyTerms * sizeof(psS32 *)); + for (i=0; i < localPolyTerms ; i++) { + polyTerms[i] = (psS32 *) psAlloc(2 * sizeof(psS32)); + } + + i=0; + // This code segment loops through each term i in the polynomial and + // 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 (num=0;num<=order;num++) { + polyTerms[i][0] = order-num; + polyTerms[i][1] = num; + i++; + } + } + + return(polyTerms); +} + +#define PS_TWENTY 20 + +// 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 + */ +void buildSums(psF64 x, + psF64 y, + psF64 sums[PS_TWENTY][PS_TWENTY], + psS32 polyOrder) +{ + psS32 i = 0; + psS32 j = 0; + psF64 xSum = 0.0; + psF64 ySum = 0.0; + + xSum = 1.0; + ySum = 1.0; + for(i=0;i<=polyOrder;i++) { + ySum = xSum; + for(j=0;j<=polyOrder;j++) { + sums[i][j] = ySum; + ySum*= y; + } + xSum*= x; + } +} + +/****************************************************************************** + + *****************************************************************************/ +psPolynomial2D *ImageFitPolynomial(psPolynomial2D *myPoly, + psImage *binnedImage, + psImage *maskImage) +{ + 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); + PS_IMAGE_CHECK_TYPE(binnedImage, PS_TYPE_F32, NULL); + PS_IMAGE_CHECK_NULL(maskImage, NULL); + PS_IMAGE_CHECK_EMPTY(maskImage, NULL); + PS_IMAGE_CHECK_TYPE(maskImage, PS_TYPE_U8, NULL); + PS_IMAGE_CHECK_SIZE_EQUAL(binnedImage, maskImage, NULL); + psS32 i; + psS32 j; + psS32 x; + psS32 y; + 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; + + for(i=0;inumRows;i++) { + for(j=0;jnumCols;j++) { + A->data.F64[i][j] = 0.0; + } + } + for(i=0;in;i++) { + B->data.F64[i] = 0.0; + } + + for (x=0;xnumRows;x++) { + for (y=0;ynumCols;y++) { + if (maskImage->data.U8[x][y] != 0) { + buildSums((psF64) x, (psF64) y, sums, polyOrder); + + /************************************************************ + Equation (7) from the 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 + that polynomial. The i-th equation is stored in row i of + matrix A[][] (matrix A[][] has origin (1,1), not (0,0)). To + compute A[i][j] we simply multiply the j-th term of the Sky + Background Polynomial (SBP) by the i-th term of SBP. + ************************************************************/ + for (aRow=0;aRowdata.F64[aRow][aCol]+= + (sums[ polyTerms[aCol][0] ][ polyTerms[aCol][1] ] * + sums[ polyTerms[aRow][0] ][ polyTerms[aRow][1] ]); + } + } + + // Build the B[] vector, which is the right-hand side of (7). + for (i=0;i<=localPolyTerms;i++) { + B->data.F64[i]+= binnedImage->data.F32[x][y] * + sums[ polyTerms[i][0] ][ polyTerms[i][1] ]; + } + } + } + } + psImage *ALUD = psMatrixLUD(NULL, outPerm, A); + psVector *C = psMatrixLUSolve(C, ALUD, B, outPerm); + + for (i=0;icoeff[ polyTerms[i][0] ][ polyTerms[i][1] ] = C->data.F64[i]; + } + + for (i=0;iimage; psImage *maskImage = in->mask; psImage *binnedImage = NULL; psPolynomial2D *myPoly; - psSpline1D *mySpline; psStatsOptions statOptions = stats->options; @@ -154,9 +306,14 @@ //XXX psWarning(PS_ERR_UNKNOWN,true, "Multiple statistical options have been requested.\n"); statOptions = getHighestPriorityStatOption(statOptions); + if (statOptions <= 0) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmSubtractSky(): unallowable stats->options was requested\n"); + return(in); + } } // Bin the input image according to input parameters. - if ((binFactor <= 0) || - (stats == NULL)) { + if ((binFactor <= 0) || (stats == NULL)) { + // Simply use the original image: no binning. if (binFactor <= 0) { psLogMsg(__func__, PS_LOG_WARN, @@ -169,4 +326,5 @@ binnedImage = origImage; } else { + // Call a private function to do the binning. binnedImage = binImage(origImage, binFactor, statOptions); } @@ -175,7 +333,7 @@ if (clipSD <= 0.0) { psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractSky(): clipSD is %f\n", clipSD) - ; + "WARNING: pmSubtractSky(): clipSD is %f\n", clipSD); } else { + // Determine the mean and standard deviation of the binned image. psF64 binnedMean; psF64 binnedStdev; @@ -189,4 +347,8 @@ psFree(myStats); + // 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. for (int row = 0; row < binnedImage->numRows ; row++) { for (int col = 0; col < binnedImage->numCols ; col++) { @@ -201,19 +363,16 @@ // XXX: fit the polynomial to the binned image if (fit == PM_FIT_POLYNOMIAL) { - // Fit a polynomial to the old overscan vector. myPoly = (psPolynomial2D *) fitSpec; - myPoly = psImageFitPolynomial(myPoly, origImage); + + // 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); - - } else if (fit == PM_FIT_SPLINE) { - // Fit a spline to the old overscan vector. - mySpline = (psSpline1D *) fitSpec; - // XXX: What do we do? We don't have 2-D splines. - return(in); } //Subtract the polynomially fitted image from the original image if (binFactor <= 0) { + // 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++) { @@ -231,4 +390,9 @@ } + + if (!((binFactor <= 0) || (stats == NULL))) { + psFree(binnedImage); + } + return(in); }