Faster multiplication of the quaternion vector does not work

Question

I need a faster quaternion-vector multiplication procedure for my math library. Now I use the canonical v' = qv(q^-1) , which gives the same result as multiplying a vector by a matrix made from a quaternion, so I'm sure of its correctness.

So far I have implemented 3 alternative methods "faster":

# 1, I have no idea where I got this from:

 v' = (q.xyz * 2 * dot(q.xyz, v)) + (v * (qw*qw - dot(q.xyz, q.zyx))) + (cross(q.xyz, v) * qw * w) 

Implemented as:

 vec3 rotateVector(const quat& q, const vec3& v) { vec3 u(qx, qy, qz); float s = qw; return vec3(u * 2.0f * vec3::dot(u, v)) + (v * (s*s - vec3::dot(u, u))) + (vec3::cross(u, v) * s * 2.0f); } 

# 2, courtesy of this beautiful blog

 t = 2 * cross(q.xyz, v); v' = v + qw * t + cross(q.xyz, t); 

Implemented as:

 __m128 rotateVector(__m128 q, __m128 v) { __m128 temp = _mm_mul_ps(vec4::cross(q, v), _mm_set1_ps(2.0f)); return _mm_add_ps( _mm_add_ps(v, _mm_mul_ps(_mm_shuffle_ps(q, q, _MM_SHUFFLE(3, 3, 3, 3)), temp)), vec4::cross(q, temp)); } 

And No. 3, from numerous sources,

 v' = v + 2.0 * cross(cross(v, q.xyz) + qw * v, q.xyz); 

implemented as:

 __m128 rotateVector(__m128 q, __m128 v) { //return v + 2.0 * cross(cross(v, q.xyz) + qw * v, q.xyz); return _mm_add_ps(v, _mm_mul_ps(_mm_set1_ps(2.0f), vec4::cross( _mm_add_ps( _mm_mul_ps(_mm_shuffle_ps(q, q, _MM_SHUFFLE(3, 3, 3, 3)), v), vec4::cross(v, q)), q))); } 

All 3 of them give incorrect results. However, I noticed some interesting patterns. First of all, # 1 and # 2 give the same result. # 3 gives the same result as me, getting from multiplying a vector by a derivative matrix if the specified matrix is ​​transposed (I discovered this by accident, earlier my square matrix code accepted matrices of strings that were incorrect).

Saving my quaternion data is defined as:

 union { __m128 data; struct { float x, y, z, w; }; float f[4]; }; 

Are my implementations incorrect or am I not seeing something here?

Answer 1

The main problem, if you want to rotate the 3d vector in the quaternion, you need to calculate all 9 scalars of the rotation matrix. In your examples, the calculation of the rotation matrix is ​​IMPLICIT. The calculation procedure may not be optimal.

If you are creating a 3x3 matrix from a quaternion and a multiplication vector, you should have the same number of arithmetic operations (@see code below).

What I recommend.

• Try to create a 3x3 matrix and multiply the vector, measure the speed and compare with the previous one.

• Analyze the explicit solution and try to optimize for the user architecture.

• try to implement alternative quaternion multiplication and derivative multiplication from the q * v * q 'equation.

// ============= pseudocode of alternative multiplication

  /** alternative way of quaternion multiplication, can speedup multiplication for some systems (PDA for example) http://mathforum.org/library/drmath/view/51464.html http://www.lboro.ac.uk/departments/ma/gallery/quat/src/quat.ps in provided code by url many bugs, have to be rewriten. */ inline xxquaternion mul_alt( const xxquaternion& q) const { float t0 = (xy)*(qy-qx); float t1 = (w+z)*(q.w+qz); float t2 = (wz)*(q.y+qx); float t3 = (x+y)*(qw-qz); float t4 = (xz)*(qz-qy); float t5 = (x+z)*(q.z+qy); float t6 = (w+y)*(qw-qx); float t7 = (wy)*(q.w+qx); float t8 = t5 + t6 + t7; float t9 = (t4 + t8)*0.5; return xxquaternion ( t3+t9-t6, t2+t9-t7, t1+t9-t8, t0+t9-t5 ); // 9 multiplications 27 addidtions 8 variables // but of couse we can clean 4 variables /* float r = w, i = z, j = y, k =x; float br = qw, bi = qz, bj = qy, bk =qx; float t0 = (kj)*(bj-bk); float t1 = (r+i)*(br+bi); float t2 = (ri)*(bj+bk); float t3 = (k+j)*(br-bi); float t4 = (ki)*(bi-bj); float t5 = (k+i)*(bi+bj); float t6 = (r+j)*(br-bk); float t7 = (rj)*(br+bk); float t8 = t5 + t6 + t7; float t9 = (t4 + t8)*0.5; float rr = t0+t9-t5; float ri = t1+t9-t8; float rj = t2+t9-t7; float rk = t3+t9-t6; return xxquaternion ( rk, rj, ri, rr ); */ } 

// ============= explicit vector rotation options

 /** rotate vector by quaternion */ inline vector3 rotate(const vector3& v)const{ xxquaternion q(vx * w + vz * y - vy * z, vy * w + vx * z - vz * x, vz * w + vy * x - vx * y, vx * x + vy * y + vz * z); return vector3(w * qx + x * qw + y * qz - z * qy, w * qy + y * qw + z * qx - x * qz, w * qz + z * qw + x * qy - y * qx)*( 1.0f/norm() ); // 29 multiplications, 20 addidtions, 4 variables // 5 /* // refrence implementation xxquaternion r = (*this)*xxquaternion(vx, vy, vz, 0)*this->inverted(); return vector3( rx, ry, rz ); */ /* // alternative implementation float wx, wy, wz, xx, yy, yz, xy, xz, zz, x2, y2, z2; x2 = qx + qx; y2 = qy + qy; z2 = qz + qz; xx = qx * x2; xy = qx * y2; xz = qx * z2; yy = qy * y2; yz = qy * z2; zz = qz * z2; wx = qw * x2; wy = qw * y2; wz = qw * z2; return vector3( vx - vx * (yy + zz) + vy * (xy - wz) + vz * (xz + wy), vy + vx * (xy + wz) - vy * (xx + zz) + vz * (yz - wx), vz + vx * (xz - wy) + vy * (yz + wx) - vz * (xx + yy) )*( 1.0f/norm() ); // 18 multiplications, 21 addidtions, 12 variables */ }; 

Good luck.