## slerp between three or more points

### slerp between three or more points

Dear all,

My field of research is computational mechanics. In particular, I have
been studying geometrically exact beam theory, in which the
configuration variables are, for each each cross, a vector describing
its position and an orthogonal tensor describing its orientation.

Now, in the finite element method you must define some kind of
interpolation. In finite elements for continua, one usually adopts
Lagrange polynomies, since the nodal variables are all translations
belonging to a vector space. However, in structural finite elements,
nodal variables include also rotations, belonging to the group SO(3)
which is not a linear group. Hence, interpolation of rotations becomes

So far, there have been authors who interpolate spins, which are
subsequentely used to update the orthogonal tensors via the
exponential mapping. Others, including myself, have resorted to
vector-like parametrization. Others still, have performed
interpolations which do not enforce the orthogonality constraint. All
those techniques have some problems (path-dependency, lack of
invariance under superimposed rigid body motions, etc).

Recently, there has been some work which resorts to SLERP, a technique
that I understand has been very popular in the computational graphics
comunity. SLERP is as far as I can see the perfect solution for
interpolating between two rotations. However, in the finite element
method, one is often interested in high-order elements, having more
than two nodes.

So, my question is: Are you aware of SLERP-like intepolations for
rotations between more than two rotations?

I have found some resources over the internet on this subject, but in
contrast to SLERP, it seem there is no simple and elegant solution. I
have lurked this newsgroup for some months, and it seems to be a very
warm forum frequented by very proficient people. So perhaps someone
can tell me what are the alternatives.

Thanks for your attention

Manuel Ritto-Corra

PS: I've posted this message yesterday, but I cannot find it in
google. So, assuming there is some problem in my news server, I am
reposting it through google. Sorry, if you see it in duplicate.

### slerp between three or more points

Hi Manuel

My understanding of your question is that you want to define a
spherical-linear interpolation between three (or more) orientations.  In
other words, the SO(3) equivalent of x = uA + vB + wC.  I will suppose that
an orientation is represented by R, being the rotation applied from some
canonical orientation.

The short answer is, such a beast will not behave as you would like it to.

The chief problem you face is not so much the non-linearity of SO(3) as the
non-commutative nature of rotations.  In fact, the interpolation you are
after is ill-defined (at least in a mathematical sense).  To demonstrate
this, it helps to form an algebra for rotations.  I'll introduce three
notations:

A + B        indicates rotation A followed by rotation B
uA            is equivalent to the matrix exponentiation A^u; u scales the
rotation A.
B - A        is a rotation such that A + (B - A) = B

Now, in the case of a slerp, we have R = uA + v(B-A), where u+v = 1.  It is
very tempting to rewrite this as a "conventional" interpolation, R' = u'A +
vB.  However, A and B do _not_ commute (unless A=B) and therefore R' /= R.
What, then, are we to make of the linear combination R = uA + vB + wC?
Clearly such a combination has little to do with interpolation between three
orientations, at least in the sense we might expect it to.

If you're prepared to accept a system where u, v and w do not map linearly
to angular distance (which I suspect you're not), you could try the
following:  Slerp between A and B to find D, and then slerp between D and C.
The expression for this is R = v(uA + (1-u)(B-A)) + (1-v)C = uvA + (1-u)B +
(1-u)A + (1-v)C.  In matrix form, assuming operators-on-left:

R = C^(1-v) A^(1-u) B^(1-u) A^(uv)

v=1 puts you at C; and if v=0, then u=0 puts you at A and u=1 puts you at B.
v has a linear mapping to angular distance which is parameterized on u, i.e.
phi(v, u) = v alpha(u).

You may have a bit more luck coercing alpha(u) to be a constant function by
changing variables.

I'd be interested to hear other people's opinions on this.

One theoretical path that I took about a year ago was to ask the question:
Given rotations A and B, what does one make of the expression lim as t->0
(A^t B^t)^(1/t)?  Clearly the limit must be the same if A and B are
reversed.  Thus under this operation I have a commutative algebra of
rotations, which I can form into linear combinations.  Unfortunately I
couldn't make anything of the expression :o).

Regards
David Turner

### slerp between three or more points

> I'll introduce three

Quote:> notations:

> A + B        indicates rotation A followed by rotation B

...

It is enough for me. Please, don't introduce such abominations. An operator
product is not a sum. A sum *is* commutative.

I don't know what Manuel Correa precisa verdadeiramente, he must explain what
is a 'higher order' interpolation:

1. A polynomial interpolation through many intemediate orientations, or
2. A multi-linear interpolation with several contributing orientations.

Then we can continue to speculate on a safer ground.

Jerzy Karczmarczuk

### slerp between three or more points

>  > I'll introduce three
> > notations:

> > A + B        indicates rotation A followed by rotation B

> ...

> It is enough for me. Please, don't introduce such abominations. An
operator
> product is not a sum. A sum *is* commutative.

Like it or not, the notation does clarify the situation.  If you'd read
further down you'd have noticed that I switched to a more conventional
notation after explaining why it doesn't work.  And '+' is _not_ necessarily
a commutative operator.

Quote:> 1. A polynomial interpolation through many intemediate orientations, or
> 2. A multi-linear interpolation with several contributing orientations.

I think from context it's clear that he means (2).  I may be wrong, in which
case he should have read the FAQ first; but I'm prepared to give him the
benefit of the doubt.

Regards
David Turner

### slerp between three or more points

Thanks David Turner and Jerzy Karczmarczuk for your feed-back.
It seems I have not expressed my self in a completely clear way.
David's comments seem to assume that I have a two-dimensional domain
to interpolate (instead of one-dimensional). Also, I am note sure I
understand the question correctly, but my answer to Jerzy Karczmarczuk
question seems to be:
3.A *polynomial* interpolation with *several* contributing
orientations.
Please read below, where I detail my request for help.

I have a one-dimensional domain (the length of a beam) which is
parametrized by one scalar argument s. I want to consider an
orthogonal tensor function that maps each value s to a rotation tensor
R(s) (rotation matrix if you want).

For the sake of simplicity, let us assume that my scalar argument s
can vary between 0 and 1 and that I know the value of an orthogonal
tensor R for 3 nodes, located at s=0, s=1/2 and s=1. (In fact, I want
to consider the general case in which I know the rotation at n
arbitrarily located values of s.)

So, given R(0), R(1/2) and R(1), how do I build a *smooth*
interpolation scheme for R(s)?

By smooth, I mean that I want the curvature to be continuous for
0<s<1. (The curvature K is the axial vector of the skew symmetric
tensor given by  Rt R' , where Rt is the transpose of R and R'=dR/ds.)
In other words, I am looking for the equivalent of Lagrange polynomies
for SO(3). Maybe, I am asking too much...

Obviously SLERP gives the perfect solution for n=2, since it leads to
a constant curvature. But, as far as I am aware, no simple
generalization for n>=3 seems to exist. Let me point out that the use
of two SLERPs (between s=0 and s=1/2 and between s=1/2 and s=1) is not
acceptable, since it is not smooth at all at s=1/2.

So, I hope all is now clear.
Thanks again,

Manuel Correa

(I am replying to myself, because although I have read the answers
from my news server, it seems I have to post through Google to be read
by others. And in google, I still cannot see the answers to my
question)

### slerp between three or more points

Quote:>So far, there have been authors who interpolate spins, which are
>subsequentely used to update the orthogonal tensors via the
>exponential mapping. Others, including myself, have resorted to
>vector-like parametrization. Others still, have performed
>interpolations which do not enforce the orthogonality constraint. All
>those techniques have some problems (path-dependency, lack of
>invariance under superimposed rigid body motions, etc).

>Recently, there has been some work which resorts to SLERP, a technique
>that I understand has been very popular in the computational graphics
>comunity. SLERP is as far as I can see the perfect solution for
>interpolating between two rotations. However, in the finite element
>method, one is often interested in high-order elements, having more
>than two nodes.

>So, my question is: Are you aware of SLERP-like intepolations for
>rotations between more than two rotations?

From the context of your question, my guess is you want something
other than a quaternion spline, probably more of a "spherical average"
idea. The non-commutative algebra of 3D rotations manifests itself, as
you suggest, in a non-flat space, with the consequence that many flat
space assumptions fall apart. Although it has been nearly two decades
since Shoemake's first paper, surprisingly little progress has been
made on questions like yours -- which suggests it ain't easy.

One recent paper that may be of interest appeared in TOG.

<http://math.ucsd.edu/~sbuss/ResearchWeb/spheremean/>

### slerp between three or more points

Hello

Quote:> So, given R(0), R(1/2) and R(1), how do I build a *smooth*
> interpolation scheme for R(s)?

Wonder of wonders.  Ok, Jerzy, I was thoroughly wrong ;-)

Manuel: what you want is a spline.  I suspect that the natural cubic spline
will suit you best.  Search for "Shoemake squad" (there are some references
in the c.g.a FAQ).  I have also seen Catmull-Rom splines implemented using
quaternions.  It depends on what you want.  Do some research.

Regards
David Turner

### slerp between three or more points

> Manuel: what you want is a spline.  I suspect that the natural cubic spline
> will suit you best.  Search for "Shoemake squad" (there are some references
> in the c.g.a FAQ).  I have also seen Catmull-Rom splines implemented using
> quaternions.

Actually, it seems that the analog of Catmull-Rom is more appropriate, since
the intermediate known points are nodes. Hm. Just a spline ... But preserving
the fact that the interpolated quaternion be unitary, no?

See here:

http://www.cis.uab.edu/info/faculty/jj/papers/mostnatural.pdf

http://citeseer.nj.nec.com/johnstone99rational.html
(actually, perhaps rather directly here:
http://www.cis.uab.edu/info/faculty/jj/papers/ratqspline.pdf
)

http://www.cs.technion.ac.il/Labs/Isl/Project/Projects_done/quaternio...

Jerzy Karczmarczuk

### slerp between three or more points

> One recent paper that may be of interest appeared in TOG.

>   <http://math.ucsd.edu/~sbuss/ResearchWeb/spheremean/>

Wow.  Thanks for the link.  The applications are mind-boggling.

Regards
David Turner

### slerp between three or more points

Quote:>Obviously SLERP gives the perfect solution for n=2, since it leads to
>a constant curvature. But, as far as I am aware, no simple
>generalization for n>=3 seems to exist. Let me point out that the use
>of two SLERPs (between s=0 and s=1/2 and between s=1/2 and s=1) is not
>acceptable, since it is not smooth at all at s=1/2.

The same paper that introduced Slerp for two points also went on to
construct smooth curves through multiple successive points. As have
any number of later papers. Furthermore, a partial answer is given
explicitly in the FAQ.

"Subject 5.25: What's the big deal with quaternions?"
<http://www.faqs.org/>

|         Squad(q1,a1,b2,q2;t) = Slerp(Slerp(q1,q2;t),
|                                      Slerp(a1,b2;t);
|                                      2t(1-t))

The FAQ omits the crucial detail of how to set a1 and b2 to get C1
continuity. At t=0 the curve hits q1 with first derivative

q1 (log(q1^-1 q2) + 2 log(q1^-1 a1))

Reversing to Squad(q2,b2,a1,q1;1-t) gives the other end properties. As
if by magic, everything "just works" when you let ai=bi while setting

-1               -1
log(q    q   ) + log(q    q   )
i    i+1         i    i-1
a  = q  exp(- -------------------------------)
i    i                      4

This in a sequence of interpolation points qk, k=1..n. The logarithm
of a unit quaternion [U sin(t),cos(t)], with U a unit vector, is pure
vector tU, and vice versa for the exponential. (I am quoting from Ken
Shoemake's Siggraph course notes; he never published Squad and its
tangent methods, despite their eventual popularity!)

If you want second derivative continuity then you need a more
sophisticated scheme, but those have been published also. Check the
references and methods of the paper I cited in my previous post.

### slerp between three or more points

Quote:> From the context of your question, my guess is you want something
> other than a quaternion spline, probably more of a "spherical average"
> idea. The non-commutative algebra of 3D rotations manifests itself, as
> you suggest, in a non-flat space, with the consequence that many flat
> space assumptions fall apart. Although it has been nearly two decades
> since Shoemake's first paper, surprisingly little progress has been
> made on questions like yours -- which suggests it ain't easy.

> One recent paper that may be of interest appeared in TOG.

>   <http://math.ucsd.edu/~sbuss/ResearchWeb/spheremean/>

I've been busy following the sugestions people kindly made to me. So,
first of all thanks to all of you. My comments follow.

Many of you think that what I really need ("o que eu preciso
verdadeiramente", as Jerzy might say) is a spline. Well I've some
(yes, I've read the FAQ before), but from my point of view, they look
ugly, at least compared to slerp.

In fact, although computational mechanics and computer graphics have
much in common -- specially in what concerns the description of
rotations or general movement in 3D space --, the needs are not quite
exactly the same.

I really don't know much about splines, but it seems to me that their
principal advantages are in constructing curves smooth curves through
(or controlled by) an arbitrary number of points (let us call them
nodes). For what I understand, in most applications splines are
piecewise cubic polynomials, even when the number of nodes increases.
These piecewise polynomials have all the desired transition
properties, and I think that is the main purpose behind their design.
All this applies to both splines in 3D space and in SO(3).

On the finite element method, on the other hand, you look for
interpolations based on a number of n nodes, and you generally want
that the degree of the polynomials to be equal to n-1. In this way,
you can build increasingly refined approximations. In fact, you talk
about the h refinement (in which we refine the mesh and, accordingly
the number of elements) and the p refinement (in which you keep the
same number of elements, but increase the degree of tthe polynomials,
adding internal nodes to each element). In 3D flat space, nothing
better than a (n-1) degree Lagrange polynomial to interpolate between
n nodes. Note that, if the correct solution is a polynomial of degree
p, it suffices to have n=p+1 and you don't gain anything from

So, what I really want is the generalization of Lagrange polynomials
to SO(3) -- or, even better, to SE(3). I don't readily dismiss
splines, but I really don't think they are what I am after. For
instance, in flat space, you don't need splines at all in the finite
element method. So, why should they be the answer in a curved space?

Slerp can be seen as a straigth line (a polynomial of degree one) in
SO(3). So for a two node element, slerp is all I need. In particular,
it has all the uniform properties, not depending on the choice of a
starting point or "pole" (one of the suggestions I received, a paper
by Johnstone and Williams, doesn't seem to satisfy my uniform
requirement).

Let me elaborate on why uniformity is so important.
Last year I published a paper , in which I used the rotation vector
to parameterize the rotations, and implemented a finite element based
on the interpolation of this rotation vector. Now, I know that for
this representation I may bump into singularities, but for the
applications I have in mind, this is not a serious problem. On the
other hand, as recently shown (actually in 1999), this kind of
interpolation is not invariant under superimposed rotations, and THAT,
from the point of view of a geometrically exact formulation, is a
capital sin.

For what I understand, most of the interpolation schemes devised for
SO(3) rely on a interpolation on a flat space which is then projected
back to SO(3). That is not acceptable for my needs, since it suffers
from the exactly same problems that the vector-like parametrization
have (OK, it could avoid singularity).

Now, the paper of Buss and Fillmore, refered by you (Just), seems very
promising. In fact by working directly on SO(3), or other Sd spheres,
it seems to have all the uniformity -- and simple elegance -- I need.
Although the paper seems to concentrate in applications to splines, it
seems to me that the same spherical avarages can be used to build my
desired Lagrange-like interpolation in SO(3). So, although the details
are still beyond my understanding, I see its great potential.

I am finishing my PhD thesis and I don't have much time to develop
this idea right now. However, somehow, all the currently used
implementations of geometrical exact beam theory seemed to lack that
special final touch. My intent in coming to this newsgroup, was the
hope that some one from the computer graphic field had some insight
that was absent from the computational mechanics comunity. I am not
disappointed.

Thank you all, for your valuable help. And of course, if you have
furher comments, I will be delighted to hear them.

Manuel Correa

 "On the differentiation of the Rodrigues formula and its
significance for the vector-like parameterization of Reissner-Simo
beam theory" in International Journal for Numerical Methods in
Engineering, Volume 55, Issue 9, 2002, pp 1005-1032 (is someone is
interested, ask me and I will send a pdf).

### slerp between three or more points

>So, what I really want is the generalization of Lagrange polynomials
>to SO(3) -- or, even better, to SE(3). I don't readily dismiss
>splines, but I really don't think they are what I am after. For
>instance, in flat space, you don't need splines at all in the finite
>element method. So, why should they be the answer in a curved space?

You said quite a lot; I'll just comment on a little.

Dealing with translation vectors is so simple, people are surprised to
learn how much more is required to deal with rotations. Combining both
is worse yet, because of a fascinating technical detail. The Lie group
SE(3) naturally has one-parameter subgroups, generated by iterating an
infinitesimal transform. These are "screw motions", and intrinsic to
the algebra. Previous experience leads one to expect SE(3) also has a
geometry to match. But it does not! That is, there is no way to assign
distances so that geodesics correspond to one-parameter subgroups. So
there can be no geometric form of "generalized linear interpolation"
(Glerp?!), much less polynomial functions or splines. Here's a little
light reading of relevance.

<http://www.seas.upenn.edu/~calin/interp_papers/c03301.pdf>

I suspect Kumar, who is at Penn, indirectly learned about this from
the same source I did. I first heard about this years ago from Eloise
Carlton, a mathematician working on apparent motion with Stanford
psychologist Roger Shepard; somewhat later their work was published in
two articles.

Carlton, E.H. & Shepard, R.N. (1990a). Psychologically
simple motions as geodesic paths: I. Asymmetric objects.
Journal of Mathematical Psychology, 34, 127-88.

Carlton, E.H. & Shepard, R.N. (1990b). Psychologically
simple motions as geodesic paths: II. Symmetric objects.
Journal of Mathematical Psychology, 34, 189-228.

The common source is likely Peter Freyd, a category theorist at Penn
who apparently suggested the idea to Shepard.

I will confine my further comments to SO(3), or S^3, which is already
difficult enough.

Lagrange interpolation can be viewed in two different ways. The first
is as a set of weight polynomials used to add points. Under this
interpretation you could apply the Buss method. Alternatively, it is
possible to build Lagrange interpolants using nested linear layers, a
tradition known as divided differences. For this you could use Slerp
repeatedly. I cannot vouch for the latter's continuity, only hope.

Example.
Given q1, q2, q3, and q4 at t1, t2, t3, and t4 (with t1<t2<t3<t4), we
wish to interpolate. We proceed as follows:

q12 = Slerp(q1,q2; t:(t1,t2))
q23 = Slerp(q2,q3; t:(t2,t3))
q34 = Slerp(q3,q4; t:(t3,t4))

q13 = Slerp(q12,q23; t:(t1,t3))
q24 = Slerp(q23,q34; t:(t2,t4))

q14 = Slerp(q13,q24; t:(t1,t4))

The mysterious notation t:(ti,tj) means use a parameterization for the
Slerp so that it hits the first point at ti and the last point at tj.
However in most cases t will be beyond those bounds, thus we will be
extrapolating (but using the usual Slerp formula).

This interpolates, as is easily shown. Take t=t2. Then q12 is q2, as
is q23. We don't know or care about q34. Then q13 must also be q2, as
must q24. Thus q14 can only be q2, as desired. And so on.

### slerp between three or more points

Quote:> Dealing with translation vectors is so simple, people are surprised to
> learn how much more is required to deal with rotations. Combining both
> is worse yet, because of a fascinating technical detail. The Lie group
> SE(3) naturally has one-parameter subgroups, generated by iterating an
> infinitesimal transform. These are "screw motions", and intrinsic to
> the algebra. Previous experience leads one to expect SE(3) also has a
> geometry to match. But it does not! That is, there is no way to assign
> distances so that geodesics correspond to one-parameter subgroups. So
> there can be no geometric form of "generalized linear interpolation"
> (Glerp?!), much less polynomial functions or splines. Here's a little
> light reading of relevance.

>   <http://www.seas.upenn.edu/~calin/interp_papers/c03301.pdf>

Hmmm. I couldn't read this pdf (fonts missing, maybe I should upgrade
my acrobat reader) but I found enough papers by Kumar and co-authors
to get an idea. Well, I must say, I had the idea that it was easy to
generalize from SO(3) to SE(3). I see that apparently it is not.

However, given Chasles' theorem for rigid body motions, it seems easy
to interpolate "linearly" between two positions: just find the central
axis and interpolate linearly both the dispacement along the axis and
the rotation around that axis (slerp). In fact, it seems to me that
this was already done in a  finite element implementation by Borri and
Botasso (at Politecnico di Milano). Maybe, I am not grasping all the
subtleties...

Quote:> I suspect Kumar, who is at Penn, indirectly learned about this from
> the same source I did.

It seems to me that you want to keep your true identity hidden. You
have, of course, the right to proceed as you like, but I confess I am
really curious about who you really are. :-)

Anyway, those references to Carlton and Shepard work were interesting
(although I don't think it will be easy to obtain them -- after all I
am in an engineering university!).

Quote:> I will confine my further comments to SO(3), or S^3, which is already
> difficult enough.

> Lagrange interpolation can be viewed in two different ways. The first
> is as a set of weight polynomials used to add points. Under this
> interpretation you could apply the Buss method. Alternatively, it is
> possible to build Lagrange interpolants using nested linear layers, a
> tradition known as divided differences. For this you could use Slerp
> repeatedly. I cannot vouch for the latter's continuity, only hope.

Nice insight. Well, I will consider both alternatives. I might add
that, besides obtaining a smooth and uniform interpolation, I want
also to be able to linearize it (either to obtain strain measures or
to be able to apply the Newton iterative method). In particular, If
R1, R2 and R3 are the orthogonal tensors at the nodes located at
T=t1,t2 and t3, I want an interpolation R(t), its derivative
R'(t)=dR/dt and the variations dR and dR' written in terms of the
variations dR1, dR2 and dR3. Of course, all this may be written in
terms of quaternions. If these values can be computed in close form,
it is better, but it is also acceptable to obtain them numerically.

Thanks, once more for all the interesting and knowledgable feed-back.

Manuel Correa

Wasn't there a discussion here some time ago about how to find the
center of a sphere given three points on its boudary?  Could someone
please forward the solution to me.  Thanks!

--
_________     __________________
\\        \\  \\

\\________\\  \\_______________     NTH (Norwegian Institute of Technology)

"Some people play hard to get, I play hard to want" (Ford Fairlane)