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)
定義(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() | 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;
}