Make QuatToMat faster and more accurate

The better accuracy is for specific cases (90 degree rotations around a
main axis: the matrix element for that axis is now 1 instead of
0.99999994). The speedup comes from doing fewer additions (multiply
seems to be faster than add for fp, at least in this situation).
This commit is contained in:
Bill Currie 2019-07-23 08:52:15 +09:00
parent 1f24a1408a
commit e164002050
2 changed files with 80 additions and 18 deletions

View file

@ -303,7 +303,7 @@ QuatExp (const quat_t a, quat_t b)
VISIBLE void
QuatToMatrix (const quat_t q, vec_t *m, int homogenous, int vertical)
{
vec_t xx, xy, xz, xw, yy, yz, yw, zz, zw, ww;
vec_t xx, xy, xz, xw, yy, yz, yw, zz, zw;
vec_t *_m[4] = {
m + (homogenous ? 0 : 0),
m + (homogenous ? 4 : 3),
@ -311,28 +311,26 @@ QuatToMatrix (const quat_t q, vec_t *m, int homogenous, int vertical)
m + (homogenous ? 12 : 9),
};
xx = q[0] * q[0];
xy = q[0] * q[1];
xz = q[0] * q[2];
xw = q[0] * q[3];
xx = 2 * q[0] * q[0];
xy = 2 * q[0] * q[1];
xz = 2 * q[0] * q[2];
xw = 2 * q[0] * q[3];
yy = q[1] * q[1];
yz = q[1] * q[2];
yw = q[1] * q[3];
yy = 2 * q[1] * q[1];
yz = 2 * q[1] * q[2];
yw = 2 * q[1] * q[3];
zz = q[2] * q[2];
zw = q[2] * q[3];
ww = q[3] * q[3];
zz = 2 * q[2] * q[2];
zw = 2 * q[2] * q[3];
if (vertical) {
VectorSet (ww + xx - yy - zz, 2 * (xy + zw), 2 * (xz - yw), _m[0]);
VectorSet (2 * (xy - zw), ww - xx + yy - zz, 2 * (yz + xw), _m[1]);
VectorSet (2 * (xz + yw), 2 * (yz - xw), ww - xx - yy + zz, _m[2]);
VectorSet (1.0f - yy - zz, xy + zw, xz - yw, _m[0]);
VectorSet (xy - zw, 1.0f - xx - zz, yz + xw, _m[1]);
VectorSet (xz + yw, yz - xw, 1.0f - xx - yy, _m[2]);
} else {
VectorSet (ww + xx - yy - zz, 2 * (xy - zw), 2 * (xz + yw), _m[0]);
VectorSet (2 * (xy + zw), ww - xx + yy - zz, 2 * (yz - xw), _m[1]);
VectorSet (2 * (xz - yw), 2 * (yz + xw), ww - xx - yy + zz, _m[2]);
VectorSet (1.0f - yy - zz, xy - zw, xz + yw, _m[0]);
VectorSet (xy + zw, 1.0f - xx - zz, yz - xw, _m[1]);
VectorSet (xz - yw, yz + xw, 1.0f - xx - yy, _m[2]);
}
if (homogenous) {
_m[0][3] = 0;

View file

@ -200,6 +200,62 @@ fail:
return 0;
}
#define s05 0.70710678118654757
static struct {
quat_t q;
vec_t expect[9];
} quat_mat_tests[] = {
{{0, 0, 0, 1},
{1, 0, 0,
0, 1, 0,
0, 0, 1}},
{{1, 0, 0, 0},
{1, 0, 0,
0, -1, 0,
0, 0, -1}},
{{0, 1, 0, 0},
{-1, 0, 0,
0, 1, 0,
0, 0, -1}},
{{0, 0, 1, 0},
{-1, 0, 0,
0, -1, 0,
0, 0, 1}},
{{0.5, 0.5, 0.5, 0.5},
{0, 0, 1,
1, 0, 0,
0, 1, 0}},
{{s05, 0.0, 0.0, s05},
{1, 0, 0,
0, 5.96046448e-8, -0.99999994,
0, 0.99999994, 5.96046448e-8}},
};
#define num_quat_mat_tests (sizeof (quat_mat_tests) / sizeof (quat_mat_tests[0]))
static int
test_quat_mat(const quat_t q, const quat_t expect)
{
int i;
vec_t m[9];
QuatToMatrix (q, m, 0, 0);
for (i = 0; i < 9; i++)
if (m[i] != expect[i]) // exact tests here
goto fail;
return 1;
fail:
printf ("%11.9g %11.9g %11.9g %11.9g\n", QuatExpand (q));
printf ("%11.9g %11.9g %11.9g %11.9g %11.9g %11.9g\n",
VectorExpand (m + 0), VectorExpand (expect + 0));
printf ("%11.9g %11.9g %11.9g %11.9g %11.9g %11.9g\n",
VectorExpand (m + 3), VectorExpand (expect + 3));
printf ("%11.9g %11.9g %11.9g %11.9g %11.9g %11.9g\n",
VectorExpand (m + 6), VectorExpand (expect + 6));
return 0;
}
int
main (int argc, const char **argv)
{
@ -228,5 +284,13 @@ main (int argc, const char **argv)
if (!test_rotation3 (test_angles[i]))
res = 1;
}
for (i = 0; i < num_quat_mat_tests; i ++) {
vec_t *q = quat_mat_tests[i].q;
vec_t *expect = quat_mat_tests[i].expect;
if (!test_quat_mat (q, expect))
res = 1;
}
return res;
}