Index: trunk/psLib/src/math/psMinimizePolyFit.c =================================================================== --- trunk/psLib/src/math/psMinimizePolyFit.c (revision 10778) +++ trunk/psLib/src/math/psMinimizePolyFit.c (revision 10848) @@ -10,6 +10,6 @@ * @author EAM, IfA * - * @version $Revision: 1.26 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-12-17 09:43:48 $ + * @version $Revision: 1.27 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-12-29 04:38:42 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -55,4 +55,8 @@ }\ }\ + +// free a local temporary F64 vector (TEMP) which is a copy of a non-F64 vector (ORIG) +# define PS_FREE_TEMP_F64_VECTOR(ORIG, TEMP) \ +if ((ORIG != NULL) && (ORIG->type.type != PS_TYPE_F64)) { psFree(TEMP); } /*****************************************************************************/ @@ -281,5 +285,4 @@ *****************************************************************************/ - /****************************************************************************** ****************************************************************************** @@ -288,12 +291,4 @@ *****************************************************************************/ -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 @@ -301,5 +296,5 @@ coefficients of that polynomial. *****************************************************************************/ -static psPolynomial1D *vectorFitPolynomial1DCheb( +static bool vectorFitPolynomial1DCheb( psPolynomial1D* myPoly, const psVector *mask, @@ -404,6 +399,5 @@ if (!psMatrixGJSolve(A, B)) { psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations. Returning NULL.\n"); - psFree(myPoly); - myPoly = NULL; + goto escape_GJ; } else { // the first nTerm entries in B correspond directly to the desired @@ -427,12 +421,10 @@ if (ALUD == NULL) { psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix. Returning NULL.\n"); - psFree(myPoly); - myPoly = NULL; + goto escape_LUD; } 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; + goto escape_LUD; } else { for (psS32 k = 0; k < numTerms; k++) { @@ -441,15 +433,22 @@ } } - - psFree(ALUD); psFree(coeffs); psFree(outPerm); } - psFree(A); psFree(B); - - return(myPoly); + return true; + +escape_LUD: + // XXX drop the LUD version! + psFree(A); + psFree(B); + return false; + +escape_GJ: + psFree(A); + psFree(B); + return false; } @@ -459,6 +458,9 @@ x and fErr vectors may be NULL. All non-NULL vectors must be of type PS_TYPE_F64. - *****************************************************************************/ -static psPolynomial1D* VectorFitPolynomial1DOrd( + +XXX EAM : since this is a private function, can we drop the ASSERTS? +XXX EAM : can we drop the LUD version? it does not calculate coeffErr values!! +*****************************************************************************/ +static bool VectorFitPolynomial1DOrd( psPolynomial1D* myPoly, const psVector *mask, @@ -469,18 +471,18 @@ { psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); - PS_ASSERT_POLY_NON_NULL(myPoly, NULL); - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, NULL); + PS_ASSERT_POLY_NON_NULL(myPoly, false); + PS_ASSERT_VECTOR_NON_NULL(f, false); + PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, false); if (mask) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, false); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false); } if (x) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, false); } if (fErr) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); - PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, false); + PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, false); } @@ -510,9 +512,8 @@ if (!psImageInit(A, 0.0) || !psVectorInit(B, 0.0)) { psError(PS_ERR_UNKNOWN, false, "Could initialize data structures A, B. Returning NULL.\n"); - psFree(myPoly); psFree(A); psFree(B); psTrace("psLib.math", 4, "---- %s() End ----\n", __func__); - return(NULL); + return false; } @@ -596,6 +597,5 @@ if (!psMatrixGJSolve(A, B)) { psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations. Returning NULL.\n"); - psFree(myPoly); - myPoly = NULL; + goto escape_GJ; } else { // the first nTerm entries in B correspond directly to the desired @@ -617,12 +617,10 @@ if (ALUD == NULL) { psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix. Returning NULL.\n"); - psFree(myPoly); - myPoly = NULL; + goto escape_LUD; } 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; + goto escape_LUD; } else { for (psS32 k = 0; k < nTerm; k++) { @@ -632,16 +630,23 @@ } } - psFree(ALUD); psFree(coeffs); psFree(outPerm); } - - psFree(A); psFree(B); psTrace("psLib.math", 4, "---- %s() End ----\n", __func__); - return (myPoly); + return true; + +escape_LUD: + psFree(A); + psFree(B); + return false; + +escape_GJ: + psFree(A); + psFree(B); + return false; } @@ -653,5 +658,5 @@ conversion only. *****************************************************************************/ -psPolynomial1D *psVectorFitPolynomial1D( +bool psVectorFitPolynomial1D( psPolynomial1D *poly, const psVector *mask, @@ -661,41 +666,40 @@ const psVector *x) { - // Internal pointers for possibly NULL or mis-typed vectors. + PS_ASSERT_POLY_NON_NULL(poly, false); + PS_ASSERT_INT_NONNEGATIVE(poly->nX, false); + + PS_ASSERT_VECTOR_NON_NULL(f, false); + PS_ASSERT_VECTOR_NON_EMPTY(f, false); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false); + if (mask != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, false); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false); + } + if (fErr != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, false); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, false); + } + if (x != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, false); + } + + // Convert input vectors to F64 if necessary. + psVector *f64 = (f->type.type == PS_TYPE_F64) ? (psVector *) f : psVectorCopy (NULL, f, PS_TYPE_F64); psVector *x64 = NULL; - psVector *f64 = NULL; + if (x != NULL) { + x64 = (x->type.type == PS_TYPE_F64) ? (psVector *) x : psVectorCopy (NULL, x, PS_TYPE_F64); + } psVector *fErr64 = NULL; - - PS_ASSERT_POLY_NON_NULL(poly, NULL); - //PS_ASSERT_INT_NONNEGATIVE(poly->nX, NULL); - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_NON_EMPTY(f, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL); - if (mask != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); - } - if (x != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL); - } if (fErr != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL); - } - - f64 = (f->type.type == PS_TYPE_F64) ? (psVector *) f : psVectorCopy(NULL, f, PS_TYPE_F64); - - if (x != NULL) { - x64 = (x->type.type == PS_TYPE_F64) ? (psVector *) x : psVectorCopy(NULL, x, PS_TYPE_F64); - } - - if (fErr != NULL) { - fErr64 = (fErr->type.type == PS_TYPE_F64) ? (psVector *) fErr : psVectorCopy(NULL, fErr, PS_TYPE_F64); - } + fErr64 = (fErr->type.type == PS_TYPE_F64) ? (psVector *) fErr : psVectorCopy (NULL, fErr, PS_TYPE_F64); + } + + bool result = true; switch (poly->type) { case PS_POLYNOMIAL_ORD: - poly = VectorFitPolynomial1DOrd(poly, mask, maskValue, f64, fErr64, x64); - if (poly == NULL) { + result = VectorFitPolynomial1DOrd(poly, mask, maskValue, f64, fErr64, x64); + if (!result) { psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial. Returning NULL.\n"); } @@ -703,5 +707,5 @@ case PS_POLYNOMIAL_CHEB: if (mask != NULL) { - // psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring mask and maskValue with Chebyshev polynomials.\n"); + psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring mask and maskValue with Chebyshev polynomials.\n"); } if (fErr != NULL) { @@ -713,5 +717,8 @@ } - poly = vectorFitPolynomial1DCheb(poly, mask, maskValue, f64, fErr64, x64); + result = vectorFitPolynomial1DCheb(poly, mask, maskValue, f64, fErr64, x64); + if (!result) { + psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial. Returning NULL.\n"); + } if (x == NULL) { @@ -721,26 +728,18 @@ default: psError(PS_ERR_UNKNOWN, true, "Incorrect polynomial type (%d). Returning NULL.\n", poly->type); - poly = NULL; + result = false; break; } // Free psVectors that were created for NULL arguments. - if (f->type.type != PS_TYPE_F64) { - psFree(f64); - } - - if ((x != NULL) && (x->type.type != PS_TYPE_F64)) { - psFree(x64); - } - - if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { - psFree(fErr64); - } - - return(poly); + PS_FREE_TEMP_F64_VECTOR (f, f64); + PS_FREE_TEMP_F64_VECTOR (x, x64); + PS_FREE_TEMP_F64_VECTOR (fErr, fErr64); + + return result; } // This function accepts F32 and F64 input vectors. -psPolynomial1D *psVectorClipFitPolynomial1D( +bool psVectorClipFitPolynomial1D( psPolynomial1D *poly, psStats *stats, @@ -752,14 +751,15 @@ { psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_PTR_NON_NULL(stats, NULL); - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL); - PS_ASSERT_VECTOR_NON_NULL(mask, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(mask, f, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); + PS_ASSERT_POLY_NON_NULL(poly, false); + PS_ASSERT_PTR_NON_NULL(stats, false); + PS_ASSERT_VECTOR_NON_NULL(f, false); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false); + PS_ASSERT_VECTOR_NON_NULL(mask, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(mask, f, false); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false); + if (fErr != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, f, NULL); - PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, f, false); + PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, false); } @@ -767,6 +767,6 @@ psVector *x = NULL; if (xIn != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(xIn, f, NULL); - PS_ASSERT_VECTOR_TYPE(xIn, f->type.type, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(xIn, f, false); + PS_ASSERT_VECTOR_TYPE(xIn, f->type.type, false); x = (psVector *) xIn; } else { @@ -781,6 +781,21 @@ } else { psError(PS_ERR_UNKNOWN, true, "Error, bad poly type.\n"); - return(NULL); - } + return false; + } + } + + // the user supplies one of various stats option pairs, + // determine the desired mean and stdev STATS options: + // XXX enforce consistency? + // XXX psStatsGetValue() probably has inverted precedence + psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN_V2); + psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV_V2); + if (!meanOption) { + psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected"); + return false; + } + if (!stdevOption) { + psError(PS_ERR_UNKNOWN, true, "no valid stdev stats option selected"); + return false; } @@ -798,14 +813,6 @@ minClipSigma = fabs(stats->clipSigma); } - - psVector *fit = NULL; psVector *resid = psVectorAlloc(f->n, PS_TYPE_F64); - // eventual expansion: user supplies one of various stats option pairs, - // eg (SAMPLE_MEAN | SAMPLE_STDEV) and the correct pair is used to - // evaluate the clipping sigma - // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used - stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - stats->options |= (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); psTrace("psLib.math", 4, "stats->clipIter is %d\n", stats->clipIter); psTrace("psLib.math", 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma); @@ -822,18 +829,17 @@ } } - poly = psVectorFitPolynomial1D(poly, mask, maskValue, f, fErr, x); - if (poly == NULL) { - psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial. Returning NULL.\n"); + if (!psVectorFitPolynomial1D(poly, mask, maskValue, f, fErr, x)) { + psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial. Returning false.\n"); if (xIn == NULL) { psFree(x); } - return(NULL); - } - - fit = psPolynomial1DEvalVector(poly, x); + return false; + } + + psVector *fit = psPolynomial1DEvalVector(poly, x); if (fit == NULL) { - psError(PS_ERR_UNKNOWN, false, "Could not call psPolynomial3DEvalVector(). Returning NULL.\n"); + psError(PS_ERR_UNKNOWN, false, "Could not call psPolynomial3DEvalVector(). Returning false.\n"); psFree(resid); - return(NULL); + return false; } for (psS32 i = 0 ; i < f->n ; i++) { @@ -856,23 +862,17 @@ if (!psVectorStats(stats, resid, NULL, mask, maskValue)) { - psError(PS_ERR_UNKNOWN, false, "Could not compute statistics on the resid vector. Returning NULL.\n"); - psFree(resid) - psFree(fit) - return(NULL); - } - # if (USE_ROBUST_STATS_FOR_CLIPPING) - psTrace("psLib.math", 5, "Median is %f\n", stats->robustMedian); - psTrace("psLib.math", 5, "Stdev is %f\n", stats->robustStdev); - psF32 minClipValue = -minClipSigma*stats->robustStdev; - psF32 maxClipValue = +maxClipSigma*stats->robustStdev; - psF32 clipMedian = stats->robustMedian; - # else - - psTrace("psLib.math", 5, "Median is %f\n", stats->sampleMedian); - psTrace("psLib.math", 5, "Stdev is %f\n", stats->sampleStdev); - psF32 minClipValue = -minClipSigma*stats->robustStdev; - psF32 maxClipValue = +maxClipSigma*stats->robustStdev; - psF32 clipMedian = stats->sampleMedian; - # endif + psError(PS_ERR_UNKNOWN, false, "Could not compute statistics on the resid vector. Returning false.\n"); + psFree(resid); + psFree(fit); + return false; + } + + double meanValue = psStatsGetValue (stats, meanOption); + double stdevValue = psStatsGetValue (stats, stdevOption); + + psTrace("psLib.math", 5, "Mean is %f\n", meanValue); + psTrace("psLib.math", 5, "Stdev is %f\n", stdevValue); + psF32 minClipValue = -minClipSigma*stdevValue; + psF32 maxClipValue = +maxClipSigma*stdevValue; // set mask if pts are not valid @@ -884,5 +884,5 @@ } - if ((resid->data.F64[i] - clipMedian > maxClipValue) || (resid->data.F64[i] - clipMedian < minClipValue)) { + if ((resid->data.F64[i] - meanValue > maxClipValue) || (resid->data.F64[i] - meanValue < minClipValue)) { if (f->type.type == PS_TYPE_F64) { psTrace("psLib.math", 6, "Masking element %d (%f). resid->data.F64[%d] is %f\n", @@ -906,4 +906,5 @@ // psTrace("psLib.math", 6, "keeping %d of %ld pts for fit\n", Nkeep, x->n); + stats->clippedNvalues = Nkeep; psFree(fit); } @@ -917,5 +918,5 @@ psTrace("psLib.math", 3, "---- %s() end ----\n", __func__); - return (poly); + return true; } @@ -933,5 +934,5 @@ *****************************************************************************/ -static psPolynomial2D* VectorFitPolynomial2DOrd( +static bool VectorFitPolynomial2DOrd( psPolynomial2D* myPoly, const psVector* mask, @@ -943,22 +944,22 @@ { psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); - PS_ASSERT_POLY_NON_NULL(myPoly, NULL); - 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); + PS_ASSERT_POLY_NON_NULL(myPoly, false); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, false); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, false); + PS_ASSERT_VECTOR_NON_NULL(f, false); + PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, false); if (fErr != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(y, fErr, NULL); - PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, NULL); - } - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(y, fErr, false); + PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, false); + } + PS_ASSERT_VECTOR_NON_NULL(x, false); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false); + PS_ASSERT_VECTOR_NON_NULL(y, false); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false); if (mask != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(y, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(y, mask, false); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false); } @@ -974,9 +975,8 @@ if (!psImageInit(A, 0.0) || !psVectorInit(B, 0.0)) { psError(PS_ERR_UNKNOWN, false, "Could initialize data structures A, B. Returning NULL.\n"); - psFree(myPoly); psFree(A); psFree(B); psTrace("psLib.math", 4, "---- %s() End ----\n", __func__); - return(NULL); + return false; } @@ -1040,21 +1040,21 @@ if (!psMatrixGJSolve(A, B)) { psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations. Returning NULL.\n"); - psFree(myPoly); - myPoly = NULL; - } else { - // select the appropriate solution entries - for (int i = 0; i < nTerm; i++) { - int l = i / nYterm; // x index - int m = i % nYterm; // y index - myPoly->coeff[l][m] = B->data.F64[i]; - myPoly->coeffErr[l][m] = sqrt(A->data.F64[i][i]); - } - } - + psFree(A); + psFree(B); + return false; + } + + // select the appropriate solution entries + for (int i = 0; i < nTerm; i++) { + int l = i / nYterm; // x index + int m = i % nYterm; // y index + myPoly->coeff[l][m] = B->data.F64[i]; + myPoly->coeffErr[l][m] = sqrt(A->data.F64[i][i]); + } psFree(A); psFree(B); psTrace("psLib.math", 4, "---- %s() end ----\n", __func__); - return (myPoly); + return true; } @@ -1065,5 +1065,5 @@ vector conversion only. *****************************************************************************/ -psPolynomial2D *psVectorFitPolynomial2D( +bool psVectorFitPolynomial2D( psPolynomial2D *poly, const psVector *mask, @@ -1074,132 +1074,64 @@ const psVector *y) { - // Internal pointers for possibly NULL or mis-typed vectors. - psVector *x64 = NULL; - psVector *y64 = NULL; - psVector *f64 = NULL; + PS_ASSERT_POLY_NON_NULL(poly, false); + PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, false); + + PS_ASSERT_VECTOR_NON_NULL(f, false); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false); + PS_ASSERT_VECTOR_NON_NULL(x, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false); + PS_ASSERT_VECTOR_NON_NULL(y, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false); + if (mask != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, false); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false); + } + if (fErr != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, false); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, false); + } + + // Convert input vectors to F64 if necessary. + psVector *f64 = (f->type.type == PS_TYPE_F64) ? (psVector *) f : psVectorCopy(NULL, f, PS_TYPE_F64); + psVector *x64 = (x->type.type == PS_TYPE_F64) ? (psVector *) x : psVectorCopy(NULL, x, PS_TYPE_F64); + psVector *y64 = (y->type.type == PS_TYPE_F64) ? (psVector *) y : psVectorCopy(NULL, y, PS_TYPE_F64); + psVector *fErr64 = NULL; - - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL); - - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL); - - if (mask != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); - } - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); if (fErr != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL); - } - - // - // Convert input vectors to F64 if necessary. - // - if (f->type.type != PS_TYPE_F64) { - f64 = psVectorCopy(NULL, f, PS_TYPE_F64); - } else { - f64 = (psVector *) f; - } - - if (x->type.type != PS_TYPE_F64) { - x64 = psVectorCopy(NULL, x, PS_TYPE_F64); - } else { - x64 = (psVector *) x; - } - - if (y->type.type != PS_TYPE_F64) { - y64 = psVectorCopy(NULL, y, PS_TYPE_F64); - } else { - y64 = (psVector *) y; - } - - // - // fErr - // - if (fErr != NULL) { - if (fErr->type.type != PS_TYPE_F64) { - fErr64 = psVectorCopy(NULL, fErr, PS_TYPE_F64); - } else { - fErr64 = (psVector *) fErr; - } - } - - if (poly->type == PS_POLYNOMIAL_ORD) { - poly = VectorFitPolynomial2DOrd(poly, mask, maskValue, f64, fErr64, x64, y64); - if (poly == NULL) { + fErr64 = (fErr->type.type == PS_TYPE_F64) ? (psVector *) fErr : psVectorCopy(NULL, fErr, PS_TYPE_F64); + } + + bool result = true; + + switch (poly->type) { + case PS_POLYNOMIAL_ORD: + result = VectorFitPolynomial2DOrd(poly, mask, maskValue, f64, fErr64, x64, y64); + if (!result) { psError(PS_ERR_UNKNOWN, true, "Could not fit polynomial. Returning NULL.\n"); - // Free psVectors that were created for NULL arguments. - if (f->type.type != PS_TYPE_F64) { - psFree(f64); - } - - if (x->type.type != PS_TYPE_F64) { - psFree(x64); - } - - if (y->type.type != PS_TYPE_F64) { - psFree(y64); - } - - if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { - psFree(fErr64); - } - return(NULL); - } - } else if (poly->type == PS_POLYNOMIAL_CHEB) { + } + break; + case PS_POLYNOMIAL_CHEB: if (mask != NULL) { psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring mask and maskValue with Chebyshev polynomials.\n"); } - psLogMsg(__func__, PS_LOG_WARN, "WARNING: 2-D Chebyshev polynomial vector fitting has not been implemented. Returning NULL.\n"); - psFree(poly); - poly = NULL; - } else { - // Free psVectors that were created for NULL arguments. - if (f->type.type != PS_TYPE_F64) { - psFree(f64); - } - - if (x->type.type != PS_TYPE_F64) { - psFree(x64); - } - - if (y->type.type != PS_TYPE_F64) { - psFree(y64); - } - - if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { - psFree(fErr64); - } + psError(PS_ERR_UNKNOWN, true, "2-D Chebyshev polynomial vector fitting has not been implemented. Returning NULL.\n"); + result = false; + break; + default: psError(PS_ERR_UNKNOWN, true, "Incorrect polynomial type. Returning NULL.\n"); - } - + result = false; + break; + } // Free psVectors that were created for NULL arguments. - if (f->type.type != PS_TYPE_F64) { - psFree(f64); - } - - if (x->type.type != PS_TYPE_F64) { - psFree(x64); - } - - if (y->type.type != PS_TYPE_F64) { - psFree(y64); - } - - if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { - psFree(fErr64); - } - - return(poly); + PS_FREE_TEMP_F64_VECTOR (f, f64); + PS_FREE_TEMP_F64_VECTOR (x, x64); + PS_FREE_TEMP_F64_VECTOR (y, y64); + PS_FREE_TEMP_F64_VECTOR (fErr, fErr64); + + return result; } -psPolynomial2D *psVectorClipFitPolynomial2D( +bool psVectorClipFitPolynomial2D( psPolynomial2D *poly, psStats *stats, @@ -1212,25 +1144,40 @@ { psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL); - PS_ASSERT_PTR_NON_NULL(stats, NULL); - PS_ASSERT_VECTOR_NON_NULL(mask, NULL); - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL); - if (mask != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); - } - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); - PS_ASSERT_VECTOR_TYPE(x, f->type.type, NULL); - - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); - PS_ASSERT_VECTOR_TYPE(y, f->type.type, NULL); + PS_ASSERT_POLY_NON_NULL(poly, false); + PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, false); + PS_ASSERT_PTR_NON_NULL(stats, false); + PS_ASSERT_VECTOR_NON_NULL(mask, false); + PS_ASSERT_VECTOR_NON_NULL(f, false); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false); + + PS_ASSERT_VECTOR_NON_NULL(x, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false); + PS_ASSERT_VECTOR_TYPE(x, f->type.type, false); + + PS_ASSERT_VECTOR_NON_NULL(y, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false); + PS_ASSERT_VECTOR_TYPE(y, f->type.type, false); + + PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, false); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false); if (fErr != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); - PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, false); + PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, false); + } + + // the user supplies one of various stats option pairs, + // determine the desired mean and stdev STATS options: + // XXX enforce consistency? + // XXX psStatsGetValue() probably has inverted precedence + psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN_V2); + psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV_V2); + if (!meanOption) { + psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected"); + return false; + } + if (!stdevOption) { + psError(PS_ERR_UNKNOWN, true, "no valid stdev stats option selected"); + return false; } @@ -1250,10 +1197,4 @@ psVector *resid = psVectorAlloc(f->n, PS_TYPE_F64); - // eventual expansion: user supplies one of various stats option pairs, - // eg (SAMPLE_MEAN | SAMPLE_STDEV) and the correct pair is used to - // evaluate the clipping sigma - // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used - stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - stats->options |= (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); psTrace("psLib.math", 4, "stats->clipIter is %d\n", stats->clipIter); psTrace("psLib.math", 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma); @@ -1270,9 +1211,8 @@ } - poly = psVectorFitPolynomial2D(poly, mask, maskValue, f, fErr, x, y); - if (poly == NULL) { - psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to the data. Returning NULL.\n"); + if (!psVectorFitPolynomial2D(poly, mask, maskValue, f, fErr, x, y)) { + psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to the data. Returning false.\n"); psFree(resid) - return(NULL); + return false; } @@ -1281,5 +1221,5 @@ psError(PS_ERR_UNKNOWN, false, "Could not call psPolynomial3DEvalVector(). Returning NULL.\n"); psFree(resid) - return(NULL); + return false; } @@ -1307,24 +1247,14 @@ psFree(resid) psFree(fit) - return(NULL); - } - # if (USE_ROBUST_STATS_FOR_CLIPPING) - psTrace("psLib.math", 5, "Median is %f\n", stats->robustMedian); - psTrace("psLib.math", 5, "Stdev is %f\n", stats->robustStdev); - psTrace("psLib.math", 5, "Sample Median is %f\n", stats->sampleMedian); - psTrace("psLib.math", 5, "Sample Stdev is %f\n", stats->sampleStdev); - psF32 minClipValue = -minClipSigma*stats->robustStdev; - psF32 maxClipValue = +maxClipSigma*stats->robustStdev; - psF32 clipMedian = stats->robustMedian; - # else - - psTrace("psLib.math", 5, "Median is %f\n", stats->sampleMedian); - psTrace("psLib.math", 5, "Stdev is %f\n", stats->sampleStdev); - psTrace("psLib.math", 5, "Robust Median is %f\n", stats->robustMedian); - psTrace("psLib.math", 5, "Robust Stdev is %f\n", stats->robustStdev); - psF32 minClipValue = -minClipSigma*stats->robustStdev; - psF32 maxClipValue = +maxClipSigma*stats->robustStdev; - psF32 clipMedian = stats->sampleMedian; - # endif + return false; + } + + double meanValue = psStatsGetValue (stats, meanOption); + double stdevValue = psStatsGetValue (stats, stdevOption); + + psTrace("psLib.math", 5, "Mean is %f\n", meanValue); + psTrace("psLib.math", 5, "Stdev is %f\n", stdevValue); + psF32 minClipValue = -minClipSigma*stdevValue; + psF32 maxClipValue = +maxClipSigma*stdevValue; // set mask if pts are not valid @@ -1336,5 +1266,5 @@ } - if ((resid->data.F64[i] - clipMedian > maxClipValue) || (resid->data.F64[i] - clipMedian < minClipValue)) { + if ((resid->data.F64[i] - meanValue > maxClipValue) || (resid->data.F64[i] - meanValue < minClipValue)) { if (fit->type.type == PS_TYPE_F64) { psTrace("psLib.math", 6, "Masking element %d (%f). resid->data.F64[%d] is %f\n", @@ -1352,6 +1282,6 @@ Nkeep++; } - psTrace("psLib.math", 4, "keeping %d of %ld pts for fit\n", Nkeep, x->n); + stats->clippedNvalues = Nkeep; psFree(fit); } @@ -1359,11 +1289,6 @@ psFree(resid); - if (poly == NULL) { - psError(PS_ERR_UNKNOWN, true, "Could not fit a polynomial to the data. Returning NULL.\n"); - return(NULL); - } - psTrace("psLib.math", 3, "---- %s() end ----\n", __func__); - return(poly); + return true; } @@ -1381,5 +1306,5 @@ *****************************************************************************/ -static psPolynomial3D* VectorFitPolynomial3DOrd( +static bool VectorFitPolynomial3DOrd( psPolynomial3D* myPoly, const psVector* mask, @@ -1392,29 +1317,28 @@ { psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); - PS_ASSERT_POLY_NON_NULL(myPoly, NULL); - PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL); - PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, NULL); - PS_ASSERT_INT_NONNEGATIVE(myPoly->nZ, NULL); - - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, NULL); + PS_ASSERT_POLY_NON_NULL(myPoly, false); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, false); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, false); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nZ, false); + + PS_ASSERT_VECTOR_NON_NULL(f, false); + PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, false); if (fErr != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(y, fErr, NULL); - PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, NULL); - } - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); - PS_ASSERT_VECTOR_NON_NULL(z, NULL); - PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F64, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(y, fErr, false); + PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, false); + } + PS_ASSERT_VECTOR_NON_NULL(x, false); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false); + PS_ASSERT_VECTOR_NON_NULL(y, false); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false); + PS_ASSERT_VECTOR_NON_NULL(z, false); + PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F64, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, false); if (mask != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); - } - + PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, false); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false); + } int nXterm = 1 + myPoly->nX; // Number of x terms @@ -1429,9 +1353,8 @@ if (!psImageInit(A, 0.0) || !psVectorInit(B, 0.0)) { psError(PS_ERR_UNKNOWN, false, "Could initialize data structures A, B. Returning NULL.\n"); - psFree(myPoly); psFree(A); psFree(B); psTrace("psLib.math", 4, "---- %s() End ----\n", __func__); - return(NULL); + return false; } @@ -1511,8 +1434,6 @@ // The matrices were overflowing, so I switched to LUD. if (!psMatrixGJSolve(A, B)) { - psFree(A); - psFree(B); psError(PS_ERR_UNKNOWN, false, "Failed to perform GaussJordan elimination.\n"); - return(NULL); + goto escape_GJ; } // select the appropriate solution entries @@ -1534,12 +1455,10 @@ if (ALUD == NULL) { psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix. Returning NULL.\n"); - psFree(myPoly); - myPoly = NULL; + goto escape_LUD; } 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; + goto escape_LUD; } else { // select the appropriate solution entries @@ -1565,5 +1484,16 @@ psTrace("psLib.math", 4, "---- %s() end ----\n", __func__); - return (myPoly); + return true; + +escape_LUD: + psFree(A); + psFree(B); + return false; + +escape_GJ: + + psFree(A); + psFree(B); + return false; } @@ -1574,5 +1504,5 @@ vector conversion only. *****************************************************************************/ -psPolynomial3D *psVectorFitPolynomial3D( +bool psVectorFitPolynomial3D( psPolynomial3D *poly, const psVector *mask, @@ -1584,147 +1514,68 @@ const psVector *z) { - // Internal pointers for possibly NULL or mis-typed vectors. - psVector *x64 = NULL; - psVector *y64 = NULL; - psVector *z64 = NULL; - psVector *f64 = NULL; + PS_ASSERT_POLY_NON_NULL(poly, false); + PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, false); + + PS_ASSERT_VECTOR_NON_NULL(f, false); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false); + PS_ASSERT_VECTOR_NON_NULL(x, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false); + PS_ASSERT_VECTOR_NON_NULL(y, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false); + PS_ASSERT_VECTOR_NON_NULL(z, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, false); + if (mask != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, false); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false); + } + if (fErr != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, false); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, false); + } + + // Convert input vectors to F64 if necessary. + psVector *f64 = (f->type.type == PS_TYPE_F64) ? (psVector *) f : psVectorCopy(NULL, f, PS_TYPE_F64); + psVector *x64 = (x->type.type == PS_TYPE_F64) ? (psVector *) x : psVectorCopy(NULL, x, PS_TYPE_F64); + psVector *y64 = (y->type.type == PS_TYPE_F64) ? (psVector *) y : psVectorCopy(NULL, y, PS_TYPE_F64); + psVector *z64 = (z->type.type == PS_TYPE_F64) ? (psVector *) z : psVectorCopy(NULL, z, PS_TYPE_F64); + psVector *fErr64 = NULL; - - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL); - - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL); - if (mask != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); - } - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); - PS_ASSERT_VECTOR_NON_NULL(z, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL); if (fErr != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); - // PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, NULL); - } - - // - // Convert input vectors to F64 if necessary. - // - if (f->type.type != PS_TYPE_F64) { - f64 = psVectorCopy(NULL, f, PS_TYPE_F64); - } else { - f64 = (psVector *) f; - } - if (x->type.type != PS_TYPE_F64) { - x64 = psVectorCopy(NULL, x, PS_TYPE_F64); - } else { - x64 = (psVector *) x; - } - if (y->type.type != PS_TYPE_F64) { - y64 = psVectorCopy(NULL, y, PS_TYPE_F64); - } else { - y64 = (psVector *) y; - } - - if (z->type.type != PS_TYPE_F64) { - z64 = psVectorCopy(NULL, z, PS_TYPE_F64); - } else { - z64 = (psVector *) z; - } - - if (fErr != NULL) { - if (fErr->type.type != PS_TYPE_F64) { - fErr64 = psVectorCopy(NULL, fErr, PS_TYPE_F64); - } else { - fErr64 = (psVector *) fErr; - } - } - - if (poly->type == PS_POLYNOMIAL_ORD) { - poly = VectorFitPolynomial3DOrd(poly, mask, maskValue, f64, fErr64, x64, y64, z64); - if (poly == NULL) { + fErr64 = (fErr->type.type == PS_TYPE_F64) ? (psVector *) fErr : psVectorCopy(NULL, fErr, PS_TYPE_F64); + } + + bool result = true; + + switch (poly->type) { + case PS_POLYNOMIAL_ORD: + result = VectorFitPolynomial3DOrd(poly, mask, maskValue, f64, fErr64, x64, y64, z64); + if (!result) { psError(PS_ERR_UNKNOWN, true, "Could not fit polynomial. Returning NULL.\n"); - // Free psVectors that were created for NULL arguments. - if (f->type.type != PS_TYPE_F64) { - psFree(f64); - } - - if (x->type.type != PS_TYPE_F64) { - psFree(x64); - } - - if (y->type.type != PS_TYPE_F64) { - psFree(y64); - } - - if (z->type.type != PS_TYPE_F64) { - psFree(z64); - } - - if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { - psFree(fErr64); - } - return(NULL); - } - } else if (poly->type == PS_POLYNOMIAL_CHEB) { + } + break; + case PS_POLYNOMIAL_CHEB: if (mask != NULL) { psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring mask and maskValue with Chebyshev polynomials.\n"); } - psLogMsg(__func__, PS_LOG_WARN, "WARNING: 3-D Chebyshev polynomial vector fitting has not been implemented. Returning NULL.\n"); - psFree(poly); - poly = NULL; - } else { - // Free psVectors that were created for NULL arguments. - if (f->type.type != PS_TYPE_F64) { - psFree(f64); - } - - if (x->type.type != PS_TYPE_F64) { - psFree(x64); - } - - if (y->type.type != PS_TYPE_F64) { - psFree(y64); - } - - if (z->type.type != PS_TYPE_F64) { - psFree(z64); - } - - if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { - psFree(fErr64); - } + psError(PS_ERR_UNKNOWN, true, "3-D Chebyshev polynomial vector fitting has not been implemented. Returning NULL.\n"); + result = false; + break; + default: psError(PS_ERR_UNKNOWN, true, "Incorrect polynomial type. Returning NULL.\n"); - } - + result = false; + break; + } // Free psVectors that were created for NULL arguments. - if (f->type.type != PS_TYPE_F64) { - psFree(f64); - } - - if (x->type.type != PS_TYPE_F64) { - psFree(x64); - } - - if (y->type.type != PS_TYPE_F64) { - psFree(y64); - } - - if (z->type.type != PS_TYPE_F64) { - psFree(z64); - } - - if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { - psFree(fErr64); - } - - return(poly); + PS_FREE_TEMP_F64_VECTOR (f, f64); + PS_FREE_TEMP_F64_VECTOR (x, x64); + PS_FREE_TEMP_F64_VECTOR (y, y64); + PS_FREE_TEMP_F64_VECTOR (z, z64); + PS_FREE_TEMP_F64_VECTOR (fErr, fErr64); + + return result; } -psPolynomial3D *psVectorClipFitPolynomial3D( +bool psVectorClipFitPolynomial3D( psPolynomial3D *poly, psStats *stats, @@ -1738,29 +1589,44 @@ { psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL); - PS_ASSERT_PTR_NON_NULL(stats, NULL); - PS_ASSERT_VECTOR_NON_NULL(mask, NULL); - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL); - if (mask != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); - } - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); - PS_ASSERT_VECTOR_TYPE(x, f->type.type, NULL); - - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); - PS_ASSERT_VECTOR_TYPE(y, f->type.type, NULL); - - PS_ASSERT_VECTOR_NON_NULL(z, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL); - PS_ASSERT_VECTOR_TYPE(z, f->type.type, NULL); + PS_ASSERT_POLY_NON_NULL(poly, false); + PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, false); + PS_ASSERT_PTR_NON_NULL(stats, false); + PS_ASSERT_VECTOR_NON_NULL(mask, false); + PS_ASSERT_VECTOR_NON_NULL(f, false); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false); + + PS_ASSERT_VECTOR_NON_NULL(x, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false); + PS_ASSERT_VECTOR_TYPE(x, f->type.type, false); + + PS_ASSERT_VECTOR_NON_NULL(y, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false); + PS_ASSERT_VECTOR_TYPE(y, f->type.type, false); + + PS_ASSERT_VECTOR_NON_NULL(z, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, false); + PS_ASSERT_VECTOR_TYPE(z, f->type.type, false); + + PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, false); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false); if (fErr != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); - PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, false); + PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, false); + } + + // the user supplies one of various stats option pairs, + // determine the desired mean and stdev STATS options: + // XXX enforce consistency? + // XXX psStatsGetValue() probably has inverted precedence + psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN_V2); + psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV_V2); + if (!meanOption) { + psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected"); + return false; + } + if (!stdevOption) { + psError(PS_ERR_UNKNOWN, true, "no valid stdev stats option selected"); + return false; } @@ -1778,13 +1644,6 @@ minClipSigma = fabs(stats->clipSigma); } - psVector *fit = NULL; psVector *resid = psVectorAlloc(f->n, PS_TYPE_F64); - // eventual expansion: user supplies one of various stats option pairs, - // eg (SAMPLE_MEAN | SAMPLE_STDEV) and the correct pair is used to - // evaluate the clipping sigma - // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used - stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - stats->options |= (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); psTrace("psLib.math", 4, "stats->clipIter is %d\n", stats->clipIter); psTrace("psLib.math", 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma); @@ -1801,16 +1660,14 @@ } - poly = psVectorFitPolynomial3D(poly, mask, maskValue, f, fErr, x, y, z); - if (poly == NULL) { + if (!psVectorFitPolynomial3D(poly, mask, maskValue, f, fErr, x, y, z)) { psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to the data. Returning NULL.\n"); psFree(resid) - psFree(fit) - return(NULL); - } - fit = psPolynomial3DEvalVector(poly, x, y, z); + return false; + } + psVector *fit = psPolynomial3DEvalVector(poly, x, y, z); if (fit == NULL) { psError(PS_ERR_UNKNOWN, false, "Could not call psPolynomial3DEvalVector(). Returning NULL.\n"); psFree(resid) - return(NULL); + return false; } for (psS32 i = 0 ; i < f->n ; i++) { @@ -1837,21 +1694,14 @@ psFree(resid) psFree(fit) - return(NULL); - } - - # if (USE_ROBUST_STATS_FOR_CLIPPING) - psTrace("psLib.math", 5, "Median is %f\n", stats->robustMedian); - psTrace("psLib.math", 5, "Stdev is %f\n", stats->robustStdev); - psF32 minClipValue = -minClipSigma*stats->robustStdev; - psF32 maxClipValue = +maxClipSigma*stats->robustStdev; - psF32 clipMedian = stats->robustMedian; - # else - - psTrace("psLib.math", 5, "Median is %f\n", stats->sampleMedian); - psTrace("psLib.math", 5, "Stdev is %f\n", stats->sampleStdev); - psF32 minClipValue = -minClipSigma*stats->robustStdev; - psF32 maxClipValue = +maxClipSigma*stats->robustStdev; - psF32 clipMedian = stats->sampleMedian; - # endif + return false; + } + + double meanValue = psStatsGetValue (stats, meanOption); + double stdevValue = psStatsGetValue (stats, stdevOption); + + psTrace("psLib.math", 5, "Mean is %f\n", meanValue); + psTrace("psLib.math", 5, "Stdev is %f\n", stdevValue); + psF32 minClipValue = -minClipSigma*stdevValue; + psF32 maxClipValue = +maxClipSigma*stdevValue; // set mask if pts are not valid @@ -1863,5 +1713,5 @@ } - if ((resid->data.F64[i] - clipMedian > maxClipValue) || (resid->data.F64[i] - clipMedian < minClipValue)) { + if ((resid->data.F64[i] - meanValue > maxClipValue) || (resid->data.F64[i] - meanValue < minClipValue)) { if (f->type.type == PS_TYPE_F64) { psTrace("psLib.math", 6, "Masking element %d (%f). resid->data.F64[%d] is %f\n", @@ -1879,6 +1729,6 @@ Nkeep++; } - psTrace("psLib.math", 6, "keeping %d of %ld pts for fit\n", Nkeep, x->n); + stats->clippedNvalues = Nkeep; psFree(fit); } @@ -1886,11 +1736,6 @@ psFree(resid); - if (poly == NULL) { - psError(PS_ERR_UNKNOWN, true, "Could not fit a polynomial to the data. Returning NULL.\n"); - return(NULL); - } - psTrace("psLib.math", 3, "---- %s() end ----\n", __func__); - return(poly); + return true; } @@ -1906,5 +1751,5 @@ *****************************************************************************/ -static psPolynomial4D* VectorFitPolynomial4DOrd( +static bool VectorFitPolynomial4DOrd( psPolynomial4D* myPoly, const psVector* mask, @@ -1918,30 +1763,30 @@ { psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); - PS_ASSERT_POLY_NON_NULL(myPoly, NULL); - PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL); - PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, NULL); - PS_ASSERT_INT_NONNEGATIVE(myPoly->nZ, NULL); - PS_ASSERT_INT_NONNEGATIVE(myPoly->nT, NULL); - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, NULL); + PS_ASSERT_POLY_NON_NULL(myPoly, false); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, false); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, false); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nZ, false); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nT, false); + PS_ASSERT_VECTOR_NON_NULL(f, false); + PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, false); if (fErr != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(y, fErr, NULL); - PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, NULL); - } - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); - PS_ASSERT_VECTOR_NON_NULL(z, NULL); - PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F64, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL); - PS_ASSERT_VECTOR_NON_NULL(t, NULL); - PS_ASSERT_VECTOR_TYPE(t, PS_TYPE_F64, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(y, fErr, false); + PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, false); + } + PS_ASSERT_VECTOR_NON_NULL(x, false); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false); + PS_ASSERT_VECTOR_NON_NULL(y, false); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false); + PS_ASSERT_VECTOR_NON_NULL(z, false); + PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F64, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, false); + PS_ASSERT_VECTOR_NON_NULL(t, false); + PS_ASSERT_VECTOR_TYPE(t, PS_TYPE_F64, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, false); if (mask) { - PS_ASSERT_VECTORS_SIZE_EQUAL(y, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(y, mask, false); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false); } @@ -1959,9 +1804,8 @@ if (!psImageInit(A, 0.0) || !psVectorInit(B, 0.0)) { psError(PS_ERR_UNKNOWN, false, "Could initialize data structures A, B. Returning NULL.\n"); - psFree(myPoly); psFree(A); psFree(B); psTrace("psLib.math", 4, "---- %s() End ----\n", __func__); - return(NULL); + return false; } @@ -2051,8 +1895,6 @@ // The GaussJordan version was overflowing, so I'm using LUD. if (!psMatrixGJSolve(A, B)) { - psFree(A); - psFree(B); psError(PS_ERR_UNKNOWN, false, "Failed to perform GaussJordan elimination.\n"); - return(NULL); + goto escape_GJ; } @@ -2076,12 +1918,10 @@ if (ALUD == NULL) { psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix. Returning NULL.\n"); - psFree(myPoly); - myPoly = NULL; + goto escape_LUD; } 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; + goto escape_LUD; } else { // select the appropriate solution entries @@ -2100,16 +1940,23 @@ } } - psFree(ALUD); psFree(coeffs); psFree(outPerm); - - } - + } psFree(A); psFree(B); psTrace("psLib.math", 4, "---- %s() end ----\n", __func__); - return (myPoly); + return true; + +escape_LUD: + psFree(A); + psFree(B); + return false; + +escape_GJ: + psFree(A); + psFree(B); + return false; } @@ -2120,5 +1967,5 @@ via vector conversion only. *****************************************************************************/ -psPolynomial4D *psVectorFitPolynomial4D( +bool psVectorFitPolynomial4D( psPolynomial4D *poly, const psVector *mask, @@ -2131,170 +1978,73 @@ const psVector *t) { - // Internal pointers for possibly NULL or mis-typed vectors. - psVector *x64 = NULL; - psVector *y64 = NULL; - psVector *z64 = NULL; - psVector *t64 = NULL; - psVector *f64 = NULL; + PS_ASSERT_POLY_NON_NULL(poly, false); + PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, false); + + PS_ASSERT_VECTOR_NON_NULL(f, false); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false); + PS_ASSERT_VECTOR_NON_NULL(x, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false); + PS_ASSERT_VECTOR_NON_NULL(y, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false); + PS_ASSERT_VECTOR_NON_NULL(z, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, false); + PS_ASSERT_VECTOR_NON_NULL(t, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, false); + if (mask) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, false); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false); + } + if (fErr != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, false); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, false); + } + + // Convert input vectors to F64 if necessary. + psVector *f64 = (f->type.type == PS_TYPE_F64) ? (psVector *) f : psVectorCopy(NULL, f, PS_TYPE_F64); + psVector *x64 = (x->type.type == PS_TYPE_F64) ? (psVector *) x : psVectorCopy(NULL, x, PS_TYPE_F64); + psVector *y64 = (y->type.type == PS_TYPE_F64) ? (psVector *) y : psVectorCopy(NULL, y, PS_TYPE_F64); + psVector *z64 = (z->type.type == PS_TYPE_F64) ? (psVector *) z : psVectorCopy(NULL, z, PS_TYPE_F64); + psVector *t64 = (t->type.type == PS_TYPE_F64) ? (psVector *) t : psVectorCopy(NULL, t, PS_TYPE_F64); + psVector *fErr64 = NULL; - - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL); - - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL); - if (mask) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); - } - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); - PS_ASSERT_VECTOR_NON_NULL(z, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL); - PS_ASSERT_VECTOR_NON_NULL(t, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL); if (fErr != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL); - } - - - // - // Convert input vector to F64 if necessary. - // - if (f->type.type != PS_TYPE_F64) { - f64 = psVectorCopy(NULL, f, PS_TYPE_F64); - } else { - f64 = (psVector *) f; - } - if (x->type.type != PS_TYPE_F64) { - x64 = psVectorCopy(NULL, x, PS_TYPE_F64); - } else { - x64 = (psVector *) x; - } - if (y->type.type != PS_TYPE_F64) { - y64 = psVectorCopy(NULL, y, PS_TYPE_F64); - } else { - y64 = (psVector *) y; - } - if (z->type.type != PS_TYPE_F64) { - z64 = psVectorCopy(NULL, z, PS_TYPE_F64); - } else { - z64 = (psVector *) z; - } - if (t->type.type != PS_TYPE_F64) { - t64 = psVectorCopy(NULL, t, PS_TYPE_F64); - } else { - t64 = (psVector *) t; - } - // - // fErr - // - if (fErr != NULL) { - if (fErr->type.type != PS_TYPE_F64) { - fErr64 = psVectorCopy(NULL, fErr, PS_TYPE_F64); - } else { - fErr64 = (psVector *) fErr; - } - } - - if (poly->type == PS_POLYNOMIAL_ORD) { - poly = VectorFitPolynomial4DOrd(poly, mask, maskValue, f64, fErr64, x64, y64, z64, t64); - if (poly == NULL) { + fErr64 = (fErr->type.type == PS_TYPE_F64) ? (psVector *) fErr : psVectorCopy(NULL, fErr, PS_TYPE_F64); + } + + bool result = true; + + switch (poly->type) { + case PS_POLYNOMIAL_ORD: + result = VectorFitPolynomial4DOrd(poly, mask, maskValue, f64, fErr64, x64, y64, z64, t64); + if (!result) { psError(PS_ERR_UNKNOWN, true, "Could not fit polynomial. Returning NULL.\n"); - // Free psVectors that were created for NULL arguments. - if (f->type.type != PS_TYPE_F64) { - psFree(f64); - } - - if (x->type.type != PS_TYPE_F64) { - psFree(x64); - } - - if (y->type.type != PS_TYPE_F64) { - psFree(y64); - } - - if (z->type.type != PS_TYPE_F64) { - psFree(z64); - } - - if (t->type.type != PS_TYPE_F64) { - psFree(t64); - } - - if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { - psFree(fErr64); - } - return(NULL); - } - } else if (poly->type == PS_POLYNOMIAL_CHEB) { + } + break; + case PS_POLYNOMIAL_CHEB: if (mask != NULL) { psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring mask and maskValue with Chebyshev polynomials.\n"); } - psLogMsg(__func__, PS_LOG_WARN, "WARNING: 4-D Chebyshev polynomial vector fitting has not been implemented. Returning NULL.\n"); - psFree(poly); - poly = NULL; - } else { - // Free psVectors that were created for NULL arguments. - if (f->type.type != PS_TYPE_F64) { - psFree(f64); - } - - if (x->type.type != PS_TYPE_F64) { - psFree(x64); - } - - if (y->type.type != PS_TYPE_F64) { - psFree(y64); - } - - if (z->type.type != PS_TYPE_F64) { - psFree(z64); - } - - if (t->type.type != PS_TYPE_F64) { - psFree(t64); - } - - if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { - psFree(fErr64); - } + psError(PS_ERR_UNKNOWN, true, "4-D Chebyshev polynomial vector fitting has not been implemented. Returning NULL.\n"); + result = false; + break; + default: psError(PS_ERR_UNKNOWN, true, "Incorrect polynomial type. Returning NULL.\n"); - } - + result = false; + break; + } // Free psVectors that were created for NULL arguments. - if (f->type.type != PS_TYPE_F64) { - psFree(f64); - } - - if (x->type.type != PS_TYPE_F64) { - psFree(x64); - } - - if (y->type.type != PS_TYPE_F64) { - psFree(y64); - } - - if (z->type.type != PS_TYPE_F64) { - psFree(z64); - } - - if (t->type.type != PS_TYPE_F64) { - psFree(t64); - } - - if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { - psFree(fErr64); - } - - return(poly); + PS_FREE_TEMP_F64_VECTOR (f, f64); + PS_FREE_TEMP_F64_VECTOR (x, x64); + PS_FREE_TEMP_F64_VECTOR (y, y64); + PS_FREE_TEMP_F64_VECTOR (z, z64); + PS_FREE_TEMP_F64_VECTOR (t, t64); + PS_FREE_TEMP_F64_VECTOR (fErr, fErr64); + + return result; } -psPolynomial4D *psVectorClipFitPolynomial4D( +bool psVectorClipFitPolynomial4D( psPolynomial4D *poly, psStats *stats, @@ -2309,33 +2059,48 @@ { psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL); - PS_ASSERT_PTR_NON_NULL(stats, NULL); - PS_ASSERT_VECTOR_NON_NULL(mask, NULL); - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL); - if (mask) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); - } - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); - PS_ASSERT_VECTOR_TYPE(x, f->type.type, NULL); - - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); - PS_ASSERT_VECTOR_TYPE(y, f->type.type, NULL); - - PS_ASSERT_VECTOR_NON_NULL(z, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL); - PS_ASSERT_VECTOR_TYPE(z, f->type.type, NULL); - - PS_ASSERT_VECTOR_NON_NULL(t, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL); - PS_ASSERT_VECTOR_TYPE(t, f->type.type, NULL); + PS_ASSERT_POLY_NON_NULL(poly, false); + PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, false); + PS_ASSERT_PTR_NON_NULL(stats, false); + PS_ASSERT_VECTOR_NON_NULL(mask, false); + PS_ASSERT_VECTOR_NON_NULL(f, false); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false); + + PS_ASSERT_VECTOR_NON_NULL(x, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false); + PS_ASSERT_VECTOR_TYPE(x, f->type.type, false); + + PS_ASSERT_VECTOR_NON_NULL(y, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false); + PS_ASSERT_VECTOR_TYPE(y, f->type.type, false); + + PS_ASSERT_VECTOR_NON_NULL(z, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, false); + PS_ASSERT_VECTOR_TYPE(z, f->type.type, false); + + PS_ASSERT_VECTOR_NON_NULL(t, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, false); + PS_ASSERT_VECTOR_TYPE(t, f->type.type, false); + + PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, false); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false); if (fErr != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); - PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, false); + PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, false); + } + + // the user supplies one of various stats option pairs, + // determine the desired mean and stdev STATS options: + // XXX enforce consistency? + // XXX psStatsGetValue() probably has inverted precedence + psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN_V2); + psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV_V2); + if (!meanOption) { + psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected"); + return false; + } + if (!stdevOption) { + psError(PS_ERR_UNKNOWN, true, "no valid stdev stats option selected"); + return false; } @@ -2353,13 +2118,6 @@ minClipSigma = fabs(stats->clipSigma); } - psVector *fit = NULL; psVector *resid = psVectorAlloc(f->n, PS_TYPE_F64); - // eventual expansion: user supplies one of various stats option pairs, - // eg (SAMPLE_MEAN | SAMPLE_STDEV) and the correct pair is used to - // evaluate the clipping sigma - // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used - stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - stats->options |= (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); psTrace("psLib.math", 4, "stats->clipIter is %d\n", stats->clipIter); psTrace("psLib.math", 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma); @@ -2376,17 +2134,15 @@ } - poly = psVectorFitPolynomial4D (poly, mask, maskValue, f, fErr, x, y, z, t); - if (poly == NULL) { + if (!psVectorFitPolynomial4D (poly, mask, maskValue, f, fErr, x, y, z, t)) { psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to the data. Returning NULL.\n"); psFree(resid) - psFree(fit) - return(NULL); - } - - fit = psPolynomial4DEvalVector (poly, x, y, z, t); + return false; + } + + psVector *fit = psPolynomial4DEvalVector (poly, x, y, z, t); if (fit == NULL) { psError(PS_ERR_UNKNOWN, false, "Could not call psPolynomial4DEvalVector(). Returning NULL.\n"); psFree(resid) - return(NULL); + return false; } for (psS32 i = 0 ; i < f->n ; i++) { @@ -2413,20 +2169,14 @@ psFree(resid) psFree(fit) - return(NULL); - } - # if (USE_ROBUST_STATS_FOR_CLIPPING) - psTrace("psLib.math", 5, "Median is %f\n", stats->robustMedian); - psTrace("psLib.math", 5, "Stdev is %f\n", stats->robustStdev); - psF32 minClipValue = -minClipSigma*stats->robustStdev; - psF32 maxClipValue = +maxClipSigma*stats->robustStdev; - psF32 clipMedian = stats->robustMedian; - # else - - psTrace("psLib.math", 5, "Median is %f\n", stats->sampleMedian); - psTrace("psLib.math", 5, "Stdev is %f\n", stats->sampleStdev); - psF32 minClipValue = -minClipSigma*stats->robustStdev; - psF32 maxClipValue = +maxClipSigma*stats->robustStdev; - psF32 clipMedian = stats->sampleMedian; - # endif + return false; + } + + double meanValue = psStatsGetValue (stats, meanOption); + double stdevValue = psStatsGetValue (stats, stdevOption); + + psTrace("psLib.math", 5, "Mean is %f\n", meanValue); + psTrace("psLib.math", 5, "Stdev is %f\n", stdevValue); + psF32 minClipValue = -minClipSigma*stdevValue; + psF32 maxClipValue = +maxClipSigma*stdevValue; // set mask if pts are not valid @@ -2438,5 +2188,5 @@ } - if ((resid->data.F64[i] - clipMedian > maxClipValue) || (resid->data.F64[i] - clipMedian < minClipValue)) { + if ((resid->data.F64[i] - meanValue > maxClipValue) || (resid->data.F64[i] - meanValue < minClipValue)) { if (f->type.type == PS_TYPE_F64) { psTrace("psLib.math", 6, "Masking element %d (%f). resid->data.F64[%d] is %f\n", @@ -2452,9 +2202,8 @@ continue; } - Nkeep++; } - psTrace("psLib.math", 6, "keeping %d of %ld pts for fit\n", Nkeep, x->n); + stats->clippedNvalues = Nkeep; psFree (fit); } @@ -2462,10 +2211,5 @@ psFree (resid); - if (poly == NULL) { - psError(PS_ERR_UNKNOWN, true, "Could not fit a polynomial to the data. Returning NULL.\n"); - return(NULL); - } - psTrace("psLib.math", 3, "---- %s() end ----\n", __func__); - return(poly); + return true; }