C++開源矩陣計算工具——Eigen的簡單用法（三）
阿新 • • 發佈：2019-01-07
本節主要涉及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）定義（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() | Invertible | ++ | + | |
| fullPivLu() | None | - | +++ | |
| householderQr() | None | ++ | + | |
| 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;
}