games101——作业0 C++数学库与Eigen库的基本使用

作业0主要是关于编程环境的搭建，这个在官方文档中已经说得很详细了，本篇博客主要针对今后作业中需要大量使用的C++数学库及Eigen库的基本使用进行说明，并完成作业要求：旋转并平移一个点

---

---

数学库常用函数

C++ 标准库 std 通过 <cmath> 头文件内置了丰富的数学函数，可对各种数字进行运算。

如果要使用 C++ 标准库提供的数学函数，首先需要包含 <cmath> 头文件

#include <cmath>

下表列出了 C++ 标准库提供的常用数学函数

函数	返回值
double cos(double)	返回弧度角（double 型）的余弦
double sin(double)	返回弧度角（double 型）的正弦
double tan(double)	返回弧度角（double 型）的正切
double acos(double)	返回参数的反余弦弧度值，其范围为 $[0,\pi ]$
double asin(double)	返回参数的反正弦弧度值，其范围为 $[-\pi /2,\pi /2]$
double atan(double)	返回参数的反正切弧度制，其范围为 $[-\pi /2,\pi /2]$
double atan2(double y, double x)	返回坐标(x,y)与原点连线与x轴正方向的夹角，其范围为 $[-\pi ,\pi ]$
double log(double x)	返回参数的自然对数ln(x)
double log10(double x)	返回以10为底的对数lg(x)
double exp(double x)	返回 ex
double pow(double x, double y)	返回 xy
double hypot(double, double)	返回两个参数的平方总和的平方根
double sqrt(double)	返回参数的平方根
int abs(int)	返回整数的绝对值
double fabs(double)	返回任意一个十进制数的绝对值
double floor(double)	返回一个小于或等于传入参数的最大整数
double ceil(double)	返回一个大于或等于传入参数的最小整数

这里的与三角函数有关的函数里面的角度都是弧度制的。

有两种方式使用 $\pi$

• 使用 #define PI acos(-1) 因为 -1 对应的反余弦弧度制就是 $\pi$
• 使用 GNU C 库中的 M_PI 宏，该库定义了多个常用的数学常量，如下表

常数	定义
M_E	自然对数的底数
M_LOG2E	M_E 以 2 为底的对数
M_LOG10E	M_E 以 10 为底的对数
M_LN2	2 的自然对数
M_LN10	10 的自然对数
M_PI	圆周率$\pi$
M_PI_2	$\pi /2$
M_PI_4	$\pi /4$
M_1_PI	$\pi$的倒数($1/\pi$)
M_2_PI	$\pi$倒数的 2 倍($2/\pi$)
M_2_SQRTPI	$2/\sqrt{\pi }$
M_SQRT2	$\sqrt{2}$
M_SQRT1_2	$1/\sqrt{2}$

下面用一些实例对数学库函数进行说明：

#include <iostream>
#include <cmath>
using namespace std;
#define PI acos(-1)
int main()
{
    // 三角函数
    cout<<cos(PI/3)<<endl; //0.5
    cout<<sin(PI/3)<<endl; //0.866025
    cout<<tan(PI/3)<<endl; //1.73205

    // 反三角函数
    cout<<acos(0.5)<<endl; // PI/3=1.0472
    cout<<asin(0.866025)<<endl; // 1.0472
    cout<<atan(1.73205)<<endl; // 1.0472
    cout<<atan2(1, -1)<<endl; // 3PI/4=2.35619

    // 对数
    cout<<log(2)<<endl; // ln(2)=0.693147
    cout<<log10(2)<<endl; // lg(2)=0.30103

    // 平方根
    cout<<sqrt(13)<<endl; //13的平方根=3.60555

    // 次方
    cout<<exp(2)<<endl;  // e^2=7.38906
    cout<<pow(2,3)<<endl; // 2^3=8
    cout<<hypot(2,3)<<endl; // sqrt(2^2+3^2)=3.60555

    //绝对值
    cout<<abs(-2)<<endl; // |-2|=2
    cout<<fabs(-0.5)<<endl; // |-0.5|=0.5

    //取整
    cout<<floor(1.4)<<endl; //1
    cout<<floor(0.5)<<endl; //0
    cout<<floor(-0.5)<<endl; //-1
    cout<<ceil(1.6)<<endl; //2
    cout<<ceil(0.5)<<endl; //1
    cout<<ceil(-0.5)<<endl; //-0
    return 0;
}

---

Eigen库

Eigen 是可以用来进行线性代数、矩阵、向量操作等运算的C++库。

下面主要对 Eigen 中的矩阵(Matrix)类的使用进行介绍

在 Eigen 中，所有矩阵和向量均是 Matrix 模板类的对象。

矩阵定义与访问

矩阵模板的参数

Matrix 类有六个模板参数，其中三个有默认值，因此只要学习三个参数就足够了。

/* 强制性的三参数模板的原型 (三个参数分别表示：标量的类型，编译时的行，编译时的列) */
Matrix<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime> 

/* 用typedef定义了很多模板，例如：Matrix4f 表示 4×4 的floats 矩阵 */ 
typedef Matrix<float, 4, 4> Matrix4f;

向量(Vectors)

向量是矩阵的特殊情况，也是用矩阵定义的

typedef Matrix<float, 3, 1> Vector3f;  
typedef Matrix<int, 1, 2> RowVector2i;

特殊动态值 (special value Dynamic)

Eigen的矩阵不仅能够在编译是确定大小（fixed size），也可以在运行时确定大小，就是所说的动态矩阵（dynamic size）。小矩阵(16及以下)建议使用固定大小，大矩阵建议使用动态矩阵。原因是固定矩阵是分配在栈上，而动态矩阵分配在堆上，因此固定矩阵相比于动态矩阵有更小的运行时间。而对于大矩阵，这个运行时间的差距基本可以忽略，而在栈上分配空间会有栈溢出的风险，且向量化时动态尺寸也更合适。

typedef Matrix<double, Dynamic, Dynamic> MatrixXd;  
typedef Matrix<int, Dynamic, 1> VectorXi;  

/* 也可使用‘行’固定‘列’动态的矩阵 */
Matrix<float, 3, Dynamic>

构造函数 (Constructors)

可以使用默认的构造函数，不执行动态分配内存，也没有初始化矩阵参数

Matrix3f a;   // a是3-by-3矩阵，包含未初始化的 float[9] 数组
MatrixXf b;   // b是动态矩阵，当前大小为 0-by-0, 没有为数组的系数分配内存

/* 矩阵的第一个参数表示“行”，数组只有一个参数。根据跟定的大小分配内存，但不初始化 */
MatrixXf a(10,15);    // a 是10-by-15阵，分配了内存，没有初始化
VectorXf b(30);       // b是动态矩阵，当前大小为 30, 分配了内存，没有初始化

/* 对于给定的矩阵，传递的参数无效 */
Matrix3f a(3,3); 

/* 对于维数最大为4的向量，可以直接初始化 */
Vector2d a(5.0, 6.0);  
Vector3d b(5.0, 6.0, 7.0);  
Vector4d c(5.0, 6.0, 7.0, 8.0);

系数访问

系数都是从0开始，矩阵默认按列存储

#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;

int main()
{
    MatrixXd m(2, 2);
    m(0, 0) = 3;
    m(1, 0) = 2.5;
    m(0, 1) = -1;
    m(1, 1) = m(1, 0) + m(0, 1);
    cout << "Here is the matrix m:" << endl;
    cout << m << endl;

    VectorXd v(2);
    v(0) = 4;
    v[1] = v[0] - 1;     //operator[] 在 vectors 中重载，意义和()相同
    cout << "Here is the vector v:" << endl;
    cout << v << endl;

    getchar();
    getchar();
}

初始化

Matrix3f m;
// 逗号隔开的初始化
m << 1, 2, 3,   4, 5, 6,   7, 8, 9;
cout << m;

//初始化向量
Vector2i a(1,2); // 一个列向量包含元素(1,2)

Eigen::Matrix3d m2 = Eigen::Matrix3d::Identity();//Eigen::Matrix3d::Zero();

Eigen::Matrix3d m3 = Eigen::Matrix3d::Random(); //随机初始化

Eigen::Vector3f v1 = Eigen::Vector3f::Zero();

//1. Eigen::RowVectorXd //行向量,之前的Vector都是列向量
using namespace Eigen;
RowVectorXd rv1(1,2,3);
RowVectorXd rv2(4);
rv2 << 1,2,3,4;
RowVectorXd joined_rv12(7);
joined_rv12 << rv1 , rv2;//rv1后接rv2
//joined_rv12 is:
1 2 3 1 2 3 4;

//2. 块操作
MatrixXf m4(2,2);
m4 << 1,2,3,4;
MatrixXf m5(4,4);
m5 << m4, m4 / 10, m4 * 10, m4;//将m5分了四块赋值
//m5 is :
1  2  0.1 0.2 
3  4  0.3 0.4
10 20  1   2
30 40  3   4

Matrix3f m6;
m6.row(0) << 1,2,3;
m6.block(1,0,2,2) << 4,5,6,7; //
//m6.block<2,2>(1,0) << 4,5,6,7;
m6.col(2).tail(2) << 6,9;
//m6 is: 
1 2 3
4 5 6
6 7 9

修改尺寸

获取矩阵的行数，列数，以及系数，可以分别使用 rows()，cols()，和 size()。改变动态大小的矩阵尺寸可以使用 resize() 方法。

#include <iostream>
#include <Eigen/Dense>
 
int main()
{
  Eigen::MatrixXd m(2,5);
  m.resize(4,3);
  std::cout << "The matrix m is of size "
            << m.rows() << "x" << m.cols() << std::endl; // 4x3
  std::cout << "It has " << m.size() << " coefficients" << std::endl; // 12
  Eigen::VectorXd v(2);
  v.resize(5);
  std::cout << "The vector v is of size " << v.size() << std::endl; // 5
  std::cout << "As a matrix, v is of size "
            << v.rows() << "x" << v.cols() << std::endl; //5x1

注意，改变一个固定尺寸的矩阵大小会引起断言错误，但下面这样没有问题(因为实际并没有修改)

#include <iostream>
#include <Eigen/Dense>
 
int main()
{
  Eigen::Matrix4d m;
  m.resize(4,4); // no operation
  std::cout << "The matrix m is of size "
            << m.rows() << "x" << m.cols() << std::endl; //4x4
}

可选模板参数

之前只讨论了三个参数，下面为全部参数

Matrix<typename Scalar,
       int RowsAtCompileTime,
       int ColsAtCompileTime,
       int Options = 0,
       int MaxRowsAtCompileTime = RowsAtCompileTime,
       int MaxColsAtCompileTime = ColsAtCompileTime>

• Options：默认是0，代表是列优先，即初始化时按列开始放，(0,0)->(0,1)->(0,2)->(1,0)->…，还可以选择 Options=RowMajor，代表行优先，即初始化按行开始放，(0,0)->(1,0)->(2,0)->(0,1)->…，比如下面为 3x3 行优先矩阵：

Matrix<float, 3, 3, RowMajor>

• MaxRowsAtCompileTime 和 MaxColsAtCompileTime：分别表示矩阵最大的行和最大的列，有时你可能不知道具体要用多大的行和列，但可以确定一个上界，这样做可以避免动态内存分配，例如下面为最大 3x4 的矩阵：

Matrix<float, Dynamic, Dynamic, 0, 3, 4>

方便的typedefs

Eigen 定义了如下矩阵的 typedefs：

• MatrixNt for Matrix<type, N, N> 比如 MatrixXi for Matrix<int, Dynamic, Dynamic>
• MatrixXNt for Matrix<type, Dynamic, N> 比如 MatrixX3i for Matrix<int, Dynamic, 3>
• MatrixNXt for Matrix<type, N, Dynamic> 比如 Matrix4Xd for Matrix<d, 4, Dynamic>
• VectorNt for Matrix<type, N, 1> 比如 Vector2f for Matrix<float, 2, 1>
• RowVectorNt for Matrix<type, 1, N> 比如 RowVector3d for Matrix<double, 1, 3>

这里

• N 可以是 2，3，4，或者 X(意味着 Dynamic)
• t 可以时 i(int)，f(float)，d(double)，cf(complex<float>)，cd(complex<double>)，虽然只提供五种类型的type，但不代表Matrix只支持五种类型。

---

矩阵运算

加法和减法

进行加法和减法的两个矩阵 a 和 b 必须有同样数量的行和列，且必须有同样的数据类型，Eigen 并不会做自动的类型提升，矩阵也不能和数进行加减运算，主要包含五种操作：a+b、a-b、-a、a+=b、a-=b

#include <iostream>
#include <Eigen/Dense>
 
int main()
{
  Eigen::Matrix2d a;
  a << 1, 2,
       3, 4;
  Eigen::MatrixXd b(2,2);
  b << 2, 3,
       1, 4;
  /*
  a + b =
  3 5
  4 8
  */
  std::cout << "a + b =\n" << a + b << std::endl;
  /*
  a - b =
 -1 -1
  2  0
  */ 
  std::cout << "a - b =\n" << a - b << std::endl;
  std::cout << "Doing a += b;" << std::endl;
  a += b;
  /*
  Now a =
  3 5
  4 8
  */
  std::cout << "Now a =\n" << a << std::endl;
  Eigen::Vector3d v(1,2,3);
  Eigen::Vector3d w(1,0,0);
  /*
  -v + w - v =
  -1
  -4
  -6
  */
  std::cout << "-v + w - v =\n" << -v + w - v << std::endl;
}

数乘和数除

矩阵可以乘上或除以一个常数，相当于矩阵中每个数乘上或除以一个常数，主要包含五种操作：matrix*scalar、scalar*matrix、matrix/scalar、matrix*=scalar、matrix/=scalar

#include <iostream>
#include <Eigen/Dense>
 
int main()
{
  Eigen::Matrix2d a;
  a << 1, 2,
       3, 4;
  Eigen::Vector3d v(1,2,3);
  /*
  a * 2.5 =
  2.5 5
  7.5 10
  */
  std::cout << "a * 2.5 =\n" << a * 2.5 << std::endl;
  /*
  0.1 * v =
  0.1
  0.2
  0.3
  */
  std::cout << "0.1 * v =\n" << 0.1 * v << std::endl;
  std::cout << "Doing v *= 2;" << std::endl;
  v *= 2;
  /*
  Now v = 
  2
  4
  6
  */
  std::cout << "Now v =\n" << v << std::endl;
}

转置和共轭

转置 $a^T$ 使用 transpose()，共轭(针对复数，虚部取反，对实数无影响) $\bar{a}$ 使用 conjugation()，共轭转置(即进行共轭也进行转置) $a^{\ast }$ 使用 adjoint()。

MatrixXcf a = MatrixXcf::Random(2,2);
/*
Here is the matrix a
(-0.211,0.68) (-0.605,0.823)
(0.597,0.566)  (0.536,-0.33)
*/
cout << "Here is the matrix a\n" << a << endl;
/*
Here is the matrix a^T
(-0.211,0.68)  (0.597,0.566)
(-0.605,0.823)  (0.536,-0.33)
*/
cout << "Here is the matrix a^T\n" << a.transpose() << endl;
/*
Here is the conjugate of a
(-0.211,-0.68) (-0.605,-0.823)
(0.597,-0.566)    (0.536,0.33)
*/
cout << "Here is the conjugate of a\n" << a.conjugate() << endl;
/*
Here is the matrix a^*
(-0.211,-0.68)  (0.597,-0.566)
(-0.605,-0.823)    (0.536,0.33)
*/
cout << "Here is the matrix a^*\n" << a.adjoint() << endl;

需要注意的是，transpose() 是赋值与转置同时进行的，b=a.transpose() 没有什么问题，但是 a=a.transpose() 就可能会出现问题：

Matrix2i a; a << 1, 2, 3, 4;
/*
Here is the matrix a:
1 2
3 4
*/
cout << "Here is the matrix a:\n" << a << endl;
 
a = a.transpose(); // !!! do NOT do this !!!
/*
and the result of the aliasing effect:
1 2
2 4
*/
cout << "and the result of the aliasing effect:\n" << a << endl;

这种情况，需要使用 transposeInPlace()，对于共轭转置可以使用 adjointInPlace()：

MatrixXf a(2,3); a << 1, 2, 3, 4, 5, 6;
/*
Here is the initial matrix a:
1 2 3
4 5 6
*/
cout << "Here is the initial matrix a:\n" << a << endl;
 
 
a.transposeInPlace();
/*
and after being transposed:
1 4
2 5
3 6
*/
cout << "and after being transposed:\n" << a << endl;

矩阵-矩阵乘法 矩阵-向量乘法

矩阵-矩阵乘法与矩阵-向量乘法使用运算符 *，主要有两种操作：a*b、a*=b：

#include <iostream>
#include <Eigen/Dense>
 
int main()
{
  Eigen::Matrix2d mat;
  mat << 1, 2,
         3, 4;
  Eigen::Vector2d u(-1,1), v(2,0);
  //对应位置元素相乘 -2 0 
  std::cout << u.cwiseProduct(v) << std::endl;
  /*
  Here is mat*mat:
  7 10
  15 22
  */
  std::cout << "Here is mat*mat:\n" << mat*mat << std::endl;
  /*
  Here is mat*u:
  1
  1
  */
  std::cout << "Here is mat*u:\n" << mat*u << std::endl;
  /*
  Here is u^T*mat:
  2 2
  */
  std::cout << "Here is u^T*mat:\n" << u.transpose()*mat << std::endl;
  /*
  Here is u^T*v:
  -2
  */
  std::cout << "Here is u^T*v:\n" << u.transpose()*v << std::endl;
  /*
  Here is u*v^T:
 -2 -0
  2  0
  */
  std::cout << "Here is u*v^T:\n" << u*v.transpose() << std::endl;
  std::cout << "Let's multiply mat by itself" << std::endl;
  mat = mat*mat;
  /*
  Now mat is mat:
  7 10
  15 22
  */
  std::cout << "Now mat is mat:\n" << mat << std::endl;
}

矩阵乘法在 Eigen 是被特殊对待的，其不会出现类似转置的别名问题，其编译 m=m*m 像这样：tmp=m*m;m=tmp。如果确定没有别名问题，可以使用noalias()，比如 c.noalias() += a * b;

点积和叉积

点积和叉积分别使用 dot() 和 cross() 方法，当然点积也可以由 u.adjoint()*v 获得 1x1 矩阵

#include <iostream>
#include <Eigen/Dense>
 
int main()
{
  Eigen::Vector3d v(1,2,3);
  Eigen::Vector3d w(0,1,2);
  
  /*	
  Dot product: 8
  */
  std::cout << "Dot product: " << v.dot(w) << std::endl;
  double dp = v.adjoint()*w; // automatic conversion of the inner product to a scalar
  /*
  Dot product via a matrix product: 8
  */
  std::cout << "Dot product via a matrix product: " << dp << std::endl;
  /*
  Cross product:
  1
  -2
  1
  */
  std::cout << "Cross product:\n" << v.cross(w) << std::endl;
}

需要注意的是叉积只能应用于大小为3的向量。

基础算术约简操作

算术约简操作就是将一个矩阵或向量计算出一个单一的值，比如求和 sum()、乘积 prod()、最大值 maxCoeff()、最小值 minCoeff()

#include <iostream>
#include <Eigen/Dense>
 
using namespace std;
int main()
{
  Eigen::Matrix2d mat;
  mat << 1, 2,
         3, 4;
  // Here is mat.sum():       10
  cout << "Here is mat.sum():       " << mat.sum()       << endl;
  // Here is mat.prod():      24
  cout << "Here is mat.prod():      " << mat.prod()      << endl;
  // Here is mat.mean():      2.5
  cout << "Here is mat.mean():      " << mat.mean()      << endl;
  // Here is mat.minCoeff():  1
  cout << "Here is mat.minCoeff():  " << mat.minCoeff()  << endl;
  // Here is mat.maxCoeff():  4
  cout << "Here is mat.maxCoeff():  " << mat.maxCoeff()  << endl;
  // Here is mat.trace():     5
  cout << "Here is mat.trace():     " << mat.trace()     << endl;
}

我们可以使用参数来获取 minCoeff 和 maxCoeff 对应最大值最小值对于的坐标

  Matrix3f m = Matrix3f::Random();
  std::ptrdiff_t i, j;
  float minOfM = m.minCoeff(&i,&j);
  /*
  Here is the matrix m:
  0.68  0.597  -0.33
  -0.211  0.823  0.536
  0.566 -0.605 -0.444
 */
  cout << "Here is the matrix m:\n" << m << endl;
  // Its minimum coefficient (-0.605) is at position (2,1)
  cout << "Its minimum coefficient (" << minOfM 
       << ") is at position (" << i << "," << j << ")\n\n";
 
  RowVector4i v = RowVector4i::Random();
  int maxOfV = v.maxCoeff(&i);
  // Here is the vector v:  1  0  3 -3
  cout << "Here is the vector v: " << v << endl;
  // Its maximum coefficient (3) is at position 2
  cout << "Its maximum coefficient (" << maxOfV 
       << ") is at position " << i << endl;

---

作业0

作业描述

给定一个点 P=(2,1)，将该点绕原点先逆时针旋转 $45^{\circ }$，再平移(1,2)，计算出变换后点的坐标（要求用齐次坐标进行计算）。

思路与代码

P=(2,1)，用齐次坐标表示为 $\left[\begin{array}{c}2\\ 1\\ 1\end{array}\right]$，绕原点逆时针旋转 $45^{\circ }$，其旋转矩阵为$\left[\begin{array}{cc}cos45^{\circ }& -sin45^{\circ }\\ sin45^{\circ }& cos45^{\circ }\end{array}\right]=\left[\begin{array}{cc}\frac{1}{\sqrt{2}}& -\frac{1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}\end{array}\right]$，加上平移矩阵(齐次坐标)下的变换矩阵为 $\left[\begin{array}{ccc}\frac{1}{\sqrt{2}}& -\frac{1}{\sqrt{2}}& 1\\ \frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}& 2\\ 0& 0& 1\end{array}\right]$，因此变换之后的P的齐次坐标为 $\left[\begin{array}{ccc}\frac{1}{\sqrt{2}}& -\frac{1}{\sqrt{2}}& 1\\ \frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}& 2\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}2\\ 1\\ 1\end{array}\right]=\left[\begin{array}{c}\frac{\sqrt{2}+2}{2}\\ \frac{3\sqrt{2}+4}{2}\\ 1\end{array}\right]$。

main.cpp

#include <iostream>
#include <eigen3/Eigen/Dense>
#include <cmath>

using namespace std;

#define PI acos(-1)
int main() {
    Eigen::Vector3d P(2.0, 1.0, 1.0); // define P(2,1) in homogeneous coordinate
    Eigen::Matrix3d Transformation;
    Transformation << cos(PI/4), -sin(PI/4), 1.0, 
                        sin(PI/4), cos(PI/4), 2.0,
                        0.0, 0.0, 1.0;
    Eigen::Vector3d P_Transform = Transformation * P;
    cout << P_Transform << endl;
    return 0;
}

编译

在 main.cpp 所在目录下，打开终端(命令行)，依次输入：

• mkdir build：创建名为 build 的文件夹
• cd build：移动到 build 文件夹下
• cmake ..：注意其中 .. 表示上一级目录，若为 . 则表示当前目录，cmake具体使用可以看这篇文章
• make：编译程序，错误提示会显示在终端中，make具体使用可以看这篇文章
• .Transformation：若上一步无错误，则可运行程序（这里的 Transformation 为可执行文件名，可参照 CMakeList.txt 中修改）

结果为 $\left[\begin{array}{c}1.70711\\ 4.12132\\ 1\end{array}\right]$