Index: /trunk/psLib/src/math/psConstants.h =================================================================== --- /trunk/psLib/src/math/psConstants.h (revision 6436) +++ /trunk/psLib/src/math/psConstants.h (revision 6437) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.84 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-01-26 21:10:22 $ + * @version $Revision: 1.85 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-02-17 00:56:48 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -169,5 +169,5 @@ } -// Return an error if the arg is lies outside the supplied range. +// Return an error if the arg lies outside the supplied range. #define PS_ASSERT_FLOAT_WITHIN_RANGE(NAME, LOWER, UPPER, RVAL) \ if ((NAME) < (LOWER) || (NAME) > (UPPER)) { \ @@ -186,5 +186,4 @@ } -// Return an error if the arg lies outside the supplied range #define PS_ASSERT_LONG_WITHIN_RANGE(NAME, LOWER, UPPER, RVAL) \ if ((NAME) < (LOWER) || (NAME) > (UPPER)) { \ @@ -519,27 +518,4 @@ } \ -/* XXX: This is not correct - #define PS_POLY_2D_DECLARE_ALLOC_STATIC(NAME, ORDER, TYPE) \ - static psPolynomial2D *(NAME) = NULL; \ - if ((NAME) == NULL) { \ - (NAME) = psPolynomial2DAlloc(ORDER, TYPE); \ - p_psMemSetPersistent((NAME), true); \ - } \ - - #define PS_POLY_3D_DECLARE_ALLOC_STATIC(NAME, ORDER, TYPE) \ - static psPolynomial3D *(NAME) = NULL; \ - if ((NAME) == NULL) { \ - (NAME) = psPolynomial3DAlloc(ORDER, TYPE); \ - p_psMemSetPersistent((NAME), true); \ - } \ - - #define PS_POLY_4D_DECLARE_ALLOC_STATIC(NAME, ORDER, TYPE) \ - static psPolynomial4D *(NAME) = NULL; \ - if ((NAME) == NULL) { \ - (NAME) = psPolynomial4DAlloc(ORDER, TYPE); \ - p_psMemSetPersistent((NAME), true); \ - } \ -*/ - #define PS_POLY_1D_D_DECLARE_ALLOC_STATIC(NAME, ORDER, TYPE) \ static psPolynomial1D *(NAME) = NULL; \ @@ -548,27 +524,4 @@ p_psMemSetPersistent((NAME), true); \ } \ - -/* XXX: This is not correct - #define PS_POLY_2D_D_DECLARE_ALLOC_STATIC(NAME, ORDER, TYPE) \ - static psPolynomial2D *(NAME) = NULL; \ - if ((NAME) == NULL) { \ - (NAME) = psPolynomial2DAlloc(ORDER, TYPE); \ - p_psMemSetPersistent((NAME), true); \ - } \ - - #define PS_POLY_3D_D_DECLARE_ALLOC_STATIC(NAME, ORDER, TYPE) \ - static psPolynomial3D *(NAME) = NULL; \ - if ((NAME) == NULL) { \ - (NAME) = psPolynomial3DAlloc(ORDER, TYPE); \ - p_psMemSetPersistent((NAME), true); \ - } \ - - #define PS_POLY_4D_D_DECLARE_ALLOC_STATIC(NAME, ORDER, TYPE) \ - static psPolynomial4D *(NAME) = NULL; \ - if ((NAME) == NULL) { \ - (NAME) = psPolynomial4DAlloc(ORDER, TYPE); \ - p_psMemSetPersistent((NAME), true); \ - } \ -*/ #define PS_POLY_PRINT_1D(NAME) \ Index: /trunk/psLib/src/math/psMathUtils.h =================================================================== --- /trunk/psLib/src/math/psMathUtils.h (revision 6436) +++ /trunk/psLib/src/math/psMathUtils.h (revision 6437) @@ -7,6 +7,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-02-02 21:09:07 $ + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-02-17 00:56:48 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -51,7 +51,4 @@ ); -// XXX: Create a single, generic, version of the vector normalize function. -// XXX: Ask IfA for a public psLib function. - psS32 p_psNormalizeVectorRange(psVector* myData, psF64 outLow, Index: /trunk/psLib/src/math/psMinimizeLMM.h =================================================================== --- /trunk/psLib/src/math/psMinimizeLMM.h (revision 6436) +++ /trunk/psLib/src/math/psMinimizeLMM.h (revision 6437) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-01-27 20:08:58 $ + * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-02-17 00:56:48 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -115,5 +115,5 @@ const psArray *x, ///< Measurement ordinates of multiple vectors const psVector *y, ///< Measurement coordinates - const psVector *yErr, ///< Errors in the measurement coordinates + const psVector *yWt, ///< Errors in the measurement coordinates psMinimizeLMChi2Func func ///< Specified function ); Index: /trunk/psLib/src/math/psPolynomial.c =================================================================== --- /trunk/psLib/src/math/psPolynomial.c (revision 6436) +++ /trunk/psLib/src/math/psPolynomial.c (revision 6437) @@ -4,9 +4,9 @@ * routines. * -* This file will hold the functions for allocated, freeing, and evaluating +* This file will hold the routiness for allocating, freeing, and evaluating * polynomials. It also contains a Gaussian functions. * -* @version $Revision: 1.142 $ $Name: not supported by cvs2svn $ -* @date $Date: 2006-02-07 23:36:15 $ +* @version $Revision: 1.143 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-02-17 00:56:48 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -159,5 +159,4 @@ used frequently and the data structure created here does not contain the outer coefficients of the Chebyshev polynomials. - *****************************************************************************/ psPolynomial1D **p_psCreateChebyshevPolys(psS32 numPolys) @@ -165,40 +164,31 @@ PS_ASSERT_INT_LARGER_THAN_OR_EQUAL(numPolys, 1, NULL); - if (0) { - psPolynomial1D **chebPolys = (psPolynomial1D **) psAlloc(numPolys * sizeof(psPolynomial1D *)); - for (psS32 i = 0; i < numPolys; i++) { - chebPolys[i] = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, i); - chebPolys[i]->coeff[i] = 1; - } - return (chebPolys); - } else { - psPolynomial1D **chebPolys = (psPolynomial1D **) psAlloc(numPolys * sizeof(psPolynomial1D *)); - for (psS32 i = 0; i < numPolys; i++) { - chebPolys[i] = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, i); - } - - // Create the Chebyshev polynomials. - // Polynomial i has i-th order. - chebPolys[0]->coeff[0] = 1.0; - if (numPolys >= 2) { - chebPolys[1]->coeff[1] = 1.0; - - for (psS32 i = 2; i < numPolys; i++) { - for (psS32 j = 0; j < chebPolys[i - 1]->nX+1; j++) { - chebPolys[i]->coeff[j + 1] = 2.0 * chebPolys[i - 1]->coeff[j]; - } - for (psS32 j = 0; j < chebPolys[i - 2]->nX+1; j++) { - chebPolys[i]->coeff[j] -= chebPolys[i - 2]->coeff[j]; - } - } - } - - if (psTraceGetLevel(__func__) >= 6) { - for (psS32 j = 0; j < numPolys; j++) { - PS_POLY_PRINT_1D(chebPolys[j]); - } - } - return (chebPolys); - } + psPolynomial1D **chebPolys = (psPolynomial1D **) psAlloc(numPolys * sizeof(psPolynomial1D *)); + for (psS32 i = 0; i < numPolys; i++) { + chebPolys[i] = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, i); + } + + // Create the Chebyshev polynomials. + // Polynomial i has i-th order. + chebPolys[0]->coeff[0] = 1.0; + if (numPolys >= 2) { + chebPolys[1]->coeff[1] = 1.0; + + for (psS32 i = 2; i < numPolys; i++) { + for (psS32 j = 0; j < chebPolys[i - 1]->nX+1; j++) { + chebPolys[i]->coeff[j + 1] = 2.0 * chebPolys[i - 1]->coeff[j]; + } + for (psS32 j = 0; j < chebPolys[i - 2]->nX+1; j++) { + chebPolys[i]->coeff[j] -= chebPolys[i - 2]->coeff[j]; + } + } + } + + if (psTraceGetLevel(__func__) >= 6) { + for (psS32 j = 0; j < numPolys; j++) { + PS_POLY_PRINT_1D(chebPolys[j]); + } + } + return (chebPolys); } @@ -223,6 +213,4 @@ for (loop_x = 0; loop_x < poly->nX+1; loop_x++) { if (poly->mask[loop_x] == 0) { - // XXX: If you set the tracelevel to 10 here, and later set the tracelevel to - // 2 or higher in the test code, you get seg faults. psTrace(__func__, 8, "polysum+= sum*coeff [%lf+= (%lf * %lf)\n", polySum, xSum, poly->coeff[loop_x]); @@ -238,5 +226,5 @@ // XXX: You can do this without having to psAlloc() vector d. // XXX: How does the mask vector effect Crenshaw's formula? -// XXX: We assume that x is scaled between -1.0 and 1.0; +// NOTE: We assume that x is scaled between -1.0 and 1.0; // XXX: Create a faster version for low-order Chebyshevs. static psF64 chebPolynomial1DEval( @@ -298,5 +286,5 @@ psFree(d); } else { - // This is old code that does not use Clenshaw's formula. Get rid of it. + // XXX: This is old code that does not use Clenshaw's formula. Get rid of it. psPolynomial1D **chebPolys = p_psCreateChebyshevPolys(1 + poly->nX); Index: /trunk/psLib/src/math/psSpline.c =================================================================== --- /trunk/psLib/src/math/psSpline.c (revision 6436) +++ /trunk/psLib/src/math/psSpline.c (revision 6437) @@ -6,6 +6,6 @@ * This file contains the routines that allocate, free, and evaluate splines. * -* @version $Revision: 1.135 $ $Name: not supported by cvs2svn $ -* @date $Date: 2006-02-02 21:09:08 $ +* @version $Revision: 1.136 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-02-17 00:56:48 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -101,5 +101,5 @@ The first and second derivatives at the endpoints were previously undefined in -the SDR. From bugzilla #???, they are required to be 0.0, impementing natural +the SDR. From bugzilla #???, they are required to be 0.0, implementing natural splines. @@ -228,5 +228,5 @@ // The following code ensures that xPtr and yPtr points to a psF32 psVector. // - // XXX: When you debug, use the vector copy and create routines here: + // XXX: Use the vector copy and create routines here: // @@ -345,5 +345,5 @@ /***************************************************************************** psSpline1DEval(): this routine takes an existing spline of arbitrary order -and an independent x value. Each determines which spline that x corresponds +and an independent x value. It determines which spline that x corresponds to by doing a bracket disection on the knots of the spline data structure (vectorBinDisectF32()). Then it evaluates the spline at that x location @@ -360,5 +360,5 @@ float x) { - psTrace(__func__, 3, "---- %s(%f) begin ----\n", __func__, x); + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); PS_ASSERT_PTR_NON_NULL(spline, NAN); PS_ASSERT_INT_NONNEGATIVE(spline->n, NAN); @@ -366,32 +366,27 @@ psS32 n = spline->n; + if ((x < spline->knots->data.F32[0]) || (x > spline->knots->data.F32[n-1])) { + // If x is outside the range of spline->knots, generate a warning + // message, then return the left, or right, endpoint. + psLogMsg(__func__, PS_LOG_WARN, + "psSpline1DEval(): x ordinate (%f) is outside the spline range (%f - %f) (%d).", + x, spline->knots->data.F32[0], spline->knots->data.F32[n-1], n); + + psS32 binNum = (x < spline->knots->data.F32[0]) ? 0 : n-1; + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(psPolynomial1DEval(spline->spline[binNum], x)); + } + psScalar tmpScalar; tmpScalar.type.type = PS_TYPE_F32; tmpScalar.data.F32 = x; psS32 binNum = p_psVectorBinDisect(spline->knots, &tmpScalar); - if (binNum < 0) { - // - // If x is outside the range of spline->knots, generate a warning - // message, then return the left, or right, endpoint. - // - psLogMsg(__func__, PS_LOG_WARN, - "psSpline1DEval(): x ordinate (%f) is outside the spline range (%f - %f) (%d).", - x, spline->knots->data.F32[0], spline->knots->data.F32[n-1], n); - - if (x < spline->knots->data.F32[0]) { - psF32 rcF32 = psPolynomial1DEval(spline->spline[0], x); - psTrace(__func__, 3, "---- %s(%f) end ----\n", __func__, rcF32); - return(rcF32); - } else if (x > spline->knots->data.F32[n-1]) { - psF32 rcF32 = psPolynomial1DEval(spline->spline[n-1], x); - psTrace(__func__, 3, "---- %s(%f) end ----\n", __func__, rcF32); - return(rcF32); - } - } - - psF32 rcF32 = psPolynomial1DEval(spline->spline[binNum], x); - psTrace(__func__, 3, "---- %s(%f) end ----\n", __func__, rcF32); - return(rcF32); + psError(PS_ERR_UNKNOWN, false, "Could not perform bin dissection on spline->knots.\n"); + return(NAN); + } + + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(psPolynomial1DEval(spline->spline[binNum], x)); } Index: /trunk/psLib/src/math/psStats.c =================================================================== --- /trunk/psLib/src/math/psStats.c (revision 6436) +++ /trunk/psLib/src/math/psStats.c (revision 6437) @@ -16,6 +16,6 @@ * use ->min and ->max (PS_STAT_USE_RANGE) * - * @version $Revision: 1.167 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-02-03 22:05:22 $ + * @version $Revision: 1.168 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-02-17 00:56:48 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -440,11 +440,10 @@ unmasked element within the specified min/max range). Otherwise, return "false". - -XXX: Can you use psVectorCountPixelMask here? *****************************************************************************/ -bool p_psVectorCheckNonEmpty(const psVector* myVector, - const psVector* maskVector, - psU32 maskVal, - psStats* stats) +bool p_psVectorCheckNonEmpty( + const psVector* myVector, + const psVector* maskVector, + psU32 maskVal, + psStats* stats) { psTrace(__func__, 4, "---- %s() begin ----\n", __func__); @@ -496,6 +495,4 @@ number of non-masked pixels in the vector that fall within the min/max range, if specified. - -XXX: Can you use psVectorCountPixelMask here? *****************************************************************************/ psS32 p_psVectorNValues(const psVector* myVector, @@ -644,9 +641,4 @@ robustHistogram with a Gaussian of width sigma. It returns a psVector of the smoothed data. - -XXX: Only PS_TYPE_F32 is supported. - -XXX: Write a general routine which smoothes a psVector. This routine should -call that. Is that possible? *****************************************************************************/ psVector *p_psVectorSmoothHistGaussian( @@ -985,6 +977,4 @@ -1: error -2: warning - -XXX: Use static vectors for tmpMask. *****************************************************************************/ psS32 p_psVectorClippedStats(const psVector* myVector, @@ -1224,8 +1214,5 @@ // - // Ensure that the y values are monotonic. - // - // XXX: This routine should probably be rewritten in a more general fashion - // so that the following checks are not necessary. + // Ensure that the y value lies within range of the y values. // if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2])) || @@ -1236,4 +1223,7 @@ } + // + // Ensure that the y values are monotonic. + // if (((y->data.F64[0] < y->data.F64[1]) && !(y->data.F64[1] <= y->data.F64[2])) || ((y->data.F64[0] > y->data.F64[1]) && !(y->data.F64[1] >= y->data.F64[2]))) { @@ -1289,6 +1279,5 @@ } else if (binNum == (xVec->n - 2)) { // XXX: Is this right? - tmpFloat = 0.5 * (xVec->data.F32[binNum] + - xVec->data.F32[binNum + 1]); + tmpFloat = 0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum + 1]); } } @@ -1677,5 +1666,5 @@ tmpScalar.data.F32 = totalDataPoints * 0.25f; if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) { - // XXX: Special case where its in first bin. Must code last bin. + // Special case where its in first bin. The last bin should be handled automatically. binLo25 = 0; } else { @@ -1684,5 +1673,5 @@ tmpScalar.data.F32 = totalDataPoints * 0.75f; if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) { - // XXX: Special case where its in first bin. Must code last bin. + // Special case where its in first bin. The last bin should be handled automatically. binHi25 = 0; } else { @@ -1704,5 +1693,4 @@ // Interpolate to find these two positions exactly: these are the upper // and lower quartile positions. - // XXX: Check for errors. // psF32 binLo25F32 = fitQuadraticSearchForYThenReturnX( @@ -1711,4 +1699,5 @@ binLo25, totalDataPoints * 0.25f); + psF32 binHi25F32 = fitQuadraticSearchForYThenReturnX( *(psVector* *)&cumulativeRobustHistogram->bounds, @@ -1716,8 +1705,18 @@ binHi25, totalDataPoints * 0.75f); + + if (isnan(binLo25F32) || isnan(binHi25F32)) { + psError(PS_ERR_UNKNOWN, false, "could not determine the robustUQ: fitQuadraticSearchForYThenReturnX() returned a NAN.\n"); + psFree(tmpStatsMinMax); + psFree(robustHistogram); + psFree(cumulativeRobustHistogram); + psFree(tmpMaskVec); + psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); + return(1); + } + stats->robustLQ = binLo25F32; stats->robustUQ = binHi25F32; psTrace(__func__, 6, "The 25 and 75 percent data point exact positions are (%f, %f).\n", binLo25F32, binHi25F32); - psS32 N50 = 0; for (psS32 i = 0 ; i < myVector->n ; i++) { @@ -1863,15 +1862,16 @@ } - // XXX: Use the min/max routines for this - psF32 minY = FLT_MAX; - psF32 maxY = -FLT_MAX; - for (psS32 i = 0 ; i < y->n ; i++) { - if (y->data.F32[i] > maxY) { - maxY = y->data.F32[i]; - } - if (y->data.F32[i] < minY) { - minY = y->data.F32[i]; - } - } + rc = p_psVectorMin(y, NULL, 0, tmpStatsMinMax); + rc|= p_psVectorMax(y, NULL, 0, tmpStatsMinMax); + if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n"); + psFree(tmpMaskVec); + psFree(robustHistogram); + psFree(cumulativeRobustHistogram); + psFree(tmpStatsMinMax); + psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); + return(1); + } + // // Normalize y to [0.0:1.0] (since the psMinimizeLMChi2Gauss1D() functions is [0.0:1.0]) @@ -1879,5 +1879,5 @@ // for (psS32 i = 0 ; i < y->n ; i++) { - y->data.F32[i]= (y->data.F32[i] - minY) / (maxY - minY); + y->data.F32[i]= (y->data.F32[i] - tmpStatsMinMax->min) / (tmpStatsMinMax->max - tmpStatsMinMax->min); } @@ -2295,5 +2295,5 @@ } } else { - // XXX: This if-statement really shouldn't be necessary. + // This if-statement really shouldn't be necessary. // However, due to numerical lack of precision, we // occasionally produce a binNum outside the range. @@ -2435,9 +2435,10 @@ XXX: Should we free stats if the asserts fail? *****************************************************************************/ -psStats* psVectorStats(psStats* stats, - const psVector* in, - const psVector* errors, - const psVector* mask, - psMaskType maskVal) +psStats* psVectorStats( + psStats* stats, + const psVector* in, + const psVector* errors, + const psVector* mask, + psMaskType maskVal) { psTrace(__func__, 3,"---- %s() begin ----\n", __func__);