Transforming Angular Velocity

Started by
13 comments, last by JoeJ 6 years, 3 months ago

My question: is it possible to transform multiple angular velocities so that they can be reinserted as one? My research is below:


// This works
quat quaternion1 = GEQuaternionFromAngleRadians(angleRadiansVector1);
quat quaternion2 = GEMultiplyQuaternions(quaternion1, GEQuaternionFromAngleRadians(angleRadiansVector2));
quat quaternion3 = GEMultiplyQuaternions(quaternion2, GEQuaternionFromAngleRadians(angleRadiansVector3));

glMultMatrixf(GEMat4FromQuaternion(quaternion3).array);

// The first two work fine but not the third. Why?
quat quaternion1 = GEQuaternionFromAngleRadians(angleRadiansVector1);
vec3 vector1 = GETransformQuaternionAndVector(quaternion1, angularVelocity1);
	
quat quaternion2 = GEQuaternionFromAngleRadians(angleRadiansVector2);
vec3 vector2 = GETransformQuaternionAndVector(quaternion2, angularVelocity2);

// This doesn't work
//quat quaternion3 = GEQuaternionFromAngleRadians(angleRadiansVector3);
//vec3 vector3 = GETransformQuaternionAndVector(quaternion3, angularVelocity3);

vec3 angleVelocity = GEAddVectors(vector1, vector2);
// Does not work: vec3 angleVelocity = GEAddVectors(vector1, GEAddVectors(vector2, vector3));

static vec3 angleRadiansVector;
vec3 angularAcceleration = GESetVector(0.0, 0.0, 0.0);

// Sending it through one angular velocity later in my motion engine
angleVelocity = GEAddVectors(angleVelocity, GEMultiplyVectorAndScalar(angularAcceleration, timeStep));
angleRadiansVector = GEAddVectors(angleRadiansVector, GEMultiplyVectorAndScalar(angleVelocity, timeStep));

glMultMatrixf(GEMat4FromEulerAngle(angleRadiansVector).array);

Also how do I combine multiple angularAcceleration variables? Is there an easier way to transform the angular values?

Advertisement
39 minutes ago, Jonathan2006 said:

My question: is it possible to transform multiple angular velocities so that they can be reinserted as one?

Do you mean finding the average?

Could you maybe explain or sketch an idea of what your expecting?

That was one I didn't think of. Thanks for the averaging idea. I still cannot get it working. I am trying to mimic OpenGL's glRotate(); except I need to know the final angular velocity and acceleration vectors. The physics engine I am working on only accepts one angular velocity and angular acceleration vector. Below I am posting more of the code:


// The Angle Radians Test (code used both by the matrix and quaternions)
	angularVelocity1 = GEAddVectors(angularVelocity1, GEMultiplyVectorAndScalar(angularAcceleration1, timeStep));
	angleRadiansVector1 = GEAddVectors(angleRadiansVector1, GEMultiplyVectorAndScalar(angularVelocity1, timeStep));

	angularVelocity2 = GEAddVectors(angularVelocity2, GEMultiplyVectorAndScalar(angularAcceleration2, timeStep));
	angleRadiansVector2 = GEAddVectors(angleRadiansVector2, GEMultiplyVectorAndScalar(angularVelocity2, timeStep));
	
	angularVelocity3 = GEAddVectors(angularVelocity3, GEMultiplyVectorAndScalar(angularAcceleration3, timeStep));
	angleRadiansVector3 = GEAddVectors(angleRadiansVector3, GEMultiplyVectorAndScalar(angularVelocity3, timeStep));

// Matrix (This works and is used to test the accuracy of the quaternions)
	mat4 rotateMatrix = GEMultiplyMat4(GEMat4FromRotate(angleRadiansVector1.x, 1.0, 0.0, 0.0), GEIdentityMat4());
	rotateMatrix = GEMultiplyMat4(GEMat4FromRotate(angleRadiansVector2.z, 0.0, 0.0, 1.0), rotateMatrix);
	rotateMatrix = GEMultiplyMat4(GEMat4FromRotate(angleRadiansVector3.x, 1.0, 0.0, 0.0), rotateMatrix);

glMultMatrixf(rotateMatrix.array);
glutWireSphere(2.0, 7.0, 7.0);

The first few lines are extra code that I was hoping I could leave out in the first post. The rotateMatrix works just like OpenGL's glRotate(); and is the desired output. I am trying to transform the angular vectors the same way OpenGL does when multipling two matrices.

Either way you combine forces that cause object to rotate force1+force2+force

And apply that to final acceleration from torque or correct me if im wrong (bit tired)  you just add acceleration since its a vector 

You'll need to transport them too if they don't apply to the same point (see Varignon theorem).

Yes, usually you represent angular velocitiy as a 3D vector. This representation is like axis and angle rotation, but the unit axis is scaled by the angle, so vector length = angle, and it still works with negative angles. As long as those vectors are in the same space, you can add / subtract them like any vector. Otherwise transform them to mach space, add, transform to another space, whatever.

The same is true for torque.

Due to the axis / angle relationship, you can also convert angular velocity to a quaternion, e.g. to integrate the angular velocity to body orientation, but be aware that quaternions have limitatins about max angle they can represent - they can mot represent multiple spins, while angular velocity as a vector can.

It's not clear to me what you try to achieve, but if you want to convert multiple rotations to a single angular velocity,

you do all those rotations in order,

calculate the rotation form original orientation to resulting orientation,

convert this final rotation to axis / angle and finally to angular velocity (where your desired speed affects only the length of the vector)

 

Edit: Looking at your code again it seems you know the things i wrote already. I assume your problem comes from flipping cases, similar to the problem of averaging 3 or more angles, but i'm not sure (you can test by using different values for your 3 rotations, in some cases it might work as intended).

What is this: GEQuaternionFromAngleRadians(angleRadiansVector3); ?

Do you mean by 'angleRadiansVector' a axis / angle rotation where the vector length = angle?

 

Oh, wait - all wrong. I think your mistake is this: On top of your snippet you do 3 rotations in a given order. Later you try to add those rotations in a order independent way and expect the same result, but this is wrong if you think about it.

So the solution should be indeed to use the rotation between initial orientation and resulting orientation after 3 rotations has been applied, as explained above.

 

 

Thanks for the posts and ideas. I already know how to get the final values for position, linear velocity, linear acceleration, and Euler angle.

I have been researching quaternions some more. if it’s even possible I would like to know the instant angular velocity and angular acceleration which I tried to find in the code below. If it’s not possible to find the instant value I am confused as to what should be a vector and what should be a quaternion in your post above.

If it matters I am trying to use instant quaternions. I had problems with oscillating after a perfect collision using just the original and resulting vectors (I have not tried it with quaternions though).

The C code below is as close as I can get to the real motion engine:


void rotateShapeX(mat4 *m, quat *q, vec3 *angleRadians, vec3 *angularVelocity, vec3 *angularAcceleration)
{
	*angularVelocity = GEAddVectors(*angularVelocity, GEMultiplyVectorAndScalar(*angularAcceleration, globalTimeStep));
	*angleRadians = GEAddVectors(*angleRadians, GEMultiplyVectorAndScalar(*angularVelocity, globalTimeStep));
	
	*q = GEQuaternionFromAxisAndAngle(*angularVelocity, angleRadians->x);
	*m = GEMultiplyMat4(GEMat4FromRotate(angleRadians->x, 1.0, 0.0, 0.0), *m);
}

void rotateShapeZ(mat4 *m, quat *q, vec3 *angleRadians, vec3 *angularVelocity, vec3 *angularAcceleration)
{
	*angularVelocity = GEAddVectors(*angularVelocity, GEMultiplyVectorAndScalar(*angularAcceleration, globalTimeStep));
	*angleRadians = GEAddVectors(*angleRadians, GEMultiplyVectorAndScalar(*angularVelocity, globalTimeStep));
	
	*q = GEQuaternionFromAxisAndAngle(*angularVelocity, angleRadians->z);
	*m = GEMultiplyMat4(GEMat4FromRotate(angleRadians->z, 0.0, 0.0, 1.0), *m);
}

void shape(void)
{
	mat4 rotateMatrix, rotateMatrix1, rotateMatrix2, rotateMatrix3;
	rotateMatrix = rotateMatrix1 = rotateMatrix2 = rotateMatrix3 = GEIdentityMat4();
	// These will be static once it works
	vec3 angularVelocityVector1 = GERadiansFromDegreesVector(GESetVector(20.0, 0.0, 0.0));
	vec3 angularVelocityVector2 = GERadiansFromDegreesVector(GESetVector(0.0, 0.0, 20.0));
	vec3 angularVelocityVector3 = GERadiansFromDegreesVector(GESetVector(20.0, 0.0, 0.0));
	// This is here so I can test only the angular velocity
	vec3 angularAccelerationVector = GESetVector(0.0, 0.0, 0.0);
	
	static vec3 angleEulerRadians, angleEulerRadians1, angleEulerRadians2, angleEulerRadians3;
	static quat quaternion, quaternion1, quaternion2, quaternion3;
	static char initialze = 0;
	
	if (initialze == 0)
	{
		angleEulerRadians = angleEulerRadians1 = angleEulerRadians2 = angleEulerRadians3 = GESetVector(0.0, 0.0, 0.0);
		quaternion = GESetQuaternion(0.0, 0.0, 0.0, 1.0);
		initialze = 1;
	}
	
	// Rotate
	rotateShapeX(&rotateMatrix1, &quaternion1, &angleEulerRadians1, &angularVelocityVector1, &angularAccelerationVector);
	rotateShapeZ(&rotateMatrix2, &quaternion2, &angleEulerRadians2, &angularVelocityVector2, &angularAccelerationVector);
	rotateShapeX(&rotateMatrix3, &quaternion3, &angleEulerRadians3, &angularVelocityVector3, &angularAccelerationVector);
	
	// Multiply matrices (works)
	rotateMatrix = GEMultiplyMat4(rotateMatrix1, rotateMatrix);
	rotateMatrix = GEMultiplyMat4(rotateMatrix2, rotateMatrix);
	rotateMatrix = GEMultiplyMat4(rotateMatrix3, rotateMatrix);
	
	// Multiply Quaternions? (does not work)
	quaternion = GETransformQuaternions(quaternion, quaternion1);
	quaternion = GETransformQuaternions(quaternion, quaternion2);
	quaternion = GETransformQuaternions(quaternion, quaternion3);
	
	vec3 angularVelocityVector = GEAxisFromQuaternion(quaternion);
	// Not sure if this is correct but I took out multipling globalTimeStep since it is already used in rotateShape*();
	angularVelocityVector = GEAddVectors(angularVelocityVector, angularAccelerationVector);
	angleEulerRadians = GEAddVectors(angleEulerRadians, angularVelocityVector);
	
	glColor3f(0.0, 1.0, 1.0);
	glPushMatrix();
		glTranslatef( 2.0, 0.0, 0.0);
		// This spins really fast
		glMultMatrixf(GEMat4FromEulerAngle(angleEulerRadians).array);
		glutWireSphere(2.0, 7.0, 7.0);
	glPopMatrix();
	
	glColor3f(1.0, 0.0, 0.0);
	glPushMatrix();
		glTranslatef(-2.0, 0.0, 0.0);
		glMultMatrixf(rotateMatrix.array); // (Works)
		glutWireSphere(2.0, 7.0, 7.0);
	glPopMatrix();
}

 

2 hours ago, Jonathan2006 said:

I am confused as to what should be a vector and what should be a quaternion in your post above.

Some code:


	quat curOrientation = object.worldSpaceOrientation;
	quat targetOrientation = curOrientation;
	targetOrientation *= quat::FromaxisAndAngle(vec(1,0,0), PI*0.5f);
	targetOrientation *= quat::FromaxisAndAngle(vec(0,1,0), PI*-0.3f);
	targetOrientation *= quat::FromaxisAndAngle(vec(0,0,1), PI*1.5f);
	// now calculate rotation from current to target
	quat diff;
	quat &qA = curOrientation; &quat qB = targetOrientation;
	{
	quat q;
    if (qA.Dot(qB) < 0.0f) // use dot product so we pick the shortest arc rotation, and invert either axis or angular part 
    {
        q[0] = qA[0]; q[1] = qA[1]; q[2] = qA[2];
        q[3] = -qA[3]; 
    }
    else
    {
        q[0] = -qA[0]; q[1] = -qA[1]; q[2] = -qA[2];
        q[3] = qA[3];
    }
            
    diff = qB * q;
	}
	// now convert diff to angular velocitiy
	float targetDuration = timestep;
	vec3 angVel (0,0,0);
	{
	// this is axis to angle conversation as you expect, but with some confusing 'fix close to zero angle issues' approach.
	    const float matchTolerance = 0.0001f;
    const float faultTolerance = 0.0005f;
	    quat q = diff;
    vec3 *omegaDir = (vec3*)&q;
    float sql = omegaDir.Dot(omegaDir);
    if (sql > (matchTolerance * matchTolerance)) 
    {
	    float length = sqrt (sql);
    if (q[3] < -1) q[3] = -1;
    if (q[3] >  1) q[3] =  1;
    angVel = (*omegaDir / length) * (2.0 * acos(q[3]));
    if (length < faultTolerance) angVel *= (length - matchTolerance) / (faultTolerance - matchTolerance);
	angVel *= targetDuration;
	}
	}
	

If we integrate angVel for targetDuration, we reach targetOrientation as desired, but we won't have zero velocity, so the body will keep spinning so overshooting the target. The easiest way to avoid this is to recalculate angVel every step, but then integrate only a fraction of it like angVel*0.3 This way the body will come to rest at desired orientation, but no longer we are able to tell how long it will take exactly.

I think my math lib has unusual rotation order, so when i write  diff = qB * q, you may need to write  diff = q * qB eventually.

 

 

2 hours ago, Jonathan2006 said:

vec3 angularVelocityVector = GEAxisFromQuaternion(quaternion); // Not sure if this is correct but I took out multipling globalTimeStep since it is already used in rotateShape*(); 

This seems wrong (you need to use angle and time as well), it should be:

vec3 angularVelocityVector = GEAxisFromQuaternion(quaternion) * GEAngleFromQuaternion(quaternion) * targetDuration;

which would be the replacement for my confusing quat to angular velocity code.  

By using timestep for duration as i did  you do the rotation as fast as possible - of course you can do it slower as well.

2 hours ago, Jonathan2006 said:

angularVelocityVector = GEAddVectors(angularVelocityVector, angularAccelerationVector);

Not sure what you intend with this.

2 hours ago, Jonathan2006 said:

angleEulerRadians = GEAddVectors(angleEulerRadians, angularVelocityVector);

Looks wrong too. You can not add angular velocity to euler angles. Euler angles are always a series of 3 rotations in order, while angular velocity encodes a single rotation about its unit direction, by axis and speed encoded in its length.

I guess this code is a try to integrate angular velocity? Here something simple i have, hope it helps:


	void IntegrateAgularVelocity (const sQuat &curOrn, const sVec3 &angvel, float timeStep, sQuat &predictedOrn)
{
    sVec3 axis;
    float fAngle = angvel.Length(); 
    
    if ( fAngle < float(0.00001) )
    {
        predictedOrn = curOrn;
        return;
    }
    else
    {
        axis = angvel * (sin (fAngle * timeStep * float(0.5)) / fAngle);
    }
    sQuat dorn (axis.x(), axis.y(), axis.z(), cos (fAngle * timeStep * float(0.5)));
    predictedOrn = dorn * curOrn;
    predictedOrn.Normalize();
}
	

 

 

Thank you again for all the help! I am still trying to wrap my head around quaternions. I just thought of something. Even if I use something like quatCur - quatOld should the angular oscillating stop if I normalize the quaternion like you did in your code? I think I understand your last post on quaternions but I am getting weird acceleration. I also added the time step back into the last integrator (hope that's the name). The reason I use the last angleEulerRadians and angularVelocityVector (integrator?) is to only test the final angular velocity. I created a gif of the rotations below. Also Here's the code I am using:


void integrateAngularVelocity(const quat curOrn, const vec3 angvel, const float timeStep, quat *predictedOrn)
{
    vec3 axis;
    float fAngle = GEMagnitudeVector(angvel);
    
    if (fAngle < 0.00001)
    {
        *predictedOrn = curOrn;
        return;
    }
    else
    {
        axis = GEMultiplyVectorAndScalar(angvel, sinf(fAngle * timeStep * 0.5) / fAngle);
    }
	
    quat dorn = GESetQuaternion(axis.x, axis.y, axis.z, cosf(fAngle * timeStep * 0.5));
    *predictedOrn = GEMultiplyQuaternions(dorn, curOrn);
    *predictedOrn = GENormalizeQuaternion(*predictedOrn);
}

void Shape(void)
{
	vec3 angularVelocityVector1 = GERadiansFromDegreesVector(GESetVector(20.0, 0.0, 0.0));
	vec3 angularVelocityVector2 = GERadiansFromDegreesVector(GESetVector(0.0, 0.0, 20.0));
	vec3 angularVelocityVector3 = GERadiansFromDegreesVector(GESetVector(20.0, 0.0, 0.0));

    static quat quaternion;

	integrateAngularVelocity(quaternion, angularVelocityVector1, globalTimeStep, &quaternion);
	integrateAngularVelocity(quaternion, angularVelocityVector2, globalTimeStep, &quaternion);
	integrateAngularVelocity(quaternion, angularVelocityVector3, globalTimeStep, &quaternion);
	
	vec3 angularVelocityVector = GEMultiplyVectorAndScalar(GEAxisFromQuaternion(quaternion), GEAngleFromQuaternion(quaternion));

	angularVelocityVector = GEAddVectors(angularVelocityVector, GEMultiplyVectorAndScalar(angularAccelerationVector, globalTimeStep));
	angleEulerRadians = GEAddVectors(angleEulerRadians, GEMultiplyVectorAndScalar(angularVelocityVector, globalTimeStep));
	
	glColor3f(0.0, 1.0, 1.0);
	glPushMatrix();
		glTranslatef( 2.0, 0.0, 0.0);
		// This spins really fast
		glMultMatrixf(GEMat4FromEulerAngle(angleEulerRadians).array);
		glutWireSphere(2.0, 7.0, 7.0);
	glPopMatrix();
}

It's not the best quality but the red sphere is the correct rotation and the cyan is the quaternion code above.

motion.gif.7bf43e77fb23a3db1e0357ea5f3f4118.gif

This topic is closed to new replies.

Advertisement