3D图形:矩阵、欧拉角、四元数与方位的故事

Posted 2023-03-25

参考技术A

又研究了将近两个星期的3D图形到了我最想研究的地方了,因为欧拉角与四元数的原因导致OpenGL ES的研究进度变缓,研究完这一块,我将教大家如何使用OpenGL ES做一个自转加公转的正立方体.效果如下.

在说矩阵、欧拉角与四元数三种与角位移的关系之前,我们先来说说 方向、方位与角位移的区别.

在现实生活中,我们很少区分"方向"和"方位"的区别(非路痴观点),比如一个朋友来看望你,但是他可能在某一个公交站下车了,你去接他,但是找不到他,你急忙给他来一个电话"兄弟,你在哪个方向呢?"或者说是"兄弟,你在哪个方位呢?",如果不细细品味这两句话,其实感觉差异不是太大.通过一痛电话的扯,然后你们成功的面基了,但是你们却并不会在意"方向"和"方位"的区别.那么在几何中,这两者到底有什么差异呢?

这里我就盗用一下书上的例子,比如一个向量如果沿着自己的方向选择是不会改变自身任何属性的,如下图所示,因为向量是只有方向没有方位的.

那么对于一个物体,情况却是不一样的,一个物体如果朝向某一个方向的时候,然后自转,那么这个物体是会发生空间上的改变的,如下图一个锥体的自转,那么它的空间位置是发生改变的,也就是锥体的方位发生了改变了.

上面让我们对物体的方向和方位的区别有了一个大体上的了解,那么我们在空间中如何描述一个方位呢?这就需要使用到 角位移 了.

我们先说一个类似的例子,我们该如何描述空间中一个物体的 位置 呢?必须要把物体放在特定的坐标系中(好像很生涩).比如,如果我们说在一个坐标系中,有一个点是[1,1,1],那么你会非常轻易的想到了这个点在空间中的位置.描述 空间位置 其实就是描述相对于给定参考点(坐标原点)的 位移 .

其实,描述一个物体的方位是一样的,我们是不可能凭空描述一个物体的方位,我盟需要一个已知方位的参考量,通过这个参考量的旋转得到当前方位,那么旋转的量就叫做 角位移 .通过概念我们知道,角位移就是用来描述方位的,类似于速度就是用来描述物体运动快慢的一样.当然了,这里我要声明的一点就是虽然角位移是用来描述方位的,但是两者是不同的.例如,我们可以这么说,一个物体的方位是如何如何的;一个物体是通过某个已知方位经过角位移XXX旋转得到.所以说,方位是用来描述一个单一的"状态".但是角位移是用来描述两个状态之间的差异.

那么,我们在实际中如何描述方位与角位移呢?具体而言,我们使用矩阵和四元数来表示"角位移",用欧拉角来表示方位.接下来,我们逐一介绍一下.

在3D环境中,描述坐标系中方位的方式就是列出这个坐标系的基向量,当然了,这些基向量是用其他表示的,并不是它本身的基向量,比如当前转换完成的坐标系的三个基向量 p [1,0,0] q[0,1,0] r[0,0,1] ,这是使用本身的坐标系表示,如果放在其他坐标系中表示当前的三个基向量可能就会发生改变.这是因为参照点选择的不同.至于基向量是如何改变的就需要在 3D图形:矩阵与线性变换 说过的旋转矩阵的相关知识了.这个就不过多的解释了.比如下图,由向量 p,q,r 组建的新的坐标系用原来的坐标系表示确实如图右边所示.

其实对于我们开发来说,我们只需要知道方位是可以使用3X3矩阵来表示的.矩阵表示的是转换后的基向量即可.接下来我们说一下使用矩阵来表示角位移有什么样的优势和缺点.我就直接拿书上所讲的了,各位看官莫怪莫怪.

当然了,我们使用矩阵来表示角位移只是作为了解而已,接下来,我们看一下如何使用欧拉角表示方位的.

很多人在大学中可能会接触到矩阵,但是欧拉角可能是接触的比较少,最少作为一个学物理的我是这样的.一开始觉得欧拉角比较难理解,但是看了3D图形之后,发现用欧拉角表示方位将会比矩阵更加的直观而且易于使用.下面我们就看一下欧拉角相关的知识.(下面的基本概念跟书上的差不多,因为我觉得书上写个就很好了,所以我就没有再次总结,所以只是写了一遍.)

首先,欧拉角的基本思想是 将角位移分解为绕三个互相垂直轴的三个旋转组成的序列 .那么这个三个互相垂直的轴是如何定义的呢?其实任意三个轴和任意顺序都是可以的,但是最常用的就是使用笛卡尔坐标系并且按照一定顺序组成的旋转序列.最常用的约定,就是所谓的 "heading-pitch-bank"约定 ,在这个系统中,一个方位被定义为heading角,一个pitch角,一个bank角.其中,在左手坐标系中,我们把heading角定义为绕y轴旋转量,pitch角为绕x轴旋转量,bank角为绕z轴旋转量.旋转法则遵守左手法则(具体请参考 3D图形:矩阵与线性变换 中的旋转模块).它的基本思想是让物体开始于"标准"方位,就是物体坐标轴和惯性坐标轴对齐.让物体做heading、pitch、bank旋转之后达到最终的空间方位.

例如下图一个锥体,一开始它自身坐标轴与惯性坐标轴是一致.

然后我把heading角设置为45°.根据左手法则(通常使用,但是决定每个旋转的正方向不一定要准守右手或者左手定则),它是会做顺时针旋转.

接着物体的坐标系就发生如下的改变了.锥体的自身坐标轴不再与惯性坐标轴一致,x,z轴都发生了对应的改变.当然了,物体的空间方位也发生了对应的改变.

然后接下来就是pitch、bank旋转,分别是绕x轴旋转和z轴旋转,跟heading旋转是类似的,最后得到锥体的最终的空间方位.这里需要注意的是 不管是 heading旋转、 pitch旋转还是bank旋转,旋转的坐标轴都是自身的坐标轴!不是惯性坐标轴!

上面,看完了"heading-pitch-bank"约定系统是如何做空间方位的旋转改变的,接下来,我们来瞅瞅关于欧拉角的其他约定.

上面我们对欧拉角的接下来,我们看一下欧拉角的优点和缺点.透露一点,其实欧拉角的缺点就是引起万向锁的原因.

其实是使用欧拉角会出现一个非常有趣的现象,那就是 万向锁 ,我们看一下"heading-pitch-bank"系统这个系统中,如果pitch角度为±90°,那么就出事了,会出现什么问题呢?heading角与bank角如果相同,那么你会发现物体最终的方位是一致的,这怎么可能,这就比较尴尬了,其实类似于这种旋转pitch角度为±90°中,物体是缺失一个旋转轴的.也就是说,当pitch角度为±90°,那么bank是0.只有heading一个旋转轴起作用,是不是懵圈了?没问题,下面我要分享一个视频,我觉得这个视频会比文字更加生动形象,请对照上面的文字自行研究.

看完使用矩阵和欧拉角表示方位.接下来,我们就看一下四元数, 四元数 一个新的概念出现在我的眼前的时候我在想,他否是因为有四个数才叫四元数,确实,四元数实际是一个标量分量和一个3D向量分量组成用来表示方位.四元数的两种记法如下所示:[ω, ν ],[ω,(x,y,z)].
复数,真心好久没用了.高中的时候我们就开始接触简单的复数了,现在简单说一下复数,其实我也顺道复习一下了.
首先,复数的形式为a+bi,其中i²=-1,a称作实部(实数部分),b称作虚部(虚数部分).对于复数的运算,我们主要说说 复数的模 , 复数的模 可以很好的表示2D中的旋转变换,我们先看看前面说到过的2D环境中的旋转矩阵.

然后,我们再看一下,一个示例,假设一个复数v = (x,y)旋转θ度得到v\',如下图所示.

为了完成此次的旋转,我们需要引入第二个复数 q = (cosθ,sinθ),现在旋转之后的复数v\'就可以使用复数的乘法计算出来了.计算过程如下所示.
v = x +yi
q = cosθ +isinθ
v\' = vq = (x +yi)(cosθ +isinθ) = (xcosθ-ysinθ)+(xsinθ+ycosθ)i
跟上面的2D环境中旋转矩阵效果是一样的.只是形式不相同而已.

上面说了这么一大堆,那么到底四元数和复数有着怎样的关系呢?其实一个四元数[w,(x,y,z)]定义了复数 w +xi +yj +zk ,也就是说一个四元数是包含着一个实部和三个虚部.
其实四元数的出现也是有故事的,我直接把书上搬过来,当做在枯燥的学习中的一个轻松时刻吧(实际上,然并卵😂😂😂),爱尔兰的数学家哈密尔顿其实一直想把复数复数从2D扩展到3D,一开始他认为,3D中的复数应该有一个实部和两个虚部,然后他没有创造出这种一个实部两个虚部有意义的复数.1843年,在他去演讲的路上他突然意识到应该有三个虚部而不是两个虚部.他把这种新复数类型行者的等式刻在了Broome桥上.这样四元数就诞生了.等式如下所示.
i² = j² = k² = -1
ij = k,ji = -k
jk = i,kj = -i
ki = j,ik = -j

我们已经知道了矩阵和欧拉角的情况,现在我们就看一下四元数是如何表示角位移的.在3D环境中任意的一个角位移都可以理解为绕某个轴旋转一定的角度,在 3D图形:矩阵与线性变换 这个里面曾经说过一个3D中绕任意轴旋转的公式(还记得当初那个验证过程吗,愣是搞了一天😭,具体验证过程就不说了,请查看原来的文章).公式如下所示.其中,θ代表着旋转角度, n 代表着旋转轴.因此轴-角对( n ,θ)定义了一个角位移:绕 n 指定的轴旋转θ角.

四元数的解释其实就是角位移的轴-角对方式,但是呢, n 和θ并不是直接放入到四元数中的.它们的形式如下所示.

那么问题来了,为什么不直接放入四元数中呢?这是有原因的,这个原因,我将会在下一篇四元数的相关运算中来说明一下.现在只要知道四元数的解释其实就是角位移的轴-角对方式即可.

</b>

自己写完这篇文章总算是对矩阵、欧拉角、四元数、角位移、方位有了一个大体的了解了.整体下来发现真心枯燥的,但是还是坚持了下来了,希望小伙伴也能坚持看完,不懂的或者有疑问可以与骚栋一起探讨.3D图像下一篇我将接着研究本篇的四元数,不过是与四元数的运算相关的知识.希望大家持续关注.

最后还是要附上<<3D数学基础 图形与游戏开发>>的pdf版的传送门.

3D图形学在游戏开发中的，矩阵，四元数，欧拉角之间的底层转换算法。

在游戏开发的过程中难免会遇到欧拉角和四元数直接的转换问题，如果有些过shader的朋友，肯定也遇到过四元数，欧拉角和矩阵直接的转换问题，这里我把这几种格式直接的转换算法写在这里有需要的朋友可以拿去有，别忘了，点赞关注。废话不多说，直接上代码、

四元数转矩阵的底层算法：

public Quaternion QuaternionMatrix(float w, float x, float y, float z)
    {
        Matrix4x4 matrix = new Matrix4x4();

        matrix.m00 = 1f - 2 * SetSquare(y) - 2 * SetSquare(z);
        matrix.m01 = 2f * (x * y) - 2f * (w * z);
        matrix.m02 = 2f * (x * z) + 2f * (w * y);
        matrix.m03 = 0.0f;

        matrix.m10 = 2f * (x * y) + 2f * (w * z);
        matrix.m11 = 1f - 2f * SetSquare(x) - 2f * SetSquare(z);
        matrix.m12 = 2f * (y * z) - 2f * (w * x);
        matrix.m13 = 0.0f;

        matrix.m20 = 2f * (x * z) - 2f * (w * y);
        matrix.m21 = 2f * (y * z) + 2f * (w * x);
        matrix.m22 = 1f - 2f * SetSquare(x) - 2f * SetSquare(y);
        matrix.m23 = 0.0f;

        matrix.m30 = 0.0f;
        matrix.m31 = 0.0f;
        matrix.m32 = 0.0f;
        matrix.m33 = 0.0f;

        float qw = Mathf.Sqrt(1f + matrix.m00 + matrix.m11 + matrix.m22) / 2;
        float wq = 4 * qw;
        float qx = (matrix.m21 - matrix.m12) / wq;
        float qy = (matrix.m02 - matrix.m20) / wq;
        float qz = (matrix.m10 - matrix.m01) / wq;
        return new Quaternion(qx, qy, qz, qw);
    }

矩阵转四元数的底层算法：

 public Quaternion MatrixToQuaternion(Matrix4x4 matrix)
    {
        float qw = Mathf.Sqrt(1f + matrix.m00 + matrix.m11 + matrix.m22) / 2;
        float w = 4 * qw;
        float qx = (matrix.m21 - matrix.m12) / w;
        float qy = (matrix.m02 - matrix.m20) / w;
        float qz = (matrix.m10 - matrix.m01) / w;
        return new Quaternion(qx, qy, qz, qw);
    }

四元数转欧拉角的底层算法实现：
这里的四元数转欧拉角我要特别做一下说明，这里我有四种实现方式，其中有三种是解决个别角度问题的

可以直接拿去用的算法

public Vector3 QauToE4(float x_, float y_, float z_, float w_)
    {

>         float check = 2.0f * (-y_ * z_ + w_ * x_);****
        if (check < -0.995f)
        {
            return new Vector3(
                -90.0f,
                0.0f,
                -Mathf.Atan2(2.0f * (x_ * z_ - w_ * y_), 1.0f - 2.0f * (y_ * y_ + z_ * z_)) * M_RADTODEG
            );
        }
        else if (check > 0.995f)
        {
            return new Vector3(
                90.0f,
                0.0f,
                Mathf.Atan2(2.0f * (x_ * z_ - w_ * y_), 1.0f - 2.0f * (y_ * y_ + z_ * z_)) * M_RADTODEG
            );
        }
        else
        {
            return new Vector3(
              Mathf.Asin(check) * M_RADTODEG,
                Mathf.Atan2(2.0f * (x_ * z_ + w_ * y_), 1.0f - 2.0f * (x_ * x_ + y_ * y_)) * M_RADTODEG,
                Mathf.Atan2(2.0f * (x_ * y_ + w_ * z_), 1.0f - 2.0f * (x_ * x_ + z_ * z_)) * M_RADTODEG
            );
        }
    }

解决个别角度问题的算法

public Vector3 QuaToE(float x, float y, float z, float w)
    {
        float h, p, b;

        float sp = -2.0f * (y * z + w * x);
        if (Mathf.Abs(sp) > 0.9999f)
        {
            p = 1.570796f * sp;

            h = Mathf.Atan2(-x * z - w * y, 0.5f - y * y - z * z);
            b = 0.0f;
        }
        else
        {
            p = Mathf.Asin(sp);
            h = Mathf.Atan2(x * z - w * y, 0.5f - x * x - y * y);
            b = Mathf.Atan2(x * y - w * z, 0.5f - x * x - z * z);
        }
        return new Vector3(h, p, b);
    }

    public Vector3 QuaToE2(float x, float y, float z, float w)
    {
        float h, p, b;
        float sp = -2.0f * (y * z - w * x);

        if (Mathf.Abs(sp) > 0.9999f)
        {
            p = 1.570796f * sp;
            h = Mathf.Atan2(-x * z + w * y, 0.5f - y * y - z * z);
            b = 0.0f;
        }
        else
        {
            p = Mathf.Asin(sp);
            h = Mathf.Atan2(x * z + w * y, 0.5f - x * x - y * y);
            b = Mathf.Atan2(x * y + w * z, 0.5f - x * x - z * z);
        }
        return new Vector3(h, p, b);
    }
    public Vector3 QauToE3(float x, float y, float z, float w)
    {

        Vector3 euler;
        const float Epsilon = 0.0009765625f;
        const float Threshold = 0.5f - Epsilon;

        float TEST = w * y - x * z;

        if (TEST < -Threshold || TEST > Threshold) // 奇异姿态,俯仰角为±90°
        {
            float sign = Mathf.Sign(TEST);

            euler.z = -2 * sign * (float)Mathf.Atan2(x, w); // yaw

            euler.y = sign * (float)(3.1415926f / 2.0); // pitch

            euler.x = 0; // roll

        }
        else
        {
            euler.x = Mathf.Atan2(2 * (y * z + w * x), w * w - x * x - y * y + z * z);
            euler.y = Mathf.Asin(-2 * (x * z - w * y));
            euler.z = Mathf.Atan2(2 * (x * y + w * z), w * w + x * x - y * y - z * z);
        }
        return euler;
    }

欧拉角转四元数算法：

public Quaternion E4ToQua(float x, float y, float z)
    {
        float w_, x_, y_, z_;
        x *= M_DEGTORAD_2;
        y *= M_DEGTORAD_2;
        z *= M_DEGTORAD_2;
        float sinX = Mathf.Sin(x);
        float cosX = Mathf. Cos(x);
        float sinY = Mathf.Sin(y);
        float cosY = Mathf.Cos(y);
        float sinZ = Mathf.Sin(z);
        float cosZ = Mathf.Cos(z);

        w_ = cosY * cosX * cosZ + sinY * sinX * sinZ;
        x_ = cosY * sinX * cosZ + sinY * cosX * sinZ;
        y_ = sinY * cosX * cosZ - cosY * sinX * sinZ;
        z_ = cosY * cosX * sinZ - sinY * sinX * cosZ;
        return new Quaternion(x_, y_, z_, w_);
    }

算法测试：
这里使用的unity进行的测试，不过这里提供的算法是通用的。