*** empty log message ***

3 files changed, 774 insertions, 963 deletions

diff --git a/src/bltGrElemLine.C b/src/bltGrElemLine.C
index 9cd9337..358006e 100644
--- a/src/bltGrElemLine.C
+++ b/src/bltGrElemLine.C
@@ -35,7 +35,6 @@
 #include <X11/Xlib.h>
 
 #include "bltGraph.h"
-#include "bltSpline.h"
 #include "bltGrElemLine.h"
 #include "bltGrElemOption.h"
 #include "bltGrAxis.h"
@@ -1154,9 +1153,9 @@ void LineElement::generateSpline(MapInfo *mapPtr)
   niPts = count;
   int result = 0;
   if (smooth_ == CUBIC)
-    result = Blt_NaturalSpline(origPts, nOrigPts, iPts, niPts);
+    result = naturalSpline(origPts, nOrigPts, iPts, niPts);
   else if (smooth_ == QUADRATIC)
-    result = Blt_QuadraticSpline(origPts, nOrigPts, iPts, niPts);
+    result = quadraticSpline(origPts, nOrigPts, iPts, niPts);
 
   // The spline interpolation failed.  We will fall back to the current
   // coordinates and do no smoothing (standard line segments)
@@ -1199,7 +1198,7 @@ void LineElement::generateParametricSpline(MapInfo *mapPtr)
     p = origPts[i];
     q = origPts[j];
     count++;
-    if (LineRectClip(&exts, &p, &q)) {
+    if (lineRectClip(&exts, &p, &q)) {
       count += (int)(hypot(q.x - p.x, q.y - p.y) * 0.5);
     }
   }
@@ -1237,7 +1236,7 @@ void LineElement::generateParametricSpline(MapInfo *mapPtr)
     /* Is any part of the interval (line segment) in the plotting
      * area?  */
 
-    if (LineRectClip(&exts, &p, &q)) {
+    if (lineRectClip(&exts, &p, &q)) {
       double dp, dq;
 
       /* Distance of original point to p. */
@@ -1262,10 +1261,9 @@ void LineElement::generateParametricSpline(MapInfo *mapPtr)
   niPts = count;
   result = 0;
   if (smooth_ == CUBIC)
-    result = Blt_NaturalParametricSpline(origPts, nOrigPts, &exts, 0,
-					 iPts, niPts);
+    result = naturalParametricSpline(origPts, nOrigPts, &exts, 0, iPts, niPts);
   else if (smooth_ == CATROM)
-    result = Blt_CatromParametricSpline(origPts, nOrigPts, iPts, niPts);
+    result = catromParametricSpline(origPts, nOrigPts, iPts, niPts);
 
   // The spline interpolation failed.  We will fall back to the current
   // coordinates and do no smoothing (standard line segments)
@@ -1786,7 +1784,7 @@ void LineElement::mapErrorBars(LineStyle **styleMap)
 	  q = graphPtr_->map2D(low, y, ops->xAxis, ops->yAxis);
 	  segPtr->p = p;
 	  segPtr->q = q;
-	  if (LineRectClip(&exts, &segPtr->p, &segPtr->q)) {
+	  if (lineRectClip(&exts, &segPtr->p, &segPtr->q)) {
 	    segPtr++;
 	    *indexPtr++ = i;
 	  }
@@ -1794,7 +1792,7 @@ void LineElement::mapErrorBars(LineStyle **styleMap)
 	  segPtr->p.x = segPtr->q.x = p.x;
 	  segPtr->p.y = p.y - stylePtr->errorBarCapWidth;
 	  segPtr->q.y = p.y + stylePtr->errorBarCapWidth;
-	  if (LineRectClip(&exts, &segPtr->p, &segPtr->q)) {
+	  if (lineRectClip(&exts, &segPtr->p, &segPtr->q)) {
 	    segPtr++;
 	    *indexPtr++ = i;
 	  }
@@ -1802,7 +1800,7 @@ void LineElement::mapErrorBars(LineStyle **styleMap)
 	  segPtr->p.x = segPtr->q.x = q.x;
 	  segPtr->p.y = q.y - stylePtr->errorBarCapWidth;
 	  segPtr->q.y = q.y + stylePtr->errorBarCapWidth;
-	  if (LineRectClip(&exts, &segPtr->p, &segPtr->q)) {
+	  if (lineRectClip(&exts, &segPtr->p, &segPtr->q)) {
 	    segPtr++;
 	    *indexPtr++ = i;
 	  }
@@ -1855,7 +1853,7 @@ void LineElement::mapErrorBars(LineStyle **styleMap)
 	  q = graphPtr_->map2D(x, low, ops->xAxis, ops->yAxis);
 	  segPtr->p = p;
 	  segPtr->q = q;
-	  if (LineRectClip(&exts, &segPtr->p, &segPtr->q)) {
+	  if (lineRectClip(&exts, &segPtr->p, &segPtr->q)) {
 	    segPtr++;
 	    *indexPtr++ = i;
 	  }
@@ -1863,7 +1861,7 @@ void LineElement::mapErrorBars(LineStyle **styleMap)
 	  segPtr->p.y = segPtr->q.y = p.y;
 	  segPtr->p.x = p.x - stylePtr->errorBarCapWidth;
 	  segPtr->q.x = p.x + stylePtr->errorBarCapWidth;
-	  if (LineRectClip(&exts, &segPtr->p, &segPtr->q)) {
+	  if (lineRectClip(&exts, &segPtr->p, &segPtr->q)) {
 	    segPtr++;
 	    *indexPtr++ = i;
 	  }
@@ -1871,7 +1869,7 @@ void LineElement::mapErrorBars(LineStyle **styleMap)
 	  segPtr->p.y = segPtr->q.y = q.y;
 	  segPtr->p.x = q.x - stylePtr->errorBarCapWidth;
 	  segPtr->q.x = q.x + stylePtr->errorBarCapWidth;
-	  if (LineRectClip(&exts, &segPtr->p, &segPtr->q)) {
+	  if (lineRectClip(&exts, &segPtr->p, &segPtr->q)) {
 	    segPtr++;
 	    *indexPtr++ = i;
 	  }
diff --git a/src/bltGrElemLine.h b/src/bltGrElemLine.h
index 3cdd24f..30f5c50 100644
--- a/src/bltGrElemLine.h
+++ b/src/bltGrElemLine.h
@@ -32,6 +32,7 @@
 
 #include <tk.h>
 
+#include "bltGraph.h"
 #include "bltGrElem.h"
 #include "bltGrPenLine.h"
 
@@ -154,6 +155,11 @@ namespace Blt {
     int simplify(Point2d*, int, int, double, int*);
     double findSplit(Point2d*, int, int, int*);
 
+    int naturalSpline(Point2d*, int, Point2d*, int);
+    int quadraticSpline(Point2d*, int, Point2d*, int);
+    int naturalParametricSpline(Point2d*, int, Region2d*, int, Point2d*, int);
+    int catromParametricSpline(Point2d*, int, Point2d*, int);
+
   public:
     LineElement(Graph*, const char*, Tcl_HashEntry*);
     virtual ~LineElement();
diff --git a/src/bltGrElemLineSpline.C b/src/bltGrElemLineSpline.C
index 3d75f11..75957a2 100644
--- a/src/bltGrElemLineSpline.C
+++ b/src/bltGrElemLineSpline.C
@@ -33,22 +33,20 @@
 #include <stdlib.h>
 #include <string.h>
 
-#include "bltSpline.h"
+#include "bltGrElemLine.h"
 
 using namespace Blt;
 
 typedef double TriDiagonalMatrix[3];
 typedef struct {
-    double b, c, d;
+  double b, c, d;
 } Cubic2D;
 
 typedef struct {
-    double b, c, d, e, f;
+  double b, c, d, e, f;
 } Quint2D;
 
-/*
- * Quadratic spline parameters
- */
+// Quadratic spline parameters
 #define E1	param[0]
 #define E2	param[1]
 #define V1	param[2]
@@ -74,34 +72,24 @@ typedef struct {
  *
  *---------------------------------------------------------------------------
  */
-static int
-Search(
-    Point2d points[],		/* Contains the abscissas of the data
-				 * points of interpolation. */
-    int nPoints,		/* Dimension of x. */
-    double key,			/* Value whose relative position in
-				 * x is to be located. */
-    int *foundPtr)		/* (out) Returns 1 if s is found in
-				 * x and 0 otherwise. */
+static int Search(Point2d points[], int nPoints, double key, int *foundPtr)
 {
-    int high, low, mid;
-
-    low = 0;
-    high = nPoints - 1;
-
-    while (high >= low) {
-	mid = (high + low) / 2;
-	if (key > points[mid].x) {
-	    low = mid + 1;
-	} else if (key < points[mid].x) {
-	    high = mid - 1;
-	} else {
-	    *foundPtr = 1;
-	    return mid;
-	}
+  int low = 0;
+  int high = nPoints - 1;
+
+  while (high >= low) {
+    int mid = (high + low) / 2;
+    if (key > points[mid].x)
+      low = mid + 1;
+    else if (key < points[mid].x)
+      high = mid - 1;
+    else {
+      *foundPtr = 1;
+      return mid;
     }
-    *foundPtr = 0;
-    return low;
+  }
+  *foundPtr = 0;
+  return low;
 }
 
 /*
@@ -118,243 +106,184 @@ Search(
  *
  *---------------------------------------------------------------------------
  */
-static int
-QuadChoose(
-    Point2d *p,			/* Coordinates of one of the points of
-				 * interpolation */
-    Point2d *q,			/* Coordinates of one of the points of
-				 * interpolation */
-    double m1,			/* Derivative condition at point P */
-    double m2,			/* Derivative condition at point Q */
-    double epsilon)		/* Error tolerance used to distinguish
-				 * cases when m1 or m2 is relatively
-				 * close to the slope or twice the
-				 * slope of the line segment joining
-				 * the points P and Q.  If
-				 * epsilon is not 0.0, then epsilon
-				 * should be greater than or equal to
-				 * machine epsilon.  */
+static int QuadChoose(Point2d* p, Point2d* q, double m1, double m2, 
+		      double epsilon)
 {
-    double slope;
-
-    /* Calculate the slope of the line joining P and Q. */
-    slope = (q->y - p->y) / (q->x - p->x);
-
-    if (slope != 0.0) {
-	double relerr;
-	double mref, mref1, mref2, prod1, prod2;
-
-	prod1 = slope * m1;
-	prod2 = slope * m2;
-
-	/* Find the absolute values of the slopes slope, m1, and m2. */
-	mref = fabs(slope);
-	mref1 = fabs(m1);
-	mref2 = fabs(m2);
-
-	/*
-	 * If the relative deviation of m1 or m2 from slope is less than
-	 * epsilon, then choose case 2 or case 3.
-	 */
-	relerr = epsilon * mref;
-	if ((fabs(slope - m1) > relerr) && (fabs(slope - m2) > relerr) &&
-	    (prod1 >= 0.0) && (prod2 >= 0.0)) {
-	    double prod;
-
-	    prod = (mref - mref1) * (mref - mref2);
-	    if (prod < 0.0) {
-		/*
-		 * l1, the line through (x1,y1) with slope m1, and l2,
-		 * the line through (x2,y2) with slope m2, intersect
-		 * at a point whose abscissa is between x1 and x2.
-		 * The abscissa becomes a knot of the spline.
-		 */
-		return 1;
-	    }
-	    if (mref1 > (mref * 2.0)) {
-		if (mref2 <= ((2.0 - epsilon) * mref)) {
-		    return 3;
-		}
-	    } else if (mref2 <= (mref * 2.0)) {
-		/*
-		 * Both l1 and l2 cross the line through
-		 * (x1+x2)/2.0,y1 and (x1+x2)/2.0,y2, which is the
-		 * midline of the rectangle formed by P and Q or both
-		 * m1 and m2 have signs different than the sign of
-		 * slope, or one of m1 and m2 has opposite sign from
-		 * slope and l1 and l2 intersect to the left of x1 or
-		 * to the right of x2.  The point (x1+x2)/2. is a knot
-		 * of the spline.
-		 */
-		return 2;
-	    } else if (mref1 <= ((2.0 - epsilon) * mref)) {
-		/*
-		 * In cases 3 and 4, sign(m1)=sign(m2)=sign(slope).
-		 * Either l1 or l2 crosses the midline, but not both.
-		 * Choose case 4 if mref1 is greater than
-		 * (2.-epsilon)*mref; otherwise, choose case 3.
-		 */
-		return 3;
-	    }
-	    /*
-	     * If neither l1 nor l2 crosses the midline, the spline
-	     * requires two knots between x1 and x2.
-	     */
-	    return 4;
-	} else {
-	    /*
-	     * The sign of at least one of the slopes m1 or m2 does not
-	     * agree with the sign of *slope*.
-	     */
-	    if ((prod1 < 0.0) && (prod2 < 0.0)) {
-		return 2;
-	    } else if (prod1 < 0.0) {
-		if (mref2 > ((epsilon + 1.0) * mref)) {
-		    return 1;
-		} else {
-		    return 2;
-		}
-	    } else if (mref1 > ((epsilon + 1.0) * mref)) {
-		return 1;
-	    } else {
-		return 2;
-	    }
-	}
-    } else if ((m1 * m2) >= 0.0) {
+  // Calculate the slope of the line joining P and Q
+  double slope = (q->y - p->y) / (q->x - p->x);
+
+  if (slope != 0.0) {
+    double prod1 = slope * m1;
+    double prod2 = slope * m2;
+
+    // Find the absolute values of the slopes slope, m1, and m2
+    double mref = fabs(slope);
+    double mref1 = fabs(m1);
+    double mref2 = fabs(m2);
+
+    // If the relative deviation of m1 or m2 from slope is less than
+    // epsilon, then choose case 2 or case 3.
+    double relerr = epsilon * mref;
+    if ((fabs(slope - m1) > relerr) && (fabs(slope - m2) > relerr) &&
+	(prod1 >= 0.0) && (prod2 >= 0.0)) {
+      double prod = (mref - mref1) * (mref - mref2);
+      if (prod < 0.0) {
+	// l1, the line through (x1,y1) with slope m1, and l2,
+	// the line through (x2,y2) with slope m2, intersect
+	// at a point whose abscissa is between x1 and x2.
+	// The abscissa becomes a knot of the spline.
+	return 1;
+      }
+      if (mref1 > (mref * 2.0)) {
+	if (mref2 <= ((2.0 - epsilon) * mref))
+	  return 3;
+      }
+      else if (mref2 <= (mref * 2.0)) {
+	// Both l1 and l2 cross the line through
+	// (x1+x2)/2.0,y1 and (x1+x2)/2.0,y2, which is the
+	// midline of the rectangle formed by P and Q or both
+	// m1 and m2 have signs different than the sign of
+	// slope, or one of m1 and m2 has opposite sign from
+	// slope and l1 and l2 intersect to the left of x1 or
+	// to the right of x2.  The point (x1+x2)/2. is a knot
+	// of the spline.
 	return 2;
-    } else {
+      }
+      else if (mref1 <= ((2.0 - epsilon) * mref)) {
+	// In cases 3 and 4, sign(m1)=sign(m2)=sign(slope).
+	// Either l1 or l2 crosses the midline, but not both.
+	// Choose case 4 if mref1 is greater than
+	// (2.-epsilon)*mref; otherwise, choose case 3.
+	return 3;
+      }
+      // If neither l1 nor l2 crosses the midline, the spline
+      // requires two knots between x1 and x2.
+      return 4;
+    }
+    else {
+      // The sign of at least one of the slopes m1 or m2 does not
+      // agree with the sign of *slope*.
+      if ((prod1 < 0.0) && (prod2 < 0.0)) {
+	return 2;
+      }
+      else if (prod1 < 0.0) {
+	if (mref2 > ((epsilon + 1.0) * mref))
+	  return 1;
+	else
+	  return 2;
+      }
+      else if (mref1 > ((epsilon + 1.0) * mref))
 	return 1;
+      else
+	return 2;
     }
+  }
+  else if ((m1 * m2) >= 0.0)
+    return 2;
+  else
+    return 1;
 }
 
 /*
  *---------------------------------------------------------------------------
- *
- * QuadCases --
- *
  *	Computes the knots and other parameters of the spline on the
  *	interval PQ.
- *
- *
  * On input--
- *
  *	P and Q are the coordinates of the points of interpolation.
- *
  *	m1 is the slope at P.
- *
  *	m2 is the slope at Q.
- *
  *	ncase controls the number and location of the knots.
- *
- *
  * On output--
  *
  *	(v1,v2),(w1,w2),(z1,z2), and (e1,e2) are the coordinates of
  *	the knots and other parameters of the spline on P.
  *	(e1,e2) and Q are used only if ncase=4.
- *
  *---------------------------------------------------------------------------
  */
-static void
-QuadCases(Point2d *p, Point2d *q, double m1, double m2, double param[], 
-	  int which)
+static void QuadCases(Point2d* p, Point2d* q, double m1, double m2, 
+		      double param[], int which)
 {
-    if ((which == 3) || (which == 4)) {	/* Parameters used in both 3 and 4 */
-	double mbar1, mbar2, mbar3, c1, d1, h1, j1, k1;
-
-	c1 = p->x + (q->y - p->y) / m1;
-	d1 = q->x + (p->y - q->y) / m2;
-	h1 = c1 * 2.0 - p->x;
-	j1 = d1 * 2.0 - q->x;
-	mbar1 = (q->y - p->y) / (h1 - p->x);
-	mbar2 = (p->y - q->y) / (j1 - q->x);
-
-	if (which == 4) {	/* Case 4. */
-	    Y1 = (p->x + c1) / 2.0;
-	    V1 = (p->x + Y1) / 2.0;
-	    V2 = m1 * (V1 - p->x) + p->y;
-	    Z1 = (d1 + q->x) / 2.0;
-	    W1 = (q->x + Z1) / 2.0;
-	    W2 = m2 * (W1 - q->x) + q->y;
-	    mbar3 = (W2 - V2) / (W1 - V1);
-	    Y2 = mbar3 * (Y1 - V1) + V2;
-	    Z2 = mbar3 * (Z1 - V1) + V2;
-	    E1 = (Y1 + Z1) / 2.0;
-	    E2 = mbar3 * (E1 - V1) + V2;
-	} else {		/* Case 3. */
-	    k1 = (p->y - q->y + q->x * mbar2 - p->x * mbar1) / (mbar2 - mbar1);
-	    if (fabs(m1) > fabs(m2)) {
-		Z1 = (k1 + p->x) / 2.0;
-	    } else {
-		Z1 = (k1 + q->x) / 2.0;
-	    }
-	    V1 = (p->x + Z1) / 2.0;
-	    V2 = p->y + m1 * (V1 - p->x);
-	    W1 = (q->x + Z1) / 2.0;
-	    W2 = q->y + m2 * (W1 - q->x);
-	    Z2 = V2 + (W2 - V2) / (W1 - V1) * (Z1 - V1);
-	}
-    } else if (which == 2) {	/* Case 2. */
-	Z1 = (p->x + q->x) / 2.0;
-	V1 = (p->x + Z1) / 2.0;
-	V2 = p->y + m1 * (V1 - p->x);
-	W1 = (Z1 + q->x) / 2.0;
-	W2 = q->y + m2 * (W1 - q->x);
-	Z2 = (V2 + W2) / 2.0;
-    } else {			/* Case 1. */
-	double ztwo;
-
-	Z1 = (p->y - q->y + m2 * q->x - m1 * p->x) / (m2 - m1);
-	ztwo = p->y + m1 * (Z1 - p->x);
-	V1 = (p->x + Z1) / 2.0;
-	V2 = (p->y + ztwo) / 2.0;
-	W1 = (Z1 + q->x) / 2.0;
-	W2 = (ztwo + q->y) / 2.0;
-	Z2 = V2 + (W2 - V2) / (W1 - V1) * (Z1 - V1);
+  if ((which == 3) || (which == 4)) {
+    double c1 = p->x + (q->y - p->y) / m1;
+    double d1 = q->x + (p->y - q->y) / m2;
+    double h1 = c1 * 2.0 - p->x;
+    double j1 = d1 * 2.0 - q->x;
+    double mbar1 = (q->y - p->y) / (h1 - p->x);
+    double mbar2 = (p->y - q->y) / (j1 - q->x);
+
+    if (which == 4) {
+      // Case 4
+      Y1 = (p->x + c1) / 2.0;
+      V1 = (p->x + Y1) / 2.0;
+      V2 = m1 * (V1 - p->x) + p->y;
+      Z1 = (d1 + q->x) / 2.0;
+      W1 = (q->x + Z1) / 2.0;
+      W2 = m2 * (W1 - q->x) + q->y;
+      double mbar3 = (W2 - V2) / (W1 - V1);
+      Y2 = mbar3 * (Y1 - V1) + V2;
+      Z2 = mbar3 * (Z1 - V1) + V2;
+      E1 = (Y1 + Z1) / 2.0;
+      E2 = mbar3 * (E1 - V1) + V2;
     }
+    else {
+      // Case 3
+      double k1 = (p->y - q->y + q->x * mbar2 - p->x * mbar1) / (mbar2 - mbar1);
+      if (fabs(m1) > fabs(m2)) {
+	Z1 = (k1 + p->x) / 2.0;
+      } else {
+	Z1 = (k1 + q->x) / 2.0;
+      }
+      V1 = (p->x + Z1) / 2.0;
+      V2 = p->y + m1 * (V1 - p->x);
+      W1 = (q->x + Z1) / 2.0;
+      W2 = q->y + m2 * (W1 - q->x);
+      Z2 = V2 + (W2 - V2) / (W1 - V1) * (Z1 - V1);
+    }
+  }
+  else if (which == 2) {
+    // Case 2
+    Z1 = (p->x + q->x) / 2.0;
+    V1 = (p->x + Z1) / 2.0;
+    V2 = p->y + m1 * (V1 - p->x);
+    W1 = (Z1 + q->x) / 2.0;
+    W2 = q->y + m2 * (W1 - q->x);
+    Z2 = (V2 + W2) / 2.0;
+  }
+  else {
+    // Case 1
+    Z1 = (p->y - q->y + m2 * q->x - m1 * p->x) / (m2 - m1);
+    double ztwo = p->y + m1 * (Z1 - p->x);
+    V1 = (p->x + Z1) / 2.0;
+    V2 = (p->y + ztwo) / 2.0;
+    W1 = (Z1 + q->x) / 2.0;
+    W2 = (ztwo + q->y) / 2.0;
+    Z2 = V2 + (W2 - V2) / (W1 - V1) * (Z1 - V1);
+  }
 }
 
-static int
-QuadSelect(Point2d *p, Point2d *q, double m1, double m2, double epsilon,
-	   double param[])
+static int QuadSelect(Point2d* p, Point2d* q, double m1, double m2, 
+		      double epsilon, double param[])
 {
-    int ncase;
-
-    ncase = QuadChoose(p, q, m1, m2, epsilon);
-    QuadCases(p, q, m1, m2, param, ncase);
-    return ncase;
+  int ncase = QuadChoose(p, q, m1, m2, epsilon);
+  QuadCases(p, q, m1, m2, param, ncase);
+  return ncase;
 }
 
-/*
- *---------------------------------------------------------------------------
- *
- * QuadGetImage --
- *
- *---------------------------------------------------------------------------
- */
-INLINE static double
-QuadGetImage(double p1, double p2, double p3, double x1, double x2, double x3)
+static double QuadGetImage(double p1, double p2, double p3, double x1, 
+			   double x2, double x3)
 {
-    double A, B, C;
-    double y;
-
-    A = x1 - x2;
-    B = x2 - x3;
-    C = x1 - x3;
+  double A = x1 - x2;
+  double B = x2 - x3;
+  double C = x1 - x3;
 
-    y = (p1 * (A * A) + p2 * 2.0 * B * A + p3 * (B * B)) / (C * C);
-    return y;
+  double y = (p1 * (A * A) + p2 * 2.0 * B * A + p3 * (B * B)) / (C * C);
+  return y;
 }
 
 /*
  *---------------------------------------------------------------------------
- *
- * QuadSpline --
- *
  *	Finds the image of a point in x.
- *
  *	On input
- *
  *	x	Contains the value at which the spline is evaluated.
  *	leftX, leftY
  *		Coordinates of the left-hand data point used in the
@@ -367,161 +296,115 @@ QuadGetImage(double p1, double p2, double p3, double x1, double x2, double x3)
  *	ncase	Controls the evaluation of the spline by indicating
  *		whether one or two knots were placed in the interval
  *		(xtabs,xtabs1).
- *
- * Results:
- *	The image of the spline at x.
- *
  *---------------------------------------------------------------------------
  */
-static void
-QuadSpline(
-    Point2d *intp,		/* Value at which spline is evaluated */
-    Point2d *left,		/* Point to the left of the data point to
-				 * be evaluated */
-    Point2d *right,		/* Point to the right of the data point to
-				 * be evaluated */
-    double param[],		/* Parameters of the spline */
-    int ncase)			/* Controls the evaluation of the
-				 * spline by indicating whether one or
-				 * two knots were placed in the
-				 * interval (leftX,rightX) */
-{
-    double y;
-
-    if (ncase == 4) {
-	/*
-	 * Case 4:  More than one knot was placed in the interval.
-	 */
-
-	/*
-	 * Determine the location of data point relative to the 1st knot.
-	 */
-	if (Y1 > intp->x) {
-	    y = QuadGetImage(left->y, V2, Y2, Y1, intp->x, left->x);
-	} else if (Y1 < intp->x) {
-	    /*
-	     * Determine the location of the data point relative to
-	     * the 2nd knot.
-	     */
-	    if (Z1 > intp->x) {
-		y = QuadGetImage(Y2, E2, Z2, Z1, intp->x, Y1);
-	    } else if (Z1 < intp->x) {
-		y = QuadGetImage(Z2, W2, right->y, right->x, intp->x, Z1);
-	    } else {
-		y = Z2;
-	    }
-	} else {
-	    y = Y2;
-	}
-    } else {
+static void QuadSpline(Point2d* intp, Point2d* left, Point2d* right,
+		       double param[], int ncase)
 
-	/*
-	 * Cases 1, 2, or 3:
-	 *
-	 * Determine the location of the data point relative to the
-	 * knot.
-	 */
-	if (Z1 < intp->x) {
-	    y = QuadGetImage(Z2, W2, right->y, right->x, intp->x, Z1);
-	} else if (Z1 > intp->x) {
-	    y = QuadGetImage(left->y, V2, Z2, Z1, intp->x, left->x);
-	} else {
-	    y = Z2;
-	}
+{
+  double y;
+
+  if (ncase == 4) {
+    // Case 4:  More than one knot was placed in the interval.
+    // Determine the location of data point relative to the 1st knot.
+    if (Y1 > intp->x)
+      y = QuadGetImage(left->y, V2, Y2, Y1, intp->x, left->x);
+    else if (Y1 < intp->x) {
+      // Determine the location of the data point relative to the 2nd knot.
+      if (Z1 > intp->x)
+	y = QuadGetImage(Y2, E2, Z2, Z1, intp->x, Y1);
+      else if (Z1 < intp->x)
+	y = QuadGetImage(Z2, W2, right->y, right->x, intp->x, Z1);
+      else
+	y = Z2;
     }
-    intp->y = y;
+    else
+      y = Y2;
+  }
+  else {
+    // Cases 1, 2, or 3:
+    // Determine the location of the data point relative to the knot.
+    if (Z1 < intp->x)
+      y = QuadGetImage(Z2, W2, right->y, right->x, intp->x, Z1);
+    else if (Z1 > intp->x)
+      y = QuadGetImage(left->y, V2, Z2, Z1, intp->x, left->x);
+    else
+      y = Z2;
+  }
+
+  intp->y = y;
 }
 
 /*
  *---------------------------------------------------------------------------
- *
- * QuadSlopes --
- *
  * 	Calculates the derivative at each of the data points.  The
  * 	slopes computed will insure that an osculatory quadratic
  * 	spline will have one additional knot between two adjacent
  * 	points of interpolation.  Convexity and monotonicity are
  * 	preserved wherever these conditions are compatible with the
  * 	data.
- *
- * Results:
- *	The output array "m" is filled with the derivates at each
- *	data point.
- *
  *---------------------------------------------------------------------------
  */
-static void
-QuadSlopes(Point2d *points, double *m, int nPoints)
+static void QuadSlopes(Point2d *points, double *m, int nPoints)
 {
-    double xbar, xmid, xhat, ydif1, ydif2;
-    double yxmid;
-    double m1, m2;
-    double m1s, m2s;
-    int i, n, l;
-
-    m1s = m2s = m1 = m2 = 0;
-    for (l = 0, i = 1, n = 2; i < (nPoints - 1); l++, i++, n++) {
-	/*
-	 * Calculate the slopes of the two lines joining three
-	 * consecutive data points.
-	 */
-	ydif1 = points[i].y - points[l].y;
-	ydif2 = points[n].y - points[i].y;
-	m1 = ydif1 / (points[i].x - points[l].x);
-	m2 = ydif2 / (points[n].x - points[i].x);
-	if (i == 1) {
-	    m1s = m1, m2s = m2;	/* Save slopes of starting point */
-	}
-	/*
-	 * If one of the preceding slopes is zero or if they have opposite
-	 * sign, assign the value zero to the derivative at the middle
-	 * point.
-	 */
-	if ((m1 == 0.0) || (m2 == 0.0) || ((m1 * m2) <= 0.0)) {
-	    m[i] = 0.0;
-	} else if (fabs(m1) > fabs(m2)) {
-	    /*
-	     * Calculate the slope by extending the line with slope m1.
-	     */
-	    xbar = ydif2 / m1 + points[i].x;
-	    xhat = (xbar + points[n].x) / 2.0;
-	    m[i] = ydif2 / (xhat - points[i].x);
-	} else {
-	    /*
-	     * Calculate the slope by extending the line with slope m2.
-	     */
-	    xbar = -ydif1 / m2 + points[i].x;
-	    xhat = (points[l].x + xbar) / 2.0;
-	    m[i] = ydif1 / (points[i].x - xhat);
-	}
+  double m1s =0;
+  double m2s =0;
+  double m1 =0;
+  double m2 =0;
+  int i, n, l;
+  for (l = 0, i = 1, n = 2; i < (nPoints - 1); l++, i++, n++) {
+    // Calculate the slopes of the two lines joining three
+    // consecutive data points.
+    double ydif1 = points[i].y - points[l].y;
+    double ydif2 = points[n].y - points[i].y;
+    m1 = ydif1 / (points[i].x - points[l].x);
+    m2 = ydif2 / (points[n].x - points[i].x);
+    if (i == 1) {
+      // Save slopes of starting point
+      m1s = m1;
+      m2s = m2;
     }
-
-    /* Calculate the slope at the last point, x(n). */
-    i = nPoints - 2;
-    n = nPoints - 1;
-    if ((m1 * m2) < 0.0) {
-	m[n] = m2 * 2.0;
-    } else {
-	xmid = (points[i].x + points[n].x) / 2.0;
-	yxmid = m[i] * (xmid - points[i].x) + points[i].y;
-	m[n] = (points[n].y - yxmid) / (points[n].x - xmid);
-	if ((m[n] * m2) < 0.0) {
-	    m[n] = 0.0;
-	}
+    // If one of the preceding slopes is zero or if they have opposite
+    // sign, assign the value zero to the derivative at the middle point.
+    if ((m1 == 0.0) || (m2 == 0.0) || ((m1 * m2) <= 0.0))
+      m[i] = 0.0;
+    else if (fabs(m1) > fabs(m2)) {
+      // Calculate the slope by extending the line with slope m1.
+      double xbar = ydif2 / m1 + points[i].x;
+      double xhat = (xbar + points[n].x) / 2.0;
+      m[i] = ydif2 / (xhat - points[i].x);
     }
-
-    /* Calculate the slope at the first point, x(0). */
-    if ((m1s * m2s) < 0.0) {
-	m[0] = m1s * 2.0;
-    } else {
-	xmid = (points[0].x + points[1].x) / 2.0;
-	yxmid = m[1] * (xmid - points[1].x) + points[1].y;
-	m[0] = (yxmid - points[0].y) / (xmid - points[0].x);
-	if ((m[0] * m1s) < 0.0) {
-	    m[0] = 0.0;
-	}
+    else {
+      // Calculate the slope by extending the line with slope m2.
+      double xbar = -ydif1 / m2 + points[i].x;
+      double xhat = (points[l].x + xbar) / 2.0;
+      m[i] = ydif1 / (points[i].x - xhat);
     }
-
+  }
+
+  // Calculate the slope at the last point, x(n).
+  i = nPoints - 2;
+  n = nPoints - 1;
+  if ((m1 * m2) < 0.0)
+    m[n] = m2 * 2.0;
+  else {
+    double xmid = (points[i].x + points[n].x) / 2.0;
+    double yxmid = m[i] * (xmid - points[i].x) + points[i].y;
+    m[n] = (points[n].y - yxmid) / (points[n].x - xmid);
+    if ((m[n] * m2) < 0.0)
+      m[n] = 0.0;
+  }
+
+  // Calculate the slope at the first point, x(0).
+  if ((m1s * m2s) < 0.0)
+    m[0] = m1s * 2.0;
+  else {
+    double xmid = (points[0].x + points[1].x) / 2.0;
+    double yxmid = m[1] * (xmid - points[1].x) + points[1].y;
+    m[0] = (yxmid - points[0].y) / (xmid - points[0].x);
+    if ((m[0] * m1s) < 0.0)
+      m[0] = 0.0;
+  }
 }
 
 /*
@@ -574,179 +457,149 @@ QuadSlopes(Point2d *points, double *m, int nPoints)
  *      QuadSpline
  *---------------------------------------------------------------------------
  */
-static int
-QuadEval(
-    Point2d origPts[],
-    int nOrigPts,
-    Point2d intpPts[],
-    int nIntpPts,
-    double *m,			/* Slope of the spline at each point
-				 * of interpolation. */
-    double epsilon)		/* Relative error tolerance (see choose) */
+static int QuadEval(Point2d origPts[], int nOrigPts, Point2d intpPts[],
+		    int nIntpPts, double *m, double epsilon)
 {
-    int error;
-    int i, j;
-    double param[10];
-    int ncase;
-    int start, end;
-    int l, p;
-    int n;
-    int found;
-
-    /* Initialize indices and set error result */
-    error = 0;
-    l = nOrigPts - 1;
-    p = l - 1;
-    ncase = 1;
-
-    /*
-     * Determine if abscissas of new vector are non-decreasing.
-     */
-    for (j = 1; j < nIntpPts; j++) {
-	if (intpPts[j].x < intpPts[j - 1].x) {
-	    return 2;
-	}
-    }
-    /*
-     * Determine if any of the points in xval are LESS than the
-     * abscissa of the first data point.
-     */
-    for (start = 0; start < nIntpPts; start++) {
-	if (intpPts[start].x >= origPts[0].x) {
-	    break;
-	}
-    }
-    /*
-     * Determine if any of the points in xval are GREATER than the
-     * abscissa of the l data point.
-     */
-    for (end = nIntpPts - 1; end >= 0; end--) {
-	if (intpPts[end].x <= origPts[l].x) {
-	    break;
-	}
-    }
-
-    if (start > 0) {
-	error = 1;		/* Set error value to indicate that
-				 * extrapolation has occurred. */
-	/*
-	 * Calculate the images of points of evaluation whose abscissas
-	 * are less than the abscissa of the first data point.
-	 */
-	ncase = QuadSelect(origPts, origPts + 1, m[0], m[1], epsilon, param);
-	for (j = 0; j < (start - 1); j++) {
-	    QuadSpline(intpPts + j, origPts, origPts + 1, param, ncase);
-	}
-	if (nIntpPts == 1) {
-	    return error;
-	}
-    }
-    if ((nIntpPts == 1) && (end != (nIntpPts - 1))) {
-	goto noExtrapolation;
-    }
+  double param[10];
+
+  // Initialize indices and set error result
+  int error = 0;
+  int l = nOrigPts - 1;
+  int p = l - 1;
+  int ncase = 1;
+
+  // Determine if abscissas of new vector are non-decreasing.
+  for (int jj=1; jj<nIntpPts; jj++) {
+    if (intpPts[jj].x < intpPts[jj - 1].x)
+      return 2;
+  }
+  // Determine if any of the points in xval are LESS than the
+  // abscissa of the first data point.
+  int start;
+  for (start = 0; start < nIntpPts; start++) {
+    if (intpPts[start].x >= origPts[0].x)
+      break;
+  }
+  // Determine if any of the points in xval are GREATER than the
+  // abscissa of the l data point.
+  int end;
+  for (end = nIntpPts - 1; end >= 0; end--) {
+    if (intpPts[end].x <= origPts[l].x)
+      break;
+  }
+
+  if (start > 0) {
+    // Set error value to indicate that extrapolation has occurred
+    error = 1;
+
+    // Calculate the images of points of evaluation whose abscissas
+    // are less than the abscissa of the first data point.
+    ncase = QuadSelect(origPts, origPts + 1, m[0], m[1], epsilon, param);
+    for (int jj=0; jj<(start - 1); jj++)
+      QuadSpline(intpPts + jj, origPts, origPts + 1, param, ncase);
+    if (nIntpPts == 1)
+      return error;
+  }
+  int ii;
+  int nn;
+  if ((nIntpPts == 1) && (end != (nIntpPts - 1)))
+    goto noExtrapolation;
     
-    /*
-     * Search locates the interval in which the first in-range
-     * point of evaluation lies.
-     */
-
-    i = Search(origPts, nOrigPts, intpPts[start].x, &found);
+  // Search locates the interval in which the first in-range
+  // point of evaluation lies.
+  int found;
+  ii = Search(origPts, nOrigPts, intpPts[start].x, &found);
     
-    n = i + 1;
-    if (n >= nOrigPts) {
-	n = nOrigPts - 1;
-	i = nOrigPts - 2;
-    }
-    /*
-     * If the first in-range point of evaluation is equal to one
-     * of the data points, assign the appropriate value from y.
-     * Continue until a point of evaluation is found which is not
-     * equal to a data point.
-     */
-    if (found) {
-	do {
-	    intpPts[start].y = origPts[i].y;
-	    start++;
-	    if (start >= nIntpPts) {
-		return error;
-	    }
-	} while (intpPts[start - 1].x == intpPts[start].x);
+  nn = ii + 1;
+  if (nn >= nOrigPts) {
+    nn = nOrigPts - 1;
+    ii = nOrigPts - 2;
+  }
+  /*
+   * If the first in-range point of evaluation is equal to one
+   * of the data points, assign the appropriate value from y.
+   * Continue until a point of evaluation is found which is not
+   * equal to a data point.
+   */
+  if (found) {
+    do {
+      intpPts[start].y = origPts[ii].y;
+      start++;
+      if (start >= nIntpPts) {
+	return error;
+      }
+    } while (intpPts[start - 1].x == intpPts[start].x);
 	
-	for (;;) {
-	    if (intpPts[start].x < origPts[n].x) {
-		break;	/* Break out of for-loop */
-	    }
-	    if (intpPts[start].x == origPts[n].x) {
-		do {
-		    intpPts[start].y = origPts[n].y;
-		    start++;
-		    if (start >= nIntpPts) {
-			return error;
-		    }
-		} while (intpPts[start].x == intpPts[start - 1].x);
-	    }
-	    i++;
-	    n++;
-	}
+    for (;;) {
+      if (intpPts[start].x < origPts[nn].x) {
+	break;	/* Break out of for-loop */
+      }
+      if (intpPts[start].x == origPts[nn].x) {
+	do {
+	  intpPts[start].y = origPts[nn].y;
+	  start++;
+	  if (start >= nIntpPts) {
+	    return error;
+	  }
+	} while (intpPts[start].x == intpPts[start - 1].x);
+      }
+      ii++;
+      nn++;
     }
-    /*
-     * Calculate the images of all the points which lie within
-     * range of the data.
-     */
-    if ((i > 0) || (error != 1)) {
-	ncase = QuadSelect(origPts + i, origPts + n, m[i], m[n], 
-			   epsilon, param);
+  }
+  /*
+   * Calculate the images of all the points which lie within
+   * range of the data.
+   */
+  if ((ii > 0) || (error != 1))
+    ncase = QuadSelect(origPts+ii, origPts+nn, m[ii], m[nn], epsilon, param);
+
+  for (int jj=start; jj<=end; jj++) {
+    // If xx(j) - x(n) is negative, do not recalculate
+    // the parameters for this section of the spline since
+    // they are already known.
+    if (intpPts[jj].x == origPts[nn].x) {
+      intpPts[jj].y = origPts[nn].y;
+      continue;
     }
-    for (j = start; j <= end; j++) {
-	/*
-	 * If xx(j) - x(n) is negative, do not recalculate
-	 * the parameters for this section of the spline since
-	 * they are already known.
-	 */
-	if (intpPts[j].x == origPts[n].x) {
-	    intpPts[j].y = origPts[n].y;
-	    continue;
-	} else if (intpPts[j].x > origPts[n].x) {
-	    double delta;
+    else if (intpPts[jj].x > origPts[nn].x) {
+      double delta;
 	    
-	    /* Determine that the routine is in the correct part of
-	       the spline. */
-	    do {
-		i++, n++;
-		delta = intpPts[j].x - origPts[n].x;
-	    } while (delta > 0.0);
+      // Determine that the routine is in the correct part of the spline
+      do {
+	ii++;
+	nn++;
+	delta = intpPts[jj].x - origPts[nn].x;
+      } while (delta > 0.0);
 	    
-	    if (delta < 0.0) {
-		ncase = QuadSelect(origPts + i, origPts + n, m[i], 
-			   m[n], epsilon, param);
-	    } else if (delta == 0.0) {
-		intpPts[j].y = origPts[n].y;
-		continue;
-	    }
-	}
-	QuadSpline(intpPts + j, origPts + i, origPts + n, param, ncase);
+      if (delta < 0.0)
+	ncase = QuadSelect(origPts+ii, origPts+nn, m[ii], m[nn], 
+			   epsilon, param);
+      else if (delta == 0.0) {
+	intpPts[jj].y = origPts[nn].y;
+	continue;
+      }
     }
+    QuadSpline(intpPts+jj, origPts+ii, origPts+nn, param, ncase);
+  }
     
-    if (end == (nIntpPts - 1)) {
-	return error;
-    }
-    if ((n == l) && (intpPts[end].x != origPts[l].x)) {
-	goto noExtrapolation;
-    }
+  if (end == (nIntpPts - 1))
+    return error;
 
-    error = 1;			/* Set error value to indicate that
-				 * extrapolation has occurred. */
-    ncase = QuadSelect(origPts + p, origPts + l, m[p], m[l], epsilon, param);
+  if ((nn == l) && (intpPts[end].x != origPts[l].x))
+    goto noExtrapolation;
 
-  noExtrapolation:
-    /*
-     * Calculate the images of the points of evaluation whose
-     * abscissas are greater than the abscissa of the last data point.
-     */
-    for (j = (end + 1); j < nIntpPts; j++) {
-	QuadSpline(intpPts + j, origPts + p, origPts + l, param, ncase);
-    }
-    return error;
+  // Set error value to indicate that extrapolation has occurred
+  error = 1;
+  ncase = QuadSelect(origPts + p, origPts + l, m[p], m[l], epsilon, param);
+
+ noExtrapolation:
+  // Calculate the images of the points of evaluation whose
+  // abscissas are greater than the abscissa of the last data point.
+  for (int jj=(end + 1); jj<nIntpPts; jj++)
+    QuadSpline(intpPts + jj, origPts + p, origPts + l, param, ncase);
+
+  return error;
 }
 
 /*
@@ -771,116 +624,108 @@ QuadEval(
  *   work   Contains the value of the first derivative at each data point
  *   ym     Contains the interpolated spline value at each data point
  */
-int
-Blt_QuadraticSpline(Point2d *origPts, int nOrigPts, Point2d *intpPts, 
-		    int nIntpPts)
+int LineElement::quadraticSpline(Point2d *origPts, int nOrigPts, 
+				 Point2d *intpPts, int nIntpPts)
 {
-    double* work = (double*)malloc(nOrigPts * sizeof(double));
-    double epsilon = 0.0;
-    /* allocate space for vectors used in calculation */
-    QuadSlopes(origPts, work, nOrigPts);
-    int result = QuadEval(origPts, nOrigPts, intpPts, nIntpPts, work, epsilon);
-    free(work);
-    if (result > 1) {
-	return 0;
-    }
-    return 1;
+  double* work = (double*)malloc(nOrigPts * sizeof(double));
+  double epsilon = 0.0;
+  /* allocate space for vectors used in calculation */
+  QuadSlopes(origPts, work, nOrigPts);
+  int result = QuadEval(origPts, nOrigPts, intpPts, nIntpPts, work, epsilon);
+  free(work);
+  if (result > 1) {
+    return 0;
+  }
+  return 1;
 }
 
 /*
  *---------------------------------------------------------------------------
- *
  * Reference:
  *	Numerical Analysis, R. Burden, J. Faires and A. Reynolds.
  *	Prindle, Weber & Schmidt 1981 pp 112
- *
- * Parameters:
- *	origPts - vector of points, assumed to be sorted along x.
- *	intpPts - vector of new points.
- *
  *---------------------------------------------------------------------------
  */
-int
-Blt_NaturalSpline(Point2d *origPts, int nOrigPts, Point2d *intpPts, 
-		  int nIntpPts)
+int LineElement::naturalSpline(Point2d *origPts, int nOrigPts, 
+			       Point2d *intpPts, int nIntpPts)
 {
-    Point2d *ip, *iend;
-    double x, dy, alpha;
-    int isKnot;
-    int i, j, n;
-
-    double* dx = (double*)malloc(sizeof(double) * nOrigPts);
-    /* Calculate vector of differences */
-    for (i = 0, j = 1; j < nOrigPts; i++, j++) {
-	dx[i] = origPts[j].x - origPts[i].x;
-	if (dx[i] < 0.0) {
-	    return 0;
-	}
-    }
-    n = nOrigPts - 1;		/* Number of intervals. */
-    TriDiagonalMatrix* A = 
-      (TriDiagonalMatrix*)malloc(sizeof(TriDiagonalMatrix) * nOrigPts);
-    if (!A) {
-	free(dx);
-	return 0;
-    }
-    /* Vectors to solve the tridiagonal matrix */
-    A[0][0] = A[n][0] = 1.0;
-    A[0][1] = A[n][1] = 0.0;
-    A[0][2] = A[n][2] = 0.0;
-
-    /* Calculate the intermediate results */
-    for (i = 0, j = 1; j < n; j++, i++) {
-	alpha = 3.0 * ((origPts[j + 1].y / dx[j]) - (origPts[j].y / dx[i]) - 
-		       (origPts[j].y / dx[j]) + (origPts[i].y / dx[i]));
-	A[j][0] = 2 * (dx[j] + dx[i]) - dx[i] * A[i][1];
-	A[j][1] = dx[j] / A[j][0];
-	A[j][2] = (alpha - dx[i] * A[i][2]) / A[j][0];
-    }
-
-    Cubic2D* eq = (Cubic2D*)malloc(sizeof(Cubic2D) * nOrigPts);
-    if (!eq) {
-	free(A);
-	free(dx);
-	return 0;
-    }
-    eq[0].c = eq[n].c = 0.0;
-    for (j = n, i = n - 1; i >= 0; i--, j--) {
-	eq[i].c = A[i][2] - A[i][1] * eq[j].c;
-	dy = origPts[i+1].y - origPts[i].y;
-	eq[i].b = (dy) / dx[i] - dx[i] * (eq[j].c + 2.0 * eq[i].c) / 3.0;
-	eq[i].d = (eq[j].c - eq[i].c) / (3.0 * dx[i]);
+  Point2d *ip, *iend;
+  double x, dy, alpha;
+  int isKnot;
+  int i, j, n;
+
+  double* dx = (double*)malloc(sizeof(double) * nOrigPts);
+  /* Calculate vector of differences */
+  for (i = 0, j = 1; j < nOrigPts; i++, j++) {
+    dx[i] = origPts[j].x - origPts[i].x;
+    if (dx[i] < 0.0) {
+      return 0;
     }
+  }
+  n = nOrigPts - 1;		/* Number of intervals. */
+  TriDiagonalMatrix* A = 
+    (TriDiagonalMatrix*)malloc(sizeof(TriDiagonalMatrix) * nOrigPts);
+  if (!A) {
+    free(dx);
+    return 0;
+  }
+  /* Vectors to solve the tridiagonal matrix */
+  A[0][0] = A[n][0] = 1.0;
+  A[0][1] = A[n][1] = 0.0;
+  A[0][2] = A[n][2] = 0.0;
+
+  /* Calculate the intermediate results */
+  for (i = 0, j = 1; j < n; j++, i++) {
+    alpha = 3.0 * ((origPts[j + 1].y / dx[j]) - (origPts[j].y / dx[i]) - 
+		   (origPts[j].y / dx[j]) + (origPts[i].y / dx[i]));
+    A[j][0] = 2 * (dx[j] + dx[i]) - dx[i] * A[i][1];
+    A[j][1] = dx[j] / A[j][0];
+    A[j][2] = (alpha - dx[i] * A[i][2]) / A[j][0];
+  }
+
+  Cubic2D* eq = (Cubic2D*)malloc(sizeof(Cubic2D) * nOrigPts);
+  if (!eq) {
     free(A);
     free(dx);
-
-    /* Now calculate the new values */
-    for (ip = intpPts, iend = ip + nIntpPts; ip < iend; ip++) {
-	ip->y = 0.0;
-	x = ip->x;
-
-	/* Is it outside the interval? */
-	if ((x < origPts[0].x) || (x > origPts[n].x)) {
-	    continue;
-	}
-	/* Search for the interval containing x in the point array */
-	i = Search(origPts, nOrigPts, x, &isKnot);
-	if (isKnot) {
-	    ip->y = origPts[i].y;
-	} else {
-	    i--;
-	    x -= origPts[i].x;
-	    ip->y = origPts[i].y + x * (eq[i].b + x * (eq[i].c + x * eq[i].d));
-	}
+    return 0;
+  }
+  eq[0].c = eq[n].c = 0.0;
+  for (j = n, i = n - 1; i >= 0; i--, j--) {
+    eq[i].c = A[i][2] - A[i][1] * eq[j].c;
+    dy = origPts[i+1].y - origPts[i].y;
+    eq[i].b = (dy) / dx[i] - dx[i] * (eq[j].c + 2.0 * eq[i].c) / 3.0;
+    eq[i].d = (eq[j].c - eq[i].c) / (3.0 * dx[i]);
+  }
+  free(A);
+  free(dx);
+
+  /* Now calculate the new values */
+  for (ip = intpPts, iend = ip + nIntpPts; ip < iend; ip++) {
+    ip->y = 0.0;
+    x = ip->x;
+
+    /* Is it outside the interval? */
+    if ((x < origPts[0].x) || (x > origPts[n].x)) {
+      continue;
     }
-    free(eq);
-    return 1;
+    /* Search for the interval containing x in the point array */
+    i = Search(origPts, nOrigPts, x, &isKnot);
+    if (isKnot) {
+      ip->y = origPts[i].y;
+    } else {
+      i--;
+      x -= origPts[i].x;
+      ip->y = origPts[i].y + x * (eq[i].b + x * (eq[i].c + x * eq[i].d));
+    }
+  }
+  free(eq);
+  return 1;
 }
 
 typedef struct {
-    double t;			/* Arc length of interval. */
-    double x;			/* 2nd derivative of X with respect to T */
-    double y;			/* 2nd derivative of Y with respect to T */
+  double t;			/* Arc length of interval. */
+  double x;			/* 2nd derivative of X with respect to T */
+  double y;			/* 2nd derivative of Y with respect to T */
 } CubicSpline;
 
 /*
@@ -900,42 +745,41 @@ typedef struct {
  * (C is upper triangle with unit diagonal, D is diagonal) is calculated.
  * Return TRUE if decomposition exist.
  */
-static int 
-SolveCubic1(TriDiagonalMatrix A[], int n)
+static int SolveCubic1(TriDiagonalMatrix A[], int n)
 {
-    int i;
-    double m_ij, m_n, m_nn, d;
-
-    if (n < 1) {
-	return 0;		/* Dimension should be at least 1 */
-    }
-    d = A[0][1];		/* D_{0,0} = A_{0,0} */
+  int i;
+  double m_ij, m_n, m_nn, d;
+
+  if (n < 1) {
+    return 0;		/* Dimension should be at least 1 */
+  }
+  d = A[0][1];		/* D_{0,0} = A_{0,0} */
+  if (d <= 0.0) {
+    return 0;		/* A (or D) should be positive definite */
+  }
+  m_n = A[0][0];		/*  A_{0,n-1}  */
+  m_nn = A[n - 1][1];		/* A_{n-1,n-1} */
+  for (i = 0; i < n - 2; i++) {
+    m_ij = A[i][2];		/*  A_{i,1}  */
+    A[i][2] = m_ij / d;	/* C_{i,i+1} */
+    A[i][0] = m_n / d;	/* C_{i,n-1} */
+    m_nn -= A[i][0] * m_n;	/* to get C_{n-1,n-1} */
+    m_n = -A[i][2] * m_n;	/* to get C_{i+1,n-1} */
+    d = A[i + 1][1] - A[i][2] * m_ij;	/* D_{i+1,i+1} */
     if (d <= 0.0) {
-	return 0;		/* A (or D) should be positive definite */
+      return 0;	/* Elements of D should be positive */
     }
-    m_n = A[0][0];		/*  A_{0,n-1}  */
-    m_nn = A[n - 1][1];		/* A_{n-1,n-1} */
-    for (i = 0; i < n - 2; i++) {
-	m_ij = A[i][2];		/*  A_{i,1}  */
-	A[i][2] = m_ij / d;	/* C_{i,i+1} */
-	A[i][0] = m_n / d;	/* C_{i,n-1} */
-	m_nn -= A[i][0] * m_n;	/* to get C_{n-1,n-1} */
-	m_n = -A[i][2] * m_n;	/* to get C_{i+1,n-1} */
-	d = A[i + 1][1] - A[i][2] * m_ij;	/* D_{i+1,i+1} */
-	if (d <= 0.0) {
-	    return 0;	/* Elements of D should be positive */
-	}
-	A[i + 1][1] = d;
-    }
-    if (n >= 2) {		/* Complete last column */
-	m_n += A[n - 2][2];	/* add A_{n-2,n-1} */
-	A[n - 2][0] = m_n / d;	/* C_{n-2,n-1} */
-	A[n - 1][1] = d = m_nn - A[n - 2][0] * m_n;	/* D_{n-1,n-1} */
-	if (d <= 0.0) {
-	    return 0;
-	}
+    A[i + 1][1] = d;
+  }
+  if (n >= 2) {		/* Complete last column */
+    m_n += A[n - 2][2];	/* add A_{n-2,n-1} */
+    A[n - 2][0] = m_n / d;	/* C_{n-2,n-1} */
+    A[n - 1][1] = d = m_nn - A[n - 2][0] * m_n;	/* D_{n-1,n-1} */
+    if (d <= 0.0) {
+      return 0;
     }
-    return 1;
+  }
+  return 1;
 }
 
 /*
@@ -943,49 +787,45 @@ SolveCubic1(TriDiagonalMatrix A[], int n)
  * decomposition calculated above (in m[][]) and the right side b given
  * in x[]. The solution x overwrites the right side in x[].
  */
-static void 
-SolveCubic2(TriDiagonalMatrix A[], CubicSpline spline[], int nIntervals)
+static void SolveCubic2(TriDiagonalMatrix A[], CubicSpline spline[], 
+			int nIntervals)
 {
-    int i;
-    double x, y;
-    int n, m;
-
-    n = nIntervals - 2;
-    m = nIntervals - 1;
-
-    /* Division by transpose of C : b = C^{-T} * b */
-    x = spline[m].x;
-    y = spline[m].y;
-    for (i = 0; i < n; i++) {
-	spline[i + 1].x -= A[i][2] * spline[i].x; /* C_{i,i+1} * x(i) */
-	spline[i + 1].y -= A[i][2] * spline[i].y; /* C_{i,i+1} * x(i) */
-	x -= A[i][0] * spline[i].x;	/* C_{i,n-1} * x(i) */
-	y -= A[i][0] * spline[i].y; /* C_{i,n-1} * x(i) */
-    }
-    if (n >= 0) {
-	/* C_{n-2,n-1} * x_{n-1} */
-	spline[m].x = x - A[n][0] * spline[n].x; 
-	spline[m].y = y - A[n][0] * spline[n].y; 
-    }
-    /* Division by D: b = D^{-1} * b */
-    for (i = 0; i < nIntervals; i++) {
-	spline[i].x /= A[i][1];
-	spline[i].y /= A[i][1];
-    }
-
-    /* Division by C: b = C^{-1} * b */
-    x = spline[m].x;
-    y = spline[m].y;
-    if (n >= 0) {
-	/* C_{n-2,n-1} * x_{n-1} */
-	spline[n].x -= A[n][0] * x;	
-	spline[n].y -= A[n][0] * y;	
-    }
-    for (i = (n - 1); i >= 0; i--) {
-	/* C_{i,i+1} * x_{i+1} + C_{i,n-1} * x_{n-1} */
-	spline[i].x -= A[i][2] * spline[i + 1].x + A[i][0] * x;
-	spline[i].y -= A[i][2] * spline[i + 1].y + A[i][0] * y;
-    }
+  int n = nIntervals - 2;
+  int m = nIntervals - 1;
+
+  // Division by transpose of C : b = C^{-T} * b 
+  double x = spline[m].x;
+  double y = spline[m].y;
+  for (int ii=0; ii<n; ii++) {
+    spline[ii + 1].x -= A[ii][2] * spline[ii].x; /* C_{i,i+1} * x(i) */
+    spline[ii + 1].y -= A[ii][2] * spline[ii].y; /* C_{i,i+1} * x(i) */
+    x -= A[ii][0] * spline[ii].x;	/* C_{i,n-1} * x(i) */
+    y -= A[ii][0] * spline[ii].y; /* C_{i,n-1} * x(i) */
+  }
+  if (n >= 0) {
+    // C_{n-2,n-1} * x_{n-1}
+    spline[m].x = x - A[n][0] * spline[n].x; 
+    spline[m].y = y - A[n][0] * spline[n].y; 
+  }
+  // Division by D: b = D^{-1} * b 
+  for (int ii=0; ii<nIntervals; ii++) {
+    spline[ii].x /= A[ii][1];
+    spline[ii].y /= A[ii][1];
+  }
+
+  // Division by C: b = C^{-1} * b
+  x = spline[m].x;
+  y = spline[m].y;
+  if (n >= 0) {
+    // C_{n-2,n-1} * x_{n-1}
+    spline[n].x -= A[n][0] * x;	
+    spline[n].y -= A[n][0] * y;	
+  }
+  for (int ii=(n - 1); ii>=0; ii--) {
+    // C_{i,i+1} * x_{i+1} + C_{i,n-1} * x_{n-1}
+    spline[ii].x -= A[ii][2] * spline[ii + 1].x + A[ii][0] * x;
+    spline[ii].y -= A[ii][2] * spline[ii + 1].y + A[ii][0] * y;
+  }
 }
 
 /*
@@ -994,286 +834,253 @@ SolveCubic2(TriDiagonalMatrix A[], CubicSpline spline[], int nIntervals)
  * length of the linear stroke. The number of points must be at least 3.
  * Note: For CLOSED_CONTOURs the first and last point must be equal.
  */
-static CubicSpline *
-CubicSlopes(
-    Point2d points[],
-    int nPoints,		/* Number of points (nPoints>=3) */
-    int isClosed,		/* CLOSED_CONTOUR or OPEN_CONTOUR  */
-    double unitX, 
-    double unitY)		/* Unit length in x and y (norm=1) */
+static CubicSpline* CubicSlopes(Point2d points[], int nPoints,
+				int isClosed, double unitX, double unitY)
 {
-    CubicSpline *s1, *s2;
-    int n, i;
-    double norm, dx, dy;
+  CubicSpline *s1, *s2;
+  int n, i;
+  double norm, dx, dy;
     
-    CubicSpline* spline = (CubicSpline*)malloc(sizeof(CubicSpline) * nPoints);
-    if (!spline)
-	return NULL;
-
-    TriDiagonalMatrix *A 
-      = (TriDiagonalMatrix*)malloc(sizeof(TriDiagonalMatrix) * nPoints);
-    if (!A) {
-	free(spline);
-	return NULL;
-    }
+  CubicSpline* spline = (CubicSpline*)malloc(sizeof(CubicSpline) * nPoints);
+  if (!spline)
+    return NULL;
+
+  TriDiagonalMatrix *A 
+    = (TriDiagonalMatrix*)malloc(sizeof(TriDiagonalMatrix) * nPoints);
+  if (!A) {
+    free(spline);
+    return NULL;
+  }
+  /*
+   * Calculate first differences in (dxdt2[i], y[i]) and interval lengths
+   * in dist[i]:
+   */
+  s1 = spline;
+  for (i = 0; i < nPoints - 1; i++) {
+    s1->x = points[i+1].x - points[i].x;
+    s1->y = points[i+1].y - points[i].y;
+
     /*
-     * Calculate first differences in (dxdt2[i], y[i]) and interval lengths
-     * in dist[i]:
+     * The Norm of a linear stroke is calculated in "normal coordinates"
+     * and used as interval length:
+     */
+    dx = s1->x / unitX;
+    dy = s1->y / unitY;
+    s1->t = sqrt(dx * dx + dy * dy);
+
+    s1->x /= s1->t;	/* first difference, with unit norm: */
+    s1->y /= s1->t;	/*   || (dxdt2[i], y[i]) || = 1      */
+    s1++;
+  }
+
+  /*
+   * Setup linear System:  Ax = b
+   */
+  n = nPoints - 2;		/* Without first and last point */
+  if (isClosed) {
+    /* First and last points must be equal for CLOSED_CONTOURs */
+    spline[nPoints - 1].t = spline[0].t;
+    spline[nPoints - 1].x = spline[0].x;
+    spline[nPoints - 1].y = spline[0].y;
+    n++;			/* Add last point (= first point) */
+  }
+  s1 = spline, s2 = s1 + 1;
+  for (i = 0; i < n; i++) {
+    /* Matrix A, mainly tridiagonal with cyclic second index 
+       ("j = j+n mod n") 
+    */
+    A[i][0] = s1->t;	/* Off-diagonal element A_{i,i-1} */
+    A[i][1] = 2.0 * (s1->t + s2->t);	/* A_{i,i} */
+    A[i][2] = s2->t;	/* Off-diagonal element A_{i,i+1} */
+
+    /* Right side b_x and b_y */
+    s1->x = (s2->x - s1->x) * 6.0;
+    s1->y = (s2->y - s1->y) * 6.0;
+
+    /* 
+     * If the linear stroke shows a cusp of more than 90 degree,
+     * the right side is reduced to avoid oscillations in the
+     * spline: 
      */
-    s1 = spline;
-    for (i = 0; i < nPoints - 1; i++) {
-	s1->x = points[i+1].x - points[i].x;
-	s1->y = points[i+1].y - points[i].y;
-
-	/*
-	 * The Norm of a linear stroke is calculated in "normal coordinates"
-	 * and used as interval length:
-	 */
-	dx = s1->x / unitX;
-	dy = s1->y / unitY;
-	s1->t = sqrt(dx * dx + dy * dy);
-
-	s1->x /= s1->t;	/* first difference, with unit norm: */
-	s1->y /= s1->t;	/*   || (dxdt2[i], y[i]) || = 1      */
-	s1++;
-    }
-
     /*
-     * Setup linear System:  Ax = b
+     * The Norm of a linear stroke is calculated in "normal coordinates"
+     * and used as interval length:
      */
-    n = nPoints - 2;		/* Without first and last point */
-    if (isClosed) {
-	/* First and last points must be equal for CLOSED_CONTOURs */
-	spline[nPoints - 1].t = spline[0].t;
-	spline[nPoints - 1].x = spline[0].x;
-	spline[nPoints - 1].y = spline[0].y;
-	n++;			/* Add last point (= first point) */
-    }
-    s1 = spline, s2 = s1 + 1;
-    for (i = 0; i < n; i++) {
-	/* Matrix A, mainly tridiagonal with cyclic second index 
-	   ("j = j+n mod n") 
-	*/
-	A[i][0] = s1->t;	/* Off-diagonal element A_{i,i-1} */
-	A[i][1] = 2.0 * (s1->t + s2->t);	/* A_{i,i} */
-	A[i][2] = s2->t;	/* Off-diagonal element A_{i,i+1} */
-
-	/* Right side b_x and b_y */
-	s1->x = (s2->x - s1->x) * 6.0;
-	s1->y = (s2->y - s1->y) * 6.0;
-
-	/* 
-	 * If the linear stroke shows a cusp of more than 90 degree,
-	 * the right side is reduced to avoid oscillations in the
-	 * spline: 
-	 */
-	/*
-	 * The Norm of a linear stroke is calculated in "normal coordinates"
-	 * and used as interval length:
-	 */
-	dx = s1->x / unitX;
-	dy = s1->y / unitY;
-	norm = sqrt(dx * dx + dy * dy) / 8.5;
-	if (norm > 1.0) {
-	    /* The first derivative will not be continuous */
-	    s1->x /= norm;
-	    s1->y /= norm;
-	}
-	s1++, s2++;
-    }
-
-    if (!isClosed) {
-	/* Third derivative is set to zero at both ends */
-	A[0][1] += A[0][0];	/* A_{0,0}     */
-	A[0][0] = 0.0;		/* A_{0,n-1}   */
-	A[n-1][1] += A[n-1][2]; /* A_{n-1,n-1} */
-	A[n-1][2] = 0.0;	/* A_{n-1,0}   */
-    }
-    /* Solve linear systems for dxdt2[] and y[] */
-
-    if (SolveCubic1(A, n)) {	/* Cholesky decomposition */
-	SolveCubic2(A, spline, n); /* A * dxdt2 = b_x */
-    } else {			/* Should not happen, but who knows ... */
-	free(A);
-	free(spline);
-	return NULL;
+    dx = s1->x / unitX;
+    dy = s1->y / unitY;
+    norm = sqrt(dx * dx + dy * dy) / 8.5;
+    if (norm > 1.0) {
+      /* The first derivative will not be continuous */
+      s1->x /= norm;
+      s1->y /= norm;
     }
-    /* Shift all second derivatives one place right and update the ends. */
-    s2 = spline + n, s1 = s2 - 1;
-    for (/* empty */; s2 > spline; s2--, s1--) {
-	s2->x = s1->x;
-	s2->y = s1->y;
-    }
-    if (isClosed) {
-	spline[0].x = spline[n].x;
-	spline[0].y = spline[n].y;
-    } else {
-	/* Third derivative is 0.0 for the first and last interval. */
-	spline[0].x = spline[1].x; 
-	spline[0].y = spline[1].y; 
-	spline[n + 1].x = spline[n].x;
-	spline[n + 1].y = spline[n].y;
-    }
-    free( A);
-    return spline;
+    s1++, s2++;
+  }
+
+  if (!isClosed) {
+    /* Third derivative is set to zero at both ends */
+    A[0][1] += A[0][0];	/* A_{0,0}     */
+    A[0][0] = 0.0;		/* A_{0,n-1}   */
+    A[n-1][1] += A[n-1][2]; /* A_{n-1,n-1} */
+    A[n-1][2] = 0.0;	/* A_{n-1,0}   */
+  }
+  /* Solve linear systems for dxdt2[] and y[] */
+
+  if (SolveCubic1(A, n)) {	/* Cholesky decomposition */
+    SolveCubic2(A, spline, n); /* A * dxdt2 = b_x */
+  } else {			/* Should not happen, but who knows ... */
+    free(A);
+    free(spline);
+    return NULL;
+  }
+  /* Shift all second derivatives one place right and update the ends. */
+  s2 = spline + n, s1 = s2 - 1;
+  for (/* empty */; s2 > spline; s2--, s1--) {
+    s2->x = s1->x;
+    s2->y = s1->y;
+  }
+  if (isClosed) {
+    spline[0].x = spline[n].x;
+    spline[0].y = spline[n].y;
+  } else {
+    /* Third derivative is 0.0 for the first and last interval. */
+    spline[0].x = spline[1].x; 
+    spline[0].y = spline[1].y; 
+    spline[n + 1].x = spline[n].x;
+    spline[n + 1].y = spline[n].y;
+  }
+  free( A);
+  return spline;
 }
 
-
-/*
- * Calculate interpolated values of the spline function (defined via p_cntr
- * and the second derivatives dxdt2[] and dydt2[]). The number of tabulated
- * values is n. On an equidistant grid n_intpol values are calculated.
- */
-static int
-CubicEval(Point2d *origPts, int nOrigPts, Point2d *intpPts, int nIntpPts,
-	  CubicSpline *spline)
+// Calculate interpolated values of the spline function (defined via p_cntr
+// and the second derivatives dxdt2[] and dydt2[]). The number of tabulated
+// values is n. On an equidistant grid n_intpol values are calculated.
+static int CubicEval(Point2d *origPts, int nOrigPts, Point2d *intpPts, 
+		     int nIntpPts, CubicSpline *spline)
 {
-    double t, tSkip, tMax;
-    Point2d q;
-    int i, j, count;
-
-    /* Sum the lengths of all the segments (intervals). */
-    tMax = 0.0;
-    for (i = 0; i < nOrigPts - 1; i++) {
-	tMax += spline[i].t;
-    }
+  double t, tSkip;
+  Point2d q;
+  int count;
+
+  /* Sum the lengths of all the segments (intervals). */
+  double tMax = 0.0;
+  for (int ii=0; ii<nOrigPts - 1; ii++)
+    tMax += spline[ii].t;
 
-    /* Need a better way of doing this... */
+  /* Need a better way of doing this... */
 
-    /* The distance between interpolated points */
-    tSkip = (1. - 1e-7) * tMax / (nIntpPts - 1);
+  /* The distance between interpolated points */
+  tSkip = (1. - 1e-7) * tMax / (nIntpPts - 1);
     
-    t = 0.0;			/* Spline parameter value. */
-    q = origPts[0];
-    count = 0;
+  t = 0.0;			/* Spline parameter value. */
+  q = origPts[0];
+  count = 0;
 
-    intpPts[count++] = q;	/* First point. */
-    t += tSkip;
+  intpPts[count++] = q;	/* First point. */
+  t += tSkip;
     
-    for (i = 0, j = 1; j < nOrigPts; i++, j++) {
-	Point2d p;
-	double d, hx, dx0, dx01, hy, dy0, dy01;
-	
-	d = spline[i].t;	/* Interval length */
-	p = q;
-	q = origPts[i+1];
-	hx = (q.x - p.x) / d;
-	hy = (q.y - p.y) / d;
-	dx0 = (spline[j].x + 2 * spline[i].x) / 6.0;
-	dy0 = (spline[j].y + 2 * spline[i].y) / 6.0;
-	dx01 = (spline[j].x - spline[i].x) / (6.0 * d);
-	dy01 = (spline[j].y - spline[i].y) / (6.0 * d);
-	while (t <= spline[i].t) { /* t in current interval ? */
-	    p.x += t * (hx + (t - d) * (dx0 + t * dx01));
-	    p.y += t * (hy + (t - d) * (dy0 + t * dy01));
-	    intpPts[count++] = p;
-	    t += tSkip;
-	}
-	/* Parameter t relative to start of next interval */
-	t -= spline[i].t;
+  for (int ii=0, jj=1; jj<nOrigPts; ii++, jj++) {
+    // Interval length
+    double d = spline[ii].t;
+    Point2d p = q;
+    q = origPts[ii+1];
+    double hx = (q.x - p.x) / d;
+    double hy = (q.y - p.y) / d;
+    double dx0 = (spline[jj].x + 2 * spline[ii].x) / 6.0;
+    double dy0 = (spline[jj].y + 2 * spline[ii].y) / 6.0;
+    double dx01 = (spline[jj].x - spline[ii].x) / (6.0 * d);
+    double dy01 = (spline[jj].y - spline[ii].y) / (6.0 * d);
+    while (t <= spline[ii].t) { /* t in current interval ? */
+      p.x += t * (hx + (t - d) * (dx0 + t * dx01));
+      p.y += t * (hy + (t - d) * (dy0 + t * dy01));
+      intpPts[count++] = p;
+      t += tSkip;
     }
-    return count;
-}
 
-/*
- * Generate a cubic spline curve through the points (x_i,y_i) which are
- * stored in the linked list p_cntr.
- * The spline is defined as a 2d-function s(t) = (x(t),y(t)), where the
- * parameter t is the length of the linear stroke.
- */
-int
-Blt_NaturalParametricSpline(Point2d *origPts, int nOrigPts, Region2d *extsPtr,
-			    int isClosed, Point2d *intpPts, int nIntpPts)
-{
-    double unitX, unitY;	/* To define norm (x,y)-plane */
-    CubicSpline *spline;
-    int result;
+    // Parameter t relative to start of next interval 
+    t -= spline[ii].t;
+  }
 
-    if (nOrigPts < 3) {
-	return 0;
-    }
-    if (isClosed) {
-	origPts[nOrigPts].x = origPts[0].x;
-	origPts[nOrigPts].y = origPts[0].y;
-	nOrigPts++;
-    }
-    /* Width and height of the grid is used at unit length (2d-norm) */
-    unitX = extsPtr->right - extsPtr->left;
-    unitY = extsPtr->bottom - extsPtr->top;
+  return count;
+}
 
-    if (unitX < FLT_EPSILON) {
-	unitX = FLT_EPSILON;
-    }
-    if (unitY < FLT_EPSILON) {
-	unitY = FLT_EPSILON;
-    }
-    /* Calculate parameters for cubic spline: 
-     *		t     = arc length of interval.
-     *		dxdt2 = second derivatives of x with respect to t, 
-     *		dydt2 = second derivatives of y with respect to t, 
-     */
-    spline = CubicSlopes(origPts, nOrigPts, isClosed, unitX, unitY);
-    if (spline == NULL) {
-	return 0;
-    }
-    result= CubicEval(origPts, nOrigPts, intpPts, nIntpPts, spline);
-    free(spline);
-    return result;
+int LineElement::naturalParametricSpline(Point2d *origPts, int nOrigPts, 
+					 Region2d *extsPtr, int isClosed,
+					 Point2d *intpPts, int nIntpPts)
+{
+  // Generate a cubic spline curve through the points (x_i,y_i) which are
+  // stored in the linked list p_cntr.
+  // The spline is defined as a 2d-function s(t) = (x(t),y(t)), where the
+  // parameter t is the length of the linear stroke.
+
+  if (nOrigPts < 3)
+    return 0;
+
+  if (isClosed) {
+    origPts[nOrigPts].x = origPts[0].x;
+    origPts[nOrigPts].y = origPts[0].y;
+    nOrigPts++;
+  }
+
+  // Width and height of the grid is used at unit length (2d-norm)
+  double unitX = extsPtr->right - extsPtr->left;
+  double unitY = extsPtr->bottom - extsPtr->top;
+  if (unitX < FLT_EPSILON)
+    unitX = FLT_EPSILON;
+  if (unitY < FLT_EPSILON)
+    unitY = FLT_EPSILON;
+
+  /* Calculate parameters for cubic spline: 
+   *		t     = arc length of interval.
+   *		dxdt2 = second derivatives of x with respect to t, 
+   *		dydt2 = second derivatives of y with respect to t, 
+   */
+  CubicSpline* spline = CubicSlopes(origPts, nOrigPts, isClosed, unitX, unitY);
+  if (spline == NULL)
+    return 0;
+
+  int result= CubicEval(origPts, nOrigPts, intpPts, nIntpPts, spline);
+
+  free(spline);
+  return result;
 }
 
-static INLINE void
-CatromCoeffs(Point2d *p, Point2d *a, Point2d *b, Point2d *c, Point2d *d)
+static void CatromCoeffs(Point2d* p, Point2d* a, Point2d* b, 
+			 Point2d* c, Point2d* d)
 {
-    a->x = -p[0].x + 3.0 * p[1].x - 3.0 * p[2].x + p[3].x;
-    b->x = 2.0 * p[0].x - 5.0 * p[1].x + 4.0 * p[2].x - p[3].x;
-    c->x = -p[0].x + p[2].x;
-    d->x = 2.0 * p[1].x;
-    a->y = -p[0].y + 3.0 * p[1].y - 3.0 * p[2].y + p[3].y;
-    b->y = 2.0 * p[0].y - 5.0 * p[1].y + 4.0 * p[2].y - p[3].y;
-    c->y = -p[0].y + p[2].y;
-    d->y = 2.0 * p[1].y;
+  a->x = -p[0].x + 3.0 * p[1].x - 3.0 * p[2].x + p[3].x;
+  b->x = 2.0 * p[0].x - 5.0 * p[1].x + 4.0 * p[2].x - p[3].x;
+  c->x = -p[0].x + p[2].x;
+  d->x = 2.0 * p[1].x;
+  a->y = -p[0].y + 3.0 * p[1].y - 3.0 * p[2].y + p[3].y;
+  b->y = 2.0 * p[0].y - 5.0 * p[1].y + 4.0 * p[2].y - p[3].y;
+  c->y = -p[0].y + p[2].y;
+  d->y = 2.0 * p[1].y;
 }
 
-/*
- *---------------------------------------------------------------------------
- *
- * Blt_ParametricCatromSpline --
- *
- *	Computes a spline based upon the data points, returning a new (larger)
- *	coordinate array of points.
- *
- * Results:
- *	None.
- *
- *---------------------------------------------------------------------------
- */
-int
-Blt_CatromParametricSpline(Point2d *points, int nPoints, Point2d *intpPts,
-			   int nIntpPts)
+int LineElement::catromParametricSpline(Point2d* points, int nPoints, 
+					Point2d* intpPts, int nIntpPts)
 {
-    int i;
-    double t;
-    int interval;
+  // The spline is computed in screen coordinates instead of data points so
+  // that we can select the abscissas of the interpolated points from each
+  // pixel horizontally across the plotting area.
+
+  Point2d* origPts = new Point2d[nPoints + 4];
+  memcpy(origPts + 1, points, sizeof(Point2d) * nPoints);
+
+  origPts[0] = origPts[1];
+  origPts[nPoints + 2] = origPts[nPoints + 1] = origPts[nPoints];
+
+  for (int ii=0; ii<nIntpPts; ii++) {
+    int interval = (int)intpPts[ii].x;
+    double t = intpPts[ii].y;
     Point2d a, b, c, d;
+    CatromCoeffs(origPts + interval, &a, &b, &c, &d);
+    intpPts[ii].x = (d.x + t * (c.x + t * (b.x + t * a.x))) / 2.0;
+    intpPts[ii].y = (d.y + t * (c.y + t * (b.y + t * a.y))) / 2.0;
+  }
 
-    /*
-     * The spline is computed in screen coordinates instead of data points so
-     * that we can select the abscissas of the interpolated points from each
-     * pixel horizontally across the plotting area.
-     */
-    Point2d* origPts = (Point2d*)malloc((nPoints + 4) * sizeof(Point2d));
-    memcpy(origPts + 1, points, sizeof(Point2d) * nPoints);
-
-    origPts[0] = origPts[1];
-    origPts[nPoints + 2] = origPts[nPoints + 1] = origPts[nPoints];
-
-    for (i = 0; i < nIntpPts; i++) {
-	interval = (int)intpPts[i].x;
-	t = intpPts[i].y;
-	CatromCoeffs(origPts + interval, &a, &b, &c, &d);
-	intpPts[i].x = (d.x + t * (c.x + t * (b.x + t * a.x))) / 2.0;
-	intpPts[i].y = (d.y + t * (c.y + t * (b.y + t * a.y))) / 2.0;
-    }
-    free(origPts);
-    return 1;
+  delete [] origPts;
+  return 1;
 }