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

时间:2023-01-31 07:36:52

本节主要涉及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;
}
输出结果为:

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

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;
}