Maths - Angle between vectors - minorlogic

By: nobody ( Nobody/Anonymous )
quternion from 2 vec  
2003-06-19 13:43

inline void RotationArc( quaternion& q, const vector3& from, const vector3& to )
{
vector3 c( from*to );
float d = from%to;
q.set( c.x, c.y, c.z, d + (float)sqrt( from.len_squared()*to.len_squared() ) );
q.normalize();
}
this create a quternion thet rotate by shortest arc from to (vectors ). If somebody have faster solution let me know.

By: martinbaker ( Martin Baker )
RE: quternion from 2 vec  
2003-06-19 17:02

Thank very much for this, I'll link it from the following web page if that's alright.

I assume that vector3 overloads * to be cross product and % is dot product. Do you have a derivation for this algorithm?

I think the equations from the code are:
x = cross(from,to).x
y = cross(from,to).y
z = cross(from,to).z
w = dot(from,to) + sqrt((from.x*from.x + from.y*from.y + from.z*from.z)* (to.x*to.x + to.y*to.y + to.z*to.z))

On this page:
http://www.euclideanspace.com/maths/algebra/vectors/angleBetween/index.htm
There is code which uses the following:
w = dot(from,to)
x = cross(from,to).x* sqrt(1.0 - w*w)
y = cross(from,to).y* sqrt(1.0 - w*w)
z = cross(from,to).z* sqrt(1.0 - w*w)

This normalises from and to before the calculation, wheras your version only normalises the result so I guess your version is more efficient. It would be nice to prove these two are equivalent.

Martin

By: nobody ( Nobody/Anonymous )
RE: quternion from 2 vec  
2003-06-19 17:44

yes , * - is cross, % - dot
lets assume that we have vectors from and to and thay are both unit.

inline void RotationArcUnit( quaternion& q, const vector3& uf, const vector3& ut )
{
vector3 c ( uf*ut );
q.set( c.x, c.y, c.z, (uf%ut) + 1.0f );
q.normalize(); // norming
}

your code is wrong , method that you try to use look:
vector3 c ( v0*v1 );
float d = v0%v1;
float w = (float)sqrt( 2.0f + d + d );
float s = 1.0f/w;
quaternion( c.x * s, c.y * s, c.z * s, 0.5f*w );
from this : cos2a = 1 - 2(sina)^2 = 2(cosa)^2 - 1;
/* but it give not stable result near 180 deg*/
//---
if we scale quat by "w" we get:
quaternion( c.x , c.y , c.z , 0.5f*w*w );
quaternion( c.x , c.y , c.z , d + 1 );
//-------------
ok , we know that cross and dot product scaled by
Len1*Len2 and scale one by this if vectors is not unit:

Len1*Len2 = sqrt(Len1*Len1)*sqrt(Len2 *Len2 ) = sqrt(Len1*Len1*Len2 *Len2 )

And at last we can catch the case of 180 deg or some vectors length == 0

inline void RotationArc( quaternion& q, const vector3& from, const vector3& to )
{
vector3 c( from*to );
float d = from%to;
q.set( c.x, c.y, c.z, d + (float)sqrt( from.len_squared()*to.len_squared() ) );
/* here check the norm of quat , if it near zero catch a special case */

q.normalize();
}
this create a quternion thet rotate by shortest arc from to (vectors ). If somebody have faster solution let me know.

By: nobody ( Nobody/Anonymous )
RE: quternion from 2 vec  
2003-06-20 08:12

Please let me know if all clear.

minorlogic

By: martinbaker ( Martin Baker )
RE: quternion from 2 vec  
2003-06-20 08:50

Minorlogic,

Well since you ask, I am having trouble working out how you got:
float w = (float)sqrt( 2.0f + d + d );

Here are my unsuccessful attempts so far to derive:

Let x,y,z,w be elements of quaternion, these can be expressed in terms of axis angle:
s = sin(angle/2)
x = v.x *s
y = v.y *s
z = v.z *s
w = cos(angle/2)

Using half angle formula from the following page gives:
http://www.euclideanspace.com/maths/geometry/trig/index.htm
sin(angle/2) = sqrt(0.5*(1 - cos (angle)))
cos(angle/2) = sqrt(0.5*(1 + cos (angle)))

so substituting in quaternion formula gives:
s = sqrt(0.5*(1 - cos (angle)))
x = v.x *s
y = v.y *s
z = v.z *s
w = sqrt(0.5*(1 + cos (angle)))

However we know that:
angle = arcos(v1.v2)
axis = v = (v1 x v2 )/(|v1| |v2|)

rearranging and assuming that v1 and v2 are normalised gives:
cos(angle) = v1.v2
v = v1 x v2

so substituting this in quaternion formula gives:
s = sqrt(0.5*(1 - v1.v2))
x = (v1 x v2).x *s
y = (v1 x v2).y *s
z = (v1 x v2).z *s
w = sqrt(0.5*(1 + v1.v2))

I realise that w in your equation may not be the same as this as w here (because quaternion may not be normalised and therefore multiplied by some factor) but even so I cant derive your equation.

Martin

By: nobody ( Nobody/Anonymous )
RE: quternion from 2 vec  
2003-06-20 10:51

Let's play with the second part of function , we assume that work with unit quats.

inline quaternion RotationArc(vector3& v0, vector3& v1)
{
v0.norm();
v1.norm();

vector3 c = cross( v0, v1 );
float d = dot( v0, v1);

// cos (2a) = 2(cosa)^2 -1 //basic trigonometry
// sin (2a) = 2sina cosa //basic trigonometry
// cos (a/2) = sqrt( 0.5*(cosa + 1) );
// sin (a/2) = 0.5*sin (a)/cos(a/2);
/*
cos (a/2) = sqrt( 0.5*(cosa + 1) ) = 2*sqrt( 0.5*(cosa + 1) ) /2 = sqrt( 4*0.5*(cosa + 1) )/2 = sqrt( 2*(cosa + 1) )/2; just simple math
*/
float cos_half_a = 0.5f*(float)sqrt( 2.0f + d + d );
float s = 0.5/cos_half_a; //sin (a/2) = 0.5*sin (a)/cos(a/2);
/* here we can return a such result
return quaternion( c.x * s, c.y * s, c.z * s, cos_half_a );
but our quat in near 180 degrees will not be stable ( not unit ) .. just experiment with it to see.
we can can scale qaternion by some value and than normalize it
let's scale by "1/s" "cos_half_a/0.5"
quaternion( c.x , c.y , c.z , cos_half_a^2/0.5 );
(cos_half_a^2)/0.5 = 0.5f*0.5f*( 2.0f + d + d )/0.5 = 1 + d;
quaternion( c.x , c.y , c.z , 1.0f + d);
normalize it
*/

quaternion q( c.x , c.y , c.z , 1.0f + d);
q.normalize();
return q;
}
//-----------------------------
if the input vectors is not unit the quat is scaled by their length sum. To get unit quat we can normalize it.

vector3 c = cross( v0, v1 );
float d = dot( v0, v1);
quaternion q( c.x * s, c.y * s, c.z * s, d);
q.normalize();
q.w += 1;
q.normalize();

but little faster and may be stable is scale 1 by sum of the length
q.set( c.x, c.y, c.z, d + (float)sqrt( from.len_squared()*to.len_squared() ) );

I hope this will be helpfull.
P.S. It is my own implementation, but i think the proof of math is clear.
minorlogic

By: martinbaker ( Martin Baker )
RE: quternion from 2 vec  
2003-06-20 16:59

Minorlogic,
Thanks for your patience, I think I'm nearly with you except that I keep coming up with an extra 'sin(angle)' term as follows. Can you spot my error?

s = sin(angle/2) = 0.5 sin(angle) / cos(angle/2)
x = v.x *s
y = v.y *s
z = v.z *s
w = cos(angle/2)

multiply x,y,z and w by 1/s

x = v.x
y = v.y
z = v.z
w = 2 * cos(angle/2) * cos(angle/2) / sin(angle)

now substitute cos(angle/2) = sqrt(0.5*(1 + cos (angle)))

x = v.x
y = v.y
z = v.z
w = (1 + cos (angle)) / sin(angle)

x = cross( v0, v1).x
y = cross( v0, v1).y
z = cross( v0, v1).z
w = (1 + dot( v0, v1)/ sin(angle)

By: nobody ( Nobody/Anonymous )
RE: quternion from 2 vec  
2003-06-23 07:14

Ok.

"Thanks for your patience, I think I'm nearly with you except that I keep coming up with an extra 'sin(angle)' term as follows. Can you spot my error"

Do you remember that cross product produce perpendicular vector witch length is sum of source vectors length and multiplied by SIN of angle between ?
So i think you forget that cross product yet multiplied by SIN( a ). it is from there.

cross = RotationAxe * Len1*Len2 *sin( a )
if vectors is unit it look
cross = RotationAxe * sin( a )
quaternion is a RotationAxe * sin( a/2 )
// sin (a/2) = 0.5*sin (a)/cos(a/2);
so if we know the cos(a/2) we get

quaternion.x = RotationAxe.x * sin( a/2 ) =
RotationAxe.x * 0.5*sin (a)/cos(a/2) = (RotationAxe.x * sin (a))*(0.5/cos(a/2)) = cross.x*(0.5/cos(a/2));

Thats why in my previous posts
float s = 0.5/cos(a/2);

So if we know the cos(a/2) we can calculate the vector part of quternion too.

So just rewrite your code:

//s = sin(angle/2) = 0.5 sin(angle) / cos(angle/2) -
s = sin(angle/2)/sin(angle) = 0.5 / cos(angle/2)
x = cross.x *s //x = cross.x * 0.5 / cos(angle/2) = (v.x*sin(a)) *0.5 / cos(angle/2);
y = cross.y *s
z = cross.z *s
w = cos(angle/2)

multiply x,y,z and w by 1/s

x = cross.x
y = cross.y
z = cross.z
//w = 2 * cos(angle/2) * cos(angle/2) / sin(angle)
w = 2 * cos(angle/2) * cos(angle/2)

now substitute cos(angle/2) = sqrt(0.5*(1 + cos (angle)))

x = cross.x
y = cross.y
z = cross.z
w = (1 + cos (angle))

x = cross( v0, v1).x
y = cross( v0, v1).y
z = cross( v0, v1).z
w = (1 + dot( v0, v1)

Thats all, if you agree , than let me know :).

minorlogic

By: martinbaker ( Martin Baker )
RE: quternion from 2 vec  
2003-06-23 15:57

Minorlogic,

> Thats all, if you agree , than let me know :).

Yes, I agree, sorry I was so slow to catch on and thank you for your patience.

I have updated the web page with your method, I hope that is alright?

Thanks,

Martin

By: nobody ( Nobody/Anonymous )
RE: quternion from 2 vec  
2003-06-23 17:36

Thanks for attention. Litle bit of your code, it still use 3 sqrt , but can only 2. Possible it is done to clear code, or some kind of stabillity of math ?
You realy don't need to norm input vectors , and it is a main trip of code. The spetial case of 180 degrees can be taken before quternion normalization, just check w is close to zero. It can be if dot between close to 1.0 or length of some input vector is close to 0. Any of this cases can be catched before normalizing.
inline void RotationArc( quaternion& q, const vector3& from, const vector3& to )
{
vector3 c = cross( from,to );
float sum_of_len = (float)sqrt( from.len_squared()*to.len_squared() ) ;
float d = dot(from, to);
q.set( c.x, c.y, c.z, d + sum_of_len);
if(q.w < 0.000001f )
{
// check the sum_of_len near zero and if it is not than it is close to 180 degrees
}
q.normalize();
}




//-----
if you want to be clear about math than better use an other nice method that can be shown on image.

inline void RotationArc( quaternion& q,
const vector3& From,
const vector3& To){

vector3 UnitFrom( From );
From.norm();
vector3 UnitTo( To );
To.norm();

vector3 bisector( UnitFrom + UnitTo );
bisector.norm();

float cos_half_a = dot( UnitFrom, bisector );
vector3 BCross = cross( UnitFrom, bisector );

q.set( BCross.x, BCross.y, BCross.z, cos_half_a );
}

it is simpler to understand , the bisecor vector lay between source vectors and have a half angle with "from".

and strange for me how you take a special cases , if angle = 0 than method work and produce correct quaternion, if angle is close to 180 degrees, than i recomend to use cross product with some axe like this:

float CosA = UnitFrom%UnitTo;
if( CosA < - 0.99999f ){ // angle close to PI
vector3 Cross( 0, UnitFrom.x, -UnitFrom.y ); // cross( (1, 0, 0),UnitFrom )
if( ( UnitFrom.z*UnitFrom.z ) > ( UnitFrom.y*UnitFrom.y ) ){
// if (0, 1, 0) Cross len> (1, 0, 0) Cross len
Cross.set( -UnitFrom.z, 0, UnitFrom.x ); // cross( (0, 1, 0),UnitFrom )
}
Cross.norm();
q.set( Cross.x, Cross.y, Cross.z, 0.0f );
}else .....

//----
I hope you are not boring yet .... :)

By: martinbaker ( Martin Baker )
RE: quternion from 2 vec  
2003-06-24 10:29

Minorlogic,

Thanks again, I am updating the webpage as you suggest. I tried expanding out the vector operations (see below). I was just checking to see if I could find any more optimisations but I could not find any.

I like the method which shows the half angle graphically, I will try to add some diagrams to the website to show this. This reminds me of the proof of the use of quaternions for rotations which involves reflection around the half angle.

> I hope you are not boring yet
No, It is a useful topic and I appreciate your help with improving the website.

Thanks Martin
------------------------------------------

float xfrom=from.x;
float yfrom=from.y;
float zfrom=from.z;
float xto=to.x;
float yto=to.y;
float zto=to.z;
float d = xfrom * xto + yfrom * yto + zfrom * zto; //dot(from, to);
float xq = yfrom * zto - yto * zfrom ; // vector3 c = cross( from,to );
float yq = zfrom * xto - zto * xfrom;
float zq = xfrom * yto - xto * yfrom;
float wq = (float)sqrt((xfrom*xfrom + yfrom*yfrom + zfrom*zfrom)*(xto*xto + yto*yto + zto*zto)) +d;
if (wq < 0.0001) {
xq =-zfrom; // Cross.set( -UnitFrom.z, 0, UnitFrom.x );
yq =yfrom;
zq =xfrom;
wq =0
}
float cnorm = (float)sqrt((xq*xq + yq*yq + zq*zq + wq*wq) ;
xq /= cnorm; // q.normalize();
yq /= cnorm;
zq /= cnorm;
wq /= cnorm;

example where vectors are 180 apart:
float xfrom=1
float yfrom=0
float zfrom=0
float xto=-1
float yto=0
float zto=0
float d = -1 //dot(from, to);
float xq = 0 ; // vector3 c = cross( from,to );
float yq = 0;
float zq = 0;
float wq = 1-1 = 0;
if (wq < 0.0001) {
xq =0; // Cross.set( -UnitFrom.z, 0, UnitFrom.x );
yq =0;
zq =1;
wq =0
}
float cnorm = (float)sqrt(0) ;
xq /= cnorm; // q.normalize();
yq /= cnorm;
zq /= cnorm;
wq /= cnorm;