IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 1755


Ignore:
Timestamp:
Sep 8, 2004, 5:38:07 PM (22 years ago)
Author:
gusciora
Message:

...

Location:
trunk/psLib/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/dataManip/psFunctions.c

    r1734 r1755  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-09-08 23:32:03 $
     9 *  @version $Revision: 1.30 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-09-09 03:38:07 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    110110/*****************************************************************************/
    111111
    112 // None
     112static psPolynomial1D* *CreateChebyshevPolys(int maxChebyPoly)
     113{
     114    psPolynomial1D **chebPolys = NULL;
     115    int i = 0;
     116    int j = 0;
     117
     118    chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly * sizeof(psPolynomial1D *));
     119    for (i = 0; i < maxChebyPoly; i++) {
     120        chebPolys[i] = psPolynomial1DAlloc(i + 1, PS_POLYNOMIAL_ORD);
     121    }
     122
     123    // Create the Chebyshev polynomials.
     124    // Polynomial i has i-th order.
     125    chebPolys[0]->coeff[0] = 1;
     126    chebPolys[1]->coeff[1] = 1;
     127    for (i = 2; i < maxChebyPoly; i++) {
     128        for (j = 0; j < chebPolys[i - 1]->n; j++) {
     129            chebPolys[i]->coeff[j + 1] = 2 * chebPolys[i - 1]->coeff[j];
     130        }
     131        for (j = 0; j < chebPolys[i - 2]->n; j++) {
     132            chebPolys[i]->coeff[j] -= chebPolys[i - 2]->coeff[j];
     133        }
     134    }
     135
     136    return (chebPolys);
     137}
    113138
    114139/*****************************************************************************/
    … …  
    523548float p_psChebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
    524549{
    525     // XXX: Waiting for Crenshaw extension to multi-dimensional Chebyshevs.
    526     psLogMsg(__func__, PS_LOG_WARN, "This function has not been coded.");
    527     return(0.0);
     550    int loop_x = 0;
     551    int loop_y = 0;
     552    float polySum = 0.0;
     553    psPolynomial1D* *chebPolys = NULL;
     554    int maxChebyPoly = 0;
     555
     556    // Determine how many Chebyshev polynomials
     557    // are needed, then create them.
     558    maxChebyPoly = myPoly->nX;
     559    if (myPoly->nY > maxChebyPoly) {
     560        maxChebyPoly = myPoly->nY;
     561    }
     562    chebPolys = CreateChebyshevPolys(maxChebyPoly);
     563
     564    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     565        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
     566            if (myPoly->mask[loop_x][loop_y] != 0) {
     567                polySum += myPoly->coeff[loop_x][loop_y] *
     568                           psPolynomial1DEval(x, chebPolys[loop_x]) *
     569                           psPolynomial1DEval(y, chebPolys[loop_y]);
     570            }
     571        }
     572    }
     573
     574    psFree(chebPolys);
     575    return(polySum);
    528576}
    529577
    … …  
    606654float p_psChebPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly)
    607655{
    608     // XXX: Waiting for Crenshaw extension to multi-dimensional Chebyshevs.
    609     psLogMsg(__func__, PS_LOG_WARN, "This function has not been coded.");
    610     return(0.0);
     656    int loop_x = 0;
     657    int loop_y = 0;
     658    int loop_z = 0;
     659    float polySum = 0.0;
     660    psPolynomial1D* *chebPolys = NULL;
     661    int maxChebyPoly = 0;
     662
     663    // Determine how many Chebyshev polynomials
     664    // are needed, then create them.
     665    maxChebyPoly = myPoly->nX;
     666    if (myPoly->nY > maxChebyPoly) {
     667        maxChebyPoly = myPoly->nY;
     668    }
     669    if (myPoly->nZ > maxChebyPoly) {
     670        maxChebyPoly = myPoly->nZ;
     671    }
     672    chebPolys = CreateChebyshevPolys(maxChebyPoly);
     673
     674    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     675        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
     676            for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
     677                if (myPoly->mask[loop_x][loop_y][loop_z] != 0) {
     678                    polySum += myPoly->coeff[loop_x][loop_y][loop_z] *
     679                               psPolynomial1DEval(x, chebPolys[loop_x]) *
     680                               psPolynomial1DEval(y, chebPolys[loop_y]) *
     681                               psPolynomial1DEval(z, chebPolys[loop_z]);
     682                }
     683            }
     684        }
     685    }
     686
     687    psFree(chebPolys);
     688    return(polySum);
    611689}
    612690
    … …  
    707785float p_psChebPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly)
    708786{
    709     // XXX: Waiting for Crenshaw extension to multi-dimensional Chebyshevs.
    710     psLogMsg(__func__, PS_LOG_WARN, "This function has not been coded.");
    711     return(0.0);
     787    int loop_w = 0;
     788    int loop_x = 0;
     789    int loop_y = 0;
     790    int loop_z = 0;
     791    float polySum = 0.0;
     792    psPolynomial1D* *chebPolys = NULL;
     793    int maxChebyPoly = 0;
     794
     795    // Determine how many Chebyshev polynomials
     796    // are needed, then create them.
     797    maxChebyPoly = myPoly->nW;
     798    if (myPoly->nX > maxChebyPoly) {
     799        maxChebyPoly = myPoly->nX;
     800    }
     801    if (myPoly->nY > maxChebyPoly) {
     802        maxChebyPoly = myPoly->nY;
     803    }
     804    if (myPoly->nZ > maxChebyPoly) {
     805        maxChebyPoly = myPoly->nZ;
     806    }
     807    chebPolys = CreateChebyshevPolys(maxChebyPoly);
     808
     809    for (loop_w = 0; loop_w < myPoly->nW; loop_w++) {
     810        for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     811            for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
     812                for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
     813                    if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) {
     814                        polySum += myPoly->coeff[loop_w][loop_x][loop_y][loop_z] *
     815                                   psPolynomial1DEval(w, chebPolys[loop_w]) *
     816                                   psPolynomial1DEval(x, chebPolys[loop_x]) *
     817                                   psPolynomial1DEval(y, chebPolys[loop_y]) *
     818                                   psPolynomial1DEval(z, chebPolys[loop_z]);
     819                    }
     820                }
     821            }
     822        }
     823    }
     824
     825    psFree(chebPolys);
     826    return(polySum);
    712827}
    713828
    … …  
    11081223double p_psDChebPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly)
    11091224{
    1110     // XXX: Waiting for Crenshaw extension to multi-dimensional Chebyshevs.
    1111     psLogMsg(__func__, PS_LOG_WARN, "This function has not been coded.");
    1112     return(0.0);
     1225    int loop_x = 0;
     1226    int loop_y = 0;
     1227    double polySum = 0.0;
     1228    psPolynomial1D* *chebPolys = NULL;
     1229    int maxChebyPoly = 0;
     1230
     1231    // Determine how many Chebyshev polynomials
     1232    // are needed, then create them.
     1233    maxChebyPoly = myPoly->nX;
     1234    if (myPoly->nY > maxChebyPoly) {
     1235        maxChebyPoly = myPoly->nY;
     1236    }
     1237    chebPolys = CreateChebyshevPolys(maxChebyPoly);
     1238
     1239    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     1240        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
     1241            if (myPoly->mask[loop_x][loop_y] != 0) {
     1242                polySum += myPoly->coeff[loop_x][loop_y] *
     1243                           psPolynomial1DEval(x, chebPolys[loop_x]) *
     1244                           psPolynomial1DEval(y, chebPolys[loop_y]);
     1245            }
     1246        }
     1247    }
     1248
     1249    psFree(chebPolys);
     1250    return(polySum);
    11131251}
    11141252
    … …  
    11901328double p_psDChebPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly)
    11911329{
    1192     // XXX: Waiting for Crenshaw extension to multi-dimensional Chebyshevs.
    1193     psLogMsg(__func__, PS_LOG_WARN, "This function has not been coded.");
    1194     return(0.0);
     1330    int loop_x = 0;
     1331    int loop_y = 0;
     1332    int loop_z = 0;
     1333    double polySum = 0.0;
     1334    psPolynomial1D* *chebPolys = NULL;
     1335    int maxChebyPoly = 0;
     1336
     1337    // Determine how many Chebyshev polynomials
     1338    // are needed, then create them.
     1339    maxChebyPoly = myPoly->nX;
     1340    if (myPoly->nY > maxChebyPoly) {
     1341        maxChebyPoly = myPoly->nY;
     1342    }
     1343    if (myPoly->nZ > maxChebyPoly) {
     1344        maxChebyPoly = myPoly->nZ;
     1345    }
     1346    chebPolys = CreateChebyshevPolys(maxChebyPoly);
     1347
     1348    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     1349        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
     1350            for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
     1351                if (myPoly->mask[loop_x][loop_y][loop_z] != 0) {
     1352                    polySum += myPoly->coeff[loop_x][loop_y][loop_z] *
     1353                               psPolynomial1DEval(x, chebPolys[loop_x]) *
     1354                               psPolynomial1DEval(y, chebPolys[loop_y]) *
     1355                               psPolynomial1DEval(z, chebPolys[loop_z]);
     1356                }
     1357            }
     1358        }
     1359    }
     1360
     1361    psFree(chebPolys);
     1362    return(polySum);
    11951363}
    11961364
    … …  
    12931461double p_psDChebPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly)
    12941462{
    1295     // XXX: Waiting for Crenshaw extension to multi-dimensional Chebyshevs.
    1296     psLogMsg(__func__, PS_LOG_WARN, "This function has not been coded.");
    1297     return(0.0);
     1463    int loop_w = 0;
     1464    int loop_x = 0;
     1465    int loop_y = 0;
     1466    int loop_z = 0;
     1467    double polySum = 0.0;
     1468    psPolynomial1D* *chebPolys = NULL;
     1469    int maxChebyPoly = 0;
     1470
     1471    // Determine how many Chebyshev polynomials
     1472    // are needed, then create them.
     1473    maxChebyPoly = myPoly->nW;
     1474    if (myPoly->nX > maxChebyPoly) {
     1475        maxChebyPoly = myPoly->nX;
     1476    }
     1477    if (myPoly->nY > maxChebyPoly) {
     1478        maxChebyPoly = myPoly->nY;
     1479    }
     1480    if (myPoly->nZ > maxChebyPoly) {
     1481        maxChebyPoly = myPoly->nZ;
     1482    }
     1483    chebPolys = CreateChebyshevPolys(maxChebyPoly);
     1484
     1485    for (loop_w = 0; loop_w < myPoly->nW; loop_w++) {
     1486        for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     1487            for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
     1488                for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
     1489                    if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) {
     1490                        polySum += myPoly->coeff[loop_w][loop_x][loop_y][loop_z] *
     1491                                   psPolynomial1DEval(w, chebPolys[loop_w]) *
     1492                                   psPolynomial1DEval(x, chebPolys[loop_x]) *
     1493                                   psPolynomial1DEval(y, chebPolys[loop_y]) *
     1494                                   psPolynomial1DEval(z, chebPolys[loop_z]);
     1495                    }
     1496                }
     1497            }
     1498        }
     1499    }
     1500
     1501    psFree(chebPolys);
     1502    return(polySum);
    12981503}
    12991504

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

    r1734 r1755  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-09-08 23:32:03 $
     9 *  @version $Revision: 1.30 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-09-09 03:38:07 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    110110/*****************************************************************************/
    111111
    112 // None
     112static psPolynomial1D* *CreateChebyshevPolys(int maxChebyPoly)
     113{
     114    psPolynomial1D **chebPolys = NULL;
     115    int i = 0;
     116    int j = 0;
     117
     118    chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly * sizeof(psPolynomial1D *));
     119    for (i = 0; i < maxChebyPoly; i++) {
     120        chebPolys[i] = psPolynomial1DAlloc(i + 1, PS_POLYNOMIAL_ORD);
     121    }
     122
     123    // Create the Chebyshev polynomials.
     124    // Polynomial i has i-th order.
     125    chebPolys[0]->coeff[0] = 1;
     126    chebPolys[1]->coeff[1] = 1;
     127    for (i = 2; i < maxChebyPoly; i++) {
     128        for (j = 0; j < chebPolys[i - 1]->n; j++) {
     129            chebPolys[i]->coeff[j + 1] = 2 * chebPolys[i - 1]->coeff[j];
     130        }
     131        for (j = 0; j < chebPolys[i - 2]->n; j++) {
     132            chebPolys[i]->coeff[j] -= chebPolys[i - 2]->coeff[j];
     133        }
     134    }
     135
     136    return (chebPolys);
     137}
    113138
    114139/*****************************************************************************/
    … …  
    523548float p_psChebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
    524549{
    525     // XXX: Waiting for Crenshaw extension to multi-dimensional Chebyshevs.
    526     psLogMsg(__func__, PS_LOG_WARN, "This function has not been coded.");
    527     return(0.0);
     550    int loop_x = 0;
     551    int loop_y = 0;
     552    float polySum = 0.0;
     553    psPolynomial1D* *chebPolys = NULL;
     554    int maxChebyPoly = 0;
     555
     556    // Determine how many Chebyshev polynomials
     557    // are needed, then create them.
     558    maxChebyPoly = myPoly->nX;
     559    if (myPoly->nY > maxChebyPoly) {
     560        maxChebyPoly = myPoly->nY;
     561    }
     562    chebPolys = CreateChebyshevPolys(maxChebyPoly);
     563
     564    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     565        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
     566            if (myPoly->mask[loop_x][loop_y] != 0) {
     567                polySum += myPoly->coeff[loop_x][loop_y] *
     568                           psPolynomial1DEval(x, chebPolys[loop_x]) *
     569                           psPolynomial1DEval(y, chebPolys[loop_y]);
     570            }
     571        }
     572    }
     573
     574    psFree(chebPolys);
     575    return(polySum);
    528576}
    529577
    … …  
    606654float p_psChebPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly)
    607655{
    608     // XXX: Waiting for Crenshaw extension to multi-dimensional Chebyshevs.
    609     psLogMsg(__func__, PS_LOG_WARN, "This function has not been coded.");
    610     return(0.0);
     656    int loop_x = 0;
     657    int loop_y = 0;
     658    int loop_z = 0;
     659    float polySum = 0.0;
     660    psPolynomial1D* *chebPolys = NULL;
     661    int maxChebyPoly = 0;
     662
     663    // Determine how many Chebyshev polynomials
     664    // are needed, then create them.
     665    maxChebyPoly = myPoly->nX;
     666    if (myPoly->nY > maxChebyPoly) {
     667        maxChebyPoly = myPoly->nY;
     668    }
     669    if (myPoly->nZ > maxChebyPoly) {
     670        maxChebyPoly = myPoly->nZ;
     671    }
     672    chebPolys = CreateChebyshevPolys(maxChebyPoly);
     673
     674    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     675        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
     676            for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
     677                if (myPoly->mask[loop_x][loop_y][loop_z] != 0) {
     678                    polySum += myPoly->coeff[loop_x][loop_y][loop_z] *
     679                               psPolynomial1DEval(x, chebPolys[loop_x]) *
     680                               psPolynomial1DEval(y, chebPolys[loop_y]) *
     681                               psPolynomial1DEval(z, chebPolys[loop_z]);
     682                }
     683            }
     684        }
     685    }
     686
     687    psFree(chebPolys);
     688    return(polySum);
    611689}
    612690
    … …  
    707785float p_psChebPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly)
    708786{
    709     // XXX: Waiting for Crenshaw extension to multi-dimensional Chebyshevs.
    710     psLogMsg(__func__, PS_LOG_WARN, "This function has not been coded.");
    711     return(0.0);
     787    int loop_w = 0;
     788    int loop_x = 0;
     789    int loop_y = 0;
     790    int loop_z = 0;
     791    float polySum = 0.0;
     792    psPolynomial1D* *chebPolys = NULL;
     793    int maxChebyPoly = 0;
     794
     795    // Determine how many Chebyshev polynomials
     796    // are needed, then create them.
     797    maxChebyPoly = myPoly->nW;
     798    if (myPoly->nX > maxChebyPoly) {
     799        maxChebyPoly = myPoly->nX;
     800    }
     801    if (myPoly->nY > maxChebyPoly) {
     802        maxChebyPoly = myPoly->nY;
     803    }
     804    if (myPoly->nZ > maxChebyPoly) {
     805        maxChebyPoly = myPoly->nZ;
     806    }
     807    chebPolys = CreateChebyshevPolys(maxChebyPoly);
     808
     809    for (loop_w = 0; loop_w < myPoly->nW; loop_w++) {
     810        for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     811            for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
     812                for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
     813                    if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) {
     814                        polySum += myPoly->coeff[loop_w][loop_x][loop_y][loop_z] *
     815                                   psPolynomial1DEval(w, chebPolys[loop_w]) *
     816                                   psPolynomial1DEval(x, chebPolys[loop_x]) *
     817                                   psPolynomial1DEval(y, chebPolys[loop_y]) *
     818                                   psPolynomial1DEval(z, chebPolys[loop_z]);
     819                    }
     820                }
     821            }
     822        }
     823    }
     824
     825    psFree(chebPolys);
     826    return(polySum);
    712827}
    713828
    … …  
    11081223double p_psDChebPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly)
    11091224{
    1110     // XXX: Waiting for Crenshaw extension to multi-dimensional Chebyshevs.
    1111     psLogMsg(__func__, PS_LOG_WARN, "This function has not been coded.");
    1112     return(0.0);
     1225    int loop_x = 0;
     1226    int loop_y = 0;
     1227    double polySum = 0.0;
     1228    psPolynomial1D* *chebPolys = NULL;
     1229    int maxChebyPoly = 0;
     1230
     1231    // Determine how many Chebyshev polynomials
     1232    // are needed, then create them.
     1233    maxChebyPoly = myPoly->nX;
     1234    if (myPoly->nY > maxChebyPoly) {
     1235        maxChebyPoly = myPoly->nY;
     1236    }
     1237    chebPolys = CreateChebyshevPolys(maxChebyPoly);
     1238
     1239    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     1240        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
     1241            if (myPoly->mask[loop_x][loop_y] != 0) {
     1242                polySum += myPoly->coeff[loop_x][loop_y] *
     1243                           psPolynomial1DEval(x, chebPolys[loop_x]) *
     1244                           psPolynomial1DEval(y, chebPolys[loop_y]);
     1245            }
     1246        }
     1247    }
     1248
     1249    psFree(chebPolys);
     1250    return(polySum);
    11131251}
    11141252
    … …  
    11901328double p_psDChebPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly)
    11911329{
    1192     // XXX: Waiting for Crenshaw extension to multi-dimensional Chebyshevs.
    1193     psLogMsg(__func__, PS_LOG_WARN, "This function has not been coded.");
    1194     return(0.0);
     1330    int loop_x = 0;
     1331    int loop_y = 0;
     1332    int loop_z = 0;
     1333    double polySum = 0.0;
     1334    psPolynomial1D* *chebPolys = NULL;
     1335    int maxChebyPoly = 0;
     1336
     1337    // Determine how many Chebyshev polynomials
     1338    // are needed, then create them.
     1339    maxChebyPoly = myPoly->nX;
     1340    if (myPoly->nY > maxChebyPoly) {
     1341        maxChebyPoly = myPoly->nY;
     1342    }
     1343    if (myPoly->nZ > maxChebyPoly) {
     1344        maxChebyPoly = myPoly->nZ;
     1345    }
     1346    chebPolys = CreateChebyshevPolys(maxChebyPoly);
     1347
     1348    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     1349        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
     1350            for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
     1351                if (myPoly->mask[loop_x][loop_y][loop_z] != 0) {
     1352                    polySum += myPoly->coeff[loop_x][loop_y][loop_z] *
     1353                               psPolynomial1DEval(x, chebPolys[loop_x]) *
     1354                               psPolynomial1DEval(y, chebPolys[loop_y]) *
     1355                               psPolynomial1DEval(z, chebPolys[loop_z]);
     1356                }
     1357            }
     1358        }
     1359    }
     1360
     1361    psFree(chebPolys);
     1362    return(polySum);
    11951363}
    11961364
    … …  
    12931461double p_psDChebPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly)
    12941462{
    1295     // XXX: Waiting for Crenshaw extension to multi-dimensional Chebyshevs.
    1296     psLogMsg(__func__, PS_LOG_WARN, "This function has not been coded.");
    1297     return(0.0);
     1463    int loop_w = 0;
     1464    int loop_x = 0;
     1465    int loop_y = 0;
     1466    int loop_z = 0;
     1467    double polySum = 0.0;
     1468    psPolynomial1D* *chebPolys = NULL;
     1469    int maxChebyPoly = 0;
     1470
     1471    // Determine how many Chebyshev polynomials
     1472    // are needed, then create them.
     1473    maxChebyPoly = myPoly->nW;
     1474    if (myPoly->nX > maxChebyPoly) {
     1475        maxChebyPoly = myPoly->nX;
     1476    }
     1477    if (myPoly->nY > maxChebyPoly) {
     1478        maxChebyPoly = myPoly->nY;
     1479    }
     1480    if (myPoly->nZ > maxChebyPoly) {
     1481        maxChebyPoly = myPoly->nZ;
     1482    }
     1483    chebPolys = CreateChebyshevPolys(maxChebyPoly);
     1484
     1485    for (loop_w = 0; loop_w < myPoly->nW; loop_w++) {
     1486        for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     1487            for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
     1488                for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
     1489                    if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) {
     1490                        polySum += myPoly->coeff[loop_w][loop_x][loop_y][loop_z] *
     1491                                   psPolynomial1DEval(w, chebPolys[loop_w]) *
     1492                                   psPolynomial1DEval(x, chebPolys[loop_x]) *
     1493                                   psPolynomial1DEval(y, chebPolys[loop_y]) *
     1494                                   psPolynomial1DEval(z, chebPolys[loop_z]);
     1495                    }
     1496                }
     1497            }
     1498        }
     1499    }
     1500
     1501    psFree(chebPolys);
     1502    return(polySum);
    12981503}
    12991504

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

    r1734 r1755  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-09-08 23:32:03 $
     9 *  @version $Revision: 1.30 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-09-09 03:38:07 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    110110/*****************************************************************************/
    111111
    112 // None
     112static psPolynomial1D* *CreateChebyshevPolys(int maxChebyPoly)
     113{
     114    psPolynomial1D **chebPolys = NULL;
     115    int i = 0;
     116    int j = 0;
     117
     118    chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly * sizeof(psPolynomial1D *));
     119    for (i = 0; i < maxChebyPoly; i++) {
     120        chebPolys[i] = psPolynomial1DAlloc(i + 1, PS_POLYNOMIAL_ORD);
     121    }
     122
     123    // Create the Chebyshev polynomials.
     124    // Polynomial i has i-th order.
     125    chebPolys[0]->coeff[0] = 1;
     126    chebPolys[1]->coeff[1] = 1;
     127    for (i = 2; i < maxChebyPoly; i++) {
     128        for (j = 0; j < chebPolys[i - 1]->n; j++) {
     129            chebPolys[i]->coeff[j + 1] = 2 * chebPolys[i - 1]->coeff[j];
     130        }
     131        for (j = 0; j < chebPolys[i - 2]->n; j++) {
     132            chebPolys[i]->coeff[j] -= chebPolys[i - 2]->coeff[j];
     133        }
     134    }
     135
     136    return (chebPolys);
     137}
    113138
    114139/*****************************************************************************/
    … …  
    523548float p_psChebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
    524549{
    525     // XXX: Waiting for Crenshaw extension to multi-dimensional Chebyshevs.
    526     psLogMsg(__func__, PS_LOG_WARN, "This function has not been coded.");
    527     return(0.0);
     550    int loop_x = 0;
     551    int loop_y = 0;
     552    float polySum = 0.0;
     553    psPolynomial1D* *chebPolys = NULL;
     554    int maxChebyPoly = 0;
     555
     556    // Determine how many Chebyshev polynomials
     557    // are needed, then create them.
     558    maxChebyPoly = myPoly->nX;
     559    if (myPoly->nY > maxChebyPoly) {
     560        maxChebyPoly = myPoly->nY;
     561    }
     562    chebPolys = CreateChebyshevPolys(maxChebyPoly);
     563
     564    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     565        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
     566            if (myPoly->mask[loop_x][loop_y] != 0) {
     567                polySum += myPoly->coeff[loop_x][loop_y] *
     568                           psPolynomial1DEval(x, chebPolys[loop_x]) *
     569                           psPolynomial1DEval(y, chebPolys[loop_y]);
     570            }
     571        }
     572    }
     573
     574    psFree(chebPolys);
     575    return(polySum);
    528576}
    529577
    … …  
    606654float p_psChebPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly)
    607655{
    608     // XXX: Waiting for Crenshaw extension to multi-dimensional Chebyshevs.
    609     psLogMsg(__func__, PS_LOG_WARN, "This function has not been coded.");
    610     return(0.0);
     656    int loop_x = 0;
     657    int loop_y = 0;
     658    int loop_z = 0;
     659    float polySum = 0.0;
     660    psPolynomial1D* *chebPolys = NULL;
     661    int maxChebyPoly = 0;
     662
     663    // Determine how many Chebyshev polynomials
     664    // are needed, then create them.
     665    maxChebyPoly = myPoly->nX;
     666    if (myPoly->nY > maxChebyPoly) {
     667        maxChebyPoly = myPoly->nY;
     668    }
     669    if (myPoly->nZ > maxChebyPoly) {
     670        maxChebyPoly = myPoly->nZ;
     671    }
     672    chebPolys = CreateChebyshevPolys(maxChebyPoly);
     673
     674    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     675        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
     676            for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
     677                if (myPoly->mask[loop_x][loop_y][loop_z] != 0) {
     678                    polySum += myPoly->coeff[loop_x][loop_y][loop_z] *
     679                               psPolynomial1DEval(x, chebPolys[loop_x]) *
     680                               psPolynomial1DEval(y, chebPolys[loop_y]) *
     681                               psPolynomial1DEval(z, chebPolys[loop_z]);
     682                }
     683            }
     684        }
     685    }
     686
     687    psFree(chebPolys);
     688    return(polySum);
    611689}
    612690
    … …  
    707785float p_psChebPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly)
    708786{
    709     // XXX: Waiting for Crenshaw extension to multi-dimensional Chebyshevs.
    710     psLogMsg(__func__, PS_LOG_WARN, "This function has not been coded.");
    711     return(0.0);
     787    int loop_w = 0;
     788    int loop_x = 0;
     789    int loop_y = 0;
     790    int loop_z = 0;
     791    float polySum = 0.0;
     792    psPolynomial1D* *chebPolys = NULL;
     793    int maxChebyPoly = 0;
     794
     795    // Determine how many Chebyshev polynomials
     796    // are needed, then create them.
     797    maxChebyPoly = myPoly->nW;
     798    if (myPoly->nX > maxChebyPoly) {
     799        maxChebyPoly = myPoly->nX;
     800    }
     801    if (myPoly->nY > maxChebyPoly) {
     802        maxChebyPoly = myPoly->nY;
     803    }
     804    if (myPoly->nZ > maxChebyPoly) {
     805        maxChebyPoly = myPoly->nZ;
     806    }
     807    chebPolys = CreateChebyshevPolys(maxChebyPoly);
     808
     809    for (loop_w = 0; loop_w < myPoly->nW; loop_w++) {
     810        for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     811            for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
     812                for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
     813                    if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) {
     814                        polySum += myPoly->coeff[loop_w][loop_x][loop_y][loop_z] *
     815                                   psPolynomial1DEval(w, chebPolys[loop_w]) *
     816                                   psPolynomial1DEval(x, chebPolys[loop_x]) *
     817                                   psPolynomial1DEval(y, chebPolys[loop_y]) *
     818                                   psPolynomial1DEval(z, chebPolys[loop_z]);
     819                    }
     820                }
     821            }
     822        }
     823    }
     824
     825    psFree(chebPolys);
     826    return(polySum);
    712827}
    713828
    … …  
    11081223double p_psDChebPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly)
    11091224{
    1110     // XXX: Waiting for Crenshaw extension to multi-dimensional Chebyshevs.
    1111     psLogMsg(__func__, PS_LOG_WARN, "This function has not been coded.");
    1112     return(0.0);
     1225    int loop_x = 0;
     1226    int loop_y = 0;
     1227    double polySum = 0.0;
     1228    psPolynomial1D* *chebPolys = NULL;
     1229    int maxChebyPoly = 0;
     1230
     1231    // Determine how many Chebyshev polynomials
     1232    // are needed, then create them.
     1233    maxChebyPoly = myPoly->nX;
     1234    if (myPoly->nY > maxChebyPoly) {
     1235        maxChebyPoly = myPoly->nY;
     1236    }
     1237    chebPolys = CreateChebyshevPolys(maxChebyPoly);
     1238
     1239    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     1240        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
     1241            if (myPoly->mask[loop_x][loop_y] != 0) {
     1242                polySum += myPoly->coeff[loop_x][loop_y] *
     1243                           psPolynomial1DEval(x, chebPolys[loop_x]) *
     1244                           psPolynomial1DEval(y, chebPolys[loop_y]);
     1245            }
     1246        }
     1247    }
     1248
     1249    psFree(chebPolys);
     1250    return(polySum);
    11131251}
    11141252
    … …  
    11901328double p_psDChebPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly)
    11911329{
    1192     // XXX: Waiting for Crenshaw extension to multi-dimensional Chebyshevs.
    1193     psLogMsg(__func__, PS_LOG_WARN, "This function has not been coded.");
    1194     return(0.0);
     1330    int loop_x = 0;
     1331    int loop_y = 0;
     1332    int loop_z = 0;
     1333    double polySum = 0.0;
     1334    psPolynomial1D* *chebPolys = NULL;
     1335    int maxChebyPoly = 0;
     1336
     1337    // Determine how many Chebyshev polynomials
     1338    // are needed, then create them.
     1339    maxChebyPoly = myPoly->nX;
     1340    if (myPoly->nY > maxChebyPoly) {
     1341        maxChebyPoly = myPoly->nY;
     1342    }
     1343    if (myPoly->nZ > maxChebyPoly) {
     1344        maxChebyPoly = myPoly->nZ;
     1345    }
     1346    chebPolys = CreateChebyshevPolys(maxChebyPoly);
     1347
     1348    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     1349        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
     1350            for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
     1351                if (myPoly->mask[loop_x][loop_y][loop_z] != 0) {
     1352                    polySum += myPoly->coeff[loop_x][loop_y][loop_z] *
     1353                               psPolynomial1DEval(x, chebPolys[loop_x]) *
     1354                               psPolynomial1DEval(y, chebPolys[loop_y]) *
     1355                               psPolynomial1DEval(z, chebPolys[loop_z]);
     1356                }
     1357            }
     1358        }
     1359    }
     1360
     1361    psFree(chebPolys);
     1362    return(polySum);
    11951363}
    11961364
    … …  
    12931461double p_psDChebPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly)
    12941462{
    1295     // XXX: Waiting for Crenshaw extension to multi-dimensional Chebyshevs.
    1296     psLogMsg(__func__, PS_LOG_WARN, "This function has not been coded.");
    1297     return(0.0);
     1463    int loop_w = 0;
     1464    int loop_x = 0;
     1465    int loop_y = 0;
     1466    int loop_z = 0;
     1467    double polySum = 0.0;
     1468    psPolynomial1D* *chebPolys = NULL;
     1469    int maxChebyPoly = 0;
     1470
     1471    // Determine how many Chebyshev polynomials
     1472    // are needed, then create them.
     1473    maxChebyPoly = myPoly->nW;
     1474    if (myPoly->nX > maxChebyPoly) {
     1475        maxChebyPoly = myPoly->nX;
     1476    }
     1477    if (myPoly->nY > maxChebyPoly) {
     1478        maxChebyPoly = myPoly->nY;
     1479    }
     1480    if (myPoly->nZ > maxChebyPoly) {
     1481        maxChebyPoly = myPoly->nZ;
     1482    }
     1483    chebPolys = CreateChebyshevPolys(maxChebyPoly);
     1484
     1485    for (loop_w = 0; loop_w < myPoly->nW; loop_w++) {
     1486        for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     1487            for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
     1488                for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
     1489                    if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) {
     1490                        polySum += myPoly->coeff[loop_w][loop_x][loop_y][loop_z] *
     1491                                   psPolynomial1DEval(w, chebPolys[loop_w]) *
     1492                                   psPolynomial1DEval(x, chebPolys[loop_x]) *
     1493                                   psPolynomial1DEval(y, chebPolys[loop_y]) *
     1494                                   psPolynomial1DEval(z, chebPolys[loop_z]);
     1495                    }
     1496                }
     1497            }
     1498        }
     1499    }
     1500
     1501    psFree(chebPolys);
     1502    return(polySum);
    12981503}
    12991504
Note: See TracChangeset for help on using the changeset viewer.