# Fundamentals of Linear Algebra and Optimization Project 1

The purpose of this project is to solve a curve interpolation problem using cubic splines.

This type of problem arises frequently in computer graphics and in robotics (path planning).

Recall from Homework 1 that polynomial curves of degree m can be dened in terms of control points and the Bezier polynomials. Cubic Bezier curves are often used because they are cheap to implement and give more exibility than quadratic Bezier curves.

A cubic Bezier curve C(t) (in R2 or R3) is specied by a list of four control points

(b0; b2; b2; b3) and is given parametrically by the equation

C(t) = (1 t)3 b0 + 3(1 t)2t b1 + 3(1 t)t2 b2 + t3 b3:

Clearly, C(0) = b0, C(1) = b3, and for t 2 [0; 1], the point C(t) belongs to the convex hull of the control points b0; b1; b2; b3.

Interpolation problems require nding curves passing through some given data points and possibly satisfying some extra constraints.

It is known that Lagrange interpolation is not very satisfactory when N is greater than 5, since Lagrange interpolants tend to oscillate in an undesirable manner. Thus, we turn to Bezier spline curves.

A Bezier spline curve F is a curve which is made up of curve segments which are Bezier curves, say C1; : : : ;Cm (m 2). We will assume that F dened on [0;m], so that for

i = 1; : : : ;m,

F(t) = Ci(t i + 1); i 1 t i:

Typically, some smoothness is required between any two junction points, that is, between any two points Ci(1) and Ci+1(0), for i = 1; : : : ;m 1. We require that Ci(1) = Ci+1(0) (C0-continuity), and typically that the derivatives of Ci at 1 and of Ci+1 at 0 agree up to second order derivatives. This is called C2-continuity, and it ensures that the tangents agree as well as the curvatures.

There are a number of interpolation problems, and we consider one of the most common problems which can be stated as follows:

1

x0 = d−1

x1

x2

x3

x4

x5

x6

x7 = d8

d0

d1

d2

d3

d4

d5

d6

d7

Figure 1: A C2 cubic interpolation spline curve passing through the points x0; x1, x2; x3,

x4; x5, x6; x7

Problem: Given N +1 data points x0; : : : ; xN, nd a C2 cubic spline curve F, such that F(i) = xi, for all i, 0 i N (N 2).

A way to solve this problem is to nd N +3 auxiliary points d1; : : : ; dN+1 called de Boor control points from which N Bezier curves can be found. Actually,

d1 = x0 and dN+1 = xN

so we only need to nd N + 1 points d0; : : : ; dN.

It turns out that the C2-continuity constraints on the N Bezier curves yield only N 1

equations, so d0 and dN can be chosen arbitrarily. In practice, d0 and dN are chosen according to various \end conditions," such as prescribed velocities at x0 and xN. For the time being, we will assume that d0 and dN are given.

Figure 1 illustrates an interpolation problem involving N + 1 = 7 + 1 = 8 data points.

The control points d0 and d7 were chosen arbitrarily.

It can be shown that d1; : : : ; dN1 are given by the linear system

2

0

[email protected]

7

2 1

1 4 1 0

. . . . . . . . .

0 1 4 1

1 7

2

1

CCCCCA

0

[email protected]

d1

d2

...

dN2

dN1

1

CCCCCA

=

0

[email protected]

6x1 3

2d0

6x2

...

6xN2

6xN1 3

2dN

1

CCCCCA

:

The derivation of the above system assumes that N 4. If N = 3, this system reduces

to

7

2 1

1 7

2

!

d1

d2

!

=

6x1 3

2d0

6x2 3

2d3

!

When N = 2, it can be shown that d1 is given by

d1 = 2x1

1

2

d0

1

2

d2:

Observe that when N 3, the above matrix is strictly diagonally dominant, so it is

invertible. Actually, the above system needs to be solved for the x-coordinates and for the y-coordinates of the dis (and also for the z-coordinates, if the points are in R3). Once the above system is solved, the Bezier cubics C1; : : : ;CN are determined as follows (we assume N 2): For 2 i N 1, the control points (bi

0; bi

1; bi

2; bi

3) of Ci are given by

bi

0 = xi1

bi

1 =

2

3

di1 +

1

3

di

bi

2 =

1

3

di1 +

2

3

di

bi

3 = xi:

The control points (b10

; b11

; b12

; b13

) of C1 are given by

b10

= x0

b11= d0

b12

=

1

2

d0 +

1

2

d1

b13

= x1;

and the control points (bN0; bN1

; bN2

; bN3

) of CN are given by

bN0

= xN1

bN1

=

1

2

dN1 +

1

2

dN

bN2

= dN

bN3

= xN:

3

Prove that the tangent vectors m0 at x0 and mN at xN are given by

m0 = 3(d0 x0)

mN = 3(xN dN):

End Conditions.

One method to determine the points d0 and dN is the natural end condition, which

consists in setting the second derivatives at x0 and at xN to be zero, that is,

C00 1 (0) = 0; C00 N(1) = 0:

(1) Prove that the second derivative at b0 of a Bezier cubic specied by the control points (b0; b1; b2; b3) is

6(b0 2b1 + b2); and the second derivative at b3 is

6(b1 2b2 + b3):

Prove that our system becomes

0

[email protected]

4 1

1 4 1 0

. . . . . . . . .

0 1 4 1

1 4

1

CCCCCA

0

[email protected]

d1

d2

...

dN2

dN1

1

CCCCCA

=

0

[email protected]

6x1 x0

6x2

...

6xN2

6xN1 xN

1

CCCCCA

;

where d0; dN are given by

d0 =

2

3

x0 +

1

3

d1

dN =

1

3

dN1 +

2

3

xN:

Note that d0 is on the line segment (x0; d1) (1=3 of the way from x0) and dN is on the line segment (dN1; xN) (1=3 of the way from xN).

In the above derivation we assumed that N 4. If N = 3, show that this system reduces

to

4 1

1 4

d1

d2

=

6x1 x0

6x2 x3

:

Another method to determine the points d0 and dN is the quadratic end condition, which consists in requiring that the second derivatives at x0 and at x1 agree, and similarly for the second derivatives at xN1 and xN; this means that

C00 1 (0) = C00 1 (1) and C00 N(0) = C00 N(1):

4

(2) Prove that our system becomes

0

[email protected]

5 1

1 4 1 0

. . . . . . . . .

0 1 4 1

1 5

1

CCCCCA

0

[email protected]

d1

d2

...

dN2

dN1

1

CCCCCA

=

0

[email protected]

7x1 x0

6x2

...

6xN2

7xN1 xN

1

CCCCCA

;

where d0; dN are given by

d0 = d1 +

2

3

x0

2

3

x1

dN = dN1 +

2

3

xN

2

3

xN1:

Geometrically, the vector d0d1 is equal to 2

3 (x0x1) and similarly the vector dN dN1

is equal to 2 3 (xN xN1); in particular, the line segments (d0; d1) and (x0; x1) are parallel, and so are (dN; dN1) and (xN; xN1).

In the above derivation we assumed that N 4. If N = 3, show that this system reduces

to

5 1

1 5

d1

d2

=

7x1 x0

7x2 x3

:

Yet another method known as Bessel end condition is to require the rst derivative C01 (0) at x0 to be equal to the rst derivative of the unique parabola passing through x0; x1; x2, for t = 0; 1; 2, and to require the rst derivative C0N (1) at xN to be equal to the rst derivative of the unique parabola passing through xN2; xN1; xN, for t = 0; 1; 2.

A parabola passing through x0 and x2 as above is given by

C(t) =

(2 t)2

4

x0 +

(2 t)t

2

b1 +

t2

4

x2;

so to require that C(1) = x1 means that

x1 =

1

4

x0 +

1

2

b1 +

1

4

x2;

which yields

b1 =

1

2

x0 + 2x1

1

2

x2:

(3) Show that

m0 = C0(0) = b1 x0 =

3

2

x0 + 2x1

1

2

x2;

5

and using the fact that d0 = x0 + 1

3m0, show that

d0 =

1

2

x0 +

1

2

x1 +

1

3

(x1 x2)

:

Geometrically, d0 is the midpoint of the line segment from x0 to a point obtained by extrap-

olation from x1 and x2 (x1 + 1

3 (x1 x2)).

Similarly, for the parabola interpolating xN2; xN1 and xN, we get

bN =

1

2

xN2 + 2xN1

1

2

xN;

and show that

mN = C0(1) = xN bN =

1

2

xN2 2xN1 +

3

2

xN;

and using the fact that dN = xN 1

3mN, show that

dN =

1

2

xN1 +

1

3

(xN1 xN2)

+

1

2

xN:

Geometrically, dN is the midpoint of the line segment from xN to a point obtained by extrapolation from xN1 and xN2 (xN1+ 1

3 (xN1xN2)). Note that the above derivation

is correct for N 2. Finally, there is the not a knot end condition, which consists in forcing the rst two Bezier segments C1 and C2 to belong to the same cubic curve, and similarly for the last two Bezier segments CN1 and CN. This amounts to require that C000 1 (1) = C000 2 (0) and C000 N1(1) = C000 N(0).

(4) Prove that the third derivative at b0 and at b3 of a Bezier cubic specied by the

control points (b0; b1; b2; b3) is

6(b0 + 3b1 3b2 + b3):

Prove that if N = 3, then

d0 =

7

18

x0 +

8

9

x1 +

7

18

x2

2

3

d2

d1 =

1

6

x0 +

4

3

x1

1

6

x2

d2 =

1

6

x1 +

4

3

x2

1

6

x3

d3 =

7

18

x1 +

8

9

x2 +

7

18

x3

2

3

d1

are already computed in terms of x0; : : : ; x3, and there is no need to solve any linear system.

6

Prove that if N = 4, then

d1 =

1

6

x0 +

4

3

x1

1

6

x2

d3 =

1

6

x2 +

4

3

x3

1

6

x4

d2 =

3

2

x2

1

4

d1

1

4

d3

d0 =

7

18

x0 +

8

9

x1 +

7

18

x2

2

3

d2

d4 =

7

18

x2 +

8

9

x3 +

7

18

x4

2

3

d2:

Prove that if N 5, our linear system becomes the (N 3) (N 3) system

0

[email protected]

4 1 0 0

1 4 1 0

. . . . . . . . .

0 1 4 1

0 0 1 4

1

CCCCCA

0

[email protected]

d2

d3

...

dN3

dN2

1

CCCCCA

=

0

[email protected]

6x2 + 1

6x0 4

3x1 + 1

6x2

6x3

...

6xN3

6xN2 + 1

6xN2 4

3xN1 + 1

6xN

1

CCCCCA;

and d0; d1; dN1; dN are given by

d0 =

7

18

x0 +

8

9

x1 +

7

18

x2

2

3

d2

d1 =

1

6

x0 +

4

3

x1

1

6

x2

dN1 =

1

6

xN2 +

4

3

xN1

1

6

xN

dN =

7

18

xN2 +

8

9

xN1 +

7

18

xN

2

3

dN2:

If N = 5, this system reduces to

4 1

1 4

!

d2

d3

!

=

6x2 + 1

6x0 4

3x1 + 1

6x2

6x3 + 1

6x3 4

3x4 + 1

6x5

!

:

(5) Implement the Gaussian elimination method with partial pivoting as well as the

method for solving a triangular system by back-substitution.

Use your program to solve several instances of the interpolation problem. Verify that no pivoting is needed.

(6) Implement the LU-factorization method for tridiagonal matrices and test it on the same interpolation problems as in (5).

Do you notice any improvement over Gaussian elimination (running time, numerical precision)?

(7) After computing d1; : : : ; dN1, compute the control points for the Bezier curves C1; : : :CN and write a program to plot these Bezier segments (for t 2 [0; 1]) to visualize the interpolating spline. Experiment with the choice of end conditions.

TOTAL: 300 points.

This type of problem arises frequently in computer graphics and in robotics (path planning).

Recall from Homework 1 that polynomial curves of degree m can be dened in terms of control points and the Bezier polynomials. Cubic Bezier curves are often used because they are cheap to implement and give more exibility than quadratic Bezier curves.

A cubic Bezier curve C(t) (in R2 or R3) is specied by a list of four control points

(b0; b2; b2; b3) and is given parametrically by the equation

C(t) = (1 t)3 b0 + 3(1 t)2t b1 + 3(1 t)t2 b2 + t3 b3:

Clearly, C(0) = b0, C(1) = b3, and for t 2 [0; 1], the point C(t) belongs to the convex hull of the control points b0; b1; b2; b3.

Interpolation problems require nding curves passing through some given data points and possibly satisfying some extra constraints.

It is known that Lagrange interpolation is not very satisfactory when N is greater than 5, since Lagrange interpolants tend to oscillate in an undesirable manner. Thus, we turn to Bezier spline curves.

A Bezier spline curve F is a curve which is made up of curve segments which are Bezier curves, say C1; : : : ;Cm (m 2). We will assume that F dened on [0;m], so that for

i = 1; : : : ;m,

F(t) = Ci(t i + 1); i 1 t i:

Typically, some smoothness is required between any two junction points, that is, between any two points Ci(1) and Ci+1(0), for i = 1; : : : ;m 1. We require that Ci(1) = Ci+1(0) (C0-continuity), and typically that the derivatives of Ci at 1 and of Ci+1 at 0 agree up to second order derivatives. This is called C2-continuity, and it ensures that the tangents agree as well as the curvatures.

There are a number of interpolation problems, and we consider one of the most common problems which can be stated as follows:

1

x0 = d−1

x1

x2

x3

x4

x5

x6

x7 = d8

d0

d1

d2

d3

d4

d5

d6

d7

Figure 1: A C2 cubic interpolation spline curve passing through the points x0; x1, x2; x3,

x4; x5, x6; x7

Problem: Given N +1 data points x0; : : : ; xN, nd a C2 cubic spline curve F, such that F(i) = xi, for all i, 0 i N (N 2).

A way to solve this problem is to nd N +3 auxiliary points d1; : : : ; dN+1 called de Boor control points from which N Bezier curves can be found. Actually,

d1 = x0 and dN+1 = xN

so we only need to nd N + 1 points d0; : : : ; dN.

It turns out that the C2-continuity constraints on the N Bezier curves yield only N 1

equations, so d0 and dN can be chosen arbitrarily. In practice, d0 and dN are chosen according to various \end conditions," such as prescribed velocities at x0 and xN. For the time being, we will assume that d0 and dN are given.

Figure 1 illustrates an interpolation problem involving N + 1 = 7 + 1 = 8 data points.

The control points d0 and d7 were chosen arbitrarily.

It can be shown that d1; : : : ; dN1 are given by the linear system

2

0

[email protected]

7

2 1

1 4 1 0

. . . . . . . . .

0 1 4 1

1 7

2

1

CCCCCA

0

[email protected]

d1

d2

...

dN2

dN1

1

CCCCCA

=

0

[email protected]

6x1 3

2d0

6x2

...

6xN2

6xN1 3

2dN

1

CCCCCA

:

The derivation of the above system assumes that N 4. If N = 3, this system reduces

to

7

2 1

1 7

2

!

d1

d2

!

=

6x1 3

2d0

6x2 3

2d3

!

When N = 2, it can be shown that d1 is given by

d1 = 2x1

1

2

d0

1

2

d2:

Observe that when N 3, the above matrix is strictly diagonally dominant, so it is

invertible. Actually, the above system needs to be solved for the x-coordinates and for the y-coordinates of the dis (and also for the z-coordinates, if the points are in R3). Once the above system is solved, the Bezier cubics C1; : : : ;CN are determined as follows (we assume N 2): For 2 i N 1, the control points (bi

0; bi

1; bi

2; bi

3) of Ci are given by

bi

0 = xi1

bi

1 =

2

3

di1 +

1

3

di

bi

2 =

1

3

di1 +

2

3

di

bi

3 = xi:

The control points (b10

; b11

; b12

; b13

) of C1 are given by

b10

= x0

b11= d0

b12

=

1

2

d0 +

1

2

d1

b13

= x1;

and the control points (bN0; bN1

; bN2

; bN3

) of CN are given by

bN0

= xN1

bN1

=

1

2

dN1 +

1

2

dN

bN2

= dN

bN3

= xN:

3

Prove that the tangent vectors m0 at x0 and mN at xN are given by

m0 = 3(d0 x0)

mN = 3(xN dN):

End Conditions.

One method to determine the points d0 and dN is the natural end condition, which

consists in setting the second derivatives at x0 and at xN to be zero, that is,

C00 1 (0) = 0; C00 N(1) = 0:

(1) Prove that the second derivative at b0 of a Bezier cubic specied by the control points (b0; b1; b2; b3) is

6(b0 2b1 + b2); and the second derivative at b3 is

6(b1 2b2 + b3):

Prove that our system becomes

0

[email protected]

4 1

1 4 1 0

. . . . . . . . .

0 1 4 1

1 4

1

CCCCCA

0

[email protected]

d1

d2

...

dN2

dN1

1

CCCCCA

=

0

[email protected]

6x1 x0

6x2

...

6xN2

6xN1 xN

1

CCCCCA

;

where d0; dN are given by

d0 =

2

3

x0 +

1

3

d1

dN =

1

3

dN1 +

2

3

xN:

Note that d0 is on the line segment (x0; d1) (1=3 of the way from x0) and dN is on the line segment (dN1; xN) (1=3 of the way from xN).

In the above derivation we assumed that N 4. If N = 3, show that this system reduces

to

4 1

1 4

d1

d2

=

6x1 x0

6x2 x3

:

Another method to determine the points d0 and dN is the quadratic end condition, which consists in requiring that the second derivatives at x0 and at x1 agree, and similarly for the second derivatives at xN1 and xN; this means that

C00 1 (0) = C00 1 (1) and C00 N(0) = C00 N(1):

4

(2) Prove that our system becomes

0

[email protected]

5 1

1 4 1 0

. . . . . . . . .

0 1 4 1

1 5

1

CCCCCA

0

[email protected]

d1

d2

...

dN2

dN1

1

CCCCCA

=

0

[email protected]

7x1 x0

6x2

...

6xN2

7xN1 xN

1

CCCCCA

;

where d0; dN are given by

d0 = d1 +

2

3

x0

2

3

x1

dN = dN1 +

2

3

xN

2

3

xN1:

Geometrically, the vector d0d1 is equal to 2

3 (x0x1) and similarly the vector dN dN1

is equal to 2 3 (xN xN1); in particular, the line segments (d0; d1) and (x0; x1) are parallel, and so are (dN; dN1) and (xN; xN1).

In the above derivation we assumed that N 4. If N = 3, show that this system reduces

to

5 1

1 5

d1

d2

=

7x1 x0

7x2 x3

:

Yet another method known as Bessel end condition is to require the rst derivative C01 (0) at x0 to be equal to the rst derivative of the unique parabola passing through x0; x1; x2, for t = 0; 1; 2, and to require the rst derivative C0N (1) at xN to be equal to the rst derivative of the unique parabola passing through xN2; xN1; xN, for t = 0; 1; 2.

A parabola passing through x0 and x2 as above is given by

C(t) =

(2 t)2

4

x0 +

(2 t)t

2

b1 +

t2

4

x2;

so to require that C(1) = x1 means that

x1 =

1

4

x0 +

1

2

b1 +

1

4

x2;

which yields

b1 =

1

2

x0 + 2x1

1

2

x2:

(3) Show that

m0 = C0(0) = b1 x0 =

3

2

x0 + 2x1

1

2

x2;

5

and using the fact that d0 = x0 + 1

3m0, show that

d0 =

1

2

x0 +

1

2

x1 +

1

3

(x1 x2)

:

Geometrically, d0 is the midpoint of the line segment from x0 to a point obtained by extrap-

olation from x1 and x2 (x1 + 1

3 (x1 x2)).

Similarly, for the parabola interpolating xN2; xN1 and xN, we get

bN =

1

2

xN2 + 2xN1

1

2

xN;

and show that

mN = C0(1) = xN bN =

1

2

xN2 2xN1 +

3

2

xN;

and using the fact that dN = xN 1

3mN, show that

dN =

1

2

xN1 +

1

3

(xN1 xN2)

+

1

2

xN:

Geometrically, dN is the midpoint of the line segment from xN to a point obtained by extrapolation from xN1 and xN2 (xN1+ 1

3 (xN1xN2)). Note that the above derivation

is correct for N 2. Finally, there is the not a knot end condition, which consists in forcing the rst two Bezier segments C1 and C2 to belong to the same cubic curve, and similarly for the last two Bezier segments CN1 and CN. This amounts to require that C000 1 (1) = C000 2 (0) and C000 N1(1) = C000 N(0).

(4) Prove that the third derivative at b0 and at b3 of a Bezier cubic specied by the

control points (b0; b1; b2; b3) is

6(b0 + 3b1 3b2 + b3):

Prove that if N = 3, then

d0 =

7

18

x0 +

8

9

x1 +

7

18

x2

2

3

d2

d1 =

1

6

x0 +

4

3

x1

1

6

x2

d2 =

1

6

x1 +

4

3

x2

1

6

x3

d3 =

7

18

x1 +

8

9

x2 +

7

18

x3

2

3

d1

are already computed in terms of x0; : : : ; x3, and there is no need to solve any linear system.

6

Prove that if N = 4, then

d1 =

1

6

x0 +

4

3

x1

1

6

x2

d3 =

1

6

x2 +

4

3

x3

1

6

x4

d2 =

3

2

x2

1

4

d1

1

4

d3

d0 =

7

18

x0 +

8

9

x1 +

7

18

x2

2

3

d2

d4 =

7

18

x2 +

8

9

x3 +

7

18

x4

2

3

d2:

Prove that if N 5, our linear system becomes the (N 3) (N 3) system

0

[email protected]

4 1 0 0

1 4 1 0

. . . . . . . . .

0 1 4 1

0 0 1 4

1

CCCCCA

0

[email protected]

d2

d3

...

dN3

dN2

1

CCCCCA

=

0

[email protected]

6x2 + 1

6x0 4

3x1 + 1

6x2

6x3

...

6xN3

6xN2 + 1

6xN2 4

3xN1 + 1

6xN

1

CCCCCA;

and d0; d1; dN1; dN are given by

d0 =

7

18

x0 +

8

9

x1 +

7

18

x2

2

3

d2

d1 =

1

6

x0 +

4

3

x1

1

6

x2

dN1 =

1

6

xN2 +

4

3

xN1

1

6

xN

dN =

7

18

xN2 +

8

9

xN1 +

7

18

xN

2

3

dN2:

If N = 5, this system reduces to

4 1

1 4

!

d2

d3

!

=

6x2 + 1

6x0 4

3x1 + 1

6x2

6x3 + 1

6x3 4

3x4 + 1

6x5

!

:

(5) Implement the Gaussian elimination method with partial pivoting as well as the

method for solving a triangular system by back-substitution.

Use your program to solve several instances of the interpolation problem. Verify that no pivoting is needed.

(6) Implement the LU-factorization method for tridiagonal matrices and test it on the same interpolation problems as in (5).

Do you notice any improvement over Gaussian elimination (running time, numerical precision)?

(7) After computing d1; : : : ; dN1, compute the control points for the Bezier curves C1; : : :CN and write a program to plot these Bezier segments (for t 2 [0; 1]) to visualize the interpolating spline. Experiment with the choice of end conditions.

TOTAL: 300 points.

You'll get 1 file (198.1KB)