Index: /trunk/psLib/src/dataManip/psFunctions.c =================================================================== --- /trunk/psLib/src/dataManip/psFunctions.c (revision 749) +++ /trunk/psLib/src/dataManip/psFunctions.c (revision 750) @@ -24,7 +24,24 @@ psGaussian(float x, float mean, - float stddev) -{ - return(exp(-((x-mean) * (x-mean)) / (2.0 * stddev * stddev))); + float stddev, + int normal) +{ + float tmp = 0.0; + + if (nomal == true) { + tmp = 1.0 / sqrtf(2.0 * PI * (stddev * stddev)); + return(tmp * exp(-((x-mean) * (x-mean)) / (2.0 * stddev * stddev))); + } else { + return(exp(-((x-mean) * (x-mean)) / (2.0 * stddev * stddev))); + } +} + +psVector *psGaussianDev(float mean, + float sigma, + int Npts) +{ + psVector *gauss = NULL; + + gauss = psVectorAlloc(Npts, PS_TYPE_FLOAT); } @@ -204,5 +221,5 @@ /***************************************************************************** - Polynomial coefficients will be accessed in (x,y,z) fashion. + Polynomial coefficients will be accessed in [w][x][y][z] fashion. *****************************************************************************/ float @@ -212,7 +229,9 @@ int loop_x = 0; float polySum = 0.0; + float xSum = 1.0; for (loop_x=0;loop_xn;loop_x++) { - polySum+= x + myPoly->coeff[loop_x]; + polySum+= xSum * myPoly->coeff[loop_x]; + xSum*=x; } @@ -229,9 +248,14 @@ int loop_y = 0; float polySum = 0.0; + float xSum = 1.0; + float ySum = 1.0; for (loop_x=0;loop_xnX;loop_x++) { + xSum = xSum; for (loop_y=0;loop_ynY;loop_y++) { - polySum+= x + myPoly->coeff[loop_x][loop_y]; - } + polySum+= ySum * myPoly->coeff[loop_x][loop_y]; + ySum*=y; + } + xSum*=x; } @@ -249,11 +273,19 @@ int loop_z = 0; float polySum = 0.0; + float xSum = 1.0; + float ySum = 1.0; + float zSum = 1.0; for (loop_x=0;loop_xnX;loop_x++) { + ySum = xSum; for (loop_y=0;loop_ynY;loop_y++) { + zSum = ySum; for (loop_z=0;loop_znZ;loop_z++) { - polySum+= x + myPoly->coeff[loop_x][loop_y][loop_z]; - } - } + polySum+= zSum * myPoly->coeff[loop_x][loop_y][loop_z]; + zSum*=z; + } + ySum*=y; + } + xSum*=x; } @@ -273,13 +305,25 @@ int loop_z = 0; float polySum = 0.0; + float wSum = 1.0; + float xSum = 1.0; + float ySum = 1.0; + float zSum = 1.0; for (loop_w=0;loop_wnW;loop_w++) { + xSum = wSum; for (loop_x=0;loop_xnX;loop_x++) { + ySum = xSum; for (loop_y=0;loop_ynY;loop_y++) { + zSum = ySum; for (loop_z=0;loop_znZ;loop_z++) { - polySum+= x + myPoly->coeff[loop_w][loop_x][loop_y][loop_z]; + polySum+= zSum * + myPoly->coeff[loop_w][loop_x][loop_y][loop_z]; + zSum*=z; } - } - } + ySum*=y; + } + xSum*=x; + } + wSum*=w; } @@ -461,13 +505,18 @@ } +/***************************************************************************** + Polynomial coefficients will be accessed in [w][x][y][z] fashion. + *****************************************************************************/ double -psevalDPolynomial1D(double x, - const psPolynomial1D *myPoly) +psDEvalPolynomial1D(double x, + const psDPolynomial1D *myPoly) { int loop_x = 0; double polySum = 0.0; + double xSum = 1.0; for (loop_x=0;loop_xn;loop_x++) { - polySum+= x + myPoly->coeff[loop_x]; + polySum+= xSum * myPoly->coeff[loop_x]; + xSum*=x; } @@ -477,16 +526,21 @@ double -psevalDPolynomial2D(double x, +psDEvalPolynomial2D(double x, double y, - const psPolynomial2D *myPoly) + const psDPolynomial2D *myPoly) { int loop_x = 0; int loop_y = 0; double polySum = 0.0; + double xSum = 1.0; + double ySum = 1.0; for (loop_x=0;loop_xnX;loop_x++) { + xSum = xSum; for (loop_y=0;loop_ynY;loop_y++) { - polySum+= x + myPoly->coeff[loop_x][loop_y]; - } + polySum+= ySum * myPoly->coeff[loop_x][loop_y]; + ySum*=y; + } + xSum*=x; } @@ -495,8 +549,8 @@ double -psevalDPolynomial3D(double x, +psDEvalPolynomial3D(double x, double y, double z, - const psPolynomial3D *myPoly) + const psDPolynomial3D *myPoly) { int loop_x = 0; @@ -504,11 +558,19 @@ int loop_z = 0; double polySum = 0.0; + double xSum = 1.0; + double ySum = 1.0; + double zSum = 1.0; for (loop_x=0;loop_xnX;loop_x++) { + ySum = xSum; for (loop_y=0;loop_ynY;loop_y++) { + zSum = ySum; for (loop_z=0;loop_znZ;loop_z++) { - polySum+= x + myPoly->coeff[loop_x][loop_y][loop_z]; - } - } + polySum+= zSum * myPoly->coeff[loop_x][loop_y][loop_z]; + zSum*=z; + } + ySum*=y; + } + xSum*=x; } @@ -517,9 +579,9 @@ double -psevalDPolynomial4D(double w, +psDEvalPolynomial4D(double w, double x, double y, double z, - const psPolynomial4D *myPoly) + const psDPolynomial4D *myPoly) { int loop_w = 0; @@ -528,15 +590,27 @@ int loop_z = 0; double polySum = 0.0; + double wSum = 1.0; + double xSum = 1.0; + double ySum = 1.0; + double zSum = 1.0; for (loop_w=0;loop_wnW;loop_w++) { + xSum = wSum; for (loop_x=0;loop_xnX;loop_x++) { + ySum = xSum; for (loop_y=0;loop_ynY;loop_y++) { + zSum = ySum; for (loop_z=0;loop_znZ;loop_z++) { - polySum+= x + myPoly->coeff[loop_w][loop_x][loop_y][loop_z]; + polySum+= zSum * + myPoly->coeff[loop_w][loop_x][loop_y][loop_z]; + zSum*=z; } - } - } - } - - return(polySum); -} + ySum*=y; + } + xSum*=x; + } + wSum*=w; + } + + return(polySum); +} Index: /trunk/psLib/src/dataManip/psFunctions.h =================================================================== --- /trunk/psLib/src/dataManip/psFunctions.h (revision 749) +++ /trunk/psLib/src/dataManip/psFunctions.h (revision 750) @@ -91,4 +91,5 @@ void psPolynomial4DFree(psPolynomial4D *myPoly ///< Polynomial to destroy ); + /** Evaluate 1D polynomial */ Index: /trunk/psLib/src/dataManip/psMinimize.c =================================================================== --- /trunk/psLib/src/dataManip/psMinimize.c (revision 749) +++ /trunk/psLib/src/dataManip/psMinimize.c (revision 750) @@ -13,8 +13,10 @@ #include "psFunctions.h" #include "psSort.h" - #include "float.h" #include - +// GUS: rewrite so there is no maximum order for the polynomials. +#define MAX_POLY_ORDER 10 +#define MAX_POLYNOMIAL_TERMS (((MAX_POLY_ORDER+1) * (MAX_POLY_ORDER + 2)) / 2) +int MyInfoLevel = 0; /****************************************************************************** This routine must minimize an arbitrary function. @@ -43,6 +45,247 @@ } +/** @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(double x, + double y, + /*@out@*/double sums[MAX_POLY_ORDER+1][MAX_POLY_ORDER+1], + int polyOrder) +{ + int i = 0; // loop index variable + int j = 0; // loop index variable + double xSum = 0.0; // The running sum of X terms + double ySum = 0.0; // The running sum of Y terms + + 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; + } +} + +/** @brief The coefficients of the matrix in equation (7) from the ADD will + * be very large if the x and y values are in the 0-511 range (ie: the sum y^7 + * for all 0 + *
  • height + *
  • width + * + * + * Output: + *
      + *
    • imageScalingFactors + *
    + * + * @return error status (PsError) indicating error information, or NULL on + * success. + */ +void buildImageScalingFactors(int height, + int width, + float **imageScalingFactors) +{ + int maxDim = 0; // The largest dimension of the image. + int i = 0; // loop index variable. + + // Calculate the maximum dimensional extent of the image. + if (height > width) { + maxDim = height; + } else { + maxDim = width; + } + + + // Allocate memory for the output array. + *imageScalingFactors = (float *) psAlloc((maxDim+10) * sizeof(float)); + + // This code is somewhat arbitrary. For an image with a height/width + // of 512x512, the scaling factors will be between 0.0-1.0. + for (i=0;i. + * polyTerms[poly][i][1] = the power to which Y is raised in the + * i-th term of in an poly-order sky + * background polynomial

    . + * + * NOTE: the C_0 term defined in the ADD begins at i=2 in our data + * structures (ie. the x/y powers of the i-th term in the sky model + * polynomial are actually stored at polyTerms[][i+2][]. There are two + * reasons for this. First, there is a term prior to C_0 in equation + * (7) of the ADD. Second, our linear algebra codes assume data is + * stored offset from index 1. + * + * Input: + *
      + *
    • polyTerms[][][] + *
    + * + * Output: + *
      + *
    • polyTerms[][][] + *
    + * + * @return error status (PsError) indicating error information, or NULL on + * success. + */ +void buildPolyTerms(/*@out@*/ int polyTerms[MAX_POLY_ORDER+1][(MAX_POLYNOMIAL_TERMS+2)][2]) +{ + int polyOrder=0; // loop index variable. + int i=0; // loop index variable. + int term = 0; // loop index variable. + int num=0; // loop index variable. + + for(polyOrder=0;polyOrder<=MAX_POLY_ORDER;polyOrder++) { + // The following 4 terms should not be used in any of the subsequent + // computation. We initialize them to zero in order to produce stable + // results for debugging purposes should they mistakenly be used. + polyTerms[polyOrder][0][0] = 0; + polyTerms[polyOrder][0][1] = 0; + polyTerms[polyOrder][1][0] = 0; + polyTerms[polyOrder][1][1] = 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. + i=2; + for (term=0;term<=polyOrder;term++) { + for (num=0;num<=term;num++) { + polyTerms[polyOrder][i][0] = term-num; + polyTerms[polyOrder][i][1] = num; + if (MyInfoLevel > 2) { + printf("%d-th order Sky polynomial term %d is x^%d y^%d\n", + polyOrder, i, + polyTerms[polyOrder][i][0], polyTerms[polyOrder][i][1]); + } + i++; + } + } + } +} + + +/** @brief This routine checks if all polyOrder-th terms in the polyOrder-th + * order sky background polynomial defined by the coefficients in the array B[] + * are consistent with zero. If true, then *flag is set to 1. Otherwise, + * *flag is set to 0. The matrix inversion code in the middle of this + * procedure draws from Numerical Recipes in C page 48. + * + * Input: + *
      + *
    • A This is the LUD decomposition of the original matrix A. + *
    • N The size of the matrix (plus 1, actually, since offset 1). + *
    • indx misc Numerical Recipes data structure. + *
    • B The coefficients of the sky polynomial. + *
    • polyOrder The degree of the sky polynomial. + *
    + * Output: + *
      + *
    • *flag Set this to 1 if we must recalculate the coefficients. + *
    + * + * @return error status (PsError) indicating error information, or NULL on + * success. + */ +void polyOrderCheck(float **A, + int N, + int *indx, + float *B, + int polyOrder, + int *flag) +{ + float **y = NULL; // This 2-D matrix will hold A^-1 + float *col = NULL; // misc NumerRecipes data structure + float *error=NULL; // will hold the sqrt() of the + // diagonal of y[][]. + int i=0; // loop-index variable + int j=0; // loop-index variable + int numPolyTerms = 0; // The number of terms in the + // polynomial. + int lastTerm = 0; // The index location of the first + // n-th order term in array B[]. + int firstTerm = 0; // Index location of last such term. + + // Allocate the necessary data structures for this procedure... + error = (float *) psAlloc((N + 1) * sizeof(float)); + col = (float *) psAlloc((N + 1) * sizeof(float)); + y = (float **) psAlloc((N + 1) * sizeof(float *)); + for(i=1;i<=N;i++) { + y[i] = (float *) psAlloc((N + 1) * sizeof(float)); + } + + // Invert the matrix A and put the result in y[][]. This code is taken + // from Numerical Recipes in C page 48. + for(j=1;j<=N;j++) { + for(i=1;i<=N;i++) { + col[i] = 0.0; + } + col[j] = 1.0; + // GUS: substitue the LUD rotine + // lubksb(A, N, indx, col); + for(i=1;i<=N;i++) { + y[i][j] = col[i]; + } + } + + // Determine where the first n-th order (in this comment, n equals + // polyOrder) polynomial term is stored in the matrix B[], and also were + // the last n-order term is stored. Then we loop over all the n-order + // terms and check if they are consistent with zero. + + numPolyTerms = (((polyOrder+1) * (polyOrder + 2)) / 2); + lastTerm = numPolyTerms + 1; + firstTerm = lastTerm - polyOrder; + *flag = 1; + for (i=firstTerm; i<=lastTerm; i++) { + error[i] = sqrtf(y[i][i]); + if (!((B[i] <= (2.0f * error[i])) && + ((-2.0f * error[i]) <= B[i]))) { + *flag = 0; + } + } + + // Free all memory allocated in this routine. + psFree(error); + psFree(col); + for(j=1;j<=N;j++) { + psFree(y[j]); + } + psFree(y); +} + + + + + + + + + /****************************************************************************** - This routine must minimize an arbotrary function. + This routine must fit a polynomial of degree myPoly to the data points + (x, y) and return the coefficients of that polynomial, as well as the + error for each data poiny (yErr). *****************************************************************************/ psPolynomial1D * Index: /trunk/psLib/src/math/psMinimize.c =================================================================== --- /trunk/psLib/src/math/psMinimize.c (revision 749) +++ /trunk/psLib/src/math/psMinimize.c (revision 750) @@ -13,8 +13,10 @@ #include "psFunctions.h" #include "psSort.h" - #include "float.h" #include - +// GUS: rewrite so there is no maximum order for the polynomials. +#define MAX_POLY_ORDER 10 +#define MAX_POLYNOMIAL_TERMS (((MAX_POLY_ORDER+1) * (MAX_POLY_ORDER + 2)) / 2) +int MyInfoLevel = 0; /****************************************************************************** This routine must minimize an arbitrary function. @@ -43,6 +45,247 @@ } +/** @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(double x, + double y, + /*@out@*/double sums[MAX_POLY_ORDER+1][MAX_POLY_ORDER+1], + int polyOrder) +{ + int i = 0; // loop index variable + int j = 0; // loop index variable + double xSum = 0.0; // The running sum of X terms + double ySum = 0.0; // The running sum of Y terms + + 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; + } +} + +/** @brief The coefficients of the matrix in equation (7) from the ADD will + * be very large if the x and y values are in the 0-511 range (ie: the sum y^7 + * for all 0 + *
  • height + *
  • width + * + * + * Output: + *
      + *
    • imageScalingFactors + *
    + * + * @return error status (PsError) indicating error information, or NULL on + * success. + */ +void buildImageScalingFactors(int height, + int width, + float **imageScalingFactors) +{ + int maxDim = 0; // The largest dimension of the image. + int i = 0; // loop index variable. + + // Calculate the maximum dimensional extent of the image. + if (height > width) { + maxDim = height; + } else { + maxDim = width; + } + + + // Allocate memory for the output array. + *imageScalingFactors = (float *) psAlloc((maxDim+10) * sizeof(float)); + + // This code is somewhat arbitrary. For an image with a height/width + // of 512x512, the scaling factors will be between 0.0-1.0. + for (i=0;i. + * polyTerms[poly][i][1] = the power to which Y is raised in the + * i-th term of in an poly-order sky + * background polynomial

    . + * + * NOTE: the C_0 term defined in the ADD begins at i=2 in our data + * structures (ie. the x/y powers of the i-th term in the sky model + * polynomial are actually stored at polyTerms[][i+2][]. There are two + * reasons for this. First, there is a term prior to C_0 in equation + * (7) of the ADD. Second, our linear algebra codes assume data is + * stored offset from index 1. + * + * Input: + *
      + *
    • polyTerms[][][] + *
    + * + * Output: + *
      + *
    • polyTerms[][][] + *
    + * + * @return error status (PsError) indicating error information, or NULL on + * success. + */ +void buildPolyTerms(/*@out@*/ int polyTerms[MAX_POLY_ORDER+1][(MAX_POLYNOMIAL_TERMS+2)][2]) +{ + int polyOrder=0; // loop index variable. + int i=0; // loop index variable. + int term = 0; // loop index variable. + int num=0; // loop index variable. + + for(polyOrder=0;polyOrder<=MAX_POLY_ORDER;polyOrder++) { + // The following 4 terms should not be used in any of the subsequent + // computation. We initialize them to zero in order to produce stable + // results for debugging purposes should they mistakenly be used. + polyTerms[polyOrder][0][0] = 0; + polyTerms[polyOrder][0][1] = 0; + polyTerms[polyOrder][1][0] = 0; + polyTerms[polyOrder][1][1] = 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. + i=2; + for (term=0;term<=polyOrder;term++) { + for (num=0;num<=term;num++) { + polyTerms[polyOrder][i][0] = term-num; + polyTerms[polyOrder][i][1] = num; + if (MyInfoLevel > 2) { + printf("%d-th order Sky polynomial term %d is x^%d y^%d\n", + polyOrder, i, + polyTerms[polyOrder][i][0], polyTerms[polyOrder][i][1]); + } + i++; + } + } + } +} + + +/** @brief This routine checks if all polyOrder-th terms in the polyOrder-th + * order sky background polynomial defined by the coefficients in the array B[] + * are consistent with zero. If true, then *flag is set to 1. Otherwise, + * *flag is set to 0. The matrix inversion code in the middle of this + * procedure draws from Numerical Recipes in C page 48. + * + * Input: + *
      + *
    • A This is the LUD decomposition of the original matrix A. + *
    • N The size of the matrix (plus 1, actually, since offset 1). + *
    • indx misc Numerical Recipes data structure. + *
    • B The coefficients of the sky polynomial. + *
    • polyOrder The degree of the sky polynomial. + *
    + * Output: + *
      + *
    • *flag Set this to 1 if we must recalculate the coefficients. + *
    + * + * @return error status (PsError) indicating error information, or NULL on + * success. + */ +void polyOrderCheck(float **A, + int N, + int *indx, + float *B, + int polyOrder, + int *flag) +{ + float **y = NULL; // This 2-D matrix will hold A^-1 + float *col = NULL; // misc NumerRecipes data structure + float *error=NULL; // will hold the sqrt() of the + // diagonal of y[][]. + int i=0; // loop-index variable + int j=0; // loop-index variable + int numPolyTerms = 0; // The number of terms in the + // polynomial. + int lastTerm = 0; // The index location of the first + // n-th order term in array B[]. + int firstTerm = 0; // Index location of last such term. + + // Allocate the necessary data structures for this procedure... + error = (float *) psAlloc((N + 1) * sizeof(float)); + col = (float *) psAlloc((N + 1) * sizeof(float)); + y = (float **) psAlloc((N + 1) * sizeof(float *)); + for(i=1;i<=N;i++) { + y[i] = (float *) psAlloc((N + 1) * sizeof(float)); + } + + // Invert the matrix A and put the result in y[][]. This code is taken + // from Numerical Recipes in C page 48. + for(j=1;j<=N;j++) { + for(i=1;i<=N;i++) { + col[i] = 0.0; + } + col[j] = 1.0; + // GUS: substitue the LUD rotine + // lubksb(A, N, indx, col); + for(i=1;i<=N;i++) { + y[i][j] = col[i]; + } + } + + // Determine where the first n-th order (in this comment, n equals + // polyOrder) polynomial term is stored in the matrix B[], and also were + // the last n-order term is stored. Then we loop over all the n-order + // terms and check if they are consistent with zero. + + numPolyTerms = (((polyOrder+1) * (polyOrder + 2)) / 2); + lastTerm = numPolyTerms + 1; + firstTerm = lastTerm - polyOrder; + *flag = 1; + for (i=firstTerm; i<=lastTerm; i++) { + error[i] = sqrtf(y[i][i]); + if (!((B[i] <= (2.0f * error[i])) && + ((-2.0f * error[i]) <= B[i]))) { + *flag = 0; + } + } + + // Free all memory allocated in this routine. + psFree(error); + psFree(col); + for(j=1;j<=N;j++) { + psFree(y[j]); + } + psFree(y); +} + + + + + + + + + /****************************************************************************** - This routine must minimize an arbotrary function. + This routine must fit a polynomial of degree myPoly to the data points + (x, y) and return the coefficients of that polynomial, as well as the + error for each data poiny (yErr). *****************************************************************************/ psPolynomial1D * Index: /trunk/psLib/src/math/psPolynomial.c =================================================================== --- /trunk/psLib/src/math/psPolynomial.c (revision 749) +++ /trunk/psLib/src/math/psPolynomial.c (revision 750) @@ -24,7 +24,24 @@ psGaussian(float x, float mean, - float stddev) -{ - return(exp(-((x-mean) * (x-mean)) / (2.0 * stddev * stddev))); + float stddev, + int normal) +{ + float tmp = 0.0; + + if (nomal == true) { + tmp = 1.0 / sqrtf(2.0 * PI * (stddev * stddev)); + return(tmp * exp(-((x-mean) * (x-mean)) / (2.0 * stddev * stddev))); + } else { + return(exp(-((x-mean) * (x-mean)) / (2.0 * stddev * stddev))); + } +} + +psVector *psGaussianDev(float mean, + float sigma, + int Npts) +{ + psVector *gauss = NULL; + + gauss = psVectorAlloc(Npts, PS_TYPE_FLOAT); } @@ -204,5 +221,5 @@ /***************************************************************************** - Polynomial coefficients will be accessed in (x,y,z) fashion. + Polynomial coefficients will be accessed in [w][x][y][z] fashion. *****************************************************************************/ float @@ -212,7 +229,9 @@ int loop_x = 0; float polySum = 0.0; + float xSum = 1.0; for (loop_x=0;loop_xn;loop_x++) { - polySum+= x + myPoly->coeff[loop_x]; + polySum+= xSum * myPoly->coeff[loop_x]; + xSum*=x; } @@ -229,9 +248,14 @@ int loop_y = 0; float polySum = 0.0; + float xSum = 1.0; + float ySum = 1.0; for (loop_x=0;loop_xnX;loop_x++) { + xSum = xSum; for (loop_y=0;loop_ynY;loop_y++) { - polySum+= x + myPoly->coeff[loop_x][loop_y]; - } + polySum+= ySum * myPoly->coeff[loop_x][loop_y]; + ySum*=y; + } + xSum*=x; } @@ -249,11 +273,19 @@ int loop_z = 0; float polySum = 0.0; + float xSum = 1.0; + float ySum = 1.0; + float zSum = 1.0; for (loop_x=0;loop_xnX;loop_x++) { + ySum = xSum; for (loop_y=0;loop_ynY;loop_y++) { + zSum = ySum; for (loop_z=0;loop_znZ;loop_z++) { - polySum+= x + myPoly->coeff[loop_x][loop_y][loop_z]; - } - } + polySum+= zSum * myPoly->coeff[loop_x][loop_y][loop_z]; + zSum*=z; + } + ySum*=y; + } + xSum*=x; } @@ -273,13 +305,25 @@ int loop_z = 0; float polySum = 0.0; + float wSum = 1.0; + float xSum = 1.0; + float ySum = 1.0; + float zSum = 1.0; for (loop_w=0;loop_wnW;loop_w++) { + xSum = wSum; for (loop_x=0;loop_xnX;loop_x++) { + ySum = xSum; for (loop_y=0;loop_ynY;loop_y++) { + zSum = ySum; for (loop_z=0;loop_znZ;loop_z++) { - polySum+= x + myPoly->coeff[loop_w][loop_x][loop_y][loop_z]; + polySum+= zSum * + myPoly->coeff[loop_w][loop_x][loop_y][loop_z]; + zSum*=z; } - } - } + ySum*=y; + } + xSum*=x; + } + wSum*=w; } @@ -461,13 +505,18 @@ } +/***************************************************************************** + Polynomial coefficients will be accessed in [w][x][y][z] fashion. + *****************************************************************************/ double -psevalDPolynomial1D(double x, - const psPolynomial1D *myPoly) +psDEvalPolynomial1D(double x, + const psDPolynomial1D *myPoly) { int loop_x = 0; double polySum = 0.0; + double xSum = 1.0; for (loop_x=0;loop_xn;loop_x++) { - polySum+= x + myPoly->coeff[loop_x]; + polySum+= xSum * myPoly->coeff[loop_x]; + xSum*=x; } @@ -477,16 +526,21 @@ double -psevalDPolynomial2D(double x, +psDEvalPolynomial2D(double x, double y, - const psPolynomial2D *myPoly) + const psDPolynomial2D *myPoly) { int loop_x = 0; int loop_y = 0; double polySum = 0.0; + double xSum = 1.0; + double ySum = 1.0; for (loop_x=0;loop_xnX;loop_x++) { + xSum = xSum; for (loop_y=0;loop_ynY;loop_y++) { - polySum+= x + myPoly->coeff[loop_x][loop_y]; - } + polySum+= ySum * myPoly->coeff[loop_x][loop_y]; + ySum*=y; + } + xSum*=x; } @@ -495,8 +549,8 @@ double -psevalDPolynomial3D(double x, +psDEvalPolynomial3D(double x, double y, double z, - const psPolynomial3D *myPoly) + const psDPolynomial3D *myPoly) { int loop_x = 0; @@ -504,11 +558,19 @@ int loop_z = 0; double polySum = 0.0; + double xSum = 1.0; + double ySum = 1.0; + double zSum = 1.0; for (loop_x=0;loop_xnX;loop_x++) { + ySum = xSum; for (loop_y=0;loop_ynY;loop_y++) { + zSum = ySum; for (loop_z=0;loop_znZ;loop_z++) { - polySum+= x + myPoly->coeff[loop_x][loop_y][loop_z]; - } - } + polySum+= zSum * myPoly->coeff[loop_x][loop_y][loop_z]; + zSum*=z; + } + ySum*=y; + } + xSum*=x; } @@ -517,9 +579,9 @@ double -psevalDPolynomial4D(double w, +psDEvalPolynomial4D(double w, double x, double y, double z, - const psPolynomial4D *myPoly) + const psDPolynomial4D *myPoly) { int loop_w = 0; @@ -528,15 +590,27 @@ int loop_z = 0; double polySum = 0.0; + double wSum = 1.0; + double xSum = 1.0; + double ySum = 1.0; + double zSum = 1.0; for (loop_w=0;loop_wnW;loop_w++) { + xSum = wSum; for (loop_x=0;loop_xnX;loop_x++) { + ySum = xSum; for (loop_y=0;loop_ynY;loop_y++) { + zSum = ySum; for (loop_z=0;loop_znZ;loop_z++) { - polySum+= x + myPoly->coeff[loop_w][loop_x][loop_y][loop_z]; + polySum+= zSum * + myPoly->coeff[loop_w][loop_x][loop_y][loop_z]; + zSum*=z; } - } - } - } - - return(polySum); -} + ySum*=y; + } + xSum*=x; + } + wSum*=w; + } + + return(polySum); +} Index: /trunk/psLib/src/math/psPolynomial.h =================================================================== --- /trunk/psLib/src/math/psPolynomial.h (revision 749) +++ /trunk/psLib/src/math/psPolynomial.h (revision 750) @@ -91,4 +91,5 @@ void psPolynomial4DFree(psPolynomial4D *myPoly ///< Polynomial to destroy ); + /** Evaluate 1D polynomial */ Index: /trunk/psLib/src/math/psSpline.c =================================================================== --- /trunk/psLib/src/math/psSpline.c (revision 749) +++ /trunk/psLib/src/math/psSpline.c (revision 750) @@ -24,7 +24,24 @@ psGaussian(float x, float mean, - float stddev) -{ - return(exp(-((x-mean) * (x-mean)) / (2.0 * stddev * stddev))); + float stddev, + int normal) +{ + float tmp = 0.0; + + if (nomal == true) { + tmp = 1.0 / sqrtf(2.0 * PI * (stddev * stddev)); + return(tmp * exp(-((x-mean) * (x-mean)) / (2.0 * stddev * stddev))); + } else { + return(exp(-((x-mean) * (x-mean)) / (2.0 * stddev * stddev))); + } +} + +psVector *psGaussianDev(float mean, + float sigma, + int Npts) +{ + psVector *gauss = NULL; + + gauss = psVectorAlloc(Npts, PS_TYPE_FLOAT); } @@ -204,5 +221,5 @@ /***************************************************************************** - Polynomial coefficients will be accessed in (x,y,z) fashion. + Polynomial coefficients will be accessed in [w][x][y][z] fashion. *****************************************************************************/ float @@ -212,7 +229,9 @@ int loop_x = 0; float polySum = 0.0; + float xSum = 1.0; for (loop_x=0;loop_xn;loop_x++) { - polySum+= x + myPoly->coeff[loop_x]; + polySum+= xSum * myPoly->coeff[loop_x]; + xSum*=x; } @@ -229,9 +248,14 @@ int loop_y = 0; float polySum = 0.0; + float xSum = 1.0; + float ySum = 1.0; for (loop_x=0;loop_xnX;loop_x++) { + xSum = xSum; for (loop_y=0;loop_ynY;loop_y++) { - polySum+= x + myPoly->coeff[loop_x][loop_y]; - } + polySum+= ySum * myPoly->coeff[loop_x][loop_y]; + ySum*=y; + } + xSum*=x; } @@ -249,11 +273,19 @@ int loop_z = 0; float polySum = 0.0; + float xSum = 1.0; + float ySum = 1.0; + float zSum = 1.0; for (loop_x=0;loop_xnX;loop_x++) { + ySum = xSum; for (loop_y=0;loop_ynY;loop_y++) { + zSum = ySum; for (loop_z=0;loop_znZ;loop_z++) { - polySum+= x + myPoly->coeff[loop_x][loop_y][loop_z]; - } - } + polySum+= zSum * myPoly->coeff[loop_x][loop_y][loop_z]; + zSum*=z; + } + ySum*=y; + } + xSum*=x; } @@ -273,13 +305,25 @@ int loop_z = 0; float polySum = 0.0; + float wSum = 1.0; + float xSum = 1.0; + float ySum = 1.0; + float zSum = 1.0; for (loop_w=0;loop_wnW;loop_w++) { + xSum = wSum; for (loop_x=0;loop_xnX;loop_x++) { + ySum = xSum; for (loop_y=0;loop_ynY;loop_y++) { + zSum = ySum; for (loop_z=0;loop_znZ;loop_z++) { - polySum+= x + myPoly->coeff[loop_w][loop_x][loop_y][loop_z]; + polySum+= zSum * + myPoly->coeff[loop_w][loop_x][loop_y][loop_z]; + zSum*=z; } - } - } + ySum*=y; + } + xSum*=x; + } + wSum*=w; } @@ -461,13 +505,18 @@ } +/***************************************************************************** + Polynomial coefficients will be accessed in [w][x][y][z] fashion. + *****************************************************************************/ double -psevalDPolynomial1D(double x, - const psPolynomial1D *myPoly) +psDEvalPolynomial1D(double x, + const psDPolynomial1D *myPoly) { int loop_x = 0; double polySum = 0.0; + double xSum = 1.0; for (loop_x=0;loop_xn;loop_x++) { - polySum+= x + myPoly->coeff[loop_x]; + polySum+= xSum * myPoly->coeff[loop_x]; + xSum*=x; } @@ -477,16 +526,21 @@ double -psevalDPolynomial2D(double x, +psDEvalPolynomial2D(double x, double y, - const psPolynomial2D *myPoly) + const psDPolynomial2D *myPoly) { int loop_x = 0; int loop_y = 0; double polySum = 0.0; + double xSum = 1.0; + double ySum = 1.0; for (loop_x=0;loop_xnX;loop_x++) { + xSum = xSum; for (loop_y=0;loop_ynY;loop_y++) { - polySum+= x + myPoly->coeff[loop_x][loop_y]; - } + polySum+= ySum * myPoly->coeff[loop_x][loop_y]; + ySum*=y; + } + xSum*=x; } @@ -495,8 +549,8 @@ double -psevalDPolynomial3D(double x, +psDEvalPolynomial3D(double x, double y, double z, - const psPolynomial3D *myPoly) + const psDPolynomial3D *myPoly) { int loop_x = 0; @@ -504,11 +558,19 @@ int loop_z = 0; double polySum = 0.0; + double xSum = 1.0; + double ySum = 1.0; + double zSum = 1.0; for (loop_x=0;loop_xnX;loop_x++) { + ySum = xSum; for (loop_y=0;loop_ynY;loop_y++) { + zSum = ySum; for (loop_z=0;loop_znZ;loop_z++) { - polySum+= x + myPoly->coeff[loop_x][loop_y][loop_z]; - } - } + polySum+= zSum * myPoly->coeff[loop_x][loop_y][loop_z]; + zSum*=z; + } + ySum*=y; + } + xSum*=x; } @@ -517,9 +579,9 @@ double -psevalDPolynomial4D(double w, +psDEvalPolynomial4D(double w, double x, double y, double z, - const psPolynomial4D *myPoly) + const psDPolynomial4D *myPoly) { int loop_w = 0; @@ -528,15 +590,27 @@ int loop_z = 0; double polySum = 0.0; + double wSum = 1.0; + double xSum = 1.0; + double ySum = 1.0; + double zSum = 1.0; for (loop_w=0;loop_wnW;loop_w++) { + xSum = wSum; for (loop_x=0;loop_xnX;loop_x++) { + ySum = xSum; for (loop_y=0;loop_ynY;loop_y++) { + zSum = ySum; for (loop_z=0;loop_znZ;loop_z++) { - polySum+= x + myPoly->coeff[loop_w][loop_x][loop_y][loop_z]; + polySum+= zSum * + myPoly->coeff[loop_w][loop_x][loop_y][loop_z]; + zSum*=z; } - } - } - } - - return(polySum); -} + ySum*=y; + } + xSum*=x; + } + wSum*=w; + } + + return(polySum); +} Index: /trunk/psLib/src/math/psSpline.h =================================================================== --- /trunk/psLib/src/math/psSpline.h (revision 749) +++ /trunk/psLib/src/math/psSpline.h (revision 750) @@ -91,4 +91,5 @@ void psPolynomial4DFree(psPolynomial4D *myPoly ///< Polynomial to destroy ); + /** Evaluate 1D polynomial */