C++开源矩阵计算工具——Eigen的简单用法（三）

本节主要涉及Eigen的块操作以及QR分解，Eigen的QR分解非常绕人，搞了很久才搞明白是怎么回事，最后是一个使用Eigen的矩阵操作完成二维高斯拟合求取光点的代码例子，关于二维高斯拟合求取光点的详细内容可参考：http://blog.csdn.net/hjx_1000/article/details/8490653

1、矩阵的块操作

        1）矩阵的块操作有两种使用方法，其定义形式为：

matrix.block(i,j,p,q);      （1）

matrix.block<p,q>(i,j);    （2）

定义（1）表示返回从矩阵的(i, j)开始，每行取p个元素，每列取q个元素所组成的临时新矩阵对象，原矩阵的元素不变。

定义（2）中block（p, q）可理解为一个p行q列的子矩阵，该定义表示从原矩阵中第(i, j)开始，获取一个p行q列的子矩阵，返回该子矩阵组成的临时 矩阵对象，原矩阵的元素不变。

详细使用情况，可参考下面的代码段：

#include <Eigen/Dense>
#include <iostream>
using namespace std;
int main()
{
Eigen::MatrixXf m(4,4);
m << 1, 2, 3, 4,
5, 6, 7, 8,
9,10,11,12,
13,14,15,16;
cout << "Block in the middle" << endl;
cout << m.block<2,2>(1,1) << endl << endl;
for (int i = 1; i <= 3; ++i)
{
cout << "Block of size " << i << "x" << i << endl;
cout << m.block(0,0,i,i) << endl << endl;
}
}

输出的结果为：

Block in the middle
 6  7
10 11

Block of size 1x1
1

Block of size 2x2
1 2
5 6

Block of size 3x3
 1  2  3
 5  6  7
 9 10 11

通过上述方式获取的子矩阵即可以作为左值也可以作为右值，也就是即可以用这个子矩阵给其他矩阵赋值，也可以给这个子矩阵对象赋值。

2）矩阵也提供了获取其指定行/列的函数，其实获取某行/列也是一种特殊的获取子块。可以通过 .col()和 .row()来完成获取指定列/行的操作，参数为列/行的索引。
注意：
（1）需与获取矩阵的行数/列数的函数（ rows(), cols() ）的进行区别，不要弄混淆。
（2）函数参数为响应行/列的索引，需注意矩阵的行列均以0开始。
下面的代码段用于演示获取矩阵的指定行列：

#include <Eigen/Dense>
#include <iostream>
using namespace std;
int main()
{
Eigen::MatrixXf m(3,3);
m << 1,2,3,
4,5,6,
7,8,9;
cout << "Here is the matrix m:" << endl << m << endl;
cout << "2nd Row: " << m.row(1) << endl;
m.col(2) += 3 * m.col(0);
cout << "After adding 3 times the first column into the third column, the matrix m is:\n";
cout << m << endl;
}

输出结果为：

Here is the matrix m:
1 2 3
4 5 6
7 8 9
2nd Row: 4 5 6
After adding 3 times the first column into the third column, the matrix m is:
 1  2  6
 4  5 18
 7  8 30

3）向量的块操作，其实向量只是一个特殊的矩阵，但是Eigen也为它单独提供了一些简化的块操作，如下三种形式：
     获取向量的前n个元素：vector.head(n);
     获取向量尾部的n个元素：vector.tail(n);
     获取从向量的第i个元素开始的n个元素：vector.segment(i,n);
     其用法可参考如下代码段：

#include <Eigen/Dense>
#include <iostream>
using namespace std;
int main()
{
Eigen::ArrayXf v(6);
v << 1, 2, 3, 4, 5, 6;
cout << "v.head(3) =" << endl << v.head(3) << endl << endl;
cout << "v.tail<3>() = " << endl << v.tail<3>() << endl << endl;
v.segment(1,4) *= 2;
cout << "after 'v.segment(1,4) *= 2', v =" << endl << v << endl;
}

输出结果为：

v.head(3) =
1
2
3

v.tail<3>() = 
4
5
6

after 'v.segment(1,4) *= 2', v =
1
4
6
8
10
6

2、QR分解
        Eigen的QR分解非常绕人，它总共提供了下面这些矩阵的分解方式：

Decomposition	Method	Requirements on the matrix	Speed	Accuracy
PartialPivLU	partialPivLu()	Invertible	++	+
FullPivLU	fullPivLu()	None	-	+++
HouseholderQR	householderQr()	None	++	+
ColPivHouseholderQR	colPivHouseholderQr()	None	+	++
FullPivHouseholderQR	fullPivHouseholderQr()	None	-	+++
LLT	llt()	Positive definite	+++	+
LDLT	ldlt()	Positive or negative semidefinite	+++	++

由于我只用到了QR分解，而且Eigen的QR分解开始使用时确实不容易入手，因此这里只提供了householderQR的分解方式的演示代码：

void QR2()
{
	Matrix3d A;
	A<<1,1,1,
		2,-1,-1,
		2,-4,5;

	HouseholderQR<Matrix3d> qr;
	qr.compute(A);
	MatrixXd R = qr.matrixQR().triangularView<Upper>();
	MatrixXd Q =  qr.householderQ();
	std::cout << "QR2(): HouseholderQR---------------------------------------------"<< std::endl;
	std::cout << "A "<< std::endl <<A << std::endl << std::endl;
	std::cout <<"qr.matrixQR()"<< std::endl << qr.matrixQR() << std::endl << std::endl;
	std::cout << "R"<< std::endl <<R << std::endl << std::endl;
	std::cout << "Q "<< std::endl <<Q << std::endl << std::endl;
	std::cout <<"Q*R" << std::endl <<Q*R << std::endl << std::endl;
}

输出结果为：

3、一个矩阵使用的例子：用矩阵操作完成二维高斯拟合，并求取光斑中心

下面的代码段是一个使用Eigen的矩阵操作完成二维高斯拟合求取光点的代码例子，关于二维高斯拟合求取光点的详细内容可参考：http://blog.csdn.net/hjx_1000/article/details/8490653

bool GetCentrePoint(float& x0,float& y0)
{
	if (m_iN<=0)
		return false;
	//QR分解
	HouseholderQR<MatrixXf> qr;
	qr.compute(m_matrix_B);
	MatrixXf R = qr.matrixQR().triangularView<Upper>();
	MatrixXf Q =  qr.householderQ();

	//块操作，获取向量或矩阵的局部
	VectorXf S;
	S = (Q.transpose()* m_Vector_A).head(5);
	MatrixXf R1;
	R1 = R.block(0,0,5,5);

	VectorXf C;
	C = R1.inverse() * S;

	x0 = -0.5 * C[1] / C[3];
	y0 = -0.5 * C[2] / C[4];

	return true;
}