本节主要涉及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 303)向量的块操作,其实向量只是一个特殊的矩阵,但是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 | +++ | ++ |
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; }