Changeset 6305 for trunk/psLib/src/math/psSpline.c
- Timestamp:
- Feb 2, 2006, 11:09:08 AM (20 years ago)
- File:
-
- 1 edited
-
trunk/psLib/src/math/psSpline.c (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/src/math/psSpline.c
r6204 r6305 4 4 * and evaluation routines. 5 5 * 6 * This file will hold the functions for allocated, freeing, and evaluating 7 * splines. 6 * This file contains the routines that allocate, free, and evaluate splines. 8 7 * 9 * @version $Revision: 1.13 4$ $Name: not supported by cvs2svn $10 * @date $Date: 2006-0 1-26 21:10:22$8 * @version $Revision: 1.135 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2006-02-02 21:09:08 $ 11 10 * 12 11 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii … … 31 30 #include "psConstants.h" 32 31 #include "psErrorText.h" 32 #include "psMathUtils.h" 33 33 34 34 /*****************************************************************************/ … … 72 72 73 73 return; 74 } 75 76 77 static void PS_POLY1D_PRINT( 78 psPolynomial1D *poly) 79 { 80 printf("-------------- PS_POLY1D_PRINT() --------------\n"); 81 printf("poly->nX is %d\n", poly->nX); 82 for (psS32 i = 0 ; i < (1 + poly->nX) ; i++) { 83 printf("poly->coeff[%d] is %f\n", i, poly->coeff[i]); 84 } 85 } 86 87 static void PS_PRINT_SPLINE2(psSpline1D *mySpline) 88 { 89 printf("-------------- PS_PRINT_SPLINE2() --------------\n"); 90 printf("mySpline->n is %d\n", mySpline->n); 91 for (psS32 i = 0 ; i < mySpline->n ; i++) { 92 PS_POLY1D_PRINT(mySpline->spline[i]); 93 } 94 PS_VECTOR_PRINT_F32(mySpline->knots); 74 95 } 75 96 … … 139 160 psTrace(__func__, 4, "---- %s() end ----\n", __func__); 140 161 return(derivs2); 141 }142 143 144 /*****************************************************************************145 vectorBinDisectTYPE(): This is a macro for a private function which takes as146 input a vector an array of data as well as a single value for that data. The147 input vector values are assumed to be non-decreasing (v[i-1] <= v[i] for all148 i). This routine does a binary disection of the vector and returns "i" such149 that (v[i] <= x <= v[i+1). If x lies outside the range of v[], then this150 routine prints a warning message and returns (-2 or -1).151 *****************************************************************************/152 #define FUNC_MACRO_VECTOR_BIN_DISECT(TYPE) \153 static psS32 vectorBinDisect##TYPE(ps##TYPE *bins, \154 psS32 numBins, \155 ps##TYPE x) \156 { \157 psS32 min; \158 psS32 max; \159 psS32 mid; \160 psTrace(__func__, 4, "---- %s() begin ----\n", __func__); \161 if (x < bins[0]) { \162 psLogMsg(__func__, PS_LOG_WARN, \163 "vectorBinDisect%s(): ordinate %f is outside vector range (%f - %f).", \164 #TYPE, x, bins[0], bins[numBins-1]); \165 return(-2); \166 } \167 if (x > bins[numBins-1]) { \168 psLogMsg(__func__, PS_LOG_WARN, \169 "vectorBinDisect%s(): ordinate %f is outside vector range (%f - %f).", \170 #TYPE, x, bins[0], bins[numBins-1]); \171 return(-1); \172 } \173 \174 min = 0; \175 max = numBins-2; \176 mid = ((max+1)-min)/2; \177 while (min != max) { \178 psTrace(__func__, 6, "(min, mid, max) is (%u, %u, %u): (x, bins) is (%f, %f)\n", min, mid, max, x, bins[mid]); \179 \180 if (x == bins[mid]) { \181 psTrace(__func__, 4, "---- %s(%d) end ----\n", __func__, mid); \182 return(mid); \183 } else if (x < bins[mid]) { \184 max = mid-1; \185 } else { \186 min = mid; \187 } \188 mid = ((max+1)+min)/2; \189 } \190 psTrace(__func__, 4, "---- %s(%d) end ----\n", __func__, min); \191 return(min); \192 } \193 194 FUNC_MACRO_VECTOR_BIN_DISECT(S8)195 FUNC_MACRO_VECTOR_BIN_DISECT(S16)196 FUNC_MACRO_VECTOR_BIN_DISECT(S32)197 FUNC_MACRO_VECTOR_BIN_DISECT(S64)198 FUNC_MACRO_VECTOR_BIN_DISECT(U8)199 FUNC_MACRO_VECTOR_BIN_DISECT(U16)200 FUNC_MACRO_VECTOR_BIN_DISECT(U32)201 FUNC_MACRO_VECTOR_BIN_DISECT(U64)202 FUNC_MACRO_VECTOR_BIN_DISECT(F32)203 FUNC_MACRO_VECTOR_BIN_DISECT(F64)204 205 /*****************************************************************************206 p_psVectorBinDisect(): A wrapper to the above p_psVectorBinDisect().207 *****************************************************************************/208 psS32 p_psVectorBinDisect(209 psVector *bins,210 psScalar *x)211 {212 PS_ASSERT_VECTOR_NON_NULL(bins, -4);213 PS_ASSERT_VECTOR_NON_EMPTY(bins, -4);214 PS_ASSERT_PTR_NON_NULL(x, -6);215 PS_ASSERT_PTR_TYPE_EQUAL(x, bins, -3);216 char* strType;217 218 switch (x->type.type) {219 case PS_TYPE_U8:220 return(vectorBinDisectU8(bins->data.U8, bins->n, x->data.U8));221 case PS_TYPE_U16:222 return(vectorBinDisectU16(bins->data.U16, bins->n, x->data.U16));223 case PS_TYPE_U32:224 return(vectorBinDisectU32(bins->data.U32, bins->n, x->data.U32));225 case PS_TYPE_U64:226 return(vectorBinDisectU64(bins->data.U64, bins->n, x->data.U64));227 case PS_TYPE_S8:228 return(vectorBinDisectS8(bins->data.S8, bins->n, x->data.S8));229 case PS_TYPE_S16:230 return(vectorBinDisectS16(bins->data.S16, bins->n, x->data.S16));231 case PS_TYPE_S32:232 return(vectorBinDisectS32(bins->data.S32, bins->n, x->data.S32));233 case PS_TYPE_S64:234 return(vectorBinDisectS64(bins->data.S64, bins->n, x->data.S64));235 case PS_TYPE_F32:236 return(vectorBinDisectF32(bins->data.F32, bins->n, x->data.F32));237 case PS_TYPE_F64:238 return(vectorBinDisectF64(bins->data.F64, bins->n, x->data.F64));239 case PS_TYPE_C32:240 PS_TYPE_NAME(strType,x->type.type);241 psError(PS_ERR_BAD_PARAMETER_TYPE,242 PS_ERRORTEXT_psSpline_TYPE_NOT_SUPPORTED,243 strType);244 return 0;245 case PS_TYPE_C64:246 PS_TYPE_NAME(strType,x->type.type);247 psError(PS_ERR_BAD_PARAMETER_TYPE,248 PS_ERRORTEXT_psSpline_TYPE_NOT_SUPPORTED,249 strType);250 return 0;251 case PS_TYPE_BOOL:252 PS_TYPE_NAME(strType,x->type.type);253 psError(PS_ERR_BAD_PARAMETER_TYPE,254 PS_ERRORTEXT_psSpline_TYPE_NOT_SUPPORTED,255 strType);256 return 0;257 }258 return(-3);259 }260 261 /*****************************************************************************262 fullInterpolate1DF32(): This routine will take as input n-element floating263 point arrays domain and range, and the x value, assumed to lie with the264 domain vector. It produces as output the (n-1)-order LaGrange interpolated265 value of x.266 267 XXX: do we error check for non-distinct domain values?268 269 XXX: This is a somewhat general function. There's no reason to put it in this270 file.271 272 XXX: Sort all these interpolation functions out. Why are there so many?273 274 XXX: The fullInterpolate1D macros are not used anywhere.275 *****************************************************************************/276 #define FUNC_MACRO_FULL_INTERPOLATE_1D(TYPE) \277 static psF32 fullInterpolate1D##TYPE(ps##TYPE *domain, \278 ps##TYPE *range, \279 unsigned int n, \280 ps##TYPE x) \281 { \282 \283 unsigned int i; \284 unsigned int m; \285 static psVector *p = NULL; \286 p = psVectorRecycle(p, n, PS_TYPE_##TYPE); \287 p_psMemSetPersistent(p, true); \288 p_psMemSetPersistent(p->data.TYPE, true); \289 \290 psTrace(__func__, 4, "---- %s() begin %u-order at x=%f) (%d data points) ----\n", __func__, n-1, x, n); \291 for (i=0;i<n;i++) { \292 psTrace(__func__, 6, "domain/range is (%f %f)\n", domain[i], range[i]); \293 } \294 \295 for (i=0;i<n;i++) { \296 p->data.TYPE[i] = range[i]; \297 psTrace(__func__, 6, "p->data.TYPE[%u] is %f\n", i, p->data.TYPE[i]); \298 \299 } \300 \301 /* From NR, during each iteration of the m loop, we are computing the \302 p_{i ... i+m} terms. \303 */ \304 for (m=1;m<n;m++) { \305 for (i=0;i<n-m;i++) { \306 /* From NR: we are computing P_{i ... i+m} \307 */ \308 p->data.TYPE[i] = (((x-domain[i+m]) * p->data.TYPE[i]) + \309 ((domain[i]-x) * p->data.TYPE[i+1])) / \310 (domain[i] - domain[i+m]); \311 psTrace(__func__, 6, "p->data.TYPE[%u] is %f\n", i, p->data.TYPE[i]); \312 } \313 } \314 psTrace(__func__, 4, "---- %s(....) end ----\n", __func__); \315 return(p->data.TYPE[0]); \316 } \317 318 /* XXX: Do this correctly.319 FUNC_MACRO_FULL_INTERPOLATE_1D(U8)320 FUNC_MACRO_FULL_INTERPOLATE_1D(U16)321 FUNC_MACRO_FULL_INTERPOLATE_1D(U32)322 FUNC_MACRO_FULL_INTERPOLATE_1D(U64)323 FUNC_MACRO_FULL_INTERPOLATE_1D(S8)324 FUNC_MACRO_FULL_INTERPOLATE_1D(S16)325 FUNC_MACRO_FULL_INTERPOLATE_1D(S32)326 FUNC_MACRO_FULL_INTERPOLATE_1D(S64)327 FUNC_MACRO_FULL_INTERPOLATE_1D(F64)328 */329 FUNC_MACRO_FULL_INTERPOLATE_1D(F32)330 331 332 /*****************************************************************************333 interpolate1DF32(): this is the base 1-D flat memory routine to perform334 LaGrange interpolation.335 *****************************************************************************/336 static psF32 interpolate1DF32(337 psF32 *domain,338 psF32 *range,339 unsigned int n,340 unsigned int order,341 psF32 x)342 {343 psTrace(__func__, 4, "---- %s() begin ----\n", __func__);344 PS_ASSERT_PTR_NON_NULL(domain, NAN)345 PS_ASSERT_PTR_NON_NULL(range, NAN)346 // XXX: Check valid values for n, order, and x?347 348 psS32 binNum;349 psS32 numIntPoints = order+1;350 psS32 origin;351 352 binNum = vectorBinDisectF32(domain, n, x);353 354 if (0 == numIntPoints%2) {355 origin = binNum - ((numIntPoints/2) - 1);356 } else {357 origin = binNum - (numIntPoints/2);358 if ((x-domain[binNum]) > (domain[binNum+1]-x)) {359 // x is closer to binNum+1.360 origin = 1 + (binNum - (numIntPoints/2));361 }362 }363 if (origin < 0) {364 origin = 0;365 }366 if ((origin + numIntPoints) > n) {367 origin = n - numIntPoints;368 }369 370 psTrace(__func__, 4, "---- %s(....) end ----\n", __func__);371 return(fullInterpolate1DF32(&domain[origin], &range[origin], order+1, x));372 }373 374 375 /*****************************************************************************376 p_psVectorInterpolate(): This routine will take as input psVectors domain and377 range, and the x value, assumed to lie with the domain vector. It produces378 as output the LaGrange interpolated value of a polynomial of the specified379 order around the point x.380 381 XXX: This stuff does not currently work with a mask.382 383 XXX: add another psScalar argument for the result (so you don't have to384 allocate it each time).385 386 XXX: The VectorCopy routines seg fault when I declare range32 as static.387 *****************************************************************************/388 psScalar *p_psVectorInterpolate(389 psVector *domain,390 psVector *range,391 psS32 order,392 psScalar *x)393 {394 psTrace(__func__, 4, "---- %s() begin ----\n", __func__);395 PS_ASSERT_VECTOR_NON_NULL(domain, NULL);396 PS_ASSERT_VECTOR_NON_NULL(range, NULL);397 PS_ASSERT_PTR_NON_NULL(x, NULL);398 PS_ASSERT_INT_NONNEGATIVE(order, NULL);399 PS_ASSERT_VECTORS_SIZE_EQUAL(domain, range, NULL);400 PS_ASSERT_PTR_TYPE_EQUAL(domain, range, NULL);401 PS_ASSERT_PTR_TYPE_EQUAL(domain, x, NULL);402 403 psVector *range32 = NULL;404 psVector *domain32 = NULL;405 if (order > (domain->n - 1)) {406 psError(PS_ERR_BAD_PARAMETER_SIZE, true,407 PS_ERRORTEXT_psSpline_NOT_ENOUGH_DATAPOINTS,408 order);409 return(NULL);410 }411 412 if (x->type.type == PS_TYPE_F32) {413 psTrace(__func__, 4, "---- %s(NULL) end ----\n", __func__);414 return(psScalarAlloc(interpolate1DF32(domain->data.F32,415 range->data.F32,416 domain->n,417 order,418 x->data.F32), PS_TYPE_F32));419 } else if (x->type.type == PS_TYPE_F64) {420 // XXX: use recycled vectors here.421 range32 = psVectorCopy(range32, range, PS_TYPE_F32);422 domain32 = psVectorCopy(domain32, domain, PS_TYPE_F32);423 424 psScalar *tmpScalar = psScalarAlloc((psF64)425 interpolate1DF32(domain32->data.F32,426 range32->data.F32,427 domain32->n,428 order,429 (psF32) x->data.F64), PS_TYPE_F64);430 psFree(range32);431 psFree(domain32);432 433 psTrace(__func__, 4, "---- %s() end ----\n", __func__);434 // XXX: Convert data type to F64?435 return(tmpScalar);436 437 } else {438 char* strType;439 PS_TYPE_NAME(strType,x->type.type);440 psError(PS_ERR_BAD_PARAMETER_TYPE,441 PS_ERRORTEXT_psSpline_TYPE_NOT_SUPPORTED,442 strType);443 }444 445 psTrace(__func__, 4, "---- %s(NULL) end ----\n", __func__);446 return(NULL);447 162 } 448 163 … … 627 342 } 628 343 629 void PS_POLY1D_PRINT(630 psPolynomial1D *poly)631 {632 printf("-------------- PS_POLY1D_PRINT() --------------\n");633 printf("poly->nX is %d\n", poly->nX);634 for (psS32 i = 0 ; i < (1 + poly->nX) ; i++) {635 printf("poly->coeff[%d] is %f\n", i, poly->coeff[i]);636 }637 }638 639 void PS_PRINT_SPLINE2(psSpline1D *mySpline)640 {641 printf("-------------- PS_PRINT_SPLINE2() --------------\n");642 printf("mySpline->n is %d\n", mySpline->n);643 for (psS32 i = 0 ; i < mySpline->n ; i++) {644 PS_POLY1D_PRINT(mySpline->spline[i]);645 }646 PS_VECTOR_PRINT_F32(mySpline->knots);647 }648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 344 685 345 /***************************************************************************** … … 706 366 707 367 psS32 n = spline->n; 708 psS32 binNum = vectorBinDisectF32(spline->knots->data.F32, (spline->n)+1, x); 368 psScalar tmpScalar; 369 tmpScalar.type.type = PS_TYPE_F32; 370 tmpScalar.data.F32 = x; 371 psS32 binNum = p_psVectorBinDisect(spline->knots, &tmpScalar); 372 709 373 if (binNum < 0) { 710 374 // … … 746 410 if (psTraceGetLevel(__func__) >= 6) { 747 411 PS_VECTOR_PRINT_F32(x); 748 // XXX: print spline412 PS_PRINT_SPLINE2((psSpline1D *) spline); 749 413 } 750 414
Note:
See TracChangeset
for help on using the changeset viewer.
