Changeset 1755 for trunk/psLib/src/math/psSpline.c

Timestamp:

Sep 8, 2004, 5:38:07 PM (22 years ago)

Author:

gusciora

Message:

...

File:

• trunk/psLib/src/math/psSpline.c (modified) (8 diffs)

trunk/psLib/src/math/psSpline.c

@@ -7 +7 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-08 23:32:03 $
+ *  @version $Revision: 1.30 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-09-09 03:38:07 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -110 +110 @@
 /*****************************************************************************/
 
-// None
+static psPolynomial1D* *CreateChebyshevPolys(int maxChebyPoly)
+{
+    psPolynomial1D **chebPolys = NULL;
+    int i = 0;
+    int j = 0;
+
+    chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly * sizeof(psPolynomial1D *));
+    for (i = 0; i < maxChebyPoly; i++) {
+        chebPolys[i] = psPolynomial1DAlloc(i + 1, PS_POLYNOMIAL_ORD);
+    }
+
+    // Create the Chebyshev polynomials.
+    // Polynomial i has i-th order.
+    chebPolys[0]->coeff[0] = 1;
+    chebPolys[1]->coeff[1] = 1;
+    for (i = 2; i < maxChebyPoly; i++) {
+        for (j = 0; j < chebPolys[i - 1]->n; j++) {
+            chebPolys[i]->coeff[j + 1] = 2 * chebPolys[i - 1]->coeff[j];
+        }
+        for (j = 0; j < chebPolys[i - 2]->n; j++) {
+            chebPolys[i]->coeff[j] -= chebPolys[i - 2]->coeff[j];
+        }
+    }
+
+    return (chebPolys);
+}
 
 /*****************************************************************************/
@@ -523 +548 @@
 float p_psChebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
 {
-    // XXX: Waiting for Crenshaw extension to multi-dimensional Chebyshevs.
-    psLogMsg(__func__, PS_LOG_WARN, "This function has not been coded.");
-    return(0.0);
+    int loop_x = 0;
+    int loop_y = 0;
+    float polySum = 0.0;
+    psPolynomial1D* *chebPolys = NULL;
+    int maxChebyPoly = 0;
+
+    // Determine how many Chebyshev polynomials
+    // are needed, then create them.
+    maxChebyPoly = myPoly->nX;
+    if (myPoly->nY > maxChebyPoly) {
+        maxChebyPoly = myPoly->nY;
+    }
+    chebPolys = CreateChebyshevPolys(maxChebyPoly);
+
+    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
+        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
+            if (myPoly->mask[loop_x][loop_y] != 0) {
+                polySum += myPoly->coeff[loop_x][loop_y] *
+                           psPolynomial1DEval(x, chebPolys[loop_x]) *
+                           psPolynomial1DEval(y, chebPolys[loop_y]);
+            }
+        }
+    }
+
+    psFree(chebPolys);
+    return(polySum);
 }
 
@@ -606 +654 @@
 float p_psChebPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly)
 {
-    // XXX: Waiting for Crenshaw extension to multi-dimensional Chebyshevs.
-    psLogMsg(__func__, PS_LOG_WARN, "This function has not been coded.");
-    return(0.0);
+    int loop_x = 0;
+    int loop_y = 0;
+    int loop_z = 0;
+    float polySum = 0.0;
+    psPolynomial1D* *chebPolys = NULL;
+    int maxChebyPoly = 0;
+
+    // Determine how many Chebyshev polynomials
+    // are needed, then create them.
+    maxChebyPoly = myPoly->nX;
+    if (myPoly->nY > maxChebyPoly) {
+        maxChebyPoly = myPoly->nY;
+    }
+    if (myPoly->nZ > maxChebyPoly) {
+        maxChebyPoly = myPoly->nZ;
+    }
+    chebPolys = CreateChebyshevPolys(maxChebyPoly);
+
+    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
+        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
+            for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
+                if (myPoly->mask[loop_x][loop_y][loop_z] != 0) {
+                    polySum += myPoly->coeff[loop_x][loop_y][loop_z] *
+                               psPolynomial1DEval(x, chebPolys[loop_x]) *
+                               psPolynomial1DEval(y, chebPolys[loop_y]) *
+                               psPolynomial1DEval(z, chebPolys[loop_z]);
+                }
+            }
+        }
+    }
+
+    psFree(chebPolys);
+    return(polySum);
 }
 
@@ -707 +785 @@
 float p_psChebPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly)
 {
-    // XXX: Waiting for Crenshaw extension to multi-dimensional Chebyshevs.
-    psLogMsg(__func__, PS_LOG_WARN, "This function has not been coded.");
-    return(0.0);
+    int loop_w = 0;
+    int loop_x = 0;
+    int loop_y = 0;
+    int loop_z = 0;
+    float polySum = 0.0;
+    psPolynomial1D* *chebPolys = NULL;
+    int maxChebyPoly = 0;
+
+    // Determine how many Chebyshev polynomials
+    // are needed, then create them.
+    maxChebyPoly = myPoly->nW;
+    if (myPoly->nX > maxChebyPoly) {
+        maxChebyPoly = myPoly->nX;
+    }
+    if (myPoly->nY > maxChebyPoly) {
+        maxChebyPoly = myPoly->nY;
+    }
+    if (myPoly->nZ > maxChebyPoly) {
+        maxChebyPoly = myPoly->nZ;
+    }
+    chebPolys = CreateChebyshevPolys(maxChebyPoly);
+
+    for (loop_w = 0; loop_w < myPoly->nW; loop_w++) {
+        for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
+            for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
+                for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
+                    if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) {
+                        polySum += myPoly->coeff[loop_w][loop_x][loop_y][loop_z] *
+                                   psPolynomial1DEval(w, chebPolys[loop_w]) *
+                                   psPolynomial1DEval(x, chebPolys[loop_x]) *
+                                   psPolynomial1DEval(y, chebPolys[loop_y]) *
+                                   psPolynomial1DEval(z, chebPolys[loop_z]);
+                    }
+                }
+            }
+        }
+    }
+
+    psFree(chebPolys);
+    return(polySum);
 }
 
@@ -1108 +1223 @@
 double p_psDChebPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly)
 {
-    // XXX: Waiting for Crenshaw extension to multi-dimensional Chebyshevs.
-    psLogMsg(__func__, PS_LOG_WARN, "This function has not been coded.");
-    return(0.0);
+    int loop_x = 0;
+    int loop_y = 0;
+    double polySum = 0.0;
+    psPolynomial1D* *chebPolys = NULL;
+    int maxChebyPoly = 0;
+
+    // Determine how many Chebyshev polynomials
+    // are needed, then create them.
+    maxChebyPoly = myPoly->nX;
+    if (myPoly->nY > maxChebyPoly) {
+        maxChebyPoly = myPoly->nY;
+    }
+    chebPolys = CreateChebyshevPolys(maxChebyPoly);
+
+    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
+        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
+            if (myPoly->mask[loop_x][loop_y] != 0) {
+                polySum += myPoly->coeff[loop_x][loop_y] *
+                           psPolynomial1DEval(x, chebPolys[loop_x]) *
+                           psPolynomial1DEval(y, chebPolys[loop_y]);
+            }
+        }
+    }
+
+    psFree(chebPolys);
+    return(polySum);
 }
 
@@ -1190 +1328 @@
 double p_psDChebPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly)
 {
-    // XXX: Waiting for Crenshaw extension to multi-dimensional Chebyshevs.
-    psLogMsg(__func__, PS_LOG_WARN, "This function has not been coded.");
-    return(0.0);
+    int loop_x = 0;
+    int loop_y = 0;
+    int loop_z = 0;
+    double polySum = 0.0;
+    psPolynomial1D* *chebPolys = NULL;
+    int maxChebyPoly = 0;
+
+    // Determine how many Chebyshev polynomials
+    // are needed, then create them.
+    maxChebyPoly = myPoly->nX;
+    if (myPoly->nY > maxChebyPoly) {
+        maxChebyPoly = myPoly->nY;
+    }
+    if (myPoly->nZ > maxChebyPoly) {
+        maxChebyPoly = myPoly->nZ;
+    }
+    chebPolys = CreateChebyshevPolys(maxChebyPoly);
+
+    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
+        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
+            for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
+                if (myPoly->mask[loop_x][loop_y][loop_z] != 0) {
+                    polySum += myPoly->coeff[loop_x][loop_y][loop_z] *
+                               psPolynomial1DEval(x, chebPolys[loop_x]) *
+                               psPolynomial1DEval(y, chebPolys[loop_y]) *
+                               psPolynomial1DEval(z, chebPolys[loop_z]);
+                }
+            }
+        }
+    }
+
+    psFree(chebPolys);
+    return(polySum);
 }
 
@@ -1293 +1461 @@
 double p_psDChebPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly)
 {
-    // XXX: Waiting for Crenshaw extension to multi-dimensional Chebyshevs.
-    psLogMsg(__func__, PS_LOG_WARN, "This function has not been coded.");
-    return(0.0);
+    int loop_w = 0;
+    int loop_x = 0;
+    int loop_y = 0;
+    int loop_z = 0;
+    double polySum = 0.0;
+    psPolynomial1D* *chebPolys = NULL;
+    int maxChebyPoly = 0;
+
+    // Determine how many Chebyshev polynomials
+    // are needed, then create them.
+    maxChebyPoly = myPoly->nW;
+    if (myPoly->nX > maxChebyPoly) {
+        maxChebyPoly = myPoly->nX;
+    }
+    if (myPoly->nY > maxChebyPoly) {
+        maxChebyPoly = myPoly->nY;
+    }
+    if (myPoly->nZ > maxChebyPoly) {
+        maxChebyPoly = myPoly->nZ;
+    }
+    chebPolys = CreateChebyshevPolys(maxChebyPoly);
+
+    for (loop_w = 0; loop_w < myPoly->nW; loop_w++) {
+        for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
+            for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
+                for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
+                    if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) {
+                        polySum += myPoly->coeff[loop_w][loop_x][loop_y][loop_z] *
+                                   psPolynomial1DEval(w, chebPolys[loop_w]) *
+                                   psPolynomial1DEval(x, chebPolys[loop_x]) *
+                                   psPolynomial1DEval(y, chebPolys[loop_y]) *
+                                   psPolynomial1DEval(z, chebPolys[loop_z]);
+                    }
+                }
+            }
+        }
+    }
+
+    psFree(chebPolys);
+    return(polySum);
 }