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

Timestamp:

Nov 15, 2005, 10:10:32 AM (21 years ago)

Author:

gusciora

Message:

Bug fix for psSpline, some new psConstant macros, ....

File:

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

trunk/psLib/src/math/psSpline.c

@@ -7 +7 @@
 *  splines.
 *
-*  @version $Revision: 1.131 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-09-29 20:06:58 $
+*  @version $Revision: 1.132 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-11-15 20:10:32 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -496 +496 @@
     PS_ASSERT_VECTOR_NON_NULL(y, NULL);
     PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL);
+    PS_ASSERT_INT_LARGER_THAN_OR_EQUAL(y->n, 2, NULL);
     psS32 numSplines = (y->n)-1;
     psTrace(__func__, 5, "numSplines is %d\n", numSplines);
+
+    //
+    // Create the psSpline1D struct.
+    //
+    psSpline1D *spline = psSpline1DAlloc();
+    spline->n = numSplines;
+    spline->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *));
+    for (psS32 i=0;i<numSplines;i++) {
+        spline->spline[i] = psPolynomial1DAlloc(3, PS_POLYNOMIAL_ORD);
+    }
+
     //
     // The following code ensures that xPtr and yPtr points to a psF32 psVector.
     //
-    psVector *xPtr = NULL;
+    // XXX: When you debug, use the vector copy and create routines here:
+    //
+
+    spline->knots = psVectorAlloc(y->n, PS_TYPE_F32);
     if (x != NULL) {
         PS_ASSERT_VECTOR_NON_NULL(x, NULL);
+        PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL);
         PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL);
-        PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL);
-        // Convert x to F32 if necessary.
-        if (PS_TYPE_F64 == x->type.type) {
-            xPtr = psVectorCopy(NULL, x, PS_TYPE_F32);
-        } else if (PS_TYPE_F32 == x->type.type) {
-            xPtr = (psVector *) x;
+        if (x->type.type == PS_TYPE_F32) {
+            for (psS32 i = 0 ; i < x->n ; i++) {
+                spline->knots->data.F32[i] = x->data.F32[i];
+            }
+        } else if (x->type.type == PS_TYPE_F64) {
+            for (psS32 i = 0 ; i < x->n ; i++) {
+                spline->knots->data.F32[i] = (psF32) x->data.F64[i];
+            }
         }
     } else {
-        // Allocate an index vector for x
-        xPtr = psVectorCreate(NULL, 0.0, (psF64) y->n, 1.0, PS_TYPE_F32);
-    }
+        for (psS32 i = 0 ; i < y->n ; i++) {
+            spline->knots->data.F32[i] = (psF32) i;
+        }
+    }
+    psVector *xPtr = spline->knots;
 
     psVector *yPtr = NULL;
@@ -526 +546 @@
 
     //
-    // Create the psSpline1D struct.
-    //
-    psSpline1D *spline = psSpline1DAlloc();
-    spline->n = numSplines;
-    spline->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *));
-    for (psS32 i=0;i<numSplines;i++) {
-        spline->spline[i] = psPolynomial1DAlloc(3, PS_POLYNOMIAL_ORD);
-    }
-    spline->knots = psVectorCopy(NULL, xPtr, PS_TYPE_F32);
-    //
     // Generate the second derivatives at each data point.
     //
@@ -548 +558 @@
         psF32 H = xPtr->data.F32[i+1] - xPtr->data.F32[i];
         if (fabs(H) <= FLT_EPSILON) {
-            psError(PS_ERR_UNKNOWN, false, "x data points are not distinct (%d %d).\n", i, i+1);
+            psError(PS_ERR_UNKNOWN, false, "x data points are not distinct (%d %d) (%f %f).\n",
+                    i, i+1, xPtr->data.F32[i], xPtr->data.F32[i+1]);
         }
         psTrace(__func__, 6, "x data (%f - %f) (%f)\n", xPtr->data.F32[i], xPtr->data.F32[i+1], H);
@@ -609 +620 @@
     }
 
-    if ((x == NULL) || (PS_TYPE_F64 == x->type.type)) {
-        psFree(xPtr);
-    }
     if (PS_TYPE_F64 == y->type.type) {
         psFree(yPtr);
@@ -618 +626 @@
     return(spline);
 }
+
+void PS_POLY1D_PRINT(
+    psPolynomial1D *poly)
+{
+    printf("-------------- PS_POLY1D_PRINT() --------------\n");
+    printf("poly->nX is %d\n", poly->nX);
+    for (psS32 i = 0 ; i < (1 + poly->nX) ; i++) {
+        printf("poly->coeff[%d] is %f\n", i, poly->coeff[i]);
+    }
+}
+
+void PS_PRINT_SPLINE2(psSpline1D *mySpline)
+{
+    printf("-------------- PS_PRINT_SPLINE2() --------------\n");
+    printf("mySpline->n is %d\n", mySpline->n);
+    for (psS32 i = 0 ; i < mySpline->n ; i++) {
+        PS_POLY1D_PRINT(mySpline->spline[i]);
+    }
+    PS_VECTOR_PRINT_F32(mySpline->knots);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 
 
@@ -650 +713 @@
         //
         psLogMsg(__func__, PS_LOG_WARN,
-                 "psSpline1DEval(): x ordinate (%f) is outside the spline range (%f - %f).",
-                 x, spline->knots->data.F32[0],
-                 spline->knots->data.F32[n-1]);
+                 "psSpline1DEval(): x ordinate (%f) is outside the spline range (%f - %f) (%d).",
+                 x, spline->knots->data.F32[0], spline->knots->data.F32[n-1], n);
 
         if (x < spline->knots->data.F32[0]) {