更快的四元数向量乘法不起作用

3
我需要一个更快的四元数-向量乘法程序来完善我的数学库。目前我正在使用标准公式v' = qv(q^-1),这与用四元数制作矩阵后将其与向量相乘得到的结果相同,因此我对其正确性有信心。
到目前为止,我已经实现了三种替代方案以提高速度:
#1,我不知道这个方法是从哪里得到的:
v' = (q.xyz * 2 * dot(q.xyz, v)) + (v * (q.w*q.w - dot(q.xyz, q.zyx))) + (cross(q.xyz, v) * q.w * w)

实现为:

vec3 rotateVector(const quat& q, const vec3& v)
{
    vec3 u(q.x, q.y, q.z);
    float s = q.w;

    return vec3(u * 2.0f * vec3::dot(u, v))
    + (v * (s*s - vec3::dot(u, u)))
    + (vec3::cross(u, v) * s * 2.0f);
}

#2,感谢这篇精彩的博客提供思路。

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

实现方式:

__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));
}

第三点,来自多个来源,

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

实现方式:

__m128 rotateVector(__m128 q, __m128 v)
{
    //return v + 2.0 * cross(cross(v, q.xyz) + q.w * 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)));
}

所有三种方法都产生了错误的结果。然而,我注意到了一些有趣的模式。首先,1号和2号方法产生相同的结果。如果将矩阵转置,则3号方法产生与将向量乘以派生矩阵得到的结果相同(我是偶然发现这个规律的,之前我的四元数转矩阵代码假设矩阵是行主序,这是不正确的)。
我的四元数数据存储定义如下:
union
{
    __m128 data;
    struct { float x, y, z, w; };
    float f[4];
};

我的实现有问题吗,还是我漏掉了什么重要的东西?


参见:https://dev59.com/EH3aa4cB1Zd3GeqPa0fU - Paul R
1
为什么你需要它更快呢?你已经证明了它是代码的瓶颈吗?你能详细描述一下你在做什么吗?如果你需要计算多个独立的四元数,那么使用SIMD(如果你使用正确)可能会更快。然而,如果你正在计算一系列依赖四元数,那么SIMD可能实际上会更慢。 - Z boson
我需要它更快,因为我在所有地方都使用四元数,包括在glsl中。我现在不关心simd vs非simd,问题是替代算法根本不起作用。 - Haydn V. Harach
1
我发现了一个问题,#1和#2有相同的名称,并且在我想要使用#2时错误地将四元数转换为使用#1。 #2实际上产生了正确的结果,并且目前是我最快的实现。现在我只需要一个在glsl中工作的版本。 - Haydn V. Harach
我也成功地通过改变叉积的顺序使#3起作用。 - Haydn V. Harach
1个回答

1
主要问题,如果您想通过四元数旋转三维向量,则需要计算旋转矩阵的所有9个标量。在您的示例中,旋转矩阵的计算是隐式的。计算顺序可能不是最优的。
如果您从四元数生成3x3矩阵并乘以向量,则应具有相同数量的算术运算(请参见底部的代码)。
我推荐以下:
1. 尝试生成3x3矩阵并乘以您的向量,测量速度并与以前进行比较。 2. 分析显式解,并尝试为自定义体系结构进行优化。 3. 尝试实现替代四元数乘法,并从方程q * v * q'派生乘法。
// ============ 替代乘法伪代码
 /**
      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's many bugs, have to be rewriten.
    */
    inline xxquaternion mul_alt( const xxquaternion& q) const {
        float t0 = (x-y)*(q.y-q.x);
        float t1 = (w+z)*(q.w+q.z);
        float t2 = (w-z)*(q.y+q.x);
        float t3 = (x+y)*(q.w-q.z);
        float t4 = (x-z)*(q.z-q.y);
        float t5 = (x+z)*(q.z+q.y);
        float t6 = (w+y)*(q.w-q.x);
        float t7 = (w-y)*(q.w+q.x);

        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 = q.w,  bi = q.z, bj = q.y, bk =q.x;
        float t0 = (k-j)*(bj-bk);
        float t1 = (r+i)*(br+bi);
        float t2 = (r-i)*(bj+bk);
        float t3 = (k+j)*(br-bi);
        float t4 = (k-i)*(bi-bj);
        float t5 = (k+i)*(bi+bj);
        float t6 = (r+j)*(br-bk);
        float t7 = (r-j)*(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 );
*/
    }

//============ 显式向量旋转变体

/**
      rotate vector by quaternion
    */
    inline vector3 rotate(const vector3& v)const{
        xxquaternion q(v.x * w + v.z * y - v.y * z,
                       v.y * w + v.x * z - v.z * x,
                       v.z * w + v.y * x - v.x * y,
                       v.x * x + v.y * y + v.z * z);

        return vector3(w * q.x + x * q.w + y * q.z - z * q.y,
                       w * q.y + y * q.w + z * q.x - x * q.z,
                       w * q.z + z * q.w + x * q.y - y * q.x)*( 1.0f/norm() );

        // 29  multiplications, 20 addidtions, 4 variables
        // 5 
       /*
        // refrence implementation  
        xxquaternion r = (*this)*xxquaternion(v.x, v.y, v.z, 0)*this->inverted();
        return vector3( r.x, r.y, r.z ); 
       */

       /*
        // alternative implementation
        float wx, wy, wz, xx, yy, yz, xy, xz, zz, x2, y2, z2;
        x2 = q.x + q.x; y2 = q.y + q.y; z2 = q.z + q.z;

        xx = q.x * x2;   xy = q.x * y2;   xz = q.x * z2;
        yy = q.y * y2;   yz = q.y * z2;   zz = q.z * z2;
        wx = q.w * x2;   wy = q.w * y2;   wz = q.w * z2;

        return vector3( v.x  - v.x * (yy + zz) + v.y * (xy - wz) + v.z * (xz + wy),
                        v.y  + v.x * (xy + wz) - v.y * (xx + zz) + v.z * (yz - wx),
                        v.z  + v.x * (xz - wy) + v.y * (yz + wx) - v.z * (xx + yy) )*( 1.0f/norm() );
        // 18 multiplications, 21 addidtions, 12 variables
       */
    };

祝你好运。


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接