|
Designing
Parametric Interpolation with C++ Templates
by Ismail Ziauddin
1. Introduction:
I describe implementation
of core algorithms for parametric interpolation of
Points, Rotations or Scalars, essentially in N
dimensional space and for N = 1 the formulation applies
to scalars. C++ template facility provides for a
function to work on different types of variables. Such
an implementation could be useful for building a
key-frame animation engine.
Parametric curves
interpolation techniques require efficient algorithms
for arclength parameterization and search for parameter
values given arclength. There are many approaches taken
for arclength paramerization, and is a subject of
ongoing research. I describe here a technique very easy
to implement that minimizes numerical error.
In Section 2 I give a
very general overview of parametric curves. A number of
excellent sources are available and the subject is well
understood. The formulation is based on [1]. Section 3
describes Arclength parameterization scheme. Section 4
gives a short description of implementing rotational
interpolation using a similar formula as in case of 3D
points, a full description and derivations are given in
[4].
2. Implicit And
Parametric Functions:
An implicit form of a
function that defines a curve in three dimensional space
is
 (1)
It is convenient in some
applications to define a curve in the form
(2)
Here X, Y and Z are
univariate polynomials in parameter u. As u varies from
0 to 1 in discrete steps, X, Y and Z trace out a curve
in 3D space. The above parametric form can be written
out as
(3)
It would be difficult to
determine constants
to achieve desired shape of curve. The shape of curve
can be intuitively outlined by a set of few control
points. Then the constants
can be replaced by points in 3D space in combination
with linearly independent basis polynomial functions in
u. This transformation can be written as
(4)
Where
are 3D points and
are basis functions.
For a cubic, that is a
curve where equation is of degree 3 in u, there will be
four control points and four linearly independent basis
functions.
Properties of parametric
curves:
Convex hull:
==> the curve lies entirely within the polygon drawn
by connecting the control points.
Affine invariance:
An affine transformation applied to the control points
results in the same curve as the curve obtained by
applying the same transformation to the original curve
of un-transformed control points.
Composite curves:
A cubic curve formed by four control points is called a
curve segment. Two curve segments
and 
can be contiguously joined by appropriately selecting
two sets of control points. The continuity will be
 in
case of bezier basis, while
for other basis such as B-Spline basis. (These basis
functions are described below)
continuity
means that derivatives
 ,
and
are the same at the joining point, i.e.,
(5)
Basis functions:
Bezier basis:
A Bezier curve of degree n can be defined in terms of a
set of points
(i=0, 1, ...n) as
(6)
(7)
Where
are binomial coefficients
(8)
For cubic case n=3,
(9)
B-Spline basis:
These are obtained by imposing three constraints of
 ,
 and
continuity at segment joints of two cubic parametric
curves.
(10)
Catmull-Rom basis:
Bezier and B-Spline basis are not interpolating, but a
class of splines known as Catmull-Rom are. These are
defined for a cubic as
(11)
NURBs:
These are 4-dimensional
rational polynomials where 4th dimension is called
weight. The projection onto 3D is obtained by
rationalizing the polynomial, that is, dividing by
weights polynomial. The curve consists of only one
segment parameterized by global parameter u which takes
values from 0 to 1.
Specifically
(12)
Where
= weight of ith control point
is the 4D control point
(13)
Blending functions
are defined recursively. Hence
(14)
Knots:
These given by
determine how and where the pieces of curves are joined.
These are increasing, non-uniformely spaced and not
necessarily distinct values. The order of the curve is
determined by number of knots (counted with multiple
values) minus number of control points, and degree of
curve is order less 1.
Weights:
With each control point there is a value attached called
weight that is considered 4th dimension. The curve
reduces to a non-rational one when all weights are
equal. Otherwise a weight controls how "hard"
the curve is pulled to the control point relative to
other points. The weights allow construction of conic
sections like circles, etc. by NURBs.
Implementation:
Here is a sample implementation of NURB basis function
// return value of basis
function for parameter u,
// t is array of knot values
REAL nb(const int i,
const int k, const REAL u, const REALARRAY& t)
{
if (k==1)
{
return ( (t[i] <= u)
&& (u < t[i+1]) ) ? 1.0 : 0.0;
}
else
{
const bool aa = (t[i+k] ==
t[i+1]);
const bool bb = (t[i+k-1] ==
t[i]);
return
(!aa && !bb) ? 0.0 :
(!bb) ? ((t[i+k]-u) / (t[i+k]-t[i+1]))
* nb(i+1, k-1, u, t) :
(!aa) ? ((u-t[i])/(t[i+k-1]-t[i]))
* nb(i, k-1, u, t) :
((u-t[i])/(t[i+k-1]-t[i]))
* nb(i, k-1, u, t) + ((t[i+k]-u) /
(t[i+k]-t[i+1])) *nb(i+1, k-1, u, t);
}
}
3. Arclength
Parameterization:
In many applications it
is useful to parameterize the curve by length of the
curve rather than by parameter u. For example in
animation, it is required to place an object at a
distance on the curve at some time t from starting
point. The distance might have been obtained from speed
of the object.
The arclength
parameterization differs from parameter u can be seen
from that points on the curve at equal parameters u are
not the same as points at equal distance parameters.
The goal is to obtain the
point on the curve given a length of curve segment.
Since a point on curve is function of u, it implies
where
is length. The solution requires computation of length,
and formulation of function  .
Computation of length:
If we consider an infinitesimal segment of the curve
with length
 ,
then it’s value can be computed from
(15)
(16)
The length of the curve
between parameter values
and  can
be obtained from
(17)
The integral can be
transposed as function of
with control points constant coefficients, and a numeric
integration technique like Romberg Integration or
Simpson’s Rule can be applied. It is also possible to
divide curve segments successively into increasing
number of subdivisions, computing length as sum of
lengths of subdivisions, and repeating the process until
difference between successive computation is small
enough. This is one way of integrating the function of
,
and I describe a method based on equation (15) in
implementation section
Search for u:
We want to find
given arclength
.
Since no closed-form solution exists for function
 ,
numeric technique needs to be applied. It is observed
that both
and
are monotonically increasing, hence a binary-search for
should converge.
Consider a segment of
cubic curve with four control points. Let it’s length
be .
Let us say we are looking for parameter value of
given length
 .
Divide the interval
([0,1] in beginning) into two intervals
and  .
Compute length of the first interval. If required length
is greater than computed length then required parameter
lies in the second interval, otherwise it lies in the
first interval. This process is repeated until length of
interval remains no greater than a small tolerance
value, or difference between lengths obtained from two
intervals is no greater than a small tolerance value.
Then midpoint of parameter interval will be required
corresponding to given length
 .
In practice a curve may
be composed of many curve segments, hence the segment of
the curve needs to be determined and then the procedure
can be applied to that segment, and an array of
arclengths of each segments is required to be stored.
NURB curves have only one segment.
4. Rotational
interpolation:
One widely used technique
for rotational interpolation is to represent rotations
by Quaternions, also known as Euler Parameters, in
contrast to representing rotations by Euler Angles that
has problems. Ken Shoemake brought Euler Parameters to
the attention of Graphics community, for example
reference [2, 3]. This scheme is well known as Spherical
Linear Interpolation, and there is less well known
Spherical Quadratic Interpolation (Squad) that is
extensive but interpolates better.
I describe a method that
is based on formulation given in [4]. Essentially, a
parametric form given by equation (4) can be rewritten
as what is called a cumulative form:
(18)
Where
are basis functions as before in case of noncumulative
form. Replacing Points in 3D space by quaternions, the
above form becomes (see [4] for complete derivation)
(19)
And (19) reduces to
(20)
For a Cubic, equation
(20) can be expanded that is suitable for direct coding
as
CLICK
HERE FOR CODE
Noting from (4) and (9)
that
Reference [4] has
descriptions for other basis functions like B-Splines,
NURBs and Hermite interpolations.
This scheme is better for
implementation in our case, because there is only an
extra function required to do rotational interpolation
in the function object that does interpolation. This
function object could be for Bezier, B-Spline or NURB
basis.
4. Implementation:
Computation of length:
//A sample
inplementation of arclength computation for one segment
of a curve.
// Template parameter T is the type, a point, a scalar
or a Euler Parameter Rotation,
// each has a function of ArcLength () defined.
// Template parameter F is structure that has a function
Evaluate(u) that returns a value of type T given u
template<typename T, typename F>
double ArcLength(const double a, const double b, const
int res, F& func)
{
double sum = 0.0;
const double du = (b - a) /
double(res-1);
double u = a;
T v1 = func.Evaluate(u);
for (int j=0; j<res-1; j++)
{
u += du;
T v2 = Evaluate(u);
sum += ::ArcLength(v1, v2);
v1 = v2;
}
return sum;
}
// actual computation of
arclength between two parameter values a and b. Calls
above ArcLength routine
// with res=number of sub-divisions required. Computes
until difference between two consecutive values is no
// greater than a small tolerance value
template<typename T, typename F>
double ArcLength(const double a, const double b, F&
func)
{
int res = KRESOLUTION;
// number of sub-divisions
double s1 = 0.0;
double s2 = ArcLength(a, b, res, func);
while (Compare(s1, s2, tol))
// return 0, -1, 1 for s1==s2, s1<s2, s1>s2 within
tol
{
s1 = s2;
res += KRESINC;
// increase resolution by a fixed number
s2 = ArcLength(a, b, res, func);
}
return s2;
}
Search for u:
// return parameter value u given arc-length l, by
binary search described in text
template<typename T, typename F>
double SearchuInSeg(const double l, F& func)
{
double la = 0.0;
double u;
double u1 = 0.0;
double u2 = 1.0;
while (true)
{
u = (u1+u2)/2.0;
double tmp = ArcLength(u1, u, func);
double lb = la+tmp;
int
cmpr = Compare(l, lb, tol);
if (cmpr == 0)
break;
else if (cmpr == -1)
u2 = u;
else
{
la += tmp;
u1 = u;
}
}
return u;
}
References:
[1] Alan Watt, Mark Watt:
Advanced Animation and Rendering Techniques,
Addison Wesley Publishing, 1992
[2] Shoemake K.:
Animating Rotations With Quaternion Curves, Computer
Graphics, SIGGRAPH Proceedings 1985.
[3] Shoemake K.:
Quaternion Calculus for Animation, SIGGRAPH Course notes
#2, 1991.
[4] M. J. Kim, M. S. Kim,
S. Y. Shin: A General Construction Scheme for Unit
Quaternion Curves With Simple High Order Derivatives,
Computer Graphics, SIGGRAPH Proceedings 1995
Contact Ismail
directly at ismailz@rediffmail.com
|
