Assignment 6 Lecture Notes

The following will act as lecture notes to help you review the material from lecture for the assignment. You can click on the links below to get directed to the appropriate section or subsection.

• Variational Time Integrators
  • Motivating Example: Simple Pendulum
    • Explicit (Forward) Euler
    • Implicit (Backward) Euler
    • Symplectic (Semi-Implicit) Euler
  • Lagrangian Mechanics
    • The Lagrangian and the Euler-Lagrange Equations
    • A Simple Application
  • The Discrete Lagrangian
    • Deriving the Discrete Euler-Lagrange (DEL) Equations
    • Deriving Symplectic Euler for the Simple Pendulum
• Keyframe Interpolation
  • Keyframes
  • Splines
  • Cardinal Splines (Catmull-Rom Splines)
  • B-splines and NURBS

Section 1: Variational Time Integrators

One common application of computer graphics animation is the simulation of physical systems such as the swinging of a simple pendulum or the flow of air over an airfoil. In these simulations, we want to graphically model the behavior of natural phenomena as they evolve dynamically in time. In order to model the behaviors correctly, we not only need to focus on the graphics that we display but also on the numerical computations that determine the state of a physical system at a particular point in time.

A state of a physical system refers to the conditions that define the system at a particular point in time. For instance, consider a ball bouncing vertically in a frictionless environment. At any given moment in time, we can specify the state of the system with the height of the ball relative to the ground and the vertical velocity of the ball. This essentially gives us a snapshot of the configuration of the system at some time $t$.

The evolution of a system from the state at time $t$ to the state at time $t+\Delta t$ is often governed by a system of differential equations that is generally hard or impossible to solve analytically. In many cases, we must resort to using numerical methods to approximate the changes in conditions between states, starting with the initial conditions for the initial state of the system. We refer to these numerical methods specifically as time integration techniques or time integrators.

As an example of a time integrator, recall Euler’s method from introductory calculus: given a position $q\left(t\right)$ and velocity $v\left(t\right)$ at time $t$ and an acceleration function $a\left(t\right)$, we compute the conditions for time $t+\Delta t$ using the equations:

$$q\left(t+\Delta t\right)=q\left(t\right)+\left(\Delta t\right)v\left(t\right)$$

$$v\left(t+\Delta t\right)=v\left(t\right)+\left(\Delta t\right)a\left(t\right)$$

In addition to Euler’s method, there are many other time integrators out there, though not all of them are equal in being “reliable” for graphical simulations. In the following subsections, we discuss one particularly appropriate family of integrators for computer graphics known as variational time integrators, which are derived from Lagrangian mechanics.

Motivating Example: Simple Pendulum

We start off by first examining a simple example of a physical system to demonstrate the different behaviors of different time integrators. Consider a simple pendulum of mass $m$ and length $l$, swinging in a frictionless medium under the influence of the gravitational acceleration $g$. The pendulum only has one degree of freedom in its angle, which we will denote $q\left(t\right)$. Figure 1 shows a visual of the system:

We know from simple physics that the equation of motion for this system is given by:

$$\ddot{q}=\frac{d^{2}q}{d{t}^2}=-\frac{g}{l}sinq$$

where we use the dot notation to denote derivatives with respect to time.

From here, we could theoretically solve for an analytical solution to the system if we were to make a small angle approximation. However, for the sake of the example, let us assume that an analytical solution is not possible (as is the case for many other systems). We would then need to approximate the time evolution of the system using a time integrator.

Let us first try to model the system using the Euler’s method we used in introductory calculus. We will refer to this method by its more technical name: the explicit Euler method. It is also often called the forward Euler method.

Explicit (Forward) Euler

The explicit Euler method is often taught in introductory calculus classes as simply “Euler’s method.” For our convenience, we will reiterate the method here: given a position $q\left(t\right)$ and velocity $v\left(t\right)$ at time $t$ and an acceleration function $a\left(t\right)$, we compute the conditions for time $t+\Delta t$ using the update rules:

$$q\left(t+\Delta t\right)=q\left(t\right)+\left(\Delta t\right)v\left(t\right)$$

$$v\left(t+\Delta t\right)=v\left(t\right)+\left(\Delta t\right)a\left(t\right)$$

To apply the explicit Euler method to the simple pendulum system, we first need to discretize the problem and introduce some notation. We break up the time between our initial time $t_0$ and our final time $t_f$ of simulation into $N$ equal time steps of length $\Delta t$. We denote the time at the $k$-th time step as $t_k=t_0+k\Delta t$ for $k=0,1,...,N-1$. The angle of the pendulum at time $t_k$ is denoted as $q_k$. Let $v=\stackrel{\cdot }{q}=\frac{dq}{dt}$ so that $v_k$ is the angular velocity of the pendulum at time $t_k$.

Using the above notation, we write the explicit Euler update rules as so:

$$q_{k+1}=q_k+\left(\Delta t\right)v_k$$

$$v_{k+1}=v_k-\left(\Delta t\right)\frac{g}{l}sin{q}_k$$

These equations are referred to as explicit equations. Given the values $q_k$ and $v_k$ at time step $k$, we can use the above update rules to explicitly compute the values $q_{k+1}$ and $v_{k+1}$ at the next time step $k+1$. Let us now examine the performance of this integrator. Consider Figure 2:

Figure 2 shows a visual display of the behavior of the explicit Euler integrator. On the left is a display of how the amplitude of the pendulum changes as time goes on. We see that the explicit Euler integrator actually causes the amplitude to increase with passing time. The phase portrait on the right conveys this information in a different manner. Here, we see that the angle $q$, represented by the x-axis, and velocity $v$, represented by the y-axis, both continuously increase from their initial values, which we mark with the dot. The outward “spiral” shown in the phase potrait tells us that the pendulum will move further and faster as time goes on, suggesting a numerical instability where the energy of the system “blows up” and goes towards infinity.

This “blowing up” behavior modeled by the explicit Euler integrator is obviously incorrect. In reality, the energy of the simple pendulum system is conserved, and the amplitude of the pendulum should remain constant. The explicit Euler integrator therefore is not a reliable integrator for simulating physical systems.

What causes this ”blowing up” behavior? To answer this question, we can examine what exactly goes on in the update rules of the explicit Euler integrator. Both the update rules for the integrator are actually first order Taylor approximations. Given $q_k$ and $v_k$, we compute $q_{k+1}$ and $v_{k+1}$ using the tangent slopes of $q$ and $v$ respectively at time $t$. Figure 3 demonstrates this with a visual:

Let the curve in Figure 3 represent the path of motion of the system in phase space. We then see that using a first order Taylor approximation for $q_{k+1}$ and $v_{k+1}$ leads us to a point that is off the actual path of the system. Because this point is off, the subsequent approximation for $q_{k+2}$ and $v_{k+2}$, which depends on the tangent at this “off” point, would then lead to another point that is even further off from the actual path of the system. As we continue making these approximations, the error between the approximated values and the true values will accumulate, eventually causing the approximated values to approach infinity. This results in the “blow up” behavior we see in Figure 2.

Knowing that the explicit Euler method is unreliable, let us try out a different time integrator: the implicit Euler method. This method is also often called the backward Euler method.

Implicit (Backward) Euler

The implicit Euler method is another popular time integrator that is a counterpart to the explicit Euler method. The implicit Euler method takes the update rules of the explicit Euler method and evaluates their right hand sides at the next time step. This causes the update rules to take the form of a set of non-linear equations.

Using the same notation as before, we write the implicit Euler update rules for the simple pendulum system as so:

$$q_{k+1}=q_k+\left(\Delta t\right)v_{k+1}$$

$$v_{k+1}=v_k-\left(\Delta t\right)\frac{g}{l}sin{q}_{k+1}$$

These equations are referred to as implicit equations as opposed to explicit equations and require the use of a non-linear solver to compute $q_{k+1}$ and $v_{k+1}$ given $q_k$ and $v_k$. We examine the behavior of implicit Euler in Figure 4:

Figure 4 shows the behavior of the implicit Euler integrator similar to how Figure 2 shows the behavior of the explicit Euler integrator. We see on the left that the implicit Euler integrator causes the amplitude of the pendulum to decrease with passing time. The phase portrait tells us that $q$ and $v$ both continuously decrease from their initial values, forming an inward “spiral” that suggests numerical dissipation. That is, the energy of the system is artificially dampened by the integrator and goes towards 0.

While the implicit Euler method is not numerically unstable like the explicit Euler method, it also does not correctly conserve the energy of the system. This is not too surprising, since the implicit Euler method is essentially a reverse explicit Euler (hence the name “backward Euler”). We can see this by rearranging the update rules for implicit Euler as so:

$$q_k=q_{k+1}-\left(\Delta t\right)v_{k+1}$$

$$v_k=v_{k+1}+\left(\Delta t\right)\frac{g}{l}sin{q}_{k+1}$$

The above equations show that the implicit Euler is just explicit Euler going backwards in time. Thus, implicit Euler also suffers from error due to first order approximations and is not a reliable integrator for simulating physical systems.

What integrator should we use then if both explicit and implicit integrators are unsatisfactory?

Symplectic (Semi-Implicit) Euler

There is an integrator known as the symplectic Euler method or semi-implicit Euler method, which involves a “strange mix” of the explicit and implicit update rules. The update rules as applied to the simple pendulum system are:

$$v_{k+1}=v_k-\left(\Delta t\right)\frac{g}{l}sin{q}_k$$

$$q_{k+1}=q_k+\left(\Delta t\right)v_{k+1}$$

The idea here is to first explicitly compute $v_{k+1}$ and then implicitly use $v_{k+1}$ to compute $q_{k+1}$. The results of this integrator are shown in Figure 5:

Similar to Figures 2 and 4, Figure 5 shows the behavior of the symplectic Euler integrator. We see from both the pendulum diagram and the phase portrait that the symplectic Euler method does correctly conserve the energy of the system.

The question now is why. Why does the symplectic Euler method produce reliable results? Also, how would we have come up with this “strange mix” of an integrator?

We will find in the following sections that the reason the symplectic Euler method correctly conserves the energy of a system is because it is derived from a principle from Lagrangian mechanics that takes into account the entire path (i.e. all the states) that the system takes to go from its initial state to its final state. In comparison, neither the forward nor implicit Euler methods consider the entire path of the system; instead, they both only consider the time evolution of the system one state at a time.

Lagrangian Mechanics

We present, in this section, a brief overview of Lagrangian mechanics as necessary for the understanding of creating reliable time integrators. Do note though that, because this class is not a physics class, our overview is not very rigorous. For those students who are interested in learning Lagrangian mechanics in depth, we recommend taking an advanced classical mechanics course. Additionally, [2] also provides a rigorous, in-depth analysis of what we cover in the following sections and more (although the work does contain a few typos).

The Lagrangian and the Euler-Lagrange Equations

At the root of Lagrangian mechanics lies the Lagrangian, which is defined to be the difference between the kinetic and potential energies of a physical system:

$$L\left(q,\stackrel{\cdot }{q}\right)=T\left(\stackrel{\cdot }{q}\right)-U\left(q\right)$$

where $L$ denotes the Lagrangian, $T$ denotes kinetic energy, and $U$ denotes potential energy. $q$ and $\stackrel{\cdot }{q}$ denote our position and velocity state vectors respectively.

We define what is called the action, denoted by $S\left(q\right)$, as:

$$S\left(q\right)={\int \; <!--nolimits-->}_{t_0}^{t_f}L\left(q\left(t\right),\stackrel{\cdot }{q}\left(t\right)\right)dt$$

To get an intuitive understanding of what the action is, let us examine for a moment the simple physical system involving a ball thrown from position $q_0$ to position $q_f$. Physics causes the ball to take a certain trajectory or path $\left\{q\left(t\right),\stackrel{\cdot }{q}\left(t\right)\right\}$ as it moves from $q_0$ to $q_f$. We can compute the Lagrangian for each point along this path, and the sum of all the computed Lagrangian values gives us the action.

The reason why we introduce this concept of “action” is because of a key principle in physics known as Hamilton’s least action principle, which states that the action of a physical system is always at a extremum (e.g. minimum). Now, before we discuss what it means for the action to be at an extremum, let us step back and consider what it means for a simple function $f\left(x\right)$ to be extremized at a point $x_t$. For $f\left(x_t\right)$ to be a extremum, we know that $f^{\prime }\left(x_t\right)=0$. From this, we might be tempted to say that for the action, $S\left(q\right)$, to be at a extremum, the derivative of the action with respect to $q$ must be 0. The issue with this statement is that, unlike $x$, $q$ does not just represent a single value. Rather, because the action is defined to be an integral over the values $q_{t_0}$ to $q_{t_f}$, $q$ is actually a path $q\left(t\right)$ between the initial and final states. This makes examining the extremum of the action a bit more complicated than what we are used to.

The key to understanding the impact of the least action principle lies in what we call variational calculus. Of course, an overview of variational calculus is beyond the scope of this course. For our purposes, we simply need to know that the analog of the derivative that we need for examining the extremum of the action is known as the variation, which we denote with $\delta$. Conceptually, a variation of a path $q$, denoted as $\delta q$, is an infinitesimal perturbation (i.e. deviation) to the path at each point except at the endpoints, where the perturbations are null. Similar to how $f^{\prime }\left(x\right)=0$ when $f\left(x\right)$ is a extremum, $\delta S\left(q\right)=0$ when $S\left(q\right)$ is a extremum.

We compute the variation of the action as so:

$$\delta S\left(q\right)=\delta {\int \; <!--nolimits-->}_{t_0}^{t_f}L\left(q,\stackrel{\cdot }{q}\right)dt={\int \; <!--nolimits-->}_{t_0}^{t_f}\left[\frac{\partial L}{\partial q}\cdot \delta q+\frac{\partial L}{\partial \stackrel{\cdot }{q}}\cdot \delta \stackrel{\cdot }{q}\right]dt$$

Using integration by parts yields:

$$\delta S\left(q\right)={\int \; <!--nolimits-->}_{t_0}^{t_f}\left[\frac{\partial L}{\partial q}-\frac{d}{dt}\left(\frac{\partial L}{\partial \stackrel{\cdot }{q}}\right)\right]\cdot \delta qdt+{\left[\frac{\partial L}{\partial \stackrel{\cdot }{q}}\cdot \delta q\right]}_{t_0}^{t_f}$$

Now, because the perturbations at the endpoints of $q$ have to be null, we have that $\delta q\left(t_0\right)=\delta q\left(t_f\right)=0$. This causes the rightmost term of the right-hand side to vanish. Following the least action principle, we then set $\delta S\left(q\right)=0$ and obtain:

$${\int \; <!--nolimits-->}_{t_0}^{t_f}\left[\frac{\partial L}{\partial q}-\frac{d}{dt}\left(\frac{\partial L}{\partial \stackrel{\cdot }{q}}\right)\right]\cdot \delta qdt=0$$

From here, we invoke what is called the fundamental lemma of calculus of variations, which tells us that, because the above integral equals 0 for all $\delta q$, the part of the integrand in the square brackets also equals 0. Applying this lemma gives us the following important formula:

$$\frac{\partial L}{\partial q}-\frac{d}{dt}\left(\frac{\partial L}{\partial \stackrel{\cdot }{q}}\right)=0$$

The above formula is known as the Euler-Lagrange equation and is essentially the core of Lagrangian mechanics. For a given Lagrangian of a physical system, we can use the Euler-Lagrange equations to derive the equations of motion of the system.

A Simple Application

As a demonstration of the power of the Euler-Lagrange equations, consider a physical system involving a particle of mass $m$ moving under the influence of a potential function $U\left(q\right)$, where we use $q\left(t\right)$ to denote the position of the particle. The Lagrangian of the system is simply:

$$L\left(q,\stackrel{\cdot }{q}\right)=\frac{1}{2}m\stackrel{\cdot }{q}{}^2-U\left(q\right)$$

Plugging $L$ into the Euler-Lagrange equations yields:

$$-\frac{\partial U}{\partial q}-m\ddot{q}=0$$

which we can express as:

$$-\nabla U\left(q\right)=F=m\ddot{q}$$

And we see that we obtain Netwon’s law of $F=ma$, which is what we would expect as the equation of motion for this system.

The Discrete Lagrangian

So what about Lagrangian mechanics allows us to create reliable integrators? The answer to this question lies in the fact that the action is a function of the entire path that the given system takes from its initial state to its final state. This fact causes the action to encode the global behavior (i.e. behavior throughout all of time) of the system. And because the Euler-Lagrange equations are derived directly from the action, they inherently also encode the global behavior of the system. As a result, any equations of motion derived from the Euler-Lagrange equations would also inherently take into account the global behavior.

Now, let us say that we manage to construct a complete, discrete equivalent of the action that – like its continuous counterpart – encodes the global behavior of the given system. The update rules that we derive from this discrete action would then also inherently take into account the global behavior. The following subsection details how we construct this discrete action and subsequently derive from it the discrete Euler-Lagrange equations.

Deriving the Discrete Euler-Lagrange (DEL) Equations

Let us first consider how we compute integrals in the discrete setting in general. Given an integral of a function $f\left(t\right)$ from $t_0$ to $t_f$, we can approximate the integral by summing up rectangle approximations of the area under the curve formed by $f$ at discrete time steps of length $\Delta t$. Let $N$ be the number of time steps. The area of the rectangle at time step $k$ for $k=0,1,...,N-1$ is given by $\left(\Delta t\right)f\left(t_k\right)$. The overall integral is then computed via:

$${\int \; <!--nolimits-->}_{t_0}^{t_f}f\left(t\right)\approx \sum _{k=0}^{N-1}\left(\Delta t\right)f\left(t_k\right)$$

We use the same thought process as above to construct a discrete equivalent of the action. We first define the discrete Lagrangian, denoted by $L_d$, as the area of a discrete rectangle given by:

$$L_d\left(q_k,q_{k+1},\Delta t\right)=\left(\Delta t\right)L\left(q_k,\frac{q_{k+1}-q_k}{\Delta t}\right)$$

where $L$ denotes the continuous Lagrangian of the system of interest. Note that we use above a discrete approximation of $\stackrel{\cdot }{q}$, which, for a discrete time step, can be written as $\stackrel{\cdot }{q}\approx \frac{q_{k+1}-q_k}{\Delta t}$. Hence, the discrete Lagrangian ends up being a function of two consecutive states, $q_k$ and $q_{k+1}$, and the size of the time step, $\Delta t$.

We then write the discrete action as:

$$S_d=\sum _{k=0}^{N-1}{L}_d\left(q_k,q_{k+1},\Delta t\right)$$

Following the continuous derivation of the Euler-Lagrange equations, we express the variation of the discrete action as:

$$\delta S_d=\sum _{k=0}^{N-1}\left[D_1{L}_d\left(q_k,q_{k+1},\Delta t\right)\cdot \delta q_k+D_2{L}_d\left(q_k,q_{k+1},\Delta t\right)\cdot \delta q_{k+1}\right]$$

where we use $D_1{L}_d$ to denote the derivative of $L_d$ with respect to the first argument, $q_k$, and $D_2{L}_d$ to denote the derivative of $L_d$ with respect to the second argument, $q_{k+1}$.

We can rearrange the sum of $S_d$ to obtain:

$$\delta S_d=\sum _{k=1}^{N-1}\left[D_2{L}_d\left(q_{k-1},q_k,\Delta t\right)+D_1{L}_d\left(q_k,q_{k+1},\Delta t\right)\right]\cdot \delta q_k$$

$$+D_1{L}_d\left(q_0,q,\Delta t\right)\cdot \delta q_0+D_2{L}_d\left(q_{N-1},q_N,\Delta t\right)\cdot \delta q_N$$

Remember that we require the variations at the boundary points, $q_0$ and $q_N$, to be 0. This eliminates the two terms outside the summation. From here, we follow the continuous derivation of the Euler-Lagrange equations and set the variation of the discrete action to 0:

$$\delta S_d=\sum _{k=1}^{N-1}\left[D_2{L}_d\left(q_{k-1},q_k,\Delta t\right)+D_1{L}_d\left(q_k,q_{k+1},\Delta t\right)\right]\cdot \delta q_k=0$$

Because the sum is 0 for all $\delta q_k$, the expression within the sum must be 0. Hence, we obtain:

We refer to the above as the discrete Euler-Lagrange (DEL) equations. Given two consecutive states (i.e. $q_{k-1}$ and $q_k$), we can use the discrete Euler-Lagrange equations to solve for the next state (i.e. $q_{k+1}$).

We can, however, go further and derive a more useful form of the DEL equations. To do so, let us first consider the DEL equations for the simple Lagrangian:

$$L\left(q,\stackrel{\cdot }{q}\right)=T\left(\stackrel{\cdot }{q}\right)-V\left(q\right)=\frac{1}{2}m\stackrel{\cdot }{q}{}^2-U\left(q\right)$$

The discrete Lagrangian for this system is:

$$L_d\left(q_k,q_{k+1},\Delta t\right)=\Delta t\left[\frac{1}{2}m{\left(\frac{q_{k+1}-q_k}{\Delta t}\right)}^2-U\left(q_k\right)\right]$$

Observe what happens when we compute:

$$D_2{L}_d\left(q_{k-1},q_k,\Delta t\right)=m\left(\frac{q_k-q_{k-1}}{\Delta t}\right)$$

Note how the right-hand side of the equation is just the momentum at time step $k$, which we will denote $p_k$. It turns out that, in general, $D_2{L}_d\left(q_{k-1},q_k,\Delta t\right)$ is indeed just $p_k$ for any Lagrangian. The full proof of this is outside the scope of this class, but we can at least see from the above example how this can be true.

Now, from the above fact and the DEL equations, we have:

$$p_k=D_2{L}_d\left(q_{k-1},q_k,\Delta t\right)=-D_1{L}_d\left(q_k,q_{k+1},\Delta t\right)$$

The above means that we can write the DEL equations in the following form:

Consider for a moment the power of the above two equations. Given the position and momentum of a physical system at some state $k$, we can solve the above two equations for the position and momentum of the next state $k+1$. Unlike the previous form of the DEL equations, which required two consecutive states to compute the next state, this form only needs one state (i.e. the state at the current time step) to compute the next state.

Equations (2) and (3) are the main DEL equations and are what we use to create reliable time integrators. Because they were derived from the action, which takes into account the entire path of the given system, the DEL equations and any time integrators derived from them inherently account for the global behavior of the system. We refer to the time integrators that we derive from the DEL equations as variational time integrators.

Deriving Symplectic Euler for the Simple Pendulum

To demonstate the power of the DEL equations, we will derive the sympletic Euler integrator for the simple pendulum system that we analyzed earlier.

We start by writing the continuous Lagrangian for the system:

$$L=T-U=\frac{1}{2}m\left(\stackrel{\cdot }{x}{}^2+\stackrel{\cdot }{y}{}^2\right)-\left(-mgy\right)$$

Using the relationships:

$$x=lsin\theta$$

$$y=lcos\theta$$

$$\stackrel{\cdot }{x}=lcos\left(\theta \right)\stackrel{\cdot }{\theta }$$

$$\stackrel{\cdot }{y}=lsin\left(\theta \right)\stackrel{\cdot }{\theta }$$

we can rewrite the Lagrangian in terms of $\theta$ and $\stackrel{\cdot }{\theta }$ as:

$$L=\frac{1}{2}m{l}^2\stackrel{\cdot }{\theta }{}^2+mglcos\theta$$

The discrete Lagrangian is then:

$$L_d\left({\theta }_k,{\theta }_{k+1},\Delta t\right)=\Delta t\left[\frac{1}{2}m{l}^2{\left(\frac{{\theta }_{k+1}-{\theta }_k}{\Delta t}\right)}^2+mglcos{\theta }_k\right]$$

Computing the respective derivatives gives us:

$$p_k=-D_1{L}_d=-\Delta t\left[-m{l}^2\left(\frac{{\theta }_{k+1}-{\theta }_k}{{\left(\Delta t\right)}^2}\right)-mglsin{\theta }_k\right]$$

$$p_{k+1}=D_2{L}_d=\Delta t\left[m{l}^2\left(\frac{{\theta }_{k+1}-{\theta }_k}{{\left(\Delta t\right)}^2}\right)\right]$$

Now, recall from classical mechanics that the angular momentum of a simple pendulum is $p=m{l}^2\omega$, where $\omega$ is the angular velocity. Substituting this definition of $p$ into the above equations and doing some algebra simplifies the above two equations into the following equations:

$${\omega }_k=\left(\frac{{\theta }_{k+1}-{\theta }_k}{\Delta t}\right)+\left(\Delta t\right)\frac{g}{l}sin{\theta }_k$$

$${\omega }_{k+1}=\frac{{\theta }_{k+1}-{\theta }_k}{\Delta t}$$

A little rearranging and substitution finally gives us:

$${\omega }_{k+1}={\omega }_k-\left(\Delta t\right)\frac{g}{l}sin{\theta }_k$$

$${\theta }_{k+1}={\theta }_k+\left(\Delta t\right){\omega }_{k+1}$$

which are exactly the equations that we introduced previously as the symplectic Euler integrator!

Section 2: Keyframe Interpolation

When we animate scenes in computer graphics, we often use a technique called keyframe interpolation to generate the animation. The idea of keyframe interpolation is best introduced with an example, so let us consider the following scenario.

We want to create a minute-long, traditional animation where each frame is hand-drawn. For the animation to appear smooth, we need to display around 30 frames per second. However, 30 frames per second for 60 seconds results in a grand total of 1800 frames to be drawn! And that is only for a minute-long animation. The amount of time and manpower it would take to generate all the frames necessary for an even longer animation would be enormous. However, with the advent of computer animation, we can now generate long animations by drawing only a small number of select frames called keyframes. We draw the keyframes, which capture the most important scenes in our animation, and generate the rest of the animation by interpolating or, in other words, using computer algorithms to “fill in the gaps” between the keyframes.

Keyframes

One particular way of specifying a keyframe is with the positions and orientations of all the points in the scene of interest. We specify the position of a point as a translation from an arbitrary chosen origin and the orientation of a point as a rotation by an angle about some axis. In addition, we also find it convenient to specify the scaling of a point as another characterization of the point’s configuration in a scene. When we interpolate between keyframes, we actually interpolate the change in translation, rotation, and scaling of each point in the scene between the frames.

Splines

We most often interpolate across keyframes by piecing together individual polynomial curves to form an overall, piecewise polynomial function called a spline. The spline goes through and connects all corresponding points across all keyframes. To see how this process works, let us consider piecing together a spline based on linear interpolation between keyframes.

Let us say that we currently have two keyframes, $K_1$ and $K_2$, that we want to linearly interpolate between. In $K_1$, we have a point that is positioned in the scene by the translation vector $p_0$. In $K_2$, that same point is now positioned in the scene by the translation vector $p_1$. The scenes in both keyframes use the same coordinate system. To linearly interpolate the change in position of the point between $K_1$ and $K_2$, we need to linearly interpolate between the translations vectors $p_0$ and $p_1$.

Our linear interpolation scheme proceeds as follows. First, we should note that we need to interpolate the components of the translation vectors one at a time. However, for simplicity in the notation, we will still use the variables, $p_0$ and $p_1$, to denote the values that we are interpolating. Now, consider a line segment connecting $p_0$ to $p_1$. We parameterize this line segment over the unit domain with the following function:

$$f\left(u\right)=\left(1-u\right)p_0+u{p}_1$$

where $u\in \left[0,1\right]$. We can then use $f\left(u\right)$ to compute the values between $p_0$ and $p_1$ to “fill in” the motion that our point takes as it moves between the two keyframes.

Now, for reasons that we will see later, we find it convenient to consider the following alternative formulation of $f\left(u\right)$:

$$f\left(u\right)=a_0+u{a}_1$$

for coefficients $a_0$ and $a_1$. We notice that the above formulation is the canonical representation of a first degree polynomial (i.e. a line). We can rewrite this formulation using vector notation. Let $u=\left[1,u\right]$ and $a=\left[a_0,a_1\right]$. Then, we have:

$$f\left(u\right)=u\cdot a$$

What exactly is $a$ though? Well, we know that $f\left(0\right)=p_0$ and $f\left(1\right)=p_1$. We also know that for $p_0$, we have $u=\left[1,0\right]$. And for $p_1$, we have $u=\left[1,1\right]$. With this, we can write the following system of equations:

$$\left[\begin{array}{c}\hfill p_0\hfill \\ \hfill p_1\hfill \end{array}\right]=\left[\begin{array}{cc}\hfill 1\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 1\hfill \end{array}\right]\left[\begin{array}{c}\hfill a_0\hfill \\ \hfill a_1\hfill \end{array}\right]$$

Let us write the above system as:

$$p=Ca$$

where we call $C$ the constraint matrix. Note how we can solve for $a$ after finding the inverse of $C$. We denote the inverse of $C$ with the variable $B$ and call it the basis matrix. We can now evaluate the curve between $p_0$ and $p_1$ as:

$$f\left(u\right)=uBp$$

We repeat the above process for each pair of values $p_i$ and $p_{i+1}$ in consecutive keyframes $K_i$ and $K_{i+1}$. In the end, we obtain multiple line segments that connect sequentially to form a spline that goes through all points $p_0,p_1,p_2,...,p_{final}$. We can use the pieces of the spline to compute the components of all translation vectors needed between keyframes to generate the animation across all keyframes.

Now, at this point, one might ask, “why did we have to make this process so complicated?” After all, our original, parameteric formulation would have been sufficient for the interpolation between two points. The answer is that while this process is unnecessary for simple linear interpolation, the whole concept of computing constraint and basis matrices becomes useful when we interpolate using more complicated curves.

Furthermore, we can imagine that simply using linear interpolation between keyframes would not produce smooth motion. Consider the shape of a spline formed from only line segments. Line segments do not transition smoothly into one another; rather, they form sharp angles at each junction. These sharp angles cause the motion of the point to result in sudden changes as we transition between keyframes. As a result, we do not obtain smooth motion from linear splines and must instead look towards more complicated curves for interpolation.

Cardinal Splines (Catmull-Rom Splines)

Before we continue, let us stop to consider the nature of a spline that would produce smooth interpolation across keyframes. We know that a spline consisting of line segments does not work because the segments form sharp angles at each junction. What is a condition that we can enforce to prevent our spline from having sharp angles?

The answer is to require the first derivatives of two connecting pieces of our spline to be equal at the junction (i.e. $C^1$ continuity). This means that we must, at the very least, use quadratics for interpolation.

It turns out that cubics are the most commonly used polynomials for keyframe interpolation. In particular, we often use a specific class of cubics called cardinal cubic curves. We call the overall resulting curve a cardinal spline. They have the nice properties of $C^1$ and $C^2$ continuity (i.e. their first and second derivatives are continuous everywhere).

We form a cardinal cubic curve in the following manner. First, let us denote the two keyframes that we are interpolating between as $K_i$ and $K_{i+1}$. Let $p_i$ and $p_{i+1}$ be the values (e.g. components of translation or scaling vectors) that we are interpolating between the two keyframes. To form a cardinal cubic curve, we need to also consider the keyframes, $K_{i-1}$ and $K_{i+2}$ and their corresponding values, $p_{i-1}$ and $p_{i+2}$. Thus, the cardinal cubic curve has four control points indexed from $i-1$ to $i+2$.

We let the derivative at the beginning of our curve segment be a scaling of the vector between the first and third control points. And we let the derivative at the end of our curve segment be a scaling of the vector between the second and fourth control points. Figure 1 below provides a visual of this.

Let us take a moment to convince ourselves that the derivatives of two consecutive cardinal curves match at the point in which the curves connect. Consider a cardinal curve $c_1$ between points $p_i$and $p_{i+1}$ connecting to a cardinal curve $c_2$ between points $p_{i+1}$ and $p_{i+2}$. The two curves connect at point $p_{i+1}$. Now, the derivative of $c_1$ at $p_{i+1}$ is given by $p_{i+2}$ - $p_i$, the difference between the second and fourth control points of $c_1$. What is the derivative of $c_2$ at $p_{i+1}$? It would the the difference between the first and third control points of $c_2$. The first and third control points would be $p_i$ and $p_{i+2}$ respectively, hence, the derivative of $c_2$ at $p_{i+1}$ would also be $p_{i+2}$ - $p_i$. Thus, the derivatives match at point $p_{i+1}$. Similar logic can be used to confirm that the derivatives at other junctions along a cardinal spline also match.

To find a mathematical representation for the cardinal curve, we undergo a process similar to the one we took for our linear interpolation example above. We first represent the cardinal curve between $p_i$ and $p_{i+1}$ with a function, $f\left(u\right)$, that is parameterized over the unit domain. Our constraints for the curve then can be expressed as:

$$f\left(0\right)=p_i$$

$$f\left(1\right)=p_{i+1}$$

$$f^{\prime }\left(0\right)=s\left(p_{i+1}-p_{i-1}\right)$$

$$f^{\prime }\left(1\right)=s\left(p_{i+2}-p_i\right)$$

where $s$ is the scaling factor for the derivative vectors. We formally define $s$ to be:

$$s=\frac{1}{2}\left(1-t\right)$$

where $t\in \left[0,1\right)$ is known as the tension parameter and is an arbitrarily set value that controls the bending of the curve. The closer we set $t$ to 1, the flatter the curve becomes around the control points. We often, however, simply set $t$ to 0. We call the set of splines for $t=0$ Catmull-Rom splines.

Since cardinal curves are cubics, we will consider the canonical representation of a third degree polynomial:

$$f\left(u\right)=a_0+u{a}_1+u^2{a}_2+u^3{a}_3=u\cdot a$$

for $u=\left[1,u,u^2,u^3\right]$ and $a=\left[a_0,a_1,a_2,a_3\right]$. We want to now write a system of equations in the form of $p=Ca$. To do so, we consider our constraints and solve for the control points:

$$p_{i-1}=f\left(1\right)-\frac{1}{s}{f}^{\prime }\left(0\right)$$

$$p_i=f\left(0\right)$$

$$p_{i+1}=f\left(1\right)$$

$$p_{i+2}=f\left(0\right)+\frac{1}{s}{f}^{\prime }\left(1\right)$$

Substituting our canonical cubic polynomial into the above set of equations yields:

$$p_{i-1}=a_0+\left(1-\frac{1}{s}\right)a_1+a_2+a_3$$

$$p_i=a_0$$

$$p_{i+1}=a_0+a_1+a_2+a_3$$

$$p_{i+2}=a_0+\frac{1}{s}{a}_1+\frac{2}{s}{a}_2+\frac{3}{s}{a}_3$$

which we can express in the form of $p=Ca$. Taking the inverse of the resulting constraint matrix yields the following basis matrix:

$$B=C^{-1}=\left[\begin{array}{cccc}\hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill -s\hfill & \hfill 0\hfill & \hfill s\hfill & \hfill 0\hfill \\ \hfill 2s\hfill & \hfill s-3\hfill & \hfill 3-2s\hfill & \hfill -s\hfill \\ \hfill -s\hfill & \hfill 2-s\hfill & \hfill s-2\hfill & \hfill s\hfill \end{array}\right]$$

For the case where $t=0$ (i.e. for Catmull-Rom splines), we have:

$$B=C^{-1}=\frac{1}{2}\left[\begin{array}{cccc}\hfill 0\hfill & \hfill 2\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill -1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 2\hfill & \hfill -5\hfill & \hfill 4\hfill & \hfill -1\hfill \\ \hfill -1\hfill & \hfill 3\hfill & \hfill -3\hfill & \hfill 1\hfill \end{array}\right]$$

Hence, our cardinal curve function can be expressed as:

$$f\left(u\right)=uBp$$

Given four control points expressed as $p=\left[p_{i-1},p_i,p_{i+1},p_{i+2}\right]$ and an appropriate $t$ parameter, we can use the above to compute the interpolated values between $p_i$ and $p_{i+1}$. This interpolation method smoothly interpolates the components of both translation and scaling vectors.

To interpolate the rotation quaternions, simply use this interpolation method on each component of the keyframe rotation quaternions, then normalize the result. Refer back to our Assignment 3 lecture notes for details about quaternions.

There are other methods for rotation interpolation, including Spherical Linear Interpolation (SLERP) [3] [4], which interpolates quaternions along the great arc of a sphere. However, SLERP does not have $C^1$ or $C^2$ continuity because it only uses the nearest two control points for interpolation. Thus, it can produce awkward motion in comparison to more sophisticated methods [5].

B-splines and NURBS

Another major class of splines is the B-spline, which is a type of non-interpolating spline. This means that a B-spline is not required to pass through its control points. An extension of B-splines called NURBS (Nonuniform rational B-splines) can have control points that are non-uniformly spaced apart. This is useful for adding control points to a NURBS curve without affecting the shape of the curve, which is extremely useful in modeling for allowing more local control of the curve. NURBS can also be used to define smooth surfaces in 3D using a technique called lofting, which combine two NURBS curves into a single surface. If you are interested in reading about the details behind NURBS and B-splines in general, feel free to check out these notes on NURBS and NURBS surfaces from previous years of CS171. For an interactive website which allows you to see how B-splines and NURBS splines change when you move the control points, click here and here.

References

[1] http://www.geometry.caltech.edu/pubs/SD06.pdf

[2] http://thesis.library.caltech.edu/2492/1/west_thesis.pdf

[3] http://dl.acm.org/citation.cfm?id=325242

[4] https://en.wikipedia.org/wiki/Slerp

[5] http://www.geometry.caltech.edu/pubs/HTZBD10.pdf