Handedness

1,334

Author

June 11, 2009 01:55 AM

Would handedness affect quaternions, vectors and matrices? I am now porting an openGL column major based application to directx 9? Having some difficulties with that...For example


x = cross(z,up);
y = cross(z,x);
 


Vec cross ( const Vec& v1, const Vec& v2 )
 {
   return Vec ( v1.y*v2.z - v1.z*v2.y,
                v1.z*v2.x - v1.x*v2.z,
                v1.x*v2.y - v1.y*v2.x  );
 }

What I can interpret is the line x = cross(z,up) must be written as x = cross(up,z) for dx. Does it sound right? The whole bunch of source code is.... would you please help me correct the ones that are wrong? I gave up after looking at it for 2 nights. Thanks a lot



# include <math.h>
# include "math.h"

const float srpi = 3.14159265f;
const float srtiny = 1.0E-6f;
# define SRNEXTZ(a,eps) ( (a)>-(eps) && (a)<(eps) )

float Vec::len () const
 {
   float f = x*x + y*y + z*z;
   if ( f>0 ) f = sqrtf ( f );
   return f;
 }

float Vec::len ( float n )
 {
   float f = sqrtf (x*x + y*y + z*z);
   if ( f>0 ) { n/=f; x*=n; y*=n; z*=n; }
   return f;
 }

Vec cross ( const Vec& v1, const Vec& v2 )
 {
   return Vec ( v1.y*v2.z - v1.z*v2.y,
                v1.z*v2.x - v1.x*v2.z,
                v1.x*v2.y - v1.y*v2.x  );
 }

void Quat::rot ( const Vec& axis, float radians )
 { 
   float f;
   
   // normalize axis:
   x=axis.x; y=axis.y; z=axis.z;
   f = x*x + y*y + z*z;
   if ( f>0 )
    { f = sqrtf ( f );
      x/=f; y/=f; z/=f;
    }
    
   // set the quaternion:
   radians/=2;
   f = sinf ( radians );
   x*=f; y*=f; z*=f;
   w = cosf ( radians );
 }

// above ok jacky

Quat operator * ( const Quat &q1, const Quat &q2 )
 {
   Quat q;

   q.w = (q1.w*q2.w) - (q1.x*q2.x + q1.y*q2.y + q1.z*q2.z); // w1*w2-dot(v1,v2)
   q.x = q1.y*q2.z - q1.z*q2.y; // cross (q1.v,q2.v)
   q.y = q1.z*q2.x - q1.x*q2.z;
   q.z = q1.x*q2.y - q1.y*q2.x;
   q.x += (q1.x*q2.w) + (q2.x*q1.w);
   q.y += (q1.y*q2.w) + (q2.y*q1.w);
   q.z += (q1.z*q2.w) + (q2.z*q1.w);

   return q;
 }

void quat_to_ts ( const Quat& q, float& tw, float& sx, float& sy )
 {
   // First test if the swing is in the singularity:
   if ( SRNEXTZ(q.z,srtiny) && SRNEXTZ(q.w,srtiny) ) { sx=sy=srpi; tw=0; return; }

   // Decompose q into twist-swing
   // by solving the equation Qtwist(t*2) * Qswing(s*2) = q
   // note: (x,y) is the normalized swing axis (x*x+y*y=1)
   // ( Ct 0 0 St ) * ( Cs xSs ySs 0 ) = ( qw qx qy qz )
   // ( CtCs  xSsCt-yStSs  xStSs+ySsCt  StCs ) = ( qw qx qy qz )  (1)
   // From (1) CtCs / StCs = qw/qz => Ct/St = qw/qz => tan(t) = qz/qw (2)
   // The swing rotation/2 s comes from:
   // From (1) (CtCs)^2 + (StCs)^2 = qw^2 + qz^2 =>  Cs = sqrt ( qw^2 + qz^2 ) (3)
   // From (1) (xSsCt-yStSs)^2 + (xStSs+ySsCt)^2 = qx^2 + qy^2 => Ss = sqrt ( qx^2 + qy^2 ) (4)
   // From (1) : |SsCt -StSs| |x| = |qx|
   //            |StSs +SsCt| |y|   |qy| (5)

   double qw, qx, qy, qz;
   if ( q.w<0 )
    { qw=-q.w; qx=-q.x; qy=-q.y; qz=-q.z; }
   else
    { qw=q.w; qx=q.x; qy=q.y; qz=q.z; }
      
   double t = atan2 ( qz, qw ); // from (2)
   double s = atan2( sqrt(qx*qx+qy*qy), sqrt(qz*qz+qw*qw) ); // from (3) and (4)

   double x=0, y=0, sins=sin(s);

   if ( !SRNEXTZ(sins,srtiny) )
    { double sint = sin(t);
      double cost = cos(t);

      // by solving the linear system in (5):
      y = (-qx*sint + qy*cost)/sins;
      x = ( qx*cost + qy*sint)/sins;
    }

   tw = float(2.0*t);
   sx = float(x*2.0*s);
   sy = float(y*2.0*s);
 }

void quat_to_st ( const Quat& q, float& sx, float& sy, float& tw )
 {
   // Decompose q into swing-twist using a similar development as in
   // function quat_to_ts
   if ( SRNEXTZ(q.z,srtiny) && SRNEXTZ(q.w,srtiny) ) { sx=sy=srpi; tw=0; }

   double qw, qx, qy, qz;
   if ( q.w<0 )
    { qw=-q.w; qx=-q.x; qy=-q.y; qz=-q.z; }
   else
    { qw=q.w; qx=q.x; qy=q.y; qz=q.z; }

   // Get the twist t:
   double t = 2.0 * atan2(qz,qw);

   double bet = atan2( sqrt(qx*qx+qy*qy), sqrt(qz*qz+qw*qw) );
   double gam = t/2.0;
   double sinc = SRNEXTZ(bet,srtiny)? 1.0 : sin(bet)/bet;
   double singam = sin(gam);
   double cosgam = cos(gam);

   sx = float( (2.0/sinc) * (cosgam*qx - singam*qy) );
   sy = float( (2.0/sinc) * (singam*qx + cosgam*qy) );
   tw = float( t );
 }

Vec operator * ( const Vec &v, const Quat &q )
 {
   Quat qv ( 0, v.x, v.y, v.z );
   qv = q * qv * q.inverse();
   return Vec ( qv.x, qv.y, qv.z );
 }

void Mat::set ( const Quat& q, const Vec& v )
 {
   // convert quaternion to rotation matrix:
   float x2  = q.x+q.x;
   float x2x = x2*q.x;
   float x2y = x2*q.y;
   float x2z = x2*q.z;
   float x2w = x2*q.w;
   float y2  = q.y+q.y;
   float y2y = y2*q.y;
   float y2z = y2*q.z;
   float y2w = y2*q.w;
   float z2  = q.z+q.z;
   float z2z = z2*q.z;
   float z2w = z2*q.w;

   e[0] = 1.0f - y2y - z2z; e[1] = x2y + z2w;        e[2]  = x2z - y2w;        e[3]=0;
   e[4] = x2y - z2w;        e[5] = 1.0f - x2x - z2z; e[6]  = y2z + x2w;        e[7]=0;
   e[8] = x2z + y2w;        e[9] = y2z - x2w;        e[10] = 1.0f - x2x - y2y; e[11]=0;

   // now update translation:
   e[12] = v.x;
   e[13] = v.y;
   e[14] = v.z;
   e[15] = 1.0f;
 }

void Mat::lookat ( const Vec& eye, const Vec& center, const Vec& up )
 {
   Vec z = eye-center; z.normalize();
   Vec x = cross ( up, z );  // Jacky for left-handed system, Vec x = cross(z,up);
   Vec y = cross ( z,  x );  // Jacky this is correct for both systems

   x.normalize();
   y.normalize();

   e[0]=x.x; e[1]=y.x; e[2]=z.x; e[3]=0;  // jacky need transpose? it's a column vector where dx is row based
   e[4]=x.y; e[5]=y.y; e[6]=z.y; e[7]=0;
   e[8]=x.z; e[9]=y.z; e[10]=z.z; e[11]=0;
   e[15]=1.0f;

   Vec v = eye*-1.0f;
   // left combine translation:
   e[12] = v.x*e[0] + v.y*e[4] + v.z*e[8];
   e[13] = v.x*e[1] + v.y*e[5] + v.z*e[9];
   e[14] = v.x*e[2] + v.y*e[6] + v.z*e[10];
 }

# define E11 e[0]
# define E12 e[1]
# define E13 e[2]
# define E14 e[3]
# define E21 e[4]
# define E22 e[5]
# define E23 e[6]
# define E24 e[7]
# define E31 e[8]
# define E32 e[9]
# define E33 e[10]
# define E34 e[11]
# define E41 e[12]
# define E42 e[13]
# define E43 e[14]
# define E44 e[15]

void Mat::inverse ( Mat &inv ) const
 {
   float m12 = E21*E32 - E22*E31;
   float m13 = E21*E33 - E23*E31;
   float m14 = E21*E34 - E24*E31;
   float m23 = E22*E33 - E23*E32;
   float m24 = E22*E34 - E24*E32;
   float m34 = E23*E34 - E24*E33;

   float d1 = E12*m34 - E13*m24 + E14*m23;
   float d2 = E11*m34 - E13*m14 + E14*m13;
   float d3 = E11*m24 - E12*m14 + E14*m12;
   float d4 = E11*m23 - E12*m13 + E13*m12;
   float d = -E41*d1 + E42*d2 - E43*d3 + E44*d4;
   if (d==0.0) return;
   d = 1.0f/d;

   inv.E11 = (E42*m34 - E43*m24 + E44*m23) * d;
   inv.E21 = (E43*m14 - E41*m34 - E44*m13) * d;
   inv.E31 = (E41*m24 - E42*m14 + E44*m12) * d;
   inv.E41 = (E42*m13 - E41*m23 - E43*m12) * d;
   inv.E14 = (E13*m24 - E12*m34 - E14*m23) * d;
   inv.E24 = (E11*m34 - E13*m14 + E14*m13) * d;
   inv.E34 = (E12*m14 - E11*m24 - E14*m12) * d;
   inv.E44 = (E11*m23 - E12*m13 + E13*m12) * d;

   m12 = E11*E42 - E12*E41;
   m13 = E11*E43 - E13*E41;
   m14 = E11*E44 - E14*E41;
   m23 = E12*E43 - E13*E42;
   m24 = E12*E44 - E14*E42;
   m34 = E13*E44 - E14*E43;
   inv.E12 = (E32*m34 - E33*m24 + E34*m23) * d;
   inv.E22 = (E33*m14 - E31*m34 - E34*m13) * d;
   inv.E32 = (E31*m24 - E32*m14 + E34*m12) * d;
   inv.E42 = (E32*m13 - E31*m23 - E33*m12) * d;
   inv.E13 = (E23*m24 - E22*m34 - E24*m23) * d;
   inv.E23 = (E21*m34 - E23*m14 + E24*m13) * d;
   inv.E33 = (E22*m14 - E21*m24 - E24*m12) * d;
   inv.E43 = (E21*m23 - E22*m13 + E23*m12) * d;
 }

Vec operator * ( const Vec &v, const Mat &m )
 {
   Vec r ( m.E11*v.x + m.E21*v.y + m.E31*v.z + m.E41,
             m.E12*v.x + m.E22*v.y + m.E32*v.z + m.E42,
             m.E13*v.x + m.E23*v.y + m.E33*v.z + m.E43  );

   float w = m.E14*v.x + m.E24*v.y + m.E34*v.z + m.E44;
   if ( w!=0.0 && w!=1.0 ) { r.x/=w; r.y/=w; r.z/=w; }
   return r;
 }

Vec operator * ( const Mat &m, const Vec &v )
 {
   Vec r ( m.E11*v.x + m.E12*v.y + m.E13*v.z + m.E14,
             m.E21*v.x + m.E22*v.y + m.E23*v.z + m.E24,
             m.E31*v.x + m.E32*v.y + m.E33*v.z + m.E34  );

   float w = m.E41*v.x + m.E42*v.y + m.E43*v.z + m.E44;
   if ( w!=0.0 && w!=1.0 ) { r.x/=w; r.y/=w; r.z/=w; }
   return r;
 }

Thanks Jack [Edited by - lucky6969b on June 11, 2009 3:37:58 AM]

lucky6969b

1,334

Author

June 11, 2009 04:46 AM

I am having difficulties understanding the inverse function done by the author.
Can I assume that he was using pivoting for a 4x4 matrix?
I am cross-checking it with a reference, but not too sure how he did that.
Thanks
Jack

RandomBystander

122

June 11, 2009 05:00 AM

Quote:Original post by lucky6969b
I am having difficulties understanding the inverse function done by the author.
Can I assume that he was using pivoting for a 4x4 matrix?
I am cross-checking it with a reference, but not too sure how he did that.
Thanks
Jack

Try to look up Cramer's_rule, that's what he seems to be doing. There is probably no pivoting going on, only computation of several determinants.

lucky6969b

1,334

Author

June 11, 2009 05:19 AM

Quote:Original post by RandomBystander
Quote:Original post by lucky6969b
I am having difficulties understanding the inverse function done by the author.
Can I assume that he was using pivoting for a 4x4 matrix?
I am cross-checking it with a reference, but not too sure how he did that.
Thanks
Jack

Try to look up Cramer's_rule, that's what he seems to be doing. There is probably no pivoting going on, only computation of several determinants.

Thanks for your reply.
I found this on the web
http://en.wikipedia.org/wiki/Cramer%27s_rule#Finding_inverse_matrix
However, still not understand it... can you explain it with more details?

Like this very first section, what does he want to achieve?

float m12 = E21*E32 - E22*E31;   float m13 = E21*E33 - E23*E31;   float m14 = E21*E34 - E24*E31;   float m23 = E22*E33 - E23*E32;   float m24 = E22*E34 - E24*E32;   float m34 = E23*E34 - E24*E33;

Thanks
Jack

Zakwayda

2,389

June 11, 2009 06:59 AM

Just a quick comment on the 'handedness' issue. Generally speaking, the (geometric) handedness of the coordinate system used should not affect your math library code (by which I mean code that constructs or manipulates vectors, matrices, quaternions, and so on). In other words (for example) a function that builds a rotation matrix or quaternion or computes the cross product of two vectors should be exactly the same regardless of the handedness of the coordinate system being used.

Projection transforms and 'view' transforms (e.g. a 'look at' transform) are exceptions in that they're generally constructed differently for left- and right-handed systems. The differences actually have to do with which direction is considered 'forward' in view space, but are also an indirect result of the coordinate system handedness and the convention that the positive x axis should point to the right in view space.

As for computing cross products, with a typical orthonormal basis the relationships between the vectors should be:

x = yXzy = zXxz = xXy

If you're seeing the arguments flipped around depending on handedness, you might be looking at a 'look at' function (where the flipping is related to which direction is to be considered forward in view space).

Even though most functions should be unaffected by handedness, you may find that the results of applying said functions change when switching handedness. For example, triangle windings may flip, models may be mirrored, and rotations may appear to go the 'wrong' way. In these cases though it's the input data that should be changed, not the functions that work with the data.

Also, keep in mind that the DirectX/D3D and OpenGL APIs differ in a few other ways as well, most notably matrix layout (row major vs. column major), vector notation (row vs. column vectors) and the near plane distance for the canonical view volume (zero vs. negative one). The issues of matrix layout, vector notation, and coordinate system handedness are often confused with each other, but in fact they are three separate and unrelated issues (although they can interact with each other in ways that can be quite confusing).

lucky6969b

1,334

Author

June 11, 2009 11:23 PM

So the conclusion is I don't need to change anything..
Is that right? :)
Thanks
Jack

Zakwayda

2,389

June 11, 2009 11:49 PM

Quote:So the conclusion is I don't need to change anything..
Is that right? :)

No, not necessarily. As I mentioned earlier, you may have to make adjustments for the differing view space and projection transform conventions, and you may also have to make adjustments to your data so that the visual results are consistent between the two APIs.

There also may be cases where you'll need to take matrix layout and/or vector notation conventions into account (note however that despite the different conventions used, D3D and OpenGL matrices are actually arranged in memory the same way in that the elements of the basis vectors are stored contiguously).

lucky6969b

1,334

Author

June 12, 2009 12:31 AM

Quote:Original post by jyk
Quote:So the conclusion is I don't need to change anything..
Is that right? :)
No, not necessarily. As I mentioned earlier, you may have to make adjustments for the differing view space and projection transform conventions, and you may also have to make adjustments to your data so that the visual results are consistent between the two APIs.

There also may be cases where you'll need to take matrix layout and/or vector notation conventions into account (note however that despite the different conventions used, D3D and OpenGL matrices are actually arranged in memory the same way in that the elements of the basis vectors are stored contiguously).

Thanks jyk,
Actually, at first blush, I am most concerned about the code that I post would work or not... :) sorry for making you repeat again.
You really gave me a lot of good info.... but would you mind checking the code for me? I just really need to know if the source code would work on Direct3D or not....
Thanks
Jack

Zakwayda

2,389

June 12, 2009 02:27 AM

Quote:You really gave me a lot of good info.... but would you mind checking the code for me? I just really need to know if the source code would work on Direct3D or not....

I'd certainly like to be able to help, but it would take me a long time to proof the code that you posted with any real accuracy. (I did spot a couple of things though - for example, it looks to me like you might indeed need to transpose the upper-left 3x3 portion of your 'look at' matrix, as the comment indicates.)

IMO, the best way to confirm the correctness of this sort of code is to gain a solid understanding of the underlying algorithms so that you can proof the implementations yourself. (You can also use reliable references, such as books, articles, or source code that you know works correctly, as a sanity check.)

Here are a couple of tips:

1. Make your code as clean and well-formated as possible. Avoid unclear and unnecessarily terse variable names (for example, there's little reason to name a variable bet instead of beta), and don't cram lots of code onto one line (for example, don't put an 'if' statement and the block of code associated with it all on the same line). The code you posted isn't bad, by the way (I've seen a lot worse), but it could still use some cleaning up.

2. Use 2-d indexing for matrix element access rather than 1-d indexing. Storage order is an implementation detail that should be kept out of the actual 'math' code. (That is, when you're writing a function to build a rotation matrix, you shouldn't be thinking about whether the matrix is row- or column-major, since that has nothing to do with the math.) Using 1-d indexing also makes the code harder to proofread (for yourself and for others) since you have to convert the indices in your head as you go, and because you have to know what the storage order is in order to know which elements are being accessed. (Actually, you're already doing this with your E** macros, but only in some places. Also, I'd recommend using an inline function rather than a macro.)

(The above is all IMO, by the way.)

If there are particular functions you're unsure of, perhaps you could clean them up and post them individually; that might make the task of proofing them seem a little less daunting :)

Steve_Segreto

2,094

June 27, 2009 10:53 PM

Whoa! That's some ugly code.

If you are porting that to DirectX 9.0 here is my recommendation:

*DO NOT* port the math functions directly, instead use D3DX.

*DO* port the data types from your system to DirectX's D3DX system.

So don't use the Vec, Quat and Mat classes methods and operators. Whenever you encounter data of these types convert them to these types:

D3DXVECTOR3, D3DXQUATERNION and D3DXMATRIX

then use the built-in D3DX math library functions (see the DirectX 9.0 SDK)

D3DXVec3Cross(), D3DXVec3Dot(), D3DXMatrixMultiply(), D3DXMatrixInverse(), D3DXMatrixTranspose(), D3DXQuaternionMultiply(), D3DXMatrixLookAtLH(), etc, etc..

Now to convert the original Vec, Quat and Mat data to DirectX formats, you will probably need the following (I'm assuming the native types were OpenGL right-handed):

D3DXVECTOR3 VecToD3DXVECTOR3( Vec& v )
{
return (D3DXVECTOR3( v.x, v.z, v.y ) );
}

D3DXQUATERNION QuatToD3DXQUATERNION( Quat &q )
{
return (D3DXQUATERNION( q.x, q.z, q.y, q.w ) );
}

D3DXMATRIX MatToD3DXMATRIX( Mat &m )
{
D3DXMATRIX outM;

//
// Right Vector
//
outM( 0, 0 ) = m.e[ 0 ];
outM( 0, 1 ) = m.e[ 2 ];
outM( 0, 2 ) = m.e[ 1 ];
outM( 0, 3 ) = m.e[ 3 ];

//
// Up Vector
//
outM( 1, 0 ) = m.e[ 8 ];
outM( 1, 1 ) = m.e[ 10 ];
outM( 1, 2 ) = m.e[ 9 ];
outM( 1, 3 ) = m.e[ 11 ];

//
// Look Vector
//
outM( 2, 0 ) = m.e[ 4 ];
outM( 2, 1 ) = m.e[ 6 ];
outM( 2, 2 ) = m.e[ 5 ];
outM( 2, 3 ) = m.e[ 7 ];

//
// Position Vector
//
outM( 3, 0 ) = m.e[ 12 ];
outM( 3, 1 ) = m.e[ 14 ];
outM( 3, 2 ) = m.e[ 13 ];
outM( 3, 3 ) = m.e[ 15 ];

return (outM);
}

Here's some explanation:

Let's say a left-handed DirectX 4x4 looks like this:

{ rx, ry, rz, 0 }
{ ux, uy, uz, 0 }
{ lx, ly, lz, 0 }
{ px, py, pz, 1 }

where the vector u is the up vector { 0, 1, 0 } and the vector l is the look vector that points into the plane of the screen and the vector r is the cross product of u x l, and of course p is the position vector in the world.

Assuming openGL 4x4's use this same naming system, here's what I did with them to get them to display correctly under DirectX:

{ rx, rz, ry, 0 }
{ lx, lz, ly, 0 }
{ ux, uz, uy, 0 }
{ px, pz, py, 1 }

In openGL the y and z coordinates are transposed (i.e. for them z is up/down and y is into the monitor), which also means the up and look vectors are transposed in the 4x4.

For quaternions I had to do the following:

DirectX:

{ x, y, z, w } or { ._41, ._42, ._43, -w } if taken from a 4x4 converted from openGL to directX and then converted to a quaternion using D3DXQuaternionRotationMatrix().

openGL:

{ x, z, y, w }

!! IMPORTANT !!

For triangles, I had to change the winding order from {v1, v2, v3} to
{v3, v2, v1}.

For (u, v) coordinates I left them the same as I had them before.

Hope this helps!

Love DAOC? Tryout my DAOC clone: https://dl.dropboxusercontent.com/u/8974528/VON_Dist.zip

Handedness

This topic is closed to new replies.

Popular Topics

Recommended Tutorials

Handedness

This topic is closed to new replies.

Popular Topics

Recommended Tutorials

Reticulating splines