Index: trunk/psLib/src/math/psSpline.c =================================================================== --- trunk/psLib/src/math/psSpline.c (revision 4405) +++ trunk/psLib/src/math/psSpline.c (revision 4422) @@ -7,6 +7,6 @@ * polynomials. It also contains a Gaussian functions. * - * @version $Revision: 1.112 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-06-28 00:53:53 $ + * @version $Revision: 1.113 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-06-29 00:43:46 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -46,12 +46,12 @@ /* TYPE DEFINITIONS */ /*****************************************************************************/ -static void polynomial1DFree(psPolynomial1D* myPoly); -static void polynomial2DFree(psPolynomial2D* myPoly); -static void polynomial3DFree(psPolynomial3D* myPoly); -static void polynomial4DFree(psPolynomial4D* myPoly); -static void dPolynomial1DFree(psDPolynomial1D* myPoly); -static void dPolynomial2DFree(psDPolynomial2D* myPoly); -static void dPolynomial3DFree(psDPolynomial3D* myPoly); -static void dPolynomial4DFree(psDPolynomial4D* myPoly); +static void polynomial1DFree(psPolynomial1D* poly); +static void polynomial2DFree(psPolynomial2D* poly); +static void polynomial3DFree(psPolynomial3D* poly); +static void polynomial4DFree(psPolynomial4D* poly); +static void dPolynomial1DFree(psDPolynomial1D* poly); +static void dPolynomial2DFree(psDPolynomial2D* poly); +static void dPolynomial3DFree(psDPolynomial3D* poly); +static void dPolynomial4DFree(psDPolynomial4D* poly); static void spline1DFree(psSpline1D *tmpSpline); static psS32 vectorBinDisectF32(psF32 *bins,psS32 numBins,psF32 x); @@ -97,138 +97,138 @@ } -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]); +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(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) -{ - for (psS32 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(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 dPolynomial1DFree(psDPolynomial1D* poly) +{ + psFree(poly->coeff); + psFree(poly->coeffErr); + psFree(poly->mask); +} + +static void dPolynomial2DFree(psDPolynomial2D* poly) +{ + for (unsigned int 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 dPolynomial3DFree(psDPolynomial3D* 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 dPolynomial4DFree(psDPolynomial4D* 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(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); + 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); } @@ -280,5 +280,5 @@ Polynomial coefficients will be accessed in [w][x][y][z] fashion. *****************************************************************************/ -static psF32 ordPolynomial1DEval(psF32 x, const psPolynomial1D* myPoly) +static psF32 ordPolynomial1DEval(psF32 x, const psPolynomial1D* poly) { psS32 loop_x = 0; @@ -289,15 +289,15 @@ "---- 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++) { + "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 %f\n", loop_x, myPoly->coeff[loop_x]); - } - - for (loop_x = 0; loop_x < myPoly->n; loop_x++) { - if (myPoly->mask[loop_x] == 0) { + "Polynomial coeff[%d] is %f\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 [%f+= (%f * %f)\n", polySum, xSum, myPoly->coeff[loop_x]); - polySum += xSum * myPoly->coeff[loop_x]; + "polysum+= sum*coeff [%f+= (%f * %f)\n", polySum, xSum, poly->coeff[loop_x]); + polySum += xSum * poly->coeff[loop_x]; } xSum *= x; @@ -310,14 +310,14 @@ // XXX: How does the mask vector effect Crenshaw's formula? // XXX: We assume that x is scaled between -1.0 and 1.0; -static psF32 chebPolynomial1DEval(psF32 x, const psPolynomial1D* myPoly) +static psF32 chebPolynomial1DEval(psF32 x, const psPolynomial1D* poly) { PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); // XXX: Create a macro for this in psConstants.h - if (myPoly->n < 1) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Error: Chebyshev polynomial as order %d.", myPoly->n); + 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 = myPoly->n; + psS32 n = poly->n; psS32 i; psF32 tmp = 0.0; @@ -325,6 +325,6 @@ // Special case where the Chebyshev poly is constant. if (n == 1) { - if (myPoly->mask[0] == 0) { - tmp += myPoly->coeff[0]; + if (poly->mask[0] == 0) { + tmp += poly->coeff[0]; } return(tmp); @@ -333,9 +333,9 @@ // Special case where the Chebyshev poly is linear. if (n == 2) { - if (myPoly->mask[0] == 0) { - tmp+= myPoly->coeff[0]; - } - if (myPoly->mask[1] == 0) { - tmp+= myPoly->coeff[1] * x; + if (poly->mask[0] == 0) { + tmp+= poly->coeff[0]; + } + if (poly->mask[1] == 0) { + tmp+= poly->coeff[1] * x; } return(tmp); @@ -344,6 +344,6 @@ // General case where the Chebyshev poly has 2 or more terms. d = psVectorAlloc(n, PS_TYPE_F32); - if(myPoly->mask[n-1] == 0) { - d->data.F32[n-1] = myPoly->coeff[n-1]; + if(poly->mask[n-1] == 0) { + d->data.F32[n-1] = poly->coeff[n-1]; } else { d->data.F32[n-1] = 0.0; @@ -351,6 +351,6 @@ d->data.F32[n-2] = (2.0 * x * d->data.F32[n-1]); - if(myPoly->mask[n-2] == 0) { - d->data.F32[n-2] += myPoly->coeff[n-2]; + if(poly->mask[n-2] == 0) { + d->data.F32[n-2] += poly->coeff[n-2]; } @@ -358,6 +358,6 @@ d->data.F32[i] = (2.0 * x * d->data.F32[i+1]) - (d->data.F32[i+2]); - if(myPoly->mask[i] == 0) { - d->data.F32[i] += myPoly->coeff[i]; + if(poly->mask[i] == 0) { + d->data.F32[i] += poly->coeff[i]; } } @@ -365,6 +365,6 @@ tmp = (x * d->data.F32[1]) - (d->data.F32[2]); - if(myPoly->mask[0] == 0) { - tmp += (0.5 * myPoly->coeff[0]); + if(poly->mask[0] == 0) { + tmp += (0.5 * poly->coeff[0]); } psFree(d); @@ -378,12 +378,12 @@ psPolynomial1D **chebPolys = NULL; - n = myPoly->n; + n = poly->n; chebPolys = createChebyshevPolys(n); tmp = 0.0; - for (i=0;in;i++) { - tmp+= (myPoly->coeff[i] * psPolynomial1DEval(x, chebPolys[i])); - } - tmp-= (myPoly->coeff[0]/2.0); + for (i=0;in;i++) { + tmp+= (poly->coeff[i] * psPolynomial1DEval(x, chebPolys[i])); + } + tmp-= (poly->coeff[0]/2.0); @@ -394,7 +394,7 @@ static psF32 ordPolynomial2DEval(psF32 x, psF32 y, - const psPolynomial2D* myPoly) -{ - PS_ASSERT_POLY_NON_NULL(myPoly, NAN); + const psPolynomial2D* poly) +{ + PS_ASSERT_POLY_NON_NULL(poly, NAN); psS32 loop_x = 0; @@ -404,9 +404,9 @@ psF32 ySum = 1.0; - for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { + for (loop_x = 0; loop_x < poly->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]; + 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; @@ -418,9 +418,9 @@ } -static psF32 chebPolynomial2DEval(psF32 x, psF32 y, const psPolynomial2D* myPoly) +static psF32 chebPolynomial2DEval(psF32 x, psF32 y, const psPolynomial2D* poly) { PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0); - PS_ASSERT_POLY_NON_NULL(myPoly, NAN); + PS_ASSERT_POLY_NON_NULL(poly, NAN); psS32 loop_x = 0; @@ -433,14 +433,14 @@ // Determine how many Chebyshev polynomials // are needed, then create them. - maxChebyPoly = myPoly->nX; - if (myPoly->nY > maxChebyPoly) { - maxChebyPoly = myPoly->nY; + maxChebyPoly = poly->nX; + if (poly->nY > maxChebyPoly) { + maxChebyPoly = poly->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] * + 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); @@ -455,5 +455,5 @@ } -static psF32 ordPolynomial3DEval(psF32 x, psF32 y, psF32 z, const psPolynomial3D* myPoly) +static psF32 ordPolynomial3DEval(psF32 x, psF32 y, psF32 z, const psPolynomial3D* poly) { psS32 loop_x = 0; @@ -465,11 +465,11 @@ psF32 zSum = 1.0; - for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { + for (loop_x = 0; loop_x < poly->nX; loop_x++) { ySum = xSum; - for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { + for (loop_y = 0; loop_y < poly->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]; + 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; @@ -483,5 +483,5 @@ } -static psF32 chebPolynomial3DEval(psF32 x, psF32 y, psF32 z, const psPolynomial3D* myPoly) +static psF32 chebPolynomial3DEval(psF32 x, psF32 y, psF32 z, const psPolynomial3D* poly) { PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); @@ -498,18 +498,18 @@ // 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; + 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 < 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] * + 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) * @@ -527,33 +527,33 @@ } -static psF32 ordPolynomial4DEval(psF32 w, psF32 x, psF32 y, psF32 z, const psPolynomial4D* myPoly) -{ - psS32 loop_w = 0; +static psF32 ordPolynomial4DEval(psF32 x, psF32 y, psF32 z, psF32 t, const psPolynomial4D* poly) +{ psS32 loop_x = 0; psS32 loop_y = 0; psS32 loop_z = 0; + psS32 loop_t = 0; psF32 polySum = 0.0; - psF32 wSum = 1.0; psF32 xSum = 1.0; psF32 ySum = 1.0; psF32 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]; + psF32 tSum = 1.0; + + for (loop_x = 0; loop_x < poly->nX; 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]; } - zSum *= z; + tSum *= t; } - ySum *= y; + zSum *= z; } - xSum *= x; - } - wSum *= w; + ySum *= y; + } + xSum *= x; } @@ -561,14 +561,14 @@ } -static psF32 chebPolynomial4DEval(psF32 w, psF32 x, psF32 y, psF32 z, const psPolynomial4D* myPoly) -{ - PS_ASSERT_FLOAT_WITHIN_RANGE(w, -1.0, 1.0, 0.0); +static psF32 chebPolynomial4DEval(psF32 x, psF32 y, psF32 z, psF32 t, const psPolynomial4D* poly) +{ PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0); PS_ASSERT_FLOAT_WITHIN_RANGE(z, -1.0, 1.0, 0.0); - psS32 loop_w = 0; + PS_ASSERT_FLOAT_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; psF32 polySum = 0.0; @@ -578,26 +578,26 @@ // 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; + 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_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(chebPolys[loop_w], w) * + 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_z], z) * + psPolynomial1DEval(chebPolys[loop_t], t); } } @@ -616,5 +616,5 @@ Polynomial coefficients will be accessed in [w][x][y][z] fashion. *****************************************************************************/ -static psF64 dOrdPolynomial1DEval(psF64 x, const psDPolynomial1D* myPoly) +static psF64 dOrdPolynomial1DEval(psF64 x, const psDPolynomial1D* poly) { psS32 loop_x = 0; @@ -622,7 +622,7 @@ psF64 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]; + for (loop_x = 0; loop_x < poly->n; loop_x++) { + if (poly->mask[loop_x] == 0) { + polySum += xSum * poly->coeff[loop_x]; } xSum *= x; @@ -634,5 +634,5 @@ // XXX: You can do this without having to psAlloc() vector d. // XXX: How does the mask vector effect Crenshaw's formula? -static psF64 dChebPolynomial1DEval(psF64 x, const psDPolynomial1D* myPoly) +static psF64 dChebPolynomial1DEval(psF64 x, const psDPolynomial1D* poly) { PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); @@ -642,20 +642,20 @@ psF64 tmp; - n = myPoly->n; + n = poly->n; d = psVectorAlloc(n, PS_TYPE_F64); - if(myPoly->mask[n-1] == 0) { - d->data.F64[n-1] = myPoly->coeff[n-1]; + 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(myPoly->mask[n-2] == 0) { - d->data.F64[n-2] += myPoly->coeff[n-2]; + 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(myPoly->mask[i] == 0) { - d->data.F64[i] += myPoly->coeff[i]; + if(poly->mask[i] == 0) { + d->data.F64[i] += poly->coeff[i]; } } @@ -663,6 +663,6 @@ tmp = (x * d->data.F64[1]) - (d->data.F64[2]); - if(myPoly->mask[0] == 0) { - tmp += (0.5 * myPoly->coeff[0]); + if(poly->mask[0] == 0) { + tmp += (0.5 * poly->coeff[0]); } @@ -673,5 +673,5 @@ static psF64 dOrdPolynomial2DEval(psF64 x, psF64 y, - const psDPolynomial2D* myPoly) + const psDPolynomial2D* poly) { psS32 loop_x = 0; @@ -681,9 +681,9 @@ psF64 ySum = 1.0; - for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { + for (loop_x = 0; loop_x < poly->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]; + 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; @@ -695,5 +695,5 @@ } -static psF64 dChebPolynomial2DEval(psF64 x, psF64 y, const psDPolynomial2D* myPoly) +static psF64 dChebPolynomial2DEval(psF64 x, psF64 y, const psDPolynomial2D* poly) { PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); @@ -708,14 +708,14 @@ // Determine how many Chebyshev polynomials // are needed, then create them. - maxChebyPoly = myPoly->nX; - if (myPoly->nY > maxChebyPoly) { - maxChebyPoly = myPoly->nY; + maxChebyPoly = poly->nX; + if (poly->nY > maxChebyPoly) { + maxChebyPoly = poly->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] * + 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); @@ -731,5 +731,5 @@ } -static psF64 dOrdPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* myPoly) +static psF64 dOrdPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* poly) { psS32 loop_x = 0; @@ -741,11 +741,11 @@ psF64 zSum = 1.0; - for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { + for (loop_x = 0; loop_x < poly->nX; loop_x++) { ySum = xSum; - for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { + for (loop_y = 0; loop_y < poly->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]; + 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; @@ -759,5 +759,5 @@ } -static psF64 dChebPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* myPoly) +static psF64 dChebPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* poly) { PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); @@ -774,18 +774,18 @@ // 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; + 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 < 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] * + 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) * @@ -803,33 +803,33 @@ } -static psF64 dOrdPolynomial4DEval(psF64 w, psF64 x, psF64 y, psF64 z, const psDPolynomial4D* myPoly) -{ - psS32 loop_w = 0; +static psF64 dOrdPolynomial4DEval(psF64 x, psF64 y, psF64 z, psF64 t, const psDPolynomial4D* poly) +{ psS32 loop_x = 0; psS32 loop_y = 0; psS32 loop_z = 0; + psS32 loop_t = 0; psF64 polySum = 0.0; - psF64 wSum = 1.0; psF64 xSum = 1.0; psF64 ySum = 1.0; psF64 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]; + psF64 tSum = 1.0; + + for (loop_x = 0; loop_x < poly->nX; 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]; } - zSum *= z; + tSum *= t; } - ySum *= y; + zSum *= z; } - xSum *= x; - } - wSum *= w; + ySum *= y; + } + xSum *= x; } @@ -837,14 +837,14 @@ } -static psF64 dChebPolynomial4DEval(psF64 w, psF64 x, psF64 y, psF64 z, const psDPolynomial4D* myPoly) -{ - PS_ASSERT_FLOAT_WITHIN_RANGE(w, -1.0, 1.0, 0.0); +static psF64 dChebPolynomial4DEval(psF64 x, psF64 y, psF64 z, psF64 t, const psDPolynomial4D* poly) +{ PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0); PS_ASSERT_FLOAT_WITHIN_RANGE(z, -1.0, 1.0, 0.0); - psS32 loop_w = 0; + PS_ASSERT_FLOAT_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; @@ -854,26 +854,26 @@ // 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; + 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_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(chebPolys[loop_w], w) * + 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_z], z) * + psPolynomial1DEval(chebPolys[loop_t], t); } } @@ -1069,10 +1069,10 @@ This routine must allocate memory for the polynomial structures. *****************************************************************************/ -psPolynomial1D* psPolynomial1DAlloc(psS32 n, +psPolynomial1D* psPolynomial1DAlloc(int n, psPolynomialType type) { PS_ASSERT_INT_POSITIVE(n, NULL); - psS32 i = 0; + int i = 0; psPolynomial1D* newPoly = NULL; @@ -1094,5 +1094,5 @@ } -psPolynomial2D* psPolynomial2DAlloc(psS32 nX, psS32 nY, +psPolynomial2D* psPolynomial2DAlloc(int nX, int nY, psPolynomialType type) { @@ -1100,6 +1100,6 @@ PS_ASSERT_INT_POSITIVE(nY, NULL); - psS32 x = 0; - psS32 y = 0; + int x = 0; + int y = 0; psPolynomial2D* newPoly = NULL; @@ -1130,5 +1130,5 @@ } -psPolynomial3D* psPolynomial3DAlloc(psS32 nX, psS32 nY, psS32 nZ, +psPolynomial3D* psPolynomial3DAlloc(int nX, int nY, int nZ, psPolynomialType type) { @@ -1176,16 +1176,16 @@ } -psPolynomial4D* psPolynomial4DAlloc(psS32 nW, psS32 nX, psS32 nY, psS32 nZ, +psPolynomial4D* psPolynomial4DAlloc(int nX, int nY, int nZ, int nT, psPolynomialType type) { - PS_ASSERT_INT_POSITIVE(nW, NULL); PS_ASSERT_INT_POSITIVE(nX, NULL); PS_ASSERT_INT_POSITIVE(nY, NULL); PS_ASSERT_INT_POSITIVE(nZ, NULL); - - psS32 w = 0; + PS_ASSERT_INT_POSITIVE(nT, NULL); + psS32 x = 0; psS32 y = 0; psS32 z = 0; + psS32 t = 0; psPolynomial4D* newPoly = NULL; @@ -1194,34 +1194,34 @@ newPoly->type = type; - newPoly->nW = nW; newPoly->nX = nX; newPoly->nY = nY; newPoly->nZ = nZ; - - newPoly->coeff = (psF32 ****)psAlloc(nW * sizeof(psF32 ***)); - newPoly->coeffErr = (psF32 ****)psAlloc(nW * sizeof(psF32 ***)); - newPoly->mask = (psU8 ****)psAlloc(nW * sizeof(psU8 ***)); - for (w = 0; w < nW; w++) { - newPoly->coeff[w] = (psF32 ***)psAlloc(nX * sizeof(psF32 **)); - newPoly->coeffErr[w] = (psF32 ***)psAlloc(nX * sizeof(psF32 **)); - newPoly->mask[w] = (psU8 ***)psAlloc(nX * sizeof(psU8 **)); - for (x = 0; x < nX; x++) { - newPoly->coeff[w][x] = (psF32 **)psAlloc(nY * sizeof(psF32 *)); - newPoly->coeffErr[w][x] = (psF32 **)psAlloc(nY * sizeof(psF32 *)); - newPoly->mask[w][x] = (psU8 **)psAlloc(nY * sizeof(psU8 *)); - for (y = 0; y < nY; y++) { - newPoly->coeff[w][x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32)); - newPoly->coeffErr[w][x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32)); - newPoly->mask[w][x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8)); + newPoly->nT = nT; + + newPoly->coeff = (psF32 ****)psAlloc(nX * sizeof(psF32 ***)); + newPoly->coeffErr = (psF32 ****)psAlloc(nX * sizeof(psF32 ***)); + newPoly->mask = (psU8 ****)psAlloc(nX * sizeof(psU8 ***)); + for (x = 0; x < nX; x++) { + newPoly->coeff[x] = (psF32 ***)psAlloc(nY * sizeof(psF32 **)); + newPoly->coeffErr[x] = (psF32 ***)psAlloc(nY * sizeof(psF32 **)); + newPoly->mask[x] = (psU8 ***)psAlloc(nY * sizeof(psU8 **)); + for (y = 0; y < nY; y++) { + newPoly->coeff[x][y] = (psF32 **)psAlloc(nZ * sizeof(psF32 *)); + newPoly->coeffErr[x][y] = (psF32 **)psAlloc(nZ * sizeof(psF32 *)); + newPoly->mask[x][y] = (psU8 **)psAlloc(nZ * sizeof(psU8 *)); + for (z = 0; z < nZ; z++) { + newPoly->coeff[x][y][z] = (psF32 *)psAlloc(nT * sizeof(psF32)); + newPoly->coeffErr[x][y][z] = (psF32 *)psAlloc(nT * sizeof(psF32)); + newPoly->mask[x][y][z] = (psU8 *)psAlloc(nT * sizeof(psU8)); } } } - for (w = 0; w < nW; w++) { - for (x = 0; x < nX; x++) { - for (y = 0; y < nY; y++) { - for (z = 0; z < nZ; z++) { - newPoly->coeff[w][x][y][z] = 0.0; - newPoly->coeffErr[w][x][y][z] = 0.0; - newPoly->mask[w][x][y][z] = 0; + 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; } } @@ -1431,10 +1431,10 @@ -psDPolynomial1D* psDPolynomial1DAlloc(psS32 n, +psDPolynomial1D* psDPolynomial1DAlloc(int n, psPolynomialType type) { PS_ASSERT_INT_POSITIVE(n, NULL); - psS32 i = 0; + unsigned int i = 0; psDPolynomial1D* newPoly = NULL; @@ -1456,5 +1456,5 @@ } -psDPolynomial2D* psDPolynomial2DAlloc(psS32 nX, psS32 nY, +psDPolynomial2D* psDPolynomial2DAlloc(int nX, int nY, psPolynomialType type) { @@ -1462,6 +1462,6 @@ PS_ASSERT_INT_POSITIVE(nY, NULL); - psS32 x = 0; - psS32 y = 0; + unsigned int x = 0; + unsigned int y = 0; psDPolynomial2D* newPoly = NULL; @@ -1492,5 +1492,5 @@ } -psDPolynomial3D* psDPolynomial3DAlloc(psS32 nX, psS32 nY, psS32 nZ, +psDPolynomial3D* psDPolynomial3DAlloc(int nX, int nY, int nZ, psPolynomialType type) { @@ -1499,7 +1499,7 @@ PS_ASSERT_INT_POSITIVE(nZ, NULL); - psS32 x = 0; - psS32 y = 0; - psS32 z = 0; + unsigned int x = 0; + unsigned int y = 0; + unsigned int z = 0; psDPolynomial3D* newPoly = NULL; @@ -1538,16 +1538,16 @@ } -psDPolynomial4D* psDPolynomial4DAlloc(psS32 nW, psS32 nX, psS32 nY, psS32 nZ, +psDPolynomial4D* psDPolynomial4DAlloc(int nX, int nY, int nZ, int nT, psPolynomialType type) { - PS_ASSERT_INT_POSITIVE(nW, NULL); PS_ASSERT_INT_POSITIVE(nX, NULL); PS_ASSERT_INT_POSITIVE(nY, NULL); PS_ASSERT_INT_POSITIVE(nZ, NULL); - - psS32 w = 0; - psS32 x = 0; - psS32 y = 0; - psS32 z = 0; + PS_ASSERT_INT_POSITIVE(nT, NULL); + + unsigned int x = 0; + unsigned int y = 0; + unsigned int z = 0; + unsigned int t = 0; psDPolynomial4D* newPoly = NULL; @@ -1556,34 +1556,34 @@ newPoly->type = type; - newPoly->nW = nW; newPoly->nX = nX; newPoly->nY = nY; newPoly->nZ = nZ; - - newPoly->coeff = (psF64 ****)psAlloc(nW * sizeof(psF64 ***)); - newPoly->coeffErr = (psF64 ****)psAlloc(nW * sizeof(psF64 ***)); - newPoly->mask = (psU8 ****)psAlloc(nW * sizeof(psU8 ***)); - for (w = 0; w < nW; w++) { - newPoly->coeff[w] = (psF64 ***)psAlloc(nX * sizeof(psF64 **)); - newPoly->coeffErr[w] = (psF64 ***)psAlloc(nX * sizeof(psF64 **)); - newPoly->mask[w] = (psU8 ***)psAlloc(nX * sizeof(psU8 **)); - for (x = 0; x < nX; x++) { - newPoly->coeff[w][x] = (psF64 **)psAlloc(nY * sizeof(psF64 *)); - newPoly->coeffErr[w][x] = (psF64 **)psAlloc(nY * sizeof(psF64 *)); - newPoly->mask[w][x] = (psU8 **)psAlloc(nY * sizeof(psU8 *)); - for (y = 0; y < nY; y++) { - newPoly->coeff[w][x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64)); - newPoly->coeffErr[w][x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64)); - newPoly->mask[w][x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8)); + newPoly->nT = nT; + + newPoly->coeff = (psF64 ****)psAlloc(nX * sizeof(psF64 ***)); + newPoly->coeffErr = (psF64 ****)psAlloc(nX * sizeof(psF64 ***)); + newPoly->mask = (psU8 ****)psAlloc(nX * sizeof(psU8 ***)); + for (x = 0; x < nX; x++) { + newPoly->coeff[x] = (psF64 ***)psAlloc(nY * sizeof(psF64 **)); + newPoly->coeffErr[x] = (psF64 ***)psAlloc(nY * sizeof(psF64 **)); + newPoly->mask[x] = (psU8 ***)psAlloc(nY * sizeof(psU8 **)); + for (y = 0; y < nY; y++) { + newPoly->coeff[x][y] = (psF64 **)psAlloc(nZ * sizeof(psF64 *)); + newPoly->coeffErr[x][y] = (psF64 **)psAlloc(nZ * sizeof(psF64 *)); + newPoly->mask[x][y] = (psU8 **)psAlloc(nZ * sizeof(psU8 *)); + for (z = 0; z < nZ; z++) { + newPoly->coeff[x][y][z] = (psF64 *)psAlloc(nT * sizeof(psF64)); + newPoly->coeffErr[x][y][z] = (psF64 *)psAlloc(nT * sizeof(psF64)); + newPoly->mask[x][y][z] = (psU8 *)psAlloc(nT * sizeof(psU8)); } } } - for (w = 0; w < nW; w++) { - for (x = 0; x < nX; x++) { - for (y = 0; y < nY; y++) { - for (z = 0; z < nZ; z++) { - newPoly->coeff[w][x][y][z] = 0.0; - newPoly->coeffErr[w][x][y][z] = 0.0; - newPoly->mask[w][x][y][z] = 0; + 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; } } @@ -1595,25 +1595,25 @@ -psF64 psDPolynomial1DEval(const psDPolynomial1D* myPoly, psF64 x) -{ - PS_ASSERT_POLY_NON_NULL(myPoly, NAN); - - if (myPoly->type == PS_POLYNOMIAL_ORD) { - return(dOrdPolynomial1DEval(x, myPoly)); - } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { - return(dChebPolynomial1DEval(x, myPoly)); +psF64 psDPolynomial1DEval(const psDPolynomial1D* poly, psF64 x) +{ + PS_ASSERT_POLY_NON_NULL(poly, NAN); + + if (poly->type == PS_POLYNOMIAL_ORD) { + return(dOrdPolynomial1DEval(x, poly)); + } else if (poly->type == PS_POLYNOMIAL_CHEB) { + return(dChebPolynomial1DEval(x, poly)); } else { psError(PS_ERR_BAD_PARAMETER_TYPE, true, PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, - myPoly->type); + poly->type); } return(NAN); } -psVector *psDPolynomial1DEvalVector(const psDPolynomial1D *myPoly, +psVector *psDPolynomial1DEvalVector(const psDPolynomial1D *poly, const psVector *x) { - PS_ASSERT_POLY_NON_NULL(myPoly, NULL); + PS_ASSERT_POLY_NON_NULL(poly, NULL); PS_ASSERT_VECTOR_NON_NULL(x, NULL); PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); @@ -1623,5 +1623,5 @@ tmp = psVectorAlloc(x->n, PS_TYPE_F64); for (psS32 i=0;in;i++) { - tmp->data.F64[i] = psDPolynomial1DEval(myPoly, + tmp->data.F64[i] = psDPolynomial1DEval(poly, x->data.F64[i]); } @@ -1631,26 +1631,26 @@ -psF64 psDPolynomial2DEval(const psDPolynomial2D* myPoly, +psF64 psDPolynomial2DEval(const psDPolynomial2D* poly, psF64 x, psF64 y) { - PS_ASSERT_POLY_NON_NULL(myPoly, NAN); - if (myPoly->type == PS_POLYNOMIAL_ORD) { - return(dOrdPolynomial2DEval(x, y, myPoly)); - } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { - return(dChebPolynomial2DEval(x, y, myPoly)); + PS_ASSERT_POLY_NON_NULL(poly, NAN); + if (poly->type == PS_POLYNOMIAL_ORD) { + return(dOrdPolynomial2DEval(x, y, poly)); + } else if (poly->type == PS_POLYNOMIAL_CHEB) { + return(dChebPolynomial2DEval(x, y, poly)); } else { psError(PS_ERR_BAD_PARAMETER_TYPE, true, PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, - myPoly->type); + poly->type); } return(NAN); } -psVector *psDPolynomial2DEvalVector(const psDPolynomial2D *myPoly, +psVector *psDPolynomial2DEvalVector(const psDPolynomial2D *poly, const psVector *x, const psVector *y) { - PS_ASSERT_POLY_NON_NULL(myPoly, NULL); + PS_ASSERT_POLY_NON_NULL(poly, NULL); PS_ASSERT_VECTOR_NON_NULL(x, NULL); PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); @@ -1671,5 +1671,5 @@ // Evaluate the polynomial for (psS32 i = 0; i < vecLen; i++) { - tmp->data.F64[i] = psDPolynomial2DEval(myPoly,x->data.F64[i],y->data.F64[i]); + tmp->data.F64[i] = psDPolynomial2DEval(poly,x->data.F64[i],y->data.F64[i]); } @@ -1679,24 +1679,24 @@ -psF64 psDPolynomial3DEval(const psDPolynomial3D* myPoly, +psF64 psDPolynomial3DEval(const psDPolynomial3D* poly, psF64 x, psF64 y, psF64 z) { - PS_ASSERT_POLY_NON_NULL(myPoly, NAN); - - if (myPoly->type == PS_POLYNOMIAL_ORD) { - return(dOrdPolynomial3DEval(x, y, z, myPoly)); - } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { - return(dChebPolynomial3DEval(x, y, z, myPoly)); + PS_ASSERT_POLY_NON_NULL(poly, NAN); + + if (poly->type == PS_POLYNOMIAL_ORD) { + return(dOrdPolynomial3DEval(x, y, z, poly)); + } else if (poly->type == PS_POLYNOMIAL_CHEB) { + return(dChebPolynomial3DEval(x, y, z, poly)); } else { psError(PS_ERR_BAD_PARAMETER_TYPE, true, PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, - myPoly->type); + poly->type); } return(NAN); } -psVector *psDPolynomial3DEvalVector(const psDPolynomial3D *myPoly, +psVector *psDPolynomial3DEvalVector(const psDPolynomial3D *poly, const psVector *x, const psVector *y, @@ -1704,5 +1704,5 @@ { - PS_ASSERT_POLY_NON_NULL(myPoly, NULL); + PS_ASSERT_POLY_NON_NULL(poly, NULL); PS_ASSERT_VECTOR_NON_NULL(x, NULL); PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); @@ -1728,5 +1728,5 @@ // Evaluate polynomial for (psS32 i = 0; i < vecLen; i++) { - tmp->data.F64[i] = psDPolynomial3DEval(myPoly, + tmp->data.F64[i] = psDPolynomial3DEval(poly, x->data.F64[i], y->data.F64[i], @@ -1738,33 +1738,31 @@ } -psF64 psDPolynomial4DEval(const psDPolynomial4D* myPoly, - psF64 w, +psF64 psDPolynomial4DEval(const psDPolynomial4D* poly, psF64 x, psF64 y, - psF64 z) -{ - PS_ASSERT_POLY_NON_NULL(myPoly, NAN); - - if (myPoly->type == PS_POLYNOMIAL_ORD) { - return(dOrdPolynomial4DEval(w,x,y,z, myPoly)); - } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { - return(dChebPolynomial4DEval(w,x,y,z, myPoly)); + psF64 z, + psF64 t) +{ + PS_ASSERT_POLY_NON_NULL(poly, NAN); + + if (poly->type == PS_POLYNOMIAL_ORD) { + return(dOrdPolynomial4DEval(x,y,z,t, poly)); + } else if (poly->type == PS_POLYNOMIAL_CHEB) { + return(dChebPolynomial4DEval(x,y,z,t, poly)); } else { psError(PS_ERR_BAD_PARAMETER_TYPE, true, PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, - myPoly->type); + poly->type); } return(NAN); } -psVector *psDPolynomial4DEvalVector(const psDPolynomial4D *myPoly, - const psVector *w, +psVector *psDPolynomial4DEvalVector(const psDPolynomial4D *poly, const psVector *x, const psVector *y, - const psVector *z) -{ - PS_ASSERT_POLY_NON_NULL(myPoly, NULL); - PS_ASSERT_VECTOR_NON_NULL(w, NULL); - PS_ASSERT_VECTOR_TYPE(w, PS_TYPE_F64, NULL); + 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); @@ -1773,17 +1771,19 @@ 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=w->n; + psS32 vecLen=x->n; // Determine the output vector size from min of input vectors + if (z->n < vecLen) { + vecLen = z->n; + } if (y->n < vecLen) { vecLen = y->n; } - if (x->n < vecLen) { - vecLen = x->n; - } - if (z->n < vecLen) { - vecLen = z->n; + if (t->n < vecLen) { + vecLen = t->n; } @@ -1793,9 +1793,9 @@ // Evaluate the polynomial for (psS32 i = 0; i < vecLen; i++) { - tmp->data.F64[i] = psDPolynomial4DEval(myPoly, - w->data.F64[i], + tmp->data.F64[i] = psDPolynomial4DEval(poly, x->data.F64[i], y->data.F64[i], - z->data.F64[i]); + z->data.F64[i], + t->data.F64[i]); }