Quaternion transformation precision

qt.js

// This code looks at precision impact of transforming a vector repeatedly by a slightly-non-unit quaternion
// Slightly-non-unit quaternions are important: they result in the process of quaternion computations naturally
// Repeated transformations are important: they may occur during simulation or complex long chains of computation

// Note that because this code runs in JS in double precision, this doesn't model floating-point roundoff.

function applyQuaternion1( q, v ) {

	const x = v.x, y = v.y, z = v.z;
	const qx = q.x, qy = q.y, qz = q.z, qw = q.w;

	// calculate quat * vector

	const ix = qw * x + qy * z - qz * y;
	const iy = qw * y + qz * x - qx * z;
	const iz = qw * z + qx * y - qy * x;
	const iw = - qx * x - qy * y - qz * z;

	// calculate result * inverse quat

	const rx = ix * qw + iw * - qx + iy * - qz - iz * - qy;
	const ry = iy * qw + iw * - qy + iz * - qx - ix * - qz;
	const rz = iz * qw + iw * - qz + ix * - qy - iy * - qx;

	return { x: rx, y: ry, z: rz };

}

function applyQuaternion2( q, v ) {

	const vx = v.x, vy = v.y, vz = v.z;
	const qx = q.x, qy = q.y, qz = q.z, qw = q.w;

	// t = 2q x v
	const tx = 2 * ( qy * vz - qz * vy );
	const ty = 2 * ( qz * vx - qx * vz );
	const tz = 2 * ( qx * vy - qy * vx );

	// v + w t + q x t
	const rx = vx + qw * tx + qy * tz - qz * ty;
	const ry = vy + qw * ty + qz * tx - qx * tz;
	const rz = vz + qw * tz + qx * ty - qy * tx;

	return { x: rx, y: ry, z: rz };
}

var axis = { x: 1 / Math.sqrt(2), y: 0, z: 1 / Math.sqrt(2) };
var angle = 0.12345;

var q = { x: axis.x * Math.sin(angle/2), y: axis.y * Math.sin(angle/2), z: axis.z * Math.sin(angle/2), w: Math.cos(angle/2) };

// Validate quaternion is unit length
console.assert(Math.abs(q.x * q.x + q.y * q.y + q.z * q.z + q.w * q.w - 1) < 1e-6);

// Copy quaternion to a unit-length version
var qu = { x: q.x, y: q.y, z: q.z, w: q.w };

// Change quaternion to non-unit-length; we choose 1+1e-5 to simulate single precision round-off error that would be common in practice outside of JS
var err = 1e-5;
q.x *= 1 + err;
q.y *= 1 + err;
q.z *= 1 + err;
q.w *= 1 + err;

// Now take a vector and multiply it repeatedly by the quaternion using both methods
// Also make a canonical copy from the unit-length quaternion
var v1 = { x: 1, y: 0, z: 0 };
var v2 = { x: 1, y: 0, z: 0 };
var vu = { x: 1, y: 0, z: 0 };

for (var i = 0; i < 1000; ++i)
{
	v1 = applyQuaternion1(q, v1);
	v2 = applyQuaternion2(q, v2);
	vu = applyQuaternion1(qu, vu); // here it doesn't matter which method we use
}

// Print all three vectors; observe that the error in the length of the second vector is much smaller
var l1 = Math.sqrt(v1.x * v1.x + v1.y * v1.y + v1.z * v1.z);
var l2 = Math.sqrt(v2.x * v2.x + v2.y * v2.y + v2.z * v2.z);
var lu = Math.sqrt(vu.x * vu.x + vu.y * vu.y + vu.z * vu.z);

console.log(v1.x, v1.y, v1.z, "length", l1);
console.log(v2.x, v2.y, v2.z, "length", l2);
console.log(vu.x, vu.y, vu.z, "length", lu);

// Print all three vectors, normalized
console.log(v1.x / l1, v1.y / l1, v1.z / l1, "deviation", (v1.x * vu.x + v1.y * vu.y + v1.z * vu.z) / l1 / lu);
console.log(v2.x / l2, v2.y / l2, v2.z / l2, "deviation", (v2.x * vu.x + v2.y * vu.y + v2.z * vu.z) / l2 / lu);
console.log(vu.x / lu, vu.y / lu, vu.z / lu, "deviation", (vu.x * vu.x + vu.y * vu.y + vu.z * vu.z) / lu / lu);