Java软光栅渲染器-四元数

目标

定义四元数,并实现四元数的运算。

源代码:

实现

关于四元数的意义、四元数与复数、四元数的共轭等概念,网上有太资料来学习,本文就不再多写了。我在实现四元数时,主要参考《3D数学基础:图形与游戏开发》第10章的内容。

四元数记法

四元数包含一个标量分量和一个3D向量分量。经常记标量分量为 w,记向量分量为单一 v 或分开的 x, y ,z。两种记法的分别如下:

  • [w, v]
  • [w, x, y, z]

在代码中,我把 w 写在 x, y, z 的后面,只要不影响计算即可。

/**
 * 四元数
 * 
 * @author yanmaoyuan
 *
 */
public class Quaternion {

    public float x, y, z, w;

    public Quaternion() {
        x = y = z = 0;
        w = 1;
    }

    public Quaternion(float x, float y, float z, float w) {
        this.x = x;
        this.y = y;
        this.z = z;
        this.w = w;
    }

    public Quaternion(Quaternion q) {
        this.x = q.x;
        this.y = q.y;
        this.z = q.z;
        this.w = q.w;
    }
}

角轴旋转

已知由单位向量 n 定义的旋转轴,以及旋转的弧度 θ,用四元数来记录这个旋转如下:

q = [cos(θ/2), sin(θ/2)n]

把单位向量n展开为 x, y, z 三个分量,记为:

q = [cos(θ/2), sin(θ/2) * x, sin(θ/2) * y, sin(θ/2) * z]

在右手坐标系中,可以使用右手定则来确定 n 和 θ 的正方向:右手握拳,让大拇指伸出朝向 n 的正方向;此时其余四指的方向(逆时针)就是 θ 的正方向。

代码实现如下:

/**
 * 使用角轴对来计算四元数
 * 
 * @param axis 由单位向量定义的旋转轴
 * @param angle 旋转角(弧度制)
 */
public Quaternion(Vector3f axis, float angle) {
    axis.normalizeLocal();

    float cosHalfAngle = (float) Math.cos(angle * 0.5);
    float sinHalfAngle = (float) Math.sin(angle * 0.5);

    w = cosHalfAngle;
    x = axis.x * sinHalfAngle;
    y = axis.y * sinHalfAngle;
    z = axis.z * sinHalfAngle;
}

/**
 * 使用角轴对来计算四元数
 * @param axis 由单位向量定义的旋转轴
 * @param angle 旋转角(弧度制)
 * @return
 */
public Quaternion fromAxisAngle(Vector3f axis, float angle) {
    axis.normalizeLocal();

    float sinHalfAngle = (float) Math.sin(angle * 0.5);
    float cosHalfAngle = (float) Math.cos(angle * 0.5);

    w = cosHalfAngle;

    x = axis.x * sinHalfAngle;
    y = axis.y * sinHalfAngle;
    z = axis.z * sinHalfAngle;

    return this;
}

如果是绕X、Y、Z轴旋转,实现还会更简单一点。

/**
 * 绕X轴旋转
 * @param angle
 * @return
 */
public Quaternion rotateX(float angle) {
    float sinHalfAngle = (float) Math.sin(angle * 0.5);
    float cosHalfAngle = (float) Math.cos(angle * 0.5);

    w = cosHalfAngle;

    x = sinHalfAngle;
    y = 0;
    z = 0;

    return this;
}

/**
 * 绕Y轴旋转
 * @param angle
 * @return
 */
public Quaternion rotateY(float angle) {
    float sinHalfAngle = (float) Math.sin(angle * 0.5);
    float cosHalfAngle = (float) Math.cos(angle * 0.5);

    w = cosHalfAngle;

    x = 0;
    y = sinHalfAngle;
    z = 0;

    return this;
}

/**
 * 绕Z轴旋转
 * @param angle
 * @return
 */
public Quaternion rotateZ(float angle) {
    float sinHalfAngle = (float) Math.sin(angle * 0.5);
    float cosHalfAngle = (float) Math.cos(angle * 0.5);

    w = cosHalfAngle;

    x = 0;
    y = 0;
    z = sinHalfAngle;

    return this;
}

记住,q 的 w 分量与 θ 有关,但它们不是同一回事。同样,v与n也有关系但不完全相同。

负四元数

四元数能求负:将每个分量都变负即可。

/**
 * 负四元数
 * 
 * @return
 */
public Quaternion negate() {
    return new Quaternion(-x, -y, -z, -w);
}

/**
 * 负四元数
 * 
 * @return
 */
public Quaternion negateLocal() {
    x = -x;
    y = -y;
    z = -z;
    w = -w;
    return this;
}

q 和 -q 代表的实际角位移是相同,很奇怪吧!如果我们将 θ 加上360°的倍数,不会改变 q 代表的角位移,但它使 q 的四个分量都变负了。因此,3D中的任意角位移都有两种不同的四元数表示方法,它们互相为负。

单位四元数

几何上,存在两个“单位”四元数,它们代表没有角位移:[1, 0] 和 [-1, 0](注意这里的 0 代表零向量)。

  • 当 θ 是360°的偶数倍时,有第一种形式,cos(θ/2)=1;
  • 当 θ 是360°的奇数倍时,有第二种形式,cos(θ/2)=-1。

这两种情况下,都有 sin(θ/2)=0,所以 n 的值无关紧要。它的意义在于:当旋转角 θ 是360°的整数倍时,方位并没有改变,并且旋转轴也是无关紧要的。

数学上,实际只有一个单位四元数:[1, 0]。用任意四元数q乘以单位四元数 [1, 0],结果任是 q。任意四元数 q 乘以另一个“几何单位”四元数 [-1, 0] 时得到 -q。几何上,因为 q 和 -q 代表的角位移相同,可认为结果是相同的。但在数学上,q 和 -q 不相等,所以 [-1, 0] 并不是“真正”的单位四元数。

代码如下:

// 单位四元数
public final static Quaternion IDENTITY = new Quaternion(0, 0, 0, 1);

/**
 * 将四元数设为 (0, 0, 0, 1)
 */
public void loadIdentity() {
    x = y = z = 0;
    w = 1;
}

/**
 * 判断是否为单位四元数
 * 
 * @return true
 */
public boolean isIdentity() {
    if (x == 0 && y == 0 && z == 0 && w == 1) {
        return true;
    } else {
        return false;
    }
}

四元数的模

四元数求模公式与向量求模是类似的。

/**
 * 求长度的平方
 * 
 * @return
 */
public float lengthSquared() {
    return x * x + y * y + z * z + w * w;
}

/**
 * 求四元数的模,或者叫长度。
 * 
 * @return
 */
public float length() {
    return (float) Math.sqrt(x * x + y * y + z * z + w * w);
}

若旋转轴n是单位向量,就有 |q| = 1。如果为了使用四元数来表示方位,我们仅使用符合这个规则的单位四元数

规范化四元数

规范化四元数的定义与向量的规范化类似。

/**
 * 规范化四元数
 * 
 * @return
 */
public Quaternion normalize() {
    float length = x * x + y * y + z * z + w * w;
    if (length != 1f && length != 0f) {
        length = (float) (1.0 / Math.sqrt(length));
        return new Quaternion(length * x, length * y, length * z, length * w);
    }
    return new Quaternion(x, y, z, w);
}

/**
 * 规范化四元数
 * 
 * @return
 */
public Quaternion normalizeLocal() {
    float length = x * x + y * y + z * z + w * w;
    if (length != 1f && length != 0f) {
        length = (float) (1.0 / Math.sqrt(length));
        x *= length;
        y *= length;
        z *= length;
        w *= length;
    }
    return this;
}

四元数的共轭

四元数的共轭记作q*,可以通过让四元数的向量部分变负来获得。

q* = [w, v]* = [w, -v]
   = [w, -x, -y, -z]

代码实现如下:

/**
 * 求共轭四元数
 * 
 * @return
 */
public Quaternion conjugate() {
    return new Quaternion(-x, -y, -z, w);
}

/**
 * 求共轭四元数
 * 
 * @return
 */
public Quaternion conjugateLocal() {
    x = -x;
    y = -y;
    z = -z;
    return this;
}

它的几何意义是绕旋转轴 n 的反方向旋转 θ 角。

四元数的逆

四元数的逆记作 q^-1,定义为四元数的共轭除以它的模。

q^-1 = q* / |q|

代码实现如下:

/**
 * 求四元数的逆。 四元数的逆 = 四元数的共轭 / 四元数的模
 * 
 * @return
 */
public Quaternion inverse() {
    float length = x * x + y * y + z * z + w * w;
    if (length != 1f && length != 0f) {
        length = (float) (1.0 / Math.sqrt(length));
        return new Quaternion(-x * length, -y * length, -z * length, w * length);
    }
    // 单位四元数的逆和共轭是相等的。
    return new Quaternion(-x, -y, -z, w);
}

/**
 * 求四元数的逆。 四元数的逆 = 四元数的共轭 / 四元数的模
 * 
 * @return
 */
public Quaternion inverseLocal() {
    float length = x * x + y * y + z * z + w * w;
    if (length != 1f && length != 0f) {
        length = (float) (1.0 / Math.sqrt(length));
        x *= length;
        y *= length;
        z *= length;
    }
    x = -x;
    y = -y;
    z = -z;
    return this;
}

一个四元数 q 乘以它的逆 q^-1,即可得到单位四元数[1, 0]。

上面的公式是四元数求逆的正式定义,但是由于我们只使用单位四元数,所以四元数的逆和共轭是相等的。

四元数乘法

四元数的乘法是根据其复数解释来相乘的,推导过程就不写了,直接上代码。

/**
 * 四元数乘法,或者叫叉乘。
 * 
 * @param q
 * @return
 */
public Quaternion mult(Quaternion q) {
    float qw = q.w, qx = q.x, qy = q.y, qz = q.z;

    float rw = w * qw - x * qx - y * qy - z * qz;
    float rx = w * qx + x * qw + y * qz - z * qy;
    float ry = w * qy + y * qw + z * qx - x * qz;
    float rz = w * qz + z * qw + x * qy - y * qx;

    return new Quaternion(rx, ry, rz, rw);
}

/**
 * 四元数乘法,或者叫叉乘。
 * 
 * @param q
 * @return
 */
public Quaternion multLocal(Quaternion q) {
    float qw = q.w, qx = q.x, qy = q.y, qz = q.z;

    float rw = w * qw - x * qx - y * qy - z * qz;
    float rx = w * qx + x * qw + y * qz - z * qy;
    float ry = w * qy + y * qw + z * qx - x * qz;
    float rz = w * qz + z * qw + x * qy - y * qx;

    x = rx;
    y = ry;
    z = rz;
    w = rw;
    return this;
}

四元数的乘法满足结合律,但不满足交换律。

a(bc) = (ab)c
ab =/= ba

两个四元数乘积的模,等于模的乘积。

|q1×q2| = |q1||q2|

这个结论非常重要,因为它保证了两个单位四元数相乘的结果还是单位四元数。

四元数乘积的逆,等于各个四元数的逆以相反的顺序相乘。

(a × b)^-1 = b^-1 x a^-1

旋转三维向量

若将一个三维空间中的点v(x, y, z) “扩展” 到四元数空间,定义为四元数 p = [0, v]。设 q 为我们讨论的旋转四元数形式 q = [cos(θ/2), sin(θ/2)*n]。执行下面的乘法,可以使三维点 p 绕 n 旋转:

p' = q x p x q^-1

由于 q 是单位四元数,可以用 q 的共轭来代替 q 的逆,写作:

p' = q x p x q*

若 p = (0, x, y, z),q = (w', x', y', z'),q* = (w', -x', -y', -z')就有:

p' = (w', x', y', z') × (0, x, y, z) × (w', -x', -y', -z')

三个四元数连续叉乘,真让人头大。。

我试着在纸上推导结果,算了一整页,才推导出来p'.w = 0。

直接上代码吧。

/**
 * 使用这个四元数来旋转一个三维向量。
 * 
 * @param v
 * @return
 */
public Vector3f mult(Vector3f v) {
    if (v.x == 0 && v.y == 0 && v.z == 0) {
        return new Vector3f(0, 0, 0);
    } else {
        float vx = v.x, vy = v.y, vz = v.z;
        float rx = w * w * vx + 2 * y * w * vz - 2 * z * w * vy + x * x * vx + 2 * y * x * vy + 2 * z * x * vz - z * z * vx - y * y * vx;
        float ry = 2 * x * y * vx + y * y * vy + 2 * z * y * vz + 2 * w * z * vx - z * z * vy + w * w * vy - 2 * x * w * vz - x * x * vy;
        float rz = 2 * x * z * vx + 2 * y * z * vy + z * z * vz - 2 * w * y * vx - y * y * vz + 2 * w * x * vy - x * x * vz + w * w * vz;
        return new Vector3f(rx, ry, rz);
    }
}

如果要进行多次旋转,情况如下:

p' = b(apa*)b*
   = (ba)p(a*b*) // 结合律
   = (ba)p(ba)*  // 四元数乘积的逆,等于各个四元数的逆以相反的顺序相乘。

注意,先进性a旋转再进行b旋转,等价于执行乘积ba代表的单一旋转。因此,四元数乘法能用来连接多次旋转,这个矩阵乘法的效果一样。

根据四元数的标准定义,这个旋转式以从右向左的顺序发生的。这非常不幸,因为它迫使我们以“由里向外”的顺序连接多次旋转,这和以矩阵形式作同样的运算是不同的(至少在使用行向量时时是不同的)。

乘法的变种

针对这个公式导致的“顺序颠倒”问题,在本书文本和代码中,我们将违背标准定义,以相反的顺序来定义四元数乘法。注意,仅仅叉乘部分收到了影响。

——《3D数学基础:图形与游戏开发》,第10.4.8节,P 149。

按照这本书上的做法,四元数乘法改成了:

  [w1, x1, y1, z1] [w2, x2, y2, z2]
= [w1, v1] [w2, v2]
= [w1w2 - v1.v2   w1v2 + w2v1 + v1 x v2] // <--注意最右叉乘的顺序
  w1w2 - x1x2 - y1y2 - z1z2
= w1x2 + x1w2 + y1z2 - z1y2
  w1y2 + y1w2 + z1x2 - x1z2
  w1z2 + z1w2 + x1y2 - y1x2

而标准定义是:

  [w1, x1, y1, z1] [w2, x2, y2, z2]
= [w1, v1] [w2, v2]
= [w1w2 - v1.v2   w1v2 + w2v1 + v2 x v1]
  w1w2 - x1x2 - y1y2 - z1z2
= w1x2 + x1w2 + z1y2 - y1z2
  w1y2 + y1w2 + x1z2 - z1x2
  w1z2 + z1w2 + y1x2 - x1y2

经过这么一改,四元数的几何性质没变,但是乘法方式改变了:

p' = q^-1 x p x q
   = q* x p x q

多个旋转连续相乘变成了自左向右,与旋转发生的顺序一致:

p' = b*(a*pa)b
   = (b*a*)p(ab)
   = (ab)* p (ab)

在我的程序中,依然按照四元数的标准方法来进行计算。

说句题外话,我似乎有点明白为什么解析某些游戏的3D动画旋转数据时总是出错了。。也许在别人的程序中,是使用这种方式计算的四元数,而jMonkeyEngine使用的是标准方式。

四元数的“差”

利用四元数的乘法和逆,就能够计算两个四元数的“差”。“差”被定义为一个方位到另一个方位的角位移。换句话说,给定方位a和b,能够计算从a旋转到b的角位移d。用四元数等式更加紧凑地表示为:ad=b。

现在来求d。

ad = b
(a^-1)ad = (a^-1)b
[1 0] d = (a^-1)b
d = a* b

这样就有了求的代表一个方位到另一个方位角位移的四元数方法。在四元数插值时会用到这种方法。

/**
 * 计算两个四元数之间的“差”。
 * 
 * 已知四元数a和b,计算从a旋转到b的角位移d。
 * a * d = b
 *     d = a的逆 * b
 * 
 * @param q
 * @return
 */
public Quaternion delta(Quaternion q) {
    return this.inverse().mult(q);
}

四元数点乘

四元数的点乘与向量点乘非常相似。

/**
 * 四元数点乘。这个值越大,说明两个四元数的旋转角度越接近。
 * 当两个四元数都是单位四元数时,返回值是它们之间夹角的余弦值。
 * 
 * @param q
 * @return
 */
public float dot(Quaternion q) {
    return x * q.x + y * q.y + z * q.z + w * q.w;
}

注意,和向量点乘一样,其结果是标量。对于单位四元数a和b,有 -1 <= a.b <= 1。通常我们只关心 a.b 的绝对值,因为 a.b = -(a.-b),所以 b 和 -b 代表相同的角位移。

四元数点乘的几何解释类似于向量点乘的几何解释。四元数点乘 a.b 的绝对值越大,a 和 b代表的角位移越“相似”。

四元数插值

我跳过了四元数的对数、指数和标量乘法运算,直接到了这里。因为slerp运算实在是太重要了,因为它可以在两个四元数间平滑插值。slerp运算避免了欧拉角插值的所有问题。

slerp 运算是一个三元操作。前两个操作数是两个四元数,将在它们中间插值。设这两个“开始”和“结束”四元数分别为q0和q1。插值参数设为变量t,t在0到1之间变化。slerp函数:slerp(q0, q1, t)

代码如下:

/**
 * 球面线性插值(Spherical Linear intERPolation)
 *
 * @param src
 * @param dest
 * @param t
 */
public Quaternion slerp(Quaternion src, Quaternion dest, float t) {
    if (src.x == dest.x && src.y == dest.y && src.z == dest.z && src.w == dest.w) {
        return new Quaternion(src);
    }

    // 用点乘计算两四元数夹角的 cos 值
    float cos = src.dot(dest);

    // 如果点乘为负,则反转一个四元数以取得短的4D“弧”
    if (cos < 0.0f) {
        cos = -cos;
        dest = dest.negate();
    }

    // 计算两个四元数的插值系数
    float srcFactor = 1 - t;
    float destFactor = t;

    // 检查它们是否过于接近以避免除零
    if (cos > 0.999f) {
        // 非常接近 -- 即线性插值
        srcFactor = 1 - t;
        destFactor = t;
    } else {
        // 计算两个四元数之间的夹角
        float angle = (float) Math.acos(cos);
        // 计算分母的倒数,这样就只需要一次除法
        float invSin = 1f / (float) Math.sin(angle);

        // 计算插值系数
        srcFactor = (float) Math.sin((1 - t) * angle) * invSin;
        destFactor = (float) Math.sin((t * angle)) * invSin;
    }

    // 插值
    float rx = (srcFactor * src.x) + (destFactor * dest.x);
    float ry = (srcFactor * src.y) + (destFactor * dest.y);
    float rz = (srcFactor * src.z) + (destFactor * dest.z);
    float rw = (srcFactor * src.w) + (destFactor * dest.w);

    // 返回插值结果
    return new Quaternion(rx, ry, rz, rw);
}

总结

目标达成。

虽然四元数很难理解,但毕竟我还是实现了四元数的常用运算。具体代码写得对不对,还得看以后的运行结果。