Index: trunk/psLib/src/math/psSpline.c =================================================================== --- trunk/psLib/src/math/psSpline.c (revision 4991) +++ trunk/psLib/src/math/psSpline.c (revision 5066) @@ -7,6 +7,6 @@ * splines. * -* @version $Revision: 1.123 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-09-11 22:18:40 $ +* @version $Revision: 1.124 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-09-19 19:53:13 $ * * @@ -43,6 +43,6 @@ /*****************************************************************************/ static void spline1DFree(psSpline1D *tmpSpline); -static psS32 vectorBinDisectF32(psF32 *bins,psS32 numBins,psF32 x); -static psS32 vectorBinDisectS32(psS32 *bins,psS32 numBins,psS32 x); +static unsigned int vectorBinDisectF32(psF32 *bins,psS32 numBins,psF32 x); +static unsigned int vectorBinDisectS32(psS32 *bins,psS32 numBins,psS32 x); /*****************************************************************************/ @@ -69,5 +69,5 @@ static void spline1DFree(psSpline1D *tmpSpline) { - psS32 i; + unsigned int i; if (tmpSpline == NULL) { @@ -101,10 +101,10 @@ static psF32 fullInterpolate1D##TYPE(ps##TYPE *domain, \ ps##TYPE *range, \ - psS32 n, \ + unsigned int n, \ ps##TYPE x) \ { \ \ - psS32 i; \ - psS32 m; \ + unsigned int i; \ + unsigned int m; \ static psVector *p = NULL; \ p = psVectorRecycle(p, n, PS_TYPE_##TYPE); \ @@ -169,6 +169,6 @@ static psF32 interpolate1DF32(psF32 *domain, psF32 *range, - psS32 n, - psS32 order, + unsigned int n, + unsigned int order, psF32 x) { @@ -178,5 +178,5 @@ psS32 binNum; - psS32 numIntPoints = order+1; + unsigned int numIntPoints = order+1; psS32 origin; @@ -229,9 +229,9 @@ "---- calculateSecondDerivs() begin ----\n"); - psS32 i; - psS32 k; + unsigned int i; + unsigned int k; psF32 sig; psF32 p; - psS32 n = y->n; + unsigned int n = y->n; psF32 *u = (psF32 *) psAlloc(n * sizeof(psF32)); psF32 *derivs2 = (psF32 *) psAlloc(n * sizeof(psF32)); @@ -253,9 +253,9 @@ psTrace(".psLib.dataManip.calculateSecondDerivs", 6, - "X[%d] is %f\n", i, X[i]); + "X[%u] is %f\n", i, X[i]); psTrace(".psLib.dataManip.calculateSecondDerivs", 6, - "Y[%d] is %f\n", i, Y[i]); + "Y[%u] is %f\n", i, Y[i]); psTrace(".psLib.dataManip.calculateSecondDerivs", 6, - "u[%d] is %f\n", i, u[i]); + "u[%u] is %f\n", i, u[i]); } @@ -268,5 +268,5 @@ psTrace(".psLib.dataManip.calculateSecondDerivs", 6, - "derivs2[%d] is %f\n", k, derivs2[k]); + "derivs2[%u] is %f\n", k, derivs2[k]); } @@ -314,5 +314,5 @@ XXX: Assumes mySpline->knots is psF32. Must add psU32 and psF64. *****************************************************************************/ -psSpline1D *psVectorFitSpline1D(psSpline1D *mySpline, ///< The spline which will be generated. +psSpline1D *psVectorFitSpline1D(psSpline1D *spline, ///< The spline which will be generated. const psVector* x, ///< Ordinates (or NULL to just use the indices) const psVector* y, ///< Coordinates @@ -321,13 +321,13 @@ PS_ASSERT_VECTOR_NON_NULL(y, NULL); PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL); - if (mySpline != NULL) { - PS_ASSERT_VECTOR_TYPE(mySpline->knots, PS_TYPE_F32, NULL); + if (spline != NULL) { + PS_ASSERT_VECTOR_TYPE(spline->knots, PS_TYPE_F32, NULL); } psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, "---- psVectorFitSpline1D() begin ----\n"); - psS32 numSplines = (y->n)-1; + unsigned int numSplines = (y->n)-1; psF32 tmp; psF32 H; - psS32 i; + unsigned int i; psF32 slope; psVector *x32 = NULL; @@ -364,17 +364,17 @@ This can not be implemented until SDR states what order spline should be created. - Should we error if mySpline is not NULL? + Should we error if spline is not NULL? Should we error if mySPline is not NULL? */ - if (mySpline == NULL) { - mySpline = psSpline1DAllocGeneric(x32, 3); - } - PS_ASSERT_PTR_NON_NULL(mySpline, NULL); - PS_ASSERT_INT_NONNEGATIVE(mySpline->n, NULL); - - if (y32->n != (1 + mySpline->n)) { + if (spline == NULL) { + spline = psSpline1DAllocGeneric(x32, 3); + } + PS_ASSERT_PTR_NON_NULL(spline, NULL); + PS_ASSERT_INT_NONNEGATIVE(spline->n, NULL); + + if (y32->n != (1 + spline->n)) { psError(PS_ERR_BAD_PARAMETER_SIZE, true, - "data size / spline size mismatch (%d %d)\n", - y32->n, mySpline->n); + "data size / spline size mismatch (%u %u)\n", + y32->n, spline->n); return(NULL); } @@ -382,5 +382,205 @@ // If these are linear splines, which means their polynomials will have // two coefficients, then we do the simple calculation. - if (2 == (mySpline->spline[0])->n) { + if (2 == (spline->spline[0])->n) { + for (i=0;in;i++) { + slope = (y32->data.F32[i+1] - y32->data.F32[i]) / + (spline->knots->data.F32[i+1] - spline->knots->data.F32[i]); + (spline->spline[i])->coeff[0] = y32->data.F32[i] - + (slope * spline->knots->data.F32[i]); + + (spline->spline[i])->coeff[1] = slope; + psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4, + "---- spline %u coeffs are (%f, %f)\n", i, + (spline->spline[i])->coeff[0], + (spline->spline[i])->coeff[1]); + } + psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4, + "---- Exiting psVectorFitSpline1D()()\n"); + return((psSpline1D *) spline); + } + + // Check if these are cubic splines (n==4). If not, psError. + if (4 != (spline->spline[0])->n) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, + "Don't know how to generate %u-order splines.", + (spline->spline[0])->n-1); + return(NULL); + } + + // If we get here, then we know these are cubic splines. We first + // generate the second derivatives at each data point. + spline->p_psDeriv2 = calculateSecondDerivs(x32, y32); + for (i=0;in;i++) + psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, + "Second deriv[%u] is %f\n", i, spline->p_psDeriv2[i]); + + // We generate the coefficients of the spline polynomials. I can't + // concisely explain how this code works. See above function comments + // and Numerical Recipes in C. + for (i=0;idata.F32[i+1] - x32->data.F32[i]; + psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, + "x data (%f - %f) (%f)\n", + x32->data.F32[i], + x32->data.F32[i+1], H); + // + // ******** Calculate 0-order term ******** + // + // From (1) + (spline->spline[i])->coeff[0] = (y32->data.F32[i] * x32->data.F32[i+1]/H); + // From (2) + ((spline->spline[i])->coeff[0])-= ((y32->data.F32[i+1] * x32->data.F32[i])/H); + // From (3) + tmp = (x32->data.F32[i+1] * x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H); + tmp-= (x32->data.F32[i+1] / H); + tmp*= (spline->p_psDeriv2)[i] * H * H / 6.0; + ((spline->spline[i])->coeff[0])+= tmp; + // From (4) + tmp = -(x32->data.F32[i] * x32->data.F32[i] * x32->data.F32[i]) / (H * H * H); + tmp+= (x32->data.F32[i] / H); + tmp*= (spline->p_psDeriv2)[i+1] * H * H / 6.0; + ((spline->spline[i])->coeff[0])+= tmp; + + // + // ******** Calculate 1-order term ******** + // + // From (1) + (spline->spline[i])->coeff[1] = -(y32->data.F32[i]) / H; + // From (2) + ((spline->spline[i])->coeff[1])+= (y32->data.F32[i+1] / H); + // From (3) + tmp = -3.0 * (x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H); + tmp+= (1.0 / H); + tmp*= ((spline->p_psDeriv2)[i]) * H * H / 6.0; + ((spline->spline[i])->coeff[1])+= tmp; + // From (4) + tmp = 3.0 * (x32->data.F32[i] * x32->data.F32[i]) / (H * H * H); + tmp-= (1.0 / H); + tmp*= ((spline->p_psDeriv2)[i+1]) * H * H / 6.0; + ((spline->spline[i])->coeff[1])+= tmp; + + // + // ******** Calculate 2-order term ******** + // + // From (3) + (spline->spline[i])->coeff[2] = ((spline->p_psDeriv2)[i]) * 3.0 * x32->data.F32[i+1] / (6.0 * H); + // From (4) + ((spline->spline[i])->coeff[2])-= (((spline->p_psDeriv2)[i+1]) * 3.0 * x32->data.F32[i] / (6.0 * H)); + + // + // ******** Calculate 3-order term ******** + // + // From (3) + (spline->spline[i])->coeff[3] = -((spline->p_psDeriv2)[i]) / (6.0 * H); + // From (4) + ((spline->spline[i])->coeff[3])+= ((spline->p_psDeriv2)[i+1]) / (6.0 * H); + + psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, + "(spline->spline[%u])->coeff[0] is %f\n", i, (spline->spline[i])->coeff[0]); + psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, + "(spline->spline[%u])->coeff[1] is %f\n", i, (spline->spline[i])->coeff[1]); + psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, + "(spline->spline[%u])->coeff[2] is %f\n", i, (spline->spline[i])->coeff[2]); + psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, + "(spline->spline[%u])->coeff[3] is %f\n", i, (spline->spline[i])->coeff[3]); + + } + + psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, + "---- psVectorFitSpline1D() end ----\n"); + return(spline); +} + + + + + + + + + +/***************************************************************************** +psVectorFitSpline1DNEW(): given a psSpline1D data structure and a set of x/y + +xF32 and yF32 are internal psVectors which are used to hold the psF32 versions +of the input data, if necessary. xPtr and yPtr are pointers to either xF32 or +the x argument. All computation is done on xPtr and yPtr. xF32 and yF32 will +simply be psFree() at the end. + +XXX: nKnots makes no sense. This number is always equal to the size of the x +an y vectors. + +XXX: How do we specify the spline order? For now, order=3. + *****************************************************************************/ +#define PS_XXX_SPLINE_ORDER 3 +psSpline1D *psVectorFitSpline1DNEW(const psVector* x, ///< Ordinates (or NULL to just use the indices) + const psVector* y, ///< Coordinates + unsigned int nKnots) +{ + psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, "---- psVectorFitSpline1DNEW() begin ----\n"); + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL); + // PS_ASSERT_INT_EQUAL(y->n, nKnots); + + // + // The following code ensures that xPtr points to a psF32 version of the + // ordinate data. + // + psVector *xF32 = NULL; + psVector *xPtr = NULL; + if (x != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL); + if (PS_TYPE_F64 == x->type.type) { + xF32 = psVectorAlloc(y->n, PS_TYPE_F32); + for (unsigned int i = 0 ; i < x->n ; i++) { + xF32->data.F32[i] = (psF32) x->data.F64[i]; + } + xPtr = xF32; + } else if (PS_TYPE_F32 == x->type.type) { + xPtr = (psVector *) x; + } else { + psError(PS_ERR_UNKNOWN, true, "psVector x is wrong type.\n"); + return(NULL); + } + } else { + // If x==NULL, create an x32 vector with x values set to (0:n). + xF32 = psVectorAlloc(y->n, PS_TYPE_F32); + for (unsigned int i = 0 ; i < x->n ; i++) { + xF32->data.F32[i] = (psF32) i; + } + xPtr = xF32; + } + + // + // If y is of type psF64, then create a new vector yF32 and convert the + // y elements. Regardless of y's type, we create a yPtr which will be + // used in the remainder of this function. + // + psVector *yF32 = NULL; + psVector *yPtr = NULL; + if (PS_TYPE_F64 == y->type.type) { + yF32 = psVectorAlloc(y->n, PS_TYPE_F32); + for (unsigned int i = 0 ; i < y->n ; i++) { + yF32->data.F32[i] = (psF32) y->data.F64[i]; + } + yPtr = yF32; + } else { + yPtr = (psVector *) y; + } + + psSpline1D *mySpline = psSpline1DAllocGeneric(xPtr, PS_XXX_SPLINE_ORDER); + + unsigned int numSplines = nKnots - 1; + psF32 tmp; + psF32 H; + unsigned int i; + psF32 slope; + // XXX: get rid of x32 and y32 (this is from old code) + psVector *x32 = xPtr; + psVector *y32 = yPtr; + + // If these are linear splines, which means their polynomials will have + // two coefficients, then we do the simple calculation. + if (1 == PS_XXX_SPLINE_ORDER) { for (i=0;in;i++) { slope = (y32->data.F32[i+1] - y32->data.F32[i]) / @@ -391,5 +591,5 @@ (mySpline->spline[i])->coeff[1] = slope; psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4, - "---- mySpline %d coeffs are (%f, %f)\n", i, + "---- mySpline %u coeffs are (%f, %f)\n", i, (mySpline->spline[i])->coeff[0], (mySpline->spline[i])->coeff[1]); @@ -400,22 +600,27 @@ } + // // Check if these are cubic splines (n==4). If not, psError. - if (4 != (mySpline->spline[0])->n) { + // + if (3 != PS_XXX_SPLINE_ORDER) { psError(PS_ERR_BAD_PARAMETER_SIZE, true, - "Don't know how to generate %d-order splines.", - (mySpline->spline[0])->n-1); + "Don't know how to generate %d-order splines.", PS_XXX_SPLINE_ORDER); return(NULL); } + // // If we get here, then we know these are cubic splines. We first // generate the second derivatives at each data point. + // mySpline->p_psDeriv2 = calculateSecondDerivs(x32, y32); for (i=0;in;i++) psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, - "Second deriv[%d] is %f\n", i, mySpline->p_psDeriv2[i]); - + "Second deriv[%u] is %f\n", i, mySpline->p_psDeriv2[i]); + + // // We generate the coefficients of the spline polynomials. I can't // concisely explain how this code works. See above function comments // and Numerical Recipes in C. + // for (i=0;idata.F32[i+1] - x32->data.F32[i]; @@ -477,216 +682,11 @@ psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, - "(mySpline->spline[%d])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]); + "(mySpline->spline[%u])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]); psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, - "(mySpline->spline[%d])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]); + "(mySpline->spline[%u])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]); psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, - "(mySpline->spline[%d])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]); + "(mySpline->spline[%u])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]); psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, - "(mySpline->spline[%d])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]); - - } - - psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, - "---- psVectorFitSpline1D() end ----\n"); - return(mySpline); -} - - - - - - - - - -/***************************************************************************** -psVectorFitSpline1DNEW(): given a psSpline1D data structure and a set of x/y - -xF32 and yF32 are internal psVectors which are used to hold the psF32 versions -of the input data, if necessary. xPtr and yPtr are pointers to either xF32 or -the x argument. All computation is done on xPtr and yPtr. xF32 and yF32 will -simply be psFree() at the end. - -XXX: nKnots makes no sense. This number is always equal to the size of the x -an y vectors. - -XXX: How do we specify the spline order? For now, order=3. - *****************************************************************************/ -#define PS_XXX_SPLINE_ORDER 3 -psSpline1D *psVectorFitSpline1DNEW(const psVector* x, ///< Ordinates (or NULL to just use the indices) - const psVector* y, ///< Coordinates - int nKnots) -{ - psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, "---- psVectorFitSpline1DNEW() begin ----\n"); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL); - // PS_ASSERT_INT_EQUAL(y->n, nKnots); - - // - // The following code ensures that xPtr points to a psF32 version of the - // ordinate data. - // - psVector *xF32 = NULL; - psVector *xPtr = NULL; - if (x != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL); - if (PS_TYPE_F64 == x->type.type) { - xF32 = psVectorAlloc(y->n, PS_TYPE_F32); - for (psS32 i = 0 ; i < x->n ; i++) { - xF32->data.F32[i] = (psF32) x->data.F64[i]; - } - xPtr = xF32; - } else if (PS_TYPE_F32 == x->type.type) { - xPtr = (psVector *) x; - } else { - psError(PS_ERR_UNKNOWN, true, "psVector x is wrong type.\n"); - return(NULL); - } - } else { - // If x==NULL, create an x32 vector with x values set to (0:n). - xF32 = psVectorAlloc(y->n, PS_TYPE_F32); - for (psS32 i = 0 ; i < x->n ; i++) { - xF32->data.F32[i] = (psF32) i; - } - xPtr = xF32; - } - - // - // If y is of type psF64, then create a new vector yF32 and convert the - // y elements. Regardless of y's type, we create a yPtr which will be - // used in the remainder of this function. - // - psVector *yF32 = NULL; - psVector *yPtr = NULL; - if (PS_TYPE_F64 == y->type.type) { - yF32 = psVectorAlloc(y->n, PS_TYPE_F32); - for (psS32 i = 0 ; i < y->n ; i++) { - yF32->data.F32[i] = (psF32) y->data.F64[i]; - } - yPtr = yF32; - } else { - yPtr = (psVector *) y; - } - - psSpline1D *mySpline = psSpline1DAllocGeneric(xPtr, PS_XXX_SPLINE_ORDER); - - psS32 numSplines = nKnots - 1; - psF32 tmp; - psF32 H; - psS32 i; - psF32 slope; - // XXX: get rid of x32 and y32 (this is from old code) - psVector *x32 = xPtr; - psVector *y32 = yPtr; - - // If these are linear splines, which means their polynomials will have - // two coefficients, then we do the simple calculation. - if (1 == PS_XXX_SPLINE_ORDER) { - for (i=0;in;i++) { - slope = (y32->data.F32[i+1] - y32->data.F32[i]) / - (mySpline->knots->data.F32[i+1] - mySpline->knots->data.F32[i]); - (mySpline->spline[i])->coeff[0] = y32->data.F32[i] - - (slope * mySpline->knots->data.F32[i]); - - (mySpline->spline[i])->coeff[1] = slope; - psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4, - "---- mySpline %d coeffs are (%f, %f)\n", i, - (mySpline->spline[i])->coeff[0], - (mySpline->spline[i])->coeff[1]); - } - psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4, - "---- Exiting psVectorFitSpline1D()()\n"); - return((psSpline1D *) mySpline); - } - - // - // Check if these are cubic splines (n==4). If not, psError. - // - if (3 != PS_XXX_SPLINE_ORDER) { - psError(PS_ERR_BAD_PARAMETER_SIZE, true, - "Don't know how to generate %d-order splines.", PS_XXX_SPLINE_ORDER); - return(NULL); - } - - // - // If we get here, then we know these are cubic splines. We first - // generate the second derivatives at each data point. - // - mySpline->p_psDeriv2 = calculateSecondDerivs(x32, y32); - for (i=0;in;i++) - psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, - "Second deriv[%d] is %f\n", i, mySpline->p_psDeriv2[i]); - - // - // We generate the coefficients of the spline polynomials. I can't - // concisely explain how this code works. See above function comments - // and Numerical Recipes in C. - // - for (i=0;idata.F32[i+1] - x32->data.F32[i]; - psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, - "x data (%f - %f) (%f)\n", - x32->data.F32[i], - x32->data.F32[i+1], H); - // - // ******** Calculate 0-order term ******** - // - // From (1) - (mySpline->spline[i])->coeff[0] = (y32->data.F32[i] * x32->data.F32[i+1]/H); - // From (2) - ((mySpline->spline[i])->coeff[0])-= ((y32->data.F32[i+1] * x32->data.F32[i])/H); - // From (3) - tmp = (x32->data.F32[i+1] * x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H); - tmp-= (x32->data.F32[i+1] / H); - tmp*= (mySpline->p_psDeriv2)[i] * H * H / 6.0; - ((mySpline->spline[i])->coeff[0])+= tmp; - // From (4) - tmp = -(x32->data.F32[i] * x32->data.F32[i] * x32->data.F32[i]) / (H * H * H); - tmp+= (x32->data.F32[i] / H); - tmp*= (mySpline->p_psDeriv2)[i+1] * H * H / 6.0; - ((mySpline->spline[i])->coeff[0])+= tmp; - - // - // ******** Calculate 1-order term ******** - // - // From (1) - (mySpline->spline[i])->coeff[1] = -(y32->data.F32[i]) / H; - // From (2) - ((mySpline->spline[i])->coeff[1])+= (y32->data.F32[i+1] / H); - // From (3) - tmp = -3.0 * (x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H); - tmp+= (1.0 / H); - tmp*= ((mySpline->p_psDeriv2)[i]) * H * H / 6.0; - ((mySpline->spline[i])->coeff[1])+= tmp; - // From (4) - tmp = 3.0 * (x32->data.F32[i] * x32->data.F32[i]) / (H * H * H); - tmp-= (1.0 / H); - tmp*= ((mySpline->p_psDeriv2)[i+1]) * H * H / 6.0; - ((mySpline->spline[i])->coeff[1])+= tmp; - - // - // ******** Calculate 2-order term ******** - // - // From (3) - (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x32->data.F32[i+1] / (6.0 * H); - // From (4) - ((mySpline->spline[i])->coeff[2])-= (((mySpline->p_psDeriv2)[i+1]) * 3.0 * x32->data.F32[i] / (6.0 * H)); - - // - // ******** Calculate 3-order term ******** - // - // From (3) - (mySpline->spline[i])->coeff[3] = -((mySpline->p_psDeriv2)[i]) / (6.0 * H); - // From (4) - ((mySpline->spline[i])->coeff[3])+= ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H); - - psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, - "(mySpline->spline[%d])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]); - psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, - "(mySpline->spline[%d])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]); - psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, - "(mySpline->spline[%d])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]); - psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, - "(mySpline->spline[%d])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]); + "(mySpline->spline[%u])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]); } @@ -726,6 +726,6 @@ XXX: What should be the default type for knots be? psF32 is assumed. *****************************************************************************/ -psSpline1D *psSpline1DAlloc(int numSplines, - int order, +psSpline1D *psSpline1DAlloc(unsigned int numSplines, + unsigned int order, float min, float max) @@ -743,5 +743,5 @@ // tmpSpline->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); - for (psS32 i=0;ispline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD); } @@ -756,5 +756,5 @@ psF32 width = (max - min) / ((psF32) numSplines); tmpSpline->knots->data.F32[0] = min; - for (psS32 i=1;iknots->data.F32[i] = min + (width * (psF32) i); } @@ -773,9 +773,9 @@ if (in->type.type == PS_TYPE_F32) { - for (psS32 i = 0 ; i < in->n ; i++) { + for (unsigned int i = 0 ; i < in->n ; i++) { out->data.F32[i] = in->data.F32[i]; } } else if (in->type.type == PS_TYPE_F64) { - for (psS32 i = 0 ; i < in->n ; i++) { + for (unsigned int i = 0 ; i < in->n ; i++) { out->data.F64[i] = in->data.F64[i]; } @@ -791,5 +791,5 @@ *****************************************************************************/ psSpline1D *psSpline1DAllocGeneric(const psVector *bounds, - int order) + unsigned int order) { PS_ASSERT_VECTOR_NON_NULL(bounds, NULL); @@ -799,5 +799,5 @@ psSpline1D *tmpSpline = (psSpline1D *) psAlloc(sizeof(psSpline1D)); - psS32 numSplines = bounds->n - 1; + unsigned int numSplines = bounds->n - 1; tmpSpline->n = numSplines; @@ -807,5 +807,5 @@ // tmpSpline->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); - for (psS32 i=0;ispline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD); } @@ -818,5 +818,5 @@ // XXX:Ensure that the knots are monotonic. // - for (psS32 i=0;in-1;i++) { + for (unsigned int i=0;in-1;i++) { if (FLT_EPSILON >= fabs(bounds->data.F32[i+1]-bounds->data.F32[i])) { psError(PS_ERR_UNKNOWN, true, "data points must be distinct ([%d] %f %f)\n", i, bounds->data.F32[i], bounds->data.F32[i+1]); @@ -839,7 +839,7 @@ *****************************************************************************/ #define FUNC_MACRO_VECTOR_BIN_DISECT(TYPE) \ -static psS32 vectorBinDisect##TYPE(ps##TYPE *bins, \ - psS32 numBins, \ - ps##TYPE x) \ +static unsigned int vectorBinDisect##TYPE(ps##TYPE *bins, \ + psS32 numBins, \ + ps##TYPE x) \ { \ psS32 min; \ @@ -906,6 +906,6 @@ XXX: Assert that the psVector and psScalar have the same type. *****************************************************************************/ -psS32 p_psVectorBinDisect(psVector *bins, - psScalar *x) +unsigned int p_psVectorBinDisect(psVector *bins, + psScalar *x) { PS_ASSERT_VECTOR_NON_NULL(bins, -4); @@ -972,5 +972,5 @@ psScalar *p_psVectorInterpolate(psVector *domain, psVector *range, - int order, + unsigned int order, psScalar *x) {