Chapter 6
Curve Modeling

6.1  Introduction

The term computer graphics was often associated with creation of realistic scenes and animated images. With the advent of both computers and applied mathematics more ambitious applications of computer graphics were sought in the last few decades. Today, one of the most important areas where computer graphics plays a central role is computer-aided design (CAD) and computer-aided manufacture (CAM). Both CAD and CAM are extensively used in a large number of areas, including aerospace, automotive engineering, marine engineering, civil engineering, and electronic engineering. The successful applications of computer graphics in engineering is largely due to the progress of computer aided geometric design (CAGD) which provides the mathematical basis for describing and processing geometric shapes and data. The term geometric modeling(which includes curve modeling, surface modeling, and solid modeling) is often used as the synonym of computer aided geometric design, although some authors argue that geometric modeling means the building up of computer representations of complex shapes from representations of similar components.

The most promising description method of geometric shape is the parametric curves and surfaces. The theory of parametric curves and surfaces are well understood in algebraic geometry and differential geometry. However, their possibilities and advantages for representing geometric entities in a CAD environment were not known until the late 1950s. The major breakthroughs in CAGD were undoubtedly the theory of Ferguson curves and patches, Coons patches, Bézier curves and surfaces, later combined with B-spline methods. Today, Bézier and B-spline representations of curves and surfaces have been the industrial standard.

In this chapter, we shall briefly review different curve representation methods and tell you why the parametric representation of curves is of mots importance. Then, we shall discuss the mathematics on Bézier and B-spline curves, which is the foundation of surface and solid modeling in the subsequent chapters. When you understand what Bézier and B-spline curves are, you will have a chance to create some applets that provide you with tools to display and manipulate these curves. These tools are usually common in most CAD systems.

6.2  Curve representation

In this section, we shall overview different curve representations, their advantages, and shortcomings. You do not have to read this section to understand subsequent sections in this chapter. However, if you are new in the area of geometric modeling, you are recommended to do so since it explains you why the parametric representation of curves, surfaces, and solids is now an established tool in computer graphics.

In the Cartesian plane, an explicit equation of a planar curve is given by y = f(x), where f(x) is a prescribed function of x. An explicit representation of curve enables us to directly compute y at any value of x. If we are asked to represent a straight line in the Cartesian coordinate by an explicit form, we will probably give the following equation:

y = kx + b,
provided that this line is not vertical to the x-axis. Otherwise, we have to represent vertical lines as x = c, where c is some constant value. Such problems inherent in explicit form is easily dealt with when solving a problem by hand. However, it is a nuisance when programming geometrical problems for a computer. Another drawback with respect to the use of explicit form is numerical stability. Referring to the above straight line, we note that the computation of y is numerically unstable if k goes to infinity, indicating the line is nearly vertical. In general, if a curve has nearly vertical tangents, we may expect overflow or rounding error problems when computing the function values. For these reasons, the use of explicit form in computer aided geometric design is very limited.

The explicit form is satisfactory when the function is single-valued and the curve has no vertical tangents. However, this precludes many curves of practical importance such as circles, ellipses and the other conic sections. An implicit equation of the form f(x,y) = 0 can avoid the difficulties of multiple values and vertical tangents inherent in the explicit form. For example, a unit circle with its center at the origin is given by x2+y2-1 = 0. If we require an explicit equation of the same circle, then it must be divided into two segments, with y = +[(r2-x2)] for the upper half and y = -[(r2-x2)] for the lower half. This kind of segmentation creates cases which are a nuisance in computer programs.

Although an implicit form of curve can overcome some limitations inherent in an explicit form, it does not enable us to compute points on a curve directly. Usually, implicitly defined curves require the solution of a non-linear equation for each point and thus numerical procedures have to be employed. Furthermore, both explicit and implicit represented curves are axis-dependent. The choice of coordinate system directly affects the ease of use.

Explicitly and implicitly defined curves are sometimes called non-parametric curves. An alternative way of describing curves is the parametric form, which uses an auxiliary parameter to represent the position of a point. For example, a unit circle with center at origin may be represented by an angle parameter u [0, 2p]:

x(u) = cosu,     y(u) = sinu.
Curves having parametric form are called parametric curves. A parametric curve that has a polynomial parameterization is called a polynomial curve, which is a standard in CAD systems for describing curves and surfaces. Since a point on a parametric curve is specified by a single value of parameter, the use of parametric techniques free us from dependence on any particular system of coordinates. Therefore, the parametric description of a curve enables coordinate transformations such as translation and rotation, required for graphical display, to be performed very simply. The parametric form also avoids problems which can arise in representing closed or multiple-valued curves and curves with vertical tangents in a fixed coordinate system. Furthermore, the parametric method lends itself to the piecewise description of curves and surfaces, which is a basic technique for the description of free form shapes. Due to these advantages, parametric curves are most commonly used in computer aided geometric design.

In Cartesian space, a point is defined by distances from the origin along the three mutually orthogonal axes x, y, and z. In vector algebra, a point is often defined by a position vector r, which is the displacement with the initial point at the origin. The path of a moving point is then described by the position vectors at successive values of the parameter, say u . Hence, the position vector r is a function of u, i.e., r = r(u). In the literature, r(u) is called the vector-valued parametric curve. Representing a parametric curve in the vector-valued form allows a uniform treatment of two-, three-, or n-dimensional space, and provides a simple yet highly expressive notation for n-dimensional problems. Thus, its use allows researchers and programmers to use simple, concise, and elegant equations to formalize their algorithms before expressing them explicitly in the Cartesian space. For these reasons, a vector-valued parametric form is used intensively for describing geometric shapes in computer graphics and computer aided geometric design.

It should be noted that the curve r(u) is said to have an infinite length if the parameter is not restricted in any specific interval, i.e., u (-,+). Conversely, the curve r(u) is said to have a finite length if u is within a closed interval, for example, u [a,b] where a, b . Given a finite parametric interval [a,b], a simple reparametrization t = (u-a)/(b-a) would normalize the parametric interval to t [0,1].

A comprehensive study on curves, including polynomial parametric curves, is beyond the scope of this chapter. Interested readers may find such information in many text books on algebraic curves. In this chapter, we discuss only two types of parametric curves, namely Bézier curves and B-spline curves. In particular, we are concerned with the representation of Bézier and B-spline curves as well as some essential geometric processing methods required for displaying and manipulating curves.

6.3  Mathematical background on Bézier Curves

As indicated by the title, this section is about the mathematics of curve, in particular, the mathematics of Bézier curves. If you are new in this area, you may find the contents of this section are quite difficult to understand. But do not be discouraged. We had the same situation as well when we were involved in this area a decade ago. It may prove useful if you read through this section, then move to the next section to design some Java applets for displaying and manipulating Bézier curves. After you have a good ``feeling'' on how the curve is displayed and related to its control polygon, you come back to read this section again. By then, I am sure that you will understand the content better. Let's study the mathematics now.

By basis conversion, a vector-valued polynomial curve of degree n can be represented in the power basis as

r(u) = n

i = 0 
ai ui,    u [0,1].
If the Cartesian coordinates of ai are given and denoted by xi, yi, and zi, we may write ai = (xi,yi,zi) and compute the Cartesian point on the curve with respect to the given u explicitly as
x(u) = n

i = 0 
xi ui,    y(u) = n

i = 0 
yi ui,    z(u) = n

i = 0 
zi ui.
When n = 3, r(u) describes a parametric cubic curve which was used by Ferguson in the early 1960s for the design of aircraft in Boeing. The usual way to design a Ferguson parametric cubic curve is to specify the positions and the first derivatives at both end points (i.e., t = 0 and t = 1) and then solve the following equations:

r(0) = a0,

r(1) = a0 +a1+a2+a3,

r(0) = a1,

r(1) = a1+2a2+3a3,

for a0, a1, a2, and a3. When the coefficients of the Ferguson curve is defined, it is very simple and efficient to evaluate points and derivatives on the curve with respect to given u since the standard Horner method can be employed. However, for interactive design of geometric shape, representing a curve in the power basis is not quite suitable as designers can hardly predict the shape of the curve while changing the coefficients ai.

P.E. Bézier of the French car company Renault introduced in 1960s a new form to represent a parametric cubic curve, which is best known as a Bézier cubic curve. In Bézier's cubic form, we write

r(u) = (1-u)3p0+3u(1-u)2p1+3u2(1-u)p2+u3p3,     u [0,1]
where, pi are called the control vertices which form the control polygon of the curve (see Figure 6.1). In some software venders, the control vertices are also called the control poles.

Figure 6.1: A Bézier cubic curve and its control polygon.

Actually, it can be shown that the control vertices pi are simply a re-arrangement of the coefficients of Ferguson's form. The important consequences of the re-arrangement are that

r(0) = p0,

r(1) = p3,

r(0) = 3(p1-p0),

r(1) = 3(p3-p2).

Therefore, the curve described by Bézier's form passes through the vertices p0 and p3. Furthermore, the curve has a tangent at p0 in the direction from p0 to p1 and at p3 has a tangent in the direction from p2 to p3. These properties result in the striking resemblance between the shape of Bézier cubic curve and the shape of its associated control polygon. Then, what is good about this resemblance? The advantage is that we can anticipate the shape of the curve while manipulating the control vertices that form the control polygon. This is obviously important for interactive design of geometric curves. For example, in order to design a cubic curve, we may proceed as follows:

Later in this section, we shall design three applets that give us an intuitive appreciation of this geometric interpretation.

We now generalize a Bézier cubic curve to a Bézier curve of degree n. In the literature, a Bézier curve of degree n (or order k = n+1) is often written as

r(u) = n

i = 0 
piBi,n(u),    u [0,1]       (5-1)
where pi are the control vertices and Bi,n(u) are the Bézier basis functions of the following form:
Bi,n(u) = Cni(1-u)n-iui,
noting that Cni are the binomial coefficients, i.e.,
Cni = n!
(n-i)!i!
.
If you are a mathematically-oriented person, you may note that the Bézier basis functions are actually the Bernstein polynomials used by Bernstein to prove the Weierstrass theorem. For u = 0, we have B0,n(0) = 1 and Bi,n(0) = 0 (i = 1,2,, n) so that r(0) = p0. Similarly, we can derive that r(1) = pn. Therefore, a Bézier curve interpolates the control vertices p0 and pn at u = 0 and t = 1 respectively. If we differentiate r(u) with respect to u, we have
r(u) = n n-1

i = 0 
(pi+1-pi)Bi,n-1(u),
which is a Bézier curve of degree n-1 with its control vertices being equal to ( pi+1-pi). It is readily seen that
r(0) = n(p1-p0),     r(1) = n(pn-p0)
indicating that the curve r(u) has a tangent at p0 in the direction from p0 to p1 and at pn has a tangent in the direction from pn-1 to pn. Therefore, the resemblance between the shape of a generic Bézier curve and the shape of its associated control polygon is still preserved although it becomes weaker when the degree of curve increases. Another important property we have to mention is that the Bézier basis functions sum identically to one, i.e., i = 0n Bi,n(u) = 1. This means that the curve lies within the convex hull formed by the control vertices pi (Figure 6.2). The convex hull property of a Bézier curve is an important criterion in solving many geometric problems such as curve-to-curve and curve-to-surface intersections by a subdivision technique.

Figure 6.2: A Bézier curve lies within the convex hull.

Before we move to the next section to create some Bézier curve applets, we need to discuss how to evaluate a Bézier curve for rendering. A simple but an inefficient evaluation method is to compute the point with respect to the given u using equation (5-1), which in turn requires us to compute the Bézier basis functions Bi,n(u). Denoting the Bézier basis functions Bi,n(u) by BFunct[i], we have the following algorithm to compute all the basis functions:

BFunct[0] = 1.0;
par = 1.0 - u;
for (i=1; i<=n; i++)
{
   saved = 0.0;
   for (k=0; k<i; k++)
   {
      temp = BFunct[k];
      BFunct[k] = saved + par * BFunct[k];
      saved = u * temp;
   }
   BFunct[i] = saved;
}

The above algorithm is more efficient than that of computing individually the basis functions Bi,n(u) = Cni(1-u)n-iui. When all the basis functions are obtained, it is straight forward to evaluate the point on the curve with respect to u. Let x[i], y[i], and z[i] denote the Cartesian coordinates of pi and rx, ry, rz the three coordinates of r(u). Then, the rx, ry, rz are computed as follows:

 
rx = ry = rz =0;
for (i=0; i<=n; i++)
{
   rx += x[i] * BFunct[i];
   ry += y[i] * BFunct[i];
   rz += z[i] * BFunct[i];
}

The above evaluation method is fine if we use it to evaluate some specific points on the curve. However, it may not be suitable for curve rendering since it is inefficient. For the purpose of curve rendering, the subdivision technique based on the de Casteljau algorithm is preferable (de Casteljau is a French engineer who coincidentally developed Bézier curves and used them for car body design at the French car company Citrön). Denoting the level of iteration by j and letting pi0 = pi, we compute pij at each iteration level as a linear combination of the vertices at the previous level, i.e.,

pij = (1- _
u
 
)pij-1+ _
u
 
pi+1j-1,    (j = 1,2,,n; i = j,j+1,,n).      (5-2)
Note that [`u] denote the given parameter where we want to evaluate the curve. After n levels of iteration, we obtain
j = 0
p00
p10
p20
p30
pn0
j = 1
p11
p21
p31
pn1
j = 2
p22
p32
pn2
·
·
j = n
pnn
The point pnn is the point on the curve with respect to the given parameter [`u]. Furthermore, the original Bézier curve is now subdivided into two at the position pnn. The left segment (i.e., the partial curve corresponding to u [0, [`u]]) is defined by the vertices p00,p11,,pnn and the right segment (i.e., the partial curve corresponding to u [[`u], 1]) is defined by the vertices pnn, pnn-1,, pn0. The method we just stated is known as the de Casteljau algorithm. Although evaluation of a point with respect to a given u using the de Casteljau algorithm is not faster than the one we discussed earlier, it is preferred in curve rendering using subdivision technique. Rather than evaluate points along a Bézier curve, the curve is approximated by a piecewise liner version obtained by subdividing the control polygon recursively. This gives a finer and finer approximation to the curve. The level of subdivision terminates when a linearity criterion is satisfied. A possible linearity criterion may be derived based on the convex hull property of Bézier curves. The subdivision technique using the de Casteljau algorithm has been proven an efficient method for rendering Bézier curves and also a powerful method in solving some important geometric problems.

The de Casteljau algorithm can also be explained geometrically: If we divide the edges of the control polygon in the ratio (1-[`u]) to [`u], connect the resulting points by straight lines, divide the new edges again in the same ratios, and repeat this process a total of n times, then the dividing point obtained in the last step is the point on the curve corresponding to [`u]. Figure 6.3 shows such a geometric interpretation for a Bézier cubic curve.

Figure 6.3: A geometric interpretation of the de Casteljau algorithm

The Bézier curves we have discussed so far are polynomial Bézier curves. In describing some geometric shapes such circles and ellipses, we cannot use polynomial Bézier curves to represent them precisely since circles and ellipses are rational curves. For example, a unit circle with its center at the origin can be represented as

x(t) = 1-t2
1+t2
,    y(t) = 2t
1+t2
,    t (-,+).
Therefore, a circle is a rational curve as it has a rational parametrization. In order to represent a rational curve precisely, a rational Bézier curve is introduced. It is given by
r(u) =
n

i = 0 
wipiBi,n(u)

n

i = 0 
wiBi,n(u)
,     u [0, 1]
where wi are called the weights. If all the weights are identically equal to, say, c, it is readily seen that the rational Bézier curve degenerates to a polynomial curve since
r(u) =
c n

i = 0 
piBi,n(u)

c n

i = 0 
Bi,n(u)
= n

i = 0 
piBi,n(u),
noting the basis functions sum to 1. Therefore, a rational Bézier curve is a generalization of polynomial Bézier curve. Beside the advantage that a rational Bézier curve can represent a circle or any other rational curve precisely, we also have additional flexibility to design the shape of curve because of the use of weights. However, it should be pointed out that negative weights are not recommended in practice as negative weights destroy the convex hull property of a Bézier curve.

Well, a rational Bézier curve is good. But how do we evaluate it? Generally speaking, we do not need additional methods to evaluate it if we extend our geometry scope to the projective geometry so as to represent the control vertices in homogeneous coordinates as

Pi = (Xi, Yi, Zi, Wi) = (wixi, wiyi, wizi, wi).
In this case, we write the curve in homogeneous coordinates as
R(u) = n

i = 0 
Pi Bi,n(u).
Accordingly, the evaluation and subdivision methods we discussed previously can all be applied directly to the above equation to compute (homogeneous) points on the curve or the new (homogeneous) vertices for subdivided curve segments. To obtain the corresponding points in the Cartesian coordinates, however, we need to divide the first three homogeneous coordinates by the fourth, i.e.,
x(u) = X(u)
W(u)
,    y(u) = Y(u)
W(u)
,    z(u) = Z(u)
W(u)
.
If you are also concerned with the derivatives of a rational curve, the following formula can be used to convert the derivatives in homogeneous coordinates to the Cartesian coordinates:
r(m)(u) =
R(m)(u)- m

i = 1 
Cmir(m-i)(u) W(i)(u)

W(u)
.

We now have some mathematical foundations about Bézier curves. Therefore, let us move to the next section to design some applets for displaying and manipulating Bézier curves.

6.4  Creating Bézier curve applets

In this section, we shall use the knowledge gained in the previous section to create three Java applets that allow you to design and manipulate Bézier curves. The designing and manipulating tools presented in this section are so common that you will find them in most CAD systems.

Before creating any Java applet, we first need to design a GBezierCurve class that has the definition of a Bézier curve of degree n and provides some member methods for data setting and geometric processing. In order to make our GBezierCurve flexible enough to deal with planar, space (i.e., 3D), polynomial and rational Bézier curves, we use an one-dimensional array to store the coordinates of the control vertices as

poles[]={x0,x1,...,xn, y0,y1,...,yn,...}

and use member variables ndim and order to indicate the dimension and rationality of the curve respectively. In the GBezierCurve, we use poles rather than vertices simply because the word is shorter. Furthermore, we use order instead of degree to simplify index expression in an array (Recalling that order=degree+1). The source code of the GBezierCurve class is given below:

The first SetPole method is used for planar curve to assign x and y to the indexed pole. The second SetPole method is used for either a planar rational curve or a space curve. The third SetPole method is used only for a space rational curve. It should be noted that the coordinates are homogeneous coordinates if the curve is rational. Since the coordinates of the control vertices are stored in a one-dimensional array, we need only two methods to get the control vertices: one for unweighted (i.e., the vertices in Cartesian coordinates) and the other for weighted. The EvalBases method evaluates the Bézier basis functions, which was discussed in the previous section. Finally, the EvalPoint method evaluates the Cartesian point on the curve with respect to the given parameter u. We save it under the name GBezierCurve.java in the directory GeomLib, which is the directory for our package GeomLib (see chapter one for detail). Let's compile it to obtain GBezierCurve.class in the same directory.

We now write the first applet ShowBezCurve1 that displays a Bézier cubic curve. The applet also allows you to use your mouse to drag any control vertex to the desired position for shape design. The source code is given below.

import java.awt.*;
import java.awt.event.*;
import java.applet.Applet;
import GeomLib.*;

public class ShowBezCurve1 extends Applet
implements MouseListener, MouseMotionListener
{
   int pole_id;
   boolean found;
   GPoint mos_pt_world = new GPoint();
   GBezierCurve bezcv;
   GService gs = new GService();

   // Initialization:
   public void init()
   {
		gs.SetOrigin(0, getSize().height);
      setBackground(Color.lightGray);
      found = false;

      // Register mouse listeners:
      addMouseListener(this);
      addMouseMotionListener(this);

      // Create curve object and set poles:
      bezcv = new GBezierCurve(4, 2);
      bezcv.SetPole(0, 20.0, 30.0);
      bezcv.SetPole(1, 50.0, 60.0);
      bezcv.SetPole(2, 85.0, 50.0);
      bezcv.SetPole(3, 110.0,25.0);
    }

   // Required to declare listener interfaces:
   public void mouseClicked(MouseEvent evt){}
   public void mouseEntered(MouseEvent evt){}
   public void mouseExited(MouseEvent evt){}
   public void mouseMoved(MouseEvent evt){}

   // Check if mouse position matches any control pole when mouse down:
   public void mousePressed(MouseEvent evt)
   {
      double dist, pole[]={0,0,0};

      for (int i=0; i< bezcv.order; i++)
      {
         bezcv.GetPole(i, pole);
         gs.DCToWorld(evt.getX(), evt.getY(), mos_pt_world);
         dist = Math.abs(pole[0]-mos_pt_world.x) + Math.abs(pole[1]-mos_pt_world.y);
         if (dist < 3.0)
         {
            pole_id = i;
            found = true;
            break;
         }
      }
   }

   // If foun=true, dynamically update selected pole to mouse position:
   public void mouseDragged(MouseEvent evt)
   {
      if (found == true)
      {
         gs.DCToWorld(evt.getX(), evt.getY(), mos_pt_world);
         bezcv.SetPole(pole_id, mos_pt_world.x, mos_pt_world.y);
         repaint();
      }
   }

   // When mouse button is released, initialize variable found:
   public void mouseReleased(MouseEvent evt)
   {
      found = false;
   }

   public void paint(Graphics g)
   {
      gs.DrawBezierCurve(g, bezcv);
   }
}

In the init() method, we register two mouse listeners and initialize the control vertices of the Bézier cubic curve. The mousePressed() method is used to edit the curve. When the mouse button is down, the mousePressed() method checks if the mouse position matches any control vertex of the curve. If so, the variable found is set to true. In turn, the mouseDragged() method dynamically updates the found control vertex to the current mouse position. Save this code under the name ShowBezCurve1.java in the directory Ch5 and compile it to obtain ShowBezCurve1.class. Then, write a web page called ShowBezCurve1.html to invoke the applet as shown in Figure 6.4. You can now do some experiments on the curve. In particular, you may want to prove the following resemblance between the shape of the curve and its control polygon:

The first applet allows you to design a Bézier cubic curve intuitively by manipulating the control vertices. However, it does not have the flexibility to design a Bézier curve of any degree. In the second applet, we shall show you how to create and display a generic Bézier curve. In particular, we shall show you how to use your mouse as an interactive tool to generate the control vertices. Let us name the second applet the ShowBezCurve2 and code it as follows.

import java.awt.*;
import java.awt.event.*;
import java.applet.*;
import GeomLib.*;

public class ShowBezCurve2 extends Applet
implements ActionListener,MouseListener,MouseMotionListener
{
   final int max_order = 16;
   int pole_id, order;
   GPoint mos_pt_world = new GPoint();
   GPoint vertex[] = new GPoint[max_order];
   boolean found, drawPolygon;
   GBezierCurve bezcv = new GBezierCurve(max_order, 2);
   GService gs = new GService();

   // Initialization:
   public void init()
   {
      setLayout(new FlowLayout());
      // Add two buttons to the applet:
      Button b = new Button("Draw Polygon");
      b.addActionListener(this);
      add(b);
      b = new Button("Draw Curve");
      b.addActionListener(this);
      add(b);

      // Register mouse listeners:
      addMouseListener(this);
      addMouseMotionListener(this);

      order = 0;
      found = false;
      drawPolygon = true;
		gs.SetOrigin(0, getSize().height);
      setBackground(Color.lightGray);
   }

   // Handle button events:
   public void actionPerformed(ActionEvent evt)
   {
      String buttonName = evt.getActionCommand();
      if (buttonName.equals("Draw Polygon"))
      {
         order = 0;
         drawPolygon = true;
      }
      else
      {
         // Since "Draw Curve" button is clicked, setup the curve:
         bezcv.order = order;
         for (int i=0; i< bezcv.order; i++)
            bezcv.SetPole(i, vertex[i].x, vertex[i].y);
         drawPolygon = false;
      }
      repaint();
   }

   // Required to declare listener interfaces:
   public void mouseClicked(MouseEvent evt){}
   public void mouseEntered(MouseEvent evt){}
   public void mouseExited(MouseEvent evt){}
   public void mouseMoved(MouseEvent evt){}

   // Either prepare data for polygon design or check if 
   // mouse position matches any existing control pole.
   public void mousePressed(MouseEvent evt)
   {
      if (drawPolygon)
      {
         if (order == 0)
         {
            vertex[order] = new GPoint();
            gs.DCToWorld(evt.getX(), evt.getY(), vertex[order++]);
         }
         if (order < max_order)
         {
            vertex[order] = new GPoint();
            gs.DCToWorld(evt.getX(), evt.getY(), vertex[order++]);
         }
      }
      else
      {
         double dist, p[]=new double[3];
         gs.DCToWorld(evt.getX(), evt.getY(), mos_pt_world);
         for (int i=0; i< bezcv.order; i++)
         {
            bezcv.GetPole(i, p);
            dist = Math.abs(p[0]-mos_pt_world.x) + Math.abs(p[1]-mos_pt_world.y);
            if (dist < 3.0)
            {
               pole_id = i;
               found = true;
               break;
            }
         }
      }
   }

   // Dynamically update control pole to current mouse position:
   public void mouseDragged(MouseEvent evt)
   {
      if (drawPolygon)
      {
         gs.DCToWorld(evt.getX(), evt.getY(), vertex[order-1]);
         repaint();
      }
      else if (found)
      {
         gs.DCToWorld(evt.getX(), evt.getY(), mos_pt_world);
         bezcv.SetPole(pole_id, mos_pt_world.x, mos_pt_world.y);
         repaint();
      }
   }

   // When mouse button is up, we initialize found=false.
   public void mouseReleased(MouseEvent evt)
   {
      found = false;
   }

   public void paint(Graphics g)
   {
      if (drawPolygon)
      {
         // Draw the polygon in the design stage:
         g.setColor(Color.red);
         for (int i=1; i< order; i++)
            gs.drawLine(g, vertex[i-1].x, vertex[i-1].y, vertex[i].x, vertex[i].y);
      }
      else
      {
         gs.DrawBezierCurve(g, bezcv);
      }
   }
}

In the init() method, we created and registered two buttons Draw polygon and Draw Curve to let you choose corresponding activity mode. Depending on which activity mode you have chosen, the actionPerformed() method either initializes the data or sets up the control vertices of a Bézier curve. If the current activity mode is Draw Polygon, the mouseReleased() method stores the current mouse position in the intermediate variable vertices[]. Otherwise, it sets the boolean variable found to false. The mousePressed() method is used to edit the curve. When the activity mode is Draw Curve and the mouse button is down, the mousePressed() method checks if the mouse position matches any control vertex of the curve. If so, the variable found is set to true. In turn, the mouseDragged() method dynamically updates the found control vertex to the current mouse position. Save the above code under the name ShowBezCurve2.java in the directory Ch5 and compile it to get ShowBezCurve2.class. Then, write a web page called ShowBezCurve2.html to invoke the applet as shown in Figure 6.5.

Our last example is to demonstrate the subdivision technique using the de Casteljau algorithm. In order to simplify the work, we assume that the curve is always subdivided at the middle point (i.e., u = 0.5). In this case, the de Casteljau algorithm (i.e., equation (5-2)) is given by

pij = 0.5(pij-1+pi+1j-1),    (j = 1,2,,n; i = j,j+1,,n)       (5-3)
noting that j denotes the level of iteration. If our goal is to subdivide the given curve once at the middle point, the work is fairly easy since only two subdivided curves are produced. On another hand, if we want to have flexibility to subdivide the curve to m level of recursion, the work is relatively involved as 2m subdivided curves will be produced. For example, 1024 subdivided curves will be generated if m = 10. Since so many curves are involved, the memory allocations and data copy have to be carefully arranged to reduce the overhead. In what follows, we shall add the subdivision method in the GBezierCurve class that can subdivide the curve to the given level of recursion. Although we do address the memory allocations and data copy issues in our implementation, the code is not fully optimized for the purpose of readability.

Now, let us open GBezierCurve.java in GeomLib directory and add in the following code fragment:

   // This is a wraper function exposed to the user:
   public void SubdivideCurve(int level, GBezierCurve bezcvs[])
   {
      int istk[] = {0};
      HalveCurve(level, 0, istk, bezcvs);
   }

   // Halving the curve recursively:
   private void HalveCurve(int level, int currentLevel, int istk[], 
                           GBezierCurve bezcvs[])
   {
      int  i, j, nd, n;
      GBezierCurve bez_lf, bez_rt;

      // Memory allocation:
      if (currentLevel + 1 == level)
      {
         bez_lf = bezcvs[istk[0]];
         istk[0] += 1;
         bez_rt = bezcvs[istk[0]];
         istk[0] += 1;
      }
      else
      {
         bez_lf = new GBezierCurve(order, ndim, rational);
         bez_rt = new GBezierCurve(order, ndim, rational);
      }

      // Copy the vertices of the current curve to the curve bez_rt
      n = (ndim + rational) * order;
      System.arraycopy(this.poles, 0, bez_rt.poles, 0, n);

      // Halving the curve by de Casteljau's algorithm:
      currentLevel++;
      for (nd=0; nd<ndim+rational; nd++)
      {
         n = nd * order;
         bez_lf.poles[n] = bez_rt.poles[n];
         for (i=1; i<order; i++)
         {
            for (j=0; j<order-i; j++)
               bez_rt.poles[n+j] = 0.5*(bez_rt.poles[n+j]+bez_rt.poles[n+j+1]);
            bez_lf.poles[n+i] = bez_rt.poles[n];
         }
      }

      // Check if halved curves need further halving:
      if (currentLevel == level)
      {
         return;
      }
      else
      {
         bez_lf.HalveCurve(level, currentLevel, istk, bezcvs);
         bez_rt.HalveCurve(level, currentLevel, istk, bezcvs);
      }
   }

As it is seen, the method HalveCurve() is declared as a private member method. Our intention is to hide the implementation detail from outsiders. Users call SubdivideCurve() to obtain subdivided curves, which has only the level of recursion as an input argument and the bezcvs[] as an output argument. The bezcvs[] is a first-in-first-out stack (one-dimensional array) that stores GBezierCurve. We use istk as the stack pointer that points to the top of the stack. Since istk is used as both an input and an output argument, it is declared as an array of size one. The bez_lf and bez_rt are used to store the left- and right-halved curves. We allocate actual memory for bez_lf and bez_rt only if the current level of recursion is less than level-1. Otherwise, we simply use them as pointers to the stack to reduce the number of memory allocations and data copy. The implementation of the de Casteljau algorithm is immediately after the memory arrangement block. Again, in order to save memory and data copy overhead, the algorithm given in (5-3) is slightly modified. The final block checks whether the current level equals the required one. If so, terminates the recursion. Otherwise, halve the left and right curves recursively.

As stated previously, the above code is not fully optimized for readability. If we introduce another stack pointer in the code, we could remove memory allocations in HalveCurve() completely. However, the code would be difficult to understand.

Save this code in GeomLib directory and compile it to update GBezierCurve.class. We now write the third applet to demonstrate the de Casteljau algorithm. In particular, we use a text field for you to input the parameter at which you want to evaluate or subdivide the curve. Let us name the applet ShowBezCurve3 and type the following code:

import java.awt.*;
import java.awt.event.*;
import java.applet.*;
import GeomLib.*;

public class ShowBezCurve3 extends Applet
implements ActionListener, MouseListener, MouseMotionListener
{
   int       level, numCurves, pole_id;
   boolean   found;
   TextField f;
   GPoint    mos_pt_world = new GPoint();
   GBezierCurve bezcv;
   GService gs = new GService();

   // Initialization:
   public void init()
   {
      // Add text field and register action listener:
      setLayout(new FlowLayout(FlowLayout.LEFT));
      add(new Label("Level of recursion (<10):"));
      f = new TextField(1);
      f.addActionListener(this);
      add(f);
      f.setText(String.valueOf(0));

      // Register mouse listeners:
      addMouseListener(this);
      addMouseMotionListener(this);

      // Create curve object and set control poles:
      bezcv = new GBezierCurve(4, 2);
      bezcv.SetPole(0, 20.0, 30.0);
      bezcv.SetPole(1, 50.0, 60.0);
      bezcv.SetPole(2, 85.0, 50.0);
      bezcv.SetPole(3, 110.0,25.0);

      found = false;
      level = 0;
      numCurves = 1;
      gs.SetOrigin(0, getSize().height);
      setBackground(Color.lightGray);
   }

   public void actionPerformed(ActionEvent evt)
   {
      // Convert text string to integer and check if input is valid:
      level = Integer.parseInt(f.getText());
      if (level > 9)
      {
         level = 9;
         f.setText(String.valueOf(level));
      }
      // Compute number of subdivided curves:
      numCurves = (int)Math.pow(2, level);
      repaint();
   }

   // Required to declare listener interfaces:
   public void mouseClicked(MouseEvent evt){}
   public void mouseEntered(MouseEvent evt){}
   public void mouseExited(MouseEvent evt){}
   public void mouseMoved(MouseEvent evt){}

   // Check if mouse position matches any control pole:
   public void mousePressed(MouseEvent evt)
   {
      double dist, pole[]=new double[3];
      gs.DCToWorld(evt.getX(), evt.getY(), mos_pt_world);
      for (int i=0; i< bezcv.order; i++)
      {
         bezcv.GetPole(i, pole);
         dist = Math.abs(pole[0]-mos_pt_world.x) + Math.abs(pole[1]-mos_pt_world.y);
         if (dist < 3.0)
         {
            pole_id = i;
            found = true;
            break;
         }
      }
   }

   // If found=true, update selected pole to mouse position:
   public void mouseDragged(MouseEvent evt)
   {
      if (found == true)
      {
         gs.DCToWorld(evt.getX(), evt.getY(), mos_pt_world);
         bezcv.SetPole(pole_id, mos_pt_world.x, mos_pt_world.y);
         repaint();
      }
   }

   // When mouse button is up, initialize found variable:
   public void mouseReleased(MouseEvent evt)
   {
      found = false;
   }

   public void paint(Graphics g)
   {
      int i, j;
      double p0[]=new double[3], p1[]=new double[3];

      // Draw the original curve:
      gs.DrawBezierCurve(g, bezcv);

      if (level > 0)
      {
         // Subdivide the curve with the level of recursion:
         GBezierCurve bezcvs[] = new GBezierCurve[numCurves];
         for (i=0; i< numCurves; i++)
            bezcvs[i] = new GBezierCurve(bezcv.order, bezcv.ndim);
         bezcv.SubdivideCurve(level, bezcvs);

         // Draw the control polygons of subdivided curves:
         g.setColor(Color.green);
         bezcvs[0].GetPole(0, p0);
         for (i=0; i< numCurves; i++)
         {
            for (j=1; j< bezcv.order; j++)
            {
               bezcvs[i].GetPole(j, p1);
               gs.drawLine(g, p0[0], p0[1], p1[0], p1[1]);
               p0[0] = p1[0]; p0[1] = p1[1];
            }
         }
      }
   }
}

The applet ShowBezCurve3 is a modified version of ShowBezCurve1. At the user interface level, we added a TextFiled to let you input the level of recursion, which is restricted to be less than 10. In the paint() method, we no longer evaluate the points along the curve. Instead, we plot only the polygon of the original and halved curves. Save this applet under the name ShowBezCurve3.java in the Ch5 directory and compile it. Then, write a web page called ShowBezCurve3.html to invoke the applet as shown in Figure 6.6.

6.5  Mathematical background on B-spline curves

In CAGD applications, a curve may have a so complicated shape that it cannot be represented by a single Bézier cubic curve since the shape of a cubic curve is not rich enough. Increasing the degree of a Bézier curve adds flexibility to the curve for shape design. However, this will significantly increase processing effort for curve evaluation and manipulation. Furthermore, a Bézier curve of high degree may cause numerical noise in computation. For these reasons, we often split the curve such that each subdivided segment can be represented by a lower degree Bézier curve. This technique is known as piecewise representation. A curve that is made of several Bézier curves is called a composite Bézier curve or a Bézier spline curve. In some area (e.g., computer data exchange), a composite Bézier cubic curve is known as the PolyBézier. If a composite Bézier curve of degree n has m Bézier curves, then the composite Bézier curve has in total m×n+1 control vertices.

A curve with complex shape may be represented by a composite Bézier curve formed by joining a number of Bézier curves with some constraints at the joints. The default constraint is that the curves are jointed smoothly. This in turn requires the continuity of the first-order derivative at the joint, which is known as the first-order parametric continuity. We may relax the constraint to require only the continuity of the tangent directions at the joint, which is known as the first-order geometric continuity. Increasing the order of continuity usually improves the smoothness of a composite Bézier curve. Although a composite Bézier curve may be used to describe a complex shape in CAGD applications, there are primarily two disadvantages associated the use of the composite Bézier curve:

  1. It is considerably involved to join Bézier curves with some order of derivatives continuity.
  2. For the reason that will become clear later, a composite Bézier curve requires more control vertices than a B-spline curve.

These disadvantages can be eliminated by working with spline curves. Originally, a spline curve was a draughtsman's aid. It was a thin elastic wooden or metal strip that was used to draw curves through certain fixed points (called nodes). The resulting curve minimizes the internal strain energy in the splines and hence is considered to be smooth. The mathematical equivalent is the cubic polynomial spline. However, conventional polynomial splines are not popular in CAD systems since they are not intuitive for iterative shape design. B-splines (sometimes, interpreted as basis splines) were investigated by a number of researchers in the 1940s. But B-splines did not gain popularity in industry until de Boor and Cox published their work in the early 1970s. Their recurrence formula to derive B-splines is still the most useful tool for computer implementation.

It is beyond the scope of this section to discuss different ways of deriving B-splines and their generic properties. Instead, we shall take you directly to the definition of a B-spline curve and then explain to you the mathematics of B-splines. Given M control vertices (or de Boor points) di (i = 0,1,,M-1), a B-spline curve of order k (or degree n = k-1) is defined as

r(u) = M-1

i = 0 
di Ni,k(u),
where, Ni,k(u) are polynomials of degree n and known as B-splines of order k. Analogous to Bézier curves, Ni,k(u) are also called the B-spline basis functions. In most practical applications, the B-spline basis functions are derived from the following knot vector (or knot sequence):
u0 = u1 = = un, uk uj uM-1, uM = uM+1 = = uM+n
where, ui are called knots. The reason to select the first and last k knots to be equal is that the control vertices d0 and dM-1 are the points on the curve. Furthermore, the tangent direction of the curve at d0 is from d0 to d1 and the tangent direction at dM-1 is from dM-2 to dM-1. Due to these properties, the shape of a B-spline curve resembles the shape of its control polygon formed by the control vertices di, although such resemblance is not as intuitive as that of a Bézier curve. In the literature, the M-k knots uk, uk+1, , uM-1 are sometimes called the interior knots. If the interior knots are all distinct, then the B-spline curve has M-k non-vanishing spans (i.e., the arc length of each span is not zero).

As we said previously, there are several methods to derive the B-spline basis functions Ni,k(u) in terms of the knot vector. We present only the recursive formula derived by de Boor and Cox as follows:

Ni,k(u) = u-ui
ui+k-1-ui
Ni,k-1(t)+ ui-u
ui+k-ui+1
Ni+1,k-1(u),
with
Ni,1 =

1, u [ui,ui+1)
0, elsewhere
These basis functions have the following properties:

The local support property is very important since it indicates that we do not need to compute all the basis functions for evaluation of a point on the curve. For the given u [uj,uj+1), we compute only k basis functions Nj-n,k(u), Nj-n+1,k(u), , Nj,k(u) because the other basis functions vanish. Denoting the basis functions by NFunct[], the recursive formula to compute the k non-vanishing basis functions is

NFunct[0] = 1.0;
for (i1=1; i1<k; i1++)
{
  left[i1] = u - knots[j+1-i1];
  right[i1] = knots[i1+j] - u;
  saved = 0.0;
  for (i2=0; i2<i1; i2++)
  {
     temp = NFunct[i2] / (right[i2+1] + left[i1-i2]);
     NFunct[i2] = saved + right[i2+1] * temp;
     saved = left[i1-i2] * temp;
  }
  NFunct[i1] = saved;
}

We note that the k non-vanishing basis functions are stored in NFunct[0], ... , NFunct[k-1] to save the memory space. Accordingly, the evaluation of a point on the curve with respect to u [uj, uj+1) is

r(u) = k-1

i = 0 
dj-n+i NFunct[i],
which is as efficient as the evaluation of a Bézier curve of the same order.

Another important technique we need to discuss is Boehm's knot insertion algorithm, which insets a knot into the existing knot vector. One reason to insert a knot into the knot vector is to gain additional flexibility to design the shape of curve. Let the control vertices, knots, and basis functions of the original B-spline curve be defined by di0, uj0, and Ni,k0(u). Then, we insert a new parameter u in [u0r,u0r+1) to form the new knot vector:

       u1j = u0j,     j = 0,1,,r

       u1r+1 = u

       u1j+1 = u0j,     j = r+1,,M+k

which defines M+1 new basis functions Ni,k1(u). Now we write the given B-spline curve in terms of the new control vertices di1 and the new basis functions, i.e.,

r(u) = M

i = 0 
N1i,k(u)d1i.
Due to the local support property of a B-spline curve, it is known that the basis functions N00,k(u),,N0r-k+1,k(u) and N0r,k(u), ,N0M-1,k(u) are not affected. Therefore, we can derive that

       d1i = d0i,     i = 0,,r-k+1

       d1i+1 = d0i,     i = r,,M-1

To compute k-1 affected poles d1r-k+2, ,d1r, we set

u = ur+11 = (1-ai)ui0+ai ui+k-10 = ui0+ai(ui+k-10-ui0).
Then, by the affine invariance, Boehm shows that
d1i = (1-ai)d0i-1+aid0i = di-10+ai(di0-di-10),     (i = r-k+2,,r)
with
ai = u - ui0
ui+k-10-ui0
.
Inserting a knot into the knot vector for multiple times can be done repeatedly using the method discussed above. However, this approach is very inefficient since we have to repeatedly update the control vertices and the knot vector after each knot insertion. Assume that we want to insert a knot u [ur0,ur+10) to the existing knot vector for p times. An improved method is summarized as follows:

  1. Copy the first part of unaffected poles:
    dpi = d0i,    i = 0,1, ,r-k+1.
  2. Compute k+p-2 affected poles dr-k+2p, , dr+p-1p:

    We use a temporary array bij to store the intermediate results. So,

        bi0 = di0,     i = r-k+1,, r

        do j = 1, p

              do i = r-k+1+j, r

                    a = (u-ui0)/(ui+k-j0-ui0)

                    bij = bi-1j-1+a(bij-1-bi-1j-1)

              enddo

        enddo

    When completed, we have the following intermediate results:

    j = 1
    b1r-k+2
    b1r-k+3
    b1r-k+4
    b1r
    j = 2
    b2r-k+3
    b2r-k+4
    b2r
    ·
    ·
    ·
    j = p
    bpr-k+p+1
    bpr
    What we need is (i) the diagonal elements br-k+21,br-k+32,,br-k+p+1p, (ii) the remaining elements at the pth row, i.e., br-k+p+2p,,brp, (iii) the remaining elements at the last column, i.e., brp-1, brp-2,,br1.
  3. Copy the last part of unaffected poles:
    dpp+i = d0i,     i = r,,m-1

It is noted that a big memory space is allocated to store the intermediate results bij. Actually, by rearranging the indices we need only an k×k array to store the intermediate results. We modify step 2 as follows:

    bi0 = dr-k+1+i0,     i = 0,1,,k-1

    do j = 1, p

          do i = j, k-1

                a = (u-ur-k+1+i0)/(ur+1+i-j0-ur-k+1+i0)

                bij = bi-1j-1+a(bij-1-bi-1j-1)

          enddo

    enddo

The intermediate results are:

j = 1
b11
b12
b13
b1k-1
j = 2
b22
b23
b2k-1
·
·
·
·
j = p
bpp
bpk-1
We then copy the intermediate results to the array of the control poles.

Actually, we can still reduce the memory size for intermediate results and improve the efficiency of code. We modify step 2 as follows:

    // First level of insertion:

    do i = 1, k-1

          a = (u-ur-k+1+i0)/(ur+i0-ur-k+1+i0)

          bi = bi-1+a(bi-bi-1)

    enddo

    dr-k+2 = b1

    dr+p-1 = bk-1

    // Loop for remaining levels of insertion:

    do j = 2, p

          p0 = bj-1

          do i = j, k-1

                p1 = bi

                a = (u-ur-k+1+i0)/(ur+1+i-j0-ur-k+1+i0)

                bi = p0+a(p1-p0)

                p0 = p1

          enddo

          dr-k+1+j = bj

          dr+p-j = bk-1

    enddo

    // Copy remaining elements at the pth row:

    dr-k+1+i = bi,    i = p+1,,k-1

As it is seen, we separate the first level of insertion from others to avoid data copying. Furthermore, we use two variables p0 and p1 to store the data so that bi can be overwritten at each level of insertion. Consequently, we need to allocate only one dimension array to store bi (i = 0,1,, k).

It is said that one reason to insert new knots into the knot vector is to gain flexibility to design the curve. We now give two more reasons: One is to subdivide the curve at the given parameter u. This is done by repeatedly inserting u until there is k-1 number of u in the knot vector. Accordingly, we can split the curve at u. The other reason is to convert a B-spline curve into a composite Bézier curve. This is achieved by repeatedly inserting knots into the knot vector until all the interior knots have the multiplicity equal to k-1. Consequently, each span of the B-spline curve is a Bézier curve of the same degree. Therefore, a composite Bézier curve is a special case of a B-spline curve. Since inserting new knots also introduces additional control vertices, a composite Bézier curve usually requires more control vertices than a B-spline curve does when they are used to represent the same geometric shape.

Up to now, we have not mentioned about how to choose the knots uj. Actually, the choice of knots depends on how we create a B-spline curve. For example, if a B-spline curve is created by first constructing the control polygon, then the curve that resembles the shape of the control polygon, we may choose the interior knots uniformly. On the other hand, if a B-spline curve is used to fit a set of points, we may then choose the knots as the ratio of the chordal length of two corresponding adjacent points to the accumulative chordal length of all the points (Detailed discussion will be given in the curve fitting section).

Finally, we briefly discuss about rational B-spline curves. Similar to a rational Bézier curve, a rational B-spline curve is used to represent a rational curve such as a circle and an ellipse precisely. The use of a rational B-spline curve also gives us more freedom to design the shape of the curve to meet our application requirements. Analogous to a rational Bézier curve, a rational B-spline curve is represented by

r(u) =
M-1

i = 0 
widi Ni,k(u)

M-1

i = 0 
wi Ni,k(u)
,
where, wi are the weights. Since the B-spline basis functions sum to one, a rational B-spline curve degenerates to a polynomial B-spline curve if all the weight are identically equal. If the knots of a rational B-spline curve is not uniformly spaced, the curve is then called the Non-Uniform Rational B-Spline curve or simply the NURBS curve. Let di = (xi, yi, zi). Then, representing the NURBS curve in homogeneous coordinates gives
R(u) = M-1

i = 0 
Ri Ni,k(u),
where, Ri = (wixi,wiyi,wizi,wi). Accordingly, we can use the evaluation and knot insertion techniques studied previously to evaluate or subdivide the curve.

So much for the math! Let us move to the next section to design an applet for displaying and manipulating B-spline curves.

6.6  Creating B-spline curve applet

In this section, we shall create an applet that allows you to manipulate the control vertices of a B-spline curve and hence the shape of the curve. It also lets you to insert a knot to the knot vector for a desired multiplicity.

Before creating the applet, we need to design the GBsplineCurve class that has the definition of a B-spline curve and member methods for data setting and geometric processing. In order to deal with planar, space, and rational B-spline curves easily, we use an one-dimensional array to store the control vertices (we call them poles here for shot). Analogous to GBezierCurve class, the GBsplineCurve class is defined as

In the GBsplineCurve class, we have two knots setting methods: One is to copy the user define knots, and the other is to set interior knots uniformly. The SetPole() methods are used to take care the need for setting planar, space, and rational curves. We also defined two methods GetWeightedPole and GetPole to access the indexed control vertex in the homogeneous and Cartesian coordinates respectively. The geometric processing methods FindSpan(), EvalBases(), EvalPoint(), and InsertKnots() are simply the Java implementation of our pseudo code given in the previous section. Let us save it under the name GBsplineCurve.java in the directory GeomLib and compile it to get GBsplineCurve.class.

We are now ready to create the applet. For simplicity, we hard-coded a B-spline cubic curve with seven control vertices (i.e., three non-vanishing spans) in the applet. When you understand the implementation detail discussed here, you may design a powerful applet that allows you to create B-spline curves in several ways as you have done with Bézier curve applets. In addition to the tool to manipulate the control vertices, we want to add two text fields to the applet. One is used for specifying the knot to be inserted and the other for the multiplicity. When you type in the desired value(s) in one or both text fields and hit the Enter key, the applet will respond to you by displaying the new control polygon and the curve with multiple knots inserted. Let's call this applet the ShowBspCurve and type the following code:

import java.awt.*;
import java.awt.event.*;
import java.applet.Applet;
import GeomLib.*;

public class ShowBspCurve extends Applet
implements ActionListener, MouseListener, MouseMotionListener
{
   int           mult, pole_id;
   double        knot_to_insert;
   boolean       found;
   TextField     knotField, multField;
   GBsplineCurve bspcv;
   GPoint mos_pt_world = new GPoint();
   GService gs = new GService();

   public void init()
   {
      knot_to_insert = 0.5;
      mult = 1;
      setLayout(new FlowLayout(FlowLayout.LEFT));
      add(new Label("The knot to insert:"));
      knotField = new TextField(2);
      knotField.addActionListener(this);
      add(knotField);
      knotField.setText(String.valueOf(knot_to_insert));

      add(new Label("Multiplicity:"));
      multField = new TextField(2);
      multField.addActionListener(this);
      add(multField);
      multField.setText(String.valueOf(mult));

      // Register mouse listeners:
      addMouseListener(this);
      addMouseMotionListener(this);
   
      // Initialize the B-spline curve:
      bspcv = new GBsplineCurve(4, 7, 2);
      bspcv.SetPole(0, 5.0, 15.0);
      bspcv.SetPole(1, 25.0, 5.0);
      bspcv.SetPole(2, 55.0, 10.0);
      bspcv.SetPole(3, 60.0, 55.0);
      bspcv.SetPole(4, 80.0, 35.0);
      bspcv.SetPole(5, 100.0, 55.0);
      bspcv.SetPole(6, 125.0, 10.0);
      bspcv.SetUniformKnots();

      found = false;
      setBackground(Color.lightGray);
      gs.SetOrigin(0, getSize().height);
   }

   public void actionPerformed(ActionEvent evt)
   {
      Double U;

      // Convert the string text to double and integer:
      U = Double.valueOf(knotField.getText());
      knot_to_insert = U.doubleValue();
      mult = Integer.parseInt(multField.getText());

      // Make sure user's input is valid:
      if (knot_to_insert <= 0.0 || knot_to_insert >= 1.0)
         knot_to_insert = 0.5;
      if (mult < 1 || mult > bspcv.order - 1)
         mult = 1;
      knotField.setText(String.valueOf(knot_to_insert));
      multField.setText(String.valueOf(mult));

      // Insert knot:
      bspcv.InsertKnots(knot_to_insert, mult);
      repaint();
   }

   // Required to declare listener interfaces:
   public void mouseClicked(MouseEvent evt){}
   public void mouseEntered(MouseEvent evt){}
   public void mouseExited(MouseEvent evt){}
   public void mouseMoved(MouseEvent evt){}

   // Check if the mouse position matches any control vertex:
   public void mousePressed(MouseEvent evt)
   {
      double dist, p[] = new double[3];
      gs.DCToWorld(evt.getX(), evt.getY(), mos_pt_world);
      for (int i=0; i< bspcv.numPoles; i++)
      {
         bspcv.GetPole(i, p);
         dist = Math.abs(p[0]-mos_pt_world.x) + Math.abs(p[1]-mos_pt_world.y);
         if (dist < 3.0)
         {
            pole_id = i;
            found = true;
            break;
         }
      }
   }

   // Update selected the control vertex position:
   public void mouseDragged(MouseEvent evt)
   {
      if (found)
      {
         gs.DCToWorld(evt.getX(), evt.getY(), mos_pt_world);
         bspcv.SetPole(pole_id, mos_pt_world.x, mos_pt_world.y);
         repaint();
      }
   }

   // Set found to false when the button is released:
   public void mouseReleased(MouseEvent evt)
   {
      found = false;
   }

   public void paint(Graphics g)
   {
      gs.DrawBsplineCurve(g, bspcv);
   }
}

We have explained the code with added comments. Save the source code under the name ShowBspCurve.java in the directory Ch5 and compile it to get ShowBspCurve.class. Then, write a HTML web page named ShowBspCurve.html to load the applet. If everything is fine, you should now see the applet shown in Figure 6.7.

6.7  Data fitting with B-spline curve

In this section, we shall take you to another exciting realm of free-form B-spline curve modeling: Representation of existing geometric shapes by B-spline curves. In practice, In particular, you will learn data fitting with B-spline curves, i.e., the construction of a B-spline curve that fits a given set of geometric data such points and derivative vectors. If the resulting curve is required to pass through every data point, then the fitting process is known as interpolation. If the curve does not necessarily passes through every point but approximates them, the fitting process is known as approximation. The most commonly used approximation method is the least squares approximation. Due to the scope of this book, we shall not go into detail about interpolation and approximation algorithms used in CAGD applications. Instead, we shall restrict ourselves to the interpolation by a B-spline cubic curve without derivative constraints.

As usual, we start with mathematics. Given a set of points qi (i = 0,1,,n) with associated parameter [`u]i, we want to construct a B-spline cubic curve having n+1 control points, i.e.,

r(u) = n

i = 0 
pi Ni,4(u)
such that
qi = r( _
u
 

i 
) = n

i = 0 
pi Ni,4( _
u
 

i 
),     (i = 0,1,,n).
It has been assumed that qi and their associated parameters [`u]i are given. We further assume that the knots sequence of a B-spline curve (i.e., uj with j = 0,1,,n+3) is also known. Then, the problem of constructing a B-spline curve becomes solving n+1 linear equations for n+1 unknowns pi. Writing this system of linear equations in matrix form gives








N0,4( _
u
 

0 
)
N1,4( _
u
 

0 
)
Nn,4( _
u
 

0 
)
N0,4( _
u
 

1 
)
N1,4( _
u
 

1 
)
Nn,4( _
u
 

1 
)
·
·
·
N0,4( _
u
 

n 
)
N1,4( _
u
 

n 
)
Nn,4( _
u
 

n 
)














p0
p1
·
pn






=





q0
q1
·
qn






.      (5-4)
Recalling from the local support property of B-splines, it is noted that there are at most four non-vanishing cubic basis functions with respect to each parameter [`u]i. Suppose that [`u]i [uj, uj+1), then the four non-vanishing basis functions are
Nj-3,4( _
u
 

i 
),    Nj-2,4( _
u
 

i 
),    Nj-1,4( _
u
 

i 
),    Nj,4( _
u
 

i 
).
Therefore, each row of the matrix in (5-4) has at most four non-vanishing elements. In other words, the matrix is a banded matrix with bandwidth 4. A linear system whose matrix has only a relatively small number of non-zero elements is called sparse. Although the method of solution to a sparse linear system is not different in principle from general methods of linear algebra, it is wasteful to use general methods because most of the arithmetic operations devoted to solving the equations involve zero operands. Furthermore, to reserve storage for zero elements wastes the memory space of your computer. For these reasons, special methods are used to solve sparse linear systems in real world programming. For simplicity, however, we shall use the Gaussian elimination method with partial pivoting to solve (5-4).

At the beginning of this section, we assumed that qi and associated parameters [`u]i as well as the knot sequence uj are all known. In practice, however, the only known information is a set of points. The choice of [`u]i and uj are up to us to determine.

To associate parameter [`u]i with qi is known as a parametrization of qi. Ideally, the parametrization of qi should be chosen such that [`u]i are proportional to the arc lengths of the curve measured from q0 to qi. This is however not feasible since the curve r(u) is not defined yet. In applications, the following methods are used to parametrize qi.

Method 1
: equally spaced parametrization, i.e.,
_
u
 

i 
= i
n
       (i = 0,1,, n).
This method is simple however may produce unwanted shapes (e.g., loops) when the data is unevenly spaced.
Method 2
: chordal parametrization. Let d = i = 1n||qi-qi-1|| and [`u]0 = 0. Then, we have
_
u
 

i 
= _
u
 

i-1 
+ ||qi-qi-1||
d
,     (i = 1,2,, n)
This is one of the most widely used methods since it gives an approximation to the arc length parametrization.

Method 3
: centripetal parametrization. Let [`u]i = 0 and d = i = 1n{||qi-qi-1||}. Then,
_
u
 

i 
= _
u
 

i-1 
+
  


||qi-qi-1||
 

d
,    (i = 1,2,, n)
This is a relatively new method that gives better results than the chordal parametrization when the data takes very sharp turns.

There are other parametrization methods as well. In this book, we use only the chordal parametrization.

We now consider how to determine the know sequence of the B-spline curve. As usual, we choose

u0 = u1 = u2 = u3   and    un = un+1 = un+2 = un+3
so that the starting and ending points of the B-spline curve will coincide with its control polygon. The simplest way to choose the remaining n-4 interior knots is again the equally spaced parametrization, i.e.,
u3+i = i
n-2
,    (i = 1,2,,n-3).
However, such simple parametrization does not take the distribution of {[`u]i} into consideration. A better parametrization is to take the average of [`u]i, i.e.,
u3+i = 1
3
( _
u
 

i 
+ _
u
 

i+1 
+ _
u
 

i+2 
),     (i = 1,2,,n-3).
This parametrization guarantees that every knot span contains at least one [`u]i.

Having enough information on mathematics, let's implement the algorithms one by one. We name method to compute chordal parametrization intuitively the ChordalParametrization and implement it as follows.

public void ChordalParametrization(int npts, int ndim, 
                        double Qpts[], double bar_u[])
{
   int     i, j, k, nd;
   double  temp, dist;

   bar_u[0] = 0.0;
   for (i=1; i<npts; i++)
   {
      j = i * ndim;
      k = j - ndim;
      dist = 0.0;
      for (nd=0; nd<ndim; nd++)
      {
         temp = Qpts[j++] - Qpts[k++];
         dist += temp * temp;
      }
      bar_u[i] = bar_u[i-1] + Math.sqrt(dist);
   }
   // Compute associated parameters:
   dist = bar_u[npts - 1];
   for (i=1; i<npts-1; i++)
      bar_u[i] /= dist;
   bar_u[npts - 1] = 1.0;
}

The method to compute the knot sequence by averaging [`u]i is named CreateKnotsByAverage and implemented as follows:

public void CreateKnotsByAverage(int ndeg, int npoles, int num_u, 
                                 double bar_u[], double knots[])
{
   int     i, j;
   double  temp;

   for (i=0; i<=ndeg; i++)
   {
       knots[i] = bar_u[0];
       knots[npoles + i] = bar_u[num_u-1];
   }
   for (j=1; j<npoles-ndeg; j++)
   {
      temp = 0.0;
      for (i=j; i<j+ndeg; i++)
         temp += bar_u[i];
      knots[ndeg+j] = temp / ndeg;
   }
}

The Gaussian elimination method can be found in most text books on numerical analysis. In our case, the following implementation is used:

public int GaussPPivot(int n_row, int n_col, double aug_mat[])
{  
   int     i, j, k, i_p=0;
   double  pivot;
   double  bidon;

   for (k=0; k<n_row-1; k++)
   {
      // Find pivot:
      pivot = 0.0;
      for (i=k; i<n_row; i++)
      {
         if (Math.abs(pivot) < Math.abs(aug_mat[i*n_col + k]))
         {
            pivot = aug_mat[i*n_col + k];
            i_p = i;
         }
      }
      if (Math.abs(pivot) < 1.0e-15)
         return 1; 

      // Check if row exchange is needed:
      if (i_p != k)
      {
         for (j=0; j<n_col; j++)
         {
            bidon = aug_mat[i_p*n_col + j];
            aug_mat[i_p*n_col + j] = aug_mat[k*n_col + j];
            aug_mat[k*n_col + j] = bidon;
         }
      }

      // Level k elimination:
      for (i=k+1; i<n_row; i++)
      {
         bidon = - aug_mat[i*n_col + k] / aug_mat[k*n_col + k];
         for (j=k; j<n_col; j++)
            aug_mat[i*n_col + j] += bidon * aug_mat[k*n_col + j];
      }

   }

   // Backward substitution:
   for (k=n_row-1; k>=0; k--)
   {
      for (i=n_row; i<n_col; i++)
      {
         for (j=k+1; j<n_row; j++)
            aug_mat[k*n_col+i] -= aug_mat[k*n_col+j] * aug_mat[j*n_col+i];
         aug_mat[k*n_col + i] /= aug_mat[k*n_col + k];
      }
   }
   return 0;
}

It is known that numerical solution of a linear system may fail for many reasons such as singular or ill-conditioned matrices. In our case the failure may occur if there are a number of coincident data points. In order to tell the calling function whether the system succeeds, we declare the return type of GaussPPivot as integer. If the system succeeds, the mthod return 0. Otherwise, it returns 1.

Finally, we implement the top level routine FitByBspCurve as follows.

public int FitByBspCurve(int npts, double Qpts[], GBsplineCurve bspcv)
{
   int     i, j, k, nd, n_row, n_col, jspn, ndeg, rc1;
   double  bar_u[], Amat[], NFunct[];

   ndeg = bspcv.order - 1;

   // Paramtrisation:
   bar_u = new double[npts];
   ChordalParametrization(npts, bspcv.ndim, Qpts, bar_u);
   CreateKnotsByAverage(ndeg, npts, npts, bar_u, bspcv.knots);

   // Compute the augument matrix:
   n_row = npts;
   n_col = n_row + bspcv.ndim;
   Amat = new double[n_row * n_col];
   NFunct = new double[bspcv.order];
   for (i=0; i<n_row; i++)
   {
      jspn = bspcv.FindSpan(bar_u[i]);
      bspcv.EvalBases(jspn, bar_u[i], NFunct);
      k = i * n_col;
      for (j=0; j<jspn-ndeg; j++)
         Amat[k++] = 0.0;
      for (j=0; j<bspcv.order; j++)
         Amat[k++] = NFunct[j];
      for (j=jspn+1; j<n_row; j++)
         Amat[k++] = 0.0;
      for (nd=0; nd<bspcv.ndim; nd++)
         Amat[k++] = Qpts[i*bspcv.ndim+nd];
   }

   // Solve the system:
   rc1 = GaussPPivot(n_row, n_col, Amat);
   if (rc1 == 0)
   {
      // Copy control vertices:
      for (i=0; i<n_row; i++)
      {
         k = i * n_col+n_row;
         bspcv.SetPole(i, Amat[k], Amat[k+1]);
      }
   }
   return rc1;
}

The four ChordalParametrization(), CreateKnotsByAverage(), GaussPPivot(), and FitByBspCurve() are all added to our GService class in the GeomLib directory.

We now design an applet ShowCurveFit to demonstrate data fitting with B-spline curves. In particular, we shall create two buttons, Generate Points and Fit with Curve, to let the user specify the task. In the process of generating points, the user clicks the mouse button to produce a number of discrete points on the applet's drawing area. These points are identified by drawing green circles. When this process is completed, the user may click the Fit with Curve button to create a B-spline curve that passes through all the points. The source code of this applet is listed below.

import java.awt.*;
import java.awt.event.*;
import java.applet.*;
import GeomLib.*;

public class ShowCurveFit extends Applet
implements ActionListener,MouseListener,MouseMotionListener
{
   final int order=4, ndim=2, max_npts = 100;
   int     error_code, npts, pole_id;
   boolean found, generatePoint;
   double  Qpts[] = new double[max_npts*ndim];
   GPoint  mos_pt_world = new GPoint();
   GBsplineCurve bspcv;
   GService gs = new GService();

   // Initialization:
   public void init()
   {
      setLayout(new FlowLayout());
      // Add two buttons to the applet:
      Button b = new Button("Generate Points");
      b.addActionListener(this);
      add(b);
      b = new Button("Fit with Curve");
      b.addActionListener(this);
      add(b);

      // Register mouse listeners:
      addMouseListener(this);
      addMouseMotionListener(this);

      npts = 0;
      error_code = 0;
      found = false;
      generatePoint = true;
      gs.SetOrigin(0, getSize().height);
      setBackground(Color.lightGray);
   }

   // Handle button events:
   public void actionPerformed(ActionEvent evt)
   {
      String buttonName = evt.getActionCommand();
      if (buttonName.equals("Generate Points"))
      {
         generatePoint = true;
         npts = 0;
      }
      else
      {
         generatePoint = false;
         bspcv = new GBsplineCurve(order, npts, 2);
         error_code = gs.FitByBspCurve(npts, Qpts, bspcv);
      }
      repaint();
   }

   // Required to declare listener interfaces:
   public void mouseClicked(MouseEvent evt){}
   public void mouseEntered(MouseEvent evt){}
   public void mouseExited(MouseEvent evt){}
   public void mouseMoved(MouseEvent evt){}

   // When mouse button is up, save the mouse position in design 
   // stage or set found variable to false in the editing stage.
   public void mouseReleased(MouseEvent evt)
   {
      if (generatePoint && npts < max_npts)
      {
         gs.DCToWorld(evt.getX(), evt.getY(), mos_pt_world);
         Qpts[npts*ndim] = mos_pt_world.x;
         Qpts[npts*ndim+1] = mos_pt_world.y;
         npts++;
         repaint();
      }
      else
      {
         found = false;
      }
   }

   // In editing stage, check if mouse position matches any control pole:
   public void mousePressed(MouseEvent evt)
   {
      if (!generatePoint)
      {
         double dist, p[]=new double[3];
         gs.DCToWorld(evt.getX(), evt.getY(), mos_pt_world);
         for (int i=0; i< bspcv.numPoles; i++)
         {
            bspcv.GetPole(i, p);
            dist = Math.abs(p[0]-mos_pt_world.x) + Math.abs(p[1]-mos_pt_world.y);
            if (dist < 3.0)
            {
               pole_id = i;
               found = true;
               break;
            }
         }
      }
   }

   // In editing stage, dynamically update selected pole to mouse position:
   public void mouseDragged(MouseEvent evt)
   {
      if (found && !generatePoint)
      {
         gs.DCToWorld(evt.getX(), evt.getY(), mos_pt_world);
         bspcv.SetPole(pole_id, mos_pt_world.x, mos_pt_world.y);
         repaint();
      }
   }

   public void paint(Graphics g)
   {
      Point pt_dc = new Point();
      // Draw green circles for data points:
      g.setColor(Color.green);
      for (int i=0; i< npts; i++)
      {
         gs.WorldToDC(Qpts[i*ndim], Qpts[i*ndim+1], pt_dc);
         g.drawOval(pt_dc.x-2, pt_dc.y-2, 4, 4);
      }
      if (!generatePoint)
      {
         if (error_code == 0)
            gs.DrawBsplineCurve(g, bspcv);
         else
            g.drawString("System failure. Please generate distinct points.", 10, 180);
      }
   }
}

In the above code, construction of a B-spline curve is done by calling a method FitByBspCurve() that is embedded in the GService class. As we mentioned earlier, a numerical solution of a linear system may fail if the system is singular. This could happen if the user clicks the mouse button repeatedly at the same potions, which in turn generates a number of coincident points. To safeguard such failures, we request the FitByBspCurve to return an error code so that we can tell the user about the failure.

Save the source code under the name ShowCurveFit.java in the directory Ch5 and compile it to get ShowCurveFit.class. Then, write a HTML web page to invoke the applet. If everything works well, you should be able to see an applet that is similar to Figure 6.8.

6.8  Summary and references

Curve modeling is the foundation of geometric modeling and plays a very important role in computer aided geometric design. Although curve modeling may sound simpler than modeling of surfaces and solids, its underlying mathematics is equivalently involved as the surface modeling does. Actually, you will find that the mathematics you have learned on curve modeling applies to surface modeling only differ in parametric dimensions. For this reason, we have put our effort in this chapter to discuss with you the essential mathematics required for designing, manipulating, and displaying free form curves (In particular, Bézier and B-spline curves). Although we try to make our topic straight forward, you may still find it difficult to understand. If this is true, please do not be discouraged since you are not alone. Many people studying NURBS (Non-Uniform Rational B-Splines) say that the acronym really means Nobody Understanding Rational B-Splines. Studying and implementing very example presented in this chapter may prove useful to understand the concepts and mathematics of curve modeling.

Due to the space limitation, we discussed in this chapter only limited algorithms in evaluation of curves. These algorithms will work in general cases. However, they may not be ideal ones. For example, it is considered inefficient to use knot insertion technique for evaluation of a B-spline curve with respect to a number of ordered parameter values. One of improved algorithms is to represent each span in the power basis by using Taylor expansion and then evaluate points and derivatives of that span by the Horner method. Interested readers may consult the references listed at the end of this section.

Finally, we want to point out that curve modeling is based on the theories of algebraic curve and geometry, differential geometry, and numerical analysis. To gain solid understanding of curve modeling, you are encouraged to explore these areas as well. We listed a few references for further readings.

  1. Clemens, C.H. (1980), A Scrapbook of Complex Curve Theory, Plenum Press.
  2. de Boor, C. (1978), A Practical Guide to B-spline, Springer-Verlag.
  3. Hildebrand, F.B. (1987), Introduction to Numerical Analysis, 2nd ed., Dover Pub.
  4. Hoschek, J. and Lasser, D. (1993), Fundamentals of Computer Aided Geometric Design, A K Peter, Ltd.
  5. Farin, G. (1988), Curves and Surfaces for Computer Aided Geometric Design , Academic Press.
  6. Piegl, L. and Tiller, W. (1997), The NURBS Book, 2nd ed., Springer-Verlag.
  7. Press, W. et al (1992), Numerical Recipes in C, 2nd ed., Cambridge University Press.
  8. Phillips, G.M and Taylor, P.J. (1973), Theory and Applications of Numerical Analysis, Academic Press.
  9. Struik, D.J. (1988), Lectures on Classical Differential Geometry, Dover Pub.
  10. Walker, R.J. (1978), Algebraic Curves, Springer-Verlag.