Index: trunk/psLib/src/math/psFunctions.c =================================================================== --- trunk/psLib/src/math/psFunctions.c (revision 4567) +++ trunk/psLib/src/math/psFunctions.c (revision 4568) @@ -7,6 +7,6 @@ * polynomials. It also contains a Gaussian functions. * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-07-15 23:19:41 $ + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-07-16 00:06:32 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -1818,88 +1818,123 @@ must exist n+1 points in "knots". -XXX: Ensure that domain[i+1] != domain[i] +XXX: Is this really needed anymore? -XXX: What should be the defualty type for knots be? psF32 is assumed. +XXX: Ensure that knots[i+1] != knots[i] + +XXX: What should be the default type for knots be? psF32 is assumed. *****************************************************************************/ -psSpline1D *psSpline1DAlloc( int numSplines, - int order, - float min, - float max) -{ +psSpline1D *psSpline1DAlloc(int numSplines, + int order, + float min, + float max) +{ + printf("HEY: psSpline1DAlloc()\n"); PS_ASSERT_INT_NONNEGATIVE(numSplines, NULL); PS_ASSERT_INT_NONNEGATIVE(order, NULL); PS_ASSERT_FLOAT_NON_EQUAL(max, min, NULL); - psSpline1D *tmp = NULL; - psS32 i; - psF32 tmpDomain; - psF32 width; - - tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D)); - tmp->n = numSplines; - - tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); - for (i=0;ispline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD); - } - - // This should be set by the psVectorFitSpline1D() - tmp->p_psDeriv2 = NULL; - - tmp->knots = psVectorAlloc(numSplines+1, PS_TYPE_F32); - width = (max - min) / ((psF32) numSplines); - - tmp->knots->data.F32[0] = min; - tmpDomain = min+width; - for (i=1;iknots->data.F32[i] = tmpDomain; - tmpDomain+= width; - } - tmp->knots->data.F32[numSplines] = max; - - psMemSetDeallocator(tmp,(psFreeFunc)spline1DFree); - return(tmp); -} - + psSpline1D *tmpSpline = (psSpline1D *) psAlloc(sizeof(psSpline1D)); + tmpSpline->n = numSplines; + + // + // XXX: We might have to allocate single or double polynomials depending on the type + // of the psVector bounds. For now, all knots and spline polynomials are 32-bit. + // + tmpSpline->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); + for (psS32 i=0;ispline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD); + } + + // This will be computed by psVectorFitSpline1D() + tmpSpline->p_psDeriv2 = NULL; + + // + // XXX:Ensure that the knots are distinct, and monotonic. + // + tmpSpline->knots = psVectorAlloc(numSplines+1, PS_TYPE_F32); + if (tmpSpline->knots == NULL) { + printf("BAD\n"); + } else { + printf("GOOD\n"); + } + psF32 width = (max - min) / ((psF32) numSplines); + tmpSpline->knots->data.F32[0] = min; + for (psS32 i=1;iknots->data.F32[i] = min + (width * (psF32) i); + } + tmpSpline->knots->data.F32[numSplines] = max; + + psMemSetDeallocator(tmpSpline, (psFreeFunc)spline1DFree); + return(tmpSpline); +} /***************************************************************************** -XXX: What should be the defualty type for knots be? psF32 is assumed. +XXX: Is there a psLib function for this? + *****************************************************************************/ +psVector *PsVectorDup(psVector *in) +{ + psVector *out = psVectorAlloc(in->n, in->type.type); + + if (in->type.type == PS_TYPE_F32) { + for (psS32 i = 0 ; i < in->n ; i++) { + out->data.F32[i] = in->data.F32[i]; + } + } else if (in->type.type == PS_TYPE_F64) { + for (psS32 i = 0 ; i < in->n ; i++) { + out->data.F64[i] = in->data.F64[i]; + } + } else { + printf("XXX: Generate an error here.\n"); + return(NULL); + } + return(out); +} + +/***************************************************************************** +XXX: What should be the default type for knots, spline polys? psF32 is assumed. *****************************************************************************/ psSpline1D *psSpline1DAllocGeneric(const psVector *bounds, int order) { + printf("HEY: psSpline1DAllocGeneric()\n"); PS_ASSERT_VECTOR_NON_NULL(bounds, NULL); PS_ASSERT_VECTOR_NON_EMPTY(bounds, NULL); + PS_ASSERT_VECTOR_TYPE(bounds, PS_TYPE_F32, NULL); PS_ASSERT_INT_NONNEGATIVE(order, NULL); - psSpline1D *tmpSpline = NULL; - unsigned int i; - unsigned int numSplines; - - tmpSpline = (psSpline1D *) psAlloc(sizeof(psSpline1D)); - - numSplines = bounds->n - 1; + psSpline1D *tmpSpline = (psSpline1D *) psAlloc(sizeof(psSpline1D)); + psS32 numSplines = bounds->n - 1; tmpSpline->n = numSplines; + // + // XXX: We might have to allocate single or double polynomials depending on the type + // of the psVector bounds. For now, all knots and spline polynomials are 32-bit. + // tmpSpline->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); - for (i=0;ispline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD); } - // This should be set by the psVectorFitSpline1D() + // This will be computed by psVectorFitSpline1D() tmpSpline->p_psDeriv2 = NULL; - tmpSpline->knots = psVectorAlloc(bounds->n, PS_TYPE_F32); - - for (i=0;in;i++) { - tmpSpline->knots->data.F32[i] = bounds->data.F32[i]; - if (i<(bounds->n-1)) { - if (FLT_EPSILON >= fabs(bounds->data.F32[i+1]-bounds->data.F32[i])) { - psError(PS_ERR_UNKNOWN, true, "data points must be distinct\n"); - } - } - } - - psMemSetDeallocator(tmpSpline,(psFreeFunc)spline1DFree); + // + // Ensure that all knots are distinct. + // XXX:Ensure that the knots are monotonic. + // + for (psS32 i=0;in-1;i++) { + if (FLT_EPSILON >= fabs(bounds->data.F32[i+1]-bounds->data.F32[i])) { + psError(PS_ERR_UNKNOWN, true, "data points must be distinct\n"); + return(NULL); + } + } + tmpSpline->knots = PsVectorDup((psVector *) bounds); + if (tmpSpline->knots == NULL) { + printf("BAD\n"); + } else { + printf("GOOD\n"); + } + + psMemSetDeallocator(tmpSpline, (psFreeFunc)spline1DFree); return(tmpSpline); }