Index: trunk/psLib/src/math/psMinimizePolyFit.c =================================================================== --- trunk/psLib/src/math/psMinimizePolyFit.c (revision 6186) +++ trunk/psLib/src/math/psMinimizePolyFit.c (revision 6193) @@ -10,8 +10,10 @@ * @author EAM, IfA * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-01-23 22:25:31 $ + * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-01-26 00:31:19 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + * + * XXX: psMatrixLUSolve() does not return error codes when the results are NANs. * * XXX: For clip-fit functions, what should we do if the mask is NULL? @@ -279,4 +281,164 @@ *****************************************************************************/ +static psPolynomial1D* VectorFitPolynomial1DOrd( + psPolynomial1D* myPoly, + const psVector *mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, + const psVector *x); + +/****************************************************************************** +vectorFitPolynomial1DCheb(): This routine will fit a Chebyshev +polynomial of degree myPoly->nX to the data points (x, y) and return the +coefficients of that polynomial. + + NOTE: We currently have implemented three algorithms. This one is + non-standard. It ignores the orthogonal properties of the Chebyshev + polys, and their known zero values. Instead, we first fit a regular + ordinary polynomial to the data. This produces the coefficients for the + various x^i terms. We then "reverse-engineer" those coefficients to + determine the coefficients of the Chebyshev polys: basically for each + c_ix^i term in the ordinary polynomial, we sum all x^i terms in the + Chebshev polys: these must be equal. This creates a set of nOrder+1 + linear equations which can be easily solved. The resulting vector is the + coefficients of the Chebyshev polys. + + This method is significantly slower than the standard NR algorithm. It + was explicitly requested that we not use the NR algorithm. + + Matrix A[[][]: + Element A[i][j] is the coefficient of the x^i term in the j-th cheby poly. + + XXX: This can be improved significantly, performance-wise. The second set + of linear equations which must be "solved" are already in upper-triangular + form. One must only perform the reverse-substitution, LUD decomposition. + + XXX: Also, we don't really need to generate the chebPolys data structure. + We simply need the matrix which corresponds to a transpose of each Cheby + polys coefficients. +*****************************************************************************/ +static psPolynomial1D *vectorFitPolynomial1DCheb( + psPolynomial1D* myPoly, + const psVector *mask, + psMaskType maskValue, + const psVector* y, + const psVector* yErr, + const psVector* x) +{ + PS_ASSERT_POLY_NON_NULL(myPoly, NULL); + PS_ASSERT_INT_LARGER_THAN_OR_EQUAL(myPoly->nX, 0, NULL); + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL); + if (yErr != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(y, yErr, NULL); + PS_ASSERT_VECTOR_TYPE(yErr, PS_TYPE_F64, NULL); + } + if (x != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(y, x, NULL); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); + } + if (mask != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(y, mask, NULL); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); + } + + psS32 polyOrder = myPoly->nX; + psS32 numPoly = polyOrder + 1; + + // + // We first fit an ordinary polynomial to the data. + // + psPolynomial1D *ordPoly = psPolynomial1DAlloc(polyOrder, PS_POLYNOMIAL_ORD); + psPolynomial1D *rc = VectorFitPolynomial1DOrd(ordPoly, mask, maskValue, y, yErr, x); + if (rc == NULL) { + psError(PS_ERR_UNKNOWN, false, "Could not fit a preliminary polynomial to the data vector. Returning NULL.\n"); + psFree(myPoly); + return(NULL); + } + + // + // Create the A-matrix and B-vector which correspond to the linear equations + // which will be solved and will then yield the Cheby poly coefficients. + // + psPolynomial1D **chebPolys = p_psCreateChebyshevPolys(numPoly); + psImage *A = psImageAlloc(numPoly, numPoly, PS_TYPE_F64); + psVector *B = psVectorAlloc(numPoly, PS_TYPE_F64); + for (psS32 i = 0 ; i < numPoly ; i++) { + for (psS32 j = 0 ; j < numPoly ; j++) { + A->data.F64[i][j] = 0.0; + } + B->data.F64[i] = ordPoly->coeff[i]; + } + + for (psS32 i = 0 ; i < numPoly ; i++) { + for (psS32 j = 0 ; j < numPoly ; j++) { + if (i <= chebPolys[j]->nX) + A->data.F64[i][j]+= chebPolys[j]->coeff[i]; + } + } + // The following statement is essential. It derives from (5.8.8) NR. + A->data.F64[0][0] = 0.5; + psFree(ordPoly); + if (psTraceGetLevel(__func__) >= 6) { + PS_IMAGE_PRINT_F64(A); + PS_VECTOR_PRINT_F64(B); + } + + if (0) { + // GaussJordan version + if (false == psGaussJordan(A, B)) { + psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations. Returning NULL.\n"); + psFree(myPoly); + myPoly = NULL; + } else { + // the first nTerm entries in B correspond directly to the desired + // polynomial coefficients. this is only true for the 1D case + for (psS32 k = 0; k < numPoly; k++) { + myPoly->coeff[k] = B->data.F64[k]; + } + } + } else { + // LUD version of the fit + psImage *ALUD = NULL; + psVector* outPerm = NULL; + psVector* coeffs = NULL; + + ALUD = psImageAlloc(numPoly, numPoly, PS_TYPE_F64); + ALUD = psMatrixLUD(ALUD, &outPerm, A); + if (ALUD == NULL) { + psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix. Returning NULL.\n"); + psFree(myPoly); + myPoly = NULL; + } else { + coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm); + if (coeffs == NULL) { + psError(PS_ERR_UNKNOWN, false, "Could not solve LUD matrix. Returning NULL.\n"); + psFree(myPoly); + myPoly = NULL; + } else { + for (psS32 k = 0; k < numPoly; k++) { + myPoly->coeff[k] = coeffs->data.F64[k]; + } + } + } + + + psFree(ALUD); + psFree(coeffs); + psFree(outPerm); + } + + psFree(A); + psFree(B); + for (psS32 i=0;in; } - // psTrace printf("vectorFitPolynomial1DChebySlow(): NUM_DATA is %d\n", NUM_DATA); psPolynomial1D **chebPolys = p_psCreateChebyshevPolys(NUM_POLY); @@ -408,10 +568,21 @@ psVector* coeffs = NULL; - // XXX: Check return codes. ALUD = psImageAlloc(NUM_POLY, NUM_POLY, PS_TYPE_F64); ALUD = psMatrixLUD(ALUD, &outPerm, A); - coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm); - for (psS32 k = 0; k < NUM_POLY; k++) { - myPoly->coeff[k] = coeffs->data.F64[k]; + if (ALUD == NULL) { + psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix. Returning NULL.\n"); + psFree(myPoly); + myPoly = NULL; + } else { + coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm); + if (coeffs == NULL) { + psError(PS_ERR_UNKNOWN, false, "Could not solve LUD matrix. Returning NULL.\n"); + psFree(myPoly); + myPoly = NULL; + } else { + for (psS32 k = 0; k < NUM_POLY; k++) { + myPoly->coeff[k] = coeffs->data.F64[k]; + } + } } @@ -437,7 +608,8 @@ polynomial. - NOTE: We currently have two algorithms. This is standard method which + NOTE: We currently have three algorithms. This is standard method which uses the orthogonal properties of the Chebyshev polys, and their known - zero values. This is significantly faster than the chi-squared approach. + zero values. This is significantly faster than the chi-squared + approaches. NOTE: This function will not work properly if the x-vector does not fully span @@ -621,5 +793,4 @@ // Build the B and A data structs. // XXX EAM : use temp pointers eg vB = B->data.F64 to save redirects - // XXX EAM : this function is only valid for data vectors of F64 for (psS32 k = 0; k < f->n; k++) { if ((mask != NULL) && (mask->data.U8[k] && maskValue)) { @@ -671,11 +842,23 @@ psVector* coeffs = NULL; - // XXX: Check return codes. ALUD = psImageAlloc(nTerm, nTerm, PS_TYPE_F64); ALUD = psMatrixLUD(ALUD, &outPerm, A); - coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm); - for (psS32 k = 0; k < nTerm; k++) { - myPoly->coeff[k] = coeffs->data.F64[k]; - } + if (ALUD == NULL) { + psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix. Returning NULL.\n"); + psFree(myPoly); + myPoly = NULL; + } else { + coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm); + if (coeffs == NULL) { + psError(PS_ERR_UNKNOWN, false, "Could not solve LUD matrix. Returning NULL.\n"); + psFree(myPoly); + myPoly = NULL; + } else { + for (psS32 k = 0; k < nTerm; k++) { + myPoly->coeff[k] = coeffs->data.F64[k]; + } + } + } + psFree(ALUD); psFree(coeffs); @@ -771,7 +954,10 @@ if (1) { - poly = vectorFitPolynomial1DChebySlow(poly, mask, maskValue, f64, fErr64, x64); + poly = vectorFitPolynomial1DCheb(poly, mask, maskValue, f64, fErr64, x64); } else { - poly = vectorFitPolynomial1DChebyFast(poly, mask, maskValue, f64, fErr64, x64); + if (0) { + poly = vectorFitPolynomial1DChebySlow(poly, mask, maskValue, f64, fErr64, x64); + poly = vectorFitPolynomial1DChebyFast(poly, mask, maskValue, f64, fErr64, x64); + } } if (x == NULL) { @@ -989,5 +1175,4 @@ PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL); PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, NULL); - PS_ASSERT_VECTOR_NON_NULL(f, NULL); PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, NULL); @@ -1560,16 +1745,26 @@ psVector* coeffs = NULL; - // XXX: Check return codes. ALUD = psImageAlloc(nTerm, nTerm, PS_TYPE_F64); ALUD = psMatrixLUD(ALUD, &outPerm, A); - coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm); - - // select the appropriate solution entries - for (psS32 ix = 0; ix < nXterm; ix++) { - for (psS32 iy = 0; iy < nYterm; iy++) { - for (psS32 iz = 0; iz < nZterm; iz++) { - psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm; - myPoly->coeff[ix][iy][iz] = coeffs->data.F64[nx]; - myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[nx][nx]); + if (ALUD == NULL) { + psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix. Returning NULL.\n"); + psFree(myPoly); + myPoly = NULL; + } else { + coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm); + if (coeffs == NULL) { + psError(PS_ERR_UNKNOWN, false, "Could not solve LUD matrix. Returning NULL.\n"); + psFree(myPoly); + myPoly = NULL; + } else { + // select the appropriate solution entries + for (psS32 ix = 0; ix < nXterm; ix++) { + for (psS32 iy = 0; iy < nYterm; iy++) { + for (psS32 iz = 0; iz < nZterm; iz++) { + psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm; + myPoly->coeff[ix][iy][iz] = coeffs->data.F64[nx]; + myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[nx][nx]); + } + } } } @@ -2108,17 +2303,27 @@ psVector* coeffs = NULL; - // XXX: Check return codes. ALUD = psImageAlloc(nTerm, nTerm, PS_TYPE_F64); ALUD = psMatrixLUD(ALUD, &outPerm, A); - coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm); - - // select the appropriate solution entries - for (psS32 ix = 0; ix < nXterm; ix++) { - for (psS32 iy = 0; iy < nYterm; iy++) { - for (psS32 iz = 0; iz < nZterm; iz++) { - for (psS32 it = 0; it < nTterm; it++) { - psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm; - myPoly->coeff[ix][iy][iz][it] = coeffs->data.F64[nx]; - myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[nx][nx]); + if (ALUD == NULL) { + psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix. Returning NULL.\n"); + psFree(myPoly); + myPoly = NULL; + } else { + coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm); + if (coeffs == NULL) { + psError(PS_ERR_UNKNOWN, false, "Could not solve LUD matrix. Returning NULL.\n"); + psFree(myPoly); + myPoly = NULL; + } else { + // select the appropriate solution entries + for (psS32 ix = 0; ix < nXterm; ix++) { + for (psS32 iy = 0; iy < nYterm; iy++) { + for (psS32 iz = 0; iz < nZterm; iz++) { + for (psS32 it = 0; it < nTterm; it++) { + psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm; + myPoly->coeff[ix][iy][iz][it] = coeffs->data.F64[nx]; + myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[nx][nx]); + } + } } }