Index: /trunk/psModules/src/pmReadoutCombine.c =================================================================== --- /trunk/psModules/src/pmReadoutCombine.c (revision 2809) +++ /trunk/psModules/src/pmReadoutCombine.c (revision 2810) @@ -5,6 +5,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-21 20:41:29 $ + * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-23 20:27:51 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii Index: /trunk/psModules/src/pmSubtractSky.c =================================================================== --- /trunk/psModules/src/pmSubtractSky.c (revision 2809) +++ /trunk/psModules/src/pmSubtractSky.c (revision 2810) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-21 20:41:29 $ + * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-23 20:27:51 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -39,5 +39,5 @@ /****************************************************************************** -getHighestPriorityStatOption(statOptions): this routine takes as in put a +getHighestPriorityStatOption(statOptions): this routine takes as input a psStats->options with multiple options set and returns one with a single option set according to the precedence set in the SDRS. @@ -45,19 +45,18 @@ psU64 getHighestPriorityStatOption(psU64 statOptions) { - switch (statOptions) { - case PS_STAT_SAMPLE_MEAN: + if (statOptions & PS_STAT_SAMPLE_MEAN) { return(PS_STAT_SAMPLE_MEAN); - case PS_STAT_SAMPLE_MEDIAN: + } else if (statOptions & PS_STAT_SAMPLE_MEDIAN) { return(PS_STAT_SAMPLE_MEDIAN); - case PS_STAT_CLIPPED_MEAN: + } else if (statOptions & PS_STAT_CLIPPED_MEAN) { return(PS_STAT_CLIPPED_MEAN); - case PS_STAT_ROBUST_MEAN: + } else if (statOptions & PS_STAT_ROBUST_MEAN) { return(PS_STAT_ROBUST_MEAN); - case PS_STAT_ROBUST_MEDIAN: + } else if (statOptions & PS_STAT_ROBUST_MEDIAN) { return(PS_STAT_ROBUST_MEDIAN); - case PS_STAT_ROBUST_MODE: + } else if (statOptions & PS_STAT_ROBUST_MODE) { return(PS_STAT_ROBUST_MODE); } - // psError(XXX); + psError(PS_ERR_UNKNOWN, true, "Unallowable option requested for statistically binning image pixels.\n"); return(-1); } @@ -129,5 +128,5 @@ /****************************************************************************** CalculatePolyTerms(xOrder, yOrder): this routine will calculate the number of -coefficients (or terms) in a polynomial of order (xOrder, yOrder). +coefficients (or terms) in a 2-D polynomial of order (xOrder, yOrder). *****************************************************************************/ psS32 CalculatePolyTerms(psS32 xOrder, psS32 yOrder) @@ -187,8 +186,8 @@ } - // 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]); + if (psTraceGetLevel(".psModule.pmSubtractSky.buildPolyTerms") >= 10) { + for (i=0; i < localPolyTerms ; i++) { + printf("x^%d * y^%d\n", polyTerms[i][0], polyTerms[i][1]); + } } @@ -198,17 +197,17 @@ /****************************************************************************** This procedure calculates various combinations of powers of x and y and stores -them in the data structure sums[][]. After it completes: +them in the data structure p_psPolySums[][]. After it completes: - sums[i][j] == x^i * y^j + p_psPolySums[i][j] == x^i * y^j -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. +XXX: Use a psImage for the p_psPolySums data structure? +XXX: Check for positive x- and y-Order. +XXX: Use variable size arrays for p_psPolySums[][]. +XXX" Must initialize p_psPolySums[][]? *****************************************************************************/ #define PS_MAX_POLYNOMIAL_ORDER 20 +psF64 p_psPolySums[PS_MAX_POLYNOMIAL_ORDER+1][PS_MAX_POLYNOMIAL_ORDER+1]; void buildSums(psF64 x, psF64 y, - psF64 sums[PS_MAX_POLYNOMIAL_ORDER][PS_MAX_POLYNOMIAL_ORDER], psS32 xOrder, psS32 yOrder) @@ -224,5 +223,5 @@ ySum = xSum; for(j=0;j<=yOrder;j++) { - sums[i][j] = ySum; + p_psPolySums[i][j] = ySum; ySum*= y; } @@ -232,5 +231,5 @@ /****************************************************************************** -ImageFitPolynomial(myPoly, binnedImage, maskImage): this private routine takes +ImageFitPolynomial(myPoly, dataImage, 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 @@ -242,5 +241,5 @@ *****************************************************************************/ psPolynomial2D *ImageFitPolynomial(psPolynomial2D *myPoly, - psImage *binnedImage, + psImage *dataImage, psImage *maskImage) { @@ -249,11 +248,11 @@ PS_POLY_CHECK_NULL(myPoly, NULL); PS_POLY_CHECK_TYPE(myPoly, PS_POLYNOMIAL_ORD, 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(dataImage, NULL); + PS_IMAGE_CHECK_EMPTY(dataImage, NULL); + PS_IMAGE_CHECK_TYPE(dataImage, 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); + PS_IMAGE_CHECK_SIZE_EQUAL(dataImage, maskImage, NULL); psS32 oldPolyX = -1; psS32 oldPolyY = -1; @@ -264,15 +263,15 @@ // polynomial if there are not enough rows/columns in the input image. - if (myPoly->nX > binnedImage->numRows) { + if (myPoly->nX > dataImage->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) { + myPoly->nX = dataImage->numRows; + } + if (myPoly->nY > dataImage->numCols) { psLogMsg(__func__, PS_LOG_WARN, "WARNING: ImageFitPolynomial(): Reducing polynomial complexity in y-dimension.\n"); oldPolyY = myPoly->nY; - myPoly->nY = binnedImage->numCols; + myPoly->nY = dataImage->numCols; } psS32 i; @@ -282,9 +281,10 @@ psS32 aRow; psS32 aCol; - // XXX: Document this. The myPoly->nX and ->nY terms are actually 1 larger + // 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); + // We determine how many coefficients will be in the polynomial that we + // are fitting to this image. 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); @@ -295,23 +295,20 @@ // Initialize A matrix and B vector. // - 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; - } - + PS_IMAGE_SET_F64(A, 0.0); + PS_VECTOR_SET_F64(B, 0.0); // // We build the A matrix and B vector. // - for (x=0;xnumRows;x++) { - for (y=0;ynumCols;y++) { + for (x=0;xnumRows;x++) { + for (y=0;ynumCols;y++) { if (maskImage->data.U8[x][y] == 0) { - buildSums((psF64) x, (psF64) y, sums, myPoly->nX-1, myPoly->nY-1); + buildSums((psF64) x, (psF64) y, myPoly->nX-1, myPoly->nY-1); /************************************************************ + This code dervies from equation (7) of the pilot ADD. However, + it is not exactly the same in that the order of the polynomial + may be different in X And Y. + Equation (7) from the pilot ADD describes 16 linear equations. The i-th equation is simply the partial derivative of the @@ -325,12 +322,12 @@ for (aCol=0;aColdata.F64[aRow][aCol]+= - (sums[ polyTerms[aCol][0] ][ polyTerms[aCol][1] ] * - sums[ polyTerms[aRow][0] ][ polyTerms[aRow][1] ]); + (p_psPolySums[ polyTerms[aCol][0] ][ polyTerms[aCol][1] ] * + p_psPolySums[ polyTerms[aRow][0] ][ polyTerms[aRow][1] ]); } } // Build the B[] vector, which is the right-hand side of (7). for (i=0;idata.F64[i]+= binnedImage->data.F32[x][y] * - sums[ polyTerms[i][0] ][ polyTerms[i][1] ]; + B->data.F64[i]+= dataImage->data.F32[x][y] * + p_psPolySums[ polyTerms[i][0] ][ polyTerms[i][1] ]; } } @@ -338,21 +335,24 @@ } - // 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]); + if (psTraceGetLevel(".psModule.pmSubtractSky.ImageFitPolynomial") >= 8) { + for (aRow=0;aRowdata.F64[aRow][aCol]); + } + } + + for (i=0;i<=localPolyTerms;i++) { + printf("B[%d] is %f\n", i, B->data.F64[i]); + } } // // Solve the matrix equations for the polynomial coefficients C. + // XXX: How do we know if these matrix operations were successful? // Aout = psMatrixLUD(Aout, outPerm, A); + PS_IMAGE_CHECK_NULL(Aout, NULL); + PS_IMAGE_CHECK_EMPTY(Aout, NULL); psVector *C = psVectorAlloc(localPolyTerms, PS_TYPE_F64); psMatrixLUSolve(C, Aout, B, outPerm); @@ -382,5 +382,5 @@ // // We restore the original size of the polynomial and set remaining - // coefficients to 0.0. + // coefficients to 0.0, if necessary. // if (oldPolyX != -1) { @@ -424,4 +424,5 @@ psTrace(".psModule.pmSubtractSky", 4, "---- pmSubtractSky() begin ----\n"); + psStatsOptions statOptions = 0; // Return the original input readout if the fit specs are poorly defined. @@ -439,16 +440,18 @@ // Determine which statistic to use when binning pixels, if any. // - psStatsOptions statOptions = stats->options; - 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 (0 >= p_psDetermineNumBits(statOptions)) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractSky(): no stats->options was requested\n"); - return(in); + if (stats != NULL) { + statOptions = stats->options; + if (1 < p_psDetermineNumBits(statOptions)) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: Multiple statistical options have been requested.\n"); + statOptions = getHighestPriorityStatOption(statOptions); + + // Save old input "stats" parameter. + oldStatOptions = stats->options; + stats->options = statOptions; + } + if (0 == p_psDetermineNumBits(statOptions)) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmSubtractSky(): no stats->options was requested\n"); + } } @@ -456,8 +459,4 @@ // 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, @@ -558,5 +557,6 @@ } } else { - //XXX: error + // We shouldn't get here since we check this above. + psError(PS_ERR_UNKNOWN, true, "Unallowable fit type."); } Index: /trunk/psModules/test/tst_pmSubtractSky.c =================================================================== --- /trunk/psModules/test/tst_pmSubtractSky.c (revision 2809) +++ /trunk/psModules/test/tst_pmSubtractSky.c (revision 2810) @@ -7,6 +7,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-21 20:41:29 $ + * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-23 20:27:51 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -22,6 +22,8 @@ #define OBJECT_INTENSITY 2000.0 static int test00(void); +static int test01(void); testDescription tests[] = { {test00, 000, "pmSubtractSky", 0, false}, + {test01, 000, "pmSubtractSky: warning, error messages", 0, false}, {NULL} }; @@ -45,6 +47,6 @@ int j; int testStatus = 0; - psImage *tmpImage1 = psImageAlloc(numCols, numRows, PS_TYPE_F32); - psReadout *myReadout = psReadoutAlloc(numCols, numRows, tmpImage1); + psImage *tmpImageF32 = psImageAlloc(numCols, numRows, PS_TYPE_F32); + psReadout *myReadout = psReadoutAlloc(numCols, numRows, tmpImageF32); psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); psPolynomial2D *myPoly = psPolynomial2DAlloc(POLY_X_ORDER, POLY_Y_ORDER, PS_POLYNOMIAL_ORD); @@ -79,6 +81,6 @@ int j; int testStatus = 0; - psImage *tmpImage1 = psImageAlloc(numCols, numRows, PS_TYPE_F32); - psReadout *myReadout = psReadoutAlloc(numCols, numRows, tmpImage1); + psImage *tmpImageF32 = psImageAlloc(numCols, numRows, PS_TYPE_F32); + psReadout *myReadout = psReadoutAlloc(numCols, numRows, tmpImageF32); psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); psPolynomial2D *myPoly = psPolynomial2DAlloc(POLY_X_ORDER, POLY_Y_ORDER, PS_POLYNOMIAL_ORD); @@ -169,2 +171,171 @@ return(testStatus); } + +#define NUM_ROWS_SMALL 16 +#define NUM_COLS_SMALL 16 +int test01( void ) +{ + int testStatus = 0; + psS32 i; + psS32 j; + psImage *tmpImageF32 = psImageAlloc(NUM_COLS_SMALL, NUM_ROWS_SMALL, PS_TYPE_F32); + psImage *tmpImageF64 = psImageAlloc(NUM_COLS_SMALL, NUM_ROWS_SMALL, PS_TYPE_F64); + psReadout *myReadout = psReadoutAlloc(NUM_COLS_SMALL, NUM_ROWS_SMALL, tmpImageF32); + psReadout *rc = NULL; + psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); + psPolynomial2D *myPoly = psPolynomial2DAlloc(POLY_X_ORDER, POLY_Y_ORDER, PS_POLYNOMIAL_ORD); + for (i=0;iimage->data.F32[i][j] = func(i, j); + } + } + + printf("----------------------------------------------------------------\n"); + printf("Calling pmSubtractSky() with NULL psReadout. Should error.\n\n"); + rc = pmSubtractSky(NULL, (void *) myPoly, PM_FIT_POLYNOMIAL, 1, myStats, 2.0); + if (rc != NULL) { + printf("ERROR: pmSubtractSky() returned a non-NULL psReadout\n"); + testStatus = false; + } + + printf("----------------------------------------------------------------\n"); + printf("Calling pmSubtractSky() with NULL psReadout->image. Should error.\n\n"); + myReadout->image = NULL; + rc = pmSubtractSky(myReadout, (void *) myPoly, PM_FIT_POLYNOMIAL, 1, myStats, 2.0); + if (rc != NULL) { + printf("ERROR: pmSubtractSky() returned a non-NULL psReadout\n"); + testStatus = false; + } + myReadout->image = tmpImageF32; + + printf("----------------------------------------------------------------\n"); + printf("Calling pmSubtractSky() with PS_TYPE_F64 psReadout->image. Should error.\n\n"); + myReadout->image = tmpImageF64; + rc = pmSubtractSky(myReadout, (void *) myPoly, PM_FIT_POLYNOMIAL, 1, myStats, 2.0); + if (rc != NULL) { + printf("ERROR: pmSubtractSky() returned a non-NULL psReadout\n"); + testStatus = false; + } + myReadout->image = tmpImageF32; + + printf("----------------------------------------------------------------\n"); + printf("Calling pmSubtractSky() with NULL fitSpec. Should return image, no ERROR, no WARNING.\n\n"); + rc = pmSubtractSky(myReadout, NULL, PM_FIT_POLYNOMIAL, 1, myStats, 2.0); + if (rc != myReadout) { + printf("ERROR: pmSubtractSky() returned something other than psReadout\n"); + testStatus = false; + } + + printf("----------------------------------------------------------------\n"); + printf("Calling pmSubtractSky() with PM_FIT_NONE fit. Should return image, no ERROR, no WARNING.\n\n"); + rc = pmSubtractSky(myReadout, (void *) myPoly, PM_FIT_NONE, 1, myStats, 2.0); + if (rc != myReadout) { + printf("ERROR: pmSubtractSky() returned something other than psReadout\n"); + testStatus = false; + } + + printf("----------------------------------------------------------------\n"); + printf("Calling pmSubtractSky() with PM_FIT_SPLINE fit. Should return image, no ERROR, no WARNING.\n\n"); + rc = pmSubtractSky(myReadout, (void *) myPoly, PM_FIT_SPLINE, 1, myStats, 2.0); + if (rc != myReadout) { + printf("ERROR: pmSubtractSky() returned something other than psReadout\n"); + testStatus = false; + } + + printf("----------------------------------------------------------------\n"); + printf("Calling pmSubtractSky() with NULL myStats. Should fit entire image, generate WARNING.\n\n"); + rc = pmSubtractSky(myReadout, (void *) myPoly, PM_FIT_POLYNOMIAL, 1, NULL, 2.0); + for (i=0;iimage->data.F32[i][j])) { + printf("ERROR: image[%d][%d] is %f, should be 0.0\n", i, j, + rc->image->data.F32[i][j]); + testStatus = false; + } + } + } + + printf("----------------------------------------------------------------\n"); + psU64 oldOptions = myStats->options; + myStats->options = 0; + printf("Calling pmSubtractSky() with no myStats->options specified. Should fit entire image, generate WARNING.\n\n"); + rc = pmSubtractSky(myReadout, (void *) myPoly, PM_FIT_POLYNOMIAL, 1, myStats, 2.0); + for (i=0;iimage->data.F32[i][j])) { + printf("ERROR: image[%d][%d] is %f, should be 0.0\n", i, j, + rc->image->data.F32[i][j]); + testStatus = false; + } + } + } + myStats->options = oldOptions; + + printf("----------------------------------------------------------------\n"); + oldOptions = myStats->options; + myStats->options = PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN; + printf("Calling pmSubtractSky() with multiple myStats->options specified. Should fit entire image, generate WARNING.\n\n"); + rc = pmSubtractSky(myReadout, (void *) myPoly, PM_FIT_POLYNOMIAL, 1, myStats, 2.0); + for (i=0;iimage->data.F32[i][j])) { + printf("ERROR: image[%d][%d] is %f, should be 0.0\n", i, j, + rc->image->data.F32[i][j]); + testStatus = false; + } + } + } + myStats->options = oldOptions; + + printf("----------------------------------------------------------------\n"); + printf("Calling pmSubtractSky() with 0 binFactor. Should fit entire image, generate WARNING.\n\n"); + rc = pmSubtractSky(myReadout, (void *) myPoly, PM_FIT_POLYNOMIAL, 0, myStats, 2.0); + for (i=0;iimage->data.F32[i][j])) { + printf("ERROR: image[%d][%d] is %f, should be 0.0\n", i, j, + rc->image->data.F32[i][j]); + testStatus = false; + } + } + } + + printf("----------------------------------------------------------------\n"); + printf("Calling pmSubtractSky() with -1 binFactor. Should fit entire image, generate WARNING.\n\n"); + rc = pmSubtractSky(myReadout, (void *) myPoly, PM_FIT_POLYNOMIAL, -1, myStats, 2.0); + for (i=0;iimage->data.F32[i][j])) { + printf("ERROR: image[%d][%d] is %f, should be 0.0\n", i, j, + rc->image->data.F32[i][j]); + testStatus = false; + } + } + } + + printf("----------------------------------------------------------------\n"); + printf("Calling pmSubtractSky() with -1.0 clipSD. Should fit entire image, generate WARNING.\n\n"); + rc = pmSubtractSky(myReadout, (void *) myPoly, PM_FIT_POLYNOMIAL, 1, myStats, -1.0); + for (i=0;iimage->data.F32[i][j])) { + printf("ERROR: image[%d][%d] is %f, should be 0.0\n", i, j, + rc->image->data.F32[i][j]); + testStatus = false; + } + } + } + + + // myReadout = pmSubtractSky(myReadout, (void *) myPoly, PM_FIT_POLYNOMIAL, + // 1, myStats, 2.0); + + + printf("----------------------------------------------------------------\n"); + psFree(myReadout); + psFree(myStats); + psFree(myPoly); + psFree(tmpImageF64); + return(testStatus); +} +// This code will