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 Bspline methods. Today, Bézier and Bspline 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 Bspline curves, which is the foundation of surface and solid modeling in the subsequent chapters. When you understand what Bézier and Bspline 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.
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:

The explicit form is satisfactory when the function is singlevalued 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 x^{2}+y^{2}1 = 0. If we require an explicit equation of the same circle, then it must be divided into two segments, with y = +Ö[(r^{2}x^{2})] for the upper half and y = Ö[(r^{2}x^{2})] 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 nonlinear equation for each point and thus numerical procedures have to be employed. Furthermore, both explicit and implicit represented curves are axisdependent. The choice of coordinate system directly affects the ease of use.
Explicitly and implicitly defined curves are sometimes called nonparametric 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]:

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 vectorvalued parametric curve. Representing a parametric curve in the vectorvalued form allows a uniform treatment of two, three, or ndimensional space, and provides a simple yet highly expressive notation for ndimensional 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 vectorvalued 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 = (ua)/(ba) 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 Bspline curves. In particular, we are concerned with the representation of Bézier and Bspline curves as well as some essential geometric processing methods required for displaying and manipulating curves.
By basis conversion, a vectorvalued polynomial curve of degree n can be represented in the power basis as


r(0) = a_{0},
r(1) = a_{0} +a_{1}+a_{2}+a_{3},
r^{¢}(0) = a_{1},
r^{¢}(1) = a_{1}+2a_{2}+3a_{3},
for a_{0}, a_{1}, a_{2}, and a_{3}. 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 a_{i}.
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

Actually, it can be shown that the control vertices p_{i} are simply a rearrangement of the coefficients of Ferguson's form. The important consequences of the rearrangement are that
r(0) = p_{0},
r(1) = p_{3},
r^{¢}(0) = 3(p_{1}p_{0}),
r^{¢}(1) = 3(p_{3}p_{2}).
Therefore, the curve described by Bézier's form passes through the vertices p_{0} and p_{3}. Furthermore, the curve has a tangent at p_{0} in the direction from p_{0} to p_{1} and at p_{3} has a tangent in the direction from p_{2} to p_{3}. 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





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 (51), which in turn requires us to compute the Bézier basis functions B_{i,n}(u). Denoting the Bézier basis functions B_{i,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 B_{i,n}(u) = C_{n}^{i}(1u)^{ni}u^{i}. 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 p_{i} 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 p_{i}^{0} = p_{i}, we compute p_{i}^{j} at each iteration level as a linear combination of the vertices at the previous level, i.e.,


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.
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



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




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.
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 onedimensional 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 onedimensional 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[order1]); 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[i1].x, vertex[i1].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 (52)) is given by

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<orderi; 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 firstinfirstout stack (onedimensional 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 righthalved curves. We allocate actual memory for bez_lf and bez_rt only if the current level of recursion is less than level1. 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 (53) 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.
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 firstorder derivative at the joint, which is known as the firstorder parametric continuity. We may relax the constraint to require only the continuity of the tangent directions at the joint, which is known as the firstorder 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:
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. Bsplines (sometimes, interpreted as basis splines) were investigated by a number of researchers in the 1940s. But Bsplines did not gain popularity in industry until de Boor and Cox published their work in the early 1970s. Their recurrence formula to derive Bsplines is still the most useful tool for computer implementation.
It is beyond the scope of this section to discuss different ways of deriving Bsplines and their generic properties. Instead, we shall take you directly to the definition of a Bspline curve and then explain to you the mathematics of Bsplines. Given M control vertices (or de Boor points) d_{i} (i = 0,1,¼,M1), a Bspline curve of order k (or degree n = k1) is defined as


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


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 Î [u_{j},u_{j+1}), we compute only k basis
functions N_{jn,k}(u), N_{jn+1,k}(u), ¼, N_{j,k}(u) because the other basis functions vanish. Denoting the basis
functions by NFunct[], the recursive formula to compute the k nonvanishing
basis functions is
NFunct[0] = 1.0; for (i1=1; i1<k; i1++) { left[i1] = u  knots[j+1i1]; right[i1] = knots[i1+j]  u; saved = 0.0; for (i2=0; i2<i1; i2++) { temp = NFunct[i2] / (right[i2+1] + left[i1i2]); NFunct[i2] = saved + right[i2+1] * temp; saved = left[i1i2] * temp; } NFunct[i1] = saved; }
We note that the k nonvanishing basis functions are stored in NFunct[0], ... , NFunct[k1] to save the memory space. Accordingly, the evaluation of a point on the curve with respect to u Î [u_{j}, u_{j+1}) is

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 Bspline curve be defined by d_{i}^{0}, u_{j}^{0}, and N_{i,k}^{0}(u). Then, we insert a new parameter u in [u^{0}_{r},u^{0}_{r+1}) to form the new knot vector:
u^{1}_{j} = u^{0}_{j}, j = 0,1,¼,r
u^{1}_{r+1} = u
u^{1}_{j+1} = u^{0}_{j}, j = r+1,¼,M+k
which defines M+1 new basis functions N_{i,k}^{1}(u). Now we write the given Bspline curve in terms of the new control vertices d_{i}^{1} and the new basis functions, i.e.,

d^{1}_{i} = d^{0}_{i}, i = 0,¼,rk+1
d^{1}_{i+1} = d^{0}_{i}, i = r,¼,M1
To compute k1 affected poles d^{1}_{rk+2}, ¼,d^{1}_{r}, we set




We use a temporary array b_{i}^{j} to store the intermediate results. So,
b_{i}^{0} = d_{i}^{0}, i = rk+1,¼, r
do j = 1, p
do i = rk+1+j, r
a = (uu_{i}^{0})/(u_{i+kj}^{0}u_{i}^{0})
b_{i}^{j} = b_{i1}^{j1}+a(b_{i}^{j1}b_{i1}^{j1})
enddo
enddo
When completed, we have the following intermediate results:


It is noted that a big memory space is allocated to store the intermediate results b_{i}^{j}. Actually, by rearranging the indices we need only an k×k array to store the intermediate results. We modify step 2 as follows:
b_{i}^{0} = d_{rk+1+i}^{0}, i = 0,1,¼,k1
do j = 1, p
do i = j, k1
a = (uu_{rk+1+i}^{0})/(u_{r+1+ij}^{0}u_{rk+1+i}^{0})
b_{i}^{j} = b_{i1}^{j1}+a(b_{i}^{j1}b_{i1}^{j1})
enddo
enddo
The intermediate results are:

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, k1
a = (uu_{rk+1+i}^{0})/(u_{r+i}^{0}u_{rk+1+i}^{0})
b_{i} = b_{i1}+a(b_{i}b_{i1})
enddo
d_{rk+2} = b_{1}
d_{r+p1} = b_{k1}
// Loop for remaining levels of insertion:
do j = 2, p
p_{0} = b_{j1}
do i = j, k1
p_{1} = b_{i}
a = (uu_{rk+1+i}^{0})/(u_{r+1+ij}^{0}u_{rk+1+i}^{0})
b_{i} = p_{0}+a(p_{1}p_{0})
p_{0} = p_{1}
enddo
d_{rk+1+j} = b_{j}
d_{r+pj} = b_{k1}
enddo
// Copy remaining elements at the pth row:
d_{rk+1+i} = b_{i}, i = p+1,¼,k1
As it is seen, we separate the first level of insertion from others to avoid data copying. Furthermore, we use two variables p_{0} and p_{1} to store the data so that b_{i} can be overwritten at each level of insertion. Consequently, we need to allocate only one dimension array to store b_{i} (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 k1 number of u in the knot vector. Accordingly, we can split the curve at u. The other reason is to convert a Bspline 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 k1. Consequently, each span of the Bspline curve is a Bézier curve of the same degree. Therefore, a composite Bézier curve is a special case of a Bspline curve. Since inserting new knots also introduces additional control vertices, a composite Bézier curve usually requires more control vertices than a Bspline 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 u_{j}. Actually, the choice of knots depends on how we create a Bspline curve. For example, if a Bspline 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 Bspline 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 Bspline curves. Similar to a rational Bézier curve, a rational Bspline curve is used to represent a rational curve such as a circle and an ellipse precisely. The use of a rational Bspline 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 Bspline curve is represented by


So much for the math! Let us move to the next section to design an applet for displaying and manipulating Bspline curves.
Before creating the applet, we need to design the GBsplineCurve class that has the definition of a Bspline curve and member methods for data setting and geometric processing. In order to deal with planar, space, and rational Bspline curves easily, we use an onedimensional 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 hardcoded a Bspline cubic curve with seven control vertices (i.e., three nonvanishing spans) in the applet. When you understand the implementation detail discussed here, you may design a powerful applet that allows you to create Bspline 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 Bspline 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); } }
As usual, we start with mathematics. Given a set of points q_{i} (i = 0,1,¼,n) with associated parameter [`u]_{i}, we want to construct a Bspline cubic curve having n+1 control points, i.e.,




At the beginning of this section, we assumed that q_{i} and associated parameters [`u]_{i} as well as the knot sequence u_{j} are all known. In practice, however, the only known information is a set of points. The choice of [`u]_{i} and u_{j} are up to us to determine.
To associate parameter [`u]_{i} with q_{i} is known as a parametrization of q_{i}. Ideally, the parametrization of q_{i} should be chosen such that [`u]_{i} are proportional to the arc lengths of the curve measured from q_{0} to q_{i}. This is however not feasible since the curve r(u) is not defined yet. In applications, the following methods are used to parametrize q_{i}.



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 Bspline curve. As usual, we choose



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[i1] + Math.sqrt(dist); } // Compute associated parameters: dist = bar_u[npts  1]; for (i=1; i<npts1; 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_u1]; } for (j=1; j<npolesndeg; 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_row1; 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.0e15) 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_row1; 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 illconditioned 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<jspnndeg; 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 Bspline 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 Bspline 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.x2, pt_dc.y2, 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 Bspline 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.
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 Bspline 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.