Re: [boost] [linear algebra] - question to David Abrahams

From: Janek Kozicki (janek_listy_at_[hidden])
Date: 2006-08-30 17:16:16

David Abrahams said:     (by the date of Wed, 30 Aug 2006 15:20:33 -0400)

> Janek Kozicki <janek_listy_at_[hidden]> writes:
> 
> > Hello,
> >
> > Not so long ago you have mentioned that you are working on a linear
> > algebra library. I'm very interested in this topic, because I'm using
> > such library in my project.
> 
> I'm sorry I've been so quiet...

Hello, much better late than never :)) Many thanks for your answer.

> > I have looked into the vault but couldn't find this library. Can we have
> > a sneak peek?
> 
> I'm sorry, not yet.  We have lots of code, but nothing worth showing
> at the linear algebra level.

ok. I have found that it is MTL4, mentioned here:

http://osl.iu.edu/research/mtl/

> > Questions:
> > 1. will your library provide cross_product(), etc... ?
>
> Yes, and it will minimize the amount of temporary storage used, so it
> should be very fast.

great! :)
 
> > 2. will your library allow Boost.Quaternions to be extended in a way that
> > those two libraries together will be useful for geometrical rotations?
> 
> Probably.  Well-designed generic libraries usually use concepts that
> are not dependent on member functions, so it should be possible to
> write the traits and other components that make that work.

then I hope that upon adding MTL4 to boost such functions (described
below in more detail) will be added too.

> > 3. will your library provide a way to solve equation Ax=b ?
> 
> Yes. 

great again :)

> > The background for question 2 above:
> > Following functions would be needed:
> >
> > operator*(vector<3,T1>,quaternion<T2>)  // most important
> 
> Is that an inner product, outer product, or pointwise multiplication?

None of above. This function operator* performs a rotation of three
element vector by quaternion. It is possible that you will prefer a
function name "rotate" instead of "operator* ".

(Although operator* is very convenient for people who work with
quaternions).

The rotation is done by first building a rotation matrix from
quaternion, and then multiplying vector by this matrix. See:

http://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation

Scroll down this webpage to the bottom, in the section titled
"Quaternions versus other representations of rotations" you can see that
matrix mentioned above.

The code to do that may look like this:

vector<3,T> operator* (const quaternion<T>& lhs,const vector<3,T>& rhs)
{
   T    aa = lhs.w*lhs.w;
   T    ab = lhs.w*lhs.x;
   T    ac = lhs.w*lhs.y;
   T    ad = lhs.w*lhs.z;
   T    bb = lhs.x*lhs.x;
   T    bc = lhs.x*lhs.y;
   T    bd = lhs.x*lhs.z;
   T    cc = lhs.y*lhs.y;
   T    cd = lhs.y*lhs.z;
   T    dd = lhs.z*lhs.z;
   T    mag = aa+bb+cc+dd; // [1]

// build matrix elements (possibly building an MTL4::matrix<3,3,T> will be a better coding here?)

   T at11 = (aa+bb-cc-dd)/mag;
   T at12 = 2*(-ad+bc)/mag;
   T at13 = 2*(ac+bd)/mag;
   T at21 = 2*(ad+bc)/mag;
   T at22 = (aa-bb+cc-dd)/mag;
   T at23 = 2*(-ab+cd)/mag;
   T at31 = 2*(-ac+bd)/mag;
   T at32 = 2*(ab+cd)/mag;
   T at33 = (aa-bb-cc+dd)/mag;

// returns rotated vector

   return vector<3,T> ( rhs.x*at11 + rhs.y*at12 + rhs.z*at13
                       ,rhs.x*at21 + rhs.y*at22 + rhs.z*at23 
                       ,rhs.x*at31 + rhs.y*at32 + rhs.z*at33 );
}

[1] this is only to ensure that quaternion has magnitude==1, in that
wikipedia article this denominator is not mentioned, because it is
assumed that quaternion already has magnitude==1. This assumption could
also be made in this function, and then variable 'mag' could be removed.

And when this function will be called with quaternion of magnitude!=1,
then results will be invalid.

> > from_rotation_matrix(matrix<3,3,T>)
> 
> returning what?

I'm sorry for being not specific.

This function is supopsed to create a quaternion when given a matrix<3,3,T>.
This matrix is supposed to be a rotation matrix. Function signature:

quaternion<T> from_rotation_matrix(matrix<3,3,T>);

sample implementation starts from line 382 in file Quaternion.h from boost vault:

http://www.boost-consulting.com/vault/index.php?directory=Math%20-%20Numerics

Although I don't know why the argument to this function is a 4x4 matrix
in this concrete implementation. Perhaps I should investigate more.

> > to_rotation_matrix()
> 
> Where's the argument?  Is this supposed to be a member function?  See
> my 4th paragraph above.

Not necessarily a member function, if not then the function signature would look like:

matrix<3,3,T> to_rotation_matrix(quaternion<T>);

This function is supposed to create a rotation matrix when given a
quaternion. Sample implementation is above in this email (part of
rotation operator* code), or starting from line 327 from Quaternion.h in
boost vault (please note that there the returned matrix is 4x4 but its
last row and column are neutral => unnecessary).

 
> > from_axis_angle(vector<3,T1> axis, T2 angle)
> 
> returning what?

Function signature:

quaternion<T> from_axis_angle(vector<3,T> axis, T angle);

This function is suposed to create a quaternion when given an axis of
rotation and angle of rotation around this axis, sample implementation
may look like:

quaternion<T> from_axis_angle(vector<3,T> axis, T angle)
{
   T half_angle = angle*0.5;
   vector<3,T> sin_axis = sin(half_angle)*axis;
   return quaternion<T>( cos(half_angle)
                        ,sin_axis.x
                        ,sin_axis.y
                        ,sin_axis.z);
}

Or look at line 305 of file Quaternion.h from boost vault.

> > to_axis_angle()
> 
> Where's the argument?

This function does the opposite, returns axis and angle when given a
quaternion. Function signature:

  void to_axis_angle(vector<3,T>& axis, T& angle, quaternion<T>);

or alternatively, without taking references to existing variables:

  std::pair<vector<3,T>,T> to_axis_angle(quaternion<T>);

This function is mathematically opposite of from_axis_angle, and
involves arcus cosinus :)

 
> > to_euler_angles()
> 
> Where's the argument?

This function takes as argument a quaternion and converts it into euler
angle representation of rotation. Function signature:

vector<3,T> to_euler_angles(quaternion<T>);

Function implementation, for example from:

http://en.wikipedia.org/wiki/Conversion_between_quaternions_and_Euler_angles

> > and few others. Similar functionality is already in boost vault:
> > directory "Math - Numerics", file Quaternion.h

just a reminder:

http://www.boost-consulting.com/vault/index.php?directory=Math%20-%20Numerics

> > Without quaternion rotations this linear algebra library would be
> > pretty useless for me, so I hope it'll be there :)
> 
> OK.

I hope that it is more clear now :) I'll be glad to explain more if
necessary (and help with coding too, according to my abilities ;).

-- 
Janek Kozicki                                                         |