Index: /trunk/psLib/src/image/psImageStats.c =================================================================== --- /trunk/psLib/src/image/psImageStats.c (revision 2261) +++ /trunk/psLib/src/image/psImageStats.c (revision 2262) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.47 $ $Name: not supported by cvs2svn $ -* @date $Date: 2004-11-02 01:52:43 $ +* @version $Revision: 1.48 $ $Name: not supported by cvs2svn $ +* @date $Date: 2004-11-02 02:03:56 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -248,10 +248,18 @@ psPolynomial1DEval((float) y, chebPolys[j]); over all pixels (x,y) in the image. - -XXX: Currently we only support F32. - *****************************************************************************/ + *****************************************************************************/ psPolynomial2D* psImageFitPolynomial(psPolynomial2D* coeffs, const psImage* input) { + PS_IMAGE_CHECK_NULL(input, NULL); + PS_IMAGE_CHECK_EMPTY(input, NULL); + if ((input->type.type != PS_TYPE_S8) && + (input->type.type != PS_TYPE_U16) && + (input->type.type != PS_TYPE_F32) && + (input->type.type != PS_TYPE_F64)) { + psError(__func__, "Unalowwable image type.\n"); + } + PS_POLY_CHECK_NULL(coeffs, NULL); + PS_POLY_CHECK_TYPE(coeffs, PS_POLYNOMIAL_CHEB, NULL); psS32 x = 0; psS32 y = 0; @@ -308,8 +316,17 @@ for (x = 0; x < input->numRows; x++) { for (y = 0; y < input->numCols; y++) { - sums[i][j] += - input->data.F32[x][y] * - psPolynomial1DEval(rScalingFactors[x], chebPolys[i]) * - psPolynomial1DEval(cScalingFactors[y], chebPolys[j]); + double pixel; + if (input->type.type == PS_TYPE_S8) { + pixel = (double) input->data.S8[x][y]; + } else if (input->type.type == PS_TYPE_U16) { + pixel = (double) input->data.U16[x][y]; + } else if (input->type.type == PS_TYPE_F32) { + pixel = (double) input->data.F32[x][y]; + } else if (input->type.type == PS_TYPE_F64) { + pixel = input->data.F64[x][y]; + } + sums[i][j] += pixel * + psPolynomial1DEval(rScalingFactors[x], chebPolys[i]) * + psPolynomial1DEval(cScalingFactors[y], chebPolys[j]); } } @@ -320,5 +337,5 @@ for (j = 0; j < coeffs->nY; j++) { coeffs->coeff[i][j] = sums[i][j]; - coeffs->coeff[i][j] /= (float)(input->numRows * input->numCols); + coeffs->coeff[i][j] /= (double)(input->numRows * input->numCols); if ((i != 0) && (j != 0)) { @@ -352,7 +369,4 @@ /***************************************************************************** XXX: Use static variables for chebyshev polynomials and scaling factors. - -XXX: Currently we only support F32. - S8, U16, F32, F64 *****************************************************************************/ psImage* psImageEvalPolynomial(psImage* input, const psPolynomial2D* coeffs) @@ -373,5 +387,4 @@ psS32 i = 0; psS32 j = 0; - // float **sums = NULL; psPolynomial1D* *chebPolys = NULL; psS32 maxChebyPoly = 0; @@ -407,5 +420,14 @@ } } - input->data.F32[x][y] = polySum; + + if (input->type.type == PS_TYPE_S8) { + input->data.S8[x][y] = (char) polySum; + } else if (input->type.type == PS_TYPE_U16) { + input->data.U16[x][y] = (short int) polySum; + } else if (input->type.type == PS_TYPE_F32) { + input->data.F32[x][y] = (float) polySum; + } else if (input->type.type == PS_TYPE_F64) { + input->data.F64[x][y] = polySum; + } } } Index: /trunk/psLib/src/imageops/psImageStats.c =================================================================== --- /trunk/psLib/src/imageops/psImageStats.c (revision 2261) +++ /trunk/psLib/src/imageops/psImageStats.c (revision 2262) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.47 $ $Name: not supported by cvs2svn $ -* @date $Date: 2004-11-02 01:52:43 $ +* @version $Revision: 1.48 $ $Name: not supported by cvs2svn $ +* @date $Date: 2004-11-02 02:03:56 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -248,10 +248,18 @@ psPolynomial1DEval((float) y, chebPolys[j]); over all pixels (x,y) in the image. - -XXX: Currently we only support F32. - *****************************************************************************/ + *****************************************************************************/ psPolynomial2D* psImageFitPolynomial(psPolynomial2D* coeffs, const psImage* input) { + PS_IMAGE_CHECK_NULL(input, NULL); + PS_IMAGE_CHECK_EMPTY(input, NULL); + if ((input->type.type != PS_TYPE_S8) && + (input->type.type != PS_TYPE_U16) && + (input->type.type != PS_TYPE_F32) && + (input->type.type != PS_TYPE_F64)) { + psError(__func__, "Unalowwable image type.\n"); + } + PS_POLY_CHECK_NULL(coeffs, NULL); + PS_POLY_CHECK_TYPE(coeffs, PS_POLYNOMIAL_CHEB, NULL); psS32 x = 0; psS32 y = 0; @@ -308,8 +316,17 @@ for (x = 0; x < input->numRows; x++) { for (y = 0; y < input->numCols; y++) { - sums[i][j] += - input->data.F32[x][y] * - psPolynomial1DEval(rScalingFactors[x], chebPolys[i]) * - psPolynomial1DEval(cScalingFactors[y], chebPolys[j]); + double pixel; + if (input->type.type == PS_TYPE_S8) { + pixel = (double) input->data.S8[x][y]; + } else if (input->type.type == PS_TYPE_U16) { + pixel = (double) input->data.U16[x][y]; + } else if (input->type.type == PS_TYPE_F32) { + pixel = (double) input->data.F32[x][y]; + } else if (input->type.type == PS_TYPE_F64) { + pixel = input->data.F64[x][y]; + } + sums[i][j] += pixel * + psPolynomial1DEval(rScalingFactors[x], chebPolys[i]) * + psPolynomial1DEval(cScalingFactors[y], chebPolys[j]); } } @@ -320,5 +337,5 @@ for (j = 0; j < coeffs->nY; j++) { coeffs->coeff[i][j] = sums[i][j]; - coeffs->coeff[i][j] /= (float)(input->numRows * input->numCols); + coeffs->coeff[i][j] /= (double)(input->numRows * input->numCols); if ((i != 0) && (j != 0)) { @@ -352,7 +369,4 @@ /***************************************************************************** XXX: Use static variables for chebyshev polynomials and scaling factors. - -XXX: Currently we only support F32. - S8, U16, F32, F64 *****************************************************************************/ psImage* psImageEvalPolynomial(psImage* input, const psPolynomial2D* coeffs) @@ -373,5 +387,4 @@ psS32 i = 0; psS32 j = 0; - // float **sums = NULL; psPolynomial1D* *chebPolys = NULL; psS32 maxChebyPoly = 0; @@ -407,5 +420,14 @@ } } - input->data.F32[x][y] = polySum; + + if (input->type.type == PS_TYPE_S8) { + input->data.S8[x][y] = (char) polySum; + } else if (input->type.type == PS_TYPE_U16) { + input->data.U16[x][y] = (short int) polySum; + } else if (input->type.type == PS_TYPE_F32) { + input->data.F32[x][y] = (float) polySum; + } else if (input->type.type == PS_TYPE_F64) { + input->data.F64[x][y] = polySum; + } } } Index: /trunk/psLib/test/image/tst_psImageStats04.c =================================================================== --- /trunk/psLib/test/image/tst_psImageStats04.c (revision 2261) +++ /trunk/psLib/test/image/tst_psImageStats04.c (revision 2262) @@ -15,5 +15,5 @@ #include #include -#define IMAGE_SIZE 8 +#define IMAGE_SIZE 5 #define NUM_COLS IMAGE_SIZE #define NUM_ROWS IMAGE_SIZE @@ -29,5 +29,5 @@ ((A) + (B * X) + (C * Y) + (D * X * Y) + (E * X * X) + (F * Y * Y)) \ -psS32 FitCheby(int numCols, int numRows) +psS32 FitChebyF32(int numCols, int numRows) { psImage *tmpImage = NULL; @@ -100,10 +100,86 @@ } +psS32 FitChebyF64(int numCols, int numRows) +{ + psImage *tmpImage = NULL; + psImage *outImage = NULL; + psS32 testStatus = true; + psS32 memLeaks = 0; + psS32 i = 0; + psS32 j = 0; + psS32 currentId = psMemGetId(); + double chi2 = 0.0; + + + printf("psImageFitPolynomial(), psImageEvalPolynom(): (%d by %d)\n", numRows, numCols); + + /*************************************************************************/ + /* Allocate and initialize data structures */ + /*************************************************************************/ + + printPositiveTestHeader(stdout, + "psImageStats functions", + "psImageFitPolynomial(), psImageEvalPolynom()"); + + tmpImage = psImageAlloc(numCols, numRows, PS_TYPE_F64); + outImage = psImageAlloc(numCols, numRows, PS_TYPE_F64); + for (i=0;idata.F64[i][j] = I_POLY(((double) i), ((double) j)); + } + } + + psPolynomial2D *chebys = psPolynomial2DAlloc(5, 5, PS_POLYNOMIAL_CHEB); + + chebys = psImageFitPolynomial(chebys, tmpImage); + + outImage = psImageEvalPolynomial(outImage, chebys); + + chi2 = 0.0; + for (i=0;idata.F64[i][j]; + double actual = outImage->data.F64[i][j]; + chi2+= (actual-expect) * (actual-expect); + if (fabs((actual - expect) / expect) > ERROR_TOL) { + printf("pixel [%d][%d] is %.2f should be %.2f\n", i, j, actual, expect); + } + } + } + chi2/= ((double) (numCols * numRows)); + printf("The chi-squared per pixel is %f\n", chi2); + + psFree(tmpImage); + psFree(outImage); + psFree(chebys); + psMemCheckCorruption(1); + printFooter(stdout, + "psImageStats functions", + "psImagePixelInterpolation()", + testStatus); + + /*************************************************************************/ + /* Deallocate data structures */ + /*************************************************************************/ + psMemCheckCorruption(1); + memLeaks = psMemCheckLeaks(currentId,NULL,stderr); + if (0 != memLeaks) { + psAbort(__func__,"Memory Leaks! (%d leaks)", memLeaks); + } + + return(!testStatus); +} + int main() { - FitCheby(1, 1); - FitCheby(1, NUM_ROWS); - FitCheby(1, 1); - FitCheby(NUM_COLS, 1); - FitCheby(NUM_COLS, NUM_ROWS); + FitChebyF32(1, 1); + FitChebyF32(1, NUM_ROWS); + FitChebyF32(1, 1); + FitChebyF32(NUM_COLS, 1); + FitChebyF32(NUM_COLS, NUM_ROWS); + FitChebyF64(1, 1); + FitChebyF64(1, NUM_ROWS); + FitChebyF64(1, 1); + FitChebyF64(NUM_COLS, 1); + FitChebyF64(NUM_COLS, NUM_ROWS); }