Index: trunk/psLib/src/math/psSpline.c =================================================================== --- trunk/psLib/src/math/psSpline.c (revision 1476) +++ trunk/psLib/src/math/psSpline.c (revision 1718) @@ -1,3 +1,2 @@ - /** @file psFunctions.c * @@ -8,6 +7,6 @@ * polynomials. It also contains a Gaussian functions. * - * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-08-11 20:07:45 $ + * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-09-08 06:00:49 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -32,4 +31,5 @@ #include "psAbort.h" #include "psFunctions.h" +#include "psLogMsg.h" #include @@ -42,5 +42,43 @@ /*****************************************************************************/ -// None +#define PS_CONVERT_VECTOR_F32(IN, OUT) \ +if (IN->type.type != PS_TYPE_F32) { \ + psLogMsg(__func__, PS_LOG_WARN, \ + "Input vector has incorrect type (%d).\n", \ + IN->type); \ + OUT = psVectorAlloc(IN->n, PS_TYPE_F32); \ + \ + if (IN->type.type == PS_TYPE_F64) { \ + for (int i=0;in;i++) { \ + OUT->data.F32[i] = (float) (IN->data.F64[i]); \ + } \ + }\ + else { \ + psAbort(__func__, "Wrong type.\n"); \ + } \ +} else { \ + OUT = (psVector *) IN; \ +} + +#define PS_CONVERT_VECTOR_F64(IN, OUT) \ +if (IN->type.type != PS_TYPE_F64) { \ + psLogMsg(__func__, PS_LOG_WARN, \ + "Input vector has incorrect type (%d).\n", \ + IN->type); \ + OUT = psVectorAlloc(IN->n, PS_TYPE_F64); \ + \ + if (IN->type.type == PS_TYPE_F32) { \ + for (int i=0;in;i++) { \ + OUT->data.F64[i] = (float) (IN->data.F32[i]); \ + } \ + }\ + else { \ + psAbort(__func__, "Wrong type.\n"); \ + } \ +} else { \ + OUT = (psVector *) IN; \ +} + + /*****************************************************************************/ @@ -107,5 +145,5 @@ } - return (tmp * exp(-((x - mean) * (x - mean)) / (2.0 * sigma * sigma))); + return(tmp * exp(-((x - mean) * (x - mean)) / (2.0 * sigma * sigma))); } @@ -138,11 +176,18 @@ // NOTE: Should I free r as well? - return (gauss); -} + return(gauss); +} + +//XXX: remove this +//typedef enum { +// PS_POLYNOMIAL_ORD, ///< Ordinary Polynomial +// PS_POLYNOMIAL_CHEB ///< Chebyshev Polynomial +//} psPolynomialType; /***************************************************************************** This routine must allocate memory for the polynomial structures. *****************************************************************************/ -psPolynomial1D* psPolynomial1DAlloc(int n) +psPolynomial1D* psPolynomial1DAlloc(int n, + psPolynomialType type) { int i = 0; @@ -151,4 +196,6 @@ newPoly = (psPolynomial1D* ) psAlloc(sizeof(psPolynomial1D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial1DFree); + + newPoly->type = type; newPoly->n = n; newPoly->coeff = (float *)psAlloc(n * sizeof(float)); @@ -161,8 +208,9 @@ } - return (newPoly); -} - -psPolynomial2D* psPolynomial2DAlloc(int nX, int nY) + return(newPoly); +} + +psPolynomial2D* psPolynomial2DAlloc(int nX, int nY, + psPolynomialType type) { int x = 0; @@ -172,4 +220,6 @@ newPoly = (psPolynomial2D* ) psAlloc(sizeof(psPolynomial2D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial2DFree); + + newPoly->type = type; newPoly->nX = nX; newPoly->nY = nY; @@ -191,8 +241,9 @@ } - return (newPoly); -} - -psPolynomial3D* psPolynomial3DAlloc(int nX, int nY, int nZ) + return(newPoly); +} + +psPolynomial3D* psPolynomial3DAlloc(int nX, int nY, int nZ, + psPolynomialType type) { int x = 0; @@ -203,4 +254,6 @@ newPoly = (psPolynomial3D* ) psAlloc(sizeof(psPolynomial3D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial3DFree); + + newPoly->type = type; newPoly->nX = nX; newPoly->nY = nY; @@ -230,8 +283,9 @@ } - return (newPoly); -} - -psPolynomial4D* psPolynomial4DAlloc(int nW, int nX, int nY, int nZ) + return(newPoly); +} + +psPolynomial4D* psPolynomial4DAlloc(int nW, int nX, int nY, int nZ, + psPolynomialType type) { int w = 0; @@ -243,4 +297,6 @@ newPoly = (psPolynomial4D* ) psAlloc(sizeof(psPolynomial4D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial4DFree); + + newPoly->type = type; newPoly->nW = nW; newPoly->nX = nX; @@ -278,5 +334,5 @@ } - return (newPoly); + return(newPoly); } @@ -359,5 +415,5 @@ XXX: Should the "coeffErr[]" should be used as well? *****************************************************************************/ -float psPolynomial1DEval(float x, const psPolynomial1D* myPoly) +float p_psOrgPolynomial1DEval(float x, const psPolynomial1D* myPoly) { int loop_x = 0; @@ -371,5 +427,5 @@ // NOTE: Do we want to flag this case? if (myPoly->n == 0) { - return (1.0); + return(1.0); } @@ -383,8 +439,67 @@ } - return (polySum); -} - -float psPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) + return(polySum); +} + +// XXX: You can do this without having to psAlloc() vector d. +float p_psChebPolynomial1DEval(float x, const psPolynomial1D* myPoly) +{ + psVector *d; + int n; + int i; + float tmp; + + n = myPoly->n; + d = psVectorAlloc(n, PS_TYPE_F32); + d->data.F32[n-1] = myPoly->coeff[n-1]; + d->data.F32[n-2] = (2.0 * x * d->data.F32[n-1]) + myPoly->coeff[n-2]; + for (i=n-3;i>=1;i--) { + d->data.F32[i] = (2.0 * x * d->data.F32[i+1]) - + (d->data.F32[i+2]) + + (myPoly->coeff[i]); + } + + tmp = (x * d->data.F32[1]) - + (d->data.F32[2]) + + (0.5 * myPoly->coeff[0]); + + psFree(d); + return(tmp); +} + +float psPolynomial1DEval(float x, const psPolynomial1D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psOrgPolynomial1DEval(x, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psChebPolynomial1DEval(x, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + +psVector *psPolynomial1DEvalVector(const psVector *x, + const psPolynomial1D *myPoly) +{ + psVector *tmp; + psVector *myX; + int i; + + PS_CONVERT_VECTOR_F32(x, myX); + + tmp = psVectorAlloc(x->n, PS_TYPE_F32); + for (i=0;in;i++) { + tmp->data.F32[i] = psPolynomial1DEval(x->data.F32[i], myPoly); + } + + if (x->type.type != PS_TYPE_F32) { + psFree(myX); + } + return(tmp); +} + + +float p_psOrgPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) { int loop_x = 0; @@ -403,8 +518,61 @@ } - return (polySum); -} - -float psPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) + return(polySum); +} + +float p_psChebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) +{ + return(0.0); +} + +float psPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psOrgPolynomial2DEval(x, y, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psChebPolynomial2DEval(x, y, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + + +psVector *psPolynomial2DEvalVector(const psVector *x, + const psVector *y, + const psPolynomial2D *myPoly) +{ + psVector *tmp; + psVector *myX; + psVector *myY; + int i; + int vecLen=x->n; + + PS_CONVERT_VECTOR_F32(x, myX); + PS_CONVERT_VECTOR_F32(y, myY); + vecLen=x->n; + if (y->n < vecLen) { + vecLen = y->n; + } + + tmp = psVectorAlloc(vecLen, PS_TYPE_F32); + for (i=0;idata.F32[i] = psPolynomial2DEval(myX->data.F32[i], + myY->data.F32[i], + myPoly); + } + + if (x->type.type != PS_TYPE_F32) { + psFree(myX); + } + if (y->type.type != PS_TYPE_F32) { + psFree(myY); + } + return(tmp); +} + + + +float p_psOrgPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) { int loop_x = 0; @@ -429,8 +597,73 @@ } - return (polySum); -} - -float psPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) + return(polySum); +} + +float p_psChebPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) +{ + return(0.0); +} + +float psPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psOrgPolynomial3DEval(x, y, z, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psChebPolynomial3DEval(x, y, z, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + +psVector *psPolynomial3DEvalVector(const psVector *x, + const psVector *y, + const psVector *z, + const psPolynomial3D *myPoly) +{ + psVector *tmp; + psVector *myX; + psVector *myY; + psVector *myZ; + int i; + int vecLen=x->n; + + PS_CONVERT_VECTOR_F32(x, myX); + PS_CONVERT_VECTOR_F32(y, myY); + PS_CONVERT_VECTOR_F32(z, myZ); + vecLen=x->n; + if (y->n < vecLen) { + vecLen = y->n; + } + if (z->n < vecLen) { + vecLen = z->n; + } + + tmp = psVectorAlloc(vecLen, PS_TYPE_F32); + for (i=0;idata.F32[i] = psPolynomial3DEval(myX->data.F32[i], + myY->data.F32[i], + myZ->data.F32[i], + myPoly); + } + + if (x->type.type != PS_TYPE_F32) { + psFree(myX); + } + if (y->type.type != PS_TYPE_F32) { + psFree(myY); + } + if (z->type.type != PS_TYPE_F32) { + psFree(myZ); + } + return(tmp); +} + + + + + + +float p_psOrgPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) { int loop_w = 0; @@ -461,8 +694,83 @@ } - return (polySum); -} - -psDPolynomial1D* psDPolynomial1DAlloc(int n) + return(polySum); +} + +float p_psChebPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) +{ + return(0.0); +} + +float psPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psOrgPolynomial4DEval(w,x,y,z, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psChebPolynomial4DEval(w,x,y,z, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + +psVector *psPolynomial4DEvalVector(const psVector *w, + const psVector *x, + const psVector *y, + const psVector *z, + const psPolynomial4D *myPoly) +{ + psVector *tmp; + psVector *myW; + psVector *myX; + psVector *myY; + psVector *myZ; + int i; + int vecLen=x->n; + + PS_CONVERT_VECTOR_F32(w, myW); + PS_CONVERT_VECTOR_F32(x, myX); + PS_CONVERT_VECTOR_F32(y, myY); + PS_CONVERT_VECTOR_F32(z, myZ); + vecLen=w->n; + if (y->n < vecLen) { + vecLen = y->n; + } + if (x->n < vecLen) { + vecLen = x->n; + } + if (z->n < vecLen) { + vecLen = z->n; + } + + tmp = psVectorAlloc(vecLen, PS_TYPE_F32); + for (i=0;idata.F32[i] = psPolynomial4DEval(myW->data.F32[i], + myX->data.F32[i], + myY->data.F32[i], + myZ->data.F32[i], + myPoly); + } + + if (w->type.type != PS_TYPE_F32) { + psFree(myW); + } + if (x->type.type != PS_TYPE_F32) { + psFree(myX); + } + if (y->type.type != PS_TYPE_F32) { + psFree(myY); + } + if (z->type.type != PS_TYPE_F32) { + psFree(myZ); + } + return(tmp); +} + + + + + +psDPolynomial1D* psDPolynomial1DAlloc(int n, + psPolynomialType type) { int i = 0; @@ -471,4 +779,6 @@ newPoly = (psDPolynomial1D* ) psAlloc(sizeof(psDPolynomial1D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial1DFree); + + newPoly->type = type; newPoly->n = n; newPoly->coeff = (double *)psAlloc(n * sizeof(double)); @@ -481,8 +791,9 @@ } - return (newPoly); -} - -psDPolynomial2D* psDPolynomial2DAlloc(int nX, int nY) + return(newPoly); +} + +psDPolynomial2D* psDPolynomial2DAlloc(int nX, int nY, + psPolynomialType type) { int x = 0; @@ -492,4 +803,6 @@ newPoly = (psDPolynomial2D* ) psAlloc(sizeof(psDPolynomial2D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial2DFree); + + newPoly->type = type; newPoly->nX = nX; newPoly->nY = nY; @@ -511,8 +824,9 @@ } - return (newPoly); -} - -psDPolynomial3D* psDPolynomial3DAlloc(int nX, int nY, int nZ) + return(newPoly); +} + +psDPolynomial3D* psDPolynomial3DAlloc(int nX, int nY, int nZ, + psPolynomialType type) { int x = 0; @@ -523,4 +837,6 @@ newPoly = (psDPolynomial3D* ) psAlloc(sizeof(psDPolynomial3D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial3DFree); + + newPoly->type = type; newPoly->nX = nX; newPoly->nY = nY; @@ -550,8 +866,9 @@ } - return (newPoly); -} - -psDPolynomial4D* psDPolynomial4DAlloc(int nW, int nX, int nY, int nZ) + return(newPoly); +} + +psDPolynomial4D* psDPolynomial4DAlloc(int nW, int nX, int nY, int nZ, + psPolynomialType type) { int w = 0; @@ -563,4 +880,6 @@ newPoly = (psDPolynomial4D* ) psAlloc(sizeof(psDPolynomial4D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial4DFree); + + newPoly->type = type; newPoly->nW = nW; newPoly->nX = nX; @@ -598,5 +917,5 @@ } - return (newPoly); + return(newPoly); } @@ -673,5 +992,5 @@ Polynomial coefficients will be accessed in [w][x][y][z] fashion. *****************************************************************************/ -double psDPolynomial1DEval(double x, const psDPolynomial1D* myPoly) +double p_psDOrgPolynomial1DEval(double x, const psDPolynomial1D* myPoly) { int loop_x = 0; @@ -681,5 +1000,5 @@ // NOTE: Do we want to flag this case? if (myPoly->n == 0) { - return (1.0); + return(1.0); } @@ -689,8 +1008,47 @@ } - return (polySum); -} - -double psDPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) + return(polySum); +} + +double p_psDChebPolynomial1DEval(double x, const psDPolynomial1D* myPoly) +{ + return(0.0); +} + +double psDPolynomial1DEval(double x, const psDPolynomial1D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psDOrgPolynomial1DEval(x, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psDChebPolynomial1DEval(x, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + +psVector *psDPolynomial1DEvalVector(const psVector *x, + const psDPolynomial1D *myPoly) +{ + psVector *tmp; + psVector *myX; + int i; + + PS_CONVERT_VECTOR_F64(x, myX); + + tmp = psVectorAlloc(x->n, PS_TYPE_F64); + for (i=0;in;i++) { + tmp->data.F64[i] = psDPolynomial1DEval(x->data.F64[i], myPoly); + } + + if (x->type.type != PS_TYPE_F64) { + psFree(myX); + } + return(tmp); +} + + + +double p_psDOrgPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) { int loop_x = 0; @@ -709,8 +1067,60 @@ } - return (polySum); -} - -double psDPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) + return(polySum); +} + +double p_psDChebPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) +{ + return(0.0); +} + +double psDPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psDOrgPolynomial2DEval(x, y, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psDChebPolynomial2DEval(x, y, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + +psVector *psDPolynomial2DEvalVector(const psVector *x, + const psVector *y, + const psDPolynomial2D *myPoly) +{ + psVector *tmp; + psVector *myX; + psVector *myY; + int i; + int vecLen=x->n; + + PS_CONVERT_VECTOR_F64(x, myX); + PS_CONVERT_VECTOR_F64(y, myY); + vecLen=x->n; + if (y->n < vecLen) { + vecLen = y->n; + } + + tmp = psVectorAlloc(vecLen, PS_TYPE_F64); + for (i=0;idata.F64[i] = psDPolynomial2DEval(myX->data.F64[i], + myY->data.F64[i], + myPoly); + } + + if (x->type.type != PS_TYPE_F64) { + psFree(myX); + } + if (y->type.type != PS_TYPE_F64) { + psFree(myY); + } + return(tmp); +} + + + +double p_psDOrgPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) { int loop_x = 0; @@ -735,8 +1145,75 @@ } - return (polySum); -} - -double psDPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) + return(polySum); +} + +double p_psDChebPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) +{ + return(0.0); +} + +double psDPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psDOrgPolynomial3DEval(x, y, z, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psDChebPolynomial3DEval(x, y, z, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + +psVector *psDPolynomial3DEvalVector(const psVector *x, + const psVector *y, + const psVector *z, + const psDPolynomial3D *myPoly) +{ + psVector *tmp; + psVector *myX; + psVector *myY; + psVector *myZ; + int i; + int vecLen=x->n; + + PS_CONVERT_VECTOR_F64(x, myX); + PS_CONVERT_VECTOR_F64(y, myY); + PS_CONVERT_VECTOR_F64(z, myZ); + vecLen=x->n; + if (y->n < vecLen) { + vecLen = y->n; + } + if (z->n < vecLen) { + vecLen = z->n; + } + + tmp = psVectorAlloc(vecLen, PS_TYPE_F64); + for (i=0;idata.F64[i] = psDPolynomial3DEval(myX->data.F64[i], + myY->data.F64[i], + myZ->data.F64[i], + myPoly); + } + + if (x->type.type != PS_TYPE_F64) { + psFree(myX); + } + if (y->type.type != PS_TYPE_F64) { + psFree(myY); + } + if (z->type.type != PS_TYPE_F64) { + psFree(myZ); + } + return(tmp); +} + + + + + + + + +double p_psDOrgPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) { int loop_w = 0; @@ -767,4 +1244,222 @@ } - return (polySum); -} + return(polySum); +} + +double p_psDChebPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) +{ + return(0.0); +} + +double psDPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psDOrgPolynomial4DEval(w,x,y,z, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psDChebPolynomial4DEval(w,x,y,z, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + +psVector *psDPolynomial4DEvalVector(const psVector *w, + const psVector *x, + const psVector *y, + const psVector *z, + const psDPolynomial4D *myPoly) +{ + psVector *tmp; + psVector *myW; + psVector *myX; + psVector *myY; + psVector *myZ; + int i; + int vecLen=x->n; + + PS_CONVERT_VECTOR_F64(w, myW); + PS_CONVERT_VECTOR_F64(x, myX); + PS_CONVERT_VECTOR_F64(y, myY); + PS_CONVERT_VECTOR_F64(z, myZ); + vecLen=w->n; + if (y->n < vecLen) { + vecLen = y->n; + } + if (x->n < vecLen) { + vecLen = x->n; + } + if (z->n < vecLen) { + vecLen = z->n; + } + + tmp = psVectorAlloc(vecLen, PS_TYPE_F64); + for (i=0;idata.F64[i] = psDPolynomial4DEval(myW->data.F64[i], + myX->data.F64[i], + myY->data.F64[i], + myZ->data.F64[i], + myPoly); + } + + if (w->type.type != PS_TYPE_F64) { + psFree(myW); + } + if (x->type.type != PS_TYPE_F64) { + psFree(myX); + } + if (y->type.type != PS_TYPE_F64) { + psFree(myY); + } + if (z->type.type != PS_TYPE_F64) { + psFree(myZ); + } + return(tmp); +} + + + + +//typedef struct { +// int n; +// psPolynomial1D **spline; +// float *domains; +//} psSpline1D; + +/***************************************************************************** + + NOTE: "n" specifies the number of spline polynomials. Therefore, there + must exist n+1 points in "domains". + *****************************************************************************/ +psSpline1D *psSpline1DAlloc(int numSplines, + int order, + float min, + float max) +{ + psSpline1D *tmp = NULL; + int i; + float tmpDomain; + float width; + + tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D)); + + tmp->n = numSplines; + + tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); + for (i=0;ispline)[i] = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD); + } + + tmp->domains = (float *) psAlloc((numSplines+1) * sizeof(float)); + width = (max - min) / ((float) numSplines); + tmpDomain = min; + for (i=0;idomains)[i] = tmpDomain; + tmpDomain+= width; + } + + return(tmp); +} + +psSpline1D *psSpline1DAllocGeneric(const psVector *bounds, + int order) +{ + psSpline1D *tmp = NULL; + int i; + int numSplines; + + tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D)); + + numSplines = bounds->n - 1; + + tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); + for (i=0;ispline)[i] = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD); + } + + tmp->domains = (float *) psAlloc((bounds->n) * sizeof(float)); + for (i=0;in;i++) { + (tmp->domains)[i] = bounds->data.F32[i]; + } + + return(tmp); +} + +/***************************************************************************** +VectorBinDisect(): This is a private function which takes as input a vector +of floating point data as well as a single floating point values. The input +vector values are assumed to be non-decreasing (v[i-1] <= v[j] for all j>=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 -1. + *****************************************************************************/ +int VectorBinDisect(float *bins, + int numBins, + float x) +{ + int min; + int max; + int mid; + + if ((x < bins[0]) || + (x > bins[numBins-1])) { + psLogMsg(__func__, PS_LOG_WARN, + "VectorBinDisect(): ordinate %f is outside vector range (%f - %f).", + x, bins[0], bins[numBins-1]); + return(-1); + } + + min = 0; + max = numBins-1; + mid = (max-min)/2; + + while (min != max) { + mid = (max-min)/2; + + if (x < bins[mid]) { + max = mid; + } else { + min = mid; + } + } + + return(min); +} + +float psSpline1DEval(const psSpline1D *spline, + float x) +{ + int binNum; + int n; + + n = spline->n; + binNum = VectorBinDisect(spline->domains, (spline->n)+1, x); + if (binNum == -1) { + psLogMsg(__func__, PS_LOG_WARN, + "psSpline1DEval(): x ordinate (%f) is outside the spline range (%f - %f).", + x, (spline->domains)[0], + (spline->domains)[n-1]); + + if (x < (spline->domains)[0]) { + return(psPolynomial1DEval(x, spline->spline[0])); + } else if (x > (spline->domains)[n-1]) { + return(psPolynomial1DEval(x, spline->spline[n-1])); + } + } + + return(psPolynomial1DEval(x, spline->spline[binNum])); +} + +psVector *psSpline1DEvalVector(const psVector *x, + const psSpline1D *spline) +{ + int i; + psVector *tmpVector; + + tmpVector = psVectorAlloc(x->n, PS_TYPE_F32); + for (i=0;in;i++) { + tmpVector->data.F32[i] = psSpline1DEval(spline, x->data.F32[i]); + } + + return(tmpVector); +} +