Index: trunk/psLib/src/math/Makefile.am =================================================================== --- trunk/psLib/src/math/Makefile.am (revision 4953) +++ trunk/psLib/src/math/Makefile.am (revision 4958) @@ -9,5 +9,6 @@ psMinimize.c \ psRandom.c \ - psFunctions.c \ + psPolynomial.c \ + psSpline.c \ psStats.c \ psUnaryOp.c @@ -23,5 +24,6 @@ psMinimize.h \ psRandom.h \ - psFunctions.h \ + psPolynomial.h \ + psSpline.h \ psStats.h \ psUnaryOp.h Index: trunk/psLib/src/math/psFunctions.c =================================================================== --- trunk/psLib/src/math/psFunctions.c (revision 4953) +++ trunk/psLib/src/math/psFunctions.c (revision 4958) @@ -7,6 +7,6 @@ * polynomials. It also contains a Gaussian functions. * -* @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-08-31 22:28:35 $ +* @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-09-07 21:35:50 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -1307,4 +1307,6 @@ /***************************************************************************** p_psVectorBinDisect(): A wrapper to the above p_psVectorBinDisect(). + +XXX: Assert that the psVector and psScalar have the same type. *****************************************************************************/ psS32 p_psVectorBinDisect(psVector *bins, Index: trunk/psLib/src/math/psMinimize.c =================================================================== --- trunk/psLib/src/math/psMinimize.c (revision 4953) +++ trunk/psLib/src/math/psMinimize.c (revision 4958) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.133 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-09-02 21:32:06 $ + * @version $Revision: 1.134 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-07 21:35:50 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -29,4 +29,5 @@ #include "psImageStructManip.h" #include "psBinaryOp.h" +#include "psLogMsg.h" /*****************************************************************************/ /* DEFINE STATEMENTS */ @@ -112,9 +113,11 @@ /****************************************************************************** -p_psBuildSums1D(x, polyOrder, sums): this routine calculates the powers of -input parameter "x" between 0 and input parameter polyOrder. The result is -returned as a psVector sums. +buildSums1D(x, polyOrder, sums): this routine calculates the powers of input +parameter "x" between 0 and input parameter polyOrder. The result is returned +as a psVector sums. XXX: Use a static vector. + +XXX: should the argument be polyOrder or numTerms? *****************************************************************************/ static void buildSums1D(psF64 x, @@ -124,14 +127,14 @@ psS32 i = 0; psF64 xSum = 0.0; + psS32 numTerms = polyOrder + 1; if (sums == NULL) { - sums = psVectorAlloc(polyOrder, PS_TYPE_F64); - } - if (polyOrder > sums->n) { - sums = psVectorRealloc(sums, polyOrder); + sums = psVectorAlloc(numTerms, PS_TYPE_F64); + } else if (numTerms > sums->n) { + sums = psVectorRealloc(sums, numTerms); } xSum = 1.0; - for (i = 0; i <= polyOrder; i++) { + for (i = 0; i < numTerms; i++) { sums->data.F64[i] = xSum; xSum *= x; @@ -708,5 +711,5 @@ those coords. -// XXX EAM this is my re-implementation of MinLM +XXX EAM this is my re-implementation of MinLM XXX: This must work for both F32 and F64. F32 is currently implemented. @@ -1306,4 +1309,6 @@ XXX: type F32 is done via vector conversion only. + +XXX: Get rid of this. Use new argument list. *****************************************************************************/ psPolynomial1D* psVectorFitPolynomial1D(psPolynomial1D* poly, @@ -2143,8 +2148,13 @@ /****************************************************************************** + ****************************************************************************** + ****************************************************************************** + ****************************************************************************** EAM Code: + ****************************************************************************** + ****************************************************************************** + ****************************************************************************** *****************************************************************************/ -// XXX EAM : my alternate EAMBuildSums1D static psVector *EAMBuildSums1D( psVector* sums, @@ -2155,9 +2165,12 @@ psF64 xSum = 0.0; + // + // XXX: Why do we multiply by 2 here? Better to do it outside and have the + // definition of this function remain sensible. + // nSum = 2*nTerm; if (sums == NULL) { sums = psVectorAlloc(nSum, PS_TYPE_F64); - } - if (nSum > sums->n) { + } else if (nSum > sums->n) { sums = psVectorRealloc(sums, nSum); } @@ -2177,6 +2190,5 @@ const psVector *x, const psVector *y, - const psVector *yErr -) + const psVector *yErr) { // I think this is 1 dimension down @@ -2277,6 +2289,5 @@ const psPolynomial2D *myPoly, const psVector *x, - const psVector *y -) + const psVector *y) { PS_ASSERT_POLY_NON_NULL(myPoly, NULL); @@ -2312,6 +2323,5 @@ psF64 y, psS32 nXterm, - psS32 nYterm -) + psS32 nYterm) { psS32 nXsum = 0; @@ -2349,6 +2359,5 @@ const psVector* y, const psVector* z, - const psVector* zErr -) + const psVector* zErr) { // I think this is 1 dimension down @@ -2439,6 +2448,5 @@ const psVector* y, const psVector* z, - const psVector* dz -) + const psVector* dz) { psVector *X; @@ -2487,4 +2495,15 @@ // XXX EAM : be careful here with F32 vs F64 vectors +/* + Basically, you repetitively fit a polynomial to a set of data points, + reject the points that did not fit well, then refit the polynomial. + + Basically, simply fit the polynomial to the data. They compare the + fit (by evaluating the x data with that polynomial and subtracting + from the original f data). That's the residual. Loop through all + data and if the ((residual - mean) > 3*stDev), mask that data point, + and fit the polynomial again. Do this 3 times. +*/ + psPolynomial2D* RobustFit2D(psPolynomial2D* poly, psVector* mask, @@ -2531,3 +2550,412 @@ } - +/****************************************************************************** + ****************************************************************************** + ****************************************************************************** + ****************************************************************************** +NEW Code: + ****************************************************************************** + ****************************************************************************** + ****************************************************************************** + *****************************************************************************/ + +#define PS_VECTOR_GEN_INDEX_F32(VEC, SIZE) \ +VEC = psVectorAlloc(SIZE, PS_TYPE_F32); \ +for (psS32 i = 0 ; i < SIZE ; i++) { \ + VEC->data.F32[i] = (psF32) i; \ +} \ + +psPolynomial1D *psVectorFitPolynomial1D_NEW( + psPolynomial1D *poly, + const psVector *mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, + const psVector *x) +{ + // Internal pointers for possibly NULL vectors. + psVector *x32 = NULL; + psVector *fErr32 = NULL; + + PS_ASSERT_POLY_NON_NULL(poly, NULL); + PS_ASSERT_VECTOR_NON_NULL(f, NULL); + PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, 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(x, PS_TYPE_F32, NULL); + } else { + PS_VECTOR_GEN_INDEX_F32(x32, f->n); + } + + if (fErr != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); + PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL); + } else { + fErr = psVectorAlloc(f->n, PS_TYPE_F32); + PS_VECTOR_SET_F32(fErr, 1.0); + } + + psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented. Returning NULL.\n"); + if (poly->type == PS_POLYNOMIAL_ORD) { + // XXX: Call EAM code + } else if (poly->type == PS_POLYNOMIAL_CHEB) { + // XXX: Call my code + } else { + printf("XXX: ERROR: incorrect polynomial type.\n"); + } + + // Free psVectors that were created for NULL arguments. + if (x == NULL) { + psFree(x32); + } + if (fErr == NULL) { + psFree(fErr32); + } + + return(NULL); +} + + +psPolynomial2D *psVectorFitPolynomial2D_NEW( + psPolynomial1D *poly, + const psVector *mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, + const psVector *x, + const psVector *y) +{ + // Internal pointers for possibly NULL vectors. + psVector *fErr32 = 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(f, PS_TYPE_F32, 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, y, NULL); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL); + if (fErr != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL); + PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL); + } else { + fErr32 = psVectorAlloc(f->n, PS_TYPE_F32); + PS_VECTOR_SET_F32(fErr32, 1.0); + } + + psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented. Returning NULL.\n"); + if (poly->type == PS_POLYNOMIAL_ORD) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: 2-D polynomial vector fitting has not been implemented. Returning NULL.\n"); + } else if (poly->type == PS_POLYNOMIAL_CHEB) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: 2-D Chebyshev polynomial vector fitting has not been implemented. Returning NULL.\n"); + } else { + printf("XXX: ERROR: incorrect polynomial type.\n"); + } + if (fErr == NULL) { + psFree(fErr32); + } + return(NULL); +} + +psPolynomial3D *psVectorFitPolynomial3D_NEW( + psPolynomial1D *poly, + const psVector *mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, + const psVector *x, + const psVector *y, + const psVector *z) +{ + // Internal pointers for possibly NULL vectors. + psVector *fErr32 = 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(f, PS_TYPE_F32, 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, y, NULL); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_NON_NULL(z, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL); + PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, NULL); + if (fErr != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL); + PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL); + } else { + fErr32 = psVectorAlloc(f->n, PS_TYPE_F32); + PS_VECTOR_SET_F32(fErr32, 1.0); + } + + psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented. Returning NULL.\n"); + if (poly->type == PS_POLYNOMIAL_ORD) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: 3-D polynomial vector fitting has not been implemented. Returning NULL.\n"); + } else if (poly->type == PS_POLYNOMIAL_CHEB) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: 3-D Chebyshev polynomial vector fitting has not been implemented. Returning NULL.\n"); + } else { + printf("XXX: ERROR: incorrect polynomial type.\n"); + } + if (fErr == NULL) { + psFree(fErr32); + } + return(NULL); +} + +psPolynomial4D *psVectorFitPolynomial4D_NEW( + psPolynomial1D *poly, + const psVector *mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, + const psVector *x, + const psVector *y, + const psVector *z, + const psVector *t) +{ + // Internal pointers for possibly NULL vectors. + psVector *fErr32 = 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(f, PS_TYPE_F32, 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, y, NULL); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_NON_NULL(z, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL); + PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_NON_NULL(t, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL); + PS_ASSERT_VECTOR_TYPE(t, PS_TYPE_F32, NULL); + if (fErr != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL); + PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL); + } else { + fErr32 = psVectorAlloc(f->n, PS_TYPE_F32); + PS_VECTOR_SET_F32(fErr32, 1.0); + } + + if (poly->type == PS_POLYNOMIAL_ORD) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: 4-D polynomial vector fitting has not been implemented. Returning NULL.\n"); + } else if (poly->type == PS_POLYNOMIAL_CHEB) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: 4-D Chebyshev polynomial vector fitting has not been implemented. Returning NULL.\n"); + } else { + printf("XXX: ERROR: incorrect polynomial type.\n"); + } + if (fErr == NULL) { + psFree(fErr32); + } + return(NULL); +} + + +psPolynomial1D *psVectorClipFitPolynomial1D_NEW( + psPolynomial1D *poly, + psStats *stats, + const psVector *mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, + const psVector *x) +{ + // Internal pointers for possibly NULL vectors. + psVector *x32 = NULL; + psVector *fErr32 = 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(f, PS_TYPE_F32, 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(x, PS_TYPE_F32, NULL); + } + if (fErr != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); + PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL); + } else { + fErr32 = psVectorAlloc(f->n, PS_TYPE_F32); + PS_VECTOR_SET_F32(fErr32, 1.0); + } + + psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented. Returning NULL.\n"); + // Free psVectors that were created for NULL arguments. + if (x == NULL) { + psFree(x32); + } + if (fErr == NULL) { + psFree(fErr32); + } + return(NULL); +} + +psPolynomial2D *psVectorClipFitPolynomial2D_NEW( + psPolynomial1D *poly, + psStats *stats, + const psVector *mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, + const psVector *x, + const psVector *y) +{ + // Internal pointers for possibly NULL vectors. + psVector *fErr32 = 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(f, PS_TYPE_F32, 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, y, NULL); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL); + if (fErr != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL); + PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL); + } else { + fErr32 = psVectorAlloc(f->n, PS_TYPE_F32); + PS_VECTOR_SET_F32(fErr32, 1.0); + } + + psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented. Returning NULL.\n"); + if (fErr == NULL) { + psFree(fErr32); + } + return(NULL); +} + +psPolynomial3D *psVectorClipFitPolynomial3D_NEW( + psPolynomial1D *poly, + psStats *stats, + const psVector *mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, + const psVector *x, + const psVector *y, + const psVector *z) +{ + // Internal pointers for possibly NULL vectors. + psVector *fErr32 = 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(f, PS_TYPE_F32, 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, y, NULL); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_NON_NULL(z, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL); + PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, NULL); + if (fErr != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL); + PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL); + } else { + fErr32 = psVectorAlloc(f->n, PS_TYPE_F32); + PS_VECTOR_SET_F32(fErr32, 1.0); + } + + psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented. Returning NULL.\n"); + if (fErr == NULL) { + psFree(fErr32); + } + return(NULL); +} + +psPolynomial4D *psVectorClipFitPolynomial4D_NEW( + psPolynomial1D *poly, + psStats *stats, + const psVector *mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, + const psVector *x, + const psVector *y, + const psVector *z, + const psVector *t) +{ + // Internal pointers for possibly NULL vectors. + psVector *fErr32 = 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(f, PS_TYPE_F32, 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, y, NULL); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_NON_NULL(z, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL); + PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_NON_NULL(t, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL); + PS_ASSERT_VECTOR_TYPE(t, PS_TYPE_F32, NULL); + if (fErr != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL); + PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL); + } else { + fErr32 = psVectorAlloc(f->n, PS_TYPE_F32); + PS_VECTOR_SET_F32(fErr32, 1.0); + } + + psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented. Returning NULL.\n"); + if (fErr == NULL) { + psFree(fErr32); + } + return(NULL); +} + + Index: trunk/psLib/src/math/psMinimize.h =================================================================== --- trunk/psLib/src/math/psMinimize.h (revision 4953) +++ trunk/psLib/src/math/psMinimize.h (revision 4958) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.53 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-08-30 01:14:13 $ + * @version $Revision: 1.54 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-07 21:35:50 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -91,4 +91,100 @@ ); + + + + +psPolynomial1D *psVectorFitPolynomial1D_NEW( + psPolynomial1D *poly, + const psVector *mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, + const psVector *x +); + +psPolynomial2D *psVectorFitPolynomial2D_NEW( + psPolynomial1D *poly, + const psVector *mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, + const psVector *x, + const psVector *y +); + +psPolynomial3D *psVectorFitPolynomial3D_NEW( + psPolynomial1D *poly, + const psVector *mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, + const psVector *x, + const psVector *y, + const psVector *z +); + +psPolynomial4D *psVectorFitPolynomial4D_NEW( + psPolynomial1D *poly, + const psVector *mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, + const psVector *x, + const psVector *y, + const psVector *z, + const psVector *t +); + + +psPolynomial1D *psVectorClipFitPolynomial1D_NEW( + psPolynomial1D *poly, + psStats *stats, + const psVector *mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, + const psVector *x +); + +psPolynomial2D *psVectorClipFitPolynomial2D_NEW( + psPolynomial1D *poly, + psStats *stats, + const psVector *mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, + const psVector *x, + const psVector *y +); + +psPolynomial3D *psVectorClipFitPolynomial3D_NEW( + psPolynomial1D *poly, + psStats *stats, + const psVector *mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, + const psVector *x, + const psVector *y, + const psVector *z +); + +psPolynomial4D *psVectorClipFitPolynomial4D_NEW( + psPolynomial1D *poly, + psStats *stats, + const psVector *mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, + const psVector *x, + const psVector *y, + const psVector *z, + const psVector *t +); + + + + /** Derive a one-dimensional spline fit. * @@ -230,4 +326,7 @@ /* \} */// End of MathGroup Functions + + + #endif // #ifndef PS_MINIMIZE_H Index: trunk/psLib/src/math/psPolynomial.c =================================================================== --- trunk/psLib/src/math/psPolynomial.c (revision 4958) +++ trunk/psLib/src/math/psPolynomial.c (revision 4958) @@ -0,0 +1,928 @@ +/** @file psFunctions.c +* +* @brief Contains basic function allocation, deallocation, and evaluation +* routines. +* +* This file will hold the functions for allocated, freeing, and evaluating +* polynomials. It also contains a Gaussian functions. +* +* @version $Revision: 1.116 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-09-07 21:35:50 $ +* +* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii +* +* XXX: Should the "coeffErr[]" be used as well? Bug ???. Ignore coeffErr +* +* XXX: In the various polyAlloc(n) functions, n is really the order of the +* polynomial plus 1. To create a 2nd-order polynomial, n == 3. +*/ +/*****************************************************************************/ +/* INCLUDE FILES */ +/*****************************************************************************/ +#include +#include +#include +#include + +#include "psRandom.h" +#include "psMemory.h" +#include "psVector.h" +#include "psScalar.h" +#include "psTrace.h" +#include "psError.h" +#include "psLogMsg.h" +#include "psPolynomial.h" +#include "psConstants.h" +#include "psErrorText.h" + +/*****************************************************************************/ +/* DEFINE STATEMENTS */ +/*****************************************************************************/ + +/*****************************************************************************/ +/* TYPE DEFINITIONS */ +/*****************************************************************************/ +static void polynomial1DFree(psPolynomial1D* poly); +static void polynomial2DFree(psPolynomial2D* poly); +static void polynomial3DFree(psPolynomial3D* poly); +static void polynomial4DFree(psPolynomial4D* poly); + +/*****************************************************************************/ +/* GLOBAL VARIABLES */ +/*****************************************************************************/ + +// None + +/*****************************************************************************/ +/* FILE STATIC VARIABLES */ +/*****************************************************************************/ + +// None + +/*****************************************************************************/ +/* FUNCTION IMPLEMENTATION - LOCAL */ +/*****************************************************************************/ + +static void polynomial1DFree(psPolynomial1D* poly) +{ + psFree(poly->coeff); + psFree(poly->coeffErr); + psFree(poly->mask); +} + +static void polynomial2DFree(psPolynomial2D* poly) +{ + unsigned int x = 0; + + for (x = 0; x < poly->nX; x++) { + psFree(poly->coeff[x]); + psFree(poly->coeffErr[x]); + psFree(poly->mask[x]); + } + psFree(poly->coeff); + psFree(poly->coeffErr); + psFree(poly->mask); +} + +static void polynomial3DFree(psPolynomial3D* poly) +{ + unsigned int x = 0; + unsigned int y = 0; + + for (x = 0; x < poly->nX; x++) { + for (y = 0; y < poly->nY; y++) { + psFree(poly->coeff[x][y]); + psFree(poly->coeffErr[x][y]); + psFree(poly->mask[x][y]); + } + psFree(poly->coeff[x]); + psFree(poly->coeffErr[x]); + psFree(poly->mask[x]); + } + + psFree(poly->coeff); + psFree(poly->coeffErr); + psFree(poly->mask); +} + +static void polynomial4DFree(psPolynomial4D* poly) +{ + unsigned int x = 0; + unsigned int y = 0; + unsigned int z = 0; + + for (x = 0; x < poly->nX; x++) { + for (y = 0; y < poly->nY; y++) { + for (z = 0; z < poly->nZ; z++) { + psFree(poly->coeff[x][y][z]); + psFree(poly->coeffErr[x][y][z]); + psFree(poly->mask[x][y][z]); + } + psFree(poly->coeff[x][y]); + psFree(poly->coeffErr[x][y]); + psFree(poly->mask[x][y]); + } + psFree(poly->coeff[x]); + psFree(poly->coeffErr[x]); + psFree(poly->mask[x]); + } + + psFree(poly->coeff); + psFree(poly->coeffErr); + psFree(poly->mask); +} + +/***************************************************************************** +createChebyshevPolys(n): this routine takes as input the required order n, +and returns as output as a pointer to an array of n psPolynomial1D +structures, corresponding to the first n Chebyshev polynomials. + +XXX: The output should be static since the Chebyshev polynomials might be +used frequently and the data structure created here does not contain the +outer coefficients of the Chebyshev polynomials. + *****************************************************************************/ +static psPolynomial1D **createChebyshevPolys(psS32 maxChebyPoly) +{ + PS_ASSERT_INT_NONNEGATIVE(maxChebyPoly, NULL); + + psPolynomial1D **chebPolys = NULL; + + chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly * sizeof(psPolynomial1D *)); + for (psS32 i = 0; i < maxChebyPoly; i++) { + chebPolys[i] = psPolynomial1DAlloc(i + 1, PS_POLYNOMIAL_ORD); + } + + // Create the Chebyshev polynomials. + // Polynomial i has i-th order. + chebPolys[0]->coeff[0] = 1; + + // XXX: Bug 296 + if (maxChebyPoly > 1) { + chebPolys[1]->coeff[1] = 1; + + for (psS32 i = 2; i < maxChebyPoly; i++) { + for (psS32 j = 0; j < chebPolys[i - 1]->n; j++) { + chebPolys[i]->coeff[j + 1] = 2 * chebPolys[i - 1]->coeff[j]; + } + for (psS32 j = 0; j < chebPolys[i - 2]->n; j++) { + chebPolys[i]->coeff[j] -= chebPolys[i - 2]->coeff[j]; + } + } + } else { + // XXX: Code this. + printf("WARNING: %d-order chebyshev polynomials not correctly implemented.\n", maxChebyPoly); + } + + return (chebPolys); +} + +/***************************************************************************** + Polynomial coefficients will be accessed in [w][x][y][z] fashion. + *****************************************************************************/ +static psF64 ordPolynomial1DEval(psF64 x, const psPolynomial1D* poly) +{ + psS32 loop_x = 0; + psF64 polySum = 0.0; + psF64 xSum = 1.0; + + psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4, + "---- Calling ordPolynomial1DEval(%lf)\n", x); + psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4, + "Polynomial order is %d\n", poly->n); + for (loop_x = 0; loop_x < poly->n; loop_x++) { + psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4, + "Polynomial coeff[%d] is %lf\n", loop_x, poly->coeff[loop_x]); + } + + for (loop_x = 0; loop_x < poly->n; loop_x++) { + if (poly->mask[loop_x] == 0) { + psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 10, + "polysum+= sum*coeff [%lf+= (%lf * %lf)\n", polySum, xSum, poly->coeff[loop_x]); + polySum += xSum * poly->coeff[loop_x]; + } + xSum *= x; + } + + return(polySum); +} + +// 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; +static psF64 chebPolynomial1DEval(psF64 x, const psPolynomial1D* poly) +{ + PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0); + // XXX: Create a macro for this in psConstants.h + if (poly->n < 1) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Error: Chebyshev polynomial as order %d.", poly->n); + return(NAN); + } + psVector *d; + psS32 n = poly->n; + psS32 i; + psF64 tmp = 0.0; + + // Special case where the Chebyshev poly is constant. + if (n == 1) { + if (poly->mask[0] == 0) { + tmp += poly->coeff[0]; + } + return(tmp); + } + + // Special case where the Chebyshev poly is linear. + if (n == 2) { + if (poly->mask[0] == 0) { + tmp+= poly->coeff[0]; + } + if (poly->mask[1] == 0) { + tmp+= poly->coeff[1] * x; + } + return(tmp); + } + + // General case where the Chebyshev poly has 2 or more terms. + d = psVectorAlloc(n, PS_TYPE_F64); + if(poly->mask[n-1] == 0) { + d->data.F64[n-1] = poly->coeff[n-1]; + } else { + d->data.F64[n-1] = 0.0; + } + + d->data.F64[n-2] = (2.0 * x * d->data.F64[n-1]); + if(poly->mask[n-2] == 0) { + d->data.F64[n-2] += poly->coeff[n-2]; + } + + for (i=n-3;i>=1;i--) { + d->data.F64[i] = (2.0 * x * d->data.F64[i+1]) - + (d->data.F64[i+2]); + if(poly->mask[i] == 0) { + d->data.F64[i] += poly->coeff[i]; + } + } + + tmp = (x * d->data.F64[1]) - + (d->data.F64[2]); + if(poly->mask[0] == 0) { + tmp += (0.5 * poly->coeff[0]); + } + psFree(d); + return(tmp); + + /* This is old code that does not use Clenshaw's formula. Get rid of it. + + psS32 n; + psS32 i; + psF32 tmp; + psPolynomial1D **chebPolys = NULL; + + n = poly->n; + chebPolys = createChebyshevPolys(n); + + tmp = 0.0; + for (i=0;in;i++) { + tmp+= (poly->coeff[i] * psPolynomial1DEval(x, chebPolys[i])); + } + tmp-= (poly->coeff[0]/2.0); + + + return(tmp); + */ +} + +static psF64 ordPolynomial2DEval(psF64 x, + psF64 y, + const psPolynomial2D* poly) +{ + PS_ASSERT_POLY_NON_NULL(poly, NAN); + + psS32 loop_x = 0; + psS32 loop_y = 0; + psF64 polySum = 0.0; + psF64 xSum = 1.0; + psF64 ySum = 1.0; + + for (loop_x = 0; loop_x < poly->nX; loop_x++) { + ySum = xSum; + for (loop_y = 0; loop_y < poly->nY; loop_y++) { + if (poly->mask[loop_x][loop_y] == 0) { + polySum += ySum * poly->coeff[loop_x][loop_y]; + } + ySum *= y; + } + xSum *= x; + } + + return(polySum); +} + +static psF64 chebPolynomial2DEval(psF64 x, psF64 y, const psPolynomial2D* poly) +{ + PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0); + PS_ASSERT_DOUBLE_WITHIN_RANGE(y, -1.0, 1.0, 0.0); + PS_ASSERT_POLY_NON_NULL(poly, NAN); + + psS32 loop_x = 0; + psS32 loop_y = 0; + psS32 i = 0; + psF64 polySum = 0.0; + psPolynomial1D* *chebPolys = NULL; + psS32 maxChebyPoly = 0; + + // Determine how many Chebyshev polynomials + // are needed, then create them. + maxChebyPoly = poly->nX; + if (poly->nY > maxChebyPoly) { + maxChebyPoly = poly->nY; + } + chebPolys = createChebyshevPolys(maxChebyPoly); + + for (loop_x = 0; loop_x < poly->nX; loop_x++) { + for (loop_y = 0; loop_y < poly->nY; loop_y++) { + if (poly->mask[loop_x][loop_y] == 0) { + polySum += poly->coeff[loop_x][loop_y] * + psPolynomial1DEval(chebPolys[loop_x], x) * + psPolynomial1DEval(chebPolys[loop_y], y); + } + } + } + for (i=0;inX; loop_x++) { + ySum = xSum; + for (loop_y = 0; loop_y < poly->nY; loop_y++) { + zSum = ySum; + for (loop_z = 0; loop_z < poly->nZ; loop_z++) { + if (poly->mask[loop_x][loop_y][loop_z] == 0) { + polySum += zSum * poly->coeff[loop_x][loop_y][loop_z]; + } + zSum *= z; + } + ySum *= y; + } + xSum *= x; + } + + return(polySum); +} + +static psF64 chebPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psPolynomial3D* poly) +{ + PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0); + PS_ASSERT_DOUBLE_WITHIN_RANGE(y, -1.0, 1.0, 0.0); + PS_ASSERT_DOUBLE_WITHIN_RANGE(z, -1.0, 1.0, 0.0); + psS32 loop_x = 0; + psS32 loop_y = 0; + psS32 loop_z = 0; + psS32 i = 0; + psF64 polySum = 0.0; + psPolynomial1D* *chebPolys = NULL; + psS32 maxChebyPoly = 0; + + // Determine how many Chebyshev polynomials + // are needed, then create them. + maxChebyPoly = poly->nX; + if (poly->nY > maxChebyPoly) { + maxChebyPoly = poly->nY; + } + if (poly->nZ > maxChebyPoly) { + maxChebyPoly = poly->nZ; + } + chebPolys = createChebyshevPolys(maxChebyPoly); + + for (loop_x = 0; loop_x < poly->nX; loop_x++) { + for (loop_y = 0; loop_y < poly->nY; loop_y++) { + for (loop_z = 0; loop_z < poly->nZ; loop_z++) { + if (poly->mask[loop_x][loop_y][loop_z] == 0) { + polySum += poly->coeff[loop_x][loop_y][loop_z] * + psPolynomial1DEval(chebPolys[loop_x], x) * + psPolynomial1DEval(chebPolys[loop_y], y) * + psPolynomial1DEval(chebPolys[loop_z], z); + } + } + } + } + + for (i=0;inX; loop_x++) { + ySum = xSum; + for (loop_y = 0; loop_y < poly->nY; loop_y++) { + zSum = ySum; + for (loop_z = 0; loop_z < poly->nZ; loop_z++) { + tSum = zSum; + for (loop_t = 0; loop_t < poly->nT; loop_t++) { + if (poly->mask[loop_x][loop_y][loop_z][loop_t] == 0) { + polySum += tSum * poly->coeff[loop_x][loop_y][loop_z][loop_t]; + } + tSum *= t; + } + zSum *= z; + } + ySum *= y; + } + xSum *= x; + } + + return(polySum); +} + +static psF64 chebPolynomial4DEval(psF64 x, psF64 y, psF64 z, psF64 t, const psPolynomial4D* poly) +{ + PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0); + PS_ASSERT_DOUBLE_WITHIN_RANGE(y, -1.0, 1.0, 0.0); + PS_ASSERT_DOUBLE_WITHIN_RANGE(z, -1.0, 1.0, 0.0); + PS_ASSERT_DOUBLE_WITHIN_RANGE(t, -1.0, 1.0, 0.0); + psS32 loop_x = 0; + psS32 loop_y = 0; + psS32 loop_z = 0; + psS32 loop_t = 0; + psS32 i = 0; + psF64 polySum = 0.0; + psPolynomial1D* *chebPolys = NULL; + psS32 maxChebyPoly = 0; + + // Determine how many Chebyshev polynomials + // are needed, then create them. + maxChebyPoly = poly->nX; + if (poly->nY > maxChebyPoly) { + maxChebyPoly = poly->nY; + } + if (poly->nZ > maxChebyPoly) { + maxChebyPoly = poly->nZ; + } + if (poly->nT > maxChebyPoly) { + maxChebyPoly = poly->nT; + } + chebPolys = createChebyshevPolys(maxChebyPoly); + + for (loop_x = 0; loop_x < poly->nX; loop_x++) { + for (loop_y = 0; loop_y < poly->nY; loop_y++) { + for (loop_z = 0; loop_z < poly->nZ; loop_z++) { + for (loop_t = 0; loop_t < poly->nT; loop_t++) { + if (poly->mask[loop_x][loop_y][loop_z][loop_t] == 0) { + polySum += poly->coeff[loop_x][loop_y][loop_z][loop_t] * + psPolynomial1DEval(chebPolys[loop_x], x) * + psPolynomial1DEval(chebPolys[loop_y], y) * + psPolynomial1DEval(chebPolys[loop_z], z) * + psPolynomial1DEval(chebPolys[loop_t], t); + } + } + } + } + } + + for (i=0;idata.F32[i] = mean + p_psRandomGaussian(r, sigma); + } + psFree(r); + + return(gauss); +} + +/***************************************************************************** + This routine must allocate memory for the polynomial structures. + *****************************************************************************/ +psPolynomial1D* psPolynomial1DAlloc(int n, + psPolynomialType type) +{ + PS_ASSERT_INT_POSITIVE(n, NULL); + + int i = 0; + psPolynomial1D* newPoly = NULL; + + newPoly = (psPolynomial1D* ) psAlloc(sizeof(psPolynomial1D)); + psMemSetDeallocator(newPoly, (psFreeFunc) polynomial1DFree); + + newPoly->type = type; + newPoly->n = n; + newPoly->coeff = psAlloc(n * sizeof(psF64)); + newPoly->coeffErr = psAlloc(n * sizeof(psF64)); + newPoly->mask = (char *)psAlloc(n * sizeof(char)); + for (i = 0; i < n; i++) { + newPoly->coeff[i] = 0.0; + newPoly->coeffErr[i] = 0.0; + newPoly->mask[i] = 0; + } + + return(newPoly); +} + +psPolynomial2D* psPolynomial2DAlloc( int nX, int nY, + psPolynomialType type) +{ + PS_ASSERT_INT_POSITIVE(nX, NULL); + PS_ASSERT_INT_POSITIVE(nY, NULL); + + int x = 0; + int y = 0; + psPolynomial2D* newPoly = NULL; + + newPoly = (psPolynomial2D* ) psAlloc(sizeof(psPolynomial2D)); + psMemSetDeallocator(newPoly, (psFreeFunc) polynomial2DFree); + + newPoly->type = type; + newPoly->nX = nX; + newPoly->nY = nY; + + newPoly->coeff = psAlloc(nX * sizeof(psF64 *)); + newPoly->coeffErr = psAlloc(nX * sizeof(psF64 *)); + newPoly->mask = (char **)psAlloc(nX * sizeof(char *)); + for (x = 0; x < nX; x++) { + newPoly->coeff[x] = psAlloc(nY * sizeof(psF64)); + newPoly->coeffErr[x] = psAlloc(nY * sizeof(psF64)); + newPoly->mask[x] = (char *)psAlloc(nY * sizeof(char)); + } + for (x = 0; x < nX; x++) { + for (y = 0; y < nY; y++) { + newPoly->coeff[x][y] = 0.0; + newPoly->coeffErr[x][y] = 0.0; + newPoly->mask[x][y] = 0; + } + } + + return(newPoly); +} + +psPolynomial3D* psPolynomial3DAlloc( int nX, int nY, int nZ, + psPolynomialType type) +{ + PS_ASSERT_INT_POSITIVE(nX, NULL); + PS_ASSERT_INT_POSITIVE(nY, NULL); + PS_ASSERT_INT_POSITIVE(nZ, NULL); + + psS32 x = 0; + psS32 y = 0; + psS32 z = 0; + psPolynomial3D* newPoly = NULL; + + newPoly = (psPolynomial3D* ) psAlloc(sizeof(psPolynomial3D)); + psMemSetDeallocator(newPoly, (psFreeFunc) polynomial3DFree); + + newPoly->type = type; + newPoly->nX = nX; + newPoly->nY = nY; + newPoly->nZ = nZ; + + newPoly->coeff = psAlloc(nX * sizeof(psF64 **)); + newPoly->coeffErr = psAlloc(nX * sizeof(psF64 **)); + newPoly->mask = (char ***)psAlloc(nX * sizeof(char **)); + for (x = 0; x < nX; x++) { + newPoly->coeff[x] = psAlloc(nY * sizeof(psF64 *)); + newPoly->coeffErr[x] = psAlloc(nY * sizeof(psF64 *)); + newPoly->mask[x] = (char **)psAlloc(nY * sizeof(char *)); + for (y = 0; y < nY; y++) { + newPoly->coeff[x][y] = psAlloc(nZ * sizeof(psF64)); + newPoly->coeffErr[x][y] = psAlloc(nZ * sizeof(psF64)); + newPoly->mask[x][y] = (char *)psAlloc(nZ * sizeof(char)); + } + } + for (x = 0; x < nX; x++) { + for (y = 0; y < nY; y++) { + for (z = 0; z < nZ; z++) { + newPoly->coeff[x][y][z] = 0.0; + newPoly->coeffErr[x][y][z] = 0.0; + newPoly->mask[x][y][z] = 0; + } + } + } + + return(newPoly); +} + +psPolynomial4D* psPolynomial4DAlloc( int nX, int nY, int nZ, int nT, + psPolynomialType type) +{ + PS_ASSERT_INT_POSITIVE(nX, NULL); + PS_ASSERT_INT_POSITIVE(nY, NULL); + PS_ASSERT_INT_POSITIVE(nZ, NULL); + PS_ASSERT_INT_POSITIVE(nT, NULL); + + psS32 x = 0; + psS32 y = 0; + psS32 z = 0; + psS32 t = 0; + psPolynomial4D* newPoly = NULL; + + newPoly = (psPolynomial4D* ) psAlloc(sizeof(psPolynomial4D)); + psMemSetDeallocator(newPoly, (psFreeFunc) polynomial4DFree); + + newPoly->type = type; + newPoly->nX = nX; + newPoly->nY = nY; + newPoly->nZ = nZ; + newPoly->nT = nT; + + newPoly->coeff = psAlloc(nX * sizeof(psF64 ***)); + newPoly->coeffErr = psAlloc(nX * sizeof(psF64 ***)); + newPoly->mask = (char ****)psAlloc(nX * sizeof(char ***)); + for (x = 0; x < nX; x++) { + newPoly->coeff[x] = psAlloc(nY * sizeof(psF64 **)); + newPoly->coeffErr[x] = psAlloc(nY * sizeof(psF64 **)); + newPoly->mask[x] = (char ***)psAlloc(nY * sizeof(char **)); + for (y = 0; y < nY; y++) { + newPoly->coeff[x][y] = psAlloc(nZ * sizeof(psF64 *)); + newPoly->coeffErr[x][y] = psAlloc(nZ * sizeof(psF64 *)); + newPoly->mask[x][y] = (char **)psAlloc(nZ * sizeof(char *)); + for (z = 0; z < nZ; z++) { + newPoly->coeff[x][y][z] = psAlloc(nT * sizeof(psF64)); + newPoly->coeffErr[x][y][z] = psAlloc(nT * sizeof(psF64)); + newPoly->mask[x][y][z] = (char *)psAlloc(nT * sizeof(char)); + } + } + } + for (x = 0; x < nX; x++) { + for (y = 0; y < nY; y++) { + for (z = 0; z < nZ; z++) { + for (t = 0; t < nT; t++) { + newPoly->coeff[x][y][z][t] = 0.0; + newPoly->coeffErr[x][y][z][t] = 0.0; + newPoly->mask[x][y][z][t] = 0; + } + } + } + } + + return(newPoly); +} + +psF64 psPolynomial1DEval(const psPolynomial1D* poly, psF64 x) +{ + PS_ASSERT_POLY_NON_NULL(poly, NAN); + + if (poly->type == PS_POLYNOMIAL_ORD) { + return(ordPolynomial1DEval(x, poly)); + } else if (poly->type == PS_POLYNOMIAL_CHEB) { + return(chebPolynomial1DEval(x, poly)); + } else { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, + poly->type); + } + return(NAN); +} + +psVector *psPolynomial1DEvalVector(const psPolynomial1D *poly, + const psVector *x) +{ + PS_ASSERT_POLY_NON_NULL(poly, NULL); + PS_ASSERT_VECTOR_NON_NULL(x, NULL); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); + + psVector *tmp; + + tmp = psVectorAlloc(x->n, PS_TYPE_F64); + for (psS32 i=0;in;i++) { + tmp->data.F64[i] = psPolynomial1DEval(poly, x->data.F64[i]); + } + + return(tmp); +} + +psF64 psPolynomial2DEval(const psPolynomial2D* poly, psF64 x, psF64 y) +{ + PS_ASSERT_POLY_NON_NULL(poly, NAN); + + if (poly->type == PS_POLYNOMIAL_ORD) { + return(ordPolynomial2DEval(x, y, poly)); + } else if (poly->type == PS_POLYNOMIAL_CHEB) { + return(chebPolynomial2DEval(x, y, poly)); + } else { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, + poly->type); + } + return(NAN); +} + +psVector *psPolynomial2DEvalVector(const psPolynomial2D *poly, + const psVector *x, + const psVector *y) + +{ + PS_ASSERT_POLY_NON_NULL(poly, NULL); + PS_ASSERT_VECTOR_NON_NULL(x, NULL); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL); + + psVector *tmp; + psS32 vecLen=x->n; + + // Determine the length of the output vector to by the minimum of the x,y vectors + if (y->n < vecLen) { + vecLen = y->n; + } + + // Create output vector to return + tmp = psVectorAlloc(vecLen, PS_TYPE_F64); + + // Evaluate the polynomial at the specified points + for (psS32 i=0; idata.F64[i] = psPolynomial2DEval(poly,x->data.F64[i],y->data.F64[i]); + } + + // Return output vector + return(tmp); +} + +psF64 psPolynomial3DEval(const psPolynomial3D* poly, psF64 x, psF64 y, psF64 z) +{ + PS_ASSERT_POLY_NON_NULL(poly, NAN); + + if (poly->type == PS_POLYNOMIAL_ORD) { + return(ordPolynomial3DEval(x, y, z, poly)); + } else if (poly->type == PS_POLYNOMIAL_CHEB) { + return(chebPolynomial3DEval(x, y, z, poly)); + } else { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, + poly->type); + } + return(NAN); +} + +psVector *psPolynomial3DEvalVector(const psPolynomial3D *poly, + const psVector *x, + const psVector *y, + const psVector *z) + +{ + PS_ASSERT_POLY_NON_NULL(poly, NULL); + PS_ASSERT_VECTOR_NON_NULL(x, NULL); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL); + PS_ASSERT_VECTOR_NON_NULL(z, NULL); + PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F64, NULL); + + psVector *tmp; + psS32 vecLen=x->n; + + // Determine the length of output vector from min of the input vectors + if (y->n < vecLen) { + vecLen = y->n; + } + if (z->n < vecLen) { + vecLen = z->n; + } + + // Allocate output vector + tmp = psVectorAlloc(vecLen, PS_TYPE_F64); + + // Evaluate polynomial + for (psS32 i = 0; i < vecLen; i++) { + tmp->data.F64[i] = psPolynomial3DEval(poly, + x->data.F64[i], + y->data.F64[i], + z->data.F64[i]); + } + + // Return output vector + return(tmp); +} + +psF64 psPolynomial4DEval(const psPolynomial4D* poly, psF64 x, psF64 y, psF64 z, psF64 t) +{ + PS_ASSERT_POLY_NON_NULL(poly, NAN); + + if (poly->type == PS_POLYNOMIAL_ORD) { + return(ordPolynomial4DEval(x,y,z,t, poly)); + } else if (poly->type == PS_POLYNOMIAL_CHEB) { + return(chebPolynomial4DEval(x,y,z,t, poly)); + } else { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, + poly->type); + } + return(NAN); +} + +psVector *psPolynomial4DEvalVector(const psPolynomial4D *poly, + const psVector *x, + const psVector *y, + const psVector *z, + const psVector *t) +{ + PS_ASSERT_POLY_NON_NULL(poly, NULL); + PS_ASSERT_VECTOR_NON_NULL(x, NULL); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL); + PS_ASSERT_VECTOR_NON_NULL(z, NULL); + PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F64, NULL); + PS_ASSERT_VECTOR_NON_NULL(t, NULL); + PS_ASSERT_VECTOR_TYPE(t, PS_TYPE_F64, NULL); + + psVector *tmp; + psS32 vecLen=x->n; + + // Determine output vector size from min of input vectors + if (z->n < vecLen) { + vecLen = z->n; + } + if (y->n < vecLen) { + vecLen = y->n; + } + if (t->n < vecLen) { + vecLen = t->n; + } + + // Allocate output vector + tmp = psVectorAlloc(vecLen, PS_TYPE_F64); + + // Evaluate polynomial + for (psS32 i = 0; i < vecLen; i++) { + tmp->data.F64[i] = psPolynomial4DEval(poly, + x->data.F64[i], + y->data.F64[i], + z->data.F64[i], + t->data.F64[i]); + } + + // Return output vector + return(tmp); +} + Index: trunk/psLib/src/math/psPolynomial.h =================================================================== --- trunk/psLib/src/math/psPolynomial.h (revision 4958) +++ trunk/psLib/src/math/psPolynomial.h (revision 4958) @@ -0,0 +1,252 @@ +/** @file psPolynomial.h + * @brief Standard Mathematical Functions. + * @ingroup Stats + * + * This file will hold the prototypes for procedures which allocate, free, + * and evaluate various polynomials. Those polynomial structures are also + * defined here. + * + * @ingroup Stats + * + * @author GLG, MHPCC + * + * @version $Revision: 1.53 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-07 21:35:50 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ + +#ifndef PS_POLYNOMIAL_H +#define PS_POLYNOMIAL_H + +#include +#include +#include +#include + +#include "psVector.h" +#include "psScalar.h" +#include "psSpline.h" + +/** \addtogroup Stats + * \{ + */ + +/** Evaluate a non-normalized Gaussian with the given mean and sigma at the + * given coordianate. + * + * Note that this is not a Gaussian deviate. The evaluated Gaussian is: + * \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f] + * + * @return float value on the gaussian curve given the input parameters + */ +float psGaussian( + float x, ///< Value at which to evaluate + float mean, ///< Mean for the Gaussian + float sigma, ///< Standard deviation for the Gaussian + bool normal ///< Indicates whether result should be normalized +); + +/** Produce a vector of random numbers from a Gaussian distribution with + * the specified mean and sigma + * + * @return psVector* vector of random numbers + * + */ +psVector* p_psGaussianDev( + psF32 mean, ///< The mean of the Gaussian + psF32 sigma, ///< The sigma of the Gaussian + psS32 Npts ///< The size of the vector +); + +/** Polynomial Type. + * + * Enumeration for Polynomial types. + */ +typedef enum { + PS_POLYNOMIAL_ORD, ///< Ordinary Polynomial + PS_POLYNOMIAL_CHEB ///< Chebyshev Polynomial +} +psPolynomialType; + +/** One-dimensional polynomial */ +typedef struct +{ + psPolynomialType type; ///< Polynomial type + int n; ///< Number of terms + psF64 *coeff; ///< Coefficients + psF64 *coeffErr; ///< Error in coefficients + char *mask; ///< Coefficient mask +} +psPolynomial1D; + +/** Two-dimensional polynomial */ +typedef struct +{ + psPolynomialType type; ///< Polynomial type + int nX; ///< Number of terms in x + int nY; ///< Number of terms in y + psF64 **coeff; ///< Coefficients + psF64 **coeffErr; ///< Error in coefficients + char **mask; ///< Coefficients mask +} +psPolynomial2D; + +/** Three-dimensional polynomial */ +typedef struct +{ + psPolynomialType type; ///< Polynomial type + int nX; ///< Number of terms in x + int nY; ///< Number of terms in y + int nZ; ///< Number of terms in z + psF64 ***coeff; ///< Coefficients + psF64 ***coeffErr; ///< Error in coefficients + char ***mask; ///< Coefficients mask +} +psPolynomial3D; + +/** Four-dimensional polynomial */ +typedef struct +{ + psPolynomialType type; ///< Polynomial type + int nX; ///< Number of terms in x + int nY; ///< Number of terms in y + int nZ; ///< Number of terms in z + int nT; ///< Number of terms in t + psF64 ****coeff; ///< Coefficients + psF64 ****coeffErr; ///< Error in coefficients + char ****mask; ///< Coefficients mask +} +psPolynomial4D; + + +/** Allocates a psPolynomial1D structure with n terms + * + * @return psPolynomial1D* new 1-D polynomial struct + */ +psPolynomial1D* psPolynomial1DAlloc( + int n, ///< Number of terms + psPolynomialType type ///< Polynomial Type +); + +/** Allocates a 2-D polynomial structure + * + * @return psPolynomial2D* new 2-D polynomial struct + */ +psPolynomial2D* psPolynomial2DAlloc( + int nX, ///< Number of terms in x + int nY, ///< Number of terms in y + psPolynomialType type ///< Polynomial Type +); + +/** Allocates a 3-D polynomial structure + * + * @return psPolynomial3D* new 3-D polynomial struct + */ +psPolynomial3D* psPolynomial3DAlloc( + int nX, ///< Number of terms in x + int nY, ///< Number of terms in y + int nZ, ///< Number of terms in z + psPolynomialType type ///< Polynomial Type +); + +/** Allocates a 4-D polynomial structure + * + * @return psPolynomial4D* new 4-D polynomial struct + */ +psPolynomial4D* psPolynomial4DAlloc( + int nX, ///< Number of terms in x + int nY, ///< Number of terms in y + int nZ, ///< Number of terms in z + int nT, ///< Number of terms in t + psPolynomialType type ///< Polynomial Type +); + +/** Evaluates a 1-D polynomial at specific coordinates. + * + * @return psF64 result of polynomial at given location + */ +psF64 psPolynomial1DEval( + const psPolynomial1D* poly, ///< Coefficients for the polynomial + psF64 x ///< location at which to evaluate +); + +/** Evaluates a 2-D polynomial at specific coordinates. + * + * @return psF64 result of polynomial at given location + */ +psF64 psPolynomial2DEval( + const psPolynomial2D* poly, ///< Coefficients for the polynomial + psF64 x, ///< x location at which to evaluate + psF64 y ///< y location at which to evaluate +); + +/** Evaluates a 3-D polynomial at specific coordinates. + * + * @return psF64 result of polynomial at given location + */ +psF64 psPolynomial3DEval( + const psPolynomial3D* poly, ///< Coefficients for the polynomial + psF64 x, ///< x location at which to evaluate + psF64 y, ///< y location at which to evaluate + psF64 z ///< z location at which to evaluate +); + +/** Evaluates a 4-D polynomial at specific coordinates. + * + * @return psF64 result of polynomial at given location + */ +psF64 psPolynomial4DEval( + const psPolynomial4D* poly, ///< Coefficients for the polynomial + psF64 x, ///< x location at which to evaluate + psF64 y, ///< y location at which to evaluate + psF64 z, ///< z location at which to evaluate + psF64 t ///< t location at which to evaluate +); + +/** Evaluates a 1-D polynomial at specific sets of coordinates + * + * @return psVector* results of polynomials at given locations + */ +psVector *psPolynomial1DEvalVector( + const psPolynomial1D *poly, ///< Coefficients for the polynomial + const psVector *x ///< x locations at which to evaluate +); + +/** Evaluates a 2-D polynomial at specific sets of coordinates + * + * @return psVector* results of polynomial at given locations + */ +psVector *psPolynomial2DEvalVector( + const psPolynomial2D *poly, ///< Coefficients for the polynomial + const psVector *x, ///< x locations at which to evaluate + const psVector *y ///< y locations at which to evaluate +); + +/** Evaluates a 3-D polynomial at specific sets of coordinates + * + * @return psVector* results of polynomial at given locations + */ +psVector *psPolynomial3DEvalVector( + const psPolynomial3D *poly, ///< Coefficients for the polynomial + const psVector *x, ///< x locations at which to evaluate + const psVector *y, ///< y locations at which to evaluate + const psVector *z ///< z locations at which to evaluate +); + +/** Evaluates a 4-D polynomial at specific sets of coordinates + * + * @return psVector* results of polynomial at given locations + */ +psVector *psPolynomial4DEvalVector( + const psPolynomial4D *poly, ///< Coefficients for the polynomial + const psVector *x, ///< x locations at which to evaluate + const psVector *y, ///< y locations at which to evaluate + const psVector *z, ///< z locations at which to evaluate + const psVector *t ///< t locations at which to evaluate +); + +/** \} */ // End of MathGroup Functions + +#endif // #ifndef PS_POLYNOMIAL_H + Index: trunk/psLib/src/math/psSpline.c =================================================================== --- trunk/psLib/src/math/psSpline.c (revision 4958) +++ trunk/psLib/src/math/psSpline.c (revision 4958) @@ -0,0 +1,622 @@ +/** @file psFunctions.c +* +* @brief Contains basic function allocation, deallocation, and evaluation +* routines. +* +* This file will hold the functions for allocated, freeing, and evaluating +* splines. +* +* @version $Revision: 1.116 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-09-07 21:35:50 $ +* +* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii +* +*/ +/*****************************************************************************/ +/* INCLUDE FILES */ +/*****************************************************************************/ +#include +#include +#include +#include + +#include "psRandom.h" +#include "psMemory.h" +#include "psVector.h" +#include "psScalar.h" +#include "psTrace.h" +#include "psError.h" +#include "psLogMsg.h" +#include "psPolynomial.h" +#include "psConstants.h" +#include "psErrorText.h" + +/*****************************************************************************/ +/* DEFINE STATEMENTS */ +/*****************************************************************************/ + +/*****************************************************************************/ +/* TYPE DEFINITIONS */ +/*****************************************************************************/ +static void spline1DFree(psSpline1D *tmpSpline); +static psS32 vectorBinDisectF32(psF32 *bins,psS32 numBins,psF32 x); +static psS32 vectorBinDisectS32(psS32 *bins,psS32 numBins,psS32 x); + +/*****************************************************************************/ +/* GLOBAL VARIABLES */ +/*****************************************************************************/ + +// None + +/*****************************************************************************/ +/* FILE STATIC VARIABLES */ +/*****************************************************************************/ + +// None + +/*****************************************************************************/ +/* FUNCTION IMPLEMENTATION - LOCAL */ +/*****************************************************************************/ + +static void spline1DFree(psSpline1D *tmpSpline) +{ + psS32 i; + + if (tmpSpline == NULL) { + return; + } + + if (tmpSpline->spline != NULL) { + for (i=0;in;i++) { + psFree((tmpSpline->spline)[i]); + } + psFree(tmpSpline->spline); + } + + if (tmpSpline->p_psDeriv2 != NULL) { + psFree(tmpSpline->p_psDeriv2); + } + psFree(tmpSpline->knots); + + return; +} + +/***************************************************************************** +fullInterpolate1DF32(): This routine will take as input n-element floating +point arrays domain and range, and the x value, assumed to lie with the +domain vector. It produces as output the (n-1)-order LaGrange interpolated +value of x. + +XXX: do we error check for non-distinct domain values? + *****************************************************************************/ +#define FUNC_MACRO_FULL_INTERPOLATE_1D(TYPE) \ +static psF32 fullInterpolate1D##TYPE(ps##TYPE *domain, \ + ps##TYPE *range, \ + psS32 n, \ + ps##TYPE x) \ +{ \ + \ + psS32 i; \ + psS32 m; \ + static psVector *p = NULL; \ + p = psVectorRecycle(p, n, PS_TYPE_##TYPE); \ + p_psMemSetPersistent(p, true); \ + p_psMemSetPersistent(p->data.TYPE, true); \ + \ + psTrace(".psLib.dataManip.psFunctions.fullInterpolate1D##TYPE", 4, \ + "---- fullInterpolate1D##TYPE() begin (%d-order at x=%f) (%d data points)----\n", n-1, x, n); \ + \ + for (i=0;idata.TYPE[i] = range[i]; \ + psTrace(".psLib.dataManip.psFunctions.fullInterpolate1D##TYPE", 6, \ + "p->data.TYPE[%d] is %f\n", i, p->data.TYPE[i]); \ + \ + } \ + \ + /* From NR, during each iteration of the m loop, we are computing the \ + p_{i ... i+m} terms. \ + */ \ + for (m=1;mdata.TYPE[i] = (((x-domain[i+m]) * p->data.TYPE[i]) + \ + ((domain[i]-x) * p->data.TYPE[i+1])) / \ + (domain[i] - domain[i+m]); \ + /*printf("((%f-%f * %f) + (%f-%f * %f)) / (%f - %f)\n", x, domain[i+m], p->data.TYPE[i], domain[i], x, p->data.TYPE[i+1], domain[i], domain[i+m]); \ + */ \ + psTrace(".psLib.dataManip.psFunctions.fullInterpolate1D##TYPE", 6, \ + "p->data.TYPE[%d] is %f\n", i, p->data.TYPE[i]); \ + } \ + } \ + psTrace(".psLib.dataManip.psFunctions.fullInterpolate1D##TYPE", 4, \ + "---- fullInterpolate1D##TYPE() end ----\n"); \ + \ + return(p->data.TYPE[0]); \ +} \ + +/* +FUNC_MACRO_FULL_INTERPOLATE_1D(U8) +FUNC_MACRO_FULL_INTERPOLATE_1D(U16) +FUNC_MACRO_FULL_INTERPOLATE_1D(U32) +FUNC_MACRO_FULL_INTERPOLATE_1D(U64) +FUNC_MACRO_FULL_INTERPOLATE_1D(S8) +FUNC_MACRO_FULL_INTERPOLATE_1D(S16) +FUNC_MACRO_FULL_INTERPOLATE_1D(S32) +FUNC_MACRO_FULL_INTERPOLATE_1D(S64) +FUNC_MACRO_FULL_INTERPOLATE_1D(F64) +*/ +FUNC_MACRO_FULL_INTERPOLATE_1D(F32) + + +/***************************************************************************** +interpolate1DF32(): this is the base 1-D flat memory routine to perform +LaGrange interpolation. + *****************************************************************************/ +static psF32 interpolate1DF32(psF32 *domain, + psF32 *range, + psS32 n, + psS32 order, + psF32 x) +{ + PS_ASSERT_PTR_NON_NULL(domain, NAN) + PS_ASSERT_PTR_NON_NULL(range, NAN) + // XXX: Check valid values for n, order, and x? + + psS32 binNum; + psS32 numIntPoints = order+1; + psS32 origin; + + psTrace(".psLib.dataManip.psFunctions.interpolate1DF32", 4, + "---- interpolate1DF32() begin ----\n"); + + binNum = vectorBinDisectF32(domain, n, x); + + if (0 == numIntPoints%2) { + origin = binNum - ((numIntPoints/2) - 1); + } else { + origin = binNum - (numIntPoints/2); + if ((x-domain[binNum]) > (domain[binNum+1]-x)) { + // x is closer to binNum+1. + origin = 1 + (binNum - (numIntPoints/2)); + } + } + if (origin < 0) { + origin = 0; + } + if ((origin + numIntPoints) > n) { + origin = n - numIntPoints; + } + + psTrace(".psLib.dataManip.psFunctions.interpolate1DF32", 4, + "---- interpolate1DF32() end ----\n"); + return(fullInterpolate1DF32(&domain[origin], &range[origin], order+1, x)); +} + +/*****************************************************************************/ +/* FUNCTION IMPLEMENTATION - PUBLIC */ +/*****************************************************************************/ + +//typedef struct { +// psS32 n; +// psPolynomial1D **spline; +// psF32 *p_psDeriv2; +// psVector *knots; +//} psSpline1D; + +/***************************************************************************** + NOTE: "n" specifies the number of spline polynomials. Therefore, there + must exist n+1 points in "knots". + +XXX: Is this really needed anymore? + +XXX: Ensure that knots[i+1] != knots[i] + +XXX: What should be the default type for knots be? psF32 is assumed. + *****************************************************************************/ +psSpline1D *psSpline1DAlloc(int numSplines, + int order, + float min, + float max) +{ + PS_ASSERT_INT_NONNEGATIVE(numSplines, NULL); + PS_ASSERT_INT_NONNEGATIVE(order, NULL); + PS_ASSERT_FLOAT_NON_EQUAL(max, min, NULL); + + psSpline1D *tmpSpline = (psSpline1D *) psAlloc(sizeof(psSpline1D)); + tmpSpline->n = numSplines; + + // + // XXX: We might have to allocate single or double polynomials depending on the type + // of the psVector bounds. For now, all knots and spline polynomials are 32-bit. + // + tmpSpline->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); + for (psS32 i=0;ispline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD); + } + + // This will be computed by psVectorFitSpline1D() + tmpSpline->p_psDeriv2 = NULL; + + // + // XXX:Ensure that the knots are distinct, and monotonic. + // + tmpSpline->knots = psVectorAlloc(numSplines+1, PS_TYPE_F32); + psF32 width = (max - min) / ((psF32) numSplines); + tmpSpline->knots->data.F32[0] = min; + for (psS32 i=1;iknots->data.F32[i] = min + (width * (psF32) i); + } + tmpSpline->knots->data.F32[numSplines] = max; + + psMemSetDeallocator(tmpSpline, (psFreeFunc)spline1DFree); + return(tmpSpline); +} + +/***************************************************************************** +XXX: Is there a psLib function for this? + *****************************************************************************/ +psVector *PsVectorDup(psVector *in) +{ + psVector *out = psVectorAlloc(in->n, in->type.type); + + if (in->type.type == PS_TYPE_F32) { + for (psS32 i = 0 ; i < in->n ; i++) { + out->data.F32[i] = in->data.F32[i]; + } + } else if (in->type.type == PS_TYPE_F64) { + for (psS32 i = 0 ; i < in->n ; i++) { + out->data.F64[i] = in->data.F64[i]; + } + } else { + printf("XXX: Generate an error here.\n"); + return(NULL); + } + return(out); +} + +/***************************************************************************** +XXX: What should be the default type for knots, spline polys? psF32 is assumed. + *****************************************************************************/ +psSpline1D *psSpline1DAllocGeneric(const psVector *bounds, + int order) +{ + PS_ASSERT_VECTOR_NON_NULL(bounds, NULL); + PS_ASSERT_VECTOR_NON_EMPTY(bounds, NULL); + PS_ASSERT_VECTOR_TYPE(bounds, PS_TYPE_F32, NULL); + PS_ASSERT_INT_NONNEGATIVE(order, NULL); + + psSpline1D *tmpSpline = (psSpline1D *) psAlloc(sizeof(psSpline1D)); + psS32 numSplines = bounds->n - 1; + tmpSpline->n = numSplines; + + // + // XXX: We might have to allocate single or double polynomials depending on the type + // of the psVector bounds. For now, all knots and spline polynomials are 32-bit. + // + tmpSpline->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); + for (psS32 i=0;ispline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD); + } + + // This will be computed by psVectorFitSpline1D() + tmpSpline->p_psDeriv2 = NULL; + + // + // Ensure that all knots are distinct. + // XXX:Ensure that the knots are monotonic. + // + for (psS32 i=0;in-1;i++) { + if (FLT_EPSILON >= fabs(bounds->data.F32[i+1]-bounds->data.F32[i])) { + psError(PS_ERR_UNKNOWN, true, "data points must be distinct ([%d] %f %f)\n", i, bounds->data.F32[i], bounds->data.F32[i+1]); + return(NULL); + } + } + tmpSpline->knots = PsVectorDup((psVector *) bounds); + + psMemSetDeallocator(tmpSpline, (psFreeFunc)spline1DFree); + return(tmpSpline); +} + +/***************************************************************************** +vectorBinDisectF32(): This is a macro for a private function which takes as +input a vector an array of data as well as a single value for that data. The +input vector values are assumed to be non-decreasing (v[i-1] <= v[i] for all +i). This routine does a binary disection of the vector and returns "i" such +that (v[i] <= x <= v[i+1). If x lies outside the range of v[], then this +routine prints a warning message and returns (-2 or -1). + *****************************************************************************/ +#define FUNC_MACRO_VECTOR_BIN_DISECT(TYPE) \ +static psS32 vectorBinDisect##TYPE(ps##TYPE *bins, \ + psS32 numBins, \ + ps##TYPE x) \ +{ \ + psS32 min; \ + psS32 max; \ + psS32 mid; \ + \ + psTrace(".psLib.dataManip.psFunctions.vectorBinDisect##TYPE", 4, \ + "---- Calling vectorBinDisect##TYPE(%f)\n", x); \ + \ + if (x < bins[0]) { \ + psLogMsg(__func__, PS_LOG_WARN, \ + "vectorBinDisect%s(): ordinate %f is outside vector range (%f - %f).", \ + #TYPE, x, bins[0], bins[numBins-1]); \ + return(-2); \ + } \ + \ + if (x > bins[numBins-1]) { \ + psLogMsg(__func__, PS_LOG_WARN, \ + "vectorBinDisect%s(): ordinate %f is outside vector range (%f - %f).", \ + #TYPE, x, bins[0], bins[numBins-1]); \ + return(-1); \ + } \ + \ + min = 0; \ + max = numBins-2; \ + mid = ((max+1)-min)/2; \ + \ + while (min != max) { \ + psTrace(".psLib.dataManip.psFunctions.vectorBinDisect##TYPE", 4, \ + "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n", \ + min, mid, max, x, bins[mid]); \ + \ + if (x == bins[mid]) { \ + psTrace(".psLib.dataManip.psFunctions.vectorBinDisect##TYPE", 4, \ + "---- Exiting vectorBinDisect##TYPE(): bin %d\n", mid); \ + return(mid); \ + } else if (x < bins[mid]) { \ + max = mid-1; \ + } else { \ + min = mid; \ + } \ + mid = ((max+1)+min)/2; \ + } \ + \ + psTrace(".psLib.dataManip.psFunctions.vectorBinDisect##TYPE", 4, \ + "---- Exiting vectorBinDisect##TYPE(): bin %d\n", min); \ + return(min); \ +} \ + +FUNC_MACRO_VECTOR_BIN_DISECT(S8) +FUNC_MACRO_VECTOR_BIN_DISECT(S16) +FUNC_MACRO_VECTOR_BIN_DISECT(S32) +FUNC_MACRO_VECTOR_BIN_DISECT(S64) +FUNC_MACRO_VECTOR_BIN_DISECT(U8) +FUNC_MACRO_VECTOR_BIN_DISECT(U16) +FUNC_MACRO_VECTOR_BIN_DISECT(U32) +FUNC_MACRO_VECTOR_BIN_DISECT(U64) +FUNC_MACRO_VECTOR_BIN_DISECT(F32) +FUNC_MACRO_VECTOR_BIN_DISECT(F64) + +/***************************************************************************** +p_psVectorBinDisect(): A wrapper to the above p_psVectorBinDisect(). + +XXX: Assert that the psVector and psScalar have the same type. + *****************************************************************************/ +psS32 p_psVectorBinDisect(psVector *bins, + psScalar *x) +{ + PS_ASSERT_VECTOR_NON_NULL(bins, -4); + PS_ASSERT_VECTOR_NON_EMPTY(bins, -4); + PS_ASSERT_PTR_NON_NULL(x, -6); + PS_ASSERT_PTR_TYPE_EQUAL(x, bins, -3); + char* strType; + + switch (x->type.type) { + case PS_TYPE_U8: + return(vectorBinDisectU8(bins->data.U8, bins->n, x->data.U8)); + case PS_TYPE_U16: + return(vectorBinDisectU16(bins->data.U16, bins->n, x->data.U16)); + case PS_TYPE_U32: + return(vectorBinDisectU32(bins->data.U32, bins->n, x->data.U32)); + case PS_TYPE_U64: + return(vectorBinDisectU64(bins->data.U64, bins->n, x->data.U64)); + case PS_TYPE_S8: + return(vectorBinDisectS8(bins->data.S8, bins->n, x->data.S8)); + case PS_TYPE_S16: + return(vectorBinDisectS16(bins->data.S16, bins->n, x->data.S16)); + case PS_TYPE_S32: + return(vectorBinDisectS32(bins->data.S32, bins->n, x->data.S32)); + case PS_TYPE_S64: + return(vectorBinDisectS64(bins->data.S64, bins->n, x->data.S64)); + case PS_TYPE_F32: + return(vectorBinDisectF32(bins->data.F32, bins->n, x->data.F32)); + case PS_TYPE_F64: + return(vectorBinDisectF64(bins->data.F64, bins->n, x->data.F64)); + case PS_TYPE_C32: + PS_TYPE_NAME(strType,x->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, + PS_ERRORTEXT_psFunctions_TYPE_NOT_SUPPORTED, + strType); + return 0; + case PS_TYPE_C64: + PS_TYPE_NAME(strType,x->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, + PS_ERRORTEXT_psFunctions_TYPE_NOT_SUPPORTED, + strType); + return 0; + case PS_TYPE_BOOL: + PS_TYPE_NAME(strType,x->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, + PS_ERRORTEXT_psFunctions_TYPE_NOT_SUPPORTED, + strType); + return 0; + } + return(-3); +} + +/***************************************************************************** +p_psVectorInterpolate(): This routine will take as input psVectors domain and +range, and the x value, assumed to lie with the domain vector. It produces +as output the LaGrange interpolated value of a polynomial of the specified +order around the point x. + +XXX: This stuff does not currently work with a mask. + +XXX: add another psScalar argument for the result. + +XXX: The VectorCopy routines seg fault when I declare range32 as static. + *****************************************************************************/ +psScalar *p_psVectorInterpolate(psVector *domain, + psVector *range, + int order, + psScalar *x) +{ + PS_ASSERT_VECTOR_NON_NULL(domain, NULL); + PS_ASSERT_VECTOR_NON_NULL(range, NULL); + PS_ASSERT_PTR_NON_NULL(x, NULL); + PS_ASSERT_INT_NONNEGATIVE(order, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(domain, range, NULL); + PS_ASSERT_PTR_TYPE_EQUAL(domain, range, NULL); + PS_ASSERT_PTR_TYPE_EQUAL(domain, x, NULL); + + psVector *range32 = NULL; + psVector *domain32 = NULL; + psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4, + "---- p_psVectorInterpolate() begin ----\n"); + + if (order > (domain->n - 1)) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, + PS_ERRORTEXT_psFunctions_NOT_ENOUGH_DATAPOINTS, + order); + return(NULL); + } + + if (x->type.type == PS_TYPE_F32) { + psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4, + "---- p_psVectorInterpolate() end ----\n"); + return(psScalarAlloc(interpolate1DF32(domain->data.F32, + range->data.F32, + domain->n, + order, + x->data.F32), PS_TYPE_F32)); + } else if (x->type.type == PS_TYPE_F64) { + // XXX: use recycled vectors here. + range32 = psVectorCopy(range32, range, PS_TYPE_F32); + domain32 = psVectorCopy(domain32, domain, PS_TYPE_F32); + + psScalar *tmpScalar = psScalarAlloc((psF64) + interpolate1DF32(domain32->data.F32, + range32->data.F32, + domain32->n, + order, + (psF32) x->data.F64), PS_TYPE_F64); + psFree(range32); + psFree(domain32); + + psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4, + "---- p_psVectorInterpolate() end ----\n"); + // XXX: Convert data type to F64? + return(tmpScalar); + + } else { + char* strType; + PS_TYPE_NAME(strType,x->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, + PS_ERRORTEXT_psFunctions_TYPE_NOT_SUPPORTED, + strType); + } + + psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4, + "return(NULL)\n"); + psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4, + "---- p_psVectorInterpolate() end ----\n"); + + return(NULL); +} + + +/***************************************************************************** +psSpline1DEval(): this routine takes an existing spline of arbitrary order +and an independent x value. Each 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 +by a call to the 1D polynomial functions. + +XXX: The spline eval functions require input and output to be F32. however + the spline fit functions require F32 and F64. + +XXX: This only works if spline0>knots if psF32. Must add support for psU32 and +psF64. + *****************************************************************************/ +float psSpline1DEval( + const psSpline1D *spline, + float x +) +{ + PS_ASSERT_PTR_NON_NULL(spline, NAN); + PS_ASSERT_INT_NONNEGATIVE(spline->n, NAN); + PS_ASSERT_VECTOR_TYPE(spline->knots, PS_TYPE_F32, NAN); + + unsigned int binNum; + unsigned int n; + + n = spline->n; + //XXX binNum = vectorBinDisectF32(spline->domains, (spline->n)+1, x); + binNum = vectorBinDisectF32(spline->knots->data.F32, (spline->n)+1, x); + if (binNum < 0) { + psLogMsg(__func__, PS_LOG_WARN, + "psSpline1DEval(): x ordinate (%f) is outside the spline range (%f - %f).", + x, spline->knots->data.F32[0], + spline->knots->data.F32[n-1]); + + if (x < spline->knots->data.F32[0]) { + return(psPolynomial1DEval(spline->spline[0], + x)); + } else if (x > spline->knots->data.F32[n-1]) { + return(psPolynomial1DEval(spline->spline[n-1], + x)); + } + } + + return(psPolynomial1DEval(spline->spline[binNum], + x)); +} + +// XXX: The spline eval functions require input and output to be F32. +// however the spline fit functions require F32 and F64. +psVector *psSpline1DEvalVector( + const psSpline1D *spline, + const psVector *x +) +{ + PS_ASSERT_PTR_NON_NULL(spline, NULL); + PS_ASSERT_VECTOR_NON_NULL(x, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL); + PS_ASSERT_VECTOR_TYPE(spline->knots, PS_TYPE_F32, NULL); + + unsigned int i; + psVector *tmpVector; + + tmpVector = psVectorAlloc(x->n, PS_TYPE_F32); + if (x->type.type == PS_TYPE_F32) { + for (i=0;in;i++) { + tmpVector->data.F32[i] = psSpline1DEval( + spline, + x->data.F32[i] + ); + } + } else if (x->type.type == PS_TYPE_F64) { + for (i=0;in;i++) { + tmpVector->data.F32[i] = psSpline1DEval( + spline, + (psF32) x->data.F64[i] + ); + } + } else { + char* strType; + PS_TYPE_NAME(strType,x->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, + PS_ERRORTEXT_psFunctions_TYPE_NOT_SUPPORTED, + strType); + return(NULL); + } + + return(tmpVector); +} Index: trunk/psLib/src/math/psSpline.h =================================================================== --- trunk/psLib/src/math/psSpline.h (revision 4958) +++ trunk/psLib/src/math/psSpline.h (revision 4958) @@ -0,0 +1,113 @@ +/** @file psSpline.h + * @brief Standard Mathematical Functions. + * @ingroup GROUP00 + * + * This file will hold the prototypes for procedures which allocate, free, + * and evaluate splines. + * + * @ingroup GROUP00 + * + * @author GLG, MHPCC + * + * @version $Revision: 1.53 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-07 21:35:50 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ + +#ifndef PS_SPLINE_H +#define PS_SPLINE_H + +#include +#include +#include +#include + +#include "psVector.h" +#include "psScalar.h" +#include "psPolynomial.h" + +/** \addtogroup GROUP00 + * \{ + */ + +/** One-Dimensional Spline */ +typedef struct +{ + int n; ///< The number of spline pieces + psPolynomial1D **spline; ///< An array of n pointers to the spline polynomials + psVector *knots; ///< The boundaries between each spline piece. Size is n+1. + psF32 *p_psDeriv2; ///< For cubic splines, the second derivative at each domain point. Size is n+1. + psF32 *p_psDomains; ///< The boundaries between each spline piece. Size is n+1. +} +psSpline1D; + +/** Allocates a psSpline1D structure + * + * Allocator for psSpline1D where the bounds are implicitly specified through specifying + * min and max values along with the number of splines. + * + * @return psSpline1D* new 1-D spline struct + */ +psSpline1D *psSpline1DAlloc( + int n, ///< Number of spline polynomials + int order, ///< Order of spline polynomials + float min, ///< Lower boundary value of spline polynomials + float max ///< Upper boundary value of spline polynomials +); + +/** Allocates a psSpline1D structure + * + * Allocator for psSpline1D where the bounds are explicitly specified. + * + * @return psSpline1D* new 1-D spline struct + */ +psSpline1D *psSpline1DAllocGeneric( + const psVector *bounds, ///< Bounds for spline polynomials + int order ///< Order of spline polynomials +); + +/** Evaluates 1-D spline polynomials at a specific coordinate. + * + * @return float result of spline polynomials evaluated at given location + */ +float psSpline1DEval( + const psSpline1D *spline, ///< Coefficients for spline polynomials + float x ///< location at which to evaluate +); + +/** Evaluates 1-D spline polynomials at a set of specific coordinates. + * + * @return psVector* results of spline polynomials evaluated at given locations + */ +psVector *psSpline1DEvalVector( + const psSpline1D *spline, ///< Coefficients of spline polynomials + const psVector *x ///< locations at which to evaluate +); + +/** Performs a binary disection on a given vector. + * Searches through an array of data for a specified value. + * + * @return psS32 corresponding index number of specified value + */ +psS32 p_psVectorBinDisect( + psVector *bins, ///< Array of non-decreasing values + psScalar *x ///< Target value to find +); + +/** Interpolates a series of data points for evaluation at a specific coordinate. Uses a + * Lagrange interpolation method. + * + * @return psScalar* Lagrange interpolation value at given location + */ +psScalar *p_psVectorInterpolate( + psVector *domain, ///< Domain (x coords) for interpolation + psVector *range, ///< Range (y coords) for interpolation + int order, ///< Order of interpolation function + psScalar *x ///< Location at which to evaluate +); + +/** \} */ // End of MathGroup Functions + +#endif // #ifndef PS_SPLINE_H + Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 4953) +++ trunk/psLib/src/math/psStats.c (revision 4958) @@ -14,6 +14,6 @@ * stats->binsize * - * @version $Revision: 1.139 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-08-30 01:14:13 $ + * @version $Revision: 1.140 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-07 21:35:50 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -46,5 +46,6 @@ /* DEFINE STATEMENTS */ /*****************************************************************************/ -#define PS_GAUSS_WIDTH 5 // The width of the Gaussian or boxcar smoothing. +#define PS_GAUSS_WIDTH 5 // The width of the Gaussian smoothing. +// This corresponds to N in the ADD. #define PS_CLIPPED_NUM_ITER_LB 1 #define PS_CLIPPED_NUM_ITER_UB 10 @@ -292,4 +293,6 @@ max of the input vector. If there was a problem with the max calculation, this routine sets stats->max to NAN. + +XXX: Do we need to factor errors into it? *****************************************************************************/ psS32 p_psVectorMax(const psVector* myVector, @@ -614,5 +617,6 @@ /****************************************************************************** p_psVectorSmoothHistGaussian(): This routine smoothes the data in the input -robustHistogram with a Gaussian of width sigma. +robustHistogram with a Gaussian of width sigma. It returns a psVector of the +smoothed data. XXX: Only PS_TYPE_F32 is supported. @@ -621,77 +625,227 @@ call that. Is that possible? *****************************************************************************/ -psVector* p_psVectorSmoothHistGaussian(psHistogram* robustHistogram, +psVector *p_psVectorSmoothHistGaussian(psHistogram *histogram, psF32 sigma) { - PS_ASSERT_PTR_NON_NULL(robustHistogram, NULL); - PS_ASSERT_PTR_NON_NULL(robustHistogram->bounds, NULL); - - psS32 i = 0; // Loop index variable - psS32 j = 0; // Loop index variable - psF32 iMid; - psF32 jMid; - psS32 numBins = robustHistogram->nums->n; - psS32 numBounds = robustHistogram->bounds->n; - psVector* smooth = psVectorAlloc(numBins, PS_TYPE_F32); + PS_ASSERT_PTR_NON_NULL(histogram, NULL); + PS_ASSERT_PTR_NON_NULL(histogram->bounds, NULL); + PS_ASSERT_PTR_NON_NULL(histogram->nums, NULL); + + psS32 numBins = histogram->nums->n; + psS32 numBounds = histogram->bounds->n; + psVector *smooth = psVectorAlloc(numBins, PS_TYPE_F32); + psF32 firstBound = histogram->bounds->data.F32[0]; + psF32 lastBound = histogram->bounds->data.F32[numBounds-1]; + psScalar x; + x.type.type = PS_TYPE_F32; psS32 jMin = 0; psS32 jMax = 0; - psF32 firstBound = robustHistogram->bounds->data.F32[0]; - psF32 lastBound = robustHistogram->bounds->data.F32[numBounds-1]; + + if (histogram->uniform == false) { + // + // We get here if the histogram is non-uniform. + // + + for (psS32 i = 0; i < numBins; i++) { + // Determine the midpoint of bin i. + psS32 iMid = PS_BIN_MIDPOINT(histogram, i); + + // + // We determine the bin numbers (jMin:jMax) corresponding to a + // range of data values surrounding iMid. The range is of size: + // 2*PS_GAUSS_WIDTH*sigma + // + x.data.F32 = iMid - (PS_GAUSS_WIDTH * sigma); + if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { + jMin = p_psVectorBinDisect( *(psVector* *)&histogram->bounds, &x); + if (jMin < 0) { + psError(PS_ERR_UNEXPECTED_NULL, + false, + PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM); + return(NULL); + } + } else if (x.data.F32 <= firstBound) { + jMin = 0; + } else if (x.data.F32 >= lastBound) { + jMin = histogram->bounds->n - 1; + } + + x.data.F32 = iMid + (PS_GAUSS_WIDTH * sigma); + if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { + jMax = p_psVectorBinDisect( *(psVector* *)&histogram->bounds, &x); + if (jMax < 0) { + psError(PS_ERR_UNEXPECTED_NULL, + false, + PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM); + return(NULL); + } + } else if (x.data.F32 <= firstBound) { + jMax = 0; + } else if (x.data.F32 >= lastBound) { + jMax = histogram->bounds->n - 1; + } + + // + // Loop from jMin to jMax, computing the gaussian of data i. + // + smooth->data.F32[i] = 0.0; + for (psS32 j = jMin ; j <= jMax ; j++) { + psS32 jMid = PS_BIN_MIDPOINT(histogram, j); + smooth->data.F32[i] += + histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true); + } + } + } else { + // + // We get here if the histogram is uniform. + // + + for (psS32 i = 0; i < numBins; i++) { + psF32 binSize = histogram->bounds->data.F32[1] - histogram->bounds->data.F32[0]; + psS32 gaussWidth = (psS32) ((PS_GAUSS_WIDTH * sigma) / binSize); + + // + // We determine the bin numbers (jMin:jMax) corresponding to a + // range of data values surrounding iMid. The range is of size: + // 2*PS_GAUSS_WIDTH*sigma + // + psS32 jMin = i - gaussWidth; + if (jMin < 0 ) { + jMin = 0; + } + psS32 jMax = i + gaussWidth; + if (jMax > (histogram->bounds->n - 1)) { + jMax = (histogram->bounds->n - 1); + } + + // + // Loop from jMin to jMax, computing the gaussian of data i. + // + smooth->data.F32[i] = 0.0; + psS32 iMid = PS_BIN_MIDPOINT(histogram, i); + for (psS32 j = jMin ; j <= jMax ; j++) { + psS32 jMid = PS_BIN_MIDPOINT(histogram, j); + smooth->data.F32[i] += + histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true); + } + } + } + + return(smooth); +} +/****************************************************************************** +p_psVectorSmoothHistGaussianNEW(): This routine smoothes the data in the input +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_psVectorSmoothHistGaussianNEW(psHistogram *histogram, + psF32 sigma) +{ + PS_ASSERT_PTR_NON_NULL(histogram, NULL); + PS_ASSERT_PTR_NON_NULL(histogram->bounds, NULL); + PS_ASSERT_PTR_NON_NULL(histogram->nums, NULL); + + psS32 numBins = histogram->nums->n; + psS32 numBounds = histogram->bounds->n; + psVector *smooth = psVectorAlloc(numBins, PS_TYPE_F32); + psF32 firstBound = histogram->bounds->data.F32[0]; + psF32 lastBound = histogram->bounds->data.F32[numBounds-1]; psScalar x; - x.type.type = PS_TYPE_F32; - for (i = 0; i < numBins; i++) { - // Determine the midpoint of bin i. - iMid = (robustHistogram->bounds->data.F32[i] + - robustHistogram->bounds->data.F32[i+1]) / 2.0; - - - // We determine the bin numbers corresponding to a range of data - // values surrounding iMid. The ranges is of size - // s*PS_GAUSS_WIDTH*sigma - - // YYY: The p_psVectorBinDisect() routine does much of the work of - // the following conditionals, however, it also reports a warning - // message. I don't want the warning message so I reproduce the - // conditionals here. Maybe p_psVectorBinDisect() should not produce - // warnings? - - x.data.F32 = iMid - (PS_GAUSS_WIDTH * sigma); - if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { - jMin = p_psVectorBinDisect( *(psVector* *)&robustHistogram->bounds, &x); - if (jMin < 0) { - psError(PS_ERR_UNEXPECTED_NULL, - false, - PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM); - return(NULL); - } - } else if (x.data.F32 <= firstBound) { - jMin = 0; - } else if (x.data.F32 >= lastBound) { - jMin = robustHistogram->bounds->n - 1; - } - - x.data.F32 = iMid + (PS_GAUSS_WIDTH * sigma); - if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { - jMax = p_psVectorBinDisect( *(psVector* *)&robustHistogram->bounds, &x); - if (jMax < 0) { - psError(PS_ERR_UNEXPECTED_NULL, - false, - PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM); - return(NULL); - } - } else if (x.data.F32 <= firstBound) { - jMax = 0; - } else if (x.data.F32 >= lastBound) { - jMax = robustHistogram->bounds->n - 1; - } - - smooth->data.F32[i] = 0.0; - for (j = jMin ; j <= jMax ; j++) { - jMid = (robustHistogram->bounds->data.F32[j] + - robustHistogram->bounds->data.F32[j+1]) / 2.0; - smooth->data.F32[i] += - robustHistogram->nums->data.F32[j] * - psGaussian(jMid, iMid, sigma, true); + psS32 jMin = 0; + psS32 jMax = 0; + + if (histogram->uniform == false) { + // + // We get here if the histogram is non-uniform. + // + + for (psS32 i = 0; i < numBins; i++) { + // Determine the midpoint of bin i. + psS32 iMid = PS_BIN_MIDPOINT(histogram, i); + + // + // We determine the bin numbers (jMin:jMax) corresponding to a + // range of data values surrounding iMid. The range is of size: + // 2*PS_GAUSS_WIDTH*sigma + // + x.data.F32 = iMid - (PS_GAUSS_WIDTH * sigma); + if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { + jMin = p_psVectorBinDisect( *(psVector* *)&histogram->bounds, &x); + if (jMin < 0) { + psError(PS_ERR_UNEXPECTED_NULL, + false, + PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM); + return(NULL); + } + } else if (x.data.F32 <= firstBound) { + jMin = 0; + } else if (x.data.F32 >= lastBound) { + jMin = histogram->bounds->n - 1; + } + + x.data.F32 = iMid + (PS_GAUSS_WIDTH * sigma); + if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { + jMax = p_psVectorBinDisect( *(psVector* *)&histogram->bounds, &x); + if (jMax < 0) { + psError(PS_ERR_UNEXPECTED_NULL, + false, + PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM); + return(NULL); + } + } else if (x.data.F32 <= firstBound) { + jMax = 0; + } else if (x.data.F32 >= lastBound) { + jMax = histogram->bounds->n - 1; + } + + // + // Loop from jMin to jMax, computing the gaussian of data i. + // + smooth->data.F32[i] = 0.0; + for (psS32 j = jMin ; j <= jMax ; j++) { + psS32 jMid = PS_BIN_MIDPOINT(histogram, j); + smooth->data.F32[i] += + histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true); + } + } + } else { + // + // We get here if the histogram is uniform. + // + + for (psS32 i = 0; i < numBins; i++) { + psF32 binSize = histogram->bounds->data.F32[1] - histogram->bounds->data.F32[0]; + psS32 gaussWidth = (psS32) ((PS_GAUSS_WIDTH * sigma) / binSize); + + // + // We determine the bin numbers (jMin:jMax) corresponding to a + // range of data values surrounding iMid. The range is of size: + // 2*PS_GAUSS_WIDTH*sigma + // + psS32 jMin = i - gaussWidth; + if (jMin < 0 ) { + jMin = 0; + } + psS32 jMax = i + gaussWidth; + if (jMax > (histogram->bounds->n - 1)) { + jMax = (histogram->bounds->n - 1); + } + + // + // Loop from jMin to jMax, computing the gaussian of data i. + // + smooth->data.F32[i] = 0.0; + psS32 iMid = PS_BIN_MIDPOINT(histogram, i); + for (psS32 j = jMin ; j <= jMax ; j++) { + psS32 jMid = PS_BIN_MIDPOINT(histogram, j); + smooth->data.F32[i] += + histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true); + } } } @@ -1142,5 +1296,5 @@ These macros and functions define the following functions: -< p_psNormalizeVectorRange(myData, low, high) + p_psNormalizeVectorRange(myData, low, high) That assumes that the low/high arguments are PS_TYPE_F64; the vector myData @@ -1251,11 +1405,10 @@ /****************************************************************************** -p_ps1DPolyMedian(myPoly, rangeLow, rangeHigh, midpoint): This routine takes -as input a 1-D polynomial of arbitrary order (though we are using 2nd-order -polynomials here) and a range of x-values for which it is defined: -[rangeLow, rangeHigh]. It determines the x-value of that polynomial such -that f(x) == midpoint. This functions uses a binary-search algorithm on the -range and assumes that the polynomial is monotonically increasing or -decreasing within that range. +p_ps1DPolyMedian(myPoly, rangeLow, rangeHigh, getThisValue): This routine +takes as input a 1-D polynomial of arbitrary order and a range of x-values for +which it is defined: [rangeLow, rangeHigh]. It determines the x-value of +that polynomial such that f(x) == getThisValue. This function uses a +binary-search algorithm on the range and assumes that the polynomial is +monotonically increasing or decreasing within that range. XXX: Terminate when f(x)-getThisValue is within some error tolerance. @@ -1269,16 +1422,10 @@ { PS_ASSERT_POLY_NON_NULL(myPoly, NAN); - PS_FLOAT_COMPARE(rangeLow, rangeHigh, NAN); + PS_ASSERT_FLOAT_LARGER_THAN(rangeHigh, rangeLow, NAN); // We ensure that the requested f(y) value, which is getThisValue, is // falls within the range of y-values of the polynomial "myPoly" in the // specified x-range (rangeLow:rangeHigh). - psF32 fLo = psPolynomial1DEval( - myPoly, - rangeLow - ); - psF32 fHi = psPolynomial1DEval( - myPoly, - rangeHigh - ); + psF32 fLo = psPolynomial1DEval(myPoly, rangeLow); + psF32 fHi = psPolynomial1DEval(myPoly, rangeHigh); if (!((fLo <= getThisValue) && (fHi >= getThisValue))) { psError(PS_ERR_UNKNOWN, @@ -1300,8 +1447,5 @@ oldMidpoint = midpoint; - f = psPolynomial1DEval( - myPoly, - midpoint - ); + f = psPolynomial1DEval(myPoly, midpoint); if (fabs(f - getThisValue) <= FLT_EPSILON) { return (midpoint); @@ -1331,4 +1475,8 @@ XXX: the vectors do not have to be the same length. Must insert the proper tests to ensure that binNum is within acceptable ranges for both vectors. + +XXX: This currently assumes that the three points are monotonically increasing +or decreasing: so, it works for the cumulative histogram vectors, but not for +arbitrary vectors. We should probably test that condition. *****************************************************************************/ psF32 fitQuadraticSearchForYThenReturnX(psVector *xVec, @@ -1345,8 +1493,4 @@ PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (yVec->n - 1), NAN); - // PS_VECTOR_DECLARE_ALLOC_STATIC(x, 3, PS_TYPE_F64); - // PS_VECTOR_DECLARE_ALLOC_STATIC(y, 3, PS_TYPE_F64); - // PS_VECTOR_DECLARE_ALLOC_STATIC(yErr, 3, PS_TYPE_F64); - // PS_POLY_1D_DECLARE_ALLOC_STATIC(myPoly, 2, PS_POLYNOMIAL_ORD); psVector *x = psVectorAlloc(3, PS_TYPE_F64); psVector *y = psVectorAlloc(3, PS_TYPE_F64); @@ -1364,4 +1508,30 @@ y->data.F64[1] = yVec->data.F32[binNum]; y->data.F64[2] = yVec->data.F32[binNum + 1]; + + // + // Ensure that the y values are monotonic. + // + // XXX: This routine should probably be rewritten in a more general fashion + // so that the folloiwng checks are not necessary. + // + if (y->data.F64[0] < y->data.F64[1]) { + if (!(y->data.F64[1] <= y->data.F64[2])) { + psError(PS_ERR_UNKNOWN, true, "This routine must be called with montically increasing or decreasing data points.\n"); + psFree(myPoly); + psFree(x); + psFree(y); + psFree(yErr); + return(NAN); + } + } else { + if (!(y->data.F64[1] >= y->data.F64[2])) { + psError(PS_ERR_UNKNOWN, true, "This routine must be called with montically increasing or decreasing data points.\n"); + psFree(myPoly); + psFree(x); + psFree(y); + psFree(yErr); + return(NAN); + } + } // Ensure that yVal is within the range of the bins we are using. @@ -1637,5 +1807,5 @@ // code no longer produces sensible results. // XXX: Since we are no longer fitting a 1-D Gaussian, we can probably - // remove some of theabove code that calculated the initial estimate + // remove some of the above code that calculated the initial estimate // for the mean and sigma. @@ -1645,6 +1815,6 @@ int index = i - modeBinNum + dL; // XXX: Should this be the natural log? - y->data.F32[index] = robustHistogramVector->data.F32[i]; - // y->data.F32[index] = logf(robustHistogramVector->data.F32[i]); + // y->data.F32[index] = robustHistogramVector->data.F32[i]; + y->data.F32[index] = logf(robustHistogramVector->data.F32[i]); x->data.F32[index] = (psF32) index; } @@ -1742,6 +1912,363 @@ psFree(robustHistogramVector); psFree(cumulativeRobustSums); + return(0); } + + +/***************************************************************************** +XXX: Is there a psLib function for this? + *****************************************************************************/ +psVector *PsVectorDup(psVector *in) +{ + psVector *out = psVectorAlloc(in->n, in->type.type); + + if (in->type.type == PS_TYPE_F32) { + for (psS32 i = 0 ; i < in->n ; i++) { + out->data.F32[i] = in->data.F32[i]; + } + } else if (in->type.type == PS_TYPE_F64) { + for (psS32 i = 0 ; i < in->n ; i++) { + out->data.F64[i] = in->data.F64[i]; + } + } else { + printf("XXX: Generate an error here.\n"); + return(NULL); + } + return(out); +} + +/****************************************************************************** +XXX: This function need to be written. Actually, it simply needs to be +retrieved from the CVS repository, since it was written earlier, then +discarded. At present, it was deleted from the CVS repository, so we might +have to retrieve it from tape. +*****************************************************************************/ +psVector *Fit1DGaussian(psVector *x, psVector*y) +{ + printf("XXX: Generate an error here.\n"); + printf("XXX: Error: This function was previously part of psStats.c, was removed, was purged from CVS, and now needs to be retrieved from tape.\n"); + return(NULL); +} + +/****************************************************************************** + +p_psVectorRobustStatsNew(myVector, maskVector, maskVal, stats): This is the new +version of the robust stats routine. + +XXX: MUST DO: If the errors in the input values are known, then the same +approach is used, except that the histograms become probability density +functions (PDFs). In this case, the input values are spread out, so that they +do not simply contribute a single unit to the histogram, but rather contribute +a fraction of a value, equivalent to the weight. In the interests of speed, a +boxcar PDF may be used to represent each input value (as opposed to a +Gaussian), where the boxcar width is equal to 2p2 ln 2 times the error and +each input value contributes constant area. Then the robust median and +standard deviation are estimated in the same manner as above. + +XXX: Check for errors in psLib routines that we call. +*****************************************************************************/ +psS32 p_psVectorRobustStatsNew(const psVector* myVector, + const psVector* errors, + const psVector* maskVector, + psU32 maskVal, + psStats* stats) +{ + psHistogram *robustHistogram = NULL; + psHistogram *cumulativeRobustHistogram = NULL; + psS32 numBins = 0; + psScalar *tmpScalar = psScalarAlloc(0.0, PS_TYPE_F32); + tmpScalar->type.type = PS_TYPE_F32; + psS32 totalDataPoints = 0; + psS32 rc = 0; + psVector *tmpMaskVec = PsVectorDup((psVector *) maskVector); + + while (1) { + // + // Determine the bin size of the robust histogram. This is done + // by computing the total range of data values and dividing by 1000.0. + // + psStats* tmpStatsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); + rc = p_psVectorMin(myVector, tmpMaskVec, maskVal, tmpStatsMinMax); + rc|= p_psVectorMax(myVector, tmpMaskVec, maskVal, 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(tmpStatsMinMax); + psFree(tmpMaskVec); + psFree(tmpScalar); + return(1); + } + psF32 binSize = (tmpStatsMinMax->max - tmpStatsMinMax->min) / 1000.0f; + + // + // If all data points have the same value, then we set the appropiate + // members of stats and return. + // + if (fabs(tmpStatsMinMax->max - tmpStatsMinMax->min) <= FLT_EPSILON) { + if (stats->options & PS_STAT_ROBUST_MEDIAN) { + stats->robustMedian = tmpStatsMinMax->min; + } + if (stats->options & PS_STAT_ROBUST_QUARTILE) { + stats->robustUQ = tmpStatsMinMax->min; + stats->robustLQ = tmpStatsMinMax->min; + } + // XXX: Set these to the number of unmasked data points? + stats->robustNfit = 0.0; + stats->robustN50 = 0.0; + psFree(tmpStatsMinMax); + psFree(tmpMaskVec); + psFree(tmpScalar); + + return(0); + } + + // + // ADD: Step 0. + // Construct the histogram with the specified bin size. + // + // NOTE: we can not specify the bin size precisely since the argument + // to psHistogramAlloc() is the number of bins, not the binSize. + // If we get here, we know that binSize != 0.0. + // + numBins = (psS32)((tmpStatsMinMax->max - tmpStatsMinMax->min) / binSize); + robustHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins); + cumulativeRobustHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins); + + // Populate the histogram array. + psVectorHistogram(robustHistogram, myVector, errors, tmpMaskVec, maskVal); + + // + // ADD: Step 1. + // Construct the cumulative histogram from the specific histogram + // + cumulativeRobustHistogram->nums->data.F32[0] = robustHistogram->nums->data.F32[0]; + for (psS32 i = 1 ; i < robustHistogram->nums->n ; i++) { + cumulativeRobustHistogram->nums->data.F32[i] = cumulativeRobustHistogram->nums->data.F32[i-1] + + robustHistogram->nums->data.F32[i]; + } + + // + // ADD: Step 2. + // Find the bin which contains the 50% data point. + // + totalDataPoints = cumulativeRobustHistogram->nums->data.F32[numBins - 1]; + tmpScalar->data.F32 = totalDataPoints/2.0; + psS32 binMedian = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); + if (binMedian != 0) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 50% data point.\n"); + psFree(tmpStatsMinMax); + psFree(robustHistogram); + psFree(cumulativeRobustHistogram); + psFree(tmpScalar); + return(1); + } + + // + // ADD: Step 3. + // Interpolate to the exact 50% position: this is the robust histogram median. + // XXX: Check for errors here! + // + stats->robustMedian = fitQuadraticSearchForYThenReturnX( + *(psVector* *)&robustHistogram->bounds, + *(psVector* *)&robustHistogram->nums, + binMedian, + totalDataPoints/2.0); + + // + // ADD: Step 4. + // Find the bins which contains the 15.8655% and 84.1345% data points. + // + tmpScalar->data.F32 = totalDataPoints * 0.158655f; + psS32 binLo = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); + tmpScalar->data.F32 = totalDataPoints * 0.841345f; + psS32 binHi = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); + if ((binLo != 0) || (binHi != 0)) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate the15.8655% and 84.1345% data point\n"); + psFree(tmpStatsMinMax); + psFree(robustHistogram); + psFree(cumulativeRobustHistogram); + psFree(tmpScalar); + return(1); + } + + // + // ADD: Step 4b. + // Interpolate Sigma to find these two positions exactly: these are the 1sigma positions. + // + psF32 binLoF32 = fitQuadraticSearchForYThenReturnX( + *(psVector* *)&robustHistogram->bounds, + *(psVector* *)&robustHistogram->nums, + binLo, + totalDataPoints * 0.158655f); + psF32 binHiF32 = fitQuadraticSearchForYThenReturnX( + *(psVector* *)&robustHistogram->bounds, + *(psVector* *)&robustHistogram->nums, + binHi, + totalDataPoints * 0.841345f); + + // + // ADD: Step 5. + // Determine SIGMA as 1/2 of the distance between these positions. + // + psF32 sigma = (binHiF32 - binLoF32) / 2.0; + + // + // ADD: Step 6. + // If the measured SIGMA is less than 2 times the bin size, exclude + // points which are more than 25 bins from the median, + // recalculate the bin size, and perform the algorithm again. + // + if (sigma < (2 * binSize)) { + psF32 medianLo = robustHistogram->bounds->data.F32[binMedian - 25]; + psF32 medianHi = robustHistogram->bounds->data.F32[binMedian + 25]; + for (psS32 i = 0 ; i < myVector->n ; i++) { + if ((myVector->data.F32[i] < medianLo) || + (myVector->data.F32[i] > medianHi)) { + tmpMaskVec->data.U8[i] = 1; + } + } + } else { + // + // ADD: Step 7. + // Find the bins which contains the 25% and 75% data points. + // + tmpScalar->data.F32 = totalDataPoints * 0.25f; + psS32 binLo25 = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); + tmpScalar->data.F32 = totalDataPoints * 0.75f; + psS32 binHi25 = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); + if ((binLo25 != 0) || (binHi25 != 0)) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 25% and 75% data points\n"); + psFree(tmpStatsMinMax); + psFree(robustHistogram); + psFree(cumulativeRobustHistogram); + psFree(tmpScalar); + return(1); + } + + // + // ADD: Step 8. + // Interpolate to find these two positions exactly: these are the upper + // and lower quartile positions. + // XXX: Check for errors. + // + psF32 binLo25F32 = fitQuadraticSearchForYThenReturnX( + *(psVector* *)&robustHistogram->bounds, + *(psVector* *)&robustHistogram->nums, + binLo25, + totalDataPoints * 0.25f); + psF32 binHi25F32 = fitQuadraticSearchForYThenReturnX( + *(psVector* *)&robustHistogram->bounds, + *(psVector* *)&robustHistogram->nums, + binHi25, + totalDataPoints * 0.75f); + + stats->robustLQ = binLo25F32; + stats->robustUQ = binHi25F32; + // XXX: No idea how to calculate stats->stdev + + // PS_BIN_MIDPOINT(robustHistogram, modeBinNum); + + // XXX: I think sumNfit == sumN50 here. + stats->robustNfit = -1; + stats->robustN50 = -1; + + // + // Perform the Robust Histogram Statistics algorithm above + // + + // + // Smooth the resulting histogram with a Gaussian with SIGMA_s = 1 + // bin. + // + // XXX: SIGMA_s is defined nowhere in the ADD. + // + psF32 SIGMA_S = 1.0; + p_psVectorSmoothHistGaussian(robustHistogram, SIGMA_S); + + // + // Find the bin with the peak value in the range 2 SIGMA of the + // robust histogram median. + // + // XXX: SIGMA is defined nowhere in the ADD. + // + psF32 SIGMA = 2.0; + psS32 binMin = binMedian - (SIGMA * PS_GAUSS_WIDTH); + if (binMin < 0) { + binMin = 0; + } + psS32 binMax = binMedian + (2 * PS_GAUSS_WIDTH); + if (binMin > (robustHistogram->nums->n - 1)) { + binMin = (robustHistogram->nums->n - 1); + } + psS32 binNum = binNum; + psF32 binMaxNums = robustHistogram->nums->data.F32[binNum]; + for (psS32 i = binNum+1 ; i <= binMax ; i++) { + if (robustHistogram->nums->data.F32[i] > binMaxNums) { + binNum = i; + binMaxNums = robustHistogram->nums->data.F32[i]; + } + } + + // + // Fit a Gaussian to the bins in the range 2 SIGMA of the robust + // histogram median. + // + // XXX: SIGMA is defined nowhere in the ADD. + // + psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); + psVector *x = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); + psS32 j = 0; + for (psS32 i = binNum ; i <= binMax ; i++) { + y->data.F32[j] = robustHistogram->nums->data.F32[i]; + x->data.F32[j] = PS_BIN_MIDPOINT(robustHistogram, i); + } + // + // XXX: This function need to be written. Actually, it simply + // needs to be retrieved from the CVS repository, since it was + // written earlier, then discarded. At present, it was deleted + // from the CVS repository, so we might have to retrieve it from + // tape. + // + psVector *params = Fit1DGaussian(x, y); + + // + // The robust mean mean_r is derived directly from the fitted + // Gaussian mean. + // + stats->robustMean = params->data.F32[0]; + + // + // The robust standard deviation, SIGMA_r is determined by + // subtracting the smoothing scale in quadrature: SIGMA_r^2 = SIGMA^2 + // - SIGMA_s^2 + // + // XXX: SIGMA and SIGMA_s are defined nowhere in the ADD. We must figure + // out what they are. + // + stats->robustStdev = sqrt(PS_SQR(SIGMA) - PS_SQR(SIGMA_S)); + + psFree(tmpStatsMinMax); + psFree(robustHistogram); + psFree(cumulativeRobustHistogram); + psFree(tmpScalar); + psFree(params); + + return(0); + } + + psFree(tmpStatsMinMax); + psFree(robustHistogram); + psFree(cumulativeRobustHistogram); + } + return(1); +} + + + + + + + + /*****************************************************************************/