master
laobai 2022-04-19 22:39:30 +08:00
parent ae00294b52
commit 4cf21e5ed5
1 changed files with 1201 additions and 907 deletions

View File

@ -25,7 +25,6 @@
#include <thread>
using namespace std;
namespace ORB_SLAM3
{
@ -38,9 +37,20 @@ namespace ORB_SLAM3
mMaxIterations = iterations;
}
bool TwoViewReconstruction::Reconstruct(const std::vector<cv::KeyPoint>& vKeys1, const std::vector<cv::KeyPoint>& vKeys2, const vector<int> &vMatches12,
/**
* @brief rt3d
* @param vKeys1
* @param vKeys2
* @param vMatches12 vKeys1vKeys2
* @param T21
* @param vP3D
* @param vbTriangulated
*/
bool TwoViewReconstruction::Reconstruct(
const std::vector<cv::KeyPoint> &vKeys1, const std::vector<cv::KeyPoint> &vKeys2, const vector<int> &vMatches12,
Sophus::SE3f &T21, vector<cv::Point3f> &vP3D, vector<bool> &vbTriangulated)
{
// 1. 准备工作提取匹配关系及准备RANSAC
mvKeys1.clear();
mvKeys2.clear();
@ -49,8 +59,11 @@ namespace ORB_SLAM3
// Fill structures with current keypoints and matches with reference frame
// Reference Frame: 1, Current Frame: 2
// 填写mvMatches12里面存放的是vKeys1vKeys2匹配点的索引
// 存放匹配点的id
mvMatches12.clear();
mvMatches12.reserve(mvKeys2.size());
// 长度与vKeys1表示对应位置的vKeys1中的点是否有匹配关系
mvbMatched1.resize(mvKeys1.size());
for (size_t i = 0, iend = vMatches12.size(); i < iend; i++)
{
@ -70,16 +83,19 @@ namespace ORB_SLAM3
vAllIndices.reserve(N);
vector<size_t> vAvailableIndices;
// 使用vAllIndices为了保证8个点选不到同一个点
for (int i = 0; i < N; i++)
{
vAllIndices.push_back(i);
}
// Generate sets of 8 points for each RANSAC iteration
// 默认200次
mvSets = vector<vector<size_t>>(mMaxIterations, vector<size_t>(8, 0));
DUtils::Random::SeedRandOnce(0);
// 2. 先遍历把200次先取好
for (int it = 0; it < mMaxIterations; it++)
{
vAvailableIndices = vAllIndices;
@ -88,10 +104,11 @@ namespace ORB_SLAM3
for (size_t j = 0; j < 8; j++)
{
int randi = DUtils::Random::RandomInt(0, vAvailableIndices.size() - 1);
int idx = vAvailableIndices[randi];
int idx = vAvailableIndices[randi]; // 这句不多余,防止重复选择
mvSets[it][j] = idx;
// 保证选不到同一个点这么做的话可以删去vAvailableIndices已选点的索引
vAvailableIndices[randi] = vAvailableIndices.back();
vAvailableIndices.pop_back();
}
@ -102,6 +119,8 @@ namespace ORB_SLAM3
float SH, SF;
Eigen::Matrix3f H, F;
// 3. 双线程分别计算
// 加ref为了提供引用
thread threadH(&TwoViewReconstruction::FindHomography, this, ref(vbMatchesInliersH), ref(SH), ref(H));
thread threadF(&TwoViewReconstruction::FindFundamental, this, ref(vbMatchesInliersF), ref(SF), ref(F));
@ -110,12 +129,15 @@ namespace ORB_SLAM3
threadF.join();
// Compute ratio of scores
if(SH+SF == 0.f) return false;
if (SH + SF == 0.f)
return false;
float RH = SH / (SH + SF);
float minParallax = 1.0;
// Try to reconstruct from homography or fundamental depending on the ratio (0.40-0.45)
// 4. 看哪个分高用哪个分别是通过H重建与通过F重建
// ORBSLAM2这里的值是0.4 TOSEE
if (RH > 0.50) // if(RH>0.40)
{
// cout << "Initialization from Homography" << endl;
@ -128,14 +150,22 @@ namespace ORB_SLAM3
}
}
/**
* @brief H
* @param vbMatchesInliers HmvMatches12
* @param score
* @param H21 12H
*/
void TwoViewReconstruction::FindHomography(vector<bool> &vbMatchesInliers, float &score, Eigen::Matrix3f &H21)
{
// Number of putative matches
// 匹配成功的个数
const int N = mvMatches12.size();
// Normalize coordinates
vector<cv::Point2f> vPn1, vPn2;
Eigen::Matrix3f T1, T2;
// 像素坐标标准化,先去均值,再除均长
Normalize(mvKeys1, vPn1, T1);
Normalize(mvKeys2, vPn2, T2);
Eigen::Matrix3f T2inv = T2.inverse();
@ -152,6 +182,7 @@ namespace ORB_SLAM3
float currentScore;
// Perform all RANSAC iterations and save the solution with highest score
// 计算归一化后的H矩阵 p2 = H21*p1
for (int it = 0; it < mMaxIterations; it++)
{
// Select a minimum set
@ -162,12 +193,13 @@ namespace ORB_SLAM3
vPn1i[j] = vPn1[mvMatches12[idx].first];
vPn2i[j] = vPn2[mvMatches12[idx].second];
}
// 计算标准化后的H矩阵
Eigen::Matrix3f Hn = ComputeH21(vPn1i, vPn2i);
// 恢复正常H
H21i = T2inv * Hn * T1;
H12i = H21i.inverse();
currentScore = CheckHomography(H21i, H12i, vbCurrentInliers, mSigma);
currentScore = CheckHomography(H21i, H12i, vbCurrentInliers, mSigma); // mSigma默认为1
if (currentScore > score)
{
@ -178,7 +210,12 @@ namespace ORB_SLAM3
}
}
/**
* @brief F
* @param vbMatchesInliers FmvMatches12
* @param score
* @param F21 12F
*/
void TwoViewReconstruction::FindFundamental(vector<bool> &vbMatchesInliers, float &score, Eigen::Matrix3f &F21)
{
// Number of putative matches
@ -216,6 +253,7 @@ namespace ORB_SLAM3
Eigen::Matrix3f Fn = ComputeF21(vPn1i, vPn2i);
// FN得到的是基于归一化坐标的F矩阵必须乘上归一化过程矩阵才是最后的基于像素坐标的F
F21i = T2t * Fn * T1;
currentScore = CheckFundamental(F21i, vbCurrentInliers, mSigma);
@ -229,6 +267,23 @@ namespace ORB_SLAM3
}
}
/**
* @brief homographynormalized DLT
* |x'| | h1 h2 h3 ||x|
* |y'| = a | h4 h5 h6 ||y| : x' = a H x, a
* |1 | | h7 h8 h9 ||1|
* 使DLT(direct linear tranform)
* x' = a H x
* ---> (x') (H x) = 0
* ---> Ah = 0
* A = | 0 0 0 -x -y -1 xy' yy' y'| h = | h1 h2 h3 h4 h5 h6 h7 h8 h9 |
* |-x -y -1 0 0 0 xx' yx' x'|
* SVDAh = 0A'A
* @param vP1 , in reference frame
* @param vP2 , in current frame
* @return
* @see Multiple View Geometry in Computer Vision - Algorithm 4.2 p109
*/
Eigen::Matrix3f TwoViewReconstruction::ComputeH21(const vector<cv::Point2f> &vP1, const vector<cv::Point2f> &vP2)
{
const int N = vP1.size();
@ -261,7 +316,6 @@ namespace ORB_SLAM3
A(2 * i + 1, 6) = -u2 * u1;
A(2 * i + 1, 7) = -u2 * v1;
A(2 * i + 1, 8) = -u2;
}
Eigen::JacobiSVD<Eigen::MatrixXf> svd(A, Eigen::ComputeFullV);
@ -271,6 +325,16 @@ namespace ORB_SLAM3
return H;
}
/**
* @brief fundamental matrixnormalized 8
* x'Fx = 0 Af = 0
* A = | x'x x'y x' y'x y'y y' x y 1 |, f = | f1 f2 f3 f4 f5 f6 f7 f8 f9 |
* SVDAf = 0A'A
* @param vP1 , in reference frame
* @param vP2 , in current frame
* @return
* @see Multiple View Geometry in Computer Vision - Algorithm 11.1 p282 ( p191)
*/
Eigen::Matrix3f TwoViewReconstruction::ComputeF21(const vector<cv::Point2f> &vP1, const vector<cv::Point2f> &vP2)
{
const int N = vP1.size();
@ -302,11 +366,19 @@ namespace ORB_SLAM3
Eigen::JacobiSVD<Eigen::Matrix3f> svd2(Fpre, Eigen::ComputeFullU | Eigen::ComputeFullV);
Eigen::Vector3f w = svd2.singularValues();
// 这里注意计算完要强制让第三个奇异值为0
w(2) = 0;
return svd2.matrixU() * Eigen::DiagonalMatrix<float, 3>(w) * svd2.matrixV().transpose();
}
/**
* @brief
* @param H21
* @param H12
* @param vbMatchesInliers mvMatches12
* @param sigma 1
*/
float TwoViewReconstruction::CheckHomography(const Eigen::Matrix3f &H21, const Eigen::Matrix3f &H12, vector<bool> &vbMatchesInliers, float sigma)
{
const int N = mvMatches12.size();
@ -334,7 +406,7 @@ namespace ORB_SLAM3
vbMatchesInliers.resize(N);
float score = 0;
// 基于卡方检验计算出的阈值 自由度为2的卡方分布显著性水平为0.05,对应的临界阈值
const float th = 5.991;
const float invSigmaSquare = 1.0 / (sigma * sigma);
@ -354,6 +426,7 @@ namespace ORB_SLAM3
// Reprojection error in first image
// x2in1 = H12*x2
// 计算投影误差2投1 1投2这么做计算累计的卡方检验分数分数越高证明内点与误差越优这么做为了平衡误差与内点个数不是说内点个数越高越好也不是说误差越小越好
const float w2in1inv = 1.0 / (h31inv * u2 + h32inv * v2 + h33inv);
const float u2in1 = (h11inv * u2 + h12inv * v2 + h13inv) * w2in1inv;
const float v2in1 = (h21inv * u2 + h22inv * v2 + h23inv) * w2in1inv;
@ -392,6 +465,12 @@ namespace ORB_SLAM3
return score;
}
/**
* @brief
* @param F21
* @param vbMatchesInliers mvMatches12
* @param sigma 1
*/
float TwoViewReconstruction::CheckFundamental(const Eigen::Matrix3f &F21, vector<bool> &vbMatchesInliers, float sigma)
{
const int N = mvMatches12.size();
@ -410,7 +489,9 @@ namespace ORB_SLAM3
float score = 0;
// 基于卡方检验计算出的阈值 自由度为1的卡方分布显著性水平为0.05,对应的临界阈值
const float th = 3.841;
// 基于卡方检验计算出的阈值 自由度为2的卡方分布显著性水平为0.05,对应的临界阈值
const float thScore = 5.991;
const float invSigmaSquare = 1.0 / (sigma * sigma);
@ -429,17 +510,19 @@ namespace ORB_SLAM3
// Reprojection error in second image
// l2=F21x1=(a2,b2,c2)
// 计算 img1 上的点在 img2 上投影得到的极线 l2 = F21 * p1 = (a2,b2,c2)
const float a2 = f11 * u1 + f12 * v1 + f13;
const float b2 = f21 * u1 + f22 * v1 + f23;
const float c2 = f31 * u1 + f32 * v1 + f33;
// 计算误差 e = (a * p2.x + b * p2.y + c) / sqrt(a * a + b * b)
const float num2 = a2 * u2 + b2 * v2 + c2;
const float squareDist1 = num2 * num2 / (a2 * a2 + b2 * b2);
const float chiSquare1 = squareDist1 * invSigmaSquare;
// 自由度为1是因为这里的计算是点到线的距离判定分数自由度为2的原因可能是为了与H矩阵持平
if (chiSquare1 > th)
bIn = false;
else
@ -447,7 +530,7 @@ namespace ORB_SLAM3
// Reprojection error in second image
// l1 =x2tF21=(a1,b1,c1)
// 与上面相同只不过反过来了
const float a1 = f11 * u2 + f21 * v2 + f31;
const float b1 = f12 * u2 + f22 * v2 + f32;
const float c1 = f13 * u2 + f23 * v2 + f33;
@ -472,27 +555,43 @@ namespace ORB_SLAM3
return score;
}
bool TwoViewReconstruction::ReconstructF(vector<bool> &vbMatchesInliers, Eigen::Matrix3f &F21, Eigen::Matrix3f &K,
/**
* @brief
* @param vbMatchesInliers mvMatches12
* @param F21
* @param K
* @param T21
* @param vP3D
* @param vbTriangulated mvKeys1
* @param minParallax 1
* @param minTriangulated 50
*/
bool TwoViewReconstruction::ReconstructF(
vector<bool> &vbMatchesInliers, Eigen::Matrix3f &F21, Eigen::Matrix3f &K,
Sophus::SE3f &T21, vector<cv::Point3f> &vP3D, vector<bool> &vbTriangulated, float minParallax, int minTriangulated)
{
// 统计了合法的匹配,后面用于对比重建出的点数
int N = 0;
for (size_t i = 0, iend = vbMatchesInliers.size(); i < iend; i++)
if (vbMatchesInliers[i])
N++;
// Compute Essential Matrix from Fundamental Matrix
// 1. 计算本质矩阵
Eigen::Matrix3f E21 = K.transpose() * F21 * K;
Eigen::Matrix3f R1, R2;
Eigen::Vector3f t;
// Recover the 4 motion hypotheses
// 2. 分解本质矩阵得到4对rt
DecomposeE(E21, R1, R2, t);
Eigen::Vector3f t1 = t;
Eigen::Vector3f t2 = -t;
// Reconstruct with the 4 hyphoteses and check
// 3. 使用四对结果分别重建
vector<cv::Point3f> vP3D1, vP3D2, vP3D3, vP3D4;
vector<bool> vbTriangulated1, vbTriangulated2, vbTriangulated3, vbTriangulated4;
float parallax1, parallax2, parallax3, parallax4;
@ -501,11 +600,13 @@ namespace ORB_SLAM3
int nGood2 = CheckRT(R2, t1, mvKeys1, mvKeys2, mvMatches12, vbMatchesInliers, K, vP3D2, 4.0 * mSigma2, vbTriangulated2, parallax2);
int nGood3 = CheckRT(R1, t2, mvKeys1, mvKeys2, mvMatches12, vbMatchesInliers, K, vP3D3, 4.0 * mSigma2, vbTriangulated3, parallax3);
int nGood4 = CheckRT(R2, t2, mvKeys1, mvKeys2, mvMatches12, vbMatchesInliers, K, vP3D4, 4.0 * mSigma2, vbTriangulated4, parallax4);
// 统计重建出点的数量最大值
int maxGood = max(nGood1, max(nGood2, max(nGood3, nGood4)));
// 起码要重建出超过百分之90的匹配点
int nMinGood = max(static_cast<int>(0.9 * N), minTriangulated);
// 4. 看看有没有脱颖而出的结果且最大值要高于要求的最低值如果大家都差不多说明有问题因为4个结果中只有一个会正常
// 大家都很多的话反而不正常了。。。
int nsimilar = 0;
if (nGood1 > 0.7 * maxGood)
nsimilar++;
@ -523,6 +624,7 @@ namespace ORB_SLAM3
}
// If best reconstruction has enough parallax initialize
// 5. 使用最好的结果输出
if (maxGood == nGood1)
{
if (parallax1 > minParallax)
@ -533,7 +635,8 @@ namespace ORB_SLAM3
T21 = Sophus::SE3f(R1, t1);
return true;
}
}else if(maxGood==nGood2)
}
else if (maxGood == nGood2)
{
if (parallax2 > minParallax)
{
@ -543,7 +646,8 @@ namespace ORB_SLAM3
T21 = Sophus::SE3f(R2, t1);
return true;
}
}else if(maxGood==nGood3)
}
else if (maxGood == nGood3)
{
if (parallax3 > minParallax)
{
@ -553,7 +657,8 @@ namespace ORB_SLAM3
T21 = Sophus::SE3f(R1, t2);
return true;
}
}else if(maxGood==nGood4)
}
else if (maxGood == nGood4)
{
if (parallax4 > minParallax)
{
@ -568,9 +673,81 @@ namespace ORB_SLAM3
return false;
}
bool TwoViewReconstruction::ReconstructH(vector<bool> &vbMatchesInliers, Eigen::Matrix3f &H21, Eigen::Matrix3f &K,
/**
* @brief HR t
* HFaugeras SVD-based decomposition Zhang SVD-based decomposition
* Motion and structure from motion in a piecewise plannar environment
* 使Faugeras SVD-based decomposition
* @see
* - Faugeras et al, Motion and structure from motion in a piecewise planar environment. International Journal of Pattern Recognition and Artificial Intelligence, 1988.
* - Deeper understanding of the homography decomposition for vision-based control
* n = (a, b, c)^t (x, y, z)ax + by + cz = d 1/d * n^t * x = 1 d
* x' = R*x + t x x'
* λ2*x' = R*(λ1*x) + t = R*(λ1*x) + t * 1/d * n^t * (λ1*x)
* x' = λ*(R + t * 1/d * n^t) * x = H^ * x
* u' = G * u G = KH^K.inv
* H^ ~= d*R + t * n^t = UV^t = U^t * H^ * V = d*U^t * R * V + (U^t * t) * (V^t * n)^t
* s = det(U) * det(V) s 1 -1 = s^2 * d*U^t * R * V + (U^t * t) * (V^t * n)^t = (s*d) * s * U^t * R * V + (U^t * t) * (V^t * n)^t
* R' = s * U^t * R * V t' = U^t * t n' = V^t * n d' = s * d
* = d' * R' + t' * n'^t s
* = | d1, 0, 0 | e1 = (1, 0, 0)^t e2 = (0, 1, 0)^t e3 = (0, 0, 1)^t
* | 0, d2, 0 |
* | 0, 0, d3 |
* n' = (a1, 0, 0)^t + (0, b1, 0)^t + (0, 0, c1)^t = a1*e1 + b1*e2 + c1*e3
*
* = [d1*e1, d2*e2, d3*e3] = [d' * R' * e1, d' * R' * e2, d' * R' * e3] + [t' * a1, t' * b1, t' * c1]
* ==> d1*e1 = d' * R' * e1 + t' * a1
* d2*e2 = d' * R' * e2 + t' * b1
* d3*e3 = d' * R' * e3 + t' * c1
*
*
* t
* d'R'(b1e1 - a1e2) = d1b1e1 - d2a1e2
* d'R'(c1e2 - b1e3) = d2c1e1 - d3b1e3
* d'R'(a1e3 - c1e1) = d3a1e3 - d1c1e1
*
* (d'^2 - d2^2)*a1^2 + (d'^2 - d1^2)*b1^2 = 0
* (d'^2 - d3^2)*b1^2 + (d'^2 - d2^2)*c1^2 = 0 d'^2 - d1^2 = x1 d'^2 - d2^2 = x2 d'^2 - d3^2 = x3
* (d'^2 - d1^2)*c1^2 + (d'^2 - d3^2)*a1^2 = 0
*
*
* | x2 x1 0 | | a1^2 |
* | 0 x3 x2 | * | b1^2 | = 0 ===> x1 * x2 * x3 = 0 (d'^2 - d1^2) * (d'^2 - d2^2) * (d'^2 - d3^2) = 0
* | x3 0 x1 | | c1^2 |
* d1 >= d2 >= d3 d' = d2 or -d2
*
* -----------------------------------------------------------------------------------------------------------------------------------------------------------------
* d' > 0 a1^2 + b1^2 + c1^2 = 1 ??????
* a1 = ε1 * sqrt((d1^2 - d2^2) / (d1^2 - d3^2))
* b1 = 0
* c1 = ε2 * sqrt((d2^2 - d3^2) / (d1^2 - d3^2)) ε1 ε2 1
* d2*e2 = d' * R' * e2 + t' * b1 => R' * e2 = e2
* | cosθ 0 -sinθ |
* ==> R' = | 0 1 0 | n' R' d'R'(a1e3 - c1e1) = d3a1e3 - d1c1e1
* | sinθ 0 cosθ |
* | cosθ 0 -sinθ | | -c1 | | -d1c1 |
* d' * | 0 1 0 | * | 0 | = | 0 | sinθ cosθ
* | sinθ 0 cosθ | | a1 | | d3a1 |
*
* R' d1*e1 = d' * R' * e1 + t' * a1
* d2*e2 = d' * R' * e2 + t' * b1
* d3*e3 = d' * R' * e3 + t' * c1
*
* t' = (d1 - d3) * (a1, 0, c1)^t
* @param vbMatchesInliers mvMatches12
* @param H21
* @param K
* @param T21
* @param vP3D
* @param vbTriangulated vbMatchesInliers
* @param minParallax 1
* @param minTriangulated 50
*/
bool TwoViewReconstruction::ReconstructH(
vector<bool> &vbMatchesInliers, Eigen::Matrix3f &H21, Eigen::Matrix3f &K,
Sophus::SE3f &T21, vector<cv::Point3f> &vP3D, vector<bool> &vbTriangulated, float minParallax, int minTriangulated)
{
// 统计了合法的匹配,后面用于对比重建出的点数
int N = 0;
for (size_t i = 0, iend = vbMatchesInliers.size(); i < iend; i++)
if (vbMatchesInliers[i])
@ -579,6 +756,8 @@ namespace ORB_SLAM3
// We recover 8 motion hypotheses using the method of Faugeras et al.
// Motion and structure from motion in a piecewise planar environment.
// International Journal of Pattern Recognition and Artificial Intelligence, 1988
// step1SVD分解Homography
// 因为特征点是图像坐标系所以讲H矩阵由相机坐标系换算到图像坐标系
Eigen::Matrix3f invK = K.inverse();
Eigen::Matrix3f A = invK * H21 * K;
@ -594,6 +773,7 @@ namespace ORB_SLAM3
float d2 = w(1);
float d3 = w(2);
// SVD分解的正常情况是特征值降序排列
if (d1 / d2 < 1.00001 || d2 / d3 < 1.00001)
{
return false;
@ -606,17 +786,38 @@ namespace ORB_SLAM3
vn.reserve(8);
// n'=[x1 0 x3] 4 posibilities e1=e3=1, e1=1 e3=-1, e1=-1 e3=1, e1=e3=-1
// step2计算法向量
// 法向量n'= [x1 0 x3]
float aux1 = sqrt((d1 * d1 - d2 * d2) / (d1 * d1 - d3 * d3));
float aux3 = sqrt((d2 * d2 - d3 * d3) / (d1 * d1 - d3 * d3));
float x1[] = {aux1, aux1, -aux1, -aux1};
float x3[] = {aux3, -aux3, aux3, -aux3};
// case d'=d2
// step3恢复旋转矩阵
// step3.1:计算 sin(theta)和cos(theta)case d'=d2
float aux_stheta = sqrt((d1 * d1 - d2 * d2) * (d2 * d2 - d3 * d3)) / ((d1 + d3) * d2);
float ctheta = (d2 * d2 + d1 * d3) / ((d1 + d3) * d2);
float stheta[] = {aux_stheta, -aux_stheta, -aux_stheta, aux_stheta};
// step3.2计算四种旋转矩阵Rt
// 计算旋转矩阵 R计算ppt中公式18
// | ctheta 0 -aux_stheta| | aux1|
// Rp = | 0 1 0 | tp = | 0 |
// | aux_stheta 0 ctheta | |-aux3|
// | ctheta 0 aux_stheta| | aux1|
// Rp = | 0 1 0 | tp = | 0 |
// |-aux_stheta 0 ctheta | | aux3|
// | ctheta 0 aux_stheta| |-aux1|
// Rp = | 0 1 0 | tp = | 0 |
// |-aux_stheta 0 ctheta | |-aux3|
// | ctheta 0 -aux_stheta| |-aux1|
// Rp = | 0 1 0 | tp = | 0 |
// | aux_stheta 0 ctheta | | aux3|
for (int i = 0; i < 4; i++)
{
Eigen::Matrix3f Rp;
@ -636,6 +837,8 @@ namespace ORB_SLAM3
tp(2) = -x3[i];
tp *= d1 - d3;
// 这里虽然对t有归一化并没有决定单目整个SLAM过程的尺度
// 因为CreateInitialMapMonocular函数对3D点深度会缩放然后反过来对 t 有改变
Eigen::Vector3f t = U * tp;
vt.push_back(t / t.norm());
@ -651,11 +854,14 @@ namespace ORB_SLAM3
}
// case d'=-d2
// step3.3:计算 sin(theta)和cos(theta)case d'=-d2
float aux_sphi = sqrt((d1 * d1 - d2 * d2) * (d2 * d2 - d3 * d3)) / ((d1 - d3) * d2);
float cphi = (d1 * d3 - d2 * d2) / ((d1 - d3) * d2);
float sphi[] = {aux_sphi, -aux_sphi, -aux_sphi, aux_sphi};
// step3.4计算四种旋转矩阵Rt
// 计算旋转矩阵 R
for (int i = 0; i < 4; i++)
{
Eigen::Matrix3f Rp;
@ -689,7 +895,6 @@ namespace ORB_SLAM3
vn.push_back(n);
}
int bestGood = 0;
int secondBestGood = 0;
int bestSolutionIdx = -1;
@ -699,6 +904,7 @@ namespace ORB_SLAM3
// Instead of applying the visibility constraints proposed in the Faugeras' paper (which could fail for points seen with low parallax)
// We reconstruct all hypotheses and check in terms of triangulated points and parallax
// step4d'=d2和d'=-d2分别对应8组(R t)通过恢复3D点并判断是否在相机正前方的方法来确定最优解
for (size_t i = 0; i < 8; i++)
{
float parallaxi;
@ -706,6 +912,7 @@ namespace ORB_SLAM3
vector<bool> vbTriangulatedi;
int nGood = CheckRT(vR[i], vt[i], mvKeys1, mvKeys2, mvMatches12, vbMatchesInliers, K, vP3Di, 4.0 * mSigma2, vbTriangulatedi, parallaxi);
// 保留最优的和次优的
if (nGood > bestGood)
{
secondBestGood = bestGood;
@ -721,7 +928,7 @@ namespace ORB_SLAM3
}
}
// step5通过判断最优是否明显好于次优从而判断该次Homography分解是否成功
if (secondBestGood < 0.75 * bestGood && bestParallax >= minParallax && bestGood > minTriangulated && bestGood > 0.9 * N)
{
T21 = Sophus::SE3f(vR[bestSolutionIdx], vt[bestSolutionIdx]);
@ -733,7 +940,12 @@ namespace ORB_SLAM3
return false;
}
/**
* @brief T
* @param vKeys
* @param vNormalizedPoints
* @param T
*/
void TwoViewReconstruction::Normalize(const vector<cv::KeyPoint> &vKeys, vector<cv::Point2f> &vNormalizedPoints, Eigen::Matrix3f &T)
{
float meanX = 0;
@ -748,6 +960,7 @@ namespace ORB_SLAM3
meanY += vKeys[i].pt.y;
}
// 1. 求均值
meanX = meanX / N;
meanY = meanY / N;
@ -763,9 +976,11 @@ namespace ORB_SLAM3
meanDevY += fabs(vNormalizedPoints[i].y);
}
// 2. 确定新原点后计算与新原点的距离均值
meanDevX = meanDevX / N;
meanDevY = meanDevY / N;
// 3. 去均值化
float sX = 1.0 / meanDevX;
float sY = 1.0 / meanDevY;
@ -775,6 +990,7 @@ namespace ORB_SLAM3
vNormalizedPoints[i].y = vNormalizedPoints[i].y * sY;
}
// 4. 计算变化矩阵
T.setZero();
T(0, 0) = sX;
T(1, 1) = sY;
@ -783,7 +999,22 @@ namespace ORB_SLAM3
T(2, 2) = 1.f;
}
int TwoViewReconstruction::CheckRT(const Eigen::Matrix3f &R, const Eigen::Vector3f &t, const vector<cv::KeyPoint> &vKeys1, const vector<cv::KeyPoint> &vKeys2,
/**
* @brief cheirality checkF
* @param R
* @param t
* @param vKeys1
* @param vKeys2
* @param vMatches12
* @param vbMatchesInliers
* @param K
* @param vP3D
* @param th2
* @param vbGood mvKeys1
* @param parallax
*/
int TwoViewReconstruction::CheckRT(
const Eigen::Matrix3f &R, const Eigen::Vector3f &t, const vector<cv::KeyPoint> &vKeys1, const vector<cv::KeyPoint> &vKeys2,
const vector<Match> &vMatches12, vector<bool> &vbMatchesInliers,
const Eigen::Matrix3f &K, vector<cv::Point3f> &vP3D, float th2, vector<bool> &vbGood, float &parallax)
{
@ -800,6 +1031,8 @@ namespace ORB_SLAM3
vCosParallax.reserve(vKeys1.size());
// Camera 1 Projection Matrix K[I|0]
// 步骤1得到一个相机的投影矩阵
// 以第一个相机的光心作为世界坐标系
Eigen::Matrix<float, 3, 4> P1;
P1.setZero();
P1.block<3, 3>(0, 0) = K;
@ -808,11 +1041,13 @@ namespace ORB_SLAM3
O1.setZero();
// Camera 2 Projection Matrix K[R|t]
// 步骤2得到第二个相机的投影矩阵
Eigen::Matrix<float, 3, 4> P2;
P2.block<3, 3>(0, 0) = R;
P2.block<3, 1>(0, 3) = t;
P2 = K * P2;
// 第二个相机的光心在世界坐标系下的坐标
Eigen::Vector3f O2 = -R.transpose() * t;
int nGood = 0;
@ -822,6 +1057,7 @@ namespace ORB_SLAM3
if (!vbMatchesInliers[i])
continue;
// kp1和kp2是匹配特征点
const cv::KeyPoint &kp1 = vKeys1[vMatches12[i].first];
const cv::KeyPoint &kp2 = vKeys2[vMatches12[i].second];
@ -829,9 +1065,9 @@ namespace ORB_SLAM3
Eigen::Vector3f x_p1(kp1.pt.x, kp1.pt.y, 1);
Eigen::Vector3f x_p2(kp2.pt.x, kp2.pt.y, 1);
// 步骤3利用三角法恢复三维点p3dC1
GeometricTools::Triangulate(x_p1, x_p2, P1, P2, p3dC1);
if (!isfinite(p3dC1(0)) || !isfinite(p3dC1(1)) || !isfinite(p3dC1(2)))
{
vbGood[vMatches12[i].first] = false;
@ -839,6 +1075,7 @@ namespace ORB_SLAM3
}
// Check parallax
// 步骤4计算视差角余弦值
Eigen::Vector3f normal1 = p3dC1 - O1;
float dist1 = normal1.norm();
@ -847,17 +1084,22 @@ namespace ORB_SLAM3
float cosParallax = normal1.dot(normal2) / (dist1 * dist2);
// 步骤5判断3D点是否在两个摄像头前方
// Check depth in front of first camera (only if enough parallax, as "infinite" points can easily go to negative depth)
// 步骤5.13D点深度为负在第一个摄像头后方淘汰
if (p3dC1(2) <= 0 && cosParallax < 0.99998)
continue;
// Check depth in front of second camera (only if enough parallax, as "infinite" points can easily go to negative depth)
// 步骤5.23D点深度为负在第二个摄像头后方淘汰
Eigen::Vector3f p3dC2 = R * p3dC1 + t;
if (p3dC2(2) <= 0 && cosParallax < 0.99998)
continue;
// 步骤6计算重投影误差
// Check reprojection error in first image
// 计算3D点在第一个图像上的投影误差
float im1x, im1y;
float invZ1 = 1.0 / p3dC1(2);
im1x = fx * p3dC1(0) * invZ1 + cx;
@ -865,10 +1107,13 @@ namespace ORB_SLAM3
float squareError1 = (im1x - kp1.pt.x) * (im1x - kp1.pt.x) + (im1y - kp1.pt.y) * (im1y - kp1.pt.y);
// 步骤6.1:重投影误差太大,跳过淘汰
// 一般视差角比较小时重投影误差比较大
if (squareError1 > th2)
continue;
// Check reprojection error in second image
// 计算3D点在第二个图像上的投影误差
float im2x, im2y;
float invZ2 = 1.0 / p3dC2(2);
im2x = fx * p3dC2(0) * invZ2 + cx;
@ -876,9 +1121,12 @@ namespace ORB_SLAM3
float squareError2 = (im2x - kp2.pt.x) * (im2x - kp2.pt.x) + (im2y - kp2.pt.y) * (im2y - kp2.pt.y);
// 步骤6.1:重投影误差太大,跳过淘汰
// 一般视差角比较小时重投影误差比较大
if (squareError2 > th2)
continue;
// 步骤7统计经过检验的3D点个数记录3D点视差角
vCosParallax.push_back(cosParallax);
vP3D[vMatches12[i].first] = cv::Point3f(p3dC1(0), p3dC1(1), p3dC1(2));
nGood++;
@ -887,11 +1135,17 @@ namespace ORB_SLAM3
vbGood[vMatches12[i].first] = true;
}
// 7 得到3D点中较小的视差角并且转换成为角度制表示
if (nGood > 0)
{
// 从小到大排序注意vCosParallax值越大视差越小
sort(vCosParallax.begin(), vCosParallax.end());
// !排序后并没有取最小的视差角,而是取一个较小的视差角
// 作者的做法如果经过检验过后的有效3D点小于50个那么就取最后那个最小的视差角(cos值最大)
// 如果大于50个就取排名第50个的较小的视差角即可为了避免3D点太多时出现太小的视差角
size_t idx = min(50, int(vCosParallax.size() - 1));
// 将这个选中的角弧度制转换为角度制
parallax = acos(vCosParallax[idx]) * 180 / CV_PI;
}
else
@ -900,11 +1154,50 @@ namespace ORB_SLAM3
return nGood;
}
/**
* @brief Essential
* https://blog.csdn.net/weixin_44580210/article/details/90344511
* FEssentialE4
* 4[R1,t],[R1,-t],[R2,t],[R2,-t]
* ##
* 3×30
* E=[t]×R=SRS
* 1 S S=UBU^T B diag(a1Za2Z...amZ00...0) Z = [0, 1; -1, 0]
*
* S S=kUZU^ Z
* | 0, 1, 0 |
* |-1, 0, 0 |
* | 0, 0, 0 |
* Z = diag(1, 1, 0) * W W
* | 0,-1, 0 |
* | 1, 0, 0 |
* | 0, 0, 1 |
* E=SR=Udiag(1,1,0)(WU^R) E
*
* ##
* P=[I0]P E SR
* 1 S=UZU^TUXV^TUXV^Tsvd
* UVE
* Udiag(1,1,0)V^T = E = SR = (UZU^T)(UXV^) = U(ZX)V^T
* ZX = diag(1,1,0) x=W X=W^T
* E SVD Udiag(1,1,0)V^E = SR S = UZU^ R = UWVTor UW^TV^
* St=00t=1线 t = U(0,0,1)^T = u3
* U StR
* P=[UWV^T +u3] or [UWV^T u3] or [UW^TV^T +u3] or [UW^TV^T u3]
* @param E Essential Matrix
* @param R1 Rotation Matrix 1
* @param R2 Rotation Matrix 2
* @param t Translation
* @see Multiple View Geometry in Computer Vision - Result 9.19 p259
*/
void TwoViewReconstruction::DecomposeE(const Eigen::Matrix3f &E, Eigen::Matrix3f &R1, Eigen::Matrix3f &R2, Eigen::Vector3f &t)
{
Eigen::JacobiSVD<Eigen::Matrix3f> svd(E, Eigen::ComputeFullU | Eigen::ComputeFullV);
// 对 t 有归一化但是这个地方并没有决定单目整个SLAM过程的尺度
// 因为CreateInitialMapMonocular函数对3D点深度会缩放然后反过来对 t 有改变
Eigen::Matrix3f U = svd.matrixU();
Eigen::Matrix3f Vt = svd.matrixV().transpose();
@ -918,6 +1211,7 @@ namespace ORB_SLAM3
W(2, 2) = 1;
R1 = U * W * Vt;
// 旋转矩阵有行列式为1的约束
if (R1.determinant() < 0)
R1 = -R1;