Index: trunk/psLib/src/math/psPolynomial.c =================================================================== --- trunk/psLib/src/math/psPolynomial.c (revision 2224) +++ trunk/psLib/src/math/psPolynomial.c (revision 2273) @@ -7,6 +7,6 @@ * polynomials. It also contains a Gaussian functions. * - * @version $Revision: 1.58 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-28 00:22:53 $ + * @version $Revision: 1.59 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-11-04 01:04:59 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -35,4 +35,7 @@ #include "psFunctions.h" #include "psConstants.h" + +#include "psDataManipErrors.h" + /*****************************************************************************/ /* DEFINE STATEMENTS */ @@ -50,4 +53,7 @@ static void dPolynomial3DFree(psDPolynomial3D* myPoly); static void dPolynomial4DFree(psDPolynomial4D* myPoly); +static void spline1DFree(psSpline1D *tmpSpline); +static psS32 vectorBinDisectF32(float *bins,psS32 numBins,float x); +static psS32 vectorBinDisectS32(psS32 *bins,psS32 numBins,psS32 x); /*****************************************************************************/ @@ -67,6 +73,167 @@ /*****************************************************************************/ +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->domains); + + return; +} + +static void polynomial1DFree(psPolynomial1D* myPoly) +{ + psFree(myPoly->coeff); + psFree(myPoly->coeffErr); + psFree(myPoly->mask); +} + +static void polynomial2DFree(psPolynomial2D* myPoly) +{ + psS32 x = 0; + + for (x = 0; x < myPoly->nX; x++) { + psFree(myPoly->coeff[x]); + psFree(myPoly->coeffErr[x]); + psFree(myPoly->mask[x]); + } + psFree(myPoly->coeff); + psFree(myPoly->coeffErr); + psFree(myPoly->mask); +} + +static void polynomial3DFree(psPolynomial3D* myPoly) +{ + psS32 x = 0; + psS32 y = 0; + + for (x = 0; x < myPoly->nX; x++) { + for (y = 0; y < myPoly->nY; y++) { + psFree(myPoly->coeff[x][y]); + psFree(myPoly->coeffErr[x][y]); + psFree(myPoly->mask[x][y]); + } + psFree(myPoly->coeff[x]); + psFree(myPoly->coeffErr[x]); + psFree(myPoly->mask[x]); + } + + psFree(myPoly->coeff); + psFree(myPoly->coeffErr); + psFree(myPoly->mask); +} + +static void polynomial4DFree(psPolynomial4D* myPoly) +{ + psS32 w = 0; + psS32 x = 0; + psS32 y = 0; + + for (w = 0; w < myPoly->nW; w++) { + for (x = 0; x < myPoly->nX; x++) { + for (y = 0; y < myPoly->nY; y++) { + psFree(myPoly->coeff[w][x][y]); + psFree(myPoly->coeffErr[w][x][y]); + psFree(myPoly->mask[w][x][y]); + } + psFree(myPoly->coeff[w][x]); + psFree(myPoly->coeffErr[w][x]); + psFree(myPoly->mask[w][x]); + } + psFree(myPoly->coeff[w]); + psFree(myPoly->coeffErr[w]); + psFree(myPoly->mask[w]); + } + + psFree(myPoly->coeff); + psFree(myPoly->coeffErr); + psFree(myPoly->mask); +} + +static void dPolynomial1DFree(psDPolynomial1D* myPoly) +{ + psFree(myPoly->coeff); + psFree(myPoly->coeffErr); + psFree(myPoly->mask); +} + +static void dPolynomial2DFree(psDPolynomial2D* myPoly) +{ + psS32 x = 0; + + for (x = 0; x < myPoly->nX; x++) { + psFree(myPoly->coeff[x]); + psFree(myPoly->coeffErr[x]); + psFree(myPoly->mask[x]); + } + psFree(myPoly->coeff); + psFree(myPoly->coeffErr); + psFree(myPoly->mask); +} + +static void dPolynomial3DFree(psDPolynomial3D* myPoly) +{ + psS32 x = 0; + psS32 y = 0; + + for (x = 0; x < myPoly->nX; x++) { + for (y = 0; y < myPoly->nY; y++) { + psFree(myPoly->coeff[x][y]); + psFree(myPoly->coeffErr[x][y]); + psFree(myPoly->mask[x][y]); + } + psFree(myPoly->coeff[x]); + psFree(myPoly->coeffErr[x]); + psFree(myPoly->mask[x]); + } + + psFree(myPoly->coeff); + psFree(myPoly->coeffErr); + psFree(myPoly->mask); +} + +static void dPolynomial4DFree(psDPolynomial4D* myPoly) +{ + psS32 w = 0; + psS32 x = 0; + psS32 y = 0; + + for (w = 0; w < myPoly->nW; w++) { + for (x = 0; x < myPoly->nX; x++) { + for (y = 0; y < myPoly->nY; y++) { + psFree(myPoly->coeff[w][x][y]); + psFree(myPoly->coeffErr[w][x][y]); + psFree(myPoly->mask[w][x][y]); + } + psFree(myPoly->coeff[w][x]); + psFree(myPoly->coeffErr[w][x]); + psFree(myPoly->mask[w][x]); + } + psFree(myPoly->coeff[w]); + psFree(myPoly->coeffErr[w]); + psFree(myPoly->mask[w]); + } + + psFree(myPoly->coeff); + psFree(myPoly->coeffErr); + psFree(myPoly->mask); +} + /***************************************************************************** -CreateChebyshevPolys(n): this routine takes as input the required order n, +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. @@ -76,5 +243,5 @@ outer coefficients of the Chebyshev polynomials. *****************************************************************************/ -static psPolynomial1D **CreateChebyshevPolys(psS32 maxChebyPoly) +static psPolynomial1D **createChebyshevPolys(psS32 maxChebyPoly) { PS_INT_CHECK_NON_NEGATIVE(maxChebyPoly, NULL); @@ -103,4 +270,655 @@ return (chebPolys); +} + +/***************************************************************************** + Polynomial coefficients will be accessed in [w][x][y][z] fashion. + + XXX: Should the "coeffErr[]" should be used as well? + *****************************************************************************/ +static float ordPolynomial1DEval(float x, const psPolynomial1D* myPoly) +{ + psS32 loop_x = 0; + float polySum = 0.0; + float xSum = 1.0; + + psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4, + "---- Calling ordPolynomial1DEval(%f)\n", x); + psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4, + "Polynomial order is %d\n", myPoly->n); + for (loop_x = 0; loop_x < myPoly->n; loop_x++) { + psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4, + "Polynomial coeff[%d] is %f\n", loop_x, myPoly->coeff[loop_x]); + } + + for (loop_x = 0; loop_x < myPoly->n; loop_x++) { + if (myPoly->mask[loop_x] == 0) { + psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 10, + "polysum+= sum*coeff [%f+= (%f * %f)\n", polySum, xSum, myPoly->coeff[loop_x]); + polySum += xSum * myPoly->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? +static float chebPolynomial1DEval(float x, const psPolynomial1D* myPoly) +{ + psVector *d; + psS32 n; + psS32 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); + /* + + psS32 n; + psS32 i; + float tmp; + psPolynomial1D **chebPolys = NULL; + + n = myPoly->n; + chebPolys = createChebyshevPolys(n); + + tmp = 0.0; + for (i=0;in;i++) { + tmp+= (myPoly->coeff[i] * psPolynomial1DEval(x, chebPolys[i])); + // printf("HMMM: psPolynomial1DEval(%f, chebPolys[%d]) is %f\n", x, i, psPolynomial1DEval(x, chebPolys[i])); + } + tmp-= (myPoly->coeff[0]/2.0); + + + return(tmp); + */ +} + +static float ordPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) +{ + PS_POLY_CHECK_NULL(myPoly, NAN); + + psS32 loop_x = 0; + psS32 loop_y = 0; + float polySum = 0.0; + float xSum = 1.0; + float ySum = 1.0; + + for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { + ySum = xSum; + for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { + if (myPoly->mask[loop_x][loop_y] == 0) { + polySum += ySum * myPoly->coeff[loop_x][loop_y]; + ySum *= y; + } + } + xSum *= x; + } + + return(polySum); +} + +static float chebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) +{ + PS_POLY_CHECK_NULL(myPoly, NAN); + + psS32 loop_x = 0; + psS32 loop_y = 0; + psS32 i = 0; + float polySum = 0.0; + psPolynomial1D* *chebPolys = NULL; + psS32 maxChebyPoly = 0; + + // Determine how many Chebyshev polynomials + // are needed, then create them. + maxChebyPoly = myPoly->nX; + if (myPoly->nY > maxChebyPoly) { + maxChebyPoly = myPoly->nY; + } + chebPolys = createChebyshevPolys(maxChebyPoly); + + for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { + for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { + if (myPoly->mask[loop_x][loop_y] == 0) { + polySum += myPoly->coeff[loop_x][loop_y] * + psPolynomial1DEval(x, chebPolys[loop_x]) * + psPolynomial1DEval(y, chebPolys[loop_y]); + } + } + } + for (i=0;inX; loop_x++) { + ySum = xSum; + for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { + zSum = ySum; + for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { + if (myPoly->mask[loop_x][loop_y][loop_z] == 0) { + polySum += zSum * myPoly->coeff[loop_x][loop_y][loop_z]; + zSum *= z; + } + } + ySum *= y; + } + xSum *= x; + } + + return(polySum); +} + +static float chebPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) +{ + psS32 loop_x = 0; + psS32 loop_y = 0; + psS32 loop_z = 0; + psS32 i = 0; + float polySum = 0.0; + psPolynomial1D* *chebPolys = NULL; + psS32 maxChebyPoly = 0; + + // Determine how many Chebyshev polynomials + // are needed, then create them. + maxChebyPoly = myPoly->nX; + if (myPoly->nY > maxChebyPoly) { + maxChebyPoly = myPoly->nY; + } + if (myPoly->nZ > maxChebyPoly) { + maxChebyPoly = myPoly->nZ; + } + chebPolys = createChebyshevPolys(maxChebyPoly); + + for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { + for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { + for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { + if (myPoly->mask[loop_x][loop_y][loop_z] == 0) { + polySum += myPoly->coeff[loop_x][loop_y][loop_z] * + psPolynomial1DEval(x, chebPolys[loop_x]) * + psPolynomial1DEval(y, chebPolys[loop_y]) * + psPolynomial1DEval(z, chebPolys[loop_z]); + } + } + } + } + + for (i=0;inW; loop_w++) { + xSum = wSum; + for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { + ySum = xSum; + for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { + zSum = ySum; + for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { + if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) { + polySum += zSum * myPoly->coeff[loop_w][loop_x][loop_y][loop_z]; + zSum *= z; + } + } + ySum *= y; + } + xSum *= x; + } + wSum *= w; + } + + return(polySum); +} + +static float chebPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) +{ + psS32 loop_w = 0; + psS32 loop_x = 0; + psS32 loop_y = 0; + psS32 loop_z = 0; + psS32 i = 0; + float polySum = 0.0; + psPolynomial1D* *chebPolys = NULL; + psS32 maxChebyPoly = 0; + + // Determine how many Chebyshev polynomials + // are needed, then create them. + maxChebyPoly = myPoly->nW; + if (myPoly->nX > maxChebyPoly) { + maxChebyPoly = myPoly->nX; + } + if (myPoly->nY > maxChebyPoly) { + maxChebyPoly = myPoly->nY; + } + if (myPoly->nZ > maxChebyPoly) { + maxChebyPoly = myPoly->nZ; + } + chebPolys = createChebyshevPolys(maxChebyPoly); + + for (loop_w = 0; loop_w < myPoly->nW; loop_w++) { + for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { + for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { + for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { + if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) { + polySum += myPoly->coeff[loop_w][loop_x][loop_y][loop_z] * + psPolynomial1DEval(w, chebPolys[loop_w]) * + psPolynomial1DEval(x, chebPolys[loop_x]) * + psPolynomial1DEval(y, chebPolys[loop_y]) * + psPolynomial1DEval(z, chebPolys[loop_z]); + } + } + } + } + } + + for (i=0;in; loop_x++) { + if (myPoly->mask[loop_x] == 0) { + polySum += xSum * myPoly->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? +static double dChebPolynomial1DEval(double x, const psDPolynomial1D* myPoly) +{ + psVector *d; + psS32 n; + psS32 i; + double tmp; + + n = myPoly->n; + d = psVectorAlloc(n, PS_TYPE_F64); + d->data.F64[n-1] = myPoly->coeff[n-1]; + d->data.F64[n-2] = (2.0 * x * d->data.F64[n-1]) + myPoly->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]) + + (myPoly->coeff[i]); + } + + tmp = (x * d->data.F64[1]) - + (d->data.F64[2]) + + (0.5 * myPoly->coeff[0]); + + psFree(d); + return(tmp); +} + +static double dOrdPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) +{ + psS32 loop_x = 0; + psS32 loop_y = 0; + double polySum = 0.0; + double xSum = 1.0; + double ySum = 1.0; + + for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { + ySum = xSum; + for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { + if (myPoly->mask[loop_x][loop_y] == 0) { + polySum += ySum * myPoly->coeff[loop_x][loop_y]; + ySum *= y; + } + } + xSum *= x; + } + + return(polySum); +} + +static double dChebPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) +{ + psS32 loop_x = 0; + psS32 loop_y = 0; + psS32 i = 0; + double polySum = 0.0; + psPolynomial1D* *chebPolys = NULL; + psS32 maxChebyPoly = 0; + + // Determine how many Chebyshev polynomials + // are needed, then create them. + maxChebyPoly = myPoly->nX; + if (myPoly->nY > maxChebyPoly) { + maxChebyPoly = myPoly->nY; + } + chebPolys = createChebyshevPolys(maxChebyPoly); + + for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { + for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { + if (myPoly->mask[loop_x][loop_y] == 0) { + polySum += myPoly->coeff[loop_x][loop_y] * + psPolynomial1DEval(x, chebPolys[loop_x]) * + psPolynomial1DEval(y, chebPolys[loop_y]); + } + } + } + + for (i=0;inX; loop_x++) { + ySum = xSum; + for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { + zSum = ySum; + for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { + if (myPoly->mask[loop_x][loop_y][loop_z] == 0) { + polySum += zSum * myPoly->coeff[loop_x][loop_y][loop_z]; + zSum *= z; + } + } + ySum *= y; + } + xSum *= x; + } + + return(polySum); +} + +static double dChebPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) +{ + psS32 loop_x = 0; + psS32 loop_y = 0; + psS32 loop_z = 0; + psS32 i = 0; + double polySum = 0.0; + psPolynomial1D* *chebPolys = NULL; + psS32 maxChebyPoly = 0; + + // Determine how many Chebyshev polynomials + // are needed, then create them. + maxChebyPoly = myPoly->nX; + if (myPoly->nY > maxChebyPoly) { + maxChebyPoly = myPoly->nY; + } + if (myPoly->nZ > maxChebyPoly) { + maxChebyPoly = myPoly->nZ; + } + chebPolys = createChebyshevPolys(maxChebyPoly); + + for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { + for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { + for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { + if (myPoly->mask[loop_x][loop_y][loop_z] == 0) { + polySum += myPoly->coeff[loop_x][loop_y][loop_z] * + psPolynomial1DEval(x, chebPolys[loop_x]) * + psPolynomial1DEval(y, chebPolys[loop_y]) * + psPolynomial1DEval(z, chebPolys[loop_z]); + } + } + } + } + + for (i=0;inW; loop_w++) { + xSum = wSum; + for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { + ySum = xSum; + for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { + zSum = ySum; + for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { + if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) { + polySum += zSum * myPoly->coeff[loop_w][loop_x][loop_y][loop_z]; + zSum *= z; + } + } + ySum *= y; + } + xSum *= x; + } + wSum *= w; + } + + return(polySum); +} + +static double dChebPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) +{ + psS32 loop_w = 0; + psS32 loop_x = 0; + psS32 loop_y = 0; + psS32 loop_z = 0; + psS32 i = 0; + double polySum = 0.0; + psPolynomial1D* *chebPolys = NULL; + psS32 maxChebyPoly = 0; + + // Determine how many Chebyshev polynomials + // are needed, then create them. + maxChebyPoly = myPoly->nW; + if (myPoly->nX > maxChebyPoly) { + maxChebyPoly = myPoly->nX; + } + if (myPoly->nY > maxChebyPoly) { + maxChebyPoly = myPoly->nY; + } + if (myPoly->nZ > maxChebyPoly) { + maxChebyPoly = myPoly->nZ; + } + chebPolys = createChebyshevPolys(maxChebyPoly); + + for (loop_w = 0; loop_w < myPoly->nW; loop_w++) { + for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { + for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { + for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { + if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) { + polySum += myPoly->coeff[loop_w][loop_x][loop_y][loop_z] * + psPolynomial1DEval(w, chebPolys[loop_w]) * + psPolynomial1DEval(x, chebPolys[loop_x]) * + psPolynomial1DEval(y, chebPolys[loop_y]) * + psPolynomial1DEval(z, chebPolys[loop_z]); + } + } + } + } + } + + for (i=0;idata.F32, true); + /* + psVector *p = psVectorAlloc(n, PS_TYPE_F32); + float tmp; + */ + + psTrace(".psLib.dataManip.psFunctions.fullInterpolate1DF32", 4, + "---- fullInterpolate1DF32() begin (%d-order at x=%f) (%d data points)----\n", n-1, x, n); + + for (i=0;idata.F32[i] = range[i]; + psTrace(".psLib.dataManip.psFunctions.fullInterpolate1DF32", 6, + "p->data.F32[%d] is %f\n", i, p->data.F32[i]); + + } + + // From NR, during each iteration of the m loop, we are computing the + // p_{i ... i+m} terms. + for (m=1;mdata.F32[i] = (((x-domain[i+m]) * p->data.F32[i]) + + ((domain[i]-x) * p->data.F32[i+1])) / + (domain[i] - domain[i+m]); + //printf("((%f-%f * %f) + (%f-%f * %f)) / (%f - %f)\n", x, domain[i+m], p->data.F32[i], domain[i], x, p->data.F32[i+1], domain[i], domain[i+m]); + psTrace(".psLib.dataManip.psFunctions.fullInterpolate1DF32", 6, + "p->data.F32[%d] is %f\n", i, p->data.F32[i]); + } + } + psTrace(".psLib.dataManip.psFunctions.fullInterpolate1DF32", 4, + "---- fullInterpolate1DF32() end ----\n"); + + /* + tmp = p->data.F32[0]; + psFree(p); + return(tmp); + */ + return(p->data.F32[0]); +} + + +/***************************************************************************** +interpolate1DF32(): this is the base 1-D flat memory routine to perform +LaGrange interpolation. + *****************************************************************************/ +static float interpolate1DF32(float *domain, + float *range, + psS32 n, + psS32 order, + float 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)); } @@ -336,152 +1154,4 @@ } -static void polynomial1DFree(psPolynomial1D* myPoly) -{ - psFree(myPoly->coeff); - psFree(myPoly->coeffErr); - psFree(myPoly->mask); -} - -static void polynomial2DFree(psPolynomial2D* myPoly) -{ - psS32 x = 0; - - for (x = 0; x < myPoly->nX; x++) { - psFree(myPoly->coeff[x]); - psFree(myPoly->coeffErr[x]); - psFree(myPoly->mask[x]); - } - psFree(myPoly->coeff); - psFree(myPoly->coeffErr); - psFree(myPoly->mask); -} - -static void polynomial3DFree(psPolynomial3D* myPoly) -{ - psS32 x = 0; - psS32 y = 0; - - for (x = 0; x < myPoly->nX; x++) { - for (y = 0; y < myPoly->nY; y++) { - psFree(myPoly->coeff[x][y]); - psFree(myPoly->coeffErr[x][y]); - psFree(myPoly->mask[x][y]); - } - psFree(myPoly->coeff[x]); - psFree(myPoly->coeffErr[x]); - psFree(myPoly->mask[x]); - } - - psFree(myPoly->coeff); - psFree(myPoly->coeffErr); - psFree(myPoly->mask); -} - -static void polynomial4DFree(psPolynomial4D* myPoly) -{ - psS32 w = 0; - psS32 x = 0; - psS32 y = 0; - - for (w = 0; w < myPoly->nW; w++) { - for (x = 0; x < myPoly->nX; x++) { - for (y = 0; y < myPoly->nY; y++) { - psFree(myPoly->coeff[w][x][y]); - psFree(myPoly->coeffErr[w][x][y]); - psFree(myPoly->mask[w][x][y]); - } - psFree(myPoly->coeff[w][x]); - psFree(myPoly->coeffErr[w][x]); - psFree(myPoly->mask[w][x]); - } - psFree(myPoly->coeff[w]); - psFree(myPoly->coeffErr[w]); - psFree(myPoly->mask[w]); - } - - psFree(myPoly->coeff); - psFree(myPoly->coeffErr); - psFree(myPoly->mask); -} - -/***************************************************************************** - Polynomial coefficients will be accessed in [w][x][y][z] fashion. - - XXX: Should the "coeffErr[]" should be used as well? - *****************************************************************************/ -float p_psOrdPolynomial1DEval(float x, const psPolynomial1D* myPoly) -{ - psS32 loop_x = 0; - float polySum = 0.0; - float xSum = 1.0; - - psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 4, - "---- Calling p_psOrdPolynomial1DEval(%f)\n", x); - psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 4, - "Polynomial order is %d\n", myPoly->n); - for (loop_x = 0; loop_x < myPoly->n; loop_x++) { - psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 4, - "Polynomial coeff[%d] is %f\n", loop_x, myPoly->coeff[loop_x]); - } - - for (loop_x = 0; loop_x < myPoly->n; loop_x++) { - if (myPoly->mask[loop_x] == 0) { - psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 10, - "polysum+= sum*coeff [%f+= (%f * %f)\n", polySum, xSum, myPoly->coeff[loop_x]); - polySum += xSum * myPoly->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? -float p_psChebPolynomial1DEval(float x, const psPolynomial1D* myPoly) -{ - psVector *d; - psS32 n; - psS32 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); - /* - - psS32 n; - psS32 i; - float tmp; - psPolynomial1D **chebPolys = NULL; - - n = myPoly->n; - chebPolys = CreateChebyshevPolys(n); - - tmp = 0.0; - for (i=0;in;i++) { - tmp+= (myPoly->coeff[i] * psPolynomial1DEval(x, chebPolys[i])); - // printf("HMMM: psPolynomial1DEval(%f, chebPolys[%d]) is %f\n", x, i, psPolynomial1DEval(x, chebPolys[i])); - } - tmp-= (myPoly->coeff[0]/2.0); - - - return(tmp); - */ -} - float psPolynomial1DEval(float x, const psPolynomial1D* myPoly) { @@ -489,9 +1159,11 @@ if (myPoly->type == PS_POLYNOMIAL_ORD) { - return(p_psOrdPolynomial1DEval(x, myPoly)); + return(ordPolynomial1DEval(x, myPoly)); } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { - return(p_psChebPolynomial1DEval(x, myPoly)); + return(chebPolynomial1DEval(x, myPoly)); } else { - psError(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, + myPoly->type); } return(0.0); @@ -520,78 +1192,19 @@ } - -float p_psOrdPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) +float psPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) { PS_POLY_CHECK_NULL(myPoly, NAN); - psS32 loop_x = 0; - psS32 loop_y = 0; - float polySum = 0.0; - float xSum = 1.0; - float ySum = 1.0; - - for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { - ySum = xSum; - for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { - if (myPoly->mask[loop_x][loop_y] == 0) { - polySum += ySum * myPoly->coeff[loop_x][loop_y]; - ySum *= y; - } - } - xSum *= x; - } - - return(polySum); -} - -float p_psChebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) -{ - PS_POLY_CHECK_NULL(myPoly, NAN); - - psS32 loop_x = 0; - psS32 loop_y = 0; - psS32 i = 0; - float polySum = 0.0; - psPolynomial1D* *chebPolys = NULL; - psS32 maxChebyPoly = 0; - - // Determine how many Chebyshev polynomials - // are needed, then create them. - maxChebyPoly = myPoly->nX; - if (myPoly->nY > maxChebyPoly) { - maxChebyPoly = myPoly->nY; - } - chebPolys = CreateChebyshevPolys(maxChebyPoly); - - for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { - for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { - if (myPoly->mask[loop_x][loop_y] == 0) { - polySum += myPoly->coeff[loop_x][loop_y] * - psPolynomial1DEval(x, chebPolys[loop_x]) * - psPolynomial1DEval(y, chebPolys[loop_y]); - } - } - } - for (i=0;itype == PS_POLYNOMIAL_ORD) { - return(p_psOrdPolynomial2DEval(x, y, myPoly)); + return(ordPolynomial2DEval(x, y, myPoly)); } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { - return(p_psChebPolynomial2DEval(x, y, myPoly)); + return(chebPolynomial2DEval(x, y, myPoly)); } else { - psError(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, + myPoly->type); } return(0.0); } - psVector *psPolynomial2DEvalVector(const psVector *x, @@ -632,75 +1245,4 @@ } - - -float p_psOrdPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) -{ - psS32 loop_x = 0; - psS32 loop_y = 0; - psS32 loop_z = 0; - float polySum = 0.0; - float xSum = 1.0; - float ySum = 1.0; - float zSum = 1.0; - - for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { - ySum = xSum; - for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { - zSum = ySum; - for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_x][loop_y][loop_z] == 0) { - polySum += zSum * myPoly->coeff[loop_x][loop_y][loop_z]; - zSum *= z; - } - } - ySum *= y; - } - xSum *= x; - } - - return(polySum); -} - -float p_psChebPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) -{ - psS32 loop_x = 0; - psS32 loop_y = 0; - psS32 loop_z = 0; - psS32 i = 0; - float polySum = 0.0; - psPolynomial1D* *chebPolys = NULL; - psS32 maxChebyPoly = 0; - - // Determine how many Chebyshev polynomials - // are needed, then create them. - maxChebyPoly = myPoly->nX; - if (myPoly->nY > maxChebyPoly) { - maxChebyPoly = myPoly->nY; - } - if (myPoly->nZ > maxChebyPoly) { - maxChebyPoly = myPoly->nZ; - } - chebPolys = CreateChebyshevPolys(maxChebyPoly); - - for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { - for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { - for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_x][loop_y][loop_z] == 0) { - polySum += myPoly->coeff[loop_x][loop_y][loop_z] * - psPolynomial1DEval(x, chebPolys[loop_x]) * - psPolynomial1DEval(y, chebPolys[loop_y]) * - psPolynomial1DEval(z, chebPolys[loop_z]); - } - } - } - } - - for (i=0;itype == PS_POLYNOMIAL_ORD) { - return(p_psOrdPolynomial3DEval(x, y, z, myPoly)); + return(ordPolynomial3DEval(x, y, z, myPoly)); } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { - return(p_psChebPolynomial3DEval(x, y, z, myPoly)); + return(chebPolynomial3DEval(x, y, z, myPoly)); } else { - psError(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, + myPoly->type); } return(0.0); @@ -765,91 +1309,4 @@ } - - - - - -float p_psOrdPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) -{ - psS32 loop_w = 0; - psS32 loop_x = 0; - psS32 loop_y = 0; - psS32 loop_z = 0; - float polySum = 0.0; - float wSum = 1.0; - float xSum = 1.0; - float ySum = 1.0; - float zSum = 1.0; - - for (loop_w = 0; loop_w < myPoly->nW; loop_w++) { - xSum = wSum; - for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { - ySum = xSum; - for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { - zSum = ySum; - for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) { - polySum += zSum * myPoly->coeff[loop_w][loop_x][loop_y][loop_z]; - zSum *= z; - } - } - ySum *= y; - } - xSum *= x; - } - wSum *= w; - } - - return(polySum); -} - -float p_psChebPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) -{ - psS32 loop_w = 0; - psS32 loop_x = 0; - psS32 loop_y = 0; - psS32 loop_z = 0; - psS32 i = 0; - float polySum = 0.0; - psPolynomial1D* *chebPolys = NULL; - psS32 maxChebyPoly = 0; - - // Determine how many Chebyshev polynomials - // are needed, then create them. - maxChebyPoly = myPoly->nW; - if (myPoly->nX > maxChebyPoly) { - maxChebyPoly = myPoly->nX; - } - if (myPoly->nY > maxChebyPoly) { - maxChebyPoly = myPoly->nY; - } - if (myPoly->nZ > maxChebyPoly) { - maxChebyPoly = myPoly->nZ; - } - chebPolys = CreateChebyshevPolys(maxChebyPoly); - - for (loop_w = 0; loop_w < myPoly->nW; loop_w++) { - for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { - for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { - for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) { - polySum += myPoly->coeff[loop_w][loop_x][loop_y][loop_z] * - psPolynomial1DEval(w, chebPolys[loop_w]) * - psPolynomial1DEval(x, chebPolys[loop_x]) * - psPolynomial1DEval(y, chebPolys[loop_y]) * - psPolynomial1DEval(z, chebPolys[loop_z]); - } - } - } - } - } - - for (i=0;itype == PS_POLYNOMIAL_ORD) { - return(p_psOrdPolynomial4DEval(w,x,y,z, myPoly)); + return(ordPolynomial4DEval(w,x,y,z, myPoly)); } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { - return(p_psChebPolynomial4DEval(w,x,y,z, myPoly)); + return(chebPolynomial4DEval(w,x,y,z, myPoly)); } else { - psError(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, + myPoly->type); } return(0.0); @@ -924,7 +1383,4 @@ return(tmp); } - - - @@ -1092,118 +1548,4 @@ } -static void dPolynomial1DFree(psDPolynomial1D* myPoly) -{ - psFree(myPoly->coeff); - psFree(myPoly->coeffErr); - psFree(myPoly->mask); -} - -static void dPolynomial2DFree(psDPolynomial2D* myPoly) -{ - psS32 x = 0; - - for (x = 0; x < myPoly->nX; x++) { - psFree(myPoly->coeff[x]); - psFree(myPoly->coeffErr[x]); - psFree(myPoly->mask[x]); - } - psFree(myPoly->coeff); - psFree(myPoly->coeffErr); - psFree(myPoly->mask); -} - -static void dPolynomial3DFree(psDPolynomial3D* myPoly) -{ - psS32 x = 0; - psS32 y = 0; - - for (x = 0; x < myPoly->nX; x++) { - for (y = 0; y < myPoly->nY; y++) { - psFree(myPoly->coeff[x][y]); - psFree(myPoly->coeffErr[x][y]); - psFree(myPoly->mask[x][y]); - } - psFree(myPoly->coeff[x]); - psFree(myPoly->coeffErr[x]); - psFree(myPoly->mask[x]); - } - - psFree(myPoly->coeff); - psFree(myPoly->coeffErr); - psFree(myPoly->mask); -} - -static void dPolynomial4DFree(psDPolynomial4D* myPoly) -{ - psS32 w = 0; - psS32 x = 0; - psS32 y = 0; - - for (w = 0; w < myPoly->nW; w++) { - for (x = 0; x < myPoly->nX; x++) { - for (y = 0; y < myPoly->nY; y++) { - psFree(myPoly->coeff[w][x][y]); - psFree(myPoly->coeffErr[w][x][y]); - psFree(myPoly->mask[w][x][y]); - } - psFree(myPoly->coeff[w][x]); - psFree(myPoly->coeffErr[w][x]); - psFree(myPoly->mask[w][x]); - } - psFree(myPoly->coeff[w]); - psFree(myPoly->coeffErr[w]); - psFree(myPoly->mask[w]); - } - - psFree(myPoly->coeff); - psFree(myPoly->coeffErr); - psFree(myPoly->mask); -} - -/***************************************************************************** - Polynomial coefficients will be accessed in [w][x][y][z] fashion. - *****************************************************************************/ -double p_psDOrdPolynomial1DEval(double x, const psDPolynomial1D* myPoly) -{ - psS32 loop_x = 0; - double polySum = 0.0; - double xSum = 1.0; - - for (loop_x = 0; loop_x < myPoly->n; loop_x++) { - if (myPoly->mask[loop_x] == 0) { - polySum += xSum * myPoly->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? -double p_psDChebPolynomial1DEval(double x, const psDPolynomial1D* myPoly) -{ - psVector *d; - psS32 n; - psS32 i; - double tmp; - - n = myPoly->n; - d = psVectorAlloc(n, PS_TYPE_F64); - d->data.F64[n-1] = myPoly->coeff[n-1]; - d->data.F64[n-2] = (2.0 * x * d->data.F64[n-1]) + myPoly->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]) + - (myPoly->coeff[i]); - } - - tmp = (x * d->data.F64[1]) - - (d->data.F64[2]) + - (0.5 * myPoly->coeff[0]); - - psFree(d); - return(tmp); -} double psDPolynomial1DEval(double x, const psDPolynomial1D* myPoly) @@ -1212,9 +1554,11 @@ if (myPoly->type == PS_POLYNOMIAL_ORD) { - return(p_psDOrdPolynomial1DEval(x, myPoly)); + return(dOrdPolynomial1DEval(x, myPoly)); } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { - return(p_psDChebPolynomial1DEval(x, myPoly)); + return(dChebPolynomial1DEval(x, myPoly)); } else { - psError(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, + myPoly->type); } return(0.0); @@ -1244,61 +1588,4 @@ - -double p_psDOrdPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) -{ - psS32 loop_x = 0; - psS32 loop_y = 0; - double polySum = 0.0; - double xSum = 1.0; - double ySum = 1.0; - - for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { - ySum = xSum; - for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { - if (myPoly->mask[loop_x][loop_y] == 0) { - polySum += ySum * myPoly->coeff[loop_x][loop_y]; - ySum *= y; - } - } - xSum *= x; - } - - return(polySum); -} - -double p_psDChebPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) -{ - psS32 loop_x = 0; - psS32 loop_y = 0; - psS32 i = 0; - double polySum = 0.0; - psPolynomial1D* *chebPolys = NULL; - psS32 maxChebyPoly = 0; - - // Determine how many Chebyshev polynomials - // are needed, then create them. - maxChebyPoly = myPoly->nX; - if (myPoly->nY > maxChebyPoly) { - maxChebyPoly = myPoly->nY; - } - chebPolys = CreateChebyshevPolys(maxChebyPoly); - - for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { - for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { - if (myPoly->mask[loop_x][loop_y] == 0) { - polySum += myPoly->coeff[loop_x][loop_y] * - psPolynomial1DEval(x, chebPolys[loop_x]) * - psPolynomial1DEval(y, chebPolys[loop_y]); - } - } - } - - for (i=0;itype == PS_POLYNOMIAL_ORD) { - return(p_psDOrdPolynomial2DEval(x, y, myPoly)); + return(dOrdPolynomial2DEval(x, y, myPoly)); } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { - return(p_psDChebPolynomial2DEval(x, y, myPoly)); + return(dChebPolynomial2DEval(x, y, myPoly)); } else { - psError(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, + myPoly->type); } return(0.0); @@ -1353,74 +1642,4 @@ - -double p_psDOrdPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) -{ - psS32 loop_x = 0; - psS32 loop_y = 0; - psS32 loop_z = 0; - double polySum = 0.0; - double xSum = 1.0; - double ySum = 1.0; - double zSum = 1.0; - - for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { - ySum = xSum; - for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { - zSum = ySum; - for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_x][loop_y][loop_z] == 0) { - polySum += zSum * myPoly->coeff[loop_x][loop_y][loop_z]; - zSum *= z; - } - } - ySum *= y; - } - xSum *= x; - } - - return(polySum); -} - -double p_psDChebPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) -{ - psS32 loop_x = 0; - psS32 loop_y = 0; - psS32 loop_z = 0; - psS32 i = 0; - double polySum = 0.0; - psPolynomial1D* *chebPolys = NULL; - psS32 maxChebyPoly = 0; - - // Determine how many Chebyshev polynomials - // are needed, then create them. - maxChebyPoly = myPoly->nX; - if (myPoly->nY > maxChebyPoly) { - maxChebyPoly = myPoly->nY; - } - if (myPoly->nZ > maxChebyPoly) { - maxChebyPoly = myPoly->nZ; - } - chebPolys = CreateChebyshevPolys(maxChebyPoly); - - for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { - for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { - for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_x][loop_y][loop_z] == 0) { - polySum += myPoly->coeff[loop_x][loop_y][loop_z] * - psPolynomial1DEval(x, chebPolys[loop_x]) * - psPolynomial1DEval(y, chebPolys[loop_y]) * - psPolynomial1DEval(z, chebPolys[loop_z]); - } - } - } - } - - for (i=0;itype == PS_POLYNOMIAL_ORD) { - return(p_psDOrdPolynomial3DEval(x, y, z, myPoly)); + return(dOrdPolynomial3DEval(x, y, z, myPoly)); } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { - return(p_psDChebPolynomial3DEval(x, y, z, myPoly)); + return(dChebPolynomial3DEval(x, y, z, myPoly)); } else { - psError(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, + myPoly->type); } return(0.0); @@ -1485,93 +1706,4 @@ } - - - - - - - -double p_psDOrdPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) -{ - psS32 loop_w = 0; - psS32 loop_x = 0; - psS32 loop_y = 0; - psS32 loop_z = 0; - double polySum = 0.0; - double wSum = 1.0; - double xSum = 1.0; - double ySum = 1.0; - double zSum = 1.0; - - for (loop_w = 0; loop_w < myPoly->nW; loop_w++) { - xSum = wSum; - for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { - ySum = xSum; - for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { - zSum = ySum; - for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) { - polySum += zSum * myPoly->coeff[loop_w][loop_x][loop_y][loop_z]; - zSum *= z; - } - } - ySum *= y; - } - xSum *= x; - } - wSum *= w; - } - - return(polySum); -} - -double p_psDChebPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) -{ - psS32 loop_w = 0; - psS32 loop_x = 0; - psS32 loop_y = 0; - psS32 loop_z = 0; - psS32 i = 0; - double polySum = 0.0; - psPolynomial1D* *chebPolys = NULL; - psS32 maxChebyPoly = 0; - - // Determine how many Chebyshev polynomials - // are needed, then create them. - maxChebyPoly = myPoly->nW; - if (myPoly->nX > maxChebyPoly) { - maxChebyPoly = myPoly->nX; - } - if (myPoly->nY > maxChebyPoly) { - maxChebyPoly = myPoly->nY; - } - if (myPoly->nZ > maxChebyPoly) { - maxChebyPoly = myPoly->nZ; - } - chebPolys = CreateChebyshevPolys(maxChebyPoly); - - for (loop_w = 0; loop_w < myPoly->nW; loop_w++) { - for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { - for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { - for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) { - polySum += myPoly->coeff[loop_w][loop_x][loop_y][loop_z] * - psPolynomial1DEval(w, chebPolys[loop_w]) * - psPolynomial1DEval(x, chebPolys[loop_x]) * - psPolynomial1DEval(y, chebPolys[loop_y]) * - psPolynomial1DEval(z, chebPolys[loop_z]); - } - } - } - } - } - - for (i=0;itype == PS_POLYNOMIAL_ORD) { - return(p_psDOrdPolynomial4DEval(w,x,y,z, myPoly)); + return(dOrdPolynomial4DEval(w,x,y,z, myPoly)); } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { - return(p_psDChebPolynomial4DEval(w,x,y,z, myPoly)); + return(dChebPolynomial4DEval(w,x,y,z, myPoly)); } else { - psError(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, + myPoly->type); } return(0.0); @@ -1699,31 +1833,8 @@ (tmp->domains)[numSplines] = max; + p_psMemSetDeallocator(tmp,(psFreeFcn)spline1DFree); return(tmp); } -// XXX: Have Robert put the dealocator in the memory file. -psS32 p_psSpline1DFree(psSpline1D *tmpSpline) -{ - psS32 i; - - if (tmpSpline == NULL) { - return(0); - } - - 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->domains); - psFree(tmpSpline); - - return(0); -} /***************************************************************************** @@ -1765,5 +1876,5 @@ /***************************************************************************** -p_psVectorBinDisectF32(): This is a private function which takes as input a +vectorBinDisectF32(): 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 @@ -1776,7 +1887,7 @@ XXX: name since we don't take psVectors as input. *****************************************************************************/ -psS32 p_psVectorBinDisectF32(float *bins, - psS32 numBins, - float x) +static psS32 vectorBinDisectF32(float *bins, + psS32 numBins, + float x) { psS32 min; @@ -1784,10 +1895,10 @@ psS32 mid; - psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4, - "---- Calling p_psVectorBinDisectF32(%f)\n", x); + psTrace(".psLib.dataManip.psFunctions.vectorBinDisectF32", 4, + "---- Calling vectorBinDisectF32(%f)\n", x); if (x < bins[0]) { psLogMsg(__func__, PS_LOG_WARN, - "p_psVectorBinDisectF32(): ordinate %f is outside vector range (%f - %f).", + "vectorBinDisectF32(): ordinate %f is outside vector range (%f - %f).", x, bins[0], bins[numBins-1]); return(-2); @@ -1796,5 +1907,5 @@ if (x > bins[numBins-1]) { psLogMsg(__func__, PS_LOG_WARN, - "p_psVectorBinDisectF32(): ordinate %f is outside vector range (%f - %f).", + "vectorBinDisectF32(): ordinate %f is outside vector range (%f - %f).", x, bins[0], bins[numBins-1]); return(-1); @@ -1806,11 +1917,11 @@ while (min != max) { - psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4, + psTrace(".psLib.dataManip.psFunctions.vectorBinDisectF32", 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.p_psVectorBinDisectF32", 4, - "---- Exiting p_psVectorBinDisectF32(): bin %d\n", mid); + psTrace(".psLib.dataManip.psFunctions.vectorBinDisectF32", 4, + "---- Exiting vectorBinDisectF32(): bin %d\n", mid); return(mid); } else if (x < bins[mid]) { @@ -1822,15 +1933,15 @@ } - psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4, - "---- Exiting p_psVectorBinDisectF32(): bin %d\n", min); + psTrace(".psLib.dataManip.psFunctions.vectorBinDisectF32", 4, + "---- Exiting vectorBinDisectF32(): bin %d\n", min); return(min); } /***************************************************************************** -p_psVectorBinDisectS32(): integer version of above. +vectorBinDisectS32(): integer version of above. *****************************************************************************/ -psS32 p_psVectorBinDisectS32(psS32 *bins, - psS32 numBins, - psS32 x) +static psS32 vectorBinDisectS32(psS32 *bins, + psS32 numBins, + psS32 x) { psS32 min; @@ -1838,11 +1949,11 @@ psS32 mid; - psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4, - "---- Calling p_psVectorBinDisectS32(%f)\n", x); + psTrace(".psLib.dataManip.psFunctions.vectorBinDisectS32", 4, + "---- Calling vectorBinDisectS32(%f)\n", x); if ((x < bins[0]) || (x > bins[numBins-1])) { psLogMsg(__func__, PS_LOG_WARN, - "p_psVectorBinDisectS32(): ordinate %f is outside vector range (%f - %f).", + "vectorBinDisectS32(): ordinate %f is outside vector range (%f - %f).", x, bins[0], bins[numBins-1]); return(-1); @@ -1854,11 +1965,11 @@ while (min != max) { - psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4, + psTrace(".psLib.dataManip.psFunctions.vectorBinDisectS32", 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.p_psVectorBinDisectS32", 4, - "---- Exiting p_psVectorBinDisectS32(): bin %d\n", min); + psTrace(".psLib.dataManip.psFunctions.vectorBinDisectS32", 4, + "---- Exiting vectorBinDisectS32(): bin %d\n", min); return(min); } else if (x < bins[mid]) { @@ -1870,6 +1981,6 @@ } - psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4, - "---- Exiting p_psVectorBinDisectS32(): bin %d\n", min); + psTrace(".psLib.dataManip.psFunctions.vectorBinDisectS32", 4, + "---- Exiting vectorBinDisectS32(): bin %d\n", min); return(min); } @@ -1884,120 +1995,16 @@ if (x->type.type == PS_TYPE_S32) { - return(p_psVectorBinDisectS32(bins->data.S32, bins->n, x->data.S32)); + return(vectorBinDisectS32(bins->data.S32, bins->n, x->data.S32)); } else if (x->type.type == PS_TYPE_F32) { - return(p_psVectorBinDisectF32(bins->data.F32, bins->n, x->data.F32)); + return(vectorBinDisectF32(bins->data.F32, bins->n, x->data.F32)); } else { - psError(__func__, "Unallowable data type."); + char* strType; + PS_TYPE_NAME(strType,x->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, + PS_ERRORTEXT_psFunctions_TYPE_NOT_SUPPORTED, + strType); return(-2); } return(-1); -} - -/***************************************************************************** -p_psInterpolate1D(): 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? - *****************************************************************************/ -float p_ps1DFullInterpolateF32(float *domain, - float *range, - psS32 n, - float x) -{ - PS_INT_CHECK_NON_NEGATIVE(n, NAN); - PS_PTR_CHECK_NULL(domain, NAN); - PS_PTR_CHECK_NULL(range, NAN); - - psS32 i; - psS32 m; - static psVector *p = NULL; - p = psVectorRecycle(p, n, PS_TYPE_F32); - p_psMemSetPersistent(p, true); - p_psMemSetPersistent(p->data.F32, true); - /* - psVector *p = psVectorAlloc(n, PS_TYPE_F32); - float tmp; - */ - - psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 4, - "---- p_ps1DFullInterpolateF32() begin (%d-order at x=%f) (%d data points)----\n", n-1, x, n); - - for (i=0;idata.F32[i] = range[i]; - psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 6, - "p->data.F32[%d] is %f\n", i, p->data.F32[i]); - - } - - // From NR, during each iteration of the m loop, we are computing the - // p_{i ... i+m} terms. - for (m=1;mdata.F32[i] = (((x-domain[i+m]) * p->data.F32[i]) + - ((domain[i]-x) * p->data.F32[i+1])) / - (domain[i] - domain[i+m]); - //printf("((%f-%f * %f) + (%f-%f * %f)) / (%f - %f)\n", x, domain[i+m], p->data.F32[i], domain[i], x, p->data.F32[i+1], domain[i], domain[i+m]); - psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 6, - "p->data.F32[%d] is %f\n", i, p->data.F32[i]); - } - } - psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 4, - "---- p_ps1DFullInterpolateF32() end ----\n"); - - /* - tmp = p->data.F32[0]; - psFree(p); - return(tmp); - */ - return(p->data.F32[0]); -} - - -/***************************************************************************** -p_ps1DInterpolateF32(): this is the base 1-D flat memory routine to perform -LaGrange interpolation. - *****************************************************************************/ -float p_ps1DInterpolateF32(float *domain, - float *range, - psS32 n, - psS32 order, - float x) -{ - psS32 binNum; - psS32 numIntPoints = order+1; - psS32 origin; - - psTrace(".psLib.dataManip.psFunctions.p_ps1DInterpolateF32", 4, - "---- p_ps1DInterpolateF32() begin ----\n"); - - binNum = p_psVectorBinDisectF32(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.p_ps1DInterpolateF32", 4, - "---- p_ps1DInterpolateF32() end ----\n"); - return(p_ps1DFullInterpolateF32(&domain[origin], &range[origin], order+1, x)); } @@ -2035,5 +2042,7 @@ if (order > (domain->n - 1)) { - psError(__func__, "not enough data points for %d-order interpolation.\n", order); + psError(PS_ERR_BAD_PARAMETER_SIZE, true, + PS_ERRORTEXT_psFunctions_NOT_ENOUGH_DATAPOINTS, + order); return(NULL); } @@ -2042,9 +2051,9 @@ psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4, "---- p_psVectorInterpolate() end ----\n"); - return(psScalarAlloc(p_ps1DInterpolateF32(domain->data.F32, - range->data.F32, - domain->n, - order, - x->data.F32), PS_TYPE_F32)); + 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. @@ -2053,9 +2062,9 @@ psScalar *tmpScalar = psScalarAlloc((double) - p_ps1DInterpolateF32(domain32->data.F32, - range32->data.F32, - domain32->n, - order, - (float) x->data.F64), PS_TYPE_F64); + interpolate1DF32(domain32->data.F32, + range32->data.F32, + domain32->n, + order, + (float) x->data.F64), PS_TYPE_F64); psFree(range32); psFree(domain32); @@ -2067,6 +2076,9 @@ } else { - // XXX psError: type not supported - psError(__func__, "type %d not supported\n", x->type.type); + char* strType; + PS_TYPE_NAME(strType,x->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, + PS_ERRORTEXT_psFunctions_TYPE_NOT_SUPPORTED, + strType); } @@ -2084,5 +2096,5 @@ and an independent x value. Each determines which spline that x corresponds to by doing a bracket disection on the domains of the spline data structure -(p_psVectorBinDisectF32()). Then it evaluates the spline at that x location +(vectorBinDisectF32()). Then it evaluates the spline at that x location by a call to the 1D polynomial functions. @@ -2100,5 +2112,5 @@ n = spline->n; - binNum = p_psVectorBinDisectF32(spline->domains, (spline->n)+1, x); + binNum = vectorBinDisectF32(spline->domains, (spline->n)+1, x); if (binNum < 0) { psLogMsg(__func__, PS_LOG_WARN, @@ -2139,5 +2151,9 @@ } } else { - psError(__func__, "Unknown data type.\n"); + char* strType; + PS_TYPE_NAME(strType,x->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, + PS_ERRORTEXT_psFunctions_TYPE_NOT_SUPPORTED, + strType); return(NULL); }