vX = Normalize(&(*vInX - vPos));
vY = Normalize(&(*vInY - vPos));
vZ = Normalize(&Cross(&vX, &vY));Finding angles of rotation between 2x2 vectors
Author
I suck in math. Seriously. I've spent days trying to figure this out, and even after reading a dozen posts on similar subjects, I still can't seem to get this to work. Problem is the following: I have a pair of orientation vectors in 3D - the "Forward" vector (Y), and the "Right" vector (X). Using these I can calculate the Z vector using a cross product.
Well, you have AX=B where A is an unknown rotation matrix, X is a matrix whose columns are the points before rotation and B is the points in the same order after rotation. So A=BX^-1. Then the eigenvector for that matrix corresponding to one is the axis of rotation. The formula for constructing an axis angle matrix should give you the means to recover the angle given the axis.
And note that there's infinite set of unit-length vectors you can use in glRotate to turn point A to point B .
1: All such vectors lie in plane with normal (A-B) .
(well, proving is left as exercise for the reader :-)
It's not hard to find angle for some vector from that set. For instance,{ vector=(A+B)*0.5 , angle=180 degrees } will turn A to B . As well will do { vector=(A Cross B) , angle=angle_between_a_and_b }. And so-on.
You need at least 2 pairs of points to find rotation,let's name 'em
P,P',Q,Q'
- P rotated to P' and Q rotated to Q'
you just need to find good vector V that can rotate both. It's placed on intersection of set of vectors that can rotate P to P' with set of vectors that can rotate Q to Q'. That is, it's intersection of plane with normal P-P' and other plane with normal Q-Q' .(read 1: and think!)
That is, it's just perpendicular to both plane normals, and it's
V=(P-P')Cross(Q-Q')
i'm bored now and it's night and i worked all day, hope you can find angle yourself. ....but no,i will solve to the end anyway....it's just angle between ( P - V*(P Dot V)/(V Dot V) ) and same formule for P' .
It's just P and P' both projected to plane with normal V, so, after turn by that angle, their projections will match. Simple.
Author
Thanks for your replies. LilBudyWizer, you lost me on the word "eigenvector". :)
Dmytry, I got the first formula you posted to work:
Using the above I am able to find the rotation axis... however your second formula doesn't make much sense? "(V Dot V)" part in particular... Could you post the entire thing? Thanks in advance!
Dmytry, I got the first formula you posted to work:
vAxis = Normalize(&Cross(&(pLink->vX - vX), &(pLink->vZ - vZ)));Using the above I am able to find the rotation axis... however your second formula doesn't make much sense? "(V Dot V)" part in particular... Could you post the entire thing? Thanks in advance!
I'm writing (V Dot V) because it's sometimes written as
V2
V2
V.length2
|V|
||V||2
It's just squared length, equal to Vx2+Vy2+Vz2
and (V Dot V) is the only notation that makes some sence by itself.
and i meant angle between vectors
( P - V*(P Dot V)/(V Dot V) )
,
( P' - V*(P' Dot V)/(V Dot V) )
You can find angle using
AngleBetweenVectors(A,B)=atan2( A Cross B , A Dot B )
If your V is normalized, then you don't need /(V Dot V) (it's there instead of normalizing)
(i don't sure angle have right sign for glRotate)
Note that it's not guaranteed to turn any Q to Q' because not every two pairs corresponds to valid rotation.
Also your notation
&(pLink->vX - vX), &(pLink->vZ - vZ)
my V is not related to your v* things
V2
V2
V.length2
|V|
||V||2
It's just squared length, equal to Vx2+Vy2+Vz2
and (V Dot V) is the only notation that makes some sence by itself.
and i meant angle between vectors
( P - V*(P Dot V)/(V Dot V) )
,
( P' - V*(P' Dot V)/(V Dot V) )
You can find angle using
AngleBetweenVectors(A,B)=atan2( A Cross B , A Dot B )
If your V is normalized, then you don't need /(V Dot V) (it's there instead of normalizing)
(i don't sure angle have right sign for glRotate)
Note that it's not guaranteed to turn any Q to Q' because not every two pairs corresponds to valid rotation.
Also your notation
&(pLink->vX - vX), &(pLink->vZ - vZ)
my V is not related to your v* things
An eigen vector is a vector x such that Ax=dx for some scalar d. Finding those isn't a simple task so it wouldn't have been a good way to find the axis of rotation. It would be easier to extract the axis and angle from the matrix, but dmytry's method would be even easier. The only problem is if P'-P=0, Q'-Q=0 or P'-P=d(Q'-Q). If any of those are true then your cross product is going to be the zero vector.
That's pretty straightforward if you have a math lib available. For instance someone posted a small .cpp file in this math forum a bit earlier. Parse the recent threads.
You have two frames (<=> 3x3 matrices). M1 : X1,Y1,(Z1) and M2 : X2,Y2,(Z2).
(
@Dmytry.
I don't think he wants a rotation between two vectors. It's between two frames ! Thus there is only one matricial solution, and all Euler angles are related to some modulo indetermlination.
)
I suppose that X1.Y1=0 and X2.Y2=0, (perp vectors) and that you have orthonormal or at least orthogonal frames.
All you need is to 'divide' M2 by M1. It can also be expressed with quaternions. In the end all you need is Euler angles.
EA3 = Euler(M2 * tM1)
(
Note that the inverse of an orthonormal matrix is simply its transposition. My math lib contains a mul by transposed, which is very efficiently implemented in SSE and 3DNow since you simply have to multiply and accumulate the. In any case this makes transposition implicit and thus gives it for free
)
Or
EA3 = Euler(Quat(M2), Quat(M1)) done.
In 2D it's very simply expressed by Complex numbers with their dual representation, r*e^it or a+i*b. For rotations, quaternions are the complex numbers of 3D (while they loose commutativity).
As you can see it would have taken me 10 seconds to write this line of code. It's quite automatic once linear algebra becomes as obvious as elementary calculus. I can't advise you anything else but make the effort to spend a few 'days' learning it. In the end it will save you a lot of time if you plan to write some decent 3D software.
[Edited by - Charles B on September 12, 2004 7:39:30 AM]
You have two frames (<=> 3x3 matrices). M1 : X1,Y1,(Z1) and M2 : X2,Y2,(Z2).
(
@Dmytry.
I don't think he wants a rotation between two vectors. It's between two frames ! Thus there is only one matricial solution, and all Euler angles are related to some modulo indetermlination.
)
I suppose that X1.Y1=0 and X2.Y2=0, (perp vectors) and that you have orthonormal or at least orthogonal frames.
All you need is to 'divide' M2 by M1. It can also be expressed with quaternions. In the end all you need is Euler angles.
EA3 = Euler(M2 * tM1)
(
Note that the inverse of an orthonormal matrix is simply its transposition. My math lib contains a mul by transposed, which is very efficiently implemented in SSE and 3DNow since you simply have to multiply and accumulate the. In any case this makes transposition implicit and thus gives it for free
)
Or
EA3 = Euler(Quat(M2), Quat(M1)) done.
In 2D it's very simply expressed by Complex numbers with their dual representation, r*e^it or a+i*b. For rotations, quaternions are the complex numbers of 3D (while they loose commutativity).
Quote:
I suck in math. Seriously. I've spent days trying to figure this out ...
As you can see it would have taken me 10 seconds to write this line of code. It's quite automatic once linear algebra becomes as obvious as elementary calculus. I can't advise you anything else but make the effort to spend a few 'days' learning it. In the end it will save you a lot of time if you plan to write some decent 3D software.
[Edited by - Charles B on September 12, 2004 7:39:30 AM]
hmmmmm....i interpreted what he said
as that he have two sets of points, and he want to find axis/angle that rotated first set to the second. Problem also was i'll-stated because there's infinitely many vectors x,y,z that can rotate A to B (even if we consider only unit-length vectors).
Anyway , doesn't matter, my solution will work if you plug two columns of one matrix as P,Q and same columns of other matrix as P', Q'
(assuming matrices are orthonormal)
Of course it can be done with quaternions, but as he said he just simply want argument to glRotate
(!!!!! and IIRC glRotate(a,x,y,z) rotates around given axis x,y,z by a degrees, it's not related to euler angles !!!!!!)
in fact it's very simple to construct quaternion from axis-and-angle, if he want quaternion instead of arguments for glRotate.
BTW:
There's some problem with my solution, because 2 points is not always enough and ,and also division-by-zero is sometimes possible.
First, for obvious reasons, P and Q must be choosen that P is not parallel to Q, that is ,there must be no scalar a that P=a*Q .
If P-P' = 0 (*): In that cases you just have to rotate around P axis, that is, V=P.
Q-Q' = 0 : same as above,except that V=Q
If both P-P' = 0 and Q-Q'=0, you don't have to rotate at all.
If (P-P')Cross(Q-Q')= 0 , you just have to choose other pair of P,Q,P',Q' .... or i think you can just set V=P-Q .
Should be robust enough now.
*************************************************
(*):better use "if |(P-P').x|<(1E-6)" , same for Q .
edit, and unnecessary final words:
and of course if i have understood what OP said... Instead of doing that crazy(but also correct) vectors stuff,
i also would rather write
M2=Mx*M1
and just find Mx=M2*inverse(M1)
as Charles B absolutely correctly wrote.(and of course, for orthonormal, inverse=transpose.)
Must be because of my participation in threads like "Do it as hard as possible"
(at least i still believe my solution works some precents faster :-)
[grin]
[Edited by - Dmytry on September 12, 2004 7:22:39 AM]
Quote:
Now I have a second set of the same vertices, in a different position. They have been rotated somehow, around an arbitrary axis that I don't know. What I need is to figure out the angles I need to use in glRotatef(angle, x, y, z); that would translate point A to point B.
as that he have two sets of points, and he want to find axis/angle that rotated first set to the second. Problem also was i'll-stated because there's infinitely many vectors x,y,z that can rotate A to B (even if we consider only unit-length vectors).
Anyway , doesn't matter, my solution will work if you plug two columns of one matrix as P,Q and same columns of other matrix as P', Q'
(assuming matrices are orthonormal)
Of course it can be done with quaternions, but as he said he just simply want argument to glRotate
(!!!!! and IIRC glRotate(a,x,y,z) rotates around given axis x,y,z by a degrees, it's not related to euler angles !!!!!!)
in fact it's very simple to construct quaternion from axis-and-angle, if he want quaternion instead of arguments for glRotate.
BTW:
There's some problem with my solution, because 2 points is not always enough and ,and also division-by-zero is sometimes possible.
First, for obvious reasons, P and Q must be choosen that P is not parallel to Q, that is ,there must be no scalar a that P=a*Q .
If P-P' = 0 (*): In that cases you just have to rotate around P axis, that is, V=P.
Q-Q' = 0 : same as above,except that V=Q
If both P-P' = 0 and Q-Q'=0, you don't have to rotate at all.
If (P-P')Cross(Q-Q')= 0 , you just have to choose other pair of P,Q,P',Q' .... or i think you can just set V=P-Q .
Should be robust enough now.
*************************************************
(*):better use "if |(P-P').x|<(1E-6)" , same for Q .
edit, and unnecessary final words:
and of course if i have understood what OP said... Instead of doing that crazy(but also correct) vectors stuff,
i also would rather write
M2=Mx*M1
and just find Mx=M2*inverse(M1)
as Charles B absolutely correctly wrote.(and of course, for orthonormal, inverse=transpose.)
Must be because of my participation in threads like "Do it as hard as possible"
(at least i still believe my solution works some precents faster :-)
[grin]
[Edited by - Dmytry on September 12, 2004 7:22:39 AM]
Quote:Original post by Dmytry
hmmmmm....i interpreted what he said
He has 3 points per set (Pos, InX, InY) <=> 2 vectors (X,Y) + one translation <=> 1 frame with the normal. I suppose these three points have something to do with 3 triangle vertices, he probably draws them like that in GL. BTW I did not speak of the translations but that's obvious (substraction).
Quote:
You need at least 2 pairs of points to find rotation,let's name 'em
P,P',Q,Q'
- P rotated to P' and Q rotated to Q'
I think your original post, 2 points per set (P,Q) <=> 1 vector and one translation is what I would call the trackball simulation algorithm except here it's affine space while the trackball issue is in linear space, just substract the translations.
Quote:
...because there's infinitely many vectors x,y,z that can rotate A to B
The quaternion that divises the two vectors, gives you the unique minimum rotation between the two vectors.
Quote:
(!!!!! and IIRC glRotate(a,x,y,z) rotates around given axis x,y,z by a degrees, it's not related to euler angles !!!!!!)
Point for you here ;) Was to quick here. Replace :
AR3 = AxisRotation( Quat(M2)/Quat(M1) )
As you suggested AxisRotation() 'cast' is indeed very fast. It mainly requires one atan2(), (and BTW it can be very fastly computed with Taylor, I can give you a very interesting link on it)
Quote:
There's some problem with my solution, ...
First, for obvious reasons, P and Q must be choosen that P is not parallel to Q.
Sure it's supposed to give a frame. The 'triangle' is not degenerate.
Quote:
and of course ... and just find Mx=M2*inverse(M1)
...
(at least i still believe my solution works some precents faster :-)
Converting to quaternions start would give a faster 'division' and the conversion to axis and rotation quasi immediate. Alas the two matrix to quaternions would probably cost a lot.
So the solution with Matrix to Axis Rotation directly is probably the fastest way to do. And if you develop your solution I am pretty certain it gives the same atomic equations and code. (*)
BTW When I'll release my math library, the function documentations will be given with benchmarks (auto generated by the project) with tables for each compiler and processor.
Quote:
Must be because of my participation in threads like "Do it as hard as possible"
(at least i still believe my solution works some precents faster :-)
Yes maybe, but you have not developped the formula enough (*), I'll probably be able to give you the exact answer one month soon :)
well,
in code:
(will work only with orthogonal P and Q)
i think should work,but haven't checked. Angle returned in radians. Sign of angle may be not what glRotate expects.
It _could be_ faster than matrix-to-quaternion stuff, but there's too many sqrts...and some branching(that actually can make it work faster for special cases).
and of course there's only one shortest-arc rotation ({ vector=(A Cross B) , angle=angle_between_a_and_b }),but we want not shortest arc rotation but rotation that can rotate both to their positions....
edit:some fixes....
in code:
(will work only with orthogonal P and Q)
const float smallnumber=1E-6;inline bool VerySmall(vec3 A){ return abs(A.x)+abs(A.y)+abs(A.z)<smallnumber;};inline bool NotVerySmall(vec3 A){ return abs(A.x)+abs(A.y)+abs(A.z)>smallnumber;};// don't have the name for it...DoThisThing(vec3 P,vec3 Q,vec3 P_,vec3 Q_,float &angle, vec3 &axis){ vec3 DP=P-P_, DQ=Q-Q_;vec3 V;if(NotVerySmall(DP)){// for branch prediction, then part should happen more frequently. if(NotVerySmall(DQ)){ V=CrossProd(DP,DQ);if(VerySmall(V)){if(DotProd(DP,DQ)>0){V=P-Q;}else{V=(P+Q);}} }else{//DQ is very small. V=DQ; }}else{//DP is very small if(NotVerySmall(DQ)){ V=DP; }else{// both is very small ,no rotation axis=vec3(1,0,0);angle=0;return(); }}Normalize(V);P-=V*DotProd(P,V);P_-=V*DotProd(P_,V);angle=atan2(Length(CrossProd(P,P_)),DotProd(P,P_));axis=V;}i think should work,but haven't checked. Angle returned in radians. Sign of angle may be not what glRotate expects.
It _could be_ faster than matrix-to-quaternion stuff, but there's too many sqrts...and some branching(that actually can make it work faster for special cases).
and of course there's only one shortest-arc rotation ({ vector=(A Cross B) , angle=angle_between_a_and_b }),but we want not shortest arc rotation but rotation that can rotate both to their positions....
edit:some fixes....
This topic is closed to new replies.
Advertisement
Popular Topics
Advertisement
