C++中的四元数（Quaternion）与复数（Complex Number）运算：在图形学与物理模拟中的优化 - 智猿学院-前后端，数据库，人工智能，云计算等领域前沿技术讲座

C++中的四元数与复数运算：在图形学与物理模拟中的优化

大家好，今天我们来深入探讨C++中四元数和复数的运算，以及它们在图形学和物理模拟中的优化应用。四元数和复数都是强大的数学工具，理解它们的底层机制和优化策略，对于开发高性能的图形和物理引擎至关重要。

1. 复数基础与C++实现

复数，顾名思义，由实部和虚部组成，通常表示为 a + bi，其中 a 和 b 是实数，i 是虚数单位，满足 i² = -1。在C++中，我们可以使用结构体或类来表示复数。

#include <iostream>
#include <cmath> // for sqrt()

struct Complex {
    double real;
    double imag;

    Complex(double r = 0.0, double i = 0.0) : real(r), imag(i) {}

    // 加法
    Complex operator+(const Complex& other) const {
        return Complex(real + other.real, imag + other.imag);
    }

    // 减法
    Complex operator-(const Complex& other) const {
        return Complex(real - other.real, imag - other.imag);
    }

    // 乘法
    Complex operator*(const Complex& other) const {
        return Complex(real * other.real - imag * other.imag, real * other.imag + imag * other.real);
    }

    // 除法
    Complex operator/(const Complex& other) const {
        double denominator = other.real * other.real + other.imag * other.imag;
        return Complex((real * other.real + imag * other.imag) / denominator, (imag * other.real - real * other.imag) / denominator);
    }

    // 共轭复数
    Complex conjugate() const {
        return Complex(real, -imag);
    }

    // 模长
    double magnitude() const {
        return std::sqrt(real * real + imag * imag);
    }

    // 辐角
    double argument() const {
        return std::atan2(imag, real);
    }

    friend std::ostream& operator<<(std::ostream& os, const Complex& c) {
        os << c.real << " + " << c.imag << "i";
        return os;
    }
};

int main() {
    Complex a(2.0, 3.0);
    Complex b(1.0, -1.0);

    std::cout << "a = " << a << std::endl;
    std::cout << "b = " << b << std::endl;
    std::cout << "a + b = " << a + b << std::endl;
    std::cout << "a * b = " << a * b << std::endl;
    std::cout << "a / b = " << a / b << std::endl;
    std::cout << "|a| = " << a.magnitude() << std::endl;

    return 0;
}

这个简单的 Complex 类实现了基本的复数运算，包括加法、减法、乘法、除法、共轭复数、模长和辐角。 operator 重载使得代码更易读，更接近数学符号。

2. 四元数基础与C++实现

四元数是复数的扩展，由一个实部和三个虚部组成，通常表示为 q = w + xi + yj + zk，其中 w, x, y, z 是实数，i, j, k 是虚数单位，满足以下关系：

i² = j² = k² = -1
ij = k, ji = -k
jk = i, kj = -i
ki = j, ik = -j

四元数可以用来表示旋转，避免了欧拉角的万向锁问题。

#include <iostream>
#include <cmath>

struct Quaternion {
    double w, x, y, z;

    Quaternion(double w = 1.0, double x = 0.0, double y = 0.0, double z = 0.0) : w(w), x(x), y(y), z(z) {}

    // 构造函数：从轴角创建四元数
    Quaternion(double angle, double x, double y, double z) {
        double halfAngle = angle * 0.5;
        double sinHalfAngle = std::sin(halfAngle);
        w = std::cos(halfAngle);
        this->x = x * sinHalfAngle;
        this->y = y * sinHalfAngle;
        this->z = z * sinHalfAngle;
        normalize();
    }

    // 加法
    Quaternion operator+(const Quaternion& other) const {
        return Quaternion(w + other.w, x + other.x, y + other.y, z + other.z);
    }

    // 减法
    Quaternion operator-(const Quaternion& other) const {
        return Quaternion(w - other.w, x - other.x, y - other.y, z - other.z);
    }

    // 乘法 (四元数乘法)
    Quaternion operator*(const Quaternion& other) const {
        return Quaternion(
            w * other.w - x * other.x - y * other.y - z * other.z,
            w * other.x + x * other.w + y * other.z - z * other.y,
            w * other.y - x * other.z + y * other.w + z * other.x,
            w * other.z + x * other.y - y * other.x + z * other.w
        );
    }

    // 共轭
    Quaternion conjugate() const {
        return Quaternion(w, -x, -y, -z);
    }

    // 模长
    double magnitude() const {
        return std::sqrt(w * w + x * x + y * y + z * z);
    }

    // 归一化
    void normalize() {
        double mag = magnitude();
        if (mag > 0.0) {
            w /= mag;
            x /= mag;
            y /= mag;
            z /= mag;
        } else {
            // Handle the case where magnitude is zero (avoid division by zero)
            w = 1.0;
            x = 0.0;
            y = 0.0;
            z = 0.0;
        }
    }

    // 旋转向量
    Vector3 rotate(const Vector3& v) const {
        Quaternion p(0.0, v.x, v.y, v.z);
        Quaternion q_conj = conjugate();
        Quaternion rotated_q = (*this) * p * q_conj;
        return Vector3(rotated_q.x, rotated_q.y, rotated_q.z);
    }

    friend std::ostream& operator<<(std::ostream& os, const Quaternion& q) {
        os << q.w << " + " << q.x << "i + " << q.y << "j + " << q.z << "k";
        return os;
    }
};

// 简单向量结构体
struct Vector3 {
    double x, y, z;
    Vector3(double x = 0.0, double y = 0.0, double z = 0.0) : x(x), y(y), z(z) {}
};

int main() {
    // 创建一个绕Y轴旋转90度的四元数
    Quaternion q(M_PI / 2.0, 0.0, 1.0, 0.0);

    // 定义一个向量
    Vector3 v(1.0, 0.0, 0.0);

    // 旋转向量
    Vector3 rotated_v = q.rotate(v);

    std::cout << "Original vector: (" << v.x << ", " << v.y << ", " << v.z << ")" << std::endl;
    std::cout << "Rotated vector: (" << rotated_v.x << ", " << rotated_v.y << ", " << rotated_v.z << ")" << std::endl;

    return 0;
}

这个 Quaternion 类包含了加法、减法、乘法（四元数乘法）、共轭、模长、归一化和旋转向量的功能。注意四元数乘法的顺序很重要，因为它是非交换的。 normalize() 函数用于确保四元数表示的是一个旋转，避免累积误差。 rotate() 函数使用四元数来旋转一个三维向量。

3. 图形学中的应用

旋转表示： 四元数是表示旋转的理想选择，避免了欧拉角的万向锁问题。在3D游戏中，可以使用四元数来表示物体的朝向。
平滑插值： 可以使用球面线性插值（Slerp）在两个四元数之间进行平滑插值，创建平滑的旋转动画。
相机控制： 四元数可以用来控制相机的旋转，提供更自然的相机控制方式。

// 球面线性插值 (Slerp)
Quaternion slerp(const Quaternion& q1, const Quaternion& q2, double t) {
    double dot = q1.w * q2.w + q1.x * q2.x + q1.y * q2.y + q1.z * q2.z;

    // Ensure dot product is positive to take the shorter path
    if (dot < 0.0) {
        return slerp(Quaternion(-q2.w, -q2.x, -q2.y, -q2.z), q1, t); // Invert q2
    }

    if (std::abs(dot) > 0.9995) {  // Vectors are almost parallel
        // Linear interpolation is sufficient
        return Quaternion(
            q1.w + t * (q2.w - q1.w),
            q1.x + t * (q2.x - q1.x),
            q1.y + t * (q2.y - q1.y),
            q1.z + t * (q2.z - q1.z)
        ).normalize();
    }

    double angle = std::acos(dot);
    double sinAngle = std::sin(angle);

    double w1 = std::sin((1.0 - t) * angle) / sinAngle;
    double w2 = std::sin(t * angle) / sinAngle;

    return Quaternion(
        w1 * q1.w + w2 * q2.w,
        w1 * q1.x + w2 * q2.x,
        w1 * q1.y + w2 * q2.y,
        w1 * q1.z + w2 * q2.z
    ).normalize();
}

int main() {
    Quaternion q1(M_PI / 4.0, 0.0, 1.0, 0.0); // 绕Y轴旋转45度
    Quaternion q2(M_PI / 2.0, 0.0, 1.0, 0.0); // 绕Y轴旋转90度

    // 在 q1 和 q2 之间插值
    Quaternion q_interp = slerp(q1, q2, 0.5); // 插值系数为 0.5

    std::cout << "q1: " << q1 << std::endl;
    std::cout << "q2: " << q2 << std::endl;
    std::cout << "Interpolated quaternion: " << q_interp << std::endl;

    return 0;
}

slerp() 函数实现了球面线性插值。首先计算两个四元数的点积，用于确定插值路径。如果点积接近 1，则可以使用线性插值，避免计算 acos 和 sin 函数。如果点积是负的，则需要反转其中一个四元数，以选择较短的插值路径。

4. 物理模拟中的应用

刚体旋转： 四元数可以用来表示刚体的旋转，避免了欧拉角的万向锁问题。
角速度更新： 可以使用四元数来更新刚体的角速度，实现更真实的旋转运动。
碰撞检测： 四元数可以用来计算碰撞后的旋转，模拟真实的物理效果。

// 假设我们有一个表示刚体状态的结构体
struct RigidBodyState {
    Vector3 position;
    Quaternion orientation; // 使用四元数表示朝向
    Vector3 linearVelocity;
    Vector3 angularVelocity; // 世界坐标系下的角速度
};

// 根据角速度更新四元数朝向
void updateOrientation(RigidBodyState& state, double deltaTime) {
    // 将世界坐标系下的角速度转换为四元数形式
    Quaternion angularVelocityQuaternion(0.0, state.angularVelocity.x, state.angularVelocity.y, state.angularVelocity.z);

    // 计算四元数导数
    Quaternion derivative = angularVelocityQuaternion * state.orientation * 0.5;

    // 更新四元数
    state.orientation.w += derivative.w * deltaTime;
    state.orientation.x += derivative.x * deltaTime;
    state.orientation.y += derivative.y * deltaTime;
    state.orientation.z += derivative.z * deltaTime;

    // 归一化四元数
    state.orientation.normalize();
}

int main() {
    RigidBodyState body;
    body.position = Vector3(0.0, 0.0, 0.0);
    body.orientation = Quaternion(1.0, 0.0, 0.0, 0.0); // 初始朝向：没有旋转
    body.linearVelocity = Vector3(0.0, 0.0, 0.0);
    body.angularVelocity = Vector3(0.0, 1.0, 0.0); // 绕Y轴旋转的角速度

    double deltaTime = 0.01; // 时间步长

    for (int i = 0; i < 100; ++i) {
        updateOrientation(body, deltaTime);
        std::cout << "Orientation at step " << i << ": " << body.orientation << std::endl;
    }

    return 0;
}

updateOrientation() 函数展示了如何使用角速度来更新四元数朝向。首先将角速度转换为四元数形式，然后计算四元数导数，最后更新四元数。每次更新后，都需要对四元数进行归一化。

5. 性能优化

内联函数： 将频繁调用的函数声明为 inline，可以减少函数调用的开销。
SIMD指令： 使用SIMD（单指令多数据）指令可以并行执行多个浮点运算，提高计算速度。例如，可以使用Intel的SSE或AVX指令集。
数据对齐： 确保四元数和向量的数据在内存中对齐，可以提高SIMD指令的效率。
避免不必要的计算： 例如，如果只需要比较两个四元数的朝向是否相同，可以使用点积而不是计算完整的旋转矩阵。
预计算： 如果某些计算结果可以重复使用，可以将其预先计算并存储起来，避免重复计算。例如，可以预先计算旋转矩阵，避免每次都从四元数计算。
使用库： 考虑使用成熟的数学库，如Eigen或GLM，它们已经对四元数和复数运算进行了优化。

#include <immintrin.h> // 包含 AVX 指令集

// 使用 AVX 指令优化四元数乘法 (假设编译器支持 AVX)
Quaternion multiply_avx(const Quaternion& q1, const Quaternion& q2) {
    // 将四元数加载到 AVX 寄存器中
    __m256d q1_vec = _mm256_set_pd(q1.w, q1.z, q1.y, q1.x);
    __m256d q2_vec = _mm256_set_pd(q2.w, q2.z, q2.y, q2.x);

    // 执行 AVX 指令进行乘法运算 (这只是一个简化的例子，需要更复杂的逻辑来实现完整的四元数乘法)
    __m256d result_vec = _mm256_mul_pd(q1_vec, q2_vec);

    // 从 AVX 寄存器中提取结果
    double result[4];
    _mm256_storeu_pd(result, result_vec);

    return Quaternion(result[3], result[2], result[1], result[0]); // 注意顺序
}

上面的代码展示了如何使用AVX指令来优化四元数乘法。这只是一个简化的例子，完整的四元数乘法需要更复杂的逻辑。使用SIMD指令可以显著提高计算速度，尤其是在处理大量数据时。需要注意的是，SIMD指令的使用需要编译器和硬件的支持。

6. 选择合适的表示方法

在某些情况下，复数可以用来简化四元数运算。例如，如果只需要表示绕一个轴的旋转，可以使用复数来表示旋转角度，然后将其转换为四元数。此外，也可以使用矩阵来表示旋转，但在某些情况下，四元数更紧凑、更高效。

7. 调试与测试

单元测试： 编写单元测试来验证四元数和复数运算的正确性。
可视化： 使用可视化工具来观察旋转效果，例如，可以使用OpenGL或DirectX来渲染旋转后的物体。
性能分析： 使用性能分析工具来确定代码中的瓶颈，并进行优化。

8. 案例分析：优化Slerp函数

Slerp函数是图形学中常用的函数，其性能至关重要。以下是一些优化Slerp函数的技巧：

提前退出: 如果两个四元数非常接近（点积接近1），使用线性插值代替Slerp，避免昂贵的 acos 和 sin 计算。
避免负点积: 确保两个四元数的点积为正，选择较短的旋转路径，避免不必要的旋转。
查表法: 预先计算一些 sin 和 cos 值，存储在表中，在需要时直接查表，避免重复计算。
SIMD优化: 使用SIMD指令并行计算多个Slerp值。

9. 实际代码示例：使用Eigen库进行四元数运算

Eigen是一个强大的C++线性代数库，提供了高效的四元数实现。以下是一个使用Eigen库进行四元数运算的示例：

#include <iostream>
#include <Eigen/Dense>
#include <Eigen/Geometry>

int main() {
    // 创建一个Eigen四元数
    Eigen::Quaterniond q1(Eigen::AngleAxisd(M_PI / 4, Eigen::Vector3d::UnitY())); // 绕Y轴旋转45度

    // 创建一个Eigen向量
    Eigen::Vector3d v(1.0, 0.0, 0.0);

    // 旋转向量
    Eigen::Vector3d rotated_v = q1 * v;

    std::cout << "Original vector: (" << v.x() << ", " << v.y() << ", " << v.z() << ")" << std::endl;
    std::cout << "Rotated vector: (" << rotated_v.x() << ", " << rotated_v.y() << ", " << rotated_v.z() << ")" << std::endl;

    // 使用Eigen进行Slerp插值
    Eigen::Quaterniond q2(Eigen::AngleAxisd(M_PI / 2, Eigen::Vector3d::UnitY())); // 绕Y轴旋转90度
    Eigen::Quaterniond q_interp = q1.slerp(0.5, q2); // 插值系数为 0.5

    std::cout << "Interpolated quaternion: (" << q_interp.w() << ", " << q_interp.x() << ", " << q_interp.y() << ", " << q_interp.z() << ")" << std::endl;

    return 0;
}

Eigen库提供了方便的API和优化的实现，可以大大简化四元数运算的代码。

四元数与复数：理解与应用是关键

掌握四元数和复数的基础知识，并理解它们在图形学和物理模拟中的应用场景，是进行优化的前提。选择合适的表示方法，并利用各种优化技巧，可以显著提高代码的性能。

优化方向：算法、数据、硬件相结合

性能优化是一个多方面的过程，需要从算法、数据结构和硬件三个方面入手。通过选择合适的算法、优化数据结构和利用硬件特性，可以最大限度地提高代码的性能。

更多IT精英技术系列讲座，到智猿学院

C++中的四元数与复数运算：在图形学与物理模拟中的优化

发表回复 取消回复

发表回复取消回复