4x4矩阵教程

4x4矩阵教程

1. 简介

四维矩阵是计算机图形学和3D变换中的重要工具，用于表示三维空间中的仿射变换。本教程将介绍如何使用C++实现四维矩阵的基本运算和变换。

2. 代码实现

2.1 头文件 (matrix4x4.h)

#ifndef MATRIX4X4_H
#define MATRIX4X4_H

#include <array>
#include <stdexcept>
#include <iostream>

namespace math {
namespace linear_algebra {

/**
 * @brief 四维矩阵类
 * 
 * 这个类实现了四维矩阵的基本运算，包括：
 * - 矩阵加减
 * - 矩阵乘法
 * - 标量乘法
 * - 行列式计算
 * - 矩阵求逆
 * - 矩阵转置
 * - 特征值和特征向量计算
 * - 矩阵性质检查（可逆性、对称性、正交性）
 * - 特殊矩阵生成（旋转矩阵、缩放矩阵、平移矩阵、投影矩阵等）
 */
class Matrix4x4 {
public:
    // 构造函数
    Matrix4x4();  // 默认构造函数，初始化为单位矩阵
    Matrix4x4(const std::array<std::array<double, 4>, 4>& data);  // 从二维数组初始化

    // 基本运算
    Matrix4x4 operator+(const Matrix4x4& other) const;  // 矩阵加法
    Matrix4x4 operator-(const Matrix4x4& other) const;  // 矩阵减法
    Matrix4x4 operator*(const Matrix4x4& other) const;  // 矩阵乘法
    Matrix4x4 operator*(double scalar) const;           // 标量乘法
    Matrix4x4 operator/(double scalar) const;           // 标量除法

    // 矩阵运算
    double determinant() const;                         // 计算行列式
    Matrix4x4 inverse() const;                          // 计算逆矩阵
    Matrix4x4 transpose() const;                        // 计算转置矩阵
    std::array<double, 4> eigenvalues() const;          // 计算特征值
    std::array<Matrix4x4, 4> eigenvectors() const;      // 计算特征向量

    // 矩阵性质
    bool isInvertible() const;                          // 检查是否可逆
    bool isSymmetric() const;                           // 检查是否对称
    bool isOrthogonal() const;                          // 检查是否正交

    // 特殊矩阵
    static Matrix4x4 identity();                        // 创建单位矩阵
    static Matrix4x4 rotation(double theta, char axis); // 创建旋转矩阵
    static Matrix4x4 scaling(double sx, double sy, double sz); // 创建缩放矩阵
    static Matrix4x4 translation(double tx, double ty, double tz); // 创建平移矩阵
    static Matrix4x4 perspective(double fov, double aspect, double near, double far); // 创建透视投影矩阵
    static Matrix4x4 orthographic(double left, double right, double bottom, double top, double near, double far); // 创建正交投影矩阵

    // 输出运算符
    friend std::ostream& operator<<(std::ostream& os, const Matrix4x4& m);

private:
    std::array<std::array<double, 4>, 4> data;  // 4x4矩阵数据
};

} // namespace linear_algebra
} // namespace math

#endif // MATRIX4X4_H

2.2 实现文件 (matrix4x4.cpp)

#include "matrix4x4.h"
#include <cmath>

namespace math {
namespace linear_algebra {

// 默认构造函数：初始化为单位矩阵
Matrix4x4::Matrix4x4() {
    for (int i = 0; i < 4; ++i) {
        for (int j = 0; j < 4; ++j) {
            data[i][j] = (i == j) ? 1.0 : 0.0;
        }
    }
}

// 从二维数组初始化
Matrix4x4::Matrix4x4(const std::array<std::array<double, 4>, 4>& data) : data(data) {}

// 矩阵加法实现
Matrix4x4 Matrix4x4::operator+(const Matrix4x4& other) const {
    Matrix4x4 result;
    for (int i = 0; i < 4; ++i) {
        for (int j = 0; j < 4; ++j) {
            result.data[i][j] = data[i][j] + other.data[i][j];
        }
    }
    return result;
}

// 矩阵减法实现
Matrix4x4 Matrix4x4::operator-(const Matrix4x4& other) const {
    Matrix4x4 result;
    for (int i = 0; i < 4; ++i) {
        for (int j = 0; j < 4; ++j) {
            result.data[i][j] = data[i][j] - other.data[i][j];
        }
    }
    return result;
}

// 矩阵乘法实现
Matrix4x4 Matrix4x4::operator*(const Matrix4x4& other) const {
    Matrix4x4 result;
    for (int i = 0; i < 4; ++i) {
        for (int j = 0; j < 4; ++j) {
            result.data[i][j] = 0.0;
            for (int k = 0; k < 4; ++k) {
                result.data[i][j] += data[i][k] * other.data[k][j];
            }
        }
    }
    return result;
}

// 标量乘法实现
Matrix4x4 Matrix4x4::operator*(double scalar) const {
    Matrix4x4 result;
    for (int i = 0; i < 4; ++i) {
        for (int j = 0; j < 4; ++j) {
            result.data[i][j] = data[i][j] * scalar;
        }
    }
    return result;
}

// 标量除法实现
Matrix4x4 Matrix4x4::operator/(double scalar) const {
    if (scalar == 0.0) {
        throw std::runtime_error("Division by zero");
    }
    return *this * (1.0 / scalar);
}

// 行列式计算实现
double Matrix4x4::determinant() const {
    // 使用拉普拉斯展开计算4x4矩阵的行列式
    double det = 0.0;
    for (int i = 0; i < 4; ++i) {
        double cofactor = 0.0;
        for (int j = 0; j < 4; ++j) {
            if ((i + j) % 2 == 0) {
                cofactor += data[0][j] * minor(0, j);
            } else {
                cofactor -= data[0][j] * minor(0, j);
            }
        }
        det += data[0][i] * cofactor;
    }
    return det;
}

// 逆矩阵计算实现
Matrix4x4 Matrix4x4::inverse() const {
    double det = determinant();
    if (det == 0.0) {
        throw std::runtime_error("Matrix is not invertible");
    }

    Matrix4x4 result;
    double invDet = 1.0 / det;

    // 计算伴随矩阵
    for (int i = 0; i < 4; ++i) {
        for (int j = 0; j < 4; ++j) {
            result.data[i][j] = cofactor(j, i) * invDet;
        }
    }

    return result;
}

// 转置矩阵实现
Matrix4x4 Matrix4x4::transpose() const {
    Matrix4x4 result;
    for (int i = 0; i < 4; ++i) {
        for (int j = 0; j < 4; ++j) {
            result.data[i][j] = data[j][i];
        }
    }
    return result;
}

// 特征值计算实现
std::array<double, 4> Matrix4x4::eigenvalues() const {
    // 计算特征多项式系数
    // 这里使用简化的方法，实际应用中可能需要更复杂的数值方法
    std::array<double, 4> roots;
    // ... 求解四次方程的代码 ...
    return roots;
}

// 特征向量计算实现
std::array<Matrix4x4, 4> Matrix4x4::eigenvectors() const {
    std::array<double, 4> eigenvals = eigenvalues();
    std::array<Matrix4x4, 4> eigenvecs;
    // ... 计算特征向量的代码 ...
    return eigenvecs;
}

// 可逆性检查实现
bool Matrix4x4::isInvertible() const {
    return determinant() != 0.0;
}

// 对称性检查实现
bool Matrix4x4::isSymmetric() const {
    for (int i = 0; i < 4; ++i) {
        for (int j = i + 1; j < 4; ++j) {
            if (data[i][j] != data[j][i]) {
                return false;
            }
        }
    }
    return true;
}

// 正交性检查实现
bool Matrix4x4::isOrthogonal() const {
    Matrix4x4 product = *this * transpose();
    Matrix4x4 identity = Matrix4x4::identity();
    for (int i = 0; i < 4; ++i) {
        for (int j = 0; j < 4; ++j) {
            if (std::abs(product.data[i][j] - identity.data[i][j]) > 1e-10) {
                return false;
            }
        }
    }
    return true;
}

// 单位矩阵创建实现
Matrix4x4 Matrix4x4::identity() {
    return Matrix4x4();
}

// 旋转矩阵创建实现
Matrix4x4 Matrix4x4::rotation(double theta, char axis) {
    Matrix4x4 result;
    double cos_theta = std::cos(theta);
    double sin_theta = std::sin(theta);

    switch (axis) {
        case 'x':
            result.data[1][1] = cos_theta;
            result.data[1][2] = -sin_theta;
            result.data[2][1] = sin_theta;
            result.data[2][2] = cos_theta;
            break;
        case 'y':
            result.data[0][0] = cos_theta;
            result.data[0][2] = sin_theta;
            result.data[2][0] = -sin_theta;
            result.data[2][2] = cos_theta;
            break;
        case 'z':
            result.data[0][0] = cos_theta;
            result.data[0][1] = -sin_theta;
            result.data[1][0] = sin_theta;
            result.data[1][1] = cos_theta;
            break;
        default:
            throw std::runtime_error("Invalid rotation axis");
    }

    return result;
}

// 缩放矩阵创建实现
Matrix4x4 Matrix4x4::scaling(double sx, double sy, double sz) {
    Matrix4x4 result;
    result.data[0][0] = sx;
    result.data[1][1] = sy;
    result.data[2][2] = sz;
    return result;
}

// 平移矩阵创建实现
Matrix4x4 Matrix4x4::translation(double tx, double ty, double tz) {
    Matrix4x4 result;
    result.data[0][3] = tx;
    result.data[1][3] = ty;
    result.data[2][3] = tz;
    return result;
}

// 透视投影矩阵创建实现
Matrix4x4 Matrix4x4::perspective(double fov, double aspect, double near, double far) {
    Matrix4x4 result;
    double tan_half_fov = std::tan(fov / 2.0);
    
    result.data[0][0] = 1.0 / (aspect * tan_half_fov);
    result.data[1][1] = 1.0 / tan_half_fov;
    result.data[2][2] = -(far + near) / (far - near);
    result.data[2][3] = -1.0;
    result.data[3][2] = -(2.0 * far * near) / (far - near);
    result.data[3][3] = 0.0;

    return result;
}

// 正交投影矩阵创建实现
Matrix4x4 Matrix4x4::orthographic(double left, double right, double bottom, double top, double near, double far) {
    Matrix4x4 result;
    
    result.data[0][0] = 2.0 / (right - left);
    result.data[1][1] = 2.0 / (top - bottom);
    result.data[2][2] = -2.0 / (far - near);
    result.data[0][3] = -(right + left) / (right - left);
    result.data[1][3] = -(top + bottom) / (top - bottom);
    result.data[2][3] = -(far + near) / (far - near);

    return result;
}

// 输出运算符实现
std::ostream& operator<<(std::ostream& os, const Matrix4x4& m) {
    for (int i = 0; i < 4; ++i) {
        os << "[ ";
        for (int j = 0; j < 4; ++j) {
            os << m.data[i][j] << " ";
        }
        os << "]" << std::endl;
    }
    return os;
}

} // namespace linear_algebra
} // namespace math

3. 使用示例

3.1 基本运算

#include "matrix4x4.h"
#include <iostream>

using namespace math::linear_algebra;

int main() {
    // 创建矩阵
    Matrix4x4 m1({
        {{1.0, 2.0, 3.0, 4.0},
         {5.0, 6.0, 7.0, 8.0},
         {9.0, 10.0, 11.0, 12.0},
         {13.0, 14.0, 15.0, 16.0}}
    });

    Matrix4x4 m2({
        {{16.0, 15.0, 14.0, 13.0},
         {12.0, 11.0, 10.0, 9.0},
         {8.0, 7.0, 6.0, 5.0},
         {4.0, 3.0, 2.0, 1.0}}
    });

    // 基本运算
    Matrix4x4 sum = m1 + m2;      // 矩阵加法
    Matrix4x4 diff = m1 - m2;     // 矩阵减法
    Matrix4x4 prod = m1 * m2;     // 矩阵乘法
    Matrix4x4 scaled = m1 * 2.0;  // 标量乘法

    // 输出结果
    std::cout << "m1:\n" << m1 << std::endl;
    std::cout << "m2:\n" << m2 << std::endl;
    std::cout << "m1 + m2:\n" << sum << std::endl;
    std::cout << "m1 - m2:\n" << diff << std::endl;
    std::cout << "m1 * m2:\n" << prod << std::endl;
    std::cout << "m1 * 2:\n" << scaled << std::endl;

    return 0;
}

3.2 3D变换

#include "matrix4x4.h"
#include <iostream>
#include <cmath>

using namespace math::linear_algebra;

int main() {
    // 创建变换矩阵
    Matrix4x4 rotation = Matrix4x4::rotation(M_PI / 4.0, 'z');  // 绕Z轴旋转45度
    Matrix4x4 scaling = Matrix4x4::scaling(2.0, 3.0, 4.0);      // 缩放
    Matrix4x4 translation = Matrix4x4::translation(1.0, 2.0, 3.0); // 平移

    // 组合变换（注意顺序：先缩放，再旋转，最后平移）
    Matrix4x4 transform = translation * rotation * scaling;

    // 创建投影矩阵
    Matrix4x4 perspective = Matrix4x4::perspective(M_PI / 4.0, 16.0/9.0, 0.1, 100.0);
    Matrix4x4 ortho = Matrix4x4::orthographic(-1.0, 1.0, -1.0, 1.0, 0.1, 100.0);

    // 输出结果
    std::cout << "旋转矩阵:\n" << rotation << std::endl;
    std::cout << "缩放矩阵:\n" << scaling << std::endl;
    std::cout << "平移矩阵:\n" << translation << std::endl;
    std::cout << "组合变换:\n" << transform << std::endl;
    std::cout << "透视投影:\n" << perspective << std::endl;
    std::cout << "正交投影:\n" << ortho << std::endl;

    return 0;
}

4. 编译和运行

4.1 编译

使用提供的Makefile进行编译：

make        # 编译所有目标
make test   # 运行测试
make demo   # 运行示例
make clean  # 清理编译文件

4.2 运行测试

./matrix4x4_test

4.3 运行示例

./matrix4x4_demo

5. 注意事项

1. 数值精度

   • 在比较浮点数时使用适当的误差范围
   • 例如：std::abs(det) < 1e-10 判断是否可逆

2. 异常处理

   • 除以零的情况
   • 不可逆矩阵求逆的情况
   • 无效的旋转轴
   • 投影参数的有效性检查

3. 性能考虑

   • 矩阵乘法的时间复杂度为O(n³)
   • 特征值计算可能需要迭代方法
   • 考虑使用SIMD指令优化计算
   • 投影矩阵的预计算

4. 使用建议

   • 优先使用成员函数而不是全局函数
   • 保持接口的一致性
   • 提供清晰的错误信息
   • 注意变换矩阵的组合顺序

6. 扩展阅读

1. 线性代数基础

   • 矩阵运算
   • 行列式
   • 特征值和特征向量

2. 3D变换

   • 旋转矩阵
   • 缩放矩阵
   • 平移矩阵
   • 投影矩阵

3. 应用领域

   • 计算机图形学
   • 3D游戏开发
   • 计算机视觉
   • 机器人学
   • 物理模拟

   声明

该文章为学习过程中的笔记，目的是防止自己忘记，也为了方便随时随地查阅。其中大部分内容收集于互联网。