_IMPLEMENTING CURVES IN C++_
by Stephen P. Johnson and Tom McReynolds


[LISTING ONE]
// curve.h - Base class for curves
#include <stdio.h>
class Point2d {
protected:
    float   x, y, w;
public:
    // set x, y
    void set_xy(float new_x, float new_y) {
         x = new_x; y = new_y; w = 1.;
    }
    // set x, y and w
    void set_xyw(float new_x, float new_y, float new_w) {
         x = new_x; y = new_y; w = new_w;
    }
    void set_x(float new_x) { x = new_x; }
    void set_y(float new_y) { y = new_y; }
    void set_w(float new_w) { w = new_w; }
    // get x, y
    void get_xy(float *ret_x, float *ret_y) {
         *ret_x = x; *ret_y = y;
    }
    // get x, y and w
    void get_xy(float *ret_x, float *ret_y, float *ret_w) {
         *ret_x = x; *ret_y = y; *ret_w = w;
    }
    float get_x(void) { return x; }
    float get_y(void) { return y; }
    float get_w(void) { return w; }
    // print out the current x, y
    void print(void) { printf("%f %f\n", x, y); }
};
class Curve {
protected:
    int     num_geom;
    Point2d *geom;
public:
    // constructor
    Curve(void);
    // destructor
    ~Curve(void);
    // set the geometry vector which is also called the control points
    void set_geom_vector(int count, Point2d *g);
    // method to display a third degree curve
    void display_curve(int n, int color, int show_geom_pts,
    int show_convex_hull);
};
class Basis_matrix_curve : public Curve {
protected:
    Point2d t_coeff[4];
    float   basis_matrix[4][4];
public:
    // constructor
    Basis_matrix_curve(void);
    // get the coefficient matrix for the t's
    void get_t_coeff(Point2d tc[4]);
    // overwrite set geometry vector to add in matrix multiply
    void set_geom_vector(int count, Point2d *g);
    // set a user defined matrix into the basis matrix
    void set_basis_matrix(float m[4][4]);
    // get the current basis matrix
    void get_basis_matrix(float m[4][4]);
    // multiply a 4x4 by 4 points that are a 2x4 matrix
    void multiply_basis_by_geometry(void);
    // method to display a third degree curve
    void display_curve(int n, int color, int show_geom_pts,
    int show_convex_hull);
};
class Hermite_curve : public Basis_matrix_curve {
public:
    // constructor for Hermite type curve
    Hermite_curve(void);
};
class Bezier_curve : public Basis_matrix_curve {
public:
    // constructor for Bezier type curve
    Bezier_curve(void);
};
class Bspline_curve : public Basis_matrix_curve {
public:
    // constructor for Bspline type curve
    Bspline_curve(void);
};
class Catmull_Rom_curve : public Basis_matrix_curve {
public:
    // constructor for Catmull-Rom type curve
    Catmull_Rom_curve(void);
};
// uniformly shaped beta-spline
class Beta_spline_curve : public Basis_matrix_curve {
protected:
    float   bias;
    float   tension;
public:
    // constructor for uniformly shaped beta-spline
    Beta_spline_curve(void);
    // methods for setting and getting the bias and tension
    void set_bias(float new_bias) {
         bias = new_bias;
         update_basis_matrix();
         multiply_basis_by_geometry();
    }
    void set_tension(float new_tension) {
         tension = new_tension;
         update_basis_matrix();
         multiply_basis_by_geometry();
    }
    float get_bias(void) { return bias; }
    float get_tension(void) { return tension; }
    // method for updating the basis matrix
    void update_basis_matrix(void);
};
class Nub_curve : public Curve {
protected:
    int     num_knots;
    float   *knots;
public:
    // constructor for NUB curves
    Nub_curve(void);
    // destructor for NUB curves
    ~Nub_curve(void);
    // setup the knot vector with user information
    void set_knot_vector(int count, float *v);
    // method to display a third degree NUB curve
    void display_curve(int n, int color, int show_geom_pts,
    int show_convex_hull);
};
class Nurb_curve : public Nub_curve {
public:
    // method to display a third degree NURB curve
    void display_curve(int n, int color, int show_geom_pts,
    int show_convex_hull);
};



[LISTING TWO]

// curve.cc - Implementation of methods for base class for curves //
#include <iostream.h>
#include <stdio.h>
#include <stdlib.h>
#ifdef __TURBOC__
#include <alloc.h>
#else
#include <malloc.h>
#endif
#include <math.h>
#include "curve.h"

#define CROSS_SIZE  5
extern "C" void line(short, short, short, short, short);

// constructor for curve
Curve::Curve(void)
{
    num_geom = 0;
    geom = (Point2d *)NULL;
}
// destructor for curve
Curve::~Curve(void)
{
    if (num_geom != 0 || geom) {
         free(geom);
    }
}
// set the geometry vector which is also called the control points
void Curve::set_geom_vector(int count, Point2d *g)
{
    if (count != num_geom && geom) {
    free(geom);
    }
    geom = (Point2d *)malloc(count * sizeof(Point2d));
    if (geom) {
         int i;
         num_geom = count;
         for (i = 0; i < num_geom; i++) {
         geom[i] = g[i];
         }
    }
}
// constructor for basis matrix class
// get the coefficients
void Basis_matrix_curve::get_t_coeff(Point2d tc[4])
{
    tc[0] = t_coeff[0];
    tc[1] = t_coeff[1];
    tc[2] = t_coeff[2];
    tc[3] = t_coeff[3];
}
// overwrite set geometry vector to do basis matrix multiply
void Basis_matrix_curve::set_geom_vector(int count, Point2d *g)
{
    Curve::set_geom_vector(count, g);
    multiply_basis_by_geometry();
}
// constructor for basis matrix class
Basis_matrix_curve::Basis_matrix_curve(void)
{
    int     i, j;
    Point2d g[4];

    g[0].set_xy(0., 0.); g[1].set_xy(0., 0.);
    g[2].set_xy(0., 0.); g[3].set_xy(0., 0.);
    Curve::set_geom_vector(4, g);

    t_coeff[0].set_xy(0., 0.);
    t_coeff[1].set_xy(0., 0.);
    t_coeff[2].set_xy(0., 0.);
    t_coeff[3].set_xy(0., 0.);
         for (i = 0; i < 4; i++)
         for (j = 0; j < 4; j++)
             if (i == j)
             basis_matrix[i][j] = 1.;
             else
             basis_matrix[i][j] = 0.;
}
// set a user defined matrix into the basis matrix
void Basis_matrix_curve::set_basis_matrix(float m[4][4])
{
    int i, j;
    for (i = 0; i < 4; i++)
         for (j = 0; j < 4; j++)
             basis_matrix[i][j] = m[i][j];
}
// get the current basis matrix
void Basis_matrix_curve::get_basis_matrix(float m[4][4])
{
    int i, j;
    for (i = 0; i < 4; i++)
         for (j = 0; j < 4; j++)
             m[i][j] = basis_matrix[i][j];
}

// multiply a 4x4 by a 2x4 (2D geometric matrix)
void Basis_matrix_curve::multiply_basis_by_geometry()
{
    t_coeff[0].set_x(
    basis_matrix[0][0] * geom[0].get_x() +
    basis_matrix[0][1] * geom[1].get_x() +
    basis_matrix[0][2] * geom[2].get_x() +
    basis_matrix[0][3] * geom[3].get_x());

    t_coeff[1].set_x(
    basis_matrix[1][0] * geom[0].get_x() +
    basis_matrix[1][1] * geom[1].get_x() +
    basis_matrix[1][2] * geom[2].get_x() +
    basis_matrix[1][3] * geom[3].get_x());
    t_coeff[2].set_x(
    basis_matrix[2][0] * geom[0].get_x() +
    basis_matrix[2][1] * geom[1].get_x() +
    basis_matrix[2][2] * geom[2].get_x() +
    basis_matrix[2][3] * geom[3].get_x());
    t_coeff[3].set_x(
    basis_matrix[3][0] * geom[0].get_x() +
    basis_matrix[3][1] * geom[1].get_x() +
    basis_matrix[3][2] * geom[2].get_x() +
    basis_matrix[3][3] * geom[3].get_x());
    t_coeff[0].set_y(
    basis_matrix[0][0] * geom[0].get_y() +
    basis_matrix[0][1] * geom[1].get_y() +
    basis_matrix[0][2] * geom[2].get_y() +
    basis_matrix[0][3] * geom[3].get_y());
    t_coeff[1].set_y(
    basis_matrix[1][0] * geom[0].get_y() +
    basis_matrix[1][1] * geom[1].get_y() +
    basis_matrix[1][2] * geom[2].get_y() +
    basis_matrix[1][3] * geom[3].get_y());

    t_coeff[2].set_y(
    basis_matrix[2][0] * geom[0].get_y() +
    basis_matrix[2][1] * geom[1].get_y() +
    basis_matrix[2][2] * geom[2].get_y() +
    basis_matrix[2][3] * geom[3].get_y());

    t_coeff[3].set_y(
    basis_matrix[3][0] * geom[0].get_y() +
    basis_matrix[3][1] * geom[1].get_y() +
    basis_matrix[3][2] * geom[2].get_y() +
    basis_matrix[3][3] * geom[3].get_y());
}
// constructor for Hermit curve
Hermite_curve::Hermite_curve(void)
{
    float   m[4][4];
    m[0][0] = 2.; m[0][1] = -2.; m[0][2] = 1.; m[0][3] = 1.;
    m[1][0] = -3.; m[1][1] = 3.; m[1][2] = -2.; m[1][3] = -1.;
    m[2][0] = 0.; m[2][1] = 0.; m[2][2] = 1.; m[2][3] = 0.;
    m[3][0] = 1.; m[3][1] = 0.; m[3][2] = 0.; m[3][3] = 0.;
    set_basis_matrix(m);
}
// constructor for Bezier curve
Bezier_curve::Bezier_curve(void)
{
    float   m[4][4];
    m[0][0] = -1.; m[0][1] = 3.; m[0][2] = -3.; m[0][3] = 1.;
    m[1][0] = 3.; m[1][1] = -6; m[1][2] = 3.; m[1][3] = 0.;
    m[2][0] = -3.; m[2][1] = 3.; m[2][2] = 0.; m[2][3] = 0.;
    m[3][0] = 1.; m[3][1] = 0.; m[3][2] = 0.; m[3][3] = 0.;
    set_basis_matrix(m);
}
// constructor for Uniform Nonrational Bspline curve
Bspline_curve::Bspline_curve(void)
{
    float   m[4][4];
    m[0][0] = -1./6.; m[0][1] = 3./6.; m[0][2] = -3./6.; m[0][3] = 1./6.;
    m[1][0] = 3./6.; m[1][1] = -6./6.; m[1][2] = 3./6.; m[1][3] = 0.;
    m[2][0] = -3./6.; m[2][1] = 0.; m[2][2] = 3./6.; m[2][3] = 0.;
    m[3][0] = 1./6.; m[3][1] = 4./6.; m[3][2] = 1./6.; m[3][3] = 0.;
    set_basis_matrix(m);
}
Catmull_Rom_curve::Catmull_Rom_curve(void)
{
    float   m[4][4];
    m[0][0] = -1./2.; m[0][1] = 3./2.; m[0][2] = -3./2.; m[0][3] = 1./2.;
    m[1][0] = 2./2.; m[1][1] = -5./2.; m[1][2] = 4./2.; m[1][3] = -1./2.;
    m[2][0] = -1./2.; m[2][1] = 0.; m[2][2] = 1./2.; m[2][3] = 0.;
    m[3][0] = 0.; m[3][1] = 2./2.; m[3][2] = 0.; m[3][3] = 0.;
    set_basis_matrix(m);
}
Beta_spline_curve::Beta_spline_curve(void)
{
    bias = 1.;
    tension = 0.;
    update_basis_matrix();
}
void Beta_spline_curve::update_basis_matrix(void)
{
    int     i, j;
    float   m[4][4];
    float   bias2 = bias*bias;
    float   bias3 = bias2*bias;
    float   delta;

    m[0][0] = -2.*bias3; m[0][1] = 2.*(tension+bias3+bias2+bias);
    m[0][2] = -2.*(tension+bias2+bias+1.); m[0][3] = 2.;

    m[1][0] = 6.*bias3; m[1][1] = -3.*(tension+2.*bias3+2.*bias2);
    m[1][2] = 3.*(tension+2.*bias2); m[1][3] = 0.;

    m[2][0] = -6.*bias3; m[2][1] = 6.*(bias3-bias);
    m[2][2] = 6.*bias; m[2][3] = 0.;

    m[3][0] = 2.*bias3; m[3][1] = tension+4.*(bias2+bias);
    m[3][2] = 2.; m[3][3] = 0.;
    delta = tension + 2.*bias3 + 4*bias2 + 4*bias + 2.;
    for (i = 0; i < 4; i++)
         for (j = 0; j < 4; j++)
              m[i][j] /= delta;
    set_basis_matrix(m);
}
float b1(int i, float t, float *knots)
{
    if (knots[i] <= 1. && t < knots[i+1])
         return 1.;
    else
         return 0.;
}
float b2(int i, float t, float *knots)
{
    float   n, d;
    float   sum = 0.;

    n = t - knots[i];
    d = knots[i+1] - knots[i];
    if (d != 0.) {
         sum += (n / d) * b1(i, t, knots);
    }
    n = knots[i+2] - t;
    d = knots[i+2] - knots[i+1];
    if (d != 0.) {
          sum += (n / d) * b1(i+1, t, knots);
    }
    return sum;
}
float b3(int i, float t, float *knots)
{
    float   n, d;
    float   sum = 0.;

    n = t - knots[i];
    d = knots[i+2] - knots[i];
    if (d != 0.) {
          sum += (n / d) * b2(i, t, knots);
    }
    n = knots[i+3] - t;
    d = knots[i+3] - knots[i+1];
    if (d != 0.) {
          sum += (n / d) * b2(i+1, t, knots);
    }
    return sum;
}
float b4(int i, float t, float *knots)
{
    float   n, d;
    float   sum = 0.;

    n = t - knots[i];
    d = knots[i+3] - knots[i];
    if (d != 0.) {
          sum += (n / d) * b3(i, t, knots);
    }
    n = knots[i+4] - t;
    d = knots[i+4] - knots[i+1];
    if (d != 0.) {
          sum += (n / d) * b3(i+1, t, knots);
    }
    return sum;
}
// constructor for NUB curves
Nub_curve::Nub_curve(void)
{
    num_knots = 0;
    knots = (float *)NULL;
}
// destructor for curve
Nub_curve::~Nub_curve(void)
{
    if (num_knots != 0 || knots) {
    free(knots);
    }
}
// setup the knot vector with user information
void Nub_curve::set_knot_vector(int count, float *v) {
    if (count != num_knots && knots) {
         free(knots);
    }
    knots = (float *)malloc(count * sizeof(float));
    if (knots) {
         int i;
         num_knots = count;
         for (i = 0; i < num_knots; i++)
             knots[i] = v[i];
    }
}
// display method for a NUB curve
void Nub_curve::display_curve(int n, int color, int show_geom_pts,
    int show_convex_hull)
{
    int     i, j, m;
    float   t;
    float   delta;
    float   x0, y0, x, y, x1, y1;
    float   bt0, bt1, bt2, bt3;
#define NUBEVAL(X, Y)           \
    bt0 = b4(i-3, t, knots);        \
    bt1 = b4(i-2, t, knots);        \
    bt2 = b4(i-1, t, knots);        \
    bt3 = b4(i, t, knots);      \
    X = geom[i-3].get_x() * bt0 +   \
    geom[i-2].get_x() * bt1 +   \
    geom[i-1].get_x() * bt2 +   \
    geom[i].get_x() * bt3;      \
    Y = geom[i-3].get_y() * bt0 +   \
    geom[i-2].get_y() * bt1 +   \
    geom[i-1].get_y() * bt2 +   \
    geom[i].get_y() * bt3;
    delta = (knots[4] - knots[3]) / (float)n;
    i = 3;
    t = knots[3];
    NUBEVAL(x0, y0);
    for (j = 2; j <= n; j++) {
          t += delta;
          NUBEVAL(x, y);
          line(
              (short)x0, (short)(floor(y0)+0.5),
              (short)x, (short)(floor(y)+0.5),
              color);
         x0 = x; y0 = y;
    }
    if (show_geom_pts) {
         for (i = 0; i < 4; i++) {
              x = geom[i].get_x();
              y = floor(geom[i].get_y())+0.5;
              line((short)x-CROSS_SIZE, (short)y, (short)x+CROSS_SIZE, 
                                                            (short)y, color);
              line((short)x, (short)y-CROSS_SIZE, (short)x, 
                                                 (short)y+CROSS_SIZE, color);
          }
    }
    if (show_convex_hull) {
          for (i = 0; i < 3; i++) {
             x0 = geom[i].get_x();
             y0 = floor(geom[i].get_y())+0.5;
             x1 = geom[i+1].get_x();
             y1 = floor(geom[i+1].get_y())+0.5;
             line((short)x0, (short)y0, (short)x1, (short)y1, color);
         }
    }
}
// display method for a NUB curve
void Nurb_curve::display_curve(int n, int color, int show_geom_pts,
    int show_convex_hull)
{
    int     i, j;
    float   t;
    float   delta;
    float   x0, y0, w, x, y, x1, y1;
    float   bt0, bt1, bt2, bt3;
#define NURBEVAL(X, Y, W)           \
    bt0 = b4(i-3, t, knots);            \
    bt1 = b4(i-2, t, knots);            \
    bt2 = b4(i-1, t, knots);            \
    bt3 = b4(i, t, knots);          \
    X = geom[i-3].get_x() * bt0 +       \
    geom[i-2].get_x() * bt1 +       \
    geom[i-1].get_x() * bt2 +       \
    geom[i].get_x() * bt3;          \
    Y = geom[i-3].get_y() * bt0 +       \
    geom[i-2].get_y() * bt1 +       \
    geom[i-1].get_y() * bt2 +       \
    geom[i].get_y() * bt3;          \
    W = geom[i-3].get_w() * bt0 +       \
    geom[i-2].get_w() * bt1 +       \
    geom[i-1].get_w() * bt2 +       \
    geom[i].get_w() * bt3;          \
    if (W != 1. && W != 0.) {           \
    X /= W; Y /= W;             \
    }
    delta = (knots[4] - knots[3]) / (float)n;
    i = 3;
    t = knots[3];
    NURBEVAL(x0, y0, w);
    for (j = 2; j <= n; j++) {
         t += delta;
         NURBEVAL(x, y, w);
         line(
             (short)x0, (short)(floor(y0)+0.5),
             (short)x, (short)(floor(y)+0.5),
             color);
         x0 = x; y0 = y;
    }
    if (show_geom_pts) {
         for (i = 0; i < 4; i++) {
             x = geom[i].get_x();
             y = floor(geom[i].get_y())+0.5;
             line((short)x-CROSS_SIZE, (short)y, (short)x+CROSS_SIZE, 
                                                             (short)y, color);
             line((short)x, (short)y-CROSS_SIZE, (short)x, 
                                                  (short)y+CROSS_SIZE, color);
         }
    }
    if (show_convex_hull) {
          for (i = 0; i < 3; i++) {
              x0 = geom[i].get_x();
              y0 = floor(geom[i].get_y())+0.5;
              x1 = geom[i+1].get_x();
              y1 = floor(geom[i+1].get_y())+0.5;
              line((short)x0, (short)y0, (short)x1, (short)y1, color);
         }
    }
}



[LISTING THREE]

// Horner's method for curve display method
void Basis_matrix_curve::display_curve(int n, int color, int 
show_geom_pts,
    int show_convex_hull)
{
    int     i;
    float   delta;
    float   t, t2, t3;
    float   x0, y0, x, y, x1, y1;

    x0 = t_coeff[3].get_x(); y0 = t_coeff[3].get_y();
    delta = 1. / (float)n;
    t = 0.;
    for (i = 0; i <= n; i++) {
         t += delta;
         t2 = t * t;
         t3 = t2 * t;
         x = t_coeff[0].get_x() * t3 +
             t_coeff[1].get_x() * t2 +
             t_coeff[2].get_x() * t +
             t_coeff[3].get_x();
         y = t_coeff[0].get_y() * t3 +
             t_coeff[1].get_y() * t2 +
             t_coeff[2].get_y() * t +
             t_coeff[3].get_y();
         line(
             (short)x0, (short)(floor(y0)+0.5),
             (short)x, (short)(floor(y)+0.5),
             color);
         x0 = x; y0 = y;
         }
    if (show_geom_pts) {
         for (i = 0; i < 4; i++) {
             x = (short)geom[i].get_x();
             y = (short)(floor(geom[i].get_y())+0.5);
             line(x-CROSS_SIZE, y, x+CROSS_SIZE, y, color);
             line(x, y-CROSS_SIZE, x, y+CROSS_SIZE, color);
         }
    }
    if (show_convex_hull) {
         for (i = 0; i < 3; i++) {
             x0 = geom[i].get_x();
             y0 = (floor(geom[i].get_y())+0.5);
             x1 = geom[i+1].get_x();
             y1 = (floor(geom[i+1].get_y())+0.5);
             line(x0, y0, x1, y1, color);
         }
    }
}



[LISTING FOUR]

// Forward differencing for curve display method
void Basis_matrix_curve::display_curve(int n, 
    int color, 
    int show_geom_pts,
    int show_convex_hull)
{
    int     i;
    float       d, delta, delta2, delta3;
    float       x0, y0, x, y, x1, y1;
    float       dx, d2x, d3x;
    float       dy, d2y, d3y;

    delta = 1. / (float)n;
    delta2 = delta * delta;
    delta3 = delta2 * delta;

    dx = t_coeff[0].get_x() * delta3 + 
    t_coeff[1].get_x() * delta2 +
    t_coeff[2].get_x() * delta;
    d2x = 6. * t_coeff[0].get_x() * delta3 + 
    2. * t_coeff[1].get_x() * delta2;
    d3x = 6. * t_coeff[0].get_x() * delta3;
    dy = t_coeff[0].get_y() * delta3 + 
    t_coeff[1].get_y() * delta2 +
    t_coeff[2].get_y() * delta;
    d2y = 6. * t_coeff[0].get_y() * delta3 + 
    2. * t_coeff[1].get_y() * delta2;
    d3y = 6. * t_coeff[0].get_y() * delta3;

    x = x0 = t_coeff[3].get_x(); y = y0 = t_coeff[3].get_y();
    for (i = 1; i <= n; i++) {
         x += dx; dx += d2x; d2x += d3x;
         y += dy; dy += d2y; d2y += d3y;
         line(
             (short)x0, (short)(floor(y0)+0.5),
             (short)x, (short)(floor(y)+0.5),
             color);
         x0 = x; y0 = y;
         }
    if (show_geom_pts) {
         for (i = 0; i < 4; i++) {
             x = (short)geom[i].get_x();
             y = (short)(floor(geom[i].get_y())+0.5);
             line((short)x-CROSS_SIZE, (short)y,
             (short)x+CROSS_SIZE, (short)y, color);
             line((short)x, (short)y-CROSS_SIZE, 
             (short)x, (short)y+CROSS_SIZE, color);
         }
    }
    if (show_convex_hull) {
    for (i = 0; i < 3; i++) {
        x0 = geom[i].get_x();
        y0 = floor(geom[i].get_y())+0.5;
        x1 = geom[i+1].get_x();
        y1 = floor(geom[i+1].get_y())+0.5;
        line((short)x0, (short)y0, 
              (short)x1, (short)y1, color);
    }
    }
}


[LISTING FIVE]

// demo.cc - demonstrating different curves
// For Unix and X, compile with:

#include <iostream.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "curve.h"

extern "C" {
void    init_graphics_device(void);
void    close_graphics_device(void);
void    clear_window(void);
void    text_output(char *);
void    text_output_and_wait(char *);
};

main()
{
    float       knot[8];
    Point2d     geom[4];
    init_graphics_device();
    Hermite_curve   herm;
    text_output("Hermite Curves");

    geom[0].set_xy(20., 20.);
    geom[1].set_xy(300., 20.);
    geom[2].set_xy(0., 500.);
    geom[3].set_xy(0., -500.);
    herm.set_geom_vector(4, geom);
    herm.display_curve(32, 14, 0, 0);

    geom[0].set_xy(20., 20.);
    geom[1].set_xy(300., 20.);
    geom[2].set_xy(0., 200.);
    geom[3].set_xy(0., -200.);
    herm.set_geom_vector(4, geom);
    herm.display_curve(32, 13, 0, 0);

    geom[0].set_xy(20., 20.);
    geom[1].set_xy(300., 20.);
    geom[2].set_xy(0., 200.);
    geom[3].set_xy(-100., 0.);
    herm.set_geom_vector(4, geom);
    herm.display_curve(32, 7, 0, 0);

    geom[0].set_xy(20., 20.);
    geom[1].set_xy(300., 20.);
    geom[2].set_xy(200., 200.);
    geom[3].set_xy(-200., -200.);
    herm.set_geom_vector(4, geom);
    herm.display_curve(32, 6, 0, 0);

    text_output_and_wait("Press return:");
    Bezier_curve    bez;
    Nub_curve       nub;
    clear_window();
    text_output("Bezier and NUB Curves");

    geom[0].set_xy(20., 20.);
    geom[1].set_xy(50., 180.);
    geom[2].set_xy(300., 50.);
    geom[3].set_xy(100., 10.);
    bez.set_geom_vector(4, geom);
    bez.display_curve(32, 15, 1, 1);

    text_output_and_wait("Press return:");
    knot[0] = knot[1] = knot[2] = knot[3] = 0.;
    knot[4] = knot[5] = knot[6] = knot[7] = 1.;
    nub.set_knot_vector(8, knot);
    nub.set_geom_vector(4, geom);
    nub.display_curve(32, 6, 0, 0);

    knot[0] = knot[1] = knot[2] = 0.;
    knot[3] = 1.; knot[4] = 2.;
    knot[5] = knot[6] = knot[7] = 3.;
    nub.set_knot_vector(8, knot);
    nub.set_geom_vector(4, geom);
    nub.display_curve(32, 5, 0, 0);

    text_output_and_wait("Press return:");
    Catmull_Rom_curve   cm;
    clear_window();

    text_output("Catmull-Rom Curve");

    geom[0].set_xy(20., 20.);
    geom[1].set_xy(50., 180.);
    geom[2].set_xy(300., 50.);
    geom[3].set_xy(100., 10.);
    cm.set_geom_vector(4, geom);
    cm.display_curve(32, 15, 1, 1);

    text_output_and_wait("Press return:");
    Beta_spline_curve   beta;
    clear_window();
    text_output("Beta-splines");

    geom[0].set_xy(20., 20.);
    geom[1].set_xy(170., 180.);
    geom[2].set_xy(150., 180.);
    geom[3].set_xy(300., 20.);

    beta.set_geom_vector(4, geom);
    beta.display_curve(32, 7, 1, 1);

    text_output_and_wait("Press return:");
    Nurb_curve      nurb;
    clear_window();
    text_output("NURB Curves with varying Knots");

    geom[0].set_xy(20., 20.);
    geom[1].set_xy(50., 180.);
    geom[2].set_xy(300., 50.);
    geom[3].set_xy(100., 10.);
    knot[0] = knot[1] = knot[2] = knot[3] = 0.;
    knot[4] = knot[5] = knot[6] = knot[7] = 1.;
    nurb.set_knot_vector(8, knot);
    nurb.set_geom_vector(4, geom);
    nurb.display_curve(32, 15, 1, 1);

    knot[0] = knot[1] = knot[2] = 0.; knot[3] = 1.;
    knot[4] = 2.; knot[5] = knot[6] = knot[7] = 3.;
    nurb.set_knot_vector(8, knot);
    nurb.set_geom_vector(4, geom);
    nurb.display_curve(32, 14, 0, 0);

    knot[0] = 0.; knot[1] = 1.; knot[2] = 2.; knot[3] = 3.;
    knot[4] = 4.; knot[5] = 5.; knot[6] = 6.; knot[7] = 7.;
    nurb.set_knot_vector(8, knot);
    nurb.set_geom_vector(4, geom);
    nurb.display_curve(32, 13, 0, 0);

    knot[0] = 0.; knot[1] = 0.; knot[2] = 1.; knot[3] = 1.;
    knot[4] = 2.; knot[5] = 2.; knot[6] = 3.; knot[7] = 3.;
    nurb.set_knot_vector(8, knot);
    nurb.set_geom_vector(4, geom);
    nurb.display_curve(32, 12, 0, 0);
    text_output_and_wait("Press return:");
    close_graphics_device();

    return 0;
}



[LISTING SIX]

// utilTC.cc - display utilities for Turbo C
// For Borland C++, compile with:
// bcc -c -mh -P -v -Ie:\bc\include util.cc

#include <iostream.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <conio.h>
#include <alloc.h>
#include <dos.h>
#include <math.h>
#include "curve.h"

#define WIDTH       320
#define HEIGHT      200
#define ABS(x)      ((x) < 0) ? -(x) : (x)

void    init_graphics_device(void);
void    close_graphics_device(void);
void    set_pixel(short, short, short, short, short, short, short);
void    line(short, short, short, short, short);
void    clear_window(void);
void    text_output(char *);
void    text_output_and_wait(char *);
void    display_string(int, char *);

// Global variables associated with the VGA display
unsigned char far *vga_fb = (unsigned char far *)MK_FP(0xA000,0);
void init_graphics_device(void)
{
    union REGS  inregs, outregs;

    inregs.h.ah = 0;        // set video mode command
    inregs.h.al = 0x13;     // 320x200x8 mode
    int86(0x10, &inregs, &outregs);

    inregs.h.ah = 0xf;      // check video mode
    int86(0x10, &inregs, &outregs);
    if (inregs.h.al != 0x13) {
    printf("cannot get 320x200x8 VGA mode : got %x\n", 
inregs.h.al);
    exit(1);
    }
}
void close_graphics_device(void)
{
    union REGS  inregs, outregs;
    inregs.h.ah = 0;        // set video mode command
    inregs.h.al = 0x2;      // 80x25 mode
    int86(0x10, &inregs, &outregs);
}
void line(short x1, short y1, short x2, short y2, short color)
{
    short   x, y;
    short   deltax, deltay;
    short   temp;
    short   err;
    short   i;
    short   swap;
    short   s1, s2;

    x = x1; y = y1;
    deltax = (short)ABS(x2 - x1);
    deltay = (short)ABS(y2 - y1);

    if ((x2 - x1) < 0.) s1 = -1; else s1 = 1;
    if ((y2 - y1) < 0.) s2 = -1; else s2 = 1;
    if (deltay > deltax) {
         temp = deltax;
         deltax = deltay;
         deltay = temp;
         swap = 1;
    }
    else
         swap = 0;

    err = 2 * deltay - deltax;
    for (i = 1; i <= deltax; i++) {
         set_pixel(0, 0, WIDTH-1, HEIGHT-1, x, y, color);
         while (err >= 0) {
             if (swap) x += s1; else y += s2;
             err -= 2 * deltax;
         }
         if (swap) y += s2; else x += s1;
         err += 2 * deltay;
    }
}
void set_pixel(short xmin, short ymin, short xmax, short ymax,
    short x, short y, short color)
{
    unsigned char far   *pixel;
    unsigned long   offset;
    if (x >= xmin && x <= xmax && y >= ymin && y <= ymax) {
         offset = (unsigned long)(((unsigned long)WIDTH * (unsigned long)y) +
                                                            (unsigned long)x);
         pixel = (unsigned char far *)(unsigned long)vga_fb + 
                                                        (unsigned long)offset;
         *pixel = (unsigned char)color;
    }
}
void clear_window()
{
    unsigned long far   *pixel32;
    unsigned long   count;

    // Clear 4 8-bit pixels at once by writing a 32 bit value 
    // with a 32-bit plane mask.
    pixel32 = (unsigned long far *)vga_fb;
    // compute number of 8-bit pixels and then divide by 4 for words
    count = (unsigned long)WIDTH * (unsigned long)HEIGHT;
    count = count >> 2;
    while (count > 0) {
         *pixel32 = (unsigned long)0;
         pixel32++;
         count--;
    }
}
void text_output(char *str)
{
    display_string(23, str);
}
void text_output_and_wait(char *str)
{
    display_string(24, str);
    getch();
}
void display_string(int row, char *str)
{
    int     i;
    union REGS  inregs, outregs;
    for (i = 0; i < strlen(str); i++) {
         inregs.h.ah = 0x2;      // set cursor position
         inregs.h.bh = 0x0;      // display page
         inregs.h.dh = row;      // near bottom of display
         inregs.h.dl = i;        // column
         int86(0x10, &inregs, &outregs);

         inregs.h.ah = 0x9;      // write character
         inregs.h.bh = 0x0;      // display page
         inregs.h.bl = 0x7;      // black BG/white FG
         inregs.h.al = str[i];       // character
         inregs.x.cx = 1;        // repeat count
         int86(0x10, &inregs, &outregs);
    }
}


[LISTING SEVEN]

/* utilXlib.cc - utilities for Xlib */
#include <X11/Xlib.h>
#include <X11/X.h>
#include <math.h>

#define WIDTH       320
#define HEIGHT      220

void    init_graphics_device(void);
void    close_graphics_device(void);
void    line(short, short, short, short, short);
void    clear_window(void);
void    text_output(char *);
void    text_output_and_wait(char *);

Display *dpy;
Window  win;
GC  gc;

void init_graphics_device(void)
{
    XEvent  event;
    int     screen;
    int     done;

    /* Open the X display and get the screen */
    dpy = XOpenDisplay(NULL);
    if (!dpy) {
    printf("could not open display\n");
    exit(1);
    }
    screen = DefaultScreen(dpy);
    /* Create a simple window which will be a child to the root window of 
     * the display. */
    win = XCreateSimpleWindow(dpy, RootWindow(dpy, screen),
    0, 0, WIDTH, HEIGHT, 2, BlackPixel(dpy, screen),
    WhitePixel(dpy, screen));
    /* Display the window */
    XMapWindow(dpy, win);
    XSelectInput(dpy, win, ExposureMask | KeyPressMask | 
ButtonPressMask);
    done = 0;
    while(!done) {
    XNextEvent(dpy, &event);
    switch(event.type) {
        case Expose:
        done = 1;
        break;
    }
    }
    /* Create an X graphics context for rendering vectors that compromise 
     * curves. */
    gc = XCreateGC(dpy, win, 0, NULL);
    XSetForeground(dpy, gc, BlackPixel(dpy, screen));
}
void close_graphics_device(void)
{
    exit(0);
}
void line(short x1, short y1, short x2, short y2, short color)
{
    /* Draw a 2D line using X calls and then flush the pixels to the server. */
    XDrawLine(dpy, win, gc, x1, y1, x2, y2);
    XFlush(dpy);
}
void clear_window()
{
    /* Clear the X window to the background color. */
    XClearWindow(dpy, win);
    XFlush(dpy);
}
void text_output(char *str)
{
    /* Put a string into the window. */
    XDrawString(dpy, win, gc, 0, 192, str, strlen(str));
    XFlush(dpy);
}
void text_output_and_wait(char *str)
{
    /* Put a string into the window and wait for user input. */
    XDrawString(dpy, win, gc, 0, 210, str, strlen(str));
    XFlush(dpy);
    getchar();
}