Changeset 1718
- Timestamp:
- Sep 7, 2004, 8:00:49 PM (22 years ago)
- Location:
- trunk/psLib/src
- Files:
-
- 5 edited
-
dataManip/psFunctions.c (modified) (34 diffs)
-
math/psPolynomial.c (modified) (34 diffs)
-
math/psSpline.c (modified) (34 diffs)
-
sys/psTrace.c (modified) (19 diffs)
-
sysUtils/psTrace.c (modified) (19 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/src/dataManip/psFunctions.c
r1476 r1718 1 2 1 /** @file psFunctions.c 3 2 * … … 8 7 * polynomials. It also contains a Gaussian functions. 9 8 * 10 * @version $Revision: 1.2 3$ $Name: not supported by cvs2svn $11 * @date $Date: 2004-0 8-11 20:07:45$9 * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ 10 * @date $Date: 2004-09-08 06:00:49 $ 12 11 * 13 12 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 32 31 #include "psAbort.h" 33 32 #include "psFunctions.h" 33 #include "psLogMsg.h" 34 34 35 35 #include <gsl/gsl_rng.h> … … 42 42 /*****************************************************************************/ 43 43 44 // None 44 #define PS_CONVERT_VECTOR_F32(IN, OUT) \ 45 if (IN->type.type != PS_TYPE_F32) { \ 46 psLogMsg(__func__, PS_LOG_WARN, \ 47 "Input vector has incorrect type (%d).\n", \ 48 IN->type); \ 49 OUT = psVectorAlloc(IN->n, PS_TYPE_F32); \ 50 \ 51 if (IN->type.type == PS_TYPE_F64) { \ 52 for (int i=0;i<IN->n;i++) { \ 53 OUT->data.F32[i] = (float) (IN->data.F64[i]); \ 54 } \ 55 }\ 56 else { \ 57 psAbort(__func__, "Wrong type.\n"); \ 58 } \ 59 } else { \ 60 OUT = (psVector *) IN; \ 61 } 62 63 #define PS_CONVERT_VECTOR_F64(IN, OUT) \ 64 if (IN->type.type != PS_TYPE_F64) { \ 65 psLogMsg(__func__, PS_LOG_WARN, \ 66 "Input vector has incorrect type (%d).\n", \ 67 IN->type); \ 68 OUT = psVectorAlloc(IN->n, PS_TYPE_F64); \ 69 \ 70 if (IN->type.type == PS_TYPE_F32) { \ 71 for (int i=0;i<IN->n;i++) { \ 72 OUT->data.F64[i] = (float) (IN->data.F32[i]); \ 73 } \ 74 }\ 75 else { \ 76 psAbort(__func__, "Wrong type.\n"); \ 77 } \ 78 } else { \ 79 OUT = (psVector *) IN; \ 80 } 81 82 45 83 46 84 /*****************************************************************************/ … … 107 145 } 108 146 109 return 147 return(tmp * exp(-((x - mean) * (x - mean)) / (2.0 * sigma * sigma))); 110 148 } 111 149 … … 138 176 139 177 // NOTE: Should I free r as well? 140 return (gauss); 141 } 178 return(gauss); 179 } 180 181 //XXX: remove this 182 //typedef enum { 183 // PS_POLYNOMIAL_ORD, ///< Ordinary Polynomial 184 // PS_POLYNOMIAL_CHEB ///< Chebyshev Polynomial 185 //} psPolynomialType; 142 186 143 187 /***************************************************************************** 144 188 This routine must allocate memory for the polynomial structures. 145 189 *****************************************************************************/ 146 psPolynomial1D* psPolynomial1DAlloc(int n) 190 psPolynomial1D* psPolynomial1DAlloc(int n, 191 psPolynomialType type) 147 192 { 148 193 int i = 0; … … 151 196 newPoly = (psPolynomial1D* ) psAlloc(sizeof(psPolynomial1D)); 152 197 p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial1DFree); 198 199 newPoly->type = type; 153 200 newPoly->n = n; 154 201 newPoly->coeff = (float *)psAlloc(n * sizeof(float)); … … 161 208 } 162 209 163 return (newPoly); 164 } 165 166 psPolynomial2D* psPolynomial2DAlloc(int nX, int nY) 210 return(newPoly); 211 } 212 213 psPolynomial2D* psPolynomial2DAlloc(int nX, int nY, 214 psPolynomialType type) 167 215 { 168 216 int x = 0; … … 172 220 newPoly = (psPolynomial2D* ) psAlloc(sizeof(psPolynomial2D)); 173 221 p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial2DFree); 222 223 newPoly->type = type; 174 224 newPoly->nX = nX; 175 225 newPoly->nY = nY; … … 191 241 } 192 242 193 return (newPoly); 194 } 195 196 psPolynomial3D* psPolynomial3DAlloc(int nX, int nY, int nZ) 243 return(newPoly); 244 } 245 246 psPolynomial3D* psPolynomial3DAlloc(int nX, int nY, int nZ, 247 psPolynomialType type) 197 248 { 198 249 int x = 0; … … 203 254 newPoly = (psPolynomial3D* ) psAlloc(sizeof(psPolynomial3D)); 204 255 p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial3DFree); 256 257 newPoly->type = type; 205 258 newPoly->nX = nX; 206 259 newPoly->nY = nY; … … 230 283 } 231 284 232 return (newPoly); 233 } 234 235 psPolynomial4D* psPolynomial4DAlloc(int nW, int nX, int nY, int nZ) 285 return(newPoly); 286 } 287 288 psPolynomial4D* psPolynomial4DAlloc(int nW, int nX, int nY, int nZ, 289 psPolynomialType type) 236 290 { 237 291 int w = 0; … … 243 297 newPoly = (psPolynomial4D* ) psAlloc(sizeof(psPolynomial4D)); 244 298 p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial4DFree); 299 300 newPoly->type = type; 245 301 newPoly->nW = nW; 246 302 newPoly->nX = nX; … … 278 334 } 279 335 280 return 336 return(newPoly); 281 337 } 282 338 … … 359 415 XXX: Should the "coeffErr[]" should be used as well? 360 416 *****************************************************************************/ 361 float p sPolynomial1DEval(float x, const psPolynomial1D* myPoly)417 float p_psOrgPolynomial1DEval(float x, const psPolynomial1D* myPoly) 362 418 { 363 419 int loop_x = 0; … … 371 427 // NOTE: Do we want to flag this case? 372 428 if (myPoly->n == 0) { 373 return 429 return(1.0); 374 430 } 375 431 … … 383 439 } 384 440 385 return (polySum); 386 } 387 388 float psPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) 441 return(polySum); 442 } 443 444 // XXX: You can do this without having to psAlloc() vector d. 445 float p_psChebPolynomial1DEval(float x, const psPolynomial1D* myPoly) 446 { 447 psVector *d; 448 int n; 449 int i; 450 float tmp; 451 452 n = myPoly->n; 453 d = psVectorAlloc(n, PS_TYPE_F32); 454 d->data.F32[n-1] = myPoly->coeff[n-1]; 455 d->data.F32[n-2] = (2.0 * x * d->data.F32[n-1]) + myPoly->coeff[n-2]; 456 for (i=n-3;i>=1;i--) { 457 d->data.F32[i] = (2.0 * x * d->data.F32[i+1]) - 458 (d->data.F32[i+2]) + 459 (myPoly->coeff[i]); 460 } 461 462 tmp = (x * d->data.F32[1]) - 463 (d->data.F32[2]) + 464 (0.5 * myPoly->coeff[0]); 465 466 psFree(d); 467 return(tmp); 468 } 469 470 float psPolynomial1DEval(float x, const psPolynomial1D* myPoly) 471 { 472 if (myPoly->type == PS_POLYNOMIAL_ORD) { 473 return(p_psOrgPolynomial1DEval(x, myPoly)); 474 } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { 475 return(p_psChebPolynomial1DEval(x, myPoly)); 476 } else { 477 psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); 478 } 479 return(0.0); 480 } 481 482 psVector *psPolynomial1DEvalVector(const psVector *x, 483 const psPolynomial1D *myPoly) 484 { 485 psVector *tmp; 486 psVector *myX; 487 int i; 488 489 PS_CONVERT_VECTOR_F32(x, myX); 490 491 tmp = psVectorAlloc(x->n, PS_TYPE_F32); 492 for (i=0;i<x->n;i++) { 493 tmp->data.F32[i] = psPolynomial1DEval(x->data.F32[i], myPoly); 494 } 495 496 if (x->type.type != PS_TYPE_F32) { 497 psFree(myX); 498 } 499 return(tmp); 500 } 501 502 503 float p_psOrgPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) 389 504 { 390 505 int loop_x = 0; … … 403 518 } 404 519 405 return (polySum); 406 } 407 408 float psPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) 520 return(polySum); 521 } 522 523 float p_psChebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) 524 { 525 return(0.0); 526 } 527 528 float psPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) 529 { 530 if (myPoly->type == PS_POLYNOMIAL_ORD) { 531 return(p_psOrgPolynomial2DEval(x, y, myPoly)); 532 } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { 533 return(p_psChebPolynomial2DEval(x, y, myPoly)); 534 } else { 535 psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); 536 } 537 return(0.0); 538 } 539 540 541 psVector *psPolynomial2DEvalVector(const psVector *x, 542 const psVector *y, 543 const psPolynomial2D *myPoly) 544 { 545 psVector *tmp; 546 psVector *myX; 547 psVector *myY; 548 int i; 549 int vecLen=x->n; 550 551 PS_CONVERT_VECTOR_F32(x, myX); 552 PS_CONVERT_VECTOR_F32(y, myY); 553 vecLen=x->n; 554 if (y->n < vecLen) { 555 vecLen = y->n; 556 } 557 558 tmp = psVectorAlloc(vecLen, PS_TYPE_F32); 559 for (i=0;i<vecLen;i++) { 560 tmp->data.F32[i] = psPolynomial2DEval(myX->data.F32[i], 561 myY->data.F32[i], 562 myPoly); 563 } 564 565 if (x->type.type != PS_TYPE_F32) { 566 psFree(myX); 567 } 568 if (y->type.type != PS_TYPE_F32) { 569 psFree(myY); 570 } 571 return(tmp); 572 } 573 574 575 576 float p_psOrgPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) 409 577 { 410 578 int loop_x = 0; … … 429 597 } 430 598 431 return (polySum); 432 } 433 434 float psPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) 599 return(polySum); 600 } 601 602 float p_psChebPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) 603 { 604 return(0.0); 605 } 606 607 float psPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) 608 { 609 if (myPoly->type == PS_POLYNOMIAL_ORD) { 610 return(p_psOrgPolynomial3DEval(x, y, z, myPoly)); 611 } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { 612 return(p_psChebPolynomial3DEval(x, y, z, myPoly)); 613 } else { 614 psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); 615 } 616 return(0.0); 617 } 618 619 psVector *psPolynomial3DEvalVector(const psVector *x, 620 const psVector *y, 621 const psVector *z, 622 const psPolynomial3D *myPoly) 623 { 624 psVector *tmp; 625 psVector *myX; 626 psVector *myY; 627 psVector *myZ; 628 int i; 629 int vecLen=x->n; 630 631 PS_CONVERT_VECTOR_F32(x, myX); 632 PS_CONVERT_VECTOR_F32(y, myY); 633 PS_CONVERT_VECTOR_F32(z, myZ); 634 vecLen=x->n; 635 if (y->n < vecLen) { 636 vecLen = y->n; 637 } 638 if (z->n < vecLen) { 639 vecLen = z->n; 640 } 641 642 tmp = psVectorAlloc(vecLen, PS_TYPE_F32); 643 for (i=0;i<vecLen;i++) { 644 tmp->data.F32[i] = psPolynomial3DEval(myX->data.F32[i], 645 myY->data.F32[i], 646 myZ->data.F32[i], 647 myPoly); 648 } 649 650 if (x->type.type != PS_TYPE_F32) { 651 psFree(myX); 652 } 653 if (y->type.type != PS_TYPE_F32) { 654 psFree(myY); 655 } 656 if (z->type.type != PS_TYPE_F32) { 657 psFree(myZ); 658 } 659 return(tmp); 660 } 661 662 663 664 665 666 667 float p_psOrgPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) 435 668 { 436 669 int loop_w = 0; … … 461 694 } 462 695 463 return (polySum); 464 } 465 466 psDPolynomial1D* psDPolynomial1DAlloc(int n) 696 return(polySum); 697 } 698 699 float p_psChebPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) 700 { 701 return(0.0); 702 } 703 704 float psPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) 705 { 706 if (myPoly->type == PS_POLYNOMIAL_ORD) { 707 return(p_psOrgPolynomial4DEval(w,x,y,z, myPoly)); 708 } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { 709 return(p_psChebPolynomial4DEval(w,x,y,z, myPoly)); 710 } else { 711 psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); 712 } 713 return(0.0); 714 } 715 716 psVector *psPolynomial4DEvalVector(const psVector *w, 717 const psVector *x, 718 const psVector *y, 719 const psVector *z, 720 const psPolynomial4D *myPoly) 721 { 722 psVector *tmp; 723 psVector *myW; 724 psVector *myX; 725 psVector *myY; 726 psVector *myZ; 727 int i; 728 int vecLen=x->n; 729 730 PS_CONVERT_VECTOR_F32(w, myW); 731 PS_CONVERT_VECTOR_F32(x, myX); 732 PS_CONVERT_VECTOR_F32(y, myY); 733 PS_CONVERT_VECTOR_F32(z, myZ); 734 vecLen=w->n; 735 if (y->n < vecLen) { 736 vecLen = y->n; 737 } 738 if (x->n < vecLen) { 739 vecLen = x->n; 740 } 741 if (z->n < vecLen) { 742 vecLen = z->n; 743 } 744 745 tmp = psVectorAlloc(vecLen, PS_TYPE_F32); 746 for (i=0;i<vecLen;i++) { 747 tmp->data.F32[i] = psPolynomial4DEval(myW->data.F32[i], 748 myX->data.F32[i], 749 myY->data.F32[i], 750 myZ->data.F32[i], 751 myPoly); 752 } 753 754 if (w->type.type != PS_TYPE_F32) { 755 psFree(myW); 756 } 757 if (x->type.type != PS_TYPE_F32) { 758 psFree(myX); 759 } 760 if (y->type.type != PS_TYPE_F32) { 761 psFree(myY); 762 } 763 if (z->type.type != PS_TYPE_F32) { 764 psFree(myZ); 765 } 766 return(tmp); 767 } 768 769 770 771 772 773 psDPolynomial1D* psDPolynomial1DAlloc(int n, 774 psPolynomialType type) 467 775 { 468 776 int i = 0; … … 471 779 newPoly = (psDPolynomial1D* ) psAlloc(sizeof(psDPolynomial1D)); 472 780 p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial1DFree); 781 782 newPoly->type = type; 473 783 newPoly->n = n; 474 784 newPoly->coeff = (double *)psAlloc(n * sizeof(double)); … … 481 791 } 482 792 483 return (newPoly); 484 } 485 486 psDPolynomial2D* psDPolynomial2DAlloc(int nX, int nY) 793 return(newPoly); 794 } 795 796 psDPolynomial2D* psDPolynomial2DAlloc(int nX, int nY, 797 psPolynomialType type) 487 798 { 488 799 int x = 0; … … 492 803 newPoly = (psDPolynomial2D* ) psAlloc(sizeof(psDPolynomial2D)); 493 804 p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial2DFree); 805 806 newPoly->type = type; 494 807 newPoly->nX = nX; 495 808 newPoly->nY = nY; … … 511 824 } 512 825 513 return (newPoly); 514 } 515 516 psDPolynomial3D* psDPolynomial3DAlloc(int nX, int nY, int nZ) 826 return(newPoly); 827 } 828 829 psDPolynomial3D* psDPolynomial3DAlloc(int nX, int nY, int nZ, 830 psPolynomialType type) 517 831 { 518 832 int x = 0; … … 523 837 newPoly = (psDPolynomial3D* ) psAlloc(sizeof(psDPolynomial3D)); 524 838 p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial3DFree); 839 840 newPoly->type = type; 525 841 newPoly->nX = nX; 526 842 newPoly->nY = nY; … … 550 866 } 551 867 552 return (newPoly); 553 } 554 555 psDPolynomial4D* psDPolynomial4DAlloc(int nW, int nX, int nY, int nZ) 868 return(newPoly); 869 } 870 871 psDPolynomial4D* psDPolynomial4DAlloc(int nW, int nX, int nY, int nZ, 872 psPolynomialType type) 556 873 { 557 874 int w = 0; … … 563 880 newPoly = (psDPolynomial4D* ) psAlloc(sizeof(psDPolynomial4D)); 564 881 p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial4DFree); 882 883 newPoly->type = type; 565 884 newPoly->nW = nW; 566 885 newPoly->nX = nX; … … 598 917 } 599 918 600 return 919 return(newPoly); 601 920 } 602 921 … … 673 992 Polynomial coefficients will be accessed in [w][x][y][z] fashion. 674 993 *****************************************************************************/ 675 double p sDPolynomial1DEval(double x, const psDPolynomial1D* myPoly)994 double p_psDOrgPolynomial1DEval(double x, const psDPolynomial1D* myPoly) 676 995 { 677 996 int loop_x = 0; … … 681 1000 // NOTE: Do we want to flag this case? 682 1001 if (myPoly->n == 0) { 683 return 1002 return(1.0); 684 1003 } 685 1004 … … 689 1008 } 690 1009 691 return (polySum); 692 } 693 694 double psDPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) 1010 return(polySum); 1011 } 1012 1013 double p_psDChebPolynomial1DEval(double x, const psDPolynomial1D* myPoly) 1014 { 1015 return(0.0); 1016 } 1017 1018 double psDPolynomial1DEval(double x, const psDPolynomial1D* myPoly) 1019 { 1020 if (myPoly->type == PS_POLYNOMIAL_ORD) { 1021 return(p_psDOrgPolynomial1DEval(x, myPoly)); 1022 } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { 1023 return(p_psDChebPolynomial1DEval(x, myPoly)); 1024 } else { 1025 psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); 1026 } 1027 return(0.0); 1028 } 1029 1030 psVector *psDPolynomial1DEvalVector(const psVector *x, 1031 const psDPolynomial1D *myPoly) 1032 { 1033 psVector *tmp; 1034 psVector *myX; 1035 int i; 1036 1037 PS_CONVERT_VECTOR_F64(x, myX); 1038 1039 tmp = psVectorAlloc(x->n, PS_TYPE_F64); 1040 for (i=0;i<x->n;i++) { 1041 tmp->data.F64[i] = psDPolynomial1DEval(x->data.F64[i], myPoly); 1042 } 1043 1044 if (x->type.type != PS_TYPE_F64) { 1045 psFree(myX); 1046 } 1047 return(tmp); 1048 } 1049 1050 1051 1052 double p_psDOrgPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) 695 1053 { 696 1054 int loop_x = 0; … … 709 1067 } 710 1068 711 return (polySum); 712 } 713 714 double psDPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) 1069 return(polySum); 1070 } 1071 1072 double p_psDChebPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) 1073 { 1074 return(0.0); 1075 } 1076 1077 double psDPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) 1078 { 1079 if (myPoly->type == PS_POLYNOMIAL_ORD) { 1080 return(p_psDOrgPolynomial2DEval(x, y, myPoly)); 1081 } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { 1082 return(p_psDChebPolynomial2DEval(x, y, myPoly)); 1083 } else { 1084 psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); 1085 } 1086 return(0.0); 1087 } 1088 1089 psVector *psDPolynomial2DEvalVector(const psVector *x, 1090 const psVector *y, 1091 const psDPolynomial2D *myPoly) 1092 { 1093 psVector *tmp; 1094 psVector *myX; 1095 psVector *myY; 1096 int i; 1097 int vecLen=x->n; 1098 1099 PS_CONVERT_VECTOR_F64(x, myX); 1100 PS_CONVERT_VECTOR_F64(y, myY); 1101 vecLen=x->n; 1102 if (y->n < vecLen) { 1103 vecLen = y->n; 1104 } 1105 1106 tmp = psVectorAlloc(vecLen, PS_TYPE_F64); 1107 for (i=0;i<vecLen;i++) { 1108 tmp->data.F64[i] = psDPolynomial2DEval(myX->data.F64[i], 1109 myY->data.F64[i], 1110 myPoly); 1111 } 1112 1113 if (x->type.type != PS_TYPE_F64) { 1114 psFree(myX); 1115 } 1116 if (y->type.type != PS_TYPE_F64) { 1117 psFree(myY); 1118 } 1119 return(tmp); 1120 } 1121 1122 1123 1124 double p_psDOrgPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) 715 1125 { 716 1126 int loop_x = 0; … … 735 1145 } 736 1146 737 return (polySum); 738 } 739 740 double psDPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) 1147 return(polySum); 1148 } 1149 1150 double p_psDChebPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) 1151 { 1152 return(0.0); 1153 } 1154 1155 double psDPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) 1156 { 1157 if (myPoly->type == PS_POLYNOMIAL_ORD) { 1158 return(p_psDOrgPolynomial3DEval(x, y, z, myPoly)); 1159 } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { 1160 return(p_psDChebPolynomial3DEval(x, y, z, myPoly)); 1161 } else { 1162 psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); 1163 } 1164 return(0.0); 1165 } 1166 1167 psVector *psDPolynomial3DEvalVector(const psVector *x, 1168 const psVector *y, 1169 const psVector *z, 1170 const psDPolynomial3D *myPoly) 1171 { 1172 psVector *tmp; 1173 psVector *myX; 1174 psVector *myY; 1175 psVector *myZ; 1176 int i; 1177 int vecLen=x->n; 1178 1179 PS_CONVERT_VECTOR_F64(x, myX); 1180 PS_CONVERT_VECTOR_F64(y, myY); 1181 PS_CONVERT_VECTOR_F64(z, myZ); 1182 vecLen=x->n; 1183 if (y->n < vecLen) { 1184 vecLen = y->n; 1185 } 1186 if (z->n < vecLen) { 1187 vecLen = z->n; 1188 } 1189 1190 tmp = psVectorAlloc(vecLen, PS_TYPE_F64); 1191 for (i=0;i<vecLen;i++) { 1192 tmp->data.F64[i] = psDPolynomial3DEval(myX->data.F64[i], 1193 myY->data.F64[i], 1194 myZ->data.F64[i], 1195 myPoly); 1196 } 1197 1198 if (x->type.type != PS_TYPE_F64) { 1199 psFree(myX); 1200 } 1201 if (y->type.type != PS_TYPE_F64) { 1202 psFree(myY); 1203 } 1204 if (z->type.type != PS_TYPE_F64) { 1205 psFree(myZ); 1206 } 1207 return(tmp); 1208 } 1209 1210 1211 1212 1213 1214 1215 1216 1217 double p_psDOrgPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) 741 1218 { 742 1219 int loop_w = 0; … … 767 1244 } 768 1245 769 return (polySum); 770 } 1246 return(polySum); 1247 } 1248 1249 double p_psDChebPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) 1250 { 1251 return(0.0); 1252 } 1253 1254 double psDPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) 1255 { 1256 if (myPoly->type == PS_POLYNOMIAL_ORD) { 1257 return(p_psDOrgPolynomial4DEval(w,x,y,z, myPoly)); 1258 } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { 1259 return(p_psDChebPolynomial4DEval(w,x,y,z, myPoly)); 1260 } else { 1261 psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); 1262 } 1263 return(0.0); 1264 } 1265 1266 psVector *psDPolynomial4DEvalVector(const psVector *w, 1267 const psVector *x, 1268 const psVector *y, 1269 const psVector *z, 1270 const psDPolynomial4D *myPoly) 1271 { 1272 psVector *tmp; 1273 psVector *myW; 1274 psVector *myX; 1275 psVector *myY; 1276 psVector *myZ; 1277 int i; 1278 int vecLen=x->n; 1279 1280 PS_CONVERT_VECTOR_F64(w, myW); 1281 PS_CONVERT_VECTOR_F64(x, myX); 1282 PS_CONVERT_VECTOR_F64(y, myY); 1283 PS_CONVERT_VECTOR_F64(z, myZ); 1284 vecLen=w->n; 1285 if (y->n < vecLen) { 1286 vecLen = y->n; 1287 } 1288 if (x->n < vecLen) { 1289 vecLen = x->n; 1290 } 1291 if (z->n < vecLen) { 1292 vecLen = z->n; 1293 } 1294 1295 tmp = psVectorAlloc(vecLen, PS_TYPE_F64); 1296 for (i=0;i<vecLen;i++) { 1297 tmp->data.F64[i] = psDPolynomial4DEval(myW->data.F64[i], 1298 myX->data.F64[i], 1299 myY->data.F64[i], 1300 myZ->data.F64[i], 1301 myPoly); 1302 } 1303 1304 if (w->type.type != PS_TYPE_F64) { 1305 psFree(myW); 1306 } 1307 if (x->type.type != PS_TYPE_F64) { 1308 psFree(myX); 1309 } 1310 if (y->type.type != PS_TYPE_F64) { 1311 psFree(myY); 1312 } 1313 if (z->type.type != PS_TYPE_F64) { 1314 psFree(myZ); 1315 } 1316 return(tmp); 1317 } 1318 1319 1320 1321 1322 //typedef struct { 1323 // int n; 1324 // psPolynomial1D **spline; 1325 // float *domains; 1326 //} psSpline1D; 1327 1328 /***************************************************************************** 1329 1330 NOTE: "n" specifies the number of spline polynomials. Therefore, there 1331 must exist n+1 points in "domains". 1332 *****************************************************************************/ 1333 psSpline1D *psSpline1DAlloc(int numSplines, 1334 int order, 1335 float min, 1336 float max) 1337 { 1338 psSpline1D *tmp = NULL; 1339 int i; 1340 float tmpDomain; 1341 float width; 1342 1343 tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D)); 1344 1345 tmp->n = numSplines; 1346 1347 tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); 1348 for (i=0;i<numSplines;i++) { 1349 (tmp->spline)[i] = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD); 1350 } 1351 1352 tmp->domains = (float *) psAlloc((numSplines+1) * sizeof(float)); 1353 width = (max - min) / ((float) numSplines); 1354 tmpDomain = min; 1355 for (i=0;i<numSplines+1;i++) { 1356 (tmp->domains)[i] = tmpDomain; 1357 tmpDomain+= width; 1358 } 1359 1360 return(tmp); 1361 } 1362 1363 psSpline1D *psSpline1DAllocGeneric(const psVector *bounds, 1364 int order) 1365 { 1366 psSpline1D *tmp = NULL; 1367 int i; 1368 int numSplines; 1369 1370 tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D)); 1371 1372 numSplines = bounds->n - 1; 1373 1374 tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); 1375 for (i=0;i<numSplines;i++) { 1376 (tmp->spline)[i] = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD); 1377 } 1378 1379 tmp->domains = (float *) psAlloc((bounds->n) * sizeof(float)); 1380 for (i=0;i<bounds->n;i++) { 1381 (tmp->domains)[i] = bounds->data.F32[i]; 1382 } 1383 1384 return(tmp); 1385 } 1386 1387 /***************************************************************************** 1388 VectorBinDisect(): This is a private function which takes as input a vector 1389 of floating point data as well as a single floating point values. The input 1390 vector values are assumed to be non-decreasing (v[i-1] <= v[j] for all j>=i). 1391 This routine does a binary disection of the vector and returns "i" such 1392 that (v[i] <= x <= v[i+1). If x lies outside the range of v[], 1393 then this routine prints a warning message and returns -1. 1394 *****************************************************************************/ 1395 int VectorBinDisect(float *bins, 1396 int numBins, 1397 float x) 1398 { 1399 int min; 1400 int max; 1401 int mid; 1402 1403 if ((x < bins[0]) || 1404 (x > bins[numBins-1])) { 1405 psLogMsg(__func__, PS_LOG_WARN, 1406 "VectorBinDisect(): ordinate %f is outside vector range (%f - %f).", 1407 x, bins[0], bins[numBins-1]); 1408 return(-1); 1409 } 1410 1411 min = 0; 1412 max = numBins-1; 1413 mid = (max-min)/2; 1414 1415 while (min != max) { 1416 mid = (max-min)/2; 1417 1418 if (x < bins[mid]) { 1419 max = mid; 1420 } else { 1421 min = mid; 1422 } 1423 } 1424 1425 return(min); 1426 } 1427 1428 float psSpline1DEval(const psSpline1D *spline, 1429 float x) 1430 { 1431 int binNum; 1432 int n; 1433 1434 n = spline->n; 1435 binNum = VectorBinDisect(spline->domains, (spline->n)+1, x); 1436 if (binNum == -1) { 1437 psLogMsg(__func__, PS_LOG_WARN, 1438 "psSpline1DEval(): x ordinate (%f) is outside the spline range (%f - %f).", 1439 x, (spline->domains)[0], 1440 (spline->domains)[n-1]); 1441 1442 if (x < (spline->domains)[0]) { 1443 return(psPolynomial1DEval(x, spline->spline[0])); 1444 } else if (x > (spline->domains)[n-1]) { 1445 return(psPolynomial1DEval(x, spline->spline[n-1])); 1446 } 1447 } 1448 1449 return(psPolynomial1DEval(x, spline->spline[binNum])); 1450 } 1451 1452 psVector *psSpline1DEvalVector(const psVector *x, 1453 const psSpline1D *spline) 1454 { 1455 int i; 1456 psVector *tmpVector; 1457 1458 tmpVector = psVectorAlloc(x->n, PS_TYPE_F32); 1459 for (i=0;i<x->n;i++) { 1460 tmpVector->data.F32[i] = psSpline1DEval(spline, x->data.F32[i]); 1461 } 1462 1463 return(tmpVector); 1464 } 1465 -
trunk/psLib/src/math/psPolynomial.c
r1476 r1718 1 2 1 /** @file psFunctions.c 3 2 * … … 8 7 * polynomials. It also contains a Gaussian functions. 9 8 * 10 * @version $Revision: 1.2 3$ $Name: not supported by cvs2svn $11 * @date $Date: 2004-0 8-11 20:07:45$9 * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ 10 * @date $Date: 2004-09-08 06:00:49 $ 12 11 * 13 12 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 32 31 #include "psAbort.h" 33 32 #include "psFunctions.h" 33 #include "psLogMsg.h" 34 34 35 35 #include <gsl/gsl_rng.h> … … 42 42 /*****************************************************************************/ 43 43 44 // None 44 #define PS_CONVERT_VECTOR_F32(IN, OUT) \ 45 if (IN->type.type != PS_TYPE_F32) { \ 46 psLogMsg(__func__, PS_LOG_WARN, \ 47 "Input vector has incorrect type (%d).\n", \ 48 IN->type); \ 49 OUT = psVectorAlloc(IN->n, PS_TYPE_F32); \ 50 \ 51 if (IN->type.type == PS_TYPE_F64) { \ 52 for (int i=0;i<IN->n;i++) { \ 53 OUT->data.F32[i] = (float) (IN->data.F64[i]); \ 54 } \ 55 }\ 56 else { \ 57 psAbort(__func__, "Wrong type.\n"); \ 58 } \ 59 } else { \ 60 OUT = (psVector *) IN; \ 61 } 62 63 #define PS_CONVERT_VECTOR_F64(IN, OUT) \ 64 if (IN->type.type != PS_TYPE_F64) { \ 65 psLogMsg(__func__, PS_LOG_WARN, \ 66 "Input vector has incorrect type (%d).\n", \ 67 IN->type); \ 68 OUT = psVectorAlloc(IN->n, PS_TYPE_F64); \ 69 \ 70 if (IN->type.type == PS_TYPE_F32) { \ 71 for (int i=0;i<IN->n;i++) { \ 72 OUT->data.F64[i] = (float) (IN->data.F32[i]); \ 73 } \ 74 }\ 75 else { \ 76 psAbort(__func__, "Wrong type.\n"); \ 77 } \ 78 } else { \ 79 OUT = (psVector *) IN; \ 80 } 81 82 45 83 46 84 /*****************************************************************************/ … … 107 145 } 108 146 109 return 147 return(tmp * exp(-((x - mean) * (x - mean)) / (2.0 * sigma * sigma))); 110 148 } 111 149 … … 138 176 139 177 // NOTE: Should I free r as well? 140 return (gauss); 141 } 178 return(gauss); 179 } 180 181 //XXX: remove this 182 //typedef enum { 183 // PS_POLYNOMIAL_ORD, ///< Ordinary Polynomial 184 // PS_POLYNOMIAL_CHEB ///< Chebyshev Polynomial 185 //} psPolynomialType; 142 186 143 187 /***************************************************************************** 144 188 This routine must allocate memory for the polynomial structures. 145 189 *****************************************************************************/ 146 psPolynomial1D* psPolynomial1DAlloc(int n) 190 psPolynomial1D* psPolynomial1DAlloc(int n, 191 psPolynomialType type) 147 192 { 148 193 int i = 0; … … 151 196 newPoly = (psPolynomial1D* ) psAlloc(sizeof(psPolynomial1D)); 152 197 p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial1DFree); 198 199 newPoly->type = type; 153 200 newPoly->n = n; 154 201 newPoly->coeff = (float *)psAlloc(n * sizeof(float)); … … 161 208 } 162 209 163 return (newPoly); 164 } 165 166 psPolynomial2D* psPolynomial2DAlloc(int nX, int nY) 210 return(newPoly); 211 } 212 213 psPolynomial2D* psPolynomial2DAlloc(int nX, int nY, 214 psPolynomialType type) 167 215 { 168 216 int x = 0; … … 172 220 newPoly = (psPolynomial2D* ) psAlloc(sizeof(psPolynomial2D)); 173 221 p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial2DFree); 222 223 newPoly->type = type; 174 224 newPoly->nX = nX; 175 225 newPoly->nY = nY; … … 191 241 } 192 242 193 return (newPoly); 194 } 195 196 psPolynomial3D* psPolynomial3DAlloc(int nX, int nY, int nZ) 243 return(newPoly); 244 } 245 246 psPolynomial3D* psPolynomial3DAlloc(int nX, int nY, int nZ, 247 psPolynomialType type) 197 248 { 198 249 int x = 0; … … 203 254 newPoly = (psPolynomial3D* ) psAlloc(sizeof(psPolynomial3D)); 204 255 p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial3DFree); 256 257 newPoly->type = type; 205 258 newPoly->nX = nX; 206 259 newPoly->nY = nY; … … 230 283 } 231 284 232 return (newPoly); 233 } 234 235 psPolynomial4D* psPolynomial4DAlloc(int nW, int nX, int nY, int nZ) 285 return(newPoly); 286 } 287 288 psPolynomial4D* psPolynomial4DAlloc(int nW, int nX, int nY, int nZ, 289 psPolynomialType type) 236 290 { 237 291 int w = 0; … … 243 297 newPoly = (psPolynomial4D* ) psAlloc(sizeof(psPolynomial4D)); 244 298 p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial4DFree); 299 300 newPoly->type = type; 245 301 newPoly->nW = nW; 246 302 newPoly->nX = nX; … … 278 334 } 279 335 280 return 336 return(newPoly); 281 337 } 282 338 … … 359 415 XXX: Should the "coeffErr[]" should be used as well? 360 416 *****************************************************************************/ 361 float p sPolynomial1DEval(float x, const psPolynomial1D* myPoly)417 float p_psOrgPolynomial1DEval(float x, const psPolynomial1D* myPoly) 362 418 { 363 419 int loop_x = 0; … … 371 427 // NOTE: Do we want to flag this case? 372 428 if (myPoly->n == 0) { 373 return 429 return(1.0); 374 430 } 375 431 … … 383 439 } 384 440 385 return (polySum); 386 } 387 388 float psPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) 441 return(polySum); 442 } 443 444 // XXX: You can do this without having to psAlloc() vector d. 445 float p_psChebPolynomial1DEval(float x, const psPolynomial1D* myPoly) 446 { 447 psVector *d; 448 int n; 449 int i; 450 float tmp; 451 452 n = myPoly->n; 453 d = psVectorAlloc(n, PS_TYPE_F32); 454 d->data.F32[n-1] = myPoly->coeff[n-1]; 455 d->data.F32[n-2] = (2.0 * x * d->data.F32[n-1]) + myPoly->coeff[n-2]; 456 for (i=n-3;i>=1;i--) { 457 d->data.F32[i] = (2.0 * x * d->data.F32[i+1]) - 458 (d->data.F32[i+2]) + 459 (myPoly->coeff[i]); 460 } 461 462 tmp = (x * d->data.F32[1]) - 463 (d->data.F32[2]) + 464 (0.5 * myPoly->coeff[0]); 465 466 psFree(d); 467 return(tmp); 468 } 469 470 float psPolynomial1DEval(float x, const psPolynomial1D* myPoly) 471 { 472 if (myPoly->type == PS_POLYNOMIAL_ORD) { 473 return(p_psOrgPolynomial1DEval(x, myPoly)); 474 } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { 475 return(p_psChebPolynomial1DEval(x, myPoly)); 476 } else { 477 psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); 478 } 479 return(0.0); 480 } 481 482 psVector *psPolynomial1DEvalVector(const psVector *x, 483 const psPolynomial1D *myPoly) 484 { 485 psVector *tmp; 486 psVector *myX; 487 int i; 488 489 PS_CONVERT_VECTOR_F32(x, myX); 490 491 tmp = psVectorAlloc(x->n, PS_TYPE_F32); 492 for (i=0;i<x->n;i++) { 493 tmp->data.F32[i] = psPolynomial1DEval(x->data.F32[i], myPoly); 494 } 495 496 if (x->type.type != PS_TYPE_F32) { 497 psFree(myX); 498 } 499 return(tmp); 500 } 501 502 503 float p_psOrgPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) 389 504 { 390 505 int loop_x = 0; … … 403 518 } 404 519 405 return (polySum); 406 } 407 408 float psPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) 520 return(polySum); 521 } 522 523 float p_psChebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) 524 { 525 return(0.0); 526 } 527 528 float psPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) 529 { 530 if (myPoly->type == PS_POLYNOMIAL_ORD) { 531 return(p_psOrgPolynomial2DEval(x, y, myPoly)); 532 } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { 533 return(p_psChebPolynomial2DEval(x, y, myPoly)); 534 } else { 535 psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); 536 } 537 return(0.0); 538 } 539 540 541 psVector *psPolynomial2DEvalVector(const psVector *x, 542 const psVector *y, 543 const psPolynomial2D *myPoly) 544 { 545 psVector *tmp; 546 psVector *myX; 547 psVector *myY; 548 int i; 549 int vecLen=x->n; 550 551 PS_CONVERT_VECTOR_F32(x, myX); 552 PS_CONVERT_VECTOR_F32(y, myY); 553 vecLen=x->n; 554 if (y->n < vecLen) { 555 vecLen = y->n; 556 } 557 558 tmp = psVectorAlloc(vecLen, PS_TYPE_F32); 559 for (i=0;i<vecLen;i++) { 560 tmp->data.F32[i] = psPolynomial2DEval(myX->data.F32[i], 561 myY->data.F32[i], 562 myPoly); 563 } 564 565 if (x->type.type != PS_TYPE_F32) { 566 psFree(myX); 567 } 568 if (y->type.type != PS_TYPE_F32) { 569 psFree(myY); 570 } 571 return(tmp); 572 } 573 574 575 576 float p_psOrgPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) 409 577 { 410 578 int loop_x = 0; … … 429 597 } 430 598 431 return (polySum); 432 } 433 434 float psPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) 599 return(polySum); 600 } 601 602 float p_psChebPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) 603 { 604 return(0.0); 605 } 606 607 float psPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) 608 { 609 if (myPoly->type == PS_POLYNOMIAL_ORD) { 610 return(p_psOrgPolynomial3DEval(x, y, z, myPoly)); 611 } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { 612 return(p_psChebPolynomial3DEval(x, y, z, myPoly)); 613 } else { 614 psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); 615 } 616 return(0.0); 617 } 618 619 psVector *psPolynomial3DEvalVector(const psVector *x, 620 const psVector *y, 621 const psVector *z, 622 const psPolynomial3D *myPoly) 623 { 624 psVector *tmp; 625 psVector *myX; 626 psVector *myY; 627 psVector *myZ; 628 int i; 629 int vecLen=x->n; 630 631 PS_CONVERT_VECTOR_F32(x, myX); 632 PS_CONVERT_VECTOR_F32(y, myY); 633 PS_CONVERT_VECTOR_F32(z, myZ); 634 vecLen=x->n; 635 if (y->n < vecLen) { 636 vecLen = y->n; 637 } 638 if (z->n < vecLen) { 639 vecLen = z->n; 640 } 641 642 tmp = psVectorAlloc(vecLen, PS_TYPE_F32); 643 for (i=0;i<vecLen;i++) { 644 tmp->data.F32[i] = psPolynomial3DEval(myX->data.F32[i], 645 myY->data.F32[i], 646 myZ->data.F32[i], 647 myPoly); 648 } 649 650 if (x->type.type != PS_TYPE_F32) { 651 psFree(myX); 652 } 653 if (y->type.type != PS_TYPE_F32) { 654 psFree(myY); 655 } 656 if (z->type.type != PS_TYPE_F32) { 657 psFree(myZ); 658 } 659 return(tmp); 660 } 661 662 663 664 665 666 667 float p_psOrgPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) 435 668 { 436 669 int loop_w = 0; … … 461 694 } 462 695 463 return (polySum); 464 } 465 466 psDPolynomial1D* psDPolynomial1DAlloc(int n) 696 return(polySum); 697 } 698 699 float p_psChebPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) 700 { 701 return(0.0); 702 } 703 704 float psPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) 705 { 706 if (myPoly->type == PS_POLYNOMIAL_ORD) { 707 return(p_psOrgPolynomial4DEval(w,x,y,z, myPoly)); 708 } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { 709 return(p_psChebPolynomial4DEval(w,x,y,z, myPoly)); 710 } else { 711 psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); 712 } 713 return(0.0); 714 } 715 716 psVector *psPolynomial4DEvalVector(const psVector *w, 717 const psVector *x, 718 const psVector *y, 719 const psVector *z, 720 const psPolynomial4D *myPoly) 721 { 722 psVector *tmp; 723 psVector *myW; 724 psVector *myX; 725 psVector *myY; 726 psVector *myZ; 727 int i; 728 int vecLen=x->n; 729 730 PS_CONVERT_VECTOR_F32(w, myW); 731 PS_CONVERT_VECTOR_F32(x, myX); 732 PS_CONVERT_VECTOR_F32(y, myY); 733 PS_CONVERT_VECTOR_F32(z, myZ); 734 vecLen=w->n; 735 if (y->n < vecLen) { 736 vecLen = y->n; 737 } 738 if (x->n < vecLen) { 739 vecLen = x->n; 740 } 741 if (z->n < vecLen) { 742 vecLen = z->n; 743 } 744 745 tmp = psVectorAlloc(vecLen, PS_TYPE_F32); 746 for (i=0;i<vecLen;i++) { 747 tmp->data.F32[i] = psPolynomial4DEval(myW->data.F32[i], 748 myX->data.F32[i], 749 myY->data.F32[i], 750 myZ->data.F32[i], 751 myPoly); 752 } 753 754 if (w->type.type != PS_TYPE_F32) { 755 psFree(myW); 756 } 757 if (x->type.type != PS_TYPE_F32) { 758 psFree(myX); 759 } 760 if (y->type.type != PS_TYPE_F32) { 761 psFree(myY); 762 } 763 if (z->type.type != PS_TYPE_F32) { 764 psFree(myZ); 765 } 766 return(tmp); 767 } 768 769 770 771 772 773 psDPolynomial1D* psDPolynomial1DAlloc(int n, 774 psPolynomialType type) 467 775 { 468 776 int i = 0; … … 471 779 newPoly = (psDPolynomial1D* ) psAlloc(sizeof(psDPolynomial1D)); 472 780 p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial1DFree); 781 782 newPoly->type = type; 473 783 newPoly->n = n; 474 784 newPoly->coeff = (double *)psAlloc(n * sizeof(double)); … … 481 791 } 482 792 483 return (newPoly); 484 } 485 486 psDPolynomial2D* psDPolynomial2DAlloc(int nX, int nY) 793 return(newPoly); 794 } 795 796 psDPolynomial2D* psDPolynomial2DAlloc(int nX, int nY, 797 psPolynomialType type) 487 798 { 488 799 int x = 0; … … 492 803 newPoly = (psDPolynomial2D* ) psAlloc(sizeof(psDPolynomial2D)); 493 804 p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial2DFree); 805 806 newPoly->type = type; 494 807 newPoly->nX = nX; 495 808 newPoly->nY = nY; … … 511 824 } 512 825 513 return (newPoly); 514 } 515 516 psDPolynomial3D* psDPolynomial3DAlloc(int nX, int nY, int nZ) 826 return(newPoly); 827 } 828 829 psDPolynomial3D* psDPolynomial3DAlloc(int nX, int nY, int nZ, 830 psPolynomialType type) 517 831 { 518 832 int x = 0; … … 523 837 newPoly = (psDPolynomial3D* ) psAlloc(sizeof(psDPolynomial3D)); 524 838 p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial3DFree); 839 840 newPoly->type = type; 525 841 newPoly->nX = nX; 526 842 newPoly->nY = nY; … … 550 866 } 551 867 552 return (newPoly); 553 } 554 555 psDPolynomial4D* psDPolynomial4DAlloc(int nW, int nX, int nY, int nZ) 868 return(newPoly); 869 } 870 871 psDPolynomial4D* psDPolynomial4DAlloc(int nW, int nX, int nY, int nZ, 872 psPolynomialType type) 556 873 { 557 874 int w = 0; … … 563 880 newPoly = (psDPolynomial4D* ) psAlloc(sizeof(psDPolynomial4D)); 564 881 p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial4DFree); 882 883 newPoly->type = type; 565 884 newPoly->nW = nW; 566 885 newPoly->nX = nX; … … 598 917 } 599 918 600 return 919 return(newPoly); 601 920 } 602 921 … … 673 992 Polynomial coefficients will be accessed in [w][x][y][z] fashion. 674 993 *****************************************************************************/ 675 double p sDPolynomial1DEval(double x, const psDPolynomial1D* myPoly)994 double p_psDOrgPolynomial1DEval(double x, const psDPolynomial1D* myPoly) 676 995 { 677 996 int loop_x = 0; … … 681 1000 // NOTE: Do we want to flag this case? 682 1001 if (myPoly->n == 0) { 683 return 1002 return(1.0); 684 1003 } 685 1004 … … 689 1008 } 690 1009 691 return (polySum); 692 } 693 694 double psDPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) 1010 return(polySum); 1011 } 1012 1013 double p_psDChebPolynomial1DEval(double x, const psDPolynomial1D* myPoly) 1014 { 1015 return(0.0); 1016 } 1017 1018 double psDPolynomial1DEval(double x, const psDPolynomial1D* myPoly) 1019 { 1020 if (myPoly->type == PS_POLYNOMIAL_ORD) { 1021 return(p_psDOrgPolynomial1DEval(x, myPoly)); 1022 } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { 1023 return(p_psDChebPolynomial1DEval(x, myPoly)); 1024 } else { 1025 psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); 1026 } 1027 return(0.0); 1028 } 1029 1030 psVector *psDPolynomial1DEvalVector(const psVector *x, 1031 const psDPolynomial1D *myPoly) 1032 { 1033 psVector *tmp; 1034 psVector *myX; 1035 int i; 1036 1037 PS_CONVERT_VECTOR_F64(x, myX); 1038 1039 tmp = psVectorAlloc(x->n, PS_TYPE_F64); 1040 for (i=0;i<x->n;i++) { 1041 tmp->data.F64[i] = psDPolynomial1DEval(x->data.F64[i], myPoly); 1042 } 1043 1044 if (x->type.type != PS_TYPE_F64) { 1045 psFree(myX); 1046 } 1047 return(tmp); 1048 } 1049 1050 1051 1052 double p_psDOrgPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) 695 1053 { 696 1054 int loop_x = 0; … … 709 1067 } 710 1068 711 return (polySum); 712 } 713 714 double psDPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) 1069 return(polySum); 1070 } 1071 1072 double p_psDChebPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) 1073 { 1074 return(0.0); 1075 } 1076 1077 double psDPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) 1078 { 1079 if (myPoly->type == PS_POLYNOMIAL_ORD) { 1080 return(p_psDOrgPolynomial2DEval(x, y, myPoly)); 1081 } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { 1082 return(p_psDChebPolynomial2DEval(x, y, myPoly)); 1083 } else { 1084 psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); 1085 } 1086 return(0.0); 1087 } 1088 1089 psVector *psDPolynomial2DEvalVector(const psVector *x, 1090 const psVector *y, 1091 const psDPolynomial2D *myPoly) 1092 { 1093 psVector *tmp; 1094 psVector *myX; 1095 psVector *myY; 1096 int i; 1097 int vecLen=x->n; 1098 1099 PS_CONVERT_VECTOR_F64(x, myX); 1100 PS_CONVERT_VECTOR_F64(y, myY); 1101 vecLen=x->n; 1102 if (y->n < vecLen) { 1103 vecLen = y->n; 1104 } 1105 1106 tmp = psVectorAlloc(vecLen, PS_TYPE_F64); 1107 for (i=0;i<vecLen;i++) { 1108 tmp->data.F64[i] = psDPolynomial2DEval(myX->data.F64[i], 1109 myY->data.F64[i], 1110 myPoly); 1111 } 1112 1113 if (x->type.type != PS_TYPE_F64) { 1114 psFree(myX); 1115 } 1116 if (y->type.type != PS_TYPE_F64) { 1117 psFree(myY); 1118 } 1119 return(tmp); 1120 } 1121 1122 1123 1124 double p_psDOrgPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) 715 1125 { 716 1126 int loop_x = 0; … … 735 1145 } 736 1146 737 return (polySum); 738 } 739 740 double psDPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) 1147 return(polySum); 1148 } 1149 1150 double p_psDChebPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) 1151 { 1152 return(0.0); 1153 } 1154 1155 double psDPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) 1156 { 1157 if (myPoly->type == PS_POLYNOMIAL_ORD) { 1158 return(p_psDOrgPolynomial3DEval(x, y, z, myPoly)); 1159 } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { 1160 return(p_psDChebPolynomial3DEval(x, y, z, myPoly)); 1161 } else { 1162 psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); 1163 } 1164 return(0.0); 1165 } 1166 1167 psVector *psDPolynomial3DEvalVector(const psVector *x, 1168 const psVector *y, 1169 const psVector *z, 1170 const psDPolynomial3D *myPoly) 1171 { 1172 psVector *tmp; 1173 psVector *myX; 1174 psVector *myY; 1175 psVector *myZ; 1176 int i; 1177 int vecLen=x->n; 1178 1179 PS_CONVERT_VECTOR_F64(x, myX); 1180 PS_CONVERT_VECTOR_F64(y, myY); 1181 PS_CONVERT_VECTOR_F64(z, myZ); 1182 vecLen=x->n; 1183 if (y->n < vecLen) { 1184 vecLen = y->n; 1185 } 1186 if (z->n < vecLen) { 1187 vecLen = z->n; 1188 } 1189 1190 tmp = psVectorAlloc(vecLen, PS_TYPE_F64); 1191 for (i=0;i<vecLen;i++) { 1192 tmp->data.F64[i] = psDPolynomial3DEval(myX->data.F64[i], 1193 myY->data.F64[i], 1194 myZ->data.F64[i], 1195 myPoly); 1196 } 1197 1198 if (x->type.type != PS_TYPE_F64) { 1199 psFree(myX); 1200 } 1201 if (y->type.type != PS_TYPE_F64) { 1202 psFree(myY); 1203 } 1204 if (z->type.type != PS_TYPE_F64) { 1205 psFree(myZ); 1206 } 1207 return(tmp); 1208 } 1209 1210 1211 1212 1213 1214 1215 1216 1217 double p_psDOrgPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) 741 1218 { 742 1219 int loop_w = 0; … … 767 1244 } 768 1245 769 return (polySum); 770 } 1246 return(polySum); 1247 } 1248 1249 double p_psDChebPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) 1250 { 1251 return(0.0); 1252 } 1253 1254 double psDPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) 1255 { 1256 if (myPoly->type == PS_POLYNOMIAL_ORD) { 1257 return(p_psDOrgPolynomial4DEval(w,x,y,z, myPoly)); 1258 } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { 1259 return(p_psDChebPolynomial4DEval(w,x,y,z, myPoly)); 1260 } else { 1261 psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); 1262 } 1263 return(0.0); 1264 } 1265 1266 psVector *psDPolynomial4DEvalVector(const psVector *w, 1267 const psVector *x, 1268 const psVector *y, 1269 const psVector *z, 1270 const psDPolynomial4D *myPoly) 1271 { 1272 psVector *tmp; 1273 psVector *myW; 1274 psVector *myX; 1275 psVector *myY; 1276 psVector *myZ; 1277 int i; 1278 int vecLen=x->n; 1279 1280 PS_CONVERT_VECTOR_F64(w, myW); 1281 PS_CONVERT_VECTOR_F64(x, myX); 1282 PS_CONVERT_VECTOR_F64(y, myY); 1283 PS_CONVERT_VECTOR_F64(z, myZ); 1284 vecLen=w->n; 1285 if (y->n < vecLen) { 1286 vecLen = y->n; 1287 } 1288 if (x->n < vecLen) { 1289 vecLen = x->n; 1290 } 1291 if (z->n < vecLen) { 1292 vecLen = z->n; 1293 } 1294 1295 tmp = psVectorAlloc(vecLen, PS_TYPE_F64); 1296 for (i=0;i<vecLen;i++) { 1297 tmp->data.F64[i] = psDPolynomial4DEval(myW->data.F64[i], 1298 myX->data.F64[i], 1299 myY->data.F64[i], 1300 myZ->data.F64[i], 1301 myPoly); 1302 } 1303 1304 if (w->type.type != PS_TYPE_F64) { 1305 psFree(myW); 1306 } 1307 if (x->type.type != PS_TYPE_F64) { 1308 psFree(myX); 1309 } 1310 if (y->type.type != PS_TYPE_F64) { 1311 psFree(myY); 1312 } 1313 if (z->type.type != PS_TYPE_F64) { 1314 psFree(myZ); 1315 } 1316 return(tmp); 1317 } 1318 1319 1320 1321 1322 //typedef struct { 1323 // int n; 1324 // psPolynomial1D **spline; 1325 // float *domains; 1326 //} psSpline1D; 1327 1328 /***************************************************************************** 1329 1330 NOTE: "n" specifies the number of spline polynomials. Therefore, there 1331 must exist n+1 points in "domains". 1332 *****************************************************************************/ 1333 psSpline1D *psSpline1DAlloc(int numSplines, 1334 int order, 1335 float min, 1336 float max) 1337 { 1338 psSpline1D *tmp = NULL; 1339 int i; 1340 float tmpDomain; 1341 float width; 1342 1343 tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D)); 1344 1345 tmp->n = numSplines; 1346 1347 tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); 1348 for (i=0;i<numSplines;i++) { 1349 (tmp->spline)[i] = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD); 1350 } 1351 1352 tmp->domains = (float *) psAlloc((numSplines+1) * sizeof(float)); 1353 width = (max - min) / ((float) numSplines); 1354 tmpDomain = min; 1355 for (i=0;i<numSplines+1;i++) { 1356 (tmp->domains)[i] = tmpDomain; 1357 tmpDomain+= width; 1358 } 1359 1360 return(tmp); 1361 } 1362 1363 psSpline1D *psSpline1DAllocGeneric(const psVector *bounds, 1364 int order) 1365 { 1366 psSpline1D *tmp = NULL; 1367 int i; 1368 int numSplines; 1369 1370 tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D)); 1371 1372 numSplines = bounds->n - 1; 1373 1374 tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); 1375 for (i=0;i<numSplines;i++) { 1376 (tmp->spline)[i] = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD); 1377 } 1378 1379 tmp->domains = (float *) psAlloc((bounds->n) * sizeof(float)); 1380 for (i=0;i<bounds->n;i++) { 1381 (tmp->domains)[i] = bounds->data.F32[i]; 1382 } 1383 1384 return(tmp); 1385 } 1386 1387 /***************************************************************************** 1388 VectorBinDisect(): This is a private function which takes as input a vector 1389 of floating point data as well as a single floating point values. The input 1390 vector values are assumed to be non-decreasing (v[i-1] <= v[j] for all j>=i). 1391 This routine does a binary disection of the vector and returns "i" such 1392 that (v[i] <= x <= v[i+1). If x lies outside the range of v[], 1393 then this routine prints a warning message and returns -1. 1394 *****************************************************************************/ 1395 int VectorBinDisect(float *bins, 1396 int numBins, 1397 float x) 1398 { 1399 int min; 1400 int max; 1401 int mid; 1402 1403 if ((x < bins[0]) || 1404 (x > bins[numBins-1])) { 1405 psLogMsg(__func__, PS_LOG_WARN, 1406 "VectorBinDisect(): ordinate %f is outside vector range (%f - %f).", 1407 x, bins[0], bins[numBins-1]); 1408 return(-1); 1409 } 1410 1411 min = 0; 1412 max = numBins-1; 1413 mid = (max-min)/2; 1414 1415 while (min != max) { 1416 mid = (max-min)/2; 1417 1418 if (x < bins[mid]) { 1419 max = mid; 1420 } else { 1421 min = mid; 1422 } 1423 } 1424 1425 return(min); 1426 } 1427 1428 float psSpline1DEval(const psSpline1D *spline, 1429 float x) 1430 { 1431 int binNum; 1432 int n; 1433 1434 n = spline->n; 1435 binNum = VectorBinDisect(spline->domains, (spline->n)+1, x); 1436 if (binNum == -1) { 1437 psLogMsg(__func__, PS_LOG_WARN, 1438 "psSpline1DEval(): x ordinate (%f) is outside the spline range (%f - %f).", 1439 x, (spline->domains)[0], 1440 (spline->domains)[n-1]); 1441 1442 if (x < (spline->domains)[0]) { 1443 return(psPolynomial1DEval(x, spline->spline[0])); 1444 } else if (x > (spline->domains)[n-1]) { 1445 return(psPolynomial1DEval(x, spline->spline[n-1])); 1446 } 1447 } 1448 1449 return(psPolynomial1DEval(x, spline->spline[binNum])); 1450 } 1451 1452 psVector *psSpline1DEvalVector(const psVector *x, 1453 const psSpline1D *spline) 1454 { 1455 int i; 1456 psVector *tmpVector; 1457 1458 tmpVector = psVectorAlloc(x->n, PS_TYPE_F32); 1459 for (i=0;i<x->n;i++) { 1460 tmpVector->data.F32[i] = psSpline1DEval(spline, x->data.F32[i]); 1461 } 1462 1463 return(tmpVector); 1464 } 1465 -
trunk/psLib/src/math/psSpline.c
r1476 r1718 1 2 1 /** @file psFunctions.c 3 2 * … … 8 7 * polynomials. It also contains a Gaussian functions. 9 8 * 10 * @version $Revision: 1.2 3$ $Name: not supported by cvs2svn $11 * @date $Date: 2004-0 8-11 20:07:45$9 * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ 10 * @date $Date: 2004-09-08 06:00:49 $ 12 11 * 13 12 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 32 31 #include "psAbort.h" 33 32 #include "psFunctions.h" 33 #include "psLogMsg.h" 34 34 35 35 #include <gsl/gsl_rng.h> … … 42 42 /*****************************************************************************/ 43 43 44 // None 44 #define PS_CONVERT_VECTOR_F32(IN, OUT) \ 45 if (IN->type.type != PS_TYPE_F32) { \ 46 psLogMsg(__func__, PS_LOG_WARN, \ 47 "Input vector has incorrect type (%d).\n", \ 48 IN->type); \ 49 OUT = psVectorAlloc(IN->n, PS_TYPE_F32); \ 50 \ 51 if (IN->type.type == PS_TYPE_F64) { \ 52 for (int i=0;i<IN->n;i++) { \ 53 OUT->data.F32[i] = (float) (IN->data.F64[i]); \ 54 } \ 55 }\ 56 else { \ 57 psAbort(__func__, "Wrong type.\n"); \ 58 } \ 59 } else { \ 60 OUT = (psVector *) IN; \ 61 } 62 63 #define PS_CONVERT_VECTOR_F64(IN, OUT) \ 64 if (IN->type.type != PS_TYPE_F64) { \ 65 psLogMsg(__func__, PS_LOG_WARN, \ 66 "Input vector has incorrect type (%d).\n", \ 67 IN->type); \ 68 OUT = psVectorAlloc(IN->n, PS_TYPE_F64); \ 69 \ 70 if (IN->type.type == PS_TYPE_F32) { \ 71 for (int i=0;i<IN->n;i++) { \ 72 OUT->data.F64[i] = (float) (IN->data.F32[i]); \ 73 } \ 74 }\ 75 else { \ 76 psAbort(__func__, "Wrong type.\n"); \ 77 } \ 78 } else { \ 79 OUT = (psVector *) IN; \ 80 } 81 82 45 83 46 84 /*****************************************************************************/ … … 107 145 } 108 146 109 return 147 return(tmp * exp(-((x - mean) * (x - mean)) / (2.0 * sigma * sigma))); 110 148 } 111 149 … … 138 176 139 177 // NOTE: Should I free r as well? 140 return (gauss); 141 } 178 return(gauss); 179 } 180 181 //XXX: remove this 182 //typedef enum { 183 // PS_POLYNOMIAL_ORD, ///< Ordinary Polynomial 184 // PS_POLYNOMIAL_CHEB ///< Chebyshev Polynomial 185 //} psPolynomialType; 142 186 143 187 /***************************************************************************** 144 188 This routine must allocate memory for the polynomial structures. 145 189 *****************************************************************************/ 146 psPolynomial1D* psPolynomial1DAlloc(int n) 190 psPolynomial1D* psPolynomial1DAlloc(int n, 191 psPolynomialType type) 147 192 { 148 193 int i = 0; … … 151 196 newPoly = (psPolynomial1D* ) psAlloc(sizeof(psPolynomial1D)); 152 197 p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial1DFree); 198 199 newPoly->type = type; 153 200 newPoly->n = n; 154 201 newPoly->coeff = (float *)psAlloc(n * sizeof(float)); … … 161 208 } 162 209 163 return (newPoly); 164 } 165 166 psPolynomial2D* psPolynomial2DAlloc(int nX, int nY) 210 return(newPoly); 211 } 212 213 psPolynomial2D* psPolynomial2DAlloc(int nX, int nY, 214 psPolynomialType type) 167 215 { 168 216 int x = 0; … … 172 220 newPoly = (psPolynomial2D* ) psAlloc(sizeof(psPolynomial2D)); 173 221 p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial2DFree); 222 223 newPoly->type = type; 174 224 newPoly->nX = nX; 175 225 newPoly->nY = nY; … … 191 241 } 192 242 193 return (newPoly); 194 } 195 196 psPolynomial3D* psPolynomial3DAlloc(int nX, int nY, int nZ) 243 return(newPoly); 244 } 245 246 psPolynomial3D* psPolynomial3DAlloc(int nX, int nY, int nZ, 247 psPolynomialType type) 197 248 { 198 249 int x = 0; … … 203 254 newPoly = (psPolynomial3D* ) psAlloc(sizeof(psPolynomial3D)); 204 255 p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial3DFree); 256 257 newPoly->type = type; 205 258 newPoly->nX = nX; 206 259 newPoly->nY = nY; … … 230 283 } 231 284 232 return (newPoly); 233 } 234 235 psPolynomial4D* psPolynomial4DAlloc(int nW, int nX, int nY, int nZ) 285 return(newPoly); 286 } 287 288 psPolynomial4D* psPolynomial4DAlloc(int nW, int nX, int nY, int nZ, 289 psPolynomialType type) 236 290 { 237 291 int w = 0; … … 243 297 newPoly = (psPolynomial4D* ) psAlloc(sizeof(psPolynomial4D)); 244 298 p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial4DFree); 299 300 newPoly->type = type; 245 301 newPoly->nW = nW; 246 302 newPoly->nX = nX; … … 278 334 } 279 335 280 return 336 return(newPoly); 281 337 } 282 338 … … 359 415 XXX: Should the "coeffErr[]" should be used as well? 360 416 *****************************************************************************/ 361 float p sPolynomial1DEval(float x, const psPolynomial1D* myPoly)417 float p_psOrgPolynomial1DEval(float x, const psPolynomial1D* myPoly) 362 418 { 363 419 int loop_x = 0; … … 371 427 // NOTE: Do we want to flag this case? 372 428 if (myPoly->n == 0) { 373 return 429 return(1.0); 374 430 } 375 431 … … 383 439 } 384 440 385 return (polySum); 386 } 387 388 float psPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) 441 return(polySum); 442 } 443 444 // XXX: You can do this without having to psAlloc() vector d. 445 float p_psChebPolynomial1DEval(float x, const psPolynomial1D* myPoly) 446 { 447 psVector *d; 448 int n; 449 int i; 450 float tmp; 451 452 n = myPoly->n; 453 d = psVectorAlloc(n, PS_TYPE_F32); 454 d->data.F32[n-1] = myPoly->coeff[n-1]; 455 d->data.F32[n-2] = (2.0 * x * d->data.F32[n-1]) + myPoly->coeff[n-2]; 456 for (i=n-3;i>=1;i--) { 457 d->data.F32[i] = (2.0 * x * d->data.F32[i+1]) - 458 (d->data.F32[i+2]) + 459 (myPoly->coeff[i]); 460 } 461 462 tmp = (x * d->data.F32[1]) - 463 (d->data.F32[2]) + 464 (0.5 * myPoly->coeff[0]); 465 466 psFree(d); 467 return(tmp); 468 } 469 470 float psPolynomial1DEval(float x, const psPolynomial1D* myPoly) 471 { 472 if (myPoly->type == PS_POLYNOMIAL_ORD) { 473 return(p_psOrgPolynomial1DEval(x, myPoly)); 474 } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { 475 return(p_psChebPolynomial1DEval(x, myPoly)); 476 } else { 477 psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); 478 } 479 return(0.0); 480 } 481 482 psVector *psPolynomial1DEvalVector(const psVector *x, 483 const psPolynomial1D *myPoly) 484 { 485 psVector *tmp; 486 psVector *myX; 487 int i; 488 489 PS_CONVERT_VECTOR_F32(x, myX); 490 491 tmp = psVectorAlloc(x->n, PS_TYPE_F32); 492 for (i=0;i<x->n;i++) { 493 tmp->data.F32[i] = psPolynomial1DEval(x->data.F32[i], myPoly); 494 } 495 496 if (x->type.type != PS_TYPE_F32) { 497 psFree(myX); 498 } 499 return(tmp); 500 } 501 502 503 float p_psOrgPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) 389 504 { 390 505 int loop_x = 0; … … 403 518 } 404 519 405 return (polySum); 406 } 407 408 float psPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) 520 return(polySum); 521 } 522 523 float p_psChebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) 524 { 525 return(0.0); 526 } 527 528 float psPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) 529 { 530 if (myPoly->type == PS_POLYNOMIAL_ORD) { 531 return(p_psOrgPolynomial2DEval(x, y, myPoly)); 532 } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { 533 return(p_psChebPolynomial2DEval(x, y, myPoly)); 534 } else { 535 psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); 536 } 537 return(0.0); 538 } 539 540 541 psVector *psPolynomial2DEvalVector(const psVector *x, 542 const psVector *y, 543 const psPolynomial2D *myPoly) 544 { 545 psVector *tmp; 546 psVector *myX; 547 psVector *myY; 548 int i; 549 int vecLen=x->n; 550 551 PS_CONVERT_VECTOR_F32(x, myX); 552 PS_CONVERT_VECTOR_F32(y, myY); 553 vecLen=x->n; 554 if (y->n < vecLen) { 555 vecLen = y->n; 556 } 557 558 tmp = psVectorAlloc(vecLen, PS_TYPE_F32); 559 for (i=0;i<vecLen;i++) { 560 tmp->data.F32[i] = psPolynomial2DEval(myX->data.F32[i], 561 myY->data.F32[i], 562 myPoly); 563 } 564 565 if (x->type.type != PS_TYPE_F32) { 566 psFree(myX); 567 } 568 if (y->type.type != PS_TYPE_F32) { 569 psFree(myY); 570 } 571 return(tmp); 572 } 573 574 575 576 float p_psOrgPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) 409 577 { 410 578 int loop_x = 0; … … 429 597 } 430 598 431 return (polySum); 432 } 433 434 float psPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) 599 return(polySum); 600 } 601 602 float p_psChebPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) 603 { 604 return(0.0); 605 } 606 607 float psPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) 608 { 609 if (myPoly->type == PS_POLYNOMIAL_ORD) { 610 return(p_psOrgPolynomial3DEval(x, y, z, myPoly)); 611 } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { 612 return(p_psChebPolynomial3DEval(x, y, z, myPoly)); 613 } else { 614 psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); 615 } 616 return(0.0); 617 } 618 619 psVector *psPolynomial3DEvalVector(const psVector *x, 620 const psVector *y, 621 const psVector *z, 622 const psPolynomial3D *myPoly) 623 { 624 psVector *tmp; 625 psVector *myX; 626 psVector *myY; 627 psVector *myZ; 628 int i; 629 int vecLen=x->n; 630 631 PS_CONVERT_VECTOR_F32(x, myX); 632 PS_CONVERT_VECTOR_F32(y, myY); 633 PS_CONVERT_VECTOR_F32(z, myZ); 634 vecLen=x->n; 635 if (y->n < vecLen) { 636 vecLen = y->n; 637 } 638 if (z->n < vecLen) { 639 vecLen = z->n; 640 } 641 642 tmp = psVectorAlloc(vecLen, PS_TYPE_F32); 643 for (i=0;i<vecLen;i++) { 644 tmp->data.F32[i] = psPolynomial3DEval(myX->data.F32[i], 645 myY->data.F32[i], 646 myZ->data.F32[i], 647 myPoly); 648 } 649 650 if (x->type.type != PS_TYPE_F32) { 651 psFree(myX); 652 } 653 if (y->type.type != PS_TYPE_F32) { 654 psFree(myY); 655 } 656 if (z->type.type != PS_TYPE_F32) { 657 psFree(myZ); 658 } 659 return(tmp); 660 } 661 662 663 664 665 666 667 float p_psOrgPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) 435 668 { 436 669 int loop_w = 0; … … 461 694 } 462 695 463 return (polySum); 464 } 465 466 psDPolynomial1D* psDPolynomial1DAlloc(int n) 696 return(polySum); 697 } 698 699 float p_psChebPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) 700 { 701 return(0.0); 702 } 703 704 float psPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) 705 { 706 if (myPoly->type == PS_POLYNOMIAL_ORD) { 707 return(p_psOrgPolynomial4DEval(w,x,y,z, myPoly)); 708 } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { 709 return(p_psChebPolynomial4DEval(w,x,y,z, myPoly)); 710 } else { 711 psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); 712 } 713 return(0.0); 714 } 715 716 psVector *psPolynomial4DEvalVector(const psVector *w, 717 const psVector *x, 718 const psVector *y, 719 const psVector *z, 720 const psPolynomial4D *myPoly) 721 { 722 psVector *tmp; 723 psVector *myW; 724 psVector *myX; 725 psVector *myY; 726 psVector *myZ; 727 int i; 728 int vecLen=x->n; 729 730 PS_CONVERT_VECTOR_F32(w, myW); 731 PS_CONVERT_VECTOR_F32(x, myX); 732 PS_CONVERT_VECTOR_F32(y, myY); 733 PS_CONVERT_VECTOR_F32(z, myZ); 734 vecLen=w->n; 735 if (y->n < vecLen) { 736 vecLen = y->n; 737 } 738 if (x->n < vecLen) { 739 vecLen = x->n; 740 } 741 if (z->n < vecLen) { 742 vecLen = z->n; 743 } 744 745 tmp = psVectorAlloc(vecLen, PS_TYPE_F32); 746 for (i=0;i<vecLen;i++) { 747 tmp->data.F32[i] = psPolynomial4DEval(myW->data.F32[i], 748 myX->data.F32[i], 749 myY->data.F32[i], 750 myZ->data.F32[i], 751 myPoly); 752 } 753 754 if (w->type.type != PS_TYPE_F32) { 755 psFree(myW); 756 } 757 if (x->type.type != PS_TYPE_F32) { 758 psFree(myX); 759 } 760 if (y->type.type != PS_TYPE_F32) { 761 psFree(myY); 762 } 763 if (z->type.type != PS_TYPE_F32) { 764 psFree(myZ); 765 } 766 return(tmp); 767 } 768 769 770 771 772 773 psDPolynomial1D* psDPolynomial1DAlloc(int n, 774 psPolynomialType type) 467 775 { 468 776 int i = 0; … … 471 779 newPoly = (psDPolynomial1D* ) psAlloc(sizeof(psDPolynomial1D)); 472 780 p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial1DFree); 781 782 newPoly->type = type; 473 783 newPoly->n = n; 474 784 newPoly->coeff = (double *)psAlloc(n * sizeof(double)); … … 481 791 } 482 792 483 return (newPoly); 484 } 485 486 psDPolynomial2D* psDPolynomial2DAlloc(int nX, int nY) 793 return(newPoly); 794 } 795 796 psDPolynomial2D* psDPolynomial2DAlloc(int nX, int nY, 797 psPolynomialType type) 487 798 { 488 799 int x = 0; … … 492 803 newPoly = (psDPolynomial2D* ) psAlloc(sizeof(psDPolynomial2D)); 493 804 p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial2DFree); 805 806 newPoly->type = type; 494 807 newPoly->nX = nX; 495 808 newPoly->nY = nY; … … 511 824 } 512 825 513 return (newPoly); 514 } 515 516 psDPolynomial3D* psDPolynomial3DAlloc(int nX, int nY, int nZ) 826 return(newPoly); 827 } 828 829 psDPolynomial3D* psDPolynomial3DAlloc(int nX, int nY, int nZ, 830 psPolynomialType type) 517 831 { 518 832 int x = 0; … … 523 837 newPoly = (psDPolynomial3D* ) psAlloc(sizeof(psDPolynomial3D)); 524 838 p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial3DFree); 839 840 newPoly->type = type; 525 841 newPoly->nX = nX; 526 842 newPoly->nY = nY; … … 550 866 } 551 867 552 return (newPoly); 553 } 554 555 psDPolynomial4D* psDPolynomial4DAlloc(int nW, int nX, int nY, int nZ) 868 return(newPoly); 869 } 870 871 psDPolynomial4D* psDPolynomial4DAlloc(int nW, int nX, int nY, int nZ, 872 psPolynomialType type) 556 873 { 557 874 int w = 0; … … 563 880 newPoly = (psDPolynomial4D* ) psAlloc(sizeof(psDPolynomial4D)); 564 881 p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial4DFree); 882 883 newPoly->type = type; 565 884 newPoly->nW = nW; 566 885 newPoly->nX = nX; … … 598 917 } 599 918 600 return 919 return(newPoly); 601 920 } 602 921 … … 673 992 Polynomial coefficients will be accessed in [w][x][y][z] fashion. 674 993 *****************************************************************************/ 675 double p sDPolynomial1DEval(double x, const psDPolynomial1D* myPoly)994 double p_psDOrgPolynomial1DEval(double x, const psDPolynomial1D* myPoly) 676 995 { 677 996 int loop_x = 0; … … 681 1000 // NOTE: Do we want to flag this case? 682 1001 if (myPoly->n == 0) { 683 return 1002 return(1.0); 684 1003 } 685 1004 … … 689 1008 } 690 1009 691 return (polySum); 692 } 693 694 double psDPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) 1010 return(polySum); 1011 } 1012 1013 double p_psDChebPolynomial1DEval(double x, const psDPolynomial1D* myPoly) 1014 { 1015 return(0.0); 1016 } 1017 1018 double psDPolynomial1DEval(double x, const psDPolynomial1D* myPoly) 1019 { 1020 if (myPoly->type == PS_POLYNOMIAL_ORD) { 1021 return(p_psDOrgPolynomial1DEval(x, myPoly)); 1022 } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { 1023 return(p_psDChebPolynomial1DEval(x, myPoly)); 1024 } else { 1025 psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); 1026 } 1027 return(0.0); 1028 } 1029 1030 psVector *psDPolynomial1DEvalVector(const psVector *x, 1031 const psDPolynomial1D *myPoly) 1032 { 1033 psVector *tmp; 1034 psVector *myX; 1035 int i; 1036 1037 PS_CONVERT_VECTOR_F64(x, myX); 1038 1039 tmp = psVectorAlloc(x->n, PS_TYPE_F64); 1040 for (i=0;i<x->n;i++) { 1041 tmp->data.F64[i] = psDPolynomial1DEval(x->data.F64[i], myPoly); 1042 } 1043 1044 if (x->type.type != PS_TYPE_F64) { 1045 psFree(myX); 1046 } 1047 return(tmp); 1048 } 1049 1050 1051 1052 double p_psDOrgPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) 695 1053 { 696 1054 int loop_x = 0; … … 709 1067 } 710 1068 711 return (polySum); 712 } 713 714 double psDPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) 1069 return(polySum); 1070 } 1071 1072 double p_psDChebPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) 1073 { 1074 return(0.0); 1075 } 1076 1077 double psDPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) 1078 { 1079 if (myPoly->type == PS_POLYNOMIAL_ORD) { 1080 return(p_psDOrgPolynomial2DEval(x, y, myPoly)); 1081 } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { 1082 return(p_psDChebPolynomial2DEval(x, y, myPoly)); 1083 } else { 1084 psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); 1085 } 1086 return(0.0); 1087 } 1088 1089 psVector *psDPolynomial2DEvalVector(const psVector *x, 1090 const psVector *y, 1091 const psDPolynomial2D *myPoly) 1092 { 1093 psVector *tmp; 1094 psVector *myX; 1095 psVector *myY; 1096 int i; 1097 int vecLen=x->n; 1098 1099 PS_CONVERT_VECTOR_F64(x, myX); 1100 PS_CONVERT_VECTOR_F64(y, myY); 1101 vecLen=x->n; 1102 if (y->n < vecLen) { 1103 vecLen = y->n; 1104 } 1105 1106 tmp = psVectorAlloc(vecLen, PS_TYPE_F64); 1107 for (i=0;i<vecLen;i++) { 1108 tmp->data.F64[i] = psDPolynomial2DEval(myX->data.F64[i], 1109 myY->data.F64[i], 1110 myPoly); 1111 } 1112 1113 if (x->type.type != PS_TYPE_F64) { 1114 psFree(myX); 1115 } 1116 if (y->type.type != PS_TYPE_F64) { 1117 psFree(myY); 1118 } 1119 return(tmp); 1120 } 1121 1122 1123 1124 double p_psDOrgPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) 715 1125 { 716 1126 int loop_x = 0; … … 735 1145 } 736 1146 737 return (polySum); 738 } 739 740 double psDPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) 1147 return(polySum); 1148 } 1149 1150 double p_psDChebPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) 1151 { 1152 return(0.0); 1153 } 1154 1155 double psDPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) 1156 { 1157 if (myPoly->type == PS_POLYNOMIAL_ORD) { 1158 return(p_psDOrgPolynomial3DEval(x, y, z, myPoly)); 1159 } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { 1160 return(p_psDChebPolynomial3DEval(x, y, z, myPoly)); 1161 } else { 1162 psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); 1163 } 1164 return(0.0); 1165 } 1166 1167 psVector *psDPolynomial3DEvalVector(const psVector *x, 1168 const psVector *y, 1169 const psVector *z, 1170 const psDPolynomial3D *myPoly) 1171 { 1172 psVector *tmp; 1173 psVector *myX; 1174 psVector *myY; 1175 psVector *myZ; 1176 int i; 1177 int vecLen=x->n; 1178 1179 PS_CONVERT_VECTOR_F64(x, myX); 1180 PS_CONVERT_VECTOR_F64(y, myY); 1181 PS_CONVERT_VECTOR_F64(z, myZ); 1182 vecLen=x->n; 1183 if (y->n < vecLen) { 1184 vecLen = y->n; 1185 } 1186 if (z->n < vecLen) { 1187 vecLen = z->n; 1188 } 1189 1190 tmp = psVectorAlloc(vecLen, PS_TYPE_F64); 1191 for (i=0;i<vecLen;i++) { 1192 tmp->data.F64[i] = psDPolynomial3DEval(myX->data.F64[i], 1193 myY->data.F64[i], 1194 myZ->data.F64[i], 1195 myPoly); 1196 } 1197 1198 if (x->type.type != PS_TYPE_F64) { 1199 psFree(myX); 1200 } 1201 if (y->type.type != PS_TYPE_F64) { 1202 psFree(myY); 1203 } 1204 if (z->type.type != PS_TYPE_F64) { 1205 psFree(myZ); 1206 } 1207 return(tmp); 1208 } 1209 1210 1211 1212 1213 1214 1215 1216 1217 double p_psDOrgPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) 741 1218 { 742 1219 int loop_w = 0; … … 767 1244 } 768 1245 769 return (polySum); 770 } 1246 return(polySum); 1247 } 1248 1249 double p_psDChebPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) 1250 { 1251 return(0.0); 1252 } 1253 1254 double psDPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) 1255 { 1256 if (myPoly->type == PS_POLYNOMIAL_ORD) { 1257 return(p_psDOrgPolynomial4DEval(w,x,y,z, myPoly)); 1258 } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { 1259 return(p_psDChebPolynomial4DEval(w,x,y,z, myPoly)); 1260 } else { 1261 psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); 1262 } 1263 return(0.0); 1264 } 1265 1266 psVector *psDPolynomial4DEvalVector(const psVector *w, 1267 const psVector *x, 1268 const psVector *y, 1269 const psVector *z, 1270 const psDPolynomial4D *myPoly) 1271 { 1272 psVector *tmp; 1273 psVector *myW; 1274 psVector *myX; 1275 psVector *myY; 1276 psVector *myZ; 1277 int i; 1278 int vecLen=x->n; 1279 1280 PS_CONVERT_VECTOR_F64(w, myW); 1281 PS_CONVERT_VECTOR_F64(x, myX); 1282 PS_CONVERT_VECTOR_F64(y, myY); 1283 PS_CONVERT_VECTOR_F64(z, myZ); 1284 vecLen=w->n; 1285 if (y->n < vecLen) { 1286 vecLen = y->n; 1287 } 1288 if (x->n < vecLen) { 1289 vecLen = x->n; 1290 } 1291 if (z->n < vecLen) { 1292 vecLen = z->n; 1293 } 1294 1295 tmp = psVectorAlloc(vecLen, PS_TYPE_F64); 1296 for (i=0;i<vecLen;i++) { 1297 tmp->data.F64[i] = psDPolynomial4DEval(myW->data.F64[i], 1298 myX->data.F64[i], 1299 myY->data.F64[i], 1300 myZ->data.F64[i], 1301 myPoly); 1302 } 1303 1304 if (w->type.type != PS_TYPE_F64) { 1305 psFree(myW); 1306 } 1307 if (x->type.type != PS_TYPE_F64) { 1308 psFree(myX); 1309 } 1310 if (y->type.type != PS_TYPE_F64) { 1311 psFree(myY); 1312 } 1313 if (z->type.type != PS_TYPE_F64) { 1314 psFree(myZ); 1315 } 1316 return(tmp); 1317 } 1318 1319 1320 1321 1322 //typedef struct { 1323 // int n; 1324 // psPolynomial1D **spline; 1325 // float *domains; 1326 //} psSpline1D; 1327 1328 /***************************************************************************** 1329 1330 NOTE: "n" specifies the number of spline polynomials. Therefore, there 1331 must exist n+1 points in "domains". 1332 *****************************************************************************/ 1333 psSpline1D *psSpline1DAlloc(int numSplines, 1334 int order, 1335 float min, 1336 float max) 1337 { 1338 psSpline1D *tmp = NULL; 1339 int i; 1340 float tmpDomain; 1341 float width; 1342 1343 tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D)); 1344 1345 tmp->n = numSplines; 1346 1347 tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); 1348 for (i=0;i<numSplines;i++) { 1349 (tmp->spline)[i] = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD); 1350 } 1351 1352 tmp->domains = (float *) psAlloc((numSplines+1) * sizeof(float)); 1353 width = (max - min) / ((float) numSplines); 1354 tmpDomain = min; 1355 for (i=0;i<numSplines+1;i++) { 1356 (tmp->domains)[i] = tmpDomain; 1357 tmpDomain+= width; 1358 } 1359 1360 return(tmp); 1361 } 1362 1363 psSpline1D *psSpline1DAllocGeneric(const psVector *bounds, 1364 int order) 1365 { 1366 psSpline1D *tmp = NULL; 1367 int i; 1368 int numSplines; 1369 1370 tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D)); 1371 1372 numSplines = bounds->n - 1; 1373 1374 tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); 1375 for (i=0;i<numSplines;i++) { 1376 (tmp->spline)[i] = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD); 1377 } 1378 1379 tmp->domains = (float *) psAlloc((bounds->n) * sizeof(float)); 1380 for (i=0;i<bounds->n;i++) { 1381 (tmp->domains)[i] = bounds->data.F32[i]; 1382 } 1383 1384 return(tmp); 1385 } 1386 1387 /***************************************************************************** 1388 VectorBinDisect(): This is a private function which takes as input a vector 1389 of floating point data as well as a single floating point values. The input 1390 vector values are assumed to be non-decreasing (v[i-1] <= v[j] for all j>=i). 1391 This routine does a binary disection of the vector and returns "i" such 1392 that (v[i] <= x <= v[i+1). If x lies outside the range of v[], 1393 then this routine prints a warning message and returns -1. 1394 *****************************************************************************/ 1395 int VectorBinDisect(float *bins, 1396 int numBins, 1397 float x) 1398 { 1399 int min; 1400 int max; 1401 int mid; 1402 1403 if ((x < bins[0]) || 1404 (x > bins[numBins-1])) { 1405 psLogMsg(__func__, PS_LOG_WARN, 1406 "VectorBinDisect(): ordinate %f is outside vector range (%f - %f).", 1407 x, bins[0], bins[numBins-1]); 1408 return(-1); 1409 } 1410 1411 min = 0; 1412 max = numBins-1; 1413 mid = (max-min)/2; 1414 1415 while (min != max) { 1416 mid = (max-min)/2; 1417 1418 if (x < bins[mid]) { 1419 max = mid; 1420 } else { 1421 min = mid; 1422 } 1423 } 1424 1425 return(min); 1426 } 1427 1428 float psSpline1DEval(const psSpline1D *spline, 1429 float x) 1430 { 1431 int binNum; 1432 int n; 1433 1434 n = spline->n; 1435 binNum = VectorBinDisect(spline->domains, (spline->n)+1, x); 1436 if (binNum == -1) { 1437 psLogMsg(__func__, PS_LOG_WARN, 1438 "psSpline1DEval(): x ordinate (%f) is outside the spline range (%f - %f).", 1439 x, (spline->domains)[0], 1440 (spline->domains)[n-1]); 1441 1442 if (x < (spline->domains)[0]) { 1443 return(psPolynomial1DEval(x, spline->spline[0])); 1444 } else if (x > (spline->domains)[n-1]) { 1445 return(psPolynomial1DEval(x, spline->spline[n-1])); 1446 } 1447 } 1448 1449 return(psPolynomial1DEval(x, spline->spline[binNum])); 1450 } 1451 1452 psVector *psSpline1DEvalVector(const psVector *x, 1453 const psSpline1D *spline) 1454 { 1455 int i; 1456 psVector *tmpVector; 1457 1458 tmpVector = psVectorAlloc(x->n, PS_TYPE_F32); 1459 for (i=0;i<x->n;i++) { 1460 tmpVector->data.F32[i] = psSpline1DEval(spline, x->data.F32[i]); 1461 } 1462 1463 return(tmpVector); 1464 } 1465 -
trunk/psLib/src/sys/psTrace.c
r1713 r1718 9 9 * @author George Gusciora, MHPCC 10 10 * 11 * @version $Revision: 1.2 1$ $Name: not supported by cvs2svn $12 * @date $Date: 2004-09-08 0 0:09:06$11 * @version $Revision: 1.22 $ $Name: not supported by cvs2svn $ 12 * @date $Date: 2004-09-08 05:59:41 $ 13 13 * 14 14 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 27 27 not the same thing as a node's "level", which corresponds to the 28 28 trace level of that node. 29 30 I think the following is the correct behavior, but not sure: 31 PS_UNKNOWN_TRACE_LEVEL: We never set the level of a component to this 32 value. This value is only used when psTraceGetLevel is called with 33 a bad component name, or if the component root is undefined, I think. 34 35 PS_DEFAULT_TRACE_LEVEL: This should only be used when adding the 36 intermediate components of a long name. Ie. the "B" in .A.B.C 29 37 30 38 *****************************************************************************/ … … 40 48 #include "psString.h" 41 49 #include "psError.h" 42 #include "psLogMsg.h" 43 44 #include "psSysUtilsErrors.h" 45 #define ERRORNAME_PREFIX PS_ERRORNAME_DOMAIN ".psTrace." 50 #include "psAbort.h" 46 51 47 52 static p_psComponent* cRoot = NULL; // The root of the trace component … … 98 103 /***************************************************************************** 99 104 Set all trace levels to zero. 105 106 XXX: Currently, no function calls this routine. 100 107 *****************************************************************************/ 101 108 void p_psTraceReset(p_psComponent* currentNode) … … 110 117 for (i = 0; i < currentNode->n; i++) { 111 118 if (NULL == currentNode->subcomp[i]) { 112 psLogMsg(ERRORNAME_PREFIX "p_psTraceReset", PS_LOG_WARN, 113 PS_ERRORTEXT_psTrace_NULL_SUBCOMPONENT, 114 i, currentNode->name); 119 psError(__func__, 120 "Sub-component %d of node %s in the trace tree is NULL.\n", i, currentNode->name); 115 121 } else { 116 122 p_psTraceReset(currentNode->subcomp[i]); … … 125 131 void psTraceReset() 126 132 { 127 psFree(cRoot);133 componentFree(cRoot); 128 134 cRoot = NULL; 129 135 } … … 135 141 to ANSI-C. 136 142 *****************************************************************************/ 137 static boolcomponentAdd(const char *addNodeName, int level)143 static void componentAdd(const char *addNodeName, int level) 138 144 { 139 145 int i = 0; // Loop index variable. … … 146 152 // XXX: Verify that this is the correct behavior. 147 153 if (strcmp("", addNodeName) == 0) { 148 psErrorMsg(ERRORNAME_PREFIX "componentAdd",PS_ERR_BAD_PARAMETER_NULL,true, 149 PS_ERRORTEXT_psTrace_ADD_NULL_COMPONENT); 150 return false; 154 psAbort(__func__, "Failed to add null component to trace tree.\n"); 151 155 } 152 156 … … 154 158 if (strcmp(".", addNodeName) == 0) { 155 159 cRoot->level = level; 156 return true;160 return; 157 161 } 158 162 159 163 if (addNodeName[0] != '.') { 160 psErrorMsg(ERRORNAME_PREFIX "componentAdd",PS_ERR_BAD_PARAMETER_VALUE,true, 161 PS_ERRORTEXT_psTrace_MALFORMED_COMPONENT_NAME, 162 addNodeName); 163 return false; 164 printf("ERROR: failed to add %s to the root component tree.\n", addNodeName); 165 exit(1); 164 166 } 165 167 … … 186 188 currentNode->subcomp = psRealloc(currentNode->subcomp, 187 189 (currentNode->n + 1) * sizeof(p_psComponent* )); 188 currentNode->n ++;190 currentNode->n = (currentNode->n) + 1; 189 191 190 192 if (pname == NULL) { 191 currentNode->subcomp[currentNode->n - 1] = componentAlloc(firstComponent, level); 193 // This is the final component to add. 194 currentNode->subcomp[(currentNode->n) - 1] = componentAlloc(firstComponent, level); 192 195 } else { 193 currentNode->subcomp[currentNode->n - 1] = componentAlloc(firstComponent, currentNode->level); 194 } 195 currentNode = currentNode->subcomp[currentNode->n - 1]; 196 } 197 } 198 199 return true; 196 // We are adding an intermediate component. The trace level 197 // is not defined. An undefined trace level inherits the 198 // trace level of it's parent. However, we do not set that 199 // specifically here since that would inheritance to be a 200 // static, one-time, type of behavior. 201 202 currentNode->subcomp[(currentNode->n) - 1] = componentAlloc(firstComponent, PS_DEFAULT_TRACE_LEVEL); 203 } 204 currentNode = currentNode->subcomp[(currentNode->n) - 1]; 205 } 206 } 200 207 } 201 208 … … 211 218 zero 212 219 *****************************************************************************/ 213 boolpsTraceSetLevel(const char *comp, // component of interest214 220 int psTraceSetLevel(const char *comp, // component of interest 221 int level) // desired trace level 215 222 { 216 223 // If the root component tree does not exist, then initialize it. … … 219 226 } 220 227 // Add the new component to the component tree. 221 if ( !componentAdd(comp, level) ) { 222 psErrorMsg(ERRORNAME_PREFIX "psTraceSetLevel", PS_ERR_UNKNOWN, false, 223 PS_ERRORTEXT_psTrace_FAILED_TO_ADD_COMPONENT,level,comp); 224 return false; 225 } 228 componentAdd(comp, level); 229 226 230 // return 0 on success. 227 return true;231 return 0; 228 232 } 229 233 … … 249 253 p_psComponent* currentNode = cRoot; 250 254 int i = 0; 255 int defaultLevel = 0; 251 256 252 257 if (NULL == currentNode) { … … 262 267 } 263 268 269 defaultLevel = cRoot->level; 264 270 strcpy(name, aname); 265 271 pname = &name[1]; … … 268 274 for (i = 0; i < currentNode->n; i++) { 269 275 if (NULL == currentNode->subcomp[i]) { 270 psLogMsg(ERRORNAME_PREFIX "p_psTraceReset", PS_LOG_WARN, 271 PS_ERRORTEXT_psTrace_NULL_SUBCOMPONENT, 272 i, currentNode->name); 276 psError(__func__, 277 "Sub-component %d of node %s in trace tree is NULL.\n", i, currentNode->name); 273 278 } 274 279 275 280 if (strcmp(currentNode->subcomp[i]->name, firstComponent) == 0) { 276 281 currentNode = currentNode->subcomp[i]; 282 // For level inheritance purpose, we save the level of this 283 // component if it is not DEFAULT. 284 if (currentNode->level != PS_DEFAULT_TRACE_LEVEL) { 285 defaultLevel = currentNode->level; 286 } 287 // Determine if thisis the last component: 277 288 if (pname == NULL) { 278 return (currentNode->level); 289 if (currentNode->level != PS_DEFAULT_TRACE_LEVEL) { 290 return (currentNode->level); 291 } else { 292 return(PS_DEFAULT_TRACE_LEVEL); 293 } 279 294 } 280 295 } 281 296 } 282 297 } 283 return (PS_UNKNOWN_TRACE_LEVEL);284 } 285 286 /***************************************************************************** 287 ps GetTraceLevel()298 return(defaultLevel); 299 } 300 301 /***************************************************************************** 302 psTraceLevelGet() 288 303 Return a trace level of "name" in the root component tree. If the 289 304 exact string of components in "name" does not exist in the root … … 303 318 // Search the component root tree, determine the trace level. 304 319 return (doGetTraceLevel(name)); 320 } 321 322 /***************************************************************************** 323 doPrintTraceLevelsOld() 324 This function recursively searches the component tree supplied by the 325 parameter "comp" and prints the name and level of each component. 326 Inputs: 327 comp: a node in the component tree. 328 level: the level of that node 329 Outputs: 330 none 331 Returns: 332 null 333 *****************************************************************************/ 334 static void doPrintTraceLevelsOld(const p_psComponent* comp, int depth) 335 { 336 int i = 0; 337 338 if (comp->name[0] == '\0') { 339 printf("%*s%-*s %d\n", depth, "", 20 - depth, 340 "(root)", (comp->level == PS_DEFAULT_TRACE_LEVEL) ? 0 : comp->level); 341 } else { 342 if (comp->level == PS_DEFAULT_TRACE_LEVEL) { 343 printf("%*s%-*s %s\n", depth, "", 20 - depth, comp->name, "."); 344 } else { 345 printf("%*s%-*s %d\n", depth, "", 20 - depth, comp->name, 346 comp->level); 347 } 348 } 349 350 for (i = 0; i < comp->n; i++) { 351 doPrintTraceLevelsOld(comp->subcomp[i], depth + 1); 352 } 305 353 } 306 354 … … 317 365 null 318 366 *****************************************************************************/ 319 static void doPrintTraceLevels(const p_psComponent* comp, int depth) 367 static void doPrintTraceLevels(const p_psComponent* comp, 368 int depth, 369 int defLevel) 320 370 { 321 371 int i = 0; 322 372 323 373 if (comp->name[0] == '\0') { 324 printf("%*s%-*s %d\n", depth, "", 20 - depth, 325 "(root)", (comp->level == PS_UNKNOWN_TRACE_LEVEL) ? 0 : comp->level); 374 psAbort(__func__, "component name is NULL\n"); 326 375 } else { 327 if (comp->level == PS_UNKNOWN_TRACE_LEVEL) { 328 printf("%*s%-*s %s\n", depth, "", 20 - depth, comp->name, "."); 376 if (comp->level == PS_DEFAULT_TRACE_LEVEL) { 377 printf("%*s%-*s %d\n", depth, "", 20 - depth, comp->name, 378 defLevel); 329 379 } else { 330 printf("%*s%-*s %d\n", depth, "", 20 - depth, comp->name, comp->level); 380 printf("%*s%-*s %d\n", depth, "", 20 - depth, comp->name, 381 comp->level); 331 382 } 332 383 } 333 384 334 385 for (i = 0; i < comp->n; i++) { 335 doPrintTraceLevels(comp->subcomp[i], depth + 1); 336 } 337 } 386 if (comp->level == PS_DEFAULT_TRACE_LEVEL) { 387 doPrintTraceLevels(comp->subcomp[i], depth + 1, defLevel); 388 } else { 389 doPrintTraceLevels(comp->subcomp[i], depth + 1, comp->level); 390 } 391 } 392 } 393 338 394 339 395 /***************************************************************************** … … 353 409 } 354 410 355 doPrintTraceLevels(cRoot, 0 );411 doPrintTraceLevels(cRoot, 0, PS_DEFAULT_TRACE_LEVEL); 356 412 } 357 413 … … 379 435 380 436 if (NULL == comp) { 381 psErrorMsg(ERRORNAME_PREFIX "psTrace", PS_ERR_BAD_PARAMETER_NULL, true, 382 PS_ERRORTEXT_psTrace_NULL_TRACETREE, 383 __func__); 384 return; 437 psError(__func__, "p_psTrace() called on a NULL trace level tree\n"); 385 438 } 386 439 // Only display this message if it's trace level is less than the level -
trunk/psLib/src/sysUtils/psTrace.c
r1713 r1718 9 9 * @author George Gusciora, MHPCC 10 10 * 11 * @version $Revision: 1.2 1$ $Name: not supported by cvs2svn $12 * @date $Date: 2004-09-08 0 0:09:06$11 * @version $Revision: 1.22 $ $Name: not supported by cvs2svn $ 12 * @date $Date: 2004-09-08 05:59:41 $ 13 13 * 14 14 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 27 27 not the same thing as a node's "level", which corresponds to the 28 28 trace level of that node. 29 30 I think the following is the correct behavior, but not sure: 31 PS_UNKNOWN_TRACE_LEVEL: We never set the level of a component to this 32 value. This value is only used when psTraceGetLevel is called with 33 a bad component name, or if the component root is undefined, I think. 34 35 PS_DEFAULT_TRACE_LEVEL: This should only be used when adding the 36 intermediate components of a long name. Ie. the "B" in .A.B.C 29 37 30 38 *****************************************************************************/ … … 40 48 #include "psString.h" 41 49 #include "psError.h" 42 #include "psLogMsg.h" 43 44 #include "psSysUtilsErrors.h" 45 #define ERRORNAME_PREFIX PS_ERRORNAME_DOMAIN ".psTrace." 50 #include "psAbort.h" 46 51 47 52 static p_psComponent* cRoot = NULL; // The root of the trace component … … 98 103 /***************************************************************************** 99 104 Set all trace levels to zero. 105 106 XXX: Currently, no function calls this routine. 100 107 *****************************************************************************/ 101 108 void p_psTraceReset(p_psComponent* currentNode) … … 110 117 for (i = 0; i < currentNode->n; i++) { 111 118 if (NULL == currentNode->subcomp[i]) { 112 psLogMsg(ERRORNAME_PREFIX "p_psTraceReset", PS_LOG_WARN, 113 PS_ERRORTEXT_psTrace_NULL_SUBCOMPONENT, 114 i, currentNode->name); 119 psError(__func__, 120 "Sub-component %d of node %s in the trace tree is NULL.\n", i, currentNode->name); 115 121 } else { 116 122 p_psTraceReset(currentNode->subcomp[i]); … … 125 131 void psTraceReset() 126 132 { 127 psFree(cRoot);133 componentFree(cRoot); 128 134 cRoot = NULL; 129 135 } … … 135 141 to ANSI-C. 136 142 *****************************************************************************/ 137 static boolcomponentAdd(const char *addNodeName, int level)143 static void componentAdd(const char *addNodeName, int level) 138 144 { 139 145 int i = 0; // Loop index variable. … … 146 152 // XXX: Verify that this is the correct behavior. 147 153 if (strcmp("", addNodeName) == 0) { 148 psErrorMsg(ERRORNAME_PREFIX "componentAdd",PS_ERR_BAD_PARAMETER_NULL,true, 149 PS_ERRORTEXT_psTrace_ADD_NULL_COMPONENT); 150 return false; 154 psAbort(__func__, "Failed to add null component to trace tree.\n"); 151 155 } 152 156 … … 154 158 if (strcmp(".", addNodeName) == 0) { 155 159 cRoot->level = level; 156 return true;160 return; 157 161 } 158 162 159 163 if (addNodeName[0] != '.') { 160 psErrorMsg(ERRORNAME_PREFIX "componentAdd",PS_ERR_BAD_PARAMETER_VALUE,true, 161 PS_ERRORTEXT_psTrace_MALFORMED_COMPONENT_NAME, 162 addNodeName); 163 return false; 164 printf("ERROR: failed to add %s to the root component tree.\n", addNodeName); 165 exit(1); 164 166 } 165 167 … … 186 188 currentNode->subcomp = psRealloc(currentNode->subcomp, 187 189 (currentNode->n + 1) * sizeof(p_psComponent* )); 188 currentNode->n ++;190 currentNode->n = (currentNode->n) + 1; 189 191 190 192 if (pname == NULL) { 191 currentNode->subcomp[currentNode->n - 1] = componentAlloc(firstComponent, level); 193 // This is the final component to add. 194 currentNode->subcomp[(currentNode->n) - 1] = componentAlloc(firstComponent, level); 192 195 } else { 193 currentNode->subcomp[currentNode->n - 1] = componentAlloc(firstComponent, currentNode->level); 194 } 195 currentNode = currentNode->subcomp[currentNode->n - 1]; 196 } 197 } 198 199 return true; 196 // We are adding an intermediate component. The trace level 197 // is not defined. An undefined trace level inherits the 198 // trace level of it's parent. However, we do not set that 199 // specifically here since that would inheritance to be a 200 // static, one-time, type of behavior. 201 202 currentNode->subcomp[(currentNode->n) - 1] = componentAlloc(firstComponent, PS_DEFAULT_TRACE_LEVEL); 203 } 204 currentNode = currentNode->subcomp[(currentNode->n) - 1]; 205 } 206 } 200 207 } 201 208 … … 211 218 zero 212 219 *****************************************************************************/ 213 boolpsTraceSetLevel(const char *comp, // component of interest214 220 int psTraceSetLevel(const char *comp, // component of interest 221 int level) // desired trace level 215 222 { 216 223 // If the root component tree does not exist, then initialize it. … … 219 226 } 220 227 // Add the new component to the component tree. 221 if ( !componentAdd(comp, level) ) { 222 psErrorMsg(ERRORNAME_PREFIX "psTraceSetLevel", PS_ERR_UNKNOWN, false, 223 PS_ERRORTEXT_psTrace_FAILED_TO_ADD_COMPONENT,level,comp); 224 return false; 225 } 228 componentAdd(comp, level); 229 226 230 // return 0 on success. 227 return true;231 return 0; 228 232 } 229 233 … … 249 253 p_psComponent* currentNode = cRoot; 250 254 int i = 0; 255 int defaultLevel = 0; 251 256 252 257 if (NULL == currentNode) { … … 262 267 } 263 268 269 defaultLevel = cRoot->level; 264 270 strcpy(name, aname); 265 271 pname = &name[1]; … … 268 274 for (i = 0; i < currentNode->n; i++) { 269 275 if (NULL == currentNode->subcomp[i]) { 270 psLogMsg(ERRORNAME_PREFIX "p_psTraceReset", PS_LOG_WARN, 271 PS_ERRORTEXT_psTrace_NULL_SUBCOMPONENT, 272 i, currentNode->name); 276 psError(__func__, 277 "Sub-component %d of node %s in trace tree is NULL.\n", i, currentNode->name); 273 278 } 274 279 275 280 if (strcmp(currentNode->subcomp[i]->name, firstComponent) == 0) { 276 281 currentNode = currentNode->subcomp[i]; 282 // For level inheritance purpose, we save the level of this 283 // component if it is not DEFAULT. 284 if (currentNode->level != PS_DEFAULT_TRACE_LEVEL) { 285 defaultLevel = currentNode->level; 286 } 287 // Determine if thisis the last component: 277 288 if (pname == NULL) { 278 return (currentNode->level); 289 if (currentNode->level != PS_DEFAULT_TRACE_LEVEL) { 290 return (currentNode->level); 291 } else { 292 return(PS_DEFAULT_TRACE_LEVEL); 293 } 279 294 } 280 295 } 281 296 } 282 297 } 283 return (PS_UNKNOWN_TRACE_LEVEL);284 } 285 286 /***************************************************************************** 287 ps GetTraceLevel()298 return(defaultLevel); 299 } 300 301 /***************************************************************************** 302 psTraceLevelGet() 288 303 Return a trace level of "name" in the root component tree. If the 289 304 exact string of components in "name" does not exist in the root … … 303 318 // Search the component root tree, determine the trace level. 304 319 return (doGetTraceLevel(name)); 320 } 321 322 /***************************************************************************** 323 doPrintTraceLevelsOld() 324 This function recursively searches the component tree supplied by the 325 parameter "comp" and prints the name and level of each component. 326 Inputs: 327 comp: a node in the component tree. 328 level: the level of that node 329 Outputs: 330 none 331 Returns: 332 null 333 *****************************************************************************/ 334 static void doPrintTraceLevelsOld(const p_psComponent* comp, int depth) 335 { 336 int i = 0; 337 338 if (comp->name[0] == '\0') { 339 printf("%*s%-*s %d\n", depth, "", 20 - depth, 340 "(root)", (comp->level == PS_DEFAULT_TRACE_LEVEL) ? 0 : comp->level); 341 } else { 342 if (comp->level == PS_DEFAULT_TRACE_LEVEL) { 343 printf("%*s%-*s %s\n", depth, "", 20 - depth, comp->name, "."); 344 } else { 345 printf("%*s%-*s %d\n", depth, "", 20 - depth, comp->name, 346 comp->level); 347 } 348 } 349 350 for (i = 0; i < comp->n; i++) { 351 doPrintTraceLevelsOld(comp->subcomp[i], depth + 1); 352 } 305 353 } 306 354 … … 317 365 null 318 366 *****************************************************************************/ 319 static void doPrintTraceLevels(const p_psComponent* comp, int depth) 367 static void doPrintTraceLevels(const p_psComponent* comp, 368 int depth, 369 int defLevel) 320 370 { 321 371 int i = 0; 322 372 323 373 if (comp->name[0] == '\0') { 324 printf("%*s%-*s %d\n", depth, "", 20 - depth, 325 "(root)", (comp->level == PS_UNKNOWN_TRACE_LEVEL) ? 0 : comp->level); 374 psAbort(__func__, "component name is NULL\n"); 326 375 } else { 327 if (comp->level == PS_UNKNOWN_TRACE_LEVEL) { 328 printf("%*s%-*s %s\n", depth, "", 20 - depth, comp->name, "."); 376 if (comp->level == PS_DEFAULT_TRACE_LEVEL) { 377 printf("%*s%-*s %d\n", depth, "", 20 - depth, comp->name, 378 defLevel); 329 379 } else { 330 printf("%*s%-*s %d\n", depth, "", 20 - depth, comp->name, comp->level); 380 printf("%*s%-*s %d\n", depth, "", 20 - depth, comp->name, 381 comp->level); 331 382 } 332 383 } 333 384 334 385 for (i = 0; i < comp->n; i++) { 335 doPrintTraceLevels(comp->subcomp[i], depth + 1); 336 } 337 } 386 if (comp->level == PS_DEFAULT_TRACE_LEVEL) { 387 doPrintTraceLevels(comp->subcomp[i], depth + 1, defLevel); 388 } else { 389 doPrintTraceLevels(comp->subcomp[i], depth + 1, comp->level); 390 } 391 } 392 } 393 338 394 339 395 /***************************************************************************** … … 353 409 } 354 410 355 doPrintTraceLevels(cRoot, 0 );411 doPrintTraceLevels(cRoot, 0, PS_DEFAULT_TRACE_LEVEL); 356 412 } 357 413 … … 379 435 380 436 if (NULL == comp) { 381 psErrorMsg(ERRORNAME_PREFIX "psTrace", PS_ERR_BAD_PARAMETER_NULL, true, 382 PS_ERRORTEXT_psTrace_NULL_TRACETREE, 383 __func__); 384 return; 437 psError(__func__, "p_psTrace() called on a NULL trace level tree\n"); 385 438 } 386 439 // Only display this message if it's trace level is less than the level
Note:
See TracChangeset
for help on using the changeset viewer.
