Thread: rotation of vector on diagonal plane

I have a right hand 3d coordinate system, and a vector at [1,1,0]. Using rotations about the x y and z axis, I want to rotate this vector in an orbit about an axis orthogonal to the vector. So that a half rotation of the vector puts it at [-1,-1,0]. In other words, a "diagonal orbit" about the origin. However, I need to describe the path of the vector using x y and z rotations.
I'm having trouble coming up with the correct solution...

Thanks

Rick

Any rotation matrix about the origin can be expressed as the product of three rotations, one about each axis. So what you are asking is certainly possible. But I am not quite sure why you want to express the rotation like this : it is nothing like so meaningful as finding the axis and the amount of rotation, and the calculation is messy.
I presume you know how to calculate the rotation matrix using conjugation: i.e. pick a new orthonormal coordinate system with one axis the axis of rotation and one the vector. The matrix for the rotation is easy in the new coordinate system, and then you use the usual formula for change of basis to calculate the matrix in the original coordinate system. Given your vector is nice this should be an easy calculation
Then you want to decompose this matrix into rotation about the x,y,z axes. Obviously you get a different answer depending on what order you do things in. The answer is not usually unique - you can often get to the final position in more than one way.

I wrote some C code to do this calculation for a general rotation matrix. I think I decomposed on the basis of rotating first about x,then y,then z. I worked things out by calculating the product of three rotations and then used trig identities to see how I could recover the original values.

To understand the code:
transform is my 3*3 matrix, and I address it as a vector, so transform[0] is the first entry in the first row, transform[5] is the last entry on the second row.

Code:

y_rot = asin(transform[6]);
if (transform[7] || transform[8])
{
  x_rot = -atan2(transform[7],transform[8]);
  z_rot = -atan2(transform[3],transform[0]);
  if (x_rot < 0 && z_rot < 0)
  {
     y_rot = pi-y_rot;
     x_rot = atan2(transform[7],-transform[8]);
     z_rot = atan2(transform[3],-transform[0]);
  }
}
else
{
  x_rot = atan2(transform[1],transform[4]);
  if (transform[6] < 0)
     x_rot = -x_rot;
   z_rot = 0;
}

atan2 calculates an angle between -pi and pi radians. atan2(y,x) is arctan(y/x) + some multiple of half pi chosen based on the quadrant implied by the sign of y and x.

I should warn you that I do everything with left hand coordinates because that is more natural for my application. However, I don't think that will make any difference to the calculation at all.

Originally Posted by alunw

Any rotation matrix about the origin can be expressed as the product of three rotations, one about each axis. So what you are asking is certainly possible. But I am not quite sure why you want to express the rotation like this : it is nothing like so meaningful as finding the axis and the amount of rotation, and the calculation is messy.
I presume you know how to calculate the rotation matrix using conjugation: i.e. pick a new orthonormal coordinate system with one axis the axis of rotation and one the vector. The matrix for the rotation is easy in the new coordinate system, and then you use the usual formula for change of basis to calculate the matrix in the original coordinate system. Given your vector is nice this should be an easy calculation
Then you want to decompose this matrix into rotation about the x,y,z axes. Obviously you get a different answer depending on what order you do things in. The answer is not usually unique - you can often get to the final position in more than one way.

I wrote some C code to do this calculation for a general rotation matrix. I think I decomposed on the basis of rotating first about x,then y,then z. I worked things out by calculating the product of three rotations and then used trig identities to see how I could recover the original values.

To understand the code:
transform is my 3*3 matrix, and I address it as a vector, so transform[0] is the first entry in the first row, transform[5] is the last entry on the second row.

Code:

y_rot = asin(transform[6]);
if (transform[7] || transform[8])
{
  x_rot = -atan2(transform[7],transform[8]);
  z_rot = -atan2(transform[3],transform[0]);
  if (x_rot < 0 && z_rot < 0)
  {
     y_rot = pi-light->y_rot;
     x_rot = atan2(transform[7],-transform[8]);
     z_rot = atan2(transform[3],-transform[0]);
  }
}
else
{
  x_rot = atan2(transform[1],transform[4]);
  if (transform[6] < 0)
     x_rot = -x_rot;
   z_rot = 0;
}

atan2 calculates an angle between -pi and pi radians. atan2(y,x) is arctan(y/x) + some multiple of half pi chosen based on the quadrant implied by the sign of y and x.

I should warn you that I do everything with left hand coordinates because that is more natural for my application. However, I don't think that will make any difference to the calculation at all.

Thanks for the reply My limitations are because of hardware. I need to tell a robotic device to move to positions based on x y z, and the hardware is in a fixed position, but the position of the object it needs to interact with can be at an angle in relation to the hardware.

I do not know exactly (or at least don't recall) how to get the new rotation matrix based on conjugation. If you can, please explain, and I will also try to look up.

THanks

R.

Here is an example of what I mean. Your original vector is (1,1,0). Lets call that y, and your axis of rotation is orthogonal to this. Lets call that z. To make things easy let's assume it is (0,0,1). Now y cross z which I make to be (1,-1,0), but I suppose could be (-1,1,0) depending on your sign convention for the cross product is a third vector orthogonal to both the other two which we can call x. (In my programs I use an LHS coordinate system, and picked the sign of the cross product so that (1,0,0) cross (0,1,0) is (0,0,1). So I think what I'm doing would work just as well in an ordinary RHS system)
So if we normalise the vectors to length 1 by dividing the first and third by sqrt(2) we have a new orthonormal coordinate system,(A rotation by theta about z looks like , so it is very easy to express any rotation. Let us call this matrix M
Now the matrix that tells us how to get from our new coordinates to our original coordinates is simply [x,y,z] i.e. . Call that matrix N
So the rotation matrix for a rotation of theta in our original coordinate system is N * M * N^-1.
Obviously you can do this for any rotation about any axis, though this example will be particularly easy to compute. In fact here I'd be tempted not to bother dividing by at all and to use the adjugate matrix instead of the inverse, and then finally just halve all the entries in the final matrix, but you have to be careful with tricks like that as the two basis vectors in the plane being rotated must be the same length.
Once you've computed this matrix you can then apply my method from the last post to compute the x,y, and z rotations required. It might be worth checking what I have done is doing things in the order you want by computing the product

That's obviously the most general matrix you can get by multiplying three rotations about the axes. Check that my code usually recovers xi,eta and zeta. (There will certainly be times when it won't).
Conjugation isn't the only way to compute 3D rotation matrices. Some people use quaternions, but I like conjugation as it is easy to work out how to do it from scratch when you have forgotten what order the matrices should go in, and it works for other kinds of transformation as well.

Originally Posted by alunw

Here is an example of what I mean. Your original vector is (1,1,0). Lets call that y, and your axis of rotation is orthogonal to this. Lets call that z. To make things easy let's assume it is (0,0,1). Now y cross z which I make to be (1,-1,0), but I suppose could be (-1,1,0) depending on your sign convention for the cross product is a third vector orthogonal to both the other two which we can call x. (In my programs I use an LHS coordinate system, and picked the sign of the cross product so that (1,0,0) cross (0,1,0) is (0,0,1). So I think what I'm doing would work just as well in an ordinary RHS system)
So if we normalise the vectors to length 1 by dividing the first and third by sqrt(2) we have a new orthonormal coordinate system,(A rotation by theta about z looks like , so it is very easy to express any rotation. Let us call this matrix M
Now the matrix that tells us how to get from our new coordinates to our original coordinates is simply [x,y,z] i.e. . Call that matrix N
So the rotation matrix for a rotation of theta in our original coordinate system is N * M * N^-1.
Obviously you can do this for any rotation about any axis, though this example will be particularly easy to compute. In fact here I'd be tempted not to bother dividing by at all and to use the adjugate matrix instead of the inverse, and then finally just halve all the entries in the final matrix, but you have to be careful with tricks like that as the two basis vectors in the plane being rotated must be the same length.
Once you've computed this matrix you can then apply my method from the last post to compute the x,y, and z rotations required. It might be worth checking what I have done is doing things in the order you want by computing the product

That's obviously the most general matrix you can get by multiplying three rotations about the axes. Check that my code usually recovers xi,eta and zeta. (There will certainly be times when it won't).
Conjugation isn't the only way to compute 3D rotation matrices. Some people use quaternions, but I like conjugation as it is easy to work out how to do it from scratch when you have forgotten what order the matrices should go in, and it works for other kinds of transformation as well.

Thanks, I pretty much get what's going on. I use a cross product to create new axes of rotation, and create axis-angle matrix, then I can extract the Euler angles. When you say it usually recovers the angles are there known cases where it doesn't? It seems to work pretty well in my tests. For all angles +- 180 degrees.

Thanks

Rick

Originally Posted by alunw

Here is an example of what I mean. Your original vector is (1,1,0). Lets call that y, and your axis of rotation is orthogonal to this. Lets call that z. To make things easy let's assume it is (0,0,1). Now y cross z which I make to be (1,-1,0), but I suppose could be (-1,1,0) depending on your sign convention for the cross product is a third vector orthogonal to both the other two which we can call x. (In my programs I use an LHS coordinate system, and picked the sign of the cross product so that (1,0,0) cross (0,1,0) is (0,0,1). So I think what I'm doing would work just as well in an ordinary RHS system)
So if we normalise the vectors to length 1 by dividing the first and third by sqrt(2) we have a new orthonormal coordinate system,(A rotation by theta about z looks like , so it is very easy to express any rotation. Let us call this matrix M
Now the matrix that tells us how to get from our new coordinates to our original coordinates is simply [x,y,z] i.e. . Call that matrix N
So the rotation matrix for a rotation of theta in our original coordinate system is N * M * N^-1.
Obviously you can do this for any rotation about any axis, though this example will be particularly easy to compute. In fact here I'd be tempted not to bother dividing by at all and to use the adjugate matrix instead of the inverse, and then finally just halve all the entries in the final matrix, but you have to be careful with tricks like that as the two basis vectors in the plane being rotated must be the same length.
Once you've computed this matrix you can then apply my method from the last post to compute the x,y, and z rotations required. It might be worth checking what I have done is doing things in the order you want by computing the product

That's obviously the most general matrix you can get by multiplying three rotations about the axes. Check that my code usually recovers xi,eta and zeta. (There will certainly be times when it won't).
Conjugation isn't the only way to compute 3D rotation matrices. Some people use quaternions, but I like conjugation as it is easy to work out how to do it from scratch when you have forgotten what order the matrices should go in, and it works for other kinds of transformation as well.

my matrices are
X
1,0,0
0,cos,-sin,
0,sin,cos

Y
cos,0,sin
0,1,0
-sin,0,cos

Z
cos,-sin,0
sin,cos,0
0,0,1

when I multiply I get in row 1 col 3 sinY.

this should correspond to your inverse arcsin(transform[6]).
Yet in my code I put arcsin(3,1) to get y and it works. That's puzzling me at the moment. I would think I would put arcsin(1,3).....

transform[6] is the first entry on the last row (which I would most likely call a[2][0]). I think you are maybe thinking in Fortran and I am thinking in C. So I think we actually agree. But in any case it would not surprise me if for some reason my code is for the transpose of what you would get. I think once you have multiplied out the matrices you will be able to deduce a suitable formula for your set up.
I also used the opposite way round for the sin signs in my y matrix from you. Possibly I didn't think about that carefully enough in my code, or it arose naturally from the way I am viewing things.

Originally Posted by alunw

transform[6] is the first entry on the last row (which I would most likely call a[2][0]). I think you are maybe thinking in Fortran and I am thinking in C. So I think we actually agree. But in any case it would not surprise me if for some reason my code is for the transpose of what you would get. I think once you have multiplied out the matrices you will be able to deduce a suitable formula for your set up.
I also used the opposite way round for the sin signs in my y matrix from you. Possibly I didn't think about that carefully enough in my code, or it arose naturally from the way I am viewing things.

I am using a axis - angle matrix I gleaned from I think Euclideanspace.com or perhaps a textbook:

CAxisAngleMatrix::CAxisAngleRotationMatrix(float theta,float X, float Y, float Z)
{

CVector _v(X,Y,Z);

CVector v = _v.createUnitVector();

float Xs = v.getX() * v.getX();
float Ys = v.getY() * v.getY();
float Zs = v.getZ() * v.getZ();
float t = 1.0 -std::cos(theta);
float s = std::sin(theta);
float c = std::cos(theta);

e11 = (t * Xs) + std::cos(theta);

e12 = t*v.getX()*v.getY() + (s * v.getZ());
e13 = t*v.getX()*v.getZ() - (s * v.getY());

e21 = t*v.getX()*v.getY() - (s * v.getZ());
e22 = t*Ys + c;
e23 = t*v.getY()*v.getZ() + (s * v.getX());d

e31 = t*v.getX()*v.getZ() + (s * v.getY());
e32 = t*v.getY()*v.getZ() - (s * v.getX());
e33 = t*Zs + c;

}

In this code an axis is passed in in the form of X Y Z, so for instance
1,0,0 for around the X axis.

So, I'm not sure what grabbing the Euler angles out of this matrix entails.
The indices that you presented ([6] I translate to e31 and just seems to work but I'd like to understand why.

in my code, e31 (element row 3 col 1) in the axis angle matrix is

e31 = (1 - cos(theta) * v.getX() * v.getZ() + (sin(theta) * v.getY());


RIck

Originally Posted by rfrank

I am using a axis - angle matrix I gleaned from I think Euclideanspace.com or perhaps a textbook:

CAxisAngleMatrix::CAxisAngleRotationMatrix(float theta,float X, float Y, float Z)
{

CVector _v(X,Y,Z);

CVector v = _v.createUnitVector();

float Xs = v.getX() * v.getX();
float Ys = v.getY() * v.getY();
float Zs = v.getZ() * v.getZ();
float t = 1.0 -std::cos(theta);
float s = std::sin(theta);
float c = std::cos(theta);

e11 = (t * Xs) + std::cos(theta);

e12 = t*v.getX()*v.getY() + (s * v.getZ());
e13 = t*v.getX()*v.getZ() - (s * v.getY());

e21 = t*v.getX()*v.getY() - (s * v.getZ());
e22 = t*Ys + c;
e23 = t*v.getY()*v.getZ() + (s * v.getX());d

e31 = t*v.getX()*v.getZ() + (s * v.getY());
e32 = t*v.getY()*v.getZ() - (s * v.getX());
e33 = t*Zs + c;

}

In this code an axis is passed in in the form of X Y Z, so for instance
1,0,0 for around the X axis.

So, I'm not sure what grabbing the Euler angles out of this matrix entails.
The indices that you presented ([6] becomes for me 1,3 (or 0,2) just seem to work but I'd like to understand why.

in my code, e13 (element row 1 col 3) is

(1-cos(theta) X Z) - sin(theta)Y

RIck

I should add further:


I set up three of these:
CAxisAngleRotation rmX(theta,xAxis.X,xAxis.Y,xAxis.Z)
CAxisAngleRotation rmY(theta,yAxis.X,yAxis.Y,yAxis.Z)
CAxisAngleRotation rmZ(theta,zAxis.X,zAxis.Y,zAxis.Z)

CAxisAngleRotation temp = rmX * rmY * rmZ;

float eX,eY,eZ;
temp.getEuler(eX,eY,eZ);

CXRotation eXm(eX);
CYRotation eYm(eY);
CZRotation eZm(eZ);

CRotation tempR = eXm * eYm * eZm;

etc.


In this code, I set up the Axis angle matrices, extract the Euler, and
the create "standard" Euler matrices based on
CXRotation::CXRotation(float radians)
{
e11 = 1.0;
e12 = 0.0;
e13 = 0.0;

e21 = 0.0;
e22 = cos(radians);
e23 = -sin(radians);

e31 = 0.0;
e32 = sin(radians);
e33 = cos(radians);

}

CYRotation::CYRotation(float radians)
{
e11 = cos(radians);
e12 = 0.0;
e13 = sin(radians);

e21 = 0.0;
e22 = 1.0;
e23 = 0;

e31 = -sin(radians);
e32 = 0;
e33 = cos(radians);

}

CZRotation::CZRotation(float radians)
{
e11 = cos(radians);
e12 = -sin(radians);
e13 = 0.0;

e21 = sin(radians);
e22 = cos(radians);
e23 = 0;

e31 = 0.0;
e32 = 0.0;
e33 = 1.0;

}


What I notice is that the end product of the axis angle rotations looks somewhat different then the end product of the standard rotation matrices,
which leads me to wonder why. I guess I don't quite understand the axis-angle matrix derivation and how it might relate to the standard matrix derivation.


I also notice that if I transform a point with the product of the 3 axis angle matrices it doesn't fully seem to work consistently. The Z axis ends up moving in some cases.

However, things seem to work with the extracted Euler angles.

Perhaps there's a handedness problem?


Thanks

R.

My example code is based on doing first rotation about x axis, then y axis, then z axis where these are fixed (ie not intrinsic to the object being rotated), but your code looks to be creating a matrix by multiplying x*y*z rotations. Either you are using row vectors and multiplying on the right, or you are doing things in the opposite order from me because for my x*y means do y first then x, or you are using intrinsic rotations. In any case my code won't apply directly (though in the first and last cases you probably just need to transpose my matrix). Looking again at the piece of real code from which I extracted my sample code for decomposing a general rotation matrix back into three rotations I seem to have used exactly the opposite sign convention for x,y, and z to you (so that in the product of three rotations I wrote out in an earlier post the y matrix represents what I really did, but the x and z rotations are the other way round). So my code is right for an LHS coordinate system with anti-clockwise rotations positive or an RHS system with clockwise rotations positive. For the standard coordinate system (RHS, anticlockwise rotations positive) I think all the signs of my angles are probably wrong which is probably why you get so many +- 180s.

I really think the best thing to do is to form the product that you intend [X,Y,Z] to represent and then see how to recover X,Y,Z from it. by using trig identities. That's basically what I did. As for Euler angles, I don't know about them: I wrote my code because I was writing software that allowed a user to rotate an object about the three axes by clicking on various controls that did a small fixed rotation. They could click on the controls in any order any number of times, but I wanted to be able to recreate the final rotation matrix later from just three rotations. It sounded as though your original problem was much the same. From a brief look at the Wikipedia page on Euler angles they are not the same decomposition of a general rotation matrix as I have done.

Regarding what I said about the decomposition not being unique I was probably over-stating things somewhat, but there are certainly some matrices for which the decomposition is not unique, as you will surely find if you try to work things out in the way I am suggesting.