Curve Fitting - Interpolation.

Copyright (c) Susan Laflin. August 1999.

Blending Functions

When the points to be interpolated, or the final shape of the profile are known, the interpolated formulae previously discussed give a good representation of the output. However in CAD, it is often the case that the curve must be modified interactively until it takes the desired form. In this case, it is usually most efficient to select `Blending Functions', f{i}(t), and use an expression of the form

F(t) = sum f{i}(t)*P{i}

where the `control points' P{i} = [ x{i},y{i},z{i} ] may be used to calculate the coordinates of points along the curve. Moving any of the control points will alter the shape of the curve and it is easy to write software which allows the user to modify the profile by dragging one or more control points to new positions. As the variable t goes from 0 to some tmax, the whole curve is generated. Examples of such functions include

Hat Functions.
Bezier Functions.
B-Spline Functions.

and these will be discussed in the following sections. Such functions must form a basis for the function space, so that any desired function may be expressed in the form given above.

Hat Functions

The simplest example of such a set of functions are the `Hat functions' which lead to piecewise linear interpolation over the data points. This may be considered a rather long-winded method of obtaining linear interpolation, but it does also introduce this method and is easier to follow than either the Bezier or the B-spline method. The graph for one of these functions is shown in the above figure and the name `Hat function' is simply taken from the shape of this graph. The value of the point P{i} will affect the profile from P{i-1} to P{i+1}, but outside this interval the Hat function is zero, indicating that P{i} has no effect on the shape of the profile.

The equation of the Hat function is given below:

H{i}(t) 	= (t-t{i-1})/(t{i}-t{i-1}) for t in [t{i-1}, t{i}]
           	= (t{i+1}-t)/(t{i+1}-t{i}) for t in [t{i}, t{i+1}]
         	=   0           for all other t
 
and the interpolation function for the whole curve takes the form:

		F(t) = sum (H{i}(t)*P{i} )

Bezier Curves

Reference: Rogers and Adams section 5.8.

Historically, the Bezier curves were developed first. They use the Bernstein polynomials as blending functions, which means that the equations have the form given below:

J{n,i}(t) = factorial{n}/({factorial{i}*factorial{n-i})* t^{i}*(1-t)^{n-i}

Remember that the definition of factorials specifies that factorial{0} = 1.

Thus in the case where n=3, we have: 
J{30}(t) = (1-t)^{3}		J{31}(t) = 3t(1-t)^{2} 
J{32}(t) = 3t^{2}(1-t)		J{33}(t) = t^{3} 
So for the n+1 control points, numbered P0,P1,P2 and P3, we have the function

F(t) = (1-t)^{3}*P0+3t(1-t)^{2}*P1+3t^{2}(1-t)*P2+t^{3}*P3

This expression may be evaluated for each coordinate independently. The diagram in the above figure gives the variation of each of the blending functions with respect to t. You will note that the resulting curve passes through P0 and P3. At all the other points of the curve, each control point has some effect on the shape of the curve, with P1 having its maximum effect at t=1/3 and P2 having its maximum effect at t=2/3.

This is referred to as `global control' since moving one of the intermediate control points will cause a slight difference in the overall shape of the curve. The other possibility is some form of `local control' where moving one control point will affect only a limited section of the curve, but we must move to the B-spline curves to achieve this.

Some examples of 4-point Bezier curves are given in the figure above. It is of course possible to approximate to a large outline by a series of Bezier curves. The diagram in the above figure shows an example of three successive 4-point Bezier curves used to give an outline. The first two curves are joined, since P3 for the first curve is equal to P0 for the second, but there is a discontinuity in the gradients. When we come to the second and third curves, P2 and P3 of the second curve and P0 and P1 of the third are collinear and this ensures a smooth join of both curve and gradient.

B-Spline Curves

Reference: Rogers and Adams sections 5.9 to 5.12.

When we come to consider a B-spline curve for n+1 points, we may start by noting the differences between Bezier and B-spline curves. Firstly, n+1 points in a Bezier curve will always give a curve of order n. For a k-point B-spline curve, k may be chosen from the set of integers 1,2,3,...,n+1. Thus we may choose from a whole family of curves, going from k=1 which is the polygon joining the control points up to k=n+1, which turns out to be the Bezier curve for those points. With a B-spline curve, the parameter t varies from 0 to tmax (where tmax = n-k+2). With a Bezier curve, tmax is always equal to 1.

To evaluate the coefficients of the B-spline curve for any given value of t, we need to use a `knot-vector'. We note that for a k-point curve, the curve will pass through any control point which is repeated at least k times. Since we require the curve to pass through the two end control points, the knot vector must take the following form.

First the value 0 is repeated k times to ensure that the curve starts at P0. Then the integers 1,2,3,...,tmax-1 are included. Finally the value tmax is included k times, to ensure that the curve ends at the control point Pn. Thus the knot vector contains n+k+1 integer values.

We also require a real array N(i,j) in which to evaluate the factors by which the control points are multiplied. The first column of N needs to be of length n+k and a k-point curve will require k columns. Since each column contains one element less than the previous one, the final kth column will contain n+1 numbers, which will multiply the n+1 control points. The number of non-zero values will not exceed k.

First we compare the value of t with the values in the knot vector, and calculate I such that

V(I) < t < V(I+1)

Then the first column of the array N is set so that N(I,1) = 1 and N(i,1)(t) = 0 for all other values of i. The remaining columns are each evaluated from the previous column according to the following formula:


N{i,j}(t)=	(t-V(i))*N{i,j-1}(t) +	(V(i+j)-t)*N{i+1,j-1}(t)
		--------------		---------------
		(V(i+j-1)-V(i))		(V(i+j)-V(i+1))

The use of this formula may be easier to understand if we look at a specific example. Let us consider the case when n=5, so we have six points called P0, P1, P2, P3, P4 and P5. For these control points, let us choose k=3, giving a third order spline curve. This tells us that tmax = n-k+2 = 4. As an example, let us evaluate this for t=0.5.

When I=2, we have V(I) < 0.5 < V(I+1) since V(2) = 0 and V(3) = 1.
The general equation takes the following form for the second column:

 N{i,2} = N{i,1}*(t-V(i))   + 	N{i+1,1}*(V(i+2)-t)
	  ----------------	-------------------
	  (V(i+1)-V(i))		(V(i+2)-V(i+1))

and for the third column:

N{i,3} =  N{i,2}*(t-V(i))  + 	N{i+1,2}*(V(i+3)-t)
	  --------------	-----------------
	  (V(i+2)-V(i))		(V(i+3)-V(i+1))
	

i   	V(i)   	N(i,1)	N(i,2)	N(i,3)	P
0    	0       0        0      0.25	P0
1    	0       0        0.5    0.625	P1
2    	0       1        0.5    0.125	P2
3    	1       0        0      0       P3
4    	2       0        0      0       P4
5    	3       0        0      0       P5
6    	4       0        0
7    	4       0
8    	4

You should note that the method generates three non-zero values and consequently the nearest three control points determine the values along the curve.

N.U.R.Bs.

Reference: Rogers and Adams section 5.13.

`NURB' stands for `non-uniform rational B-spline' and may be explained as follows:

A polynomial B-spline, as discussed in the previous section, may be expressed in the form:

C(t) = sum (N{i,k}(t).Pi )

where the control-points Pi have coordinates in 3D space, the parameter t varies from 0 to tmax, and the N{i,k} are polynomials in t of order k (and degree k-1). Note that the points are now numbered from 1 to m instead of 0 to n.

The knot-vector vj is given by 
0 =t1 =t2 = ..... =tk < t{k+1} <= t{k+2} ....<= tm =t{m+1}=....=t{m+k} =tmax

If the tj in the centre section are equally spaced, i.e. if there is some positive real number d (in the previous section d=1) such that

t{l+1} - t{l} =d for all k <= l <= m

then tj is said to be a `uniform knot vector'. We may generalise this to a knot-vector over an interval (a,b) instead of (0,tmax) and non-integer values of d and still have a `uniform' B-spline. However, if the increments between successive values in the knot-vector are unequal, then we have a `non-uniform B-spline'. The algorithm used in the previous section may still be used to calculate the values of N{i,k} for any given t.

A `rational' B-spline is obtained when we move to general homogeneous coordinates and so instead of the 3 values (x,y,z) for each point, we have the four values (x,y,z,h). The fourth coordinate ,h, must always be greater than zero (i.e. positive) and when the coordinates are scaled so that h=1, then the values of x, y and z correspond to those of our usual 3D space. If we denote the control points in homogeneous coordinates by Ph{i}, then the curve in homogeneous coordinates is

Ch(t) = sum ( N{i,k}(t). Ph{i}

and when this is converted into the corresponding equation in 3D space, we get the `rational B-spline', which has the equation:

C(t) = sum ( N{i,k}(t) h{i} Pi)/ sum ( N{i,k}(t) h{i}

and this has an associated knot-vector vj where a <= vj <= b and the curve is defined for values of t in the interval (a,b). In the general case, the knot-vector need not be uniform. If all the h{i} are equal to 1, then the equation reduces to that of the polynomial spline, since sum (N{i,k}(t) )=1 . Many of the curves we use regularly can be expressed in this form. For example:

1. Straight line (or polyline). This is obtained by setting k=2 in the above expression.

2. Bezier Curve. This is obtained by setting k=m in the above expression.

3. Circle. There are several sets of control points which will give a circle, two of them lie on either the circumscribing equilateral triangle or on the circumscribing square (and probably on any regular polygon around the circle). Taking the case of the square, we have the control points shown below and the following values of h{i} and vj.
The values of h{i} are (1, 1, 2, 1, 1, 1, 2, 1, 1) and the knot-vector, vj is (0, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 4).