IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 7, 2004, 8:00:49 PM (22 years ago)
Author:
gusciora
Message:

...

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/math/psPolynomial.c

    r1476 r1718  
    1 
    21/** @file  psFunctions.c
    32 *
    … …  
    87 *  polynomials.  It also contains a Gaussian functions.
    98 *
    10  *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2004-08-11 20:07:45 $
     9 *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-09-08 06:00:49 $
    1211 *
    1312 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    3231#include "psAbort.h"
    3332#include "psFunctions.h"
     33#include "psLogMsg.h"
    3434
    3535#include <gsl/gsl_rng.h>
    … …  
    4242/*****************************************************************************/
    4343
    44 // None
     44#define PS_CONVERT_VECTOR_F32(IN, OUT) \
     45if (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) \
     64if (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
    4583
    4684/*****************************************************************************/
    … …  
    107145    }
    108146
    109     return (tmp * exp(-((x - mean) * (x - mean)) / (2.0 * sigma * sigma)));
     147    return(tmp * exp(-((x - mean) * (x - mean)) / (2.0 * sigma * sigma)));
    110148}
    111149
    … …  
    138176
    139177    // 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;
    142186
    143187/*****************************************************************************
    144188    This routine must allocate memory for the polynomial structures.
    145189 *****************************************************************************/
    146 psPolynomial1D* psPolynomial1DAlloc(int n)
     190psPolynomial1D* psPolynomial1DAlloc(int n,
     191                                    psPolynomialType type)
    147192{
    148193    int i = 0;
    … …  
    151196    newPoly = (psPolynomial1D* ) psAlloc(sizeof(psPolynomial1D));
    152197    p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial1DFree);
     198
     199    newPoly->type = type;
    153200    newPoly->n = n;
    154201    newPoly->coeff = (float *)psAlloc(n * sizeof(float));
    … …  
    161208    }
    162209
    163     return (newPoly);
    164 }
    165 
    166 psPolynomial2D* psPolynomial2DAlloc(int nX, int nY)
     210    return(newPoly);
     211}
     212
     213psPolynomial2D* psPolynomial2DAlloc(int nX, int nY,
     214                                    psPolynomialType type)
    167215{
    168216    int x = 0;
    … …  
    172220    newPoly = (psPolynomial2D* ) psAlloc(sizeof(psPolynomial2D));
    173221    p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial2DFree);
     222
     223    newPoly->type = type;
    174224    newPoly->nX = nX;
    175225    newPoly->nY = nY;
    … …  
    191241    }
    192242
    193     return (newPoly);
    194 }
    195 
    196 psPolynomial3D* psPolynomial3DAlloc(int nX, int nY, int nZ)
     243    return(newPoly);
     244}
     245
     246psPolynomial3D* psPolynomial3DAlloc(int nX, int nY, int nZ,
     247                                    psPolynomialType type)
    197248{
    198249    int x = 0;
    … …  
    203254    newPoly = (psPolynomial3D* ) psAlloc(sizeof(psPolynomial3D));
    204255    p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial3DFree);
     256
     257    newPoly->type = type;
    205258    newPoly->nX = nX;
    206259    newPoly->nY = nY;
    … …  
    230283    }
    231284
    232     return (newPoly);
    233 }
    234 
    235 psPolynomial4D* psPolynomial4DAlloc(int nW, int nX, int nY, int nZ)
     285    return(newPoly);
     286}
     287
     288psPolynomial4D* psPolynomial4DAlloc(int nW, int nX, int nY, int nZ,
     289                                    psPolynomialType type)
    236290{
    237291    int w = 0;
    … …  
    243297    newPoly = (psPolynomial4D* ) psAlloc(sizeof(psPolynomial4D));
    244298    p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial4DFree);
     299
     300    newPoly->type = type;
    245301    newPoly->nW = nW;
    246302    newPoly->nX = nX;
    … …  
    278334    }
    279335
    280     return (newPoly);
     336    return(newPoly);
    281337}
    282338
    … …  
    359415    XXX: Should the "coeffErr[]" should be used as well?
    360416 *****************************************************************************/
    361 float psPolynomial1DEval(float x, const psPolynomial1D* myPoly)
     417float p_psOrgPolynomial1DEval(float x, const psPolynomial1D* myPoly)
    362418{
    363419    int loop_x = 0;
    … …  
    371427    // NOTE: Do we want to flag this case?
    372428    if (myPoly->n == 0) {
    373         return (1.0);
     429        return(1.0);
    374430    }
    375431
    … …  
    383439    }
    384440
    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.
     445float 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
     470float 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
     482psVector *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
     503float p_psOrgPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
    389504{
    390505    int loop_x = 0;
    … …  
    403518    }
    404519
    405     return (polySum);
    406 }
    407 
    408 float psPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly)
     520    return(polySum);
     521}
     522
     523float p_psChebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
     524{
     525    return(0.0);
     526}
     527
     528float 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
     541psVector *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
     576float p_psOrgPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly)
    409577{
    410578    int loop_x = 0;
    … …  
    429597    }
    430598
    431     return (polySum);
    432 }
    433 
    434 float psPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly)
     599    return(polySum);
     600}
     601
     602float p_psChebPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly)
     603{
     604    return(0.0);
     605}
     606
     607float 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
     619psVector *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
     667float p_psOrgPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly)
    435668{
    436669    int loop_w = 0;
    … …  
    461694    }
    462695
    463     return (polySum);
    464 }
    465 
    466 psDPolynomial1D* psDPolynomial1DAlloc(int n)
     696    return(polySum);
     697}
     698
     699float p_psChebPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly)
     700{
     701    return(0.0);
     702}
     703
     704float 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
     716psVector *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
     773psDPolynomial1D* psDPolynomial1DAlloc(int n,
     774                                      psPolynomialType type)
    467775{
    468776    int i = 0;
    … …  
    471779    newPoly = (psDPolynomial1D* ) psAlloc(sizeof(psDPolynomial1D));
    472780    p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial1DFree);
     781
     782    newPoly->type = type;
    473783    newPoly->n = n;
    474784    newPoly->coeff = (double *)psAlloc(n * sizeof(double));
    … …  
    481791    }
    482792
    483     return (newPoly);
    484 }
    485 
    486 psDPolynomial2D* psDPolynomial2DAlloc(int nX, int nY)
     793    return(newPoly);
     794}
     795
     796psDPolynomial2D* psDPolynomial2DAlloc(int nX, int nY,
     797                                      psPolynomialType type)
    487798{
    488799    int x = 0;
    … …  
    492803    newPoly = (psDPolynomial2D* ) psAlloc(sizeof(psDPolynomial2D));
    493804    p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial2DFree);
     805
     806    newPoly->type = type;
    494807    newPoly->nX = nX;
    495808    newPoly->nY = nY;
    … …  
    511824    }
    512825
    513     return (newPoly);
    514 }
    515 
    516 psDPolynomial3D* psDPolynomial3DAlloc(int nX, int nY, int nZ)
     826    return(newPoly);
     827}
     828
     829psDPolynomial3D* psDPolynomial3DAlloc(int nX, int nY, int nZ,
     830                                      psPolynomialType type)
    517831{
    518832    int x = 0;
    … …  
    523837    newPoly = (psDPolynomial3D* ) psAlloc(sizeof(psDPolynomial3D));
    524838    p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial3DFree);
     839
     840    newPoly->type = type;
    525841    newPoly->nX = nX;
    526842    newPoly->nY = nY;
    … …  
    550866    }
    551867
    552     return (newPoly);
    553 }
    554 
    555 psDPolynomial4D* psDPolynomial4DAlloc(int nW, int nX, int nY, int nZ)
     868    return(newPoly);
     869}
     870
     871psDPolynomial4D* psDPolynomial4DAlloc(int nW, int nX, int nY, int nZ,
     872                                      psPolynomialType type)
    556873{
    557874    int w = 0;
    … …  
    563880    newPoly = (psDPolynomial4D* ) psAlloc(sizeof(psDPolynomial4D));
    564881    p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial4DFree);
     882
     883    newPoly->type = type;
    565884    newPoly->nW = nW;
    566885    newPoly->nX = nX;
    … …  
    598917    }
    599918
    600     return (newPoly);
     919    return(newPoly);
    601920}
    602921
    … …  
    673992    Polynomial coefficients will be accessed in [w][x][y][z] fashion.
    674993 *****************************************************************************/
    675 double psDPolynomial1DEval(double x, const psDPolynomial1D* myPoly)
     994double p_psDOrgPolynomial1DEval(double x, const psDPolynomial1D* myPoly)
    676995{
    677996    int loop_x = 0;
    … …  
    6811000    // NOTE: Do we want to flag this case?
    6821001    if (myPoly->n == 0) {
    683         return (1.0);
     1002        return(1.0);
    6841003    }
    6851004
    … …  
    6891008    }
    6901009
    691     return (polySum);
    692 }
    693 
    694 double psDPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly)
     1010    return(polySum);
     1011}
     1012
     1013double p_psDChebPolynomial1DEval(double x, const psDPolynomial1D* myPoly)
     1014{
     1015    return(0.0);
     1016}
     1017
     1018double 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
     1030psVector *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
     1052double p_psDOrgPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly)
    6951053{
    6961054    int loop_x = 0;
    … …  
    7091067    }
    7101068
    711     return (polySum);
    712 }
    713 
    714 double psDPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly)
     1069    return(polySum);
     1070}
     1071
     1072double p_psDChebPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly)
     1073{
     1074    return(0.0);
     1075}
     1076
     1077double 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
     1089psVector *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
     1124double p_psDOrgPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly)
    7151125{
    7161126    int loop_x = 0;
    … …  
    7351145    }
    7361146
    737     return (polySum);
    738 }
    739 
    740 double psDPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly)
     1147    return(polySum);
     1148}
     1149
     1150double p_psDChebPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly)
     1151{
     1152    return(0.0);
     1153}
     1154
     1155double 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
     1167psVector *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
     1217double p_psDOrgPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly)
    7411218{
    7421219    int loop_w = 0;
    … …  
    7671244    }
    7681245
    769     return (polySum);
    770 }
     1246    return(polySum);
     1247}
     1248
     1249double p_psDChebPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly)
     1250{
     1251    return(0.0);
     1252}
     1253
     1254double 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
     1266psVector *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 *****************************************************************************/
     1333psSpline1D *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
     1363psSpline1D *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/*****************************************************************************
     1388VectorBinDisect(): This is a private function which takes as input a vector
     1389of floating point data as well as a single floating point values.  The input
     1390vector values are assumed to be non-decreasing (v[i-1] <= v[j] for all j>=i).
     1391This routine does a binary disection of the vector and returns "i" such
     1392that (v[i] <= x <= v[i+1).  If x lies outside the range of v[],
     1393then this routine prints a warning message and returns -1.
     1394 *****************************************************************************/
     1395int 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
     1428float 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
     1452psVector *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
Note: See TracChangeset for help on using the changeset viewer.