Re: [eigen] Eigen and rigid body simulation
• To: eigen@xxxxxxxxxxxxxxxxxxx
• Subject: Re: [eigen] Eigen and rigid body simulation
• From: Gael Guennebaud <gael.guennebaud@xxxxxxxxx>
• Date: Wed, 2 Dec 2009 11:31:03 +0100
• Domainkey-signature: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:in-reply-to:references:date:message-id:subject:from:to :content-type; b=oOQ8oxwrP4W0rt+jHAeia47nFv7J8NtdoiEAK872NeD55CydVBOrI4PsVo0vWbtmJp iSAu3Rxu65NrOGBvS+oa7TpXAJAMJ+UHBUNgaqHK5L+p5PXfeU/7L02cyE3EqkElucW3 RgfSYOk9o/hpfjuZgmXBYrWCPAEAAPh/u3Veo=

I also like Maxime's idea, but I'm wondering whether the Lie<G> group class should not be a decorator to G rather than a traits class. The main raisons for that are 1) to make it less tedious to use, 2) the current design of the geometry module already goes to that direction since all our transformation classes already implement inverse, identity and composition via operator* (e.g., Translation::operator* actually performs a +).

So basically, a Lie object should implement the following interface:

template<class G>
struct Lie {
typedef some_type Algebra;
typedef some_type CoAlgebra;

static G Identity();
G inverse() const;
G operator*(const G&) const;

G exp() const;
G log() const;
};

and then, e.g., Translation<S,D> could inherit Lie<Translation<S,D> >, NonUniformScaling<S,D> would inherit Lie<NonUniformScaling<S,D> >, Quaternion<S> would inherit Lie<Quaternion<S>>, etc.

Each of these specializations of Lie<> would also provide the storage and respective ctors. This way one can still directly instanciate a Lie<Translation<S,D> > if he really wants a "pure" Lie group object, and more importantly, write generic algorithms taking only objects implementing the Lie group concept:

template<G>
Lie<G> fast_exp(const Lie<G>& x, int n)
{
if( n == 0 )
return Lie<G>::Identity();

if( n % 2 )
return x * fast_exp(g*g, n/2);
else
return fast_exp( g*g, n/2);
}

However, the exp/log functions could now become ambiguous, especially for Quaternion. So if you think the above proposal makes sens, we could use lieExp/lieLog to make it more explicit they belong to the Lie algebra.

Now I have to acknowledge that your design has so advantages too. In particular it really well separate the Lie algebra, and make it clear when you are using operations belonging to the Lie algebra. In that sense it looks a bit cleaner. Are there any other advantages that I'm missing ?

Also I don't understand why do you use functors for the log/exp functions and not static functions ?

gael.

On Tue, Dec 1, 2009 at 11:24 PM, Benoit Jacob wrote:
2009/11/30 Maxime Tournier <maxime.tournier@xxxxxxxxxxxx>:
> Hi everyone,
>
> As Mathieu said, I have some generic code that allows one to abstract
> the Lie group structure of a type, so that one can write any Lie group
> algorithm no matter the group. I told Gael about it few days ago, and
> since it could be of interest to others, here's an overview of what I
> did.
>
> The approach somehow follows c++0x's concepts proposal, in that the
> Lie group structure for a type T is encoded in a separate template
> class, Lie<T>, containing all the types and operations that define the
> Lie group structure.
>
> One can describe the Lie group structure for a type T by specializing
> the template Lie for this type, giving the associated types and
> implementing required operations.  Every specialization should conform
> to the same "interface".
>
> Here are two basic examples with only the algebraic structure:
>
> // here is S^3 with quaternion multiplication
> template<class Real>
> struct Lie< Quaternion<Real> > {
>
>  static Quaternion<Real> id() { return Quaternion<Real>::Identity(); }
>  static Quaternion<Real> inv(const Quaternion<Real>& q) { return q.inverse();
> }
>  static Quaternion<Real> comp(const Quaternion<Real>& q1,
>                               const Quaternion<Real>& q2) {
>    return Quaternion<Real>::Identity();
>  }
>
> };
>
> // here is R^3 with addition
> template<>
> struct Lie< Vector3d > {
>
>  static Vector3d id() { return Vector3d::Zero(); }
>  static Vector3d inv(const Vector3d& v) {
>    return -v;
>  }
>  static Vector3d comp(const Vector3d& a,
>                       const Vector3d& b) {
>    return a + b;
>  }
>
> };
>
> You can then write an algorithm for any kind of group like this:
>
> template<class G>
> G fast_exponentiation(const G& g, int n) {
>  if( n == 0 ) {
>    return Lie<G>::id();
>  }
>
>  if( n % 2 ) {
>    return Lie<G>::comp(  fast_exponentiation( Lie<G>::comp(g, g), n/2), g );
>  } else {
>    return fast_exponentiation( Lie<G>::comp(g, g), n/2) );
>  }
>
> }
>
> Now I know it's a bit tedious to write, but:
>
> - No modification to exisiting types is needed whatsoever. This is A
>  Good Thing. No curious inheritance pattern is introduced in
>  particular. Actually, underlying types are completely independant
>  from this new structure. This should make Eigen types developers
>  happy since we're not touching anything in existing Eigen code :-)
>
> - We can quite easily construct new groups from existing ones:
>
> template<class G1, class G2>
> struct Lie< std::pair<G1, G2> > {
>
>  // ...
>
> };
>
> ...which essentialy says that the product of two groups is a group for
> the product law. This is particularly useful when using std::tuple and
> variadic templates, since product-groups are automatically constructed
> at compilation.

I understand, your design seems very good.

> The complete Lie template specification in my code is the following:
>
> template<class G>
> struct Lie {
>  typedef some_type algebra;
>  typedef some_type coalgebra;
>
>  static G id();
>  static G inv(const G&);
>  static G comp(const G&, const G& );
>
>  // default-constructible
>  typedef some_functor exp;
>  typedef some_functor log;
>
>  // G-constructible
>
> };

OK I understand. ad is inner automorphisms? and ad_T is the inverse?

>
> The only problem with it is the identity for dynamic-sized types, for
> example VectorXd: what should be the size of the result ?

Do you mean as a compile time value or as a runtime value?

At compile time:

Throughout Eigen we use a special value, Dynamic, to mean that. For
example, VectorXd is just:
matrix<double,Dynamic,1>.
Is that what you were looking for?

At run time: in the dynamic-size case, you could require id() to take
a size argument. So id(void) would use a static assertion that the
size is fixed at compile time. See what is already done in
MatrixBase::Identity().

> So maybe it
> should accept some optional contextual hint for the dimension.

in Eigen, all types inheriting MatrixBase (hence VectorXd) have a
member enum SizeAtCompileTime telling the number of coefficients
(rows*cols)

so:
VectorXd::SizeAtCompileTime == Dynamic

Other types could be harmonized if that is useful to you, e.g. we
could set Quaternion::SizeAtCompileTime==4 if that's useful.

> Once this is done, it is quite straightforward to implement the
> tangent bundle, then differentiable functions, and then automatically
> compute differentials for complicated function compositions using
> _expression_ templates. When dealing with complicated functions on
> complicated non-flat spaces, this is a real joy to use.
>
> I have the Lie group structure implemented for quaternions, Rohit
> Garg's rigid transforms, row/column vectors, as well as any std::tuple
> combination of these. Based on this I have code that can interpolate
> (linear/~spline), compute geometric mean, and do statistical analysis
> while preserving the Lie group (hence manifold) structure.
>
> So if you need to spline interpolate homographies
> in a non degenerate way for example, all you have to do is to
> implement the Lie group structure for homographies and the algorithm
> is already coded for you :-)
>
> So my question is: would anyone be interested in me porting all/parts
> of this stuff into Eigen ? If so where/how should I start ?

It seems that this would be nice content for a new module, "Lie". Does
that sound good?

You would start by forking the eigen repository at bitbucket, and
creating a new Lie module in unsupported/.

If you prefer it's always possible to develop directly in our
repository, it's just that forking gives you that extra flexibility.

I would also like to understand how this interplays with Mathieu's
proposal. If you create a new Lie module, should Mathieu's classes go
into it? Or should only Torque/Twist/Wrench go into it while "so(3)"
is considered so important for geometry that it goes into the Geometry
module? If Mathieu's stuff gets split between Geometry and Lie, does
that mean that Lie should #include Geometry? Is that OK? Is that
something you wanted to do anyway?

Benoit

 Mail converted by MHonArc 2.6.19+ http://listengine.tuxfamily.org/