文章目录
视觉里程计的目的是根据相邻图像信息估计相机的运动(旋转R、平移t)，主要分为特征点法和直接法，本节将介绍特征点法，直接法下节介绍。
一、特征点法
1.1 特征点
- 特征点=关键点（图像的位置、大小、朝向等）+描述子（关键点周围像素信息）
- 特征点的性质：
1）可重复性：相同特征可在不同图像下找到
2）可区别性：不同特征有不同表达
3）高效性：特征点数远远小于像素数
4）本地性：特征仅与一小片图像区域相关 - 图像特征：SIFT、SURF、ORB
1.2 ORB特征
- FAST角点提取：在原FAST基础上增加了尺度和旋转
1）FAST是一种角点，用于检测局部像素灰度变化明显的地方。
2）以像素点为圆心画一个圆，选取圆上的点，当连续的点与像素点的亮度差距大于特定值时，则可被认为是特征点。
3）FAST不具备方向和尺度，因此ORB考虑添加尺度(构建图像金字塔实现)和旋转(采用灰度质心法) - BEIEF描述子：在原BRIEF基础上考虑了方向
1）BEIEF是一种二进制描述子，其描述向量由0/1构成。
2）两个像素点p和q，当p>q，取1；反之取0。
3）采用BEIEF随即选点，速度快；二进制存储方便；适用实时的图像匹配。
4）原BRIEF描述子不具备旋转不变性，ORB在其基础上考虑了方向。
1.3 特征匹配
- 特征匹配解决了SLAM的数据关联问题，即当前路标和之前路标的对应关系。但存在误匹配现象。
- 度量方法：
1）浮点类型描述子：欧式距离度量
2）二进制描述子：汉明距离度量，指不同位数的个数 - 特征匹配方法：暴力法、快速近似最邻近(FLANN)法
二、特征点的提取与匹配
2.1 使用OpenCV提取匹配ORB
2.1.1 CMakeLists.txt
cmake_minimum_required(VERSION 2.8)
project(vo1)
set(CMAKE_BUILD_TYPE "Release")
add_definitions("-DENABLE_SSE")
set(CMAKE_CXX_FLAGS "-std=c++14 -msse4")
list(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake)
find_package(OpenCV REQUIRED)
find_package(G2O REQUIRED)
find_package(Sophus REQUIRED)
include_directories(
${OpenCV_INCLUDE_DIRS}
${G2O_INCLUDE_DIRS}
${Sophus_INCLUDE_DIRS}
"/usr/include/eigen3/"
)
add_executable(orb_cv orb_cv.cpp)
target_link_libraries(orb_cv ${OpenCV_LIBS} ${FMT_LIBRARIES} fmt)
add_executable(orb_self orb_self.cpp)
target_link_libraries(orb_self ${OpenCV_LIBS} ${FMT_LIBRARIES} fmt)
# add_executable( pose_estimation_2d2d pose_estimation_2d2d.cpp extra.cpp ) # use this if in OpenCV2
add_executable(pose_estimation_2d2d pose_estimation_2d2d.cpp)
target_link_libraries(pose_estimation_2d2d ${OpenCV_LIBS} ${FMT_LIBRARIES} fmt)
# add_executable( triangulation triangulation.cpp extra.cpp) # use this if in opencv2
add_executable(triangulation triangulation.cpp)
target_link_libraries(triangulation ${OpenCV_LIBS} ${FMT_LIBRARIES} fmt)
add_executable(pose_estimation_3d2d pose_estimation_3d2d.cpp)
target_link_libraries(pose_estimation_3d2d
g2o_core g2o_stuff
${OpenCV_LIBS} ${FMT_LIBRARIES} fmt)
add_executable(pose_estimation_3d3d pose_estimation_3d3d.cpp)
target_link_libraries(pose_estimation_3d3d
g2o_core g2o_stuff
${OpenCV_LIBS} ${FMT_LIBRARIES} fmt)
2.1.2 代码展示
#include <iostream>
#include <opencv2/core/core.hpp>
#include <opencv2/features2d/features2d.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <chrono>
using namespace std;
using namespace cv;
int main(int argc, char **argv) {
// 1. 读取图像
Mat img_1 = imread("/home/robot/桌面/slambook2-master/ch7/1.png", cv::IMREAD_COLOR);
Mat img_2 = imread("/home/robot/桌面/slambook2-master/ch7/2.png", cv::IMREAD_COLOR);
assert(img_1.data != nullptr && img_2.data != nullptr);
// 2. 初始化
std::vector<KeyPoint> keypoints_1, keypoints_2;
Mat descriptors_1, descriptors_2;
Ptr<FeatureDetector> detector = ORB::create(); //关键点：里面可以添加参数设置要提取的特征点的个数
Ptr<DescriptorExtractor> descriptor = ORB::create();//描述子
Ptr<DescriptorMatcher> matcher = DescriptorMatcher::create("BruteForce-Hamming");//匹配：这里采用汉明距离
// 3. 第一步:检测 Oriented FAST 角点位置
chrono::steady_clock::time_point t1 = chrono::steady_clock::now();//计时
detector->detect(img_1, keypoints_1);
detector->detect(img_2, keypoints_2);
// 4. 第二步:根据角点位置计算 BRIEF 描述子
descriptor->compute(img_1, keypoints_1, descriptors_1);
descriptor->compute(img_2, keypoints_2, descriptors_2);
chrono::steady_clock::time_point t2 = chrono::steady_clock::now();//计时
chrono::duration<double> time_used = chrono::duration_cast<chrono::duration<double>>(t2 - t1);
cout << "extract ORB cost = " << time_used.count() << " seconds. " << endl;
// 在图中画出特征点的位置
Mat outimg1;
drawKeypoints(img_1, keypoints_1, outimg1, Scalar::all(-1), DrawMatchesFlags::DEFAULT);
imshow("ORB features", outimg1);
// 5. 第三步:对两幅图像中的BRIEF描述子进行匹配，使用 Hamming 距离
vector<DMatch> matches;
t1 = chrono::steady_clock::now();
matcher->match(descriptors_1, descriptors_2, matches);
t2 = chrono::steady_clock::now();
time_used = chrono::duration_cast<chrono::duration<double>>(t2 - t1);
cout << "match ORB cost = " << time_used.count() << " seconds. " << endl;
// 6. 第四步:匹配点对筛选
// 计算最小距离和最大距离
auto min_max = minmax_element(matches.begin(), matches.end(),
[](const DMatch &m1, const DMatch &m2) { return m1.distance < m2.distance; });
double min_dist = min_max.first->distance;
double max_dist = min_max.second->distance;
printf("-- Max dist : %f \n", max_dist);
printf("-- Min dist : %f \n", min_dist);
//当描述子之间的距离大于两倍的最小距离时,即认为匹配有误.但有时候最小距离会非常小,设置一个经验值30作为下限.
std::vector<DMatch> good_matches;
for (int i = 0; i < descriptors_1.rows; i++) {
if (matches[i].distance <= max(2 * min_dist, 30.0)) {
good_matches.push_back(matches[i]);
}
}
// 7. 第五步:绘制匹配结果
Mat img_match;
Mat img_goodmatch;
drawMatches(img_1, keypoints_1, img_2, keypoints_2, matches, img_match);
drawMatches(img_1, keypoints_1, img_2, keypoints_2, good_matches, img_goodmatch);
imshow("all matches", img_match);
imshow("good matches", img_goodmatch);
waitKey(0);
return 0;
}
2.1.3 结果
特征点位置
特征点匹配以及筛选匹配点
2.2 手写ORB特征
2.2.1 CMakeLists.txt
参见2.1.1
2.2.2 代码展示
#include <opencv2/opencv.hpp>
#include <string>
#include <nmmintrin.h>
#include <chrono>
using namespace std;
// global variables
string first_file = "/home/robot/桌面/slambook2-master/ch7/1.png";
string second_file = "/home/robot/桌面/slambook2-master/ch7/2.png";
// 32 bit unsigned int, will have 8, 8x32=256
typedef vector<uint32_t> DescType; // Descriptor type
/** ORB关键点的计算描述子
* compute descriptor of orb keypoints
* @param img input image
* @param keypoints detected fast keypoints
* @param descriptors descriptors
*
* NOTE: if a keypoint goes outside the image boundary (8 pixels), descriptors will not be computed and will be left as
* empty
*/
void ComputeORB(const cv::Mat &img, vector<cv::KeyPoint> &keypoints, vector<DescType> &descriptors);
/** 暴力匹配两组描述子
* brute-force match two sets of descriptors
* @param desc1 the first descriptor
* @param desc2 the second descriptor
* @param matches matches of two images
*/
void BfMatch(const vector<DescType> &desc1, const vector<DescType> &desc2, vector<cv::DMatch> &matches);
int main(int argc, char **argv) {
// 1. 读取图片
cv::Mat first_image = cv::imread(first_file, 0);
cv::Mat second_image = cv::imread(second_file, 0);
assert(first_image.data != nullptr && second_image.data != nullptr);
// 2. 对第一张图片处理：使用阈值=40检测FAST关键点1
chrono::steady_clock::time_point t1 = chrono::steady_clock::now();
vector<cv::KeyPoint> keypoints1;
cv::FAST(first_image, keypoints1, 40); //使用FAST算法检测拐角
vector<DescType> descriptor1;
ComputeORB(first_image, keypoints1, descriptor1);//暴力匹配两组描述子
// 3. 对第二张图片处理
vector<cv::KeyPoint> keypoints2;
vector<DescType> descriptor2;
cv::FAST(second_image, keypoints2, 40); //使用FAST算法检测拐角
ComputeORB(second_image, keypoints2, descriptor2);//暴力匹配两组描述子
chrono::steady_clock::time_point t2 = chrono::steady_clock::now();
chrono::duration<double> time_used = chrono::duration_cast<chrono::duration<double>>(t2 - t1);
cout << "extract ORB cost = " << time_used.count() << " seconds. " << endl;
// 4. 匹配
vector<cv::DMatch> matches;
t1 = chrono::steady_clock::now();
BfMatch(descriptor1, descriptor2, matches);//暴力匹配
t2 = chrono::steady_clock::now();
time_used = chrono::duration_cast<chrono::duration<double>>(t2 - t1);
cout << "match ORB cost = " << time_used.count() << " seconds. " << endl;
cout << "matches: " << matches.size() << endl;
// 5. 画图
cv::Mat image_show;
cv::drawMatches(first_image, keypoints1, second_image, keypoints2, matches, image_show);
cv::imshow("matches", image_show);
cv::imwrite("matches.png", image_show);
cv::waitKey(0);
cout << "done." << endl;
return 0;
}
// -------------------------------------------------------------------------------------------------- //
// ORB pattern
int ORB_pattern[256 * 4] = {
8, -3, 9, 5/*mean (0), correlation (0)*/,
4, 2, 7, -12/*mean (1.12461e-05), correlation (0.0437584)*/,
-11, 9, -8, 2/*mean (3.37382e-05), correlation (0.0617409)*/,
7, -12, 12, -13/*mean (5.62303e-05), correlation (0.0636977)*/,
2, -13, 2, 12/*mean (0.000134953), correlation (0.085099)*/,
1, -7, 1, 6/*mean (0.000528565), correlation (0.0857175)*/,
-2, -10, -2, -4/*mean (0.0188821), correlation (0.0985774)*/,
-13, -13, -11, -8/*mean (0.0363135), correlation (0.0899616)*/,
-13, -3, -12, -9/*mean (0.121806), correlation (0.099849)*/,
10, 4, 11, 9/*mean (0.122065), correlation (0.093285)*/,
-13, -8, -8, -9/*mean (0.162787), correlation (0.0942748)*/,
-11, 7, -9, 12/*mean (0.21561), correlation (0.0974438)*/,
7, 7, 12, 6/*mean (0.160583), correlation (0.130064)*/,
-4, -5, -3, 0/*mean (0.228171), correlation (0.132998)*/,
-13, 2, -12, -3/*mean (0.00997526), correlation (0.145926)*/,
-9, 0, -7, 5/*mean (0.198234), correlation (0.143636)*/,
12, -6, 12, -1/*mean (0.0676226), correlation (0.16689)*/,
-3, 6, -2, 12/*mean (0.166847), correlation (0.171682)*/,
-6, -13, -4, -8/*mean (0.101215), correlation (0.179716)*/,
11, -13, 12, -8/*mean (0.200641), correlation (0.192279)*/,
4, 7, 5, 1/*mean (0.205106), correlation (0.186848)*/,
5, -3, 10, -3/*mean (0.234908), correlation (0.192319)*/,
3, -7, 6, 12/*mean (0.0709964), correlation (0.210872)*/,
-8, -7, -6, -2/*mean (0.0939834), correlation (0.212589)*/,
-2, 11, -1, -10/*mean (0.127778), correlation (0.20866)*/,
-13, 12, -8, 10/*mean (0.14783), correlation (0.206356)*/,
-7, 3, -5, -3/*mean (0.182141), correlation (0.198942)*/,
-4, 2, -3, 7/*mean (0.188237), correlation (0.21384)*/,
-10, -12, -6, 11/*mean (0.14865), correlation (0.23571)*/,
5, -12, 6, -7/*mean (0.222312), correlation (0.23324)*/,
5, -6, 7, -1/*mean (0.229082), correlation (0.23389)*/,
1, 0, 4, -5/*mean (0.241577), correlation (0.215286)*/,
9, 11, 11, -13/*mean (0.00338507), correlation (0.251373)*/,
4, 7, 4, 12/*mean (0.131005), correlation (0.257622)*/,
2, -1, 4, 4/*mean (0.152755), correlation (0.255205)*/,
-4, -12, -2, 7/*mean (0.182771), correlation (0.244867)*/,
-8, -5, -7, -10/*mean (0.186898), correlation (0.23901)*/,
4, 11, 9, 12/*mean (0.226226), correlation (0.258255)*/,
0, -8, 1, -13/*mean (0.0897886), correlation (0.274827)*/,
-13, -2, -8, 2/*mean (0.148774), correlation (0.28065)*/,
-3, -2, -2, 3/*mean (0.153048), correlation (0.283063)*/,
-6, 9, -4, -9/*mean (0.169523), correlation (0.278248)*/,
8, 12, 10, 7/*mean (0.225337), correlation (0.282851)*/,
0, 9, 1, 3/*mean (0.226687), correlation (0.278734)*/,
7, -5, 11, -10/*mean (0.00693882), correlation (0.305161)*/,
-13, -6, -11, 0/*mean (0.0227283), correlation (0.300181)*/,
10, 7, 12, 1/*mean (0.125517), correlation (0.31089)*/,
-6, -3, -6, 12/*mean (0.131748), correlation (0.312779)*/,
10, -9, 12, -4/*mean (0.144827), correlation (0.292797)*/,
-13, 8, -8, -12/*mean (0.149202), correlation (0.308918)*/,
-13, 0, -8, -4/*mean (0.160909), correlation (0.310013)*/,
3, 3, 7, 8/*mean (0.177755), correlation (0.309394)*/,
5, 7, 10, -7/*mean (0.212337), correlation (0.310315)*/,
-1, 7, 1, -12/*mean (0.214429), correlation (0.311933)*/,
3, -10, 5, 6/*mean (0.235807), correlation (0.313104)*/,
2, -4, 3, -10/*mean (0.00494827), correlation (0.344948)*/,
-13, 0, -13, 5/*mean (0.0549145), correlation (0.344675)*/,
-13, -7, -12, 12/*mean (0.103385), correlation (0.342715)*/,
-13, 3, -11, 8/*mean (0.134222), correlation (0.322922)*/,
-7, 12, -4, 7/*mean (0.153284), correlation (0.337061)*/,
6, -10, 12, 8/*mean (0.154881), correlation (0.329257)*/,
-9, -1, -7, -6/*mean (0.200967), correlation (0.33312)*/,
-2, -5, 0, 12/*mean (0.201518), correlation (0.340635)*/,
-12, 5, -7, 5/*mean (0.207805), correlation (0.335631)*/,
3, -10, 8, -13/*mean (0.224438), correlation (0.34504)*/,
-7, -7, -4, 5/*mean (0.239361), correlation (0.338053)*/,
-3, -2, -1, -7/*mean (0.240744), correlation (0.344322)*/,
2, 9, 5, -11/*mean (0.242949), correlation (0.34145)*/,
-11, -13, -5, -13/*mean (0.244028), correlation (0.336861)*/,
-1, 6, 0, -1/*mean (0.247571), correlation (0.343684)*/,
5, -3, 5, 2/*mean (0.000697256), correlation (0.357265)*/,
-4, -13, -4, 12/*mean (0.00213675), correlation (0.373827)*/,
-9, -6, -9, 6/*mean (0.0126856), correlation (0.373938)*/,
-12, -10, -8, -4/*mean (0.0152497), correlation (0.364237)*/,
10, 2, 12, -3/*mean (0.0299933), correlation (0.345292)*/,
7, 12, 12, 12/*mean (0.0307242), correlation (0.366299)*/,
-7, -13, -6, 5/*mean (0.0534975), correlation (0.368357)*/,
-4, 9, -3, 4/*mean (0.099865), correlation (0.372276)*/,
7, -1, 12, 2/*mean (0.117083), correlation (0.364529)*/,
-7, 6, -5, 1/*mean (0.126125), correlation (0.369606)*/,
-13, 11, -12, 5/*mean (0.130364), correlation (0.358502)*/,
-3, 7, -2, -6/*mean (0.131691), correlation (0.375531)*/,
7, -8, 12, -7/*mean (0.160166), correlation (0.379508)*/,
-13, -7, -11, -12/*mean (0.167848), correlation (0.353343)*/,
1, -3, 12, 12/*mean (0.183378), correlation (0.371916)*/,
2, -6, 3, 0/*mean (0.228711), correlation (0.371761)*/,
-4, 3, -2, -13/*mean (0.247211), correlation (0.364063)*/,
-1, -13, 1, 9/*mean (0.249325), correlation (0.378139)*/,
7, 1, 8, -6/*mean (0.000652272), correlation (0.411682)*/,
1, -1, 3, 12/*mean (0.00248538), correlation (0.392988)*/,
9, 1, 12, 6/*mean (0.0206815), correlation (0.386106)*/,
-1, -9, -1, 3/*mean (0.0364485), correlation (0.410752)*/,
-13, -13, -10, 5/*mean (0.0376068), correlation (0.398374)*/,
7, 7, 10, 12/*mean (0.0424202), correlation (0.405663)*/,
12, -5, 12, 9/*mean (0.0942645), correlation (0.410422)*/,
6, 3, 7, 11/*mean (0.1074), correlation (0.413224)*/,
5, -13, 6, 10/*mean (0.109256), correlation (0.408646)*/,
2, -12, 2, 3/*mean (0.131691), correlation (0.416076)*/,
3, 8, 4, -6/*mean (0.165081), correlation (0.417569)*/,
2, 6, 12, -13/*mean (0.171874), correlation (0.408471)*/,
9, -12, 10, 3/*mean (0.175146), correlation (0.41296)*/,
-8, 4, -7, 9/*mean (0.183682), correlation (0.402956)*/,
-11, 12, -4, -6/*mean (0.184672), correlation (0.416125)*/,
1, 12, 2, -8/*mean (0.191487), correlation (0.386696)*/,
6, -9, 7, -4/*mean (0.192668), correlation (0.394771)*/,
2, 3, 3, -2/*mean (0.200157), correlation (0.408303)*/,
6, 3, 11, 0/*mean (0.204588), correlation (0.411762)*/,
3, -3, 8, -8/*mean (0.205904), correlation (0.416294)*/,
7, 8, 9, 3/*mean (0.213237), correlation (0.409306)*/,
-11, -5, -6, -4/*mean (0.243444), correlation (0.395069)*/,
-10, 11, -5, 10/*mean (0.247672), correlation (0.413392)*/,
-5, -8, -3, 12/*mean (0.24774), correlation (0.411416)*/,
-10, 5, -9, 0/*mean (0.00213675), correlation (0.454003)*/,
8, -1, 12, -6/*mean (0.0293635), correlation (0.455368)*/,
4, -6, 6, -11/*mean (0.0404971), correlation (0.457393)*/,
-10, 12, -8, 7/*mean (0.0481107), correlation (0.448364)*/,
4, -2, 6, 7/*mean (0.050641), correlation (0.455019)*/,
-2, 0, -2, 12/*mean (0.0525978), correlation (0.44338)*/,
-5, -8, -5, 2/*mean (0.0629667), correlation (0.457096)*/,
7, -6, 10, 12/*mean (0.0653846), correlation (0.445623)*/,
-9, -13, -8, -8/*mean (0.0858749), correlation (0.449789)*/,
-5, -13, -5, -2/*mean (0.122402), correlation (0.450201)*/,
8, -8, 9, -13/*mean (0.125416), correlation (0.453224)*/,
-9, -11, -9, 0/*mean (0.130128), correlation (0.458724)*/,
1, -8, 1, -2/*mean (0.132467), correlation (0.440133)*/,
7, -4, 9, 1/*mean (0.132692), correlation (0.454)*/,
-2, 1, -1, -4/*mean (0.135695), correlation (0.455739)*/,
11, -6, 12, -11/*mean (0.142904), correlation (0.446114)*/,
-12, -9, -6, 4/*mean (0.146165), correlation (0.451473)*/,
3, 7, 7, 12/*mean (0.147627), correlation (0.456643)*/,
5, 5, 10, 8/*mean (0.152901), correlation (0.455036)*/,
0, -4, 2, 8/*mean (0.167083), correlation (0.459315)*/,
-9, 12, -5, -13/*mean (0.173234), correlation (0.454706)*/,
0, 7, 2, 12/*mean (0.18312), correlation (0.433855)*/,
-1, 2, 1, 7/*mean (0.185504), correlation (0.443838)*/,
5, 11, 7, -9/*mean (0.185706), correlation (0.451123)*/,
3, 5, 6, -8/*mean (0.188968), correlation (0.455808)*/,
-13, -4, -8, 9/*mean (0.191667), correlation (0.459128)*/,
-5, 9, -3, -3/*mean (0.193196), correlation (0.458364)*/,
-4, -7, -3, -12/*mean (0.196536), correlation (0.455782)*/,
6, 5, 8, 0/*mean (0.1972), correlation (0.450481)*/,
-7, 6, -6, 12/*mean (0.199438), correlation (0.458156)*/,
-13, 6, -5, -2/*mean (0.211224), correlation (0.449548)*/,
1, -10, 3, 10/*mean (0.211718), correlation (0.440606)*/,
4, 1, 8, -4/*mean (0.213034), correlation (0.443177)*/,
-2, -2, 2, -13/*mean (0.234334), correlation (0.455304)*/,
2, -12, 12, 12/*mean (0.235684), correlation (0.443436)*/,
-2, -13, 0, -6/*mean (0.237674), correlation (0.452525)*/,
4, 1, 9, 3/*mean (0.23962), correlation (0.444824)*/,
-6, -10, -3, -5/*mean (0.248459), correlation (0.439621)*/,
-3, -13, -1, 1/*mean (0.249505), correlation (0.456666)*/,
7, 5, 12, -11/*mean (0.00119208), correlation (0.495466)*/,
4, -2, 5, -7/*mean (0.00372245), correlation (0.484214)*/,
-13, 9, -9, -5/*mean (0.00741116), correlation (0.499854)*/,
7, 1, 8, 6/*mean (0.0208952), correlation (0.499773)*/,
7, -8, 7, 6/*mean (0.0220085), correlation (0.501609)*/,
-7, -4, -7, 1/*mean (0.0233806), correlation (0.496568)*/,
-8, 11, -7, -8/*mean (0.0236505), correlation (0.489719)*/,
-13, 6, -12, -8/*mean (0.0268781), correlation (0.503487)*/,
2, 4, 3, 9/*mean (0.0323324), correlation (0.501938)*/,
10, -5, 12, 3/*mean (0.0399235), correlation (0.494029)*/,
-6, -5, -6, 7/*mean (0.0420153), correlation (0.486579)*/,
8, -3, 9, -8/*mean (0.0548021), correlation (0.484237)*/,
2, -12, 2, 8/*mean (0.0616622), correlation (0.496642)*/,
-11, -2, -10, 3/*mean (0.0627755), correlation (0.498563)*/,
-12, -13, -7, -9/*mean (0.0829622), correlation (0.495491)*/,
-11, 0, -10, -5/*mean (0.0843342), correlation (0.487146)*/,
5, -3, 11, 8/*mean (0.0929937), correlation (0.502315)*/,
-2, -13, -1, 12/*mean (0.113327), correlation (0.48941)*/,
-1, -8, 0, 9/*mean (0.132119), correlation (0.467268)*/,
-13, -11, -12, -5/*mean (0.136269), correlation (0.498771)*/,
-10, -2, -10, 11/*mean (0.142173), correlation (0.498714)*/,
-3, 9, -2, -13/*mean (0.144141), correlation (0.491973)*/,
2, -3, 3, 2/*mean (0.14892), correlation (0.500782)*/,
-9, -13, -4, 0/*mean (0.150371), correlation (0.498211)*/,
-4, 6, -3, -10/*mean (0.152159), correlation (0.495547)*/,
-4, 12, -2, -7/*mean (0.156152), correlation (0.496925)*/,
-6, -11, -4, 9/*mean (0.15749), correlation (0.499222)*/,
6, -3, 6, 11/*mean (0.159211), correlation (0.503821)*/,
-13, 11, -5, 5/*mean (0.162427), correlation (0.501907)*/,
11, 11, 12, 6/*mean (0.16652), correlation (0.497632)*/,
7, -5, 12, -2/*mean (0.169141), correlation (0.484474)*/,
-1, 12, 0, 7/*mean (0.169456), correlation (0.495339)*/,
-4, -8, -3, -2/*mean (0.171457), correlation (0.487251)*/,
-7, 1, -6, 7/*mean (0.175), correlation (0.500024)*/,
-13, -12, -8, -13/*mean (0.175866), correlation (0.497523)*/,
-7, -2, -6, -8/*mean (0.178273), correlation (0.501854)*/,
-8, 5, -6, -9/*mean (0.181107), correlation (0.494888)*/,
-5, -1, -4, 5/*mean (0.190227), correlation (0.482557)*/,
-13, 7, -8, 10/*mean (0.196739), correlation (0.496503)*/,
1, 5, 5, -13/*mean (0.19973), correlation (0.499759)*/,
1, 0, 10, -13/*mean (0.204465), correlation (0.49873)*/,
9, 12, 10, -1/*mean (0.209334), correlation (0.49063)*/,
5, -8, 10, -9/*mean (0.211134), correlation (0.503011)*/,
-1, 11, 1, -13/*mean (0.212), correlation (0.499414)*/,
-9, -3, -6, 2/*mean (0.212168), correlation (0.480739)*/,
-1, -10, 1, 12/*mean (0.212731), correlation (0.502523)*/,
-13, 1, -8, -10/*mean (0.21327), correlation (0.489786)*/,
8, -11, 10, -6/*mean (0.214159), correlation (0.488246)*/,
2, -13, 3, -6/*mean (0.216993), correlation (0.50287)*/,
7, -13, 12, -9/*mean (0.223639), correlation (0.470502)*/,
-10, -10, -5, -7/*mean (0.224089), correlation (0.500852)*/,
-10, -8, -8, -13/*mean (0.228666), correlation (0.502629)*/,
4, -6, 8, 5/*mean (0.22906), correlation (0.498305)*/,
3, 12, 8, -13/*mean (0.233378), correlation (0.503825)*/,
-4, 2, -3, -3/*mean (0.234323), correlation (0.476692)*/,
5, -13, 10, -12/*mean (0.236392), correlation (0.475462)*/,
4, -13, 5, -1/*mean (0.236842), correlation (0.504132)*/,
-9, 9, -4, 3/*mean (0.236977), correlation (0.497739)*/,
0, 3, 3, -9/*mean (0.24314), correlation (0.499398)*/,
-12, 1, -6, 1/*mean (0.243297), correlation (0.489447)*/,
3, 2, 4, -8/*mean (0.00155196), correlation (0.553496)*/,
-10, -10, -10, 9/*mean (0.00239541), correlation (0.54297)*/,
8, -13, 12, 12/*mean (0.0034413), correlation (0.544361)*/,
-8, -12, -6, -5/*mean (0.003565), correlation (0.551225)*/,
2, 2, 3, 7/*mean (0.00835583), correlation (0.55285)*/,
10, 6, 11, -8/*mean (0.00885065), correlation (0.540913)*/,
6, 8, 8, -12/*mean (0.0101552), correlation (0.551085)*/,
-7, 10, -6, 5/*mean (0.0102227), correlation (0.533635)*/,
-3, -9, -3, 9/*mean (0.0110211), correlation (0.543121)*/,
-1, -13, -1, 5/*mean (0.0113473), correlation (0.550173)*/,
-3, -7, -3, 4/*mean (0.0140913), correlation (0.554774)*/,
-8, -2, -8, 3/*mean (0.017049), correlation (0.55461)*/,
4, 2, 12, 12/*mean (0.01778), correlation (0.546921)*/,
2, -5, 3, 11/*mean (0.0224022), correlation (0.549667)*/,
6, -9, 11, -13/*mean (0.029161), correlation (0.546295)*/,
3, -1, 7, 12/*mean (0.0303081), correlation (0.548599)*/,
11, -1, 12, 4/*mean (0.0355151), correlation (0.523943)*/,
-3, 0, -3, 6/*mean (0.0417904), correlation (0.543395)*/,
4, -11, 4, 12/*mean (0.0487292), correlation (0.542818)*/,
2, -4, 2, 1/*mean (0.0575124), correlation (0.554888)*/,
-10, -6, -8, 1/*mean (0.0594242), correlation (0.544026)*/,
-13, 7, -11, 1/*mean (0.0597391), correlation (0.550524)*/,
-13, 12, -11, -13/*mean (0.0608974), correlation (0.55383)*/,
6, 0, 11, -13/*mean (0.065126), correlation (0.552006)*/,
0, -1, 1, 4/*mean (0.074224), correlation (0.546372)*/,
-13, 3, -9, -2/*mean (0.0808592), correlation (0.554875)*/,
-9, 8, -6, -3/*mean (0.0883378), correlation (0.551178)*/,
-13, -6, -8, -2/*mean (0.0901035), correlation (0.548446)*/,
5, -9, 8, 10/*mean (0.0949843), correlation (0.554694)*/,
2, 7, 3, -9/*mean (0.0994152), correlation (0.550979)*/,
-1, -6, -1, -1/*mean (0.10045), correlation (0.552714)*/,
9, 5, 11, -2/*mean (0.100686), correlation (0.552594)*/,
11, -3, 12, -8/*mean (0.101091), correlation (0.532394)*/,
3, 0, 3, 5/*mean (0.101147), correlation (0.525576)*/,
-1, 4, 0, 10/*mean (0.105263), correlation (0.531498)*/,
3, -6, 4, 5/*mean (0.110785), correlation (0.540491)*/,
-13, 0, -10, 5/*mean (0.112798), correlation (0.536582)*/,
5, 8, 12, 11/*mean (0.114181), correlation (0.555793)*/,
8, 9, 9, -6/*mean (0.117431), correlation (0.553763)*/,
7, -4, 8, -12/*mean (0.118522), correlation (0.553452)*/,
-10, 4, -10, 9/*mean (0.12094), correlation (0.554785)*/,
7, 3, 12, 4/*mean (0.122582), correlation (0.555825)*/,
9, -7, 10, -2/*mean (0.124978), correlation (0.549846)*/,
7, 0, 12, -2/*mean (0.127002), correlation (0.537452)*/,
-1, -6, 0, -11/*mean (0.127148), correlation (0.547401)*/
};
// 计算描述子compute the descriptor
void ComputeORB(const cv::Mat &img, vector<cv::KeyPoint> &keypoints, vector<DescType> &descriptors) {
const int half_patch_size = 8;
const int half_boundary = 16;
int bad_points = 0;
for (auto &kp: keypoints) {
if (kp.pt.x < half_boundary || kp.pt.y < half_boundary ||
kp.pt.x >= img.cols - half_boundary || kp.pt.y >= img.rows - half_boundary) {
// outside
bad_points++;
descriptors.push_back({});
continue;
}
float m01 = 0, m10 = 0;
for (int dx = -half_patch_size; dx < half_patch_size; ++dx) {
for (int dy = -half_patch_size; dy < half_patch_size; ++dy) {
uchar pixel = img.at<uchar>(kp.pt.y + dy, kp.pt.x + dx);
m10 += dx * pixel;
m01 += dy * pixel;
}
}
// angle should be arc tan(m01/m10);
float m_sqrt = sqrt(m01 * m01 + m10 * m10) + 1e-18; // avoid divide by zero
float sin_theta = m01 / m_sqrt;
float cos_theta = m10 / m_sqrt;
// compute the angle of this point
DescType desc(8, 0);
for (int i = 0; i < 8; i++) {
uint32_t d = 0;
for (int k = 0; k < 32; k++) {
int idx_pq = i * 32 + k;
cv::Point2f p(ORB_pattern[idx_pq * 4], ORB_pattern[idx_pq * 4 + 1]);
cv::Point2f q(ORB_pattern[idx_pq * 4 + 2], ORB_pattern[idx_pq * 4 + 3]);
// rotate with theta
cv::Point2f pp = cv::Point2f(cos_theta * p.x - sin_theta * p.y, sin_theta * p.x + cos_theta * p.y)
+ kp.pt;
cv::Point2f qq = cv::Point2f(cos_theta * q.x - sin_theta * q.y, sin_theta * q.x + cos_theta * q.y)
+ kp.pt;
if(pp.y < 0)
pp.y=0;
if (img.at<uchar>(pp.y, pp.x) < img.at<uchar>(qq.y, qq.x)) {
d |= 1 << k;
}
}
desc[i] = d;
}
descriptors.push_back(desc);
}
cout << "bad/total: " << bad_points << "/" << keypoints.size() << endl;
}
// 暴力匹配brute-force matching
void BfMatch(const vector<DescType> &desc1, const vector<DescType> &desc2, vector<cv::DMatch> &matches) {
const int d_max = 40;
// 计算每个特征点与所有特征点的距离，取最近
for (size_t i1 = 0; i1 < desc1.size(); ++i1) {
if (desc1[i1].empty()) continue;
cv::DMatch m{i1, 0, 256};
for (size_t i2 = 0; i2 < desc2.size(); ++i2) {
if (desc2[i2].empty()) continue;
int distance = 0;
for (int k = 0; k < 8; k++) {
distance += _mm_popcnt_u32(desc1[i1][k] ^ desc2[i2][k]);
}
if (distance < d_max && distance < m.distance) {
m.distance = distance;
m.trainIdx = i2;
}
}
if (m.distance < d_max) {
matches.push_back(m);
}
}
}
2.2.3 结果
三、根据特征点估计相机运动：旋转R 、平移t
估计相机运动的三种方法：
四、对极几何法（2D->2D）：已知二维特征点(像素坐标)估计相机位姿
4.1 理论部分
- 对极约束：
- 特征点不共面时，采用本质矩阵E或者基础矩阵F(八点法)，然后用H恢复旋转R和平移t（采用奇异值分解SVD可得到4个分解）
- 特征点共面时，采用单应矩阵H，然后用H恢复旋转R和平移t
- 关于八点法的讨论：
1）八点法应用于单目SLAM的初始化，初始化后2D->2D问题就转化成3D->2D，之后就可用PnP求解
2）具有尺度不确定性
3）不能纯旋转，还要有平移
4）多于八点时，使用最小二乘或RANSAC
4.2 CMakeLists.txt
参见2.1.1
4.3 代码展示
#include <iostream>
#include <opencv2/core/core.hpp>
#include <opencv2/features2d/features2d.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/calib3d/calib3d.hpp>
// #include "extra.h" // use this if in OpenCV2
#include "opencv2/imgcodecs/legacy/constants_c.h"
using namespace std;
using namespace cv;
/****************************************************
* 本程序演示了如何使用2D-2D的特征匹配估计相机运动
* **************************************************/
void find_feature_matches(
const Mat &img_1, const Mat &img_2,
std::vector<KeyPoint> &keypoints_1,
std::vector<KeyPoint> &keypoints_2,
std::vector<DMatch> &matches);
void pose_estimation_2d2d(
std::vector<KeyPoint> keypoints_1,
std::vector<KeyPoint> keypoints_2,
std::vector<DMatch> matches,
Mat &R, Mat &t);
// 像素坐标转相机归一化坐标
Point2d pixel2cam(const Point2d &p, const Mat &K);
int main(int argc, char **argv) {
// 1.读取图像
Mat img_1 = imread("/home/robot/桌面/slambook2-master/ch7/1.png", cv::IMREAD_COLOR);
Mat img_2 = imread("/home/robot/桌面/slambook2-master/ch7/2.png", cv::IMREAD_COLOR);
assert(img_1.data && img_2.data && "Can not load images!");
// 2. 特征点的提取与匹配
vector<KeyPoint> keypoints_1, keypoints_2;
vector<DMatch> matches;
find_feature_matches(img_1, img_2, keypoints_1, keypoints_2, matches);
cout << "一共找到了" << matches.size() << "组匹配点" << endl;
// 3. 估计两张图像间运动(2D->2D 对极几何)：根据特征点求出：旋转R、平移t
Mat R, t;
pose_estimation_2d2d(keypoints_1, keypoints_2, matches, R, t);
// 4. 验证E=t^R*scale
Mat t_x =
(Mat_<double>(3, 3) <<
0, -t.at<double>(2, 0), t.at<double>(1, 0),
t.at<double>(2, 0), 0, -t.at<double>(0, 0),
-t.at<double>(1, 0), t.at<double>(0, 0), 0);// t变成t^
cout << "t^R=" << endl << t_x * R << endl;// E=t^R
// 5. 验证对极约束
Mat K = (Mat_<double>(3, 3) << 520.9, 0, 325.1, 0, 521.0, 249.7, 0, 0, 1);
for (DMatch m: matches) {
Point2d pt1 = pixel2cam(keypoints_1[m.queryIdx].pt, K);// 像素坐标x1➗内参K->归一化坐标y1
Mat y1 = (Mat_<double>(3, 1) << pt1.x, pt1.y, 1);
Point2d pt2 = pixel2cam(keypoints_2[m.trainIdx].pt, K);// 像素坐标x2➗内参K->归一化坐标y2
Mat y2 = (Mat_<double>(3, 1) << pt2.x, pt2.y, 1);
Mat d = y2.t() * t_x * R * y1; // 验证对极约束y2^ * t^ * R * y1 约等于 0
cout << "epipolar constraint = " << d << endl;
}
return 0;
}
void find_feature_matches(const Mat &img_1, const Mat &img_2,
std::vector<KeyPoint> &keypoints_1,
std::vector<KeyPoint> &keypoints_2,
std::vector<DMatch> &matches) {
//-- 初始化
Mat descriptors_1, descriptors_2;
// used in OpenCV3
Ptr<FeatureDetector> detector = ORB::create();
Ptr<DescriptorExtractor> descriptor = ORB::create();
// use this if you are in OpenCV2
// Ptr<FeatureDetector> detector = FeatureDetector::create ( "ORB" );
// Ptr<DescriptorExtractor> descriptor = DescriptorExtractor::create ( "ORB" );
Ptr<DescriptorMatcher> matcher = DescriptorMatcher::create("BruteForce-Hamming");
//-- 第一步:检测 Oriented FAST 角点位置
detector->detect(img_1, keypoints_1);
detector->detect(img_2, keypoints_2);
//-- 第二步:根据角点位置计算 BRIEF 描述子
descriptor->compute(img_1, keypoints_1, descriptors_1);
descriptor->compute(img_2, keypoints_2, descriptors_2);
//-- 第三步:对两幅图像中的BRIEF描述子进行匹配，使用 Hamming 距离
vector<DMatch> match;
//BFMatcher matcher ( NORM_HAMMING );
matcher->match(descriptors_1, descriptors_2, match);
//-- 第四步:匹配点对筛选
double min_dist = 10000, max_dist = 0;
//找出所有匹配之间的最小距离和最大距离, 即是最相似的和最不相似的两组点之间的距离
for (int i = 0; i < descriptors_1.rows; i++) {
double dist = match[i].distance;
if (dist < min_dist) min_dist = dist;
if (dist > max_dist) max_dist = dist;
}
printf("-- Max dist : %f \n", max_dist);
printf("-- Min dist : %f \n", min_dist);
//当描述子之间的距离大于两倍的最小距离时,即认为匹配有误.但有时候最小距离会非常小,设置一个经验值30作为下限.
for (int i = 0; i < descriptors_1.rows; i++) {
if (match[i].distance <= max(2 * min_dist, 30.0)) {
matches.push_back(match[i]);
}
}
}
// 课本P99
Point2d pixel2cam(const Point2d &p, const Mat &K) {
return Point2d
(
(p.x - K.at<double>(0, 2)) / K.at<double>(0, 0),
(p.y - K.at<double>(1, 2)) / K.at<double>(1, 1)
);
}
void pose_estimation_2d2d(std::vector<KeyPoint> keypoints_1,
std::vector<KeyPoint> keypoints_2,
std::vector<DMatch> matches,
Mat &R, Mat &t) {
// 相机内参,TUM Freiburg2
Mat K = (Mat_<double>(3, 3) << 520.9, 0, 325.1, 0, 521.0, 249.7, 0, 0, 1);
//-- 把匹配点转换为vector<Point2f>的形式
vector<Point2f> points1;
vector<Point2f> points2;
for (int i = 0; i < (int) matches.size(); i++) {
points1.push_back(keypoints_1[matches[i].queryIdx].pt);
points2.push_back(keypoints_2[matches[i].trainIdx].pt);
}
//-- 计算基础矩阵
Mat fundamental_matrix;
fundamental_matrix = findFundamentalMat(points1, points2, FM_8POINT);
cout << "fundamental_matrix is " << endl << fundamental_matrix << endl;
//-- 计算本质矩阵
Point2d principal_point(325.1, 249.7); //相机光心, TUM dataset标定值
double focal_length = 521; //相机焦距, TUM dataset标定值
Mat essential_matrix;
essential_matrix = findEssentialMat(points1, points2, focal_length, principal_point);
cout << "essential_matrix is " << endl << essential_matrix << endl;
//-- 计算单应矩阵
//-- 但是本例中场景不是平面，单应矩阵意义不大
Mat homography_matrix;
homography_matrix = findHomography(points1, points2, RANSAC, 3);
cout << "homography_matrix is " << endl << homography_matrix << endl;
//-- 从本质矩阵中恢复旋转和平移信息.
// 此函数仅在Opencv3中提供
recoverPose(essential_matrix, points1, points2, R, t, focal_length, principal_point);
cout << "R is " << endl << R << endl;
cout << "t is " << endl << t << endl;
}
4.4 结果
-- Max dist : 94.000000
-- Min dist : 4.000000
一共找到了79组匹配点
fundamental_matrix is
[4.54443750398184e-06, 0.0001333855576992603, -0.01798499246479044;
-0.0001275657012964255, 2.266794804645652e-05, -0.01416678429206633;
0.01814994639971766, 0.004146055870980492, 1]
essential_matrix is
[-0.008455114492964278, 0.05451570701059781, 0.1546375809484052;
-0.008287154708445212, 0.03351311565984172, -0.6896472136971504;
-0.1153993974485718, 0.6945899967012867, 0.02159624094256633]
homography_matrix is
[0.9261214237658335, -0.1445322040802305, 33.26921164265664;
0.04535424230636757, 0.9386696658342905, 8.570980713233848;
-1.006198269424755e-05, -3.008140685985328e-05, 1]
R is
[0.9956584940813579, -0.05615340406690447, 0.07423582945816433;
0.05268846331440004, 0.9974645001566195, 0.04783823534446425;
-0.07673388428334535, -0.0437191735855581, 0.9960926386957119]
t is
[-0.9726703113454949;
-0.2153829834753195;
0.08673313009645391]
t^R=
[0.01195733758736675, -0.07709685221674556, -0.2186905642298021;
0.01171980658216709, -0.04739470268352609, 0.9753084428633267;
0.1631993929614534, -0.9822985936236425, -0.03054169683725466]
epipolar constraint = [-0.0005617285518606241]
epipolar constraint = [0.002891683190146016]
epipolar constraint = [-0.0001941259398173245]
epipolar constraint = [0.003462947761727536]
......//太多省略部分
epipolar constraint = [0.005653889777384447]
epipolar constraint = [0.0008830143247820065]
epipolar constraint = [-0.001103292290051336]
epipolar constraint = [-0.003982708195313309]
epipolar constraint = [-0.0053874915375101]
五、三角测量：根据相机运动估计特征点世界坐标
5.1 理论部分
- 三角测量：通过不同位置对同一路标观察，推断路标点的距离
- 通过R和t推出特征点在两个相机下的深度s1,s2
- 三角测量是通过平移得到的，纯旋转不行
- 平移小，不确定性大；平移大，三角测量将更加精确
- 三角化的矛盾/视差：要想提高三角化的精度？
1）提高特征点的精度，也就是提高图像分辨率，将导致图像变大，增加计算成本
2）平移量增大，导致图像外观发生明显变化
5.2 CMakeLists.txt
参见2.1.1
5.3 代码展示
#include <iostream>
#include <opencv2/opencv.hpp>
// #include "extra.h" // used in opencv2
using namespace std;
using namespace cv;
void find_feature_matches(
const Mat &img_1, const Mat &img_2,
std::vector<KeyPoint> &keypoints_1,
std::vector<KeyPoint> &keypoints_2,
std::vector<DMatch> &matches);
void pose_estimation_2d2d(
const std::vector<KeyPoint> &keypoints_1,
const std::vector<KeyPoint> &keypoints_2,
const std::vector<DMatch> &matches,
Mat &R, Mat &t);
void triangulation(
const vector<KeyPoint> &keypoint_1,
const vector<KeyPoint> &keypoint_2,
const std::vector<DMatch> &matches,
const Mat &R, const Mat &t,
vector<Point3d> &points
);
/// 作图用
inline cv::Scalar get_color(float depth) {
float up_th = 50, low_th = 10, th_range = up_th - low_th;
if (depth > up_th) depth = up_th;
if (depth < low_th) depth = low_th;
return cv::Scalar(255 * depth / th_range, 0, 255 * (1 - depth / th_range));
}
// 像素坐标转相机归一化坐标
Point2f pixel2cam(const Point2d &p, const Mat &K);
int main(int argc, char **argv) {
// 1. 读取图像
Mat img_1 = imread("/home/robot/桌面/slambook2-master/ch7/1.png", cv::IMREAD_COLOR);
Mat img_2 = imread("/home/robot/桌面/slambook2-master/ch7/2.png", cv::IMREAD_COLOR);
// 2. 特征点提取与匹配
vector<KeyPoint> keypoints_1, keypoints_2;
vector<DMatch> matches;//声明匹配点：经过筛选得到匹配点
find_feature_matches(img_1, img_2, keypoints_1, keypoints_2, matches);
cout << "一共找到了" << matches.size() << "组匹配点" << endl;
// 3. 估计两张图像间运动(2D->2D 对极几何)：根据特征点求出：旋转R、平移t
Mat R, t;
pose_estimation_2d2d(keypoints_1, keypoints_2, matches, R, t);
// 4. 三角化：根据R t估计特征点points的空间位置
vector<Point3d> points;
triangulation(keypoints_1, keypoints_2, matches, R, t, points);
// 5. 验证三角化点与特征点的重投影关系,并且将特征点在图中画出
Mat K = (Mat_<double>(3, 3) << 520.9, 0, 325.1, 0, 521.0, 249.7, 0, 0, 1);//内参
Mat img1_plot = img_1.clone();
Mat img2_plot = img_2.clone();
for (int i = 0; i < matches.size(); i++) {
// 第一个图
float depth1 = points[i].z;
cout << "depth: " << depth1 << endl;
Point2d pt1_cam = pixel2cam(keypoints_1[matches[i].queryIdx].pt, K);//像素坐标->归一化坐标
cv::circle(img1_plot, keypoints_1[matches[i].queryIdx].pt, 2, get_color(depth1), 2);//画出
// 第二个图：世界坐标旋转平移
Mat pt2_trans = R * (Mat_<double>(3, 1) << points[i].x, points[i].y, points[i].z) + t;
float depth2 = pt2_trans.at<double>(2, 0);
cv::circle(img2_plot, keypoints_2[matches[i].trainIdx].pt, 2, get_color(depth2), 2);
}
cv::imshow("img 1", img1_plot);
cv::imshow("img 2", img2_plot);
cv::waitKey();
return 0;
}
void find_feature_matches(const Mat &img_1, const Mat &img_2,
std::vector<KeyPoint> &keypoints_1,
std::vector<KeyPoint> &keypoints_2,
std::vector<DMatch> &matches) {
//-- 初始化
Mat descriptors_1, descriptors_2;
// used in OpenCV3
Ptr<FeatureDetector> detector = ORB::create();
Ptr<DescriptorExtractor> descriptor = ORB::create();
// use this if you are in OpenCV2
// Ptr<FeatureDetector> detector = FeatureDetector::create ( "ORB" );
// Ptr<DescriptorExtractor> descriptor = DescriptorExtractor::create ( "ORB" );
Ptr<DescriptorMatcher> matcher = DescriptorMatcher::create("BruteForce-Hamming");
//-- 第一步:检测 Oriented FAST 角点位置
detector->detect(img_1, keypoints_1);
detector->detect(img_2, keypoints_2);
//-- 第二步:根据角点位置计算 BRIEF 描述子
descriptor->compute(img_1, keypoints_1, descriptors_1);
descriptor->compute(img_2, keypoints_2, descriptors_2);
//-- 第三步:对两幅图像中的BRIEF描述子进行匹配，使用 Hamming 距离
vector<DMatch> match;
// BFMatcher matcher ( NORM_HAMMING );
matcher->match(descriptors_1, descriptors_2, match);
//-- 第四步:匹配点对筛选
double min_dist = 10000, max_dist = 0;
//找出所有匹配之间的最小距离和最大距离, 即是最相似的和最不相似的两组点之间的距离
for (int i = 0; i < descriptors_1.rows; i++) {
double dist = match[i].distance;
if (dist < min_dist) min_dist = dist;
if (dist > max_dist) max_dist = dist;
}
printf("-- Max dist : %f \n", max_dist);
printf("-- Min dist : %f \n", min_dist);
//当描述子之间的距离大于两倍的最小距离时,即认为匹配有误.但有时候最小距离会非常小,设置一个经验值30作为下限.
for (int i = 0; i < descriptors_1.rows; i++) {
if (match[i].distance <= max(2 * min_dist, 30.0)) {
matches.push_back(match[i]);
}
}
}
void pose_estimation_2d2d(
const std::vector<KeyPoint> &keypoints_1,
const std::vector<KeyPoint> &keypoints_2,
const std::vector<DMatch> &matches,
Mat &R, Mat &t) {
// 相机内参,TUM Freiburg2
Mat K = (Mat_<double>(3, 3) << 520.9, 0, 325.1, 0, 521.0, 249.7, 0, 0, 1);
//-- 把匹配点转换为vector<Point2f>的形式
vector<Point2f> points1;
vector<Point2f> points2;
for (int i = 0; i < (int) matches.size(); i++) {
points1.push_back(keypoints_1[matches[i].queryIdx].pt);
points2.push_back(keypoints_2[matches[i].trainIdx].pt);
}
//-- 计算本质矩阵
Point2d principal_point(325.1, 249.7); //相机主点, TUM dataset标定值
int focal_length = 521; //相机焦距, TUM dataset标定值
Mat essential_matrix;
essential_matrix = findEssentialMat(points1, points2, focal_length, principal_point);
//-- 从本质矩阵中恢复旋转和平移信息.
recoverPose(essential_matrix, points1, points2, R, t, focal_length, principal_point);
}
void triangulation(
const vector<KeyPoint> &keypoint_1,
const vector<KeyPoint> &keypoint_2,
const std::vector<DMatch> &matches,
const Mat &R, const Mat &t,
vector<Point3d> &points) {
Mat T1 = (Mat_<float>(3, 4) <<
1, 0, 0, 0,
0, 1, 0, 0,
0, 0, 1, 0);
//外参
Mat T2 = (Mat_<float>(3, 4) <<
R.at<double>(0, 0), R.at<double>(0, 1), R.at<double>(0, 2), t.at<double>(0, 0),
R.at<double>(1, 0), R.at<double>(1, 1), R.at<double>(1, 2), t.at<double>(1, 0),
R.at<double>(2, 0), R.at<double>(2, 1), R.at<double>(2, 2), t.at<double>(2, 0)
);
Mat K = (Mat_<double>(3, 3) << 520.9, 0, 325.1, 0, 521.0, 249.7, 0, 0, 1);//内参
vector<Point2f> pts_1, pts_2;
for (DMatch m:matches) {
// 将像素坐标转换至相机坐标
pts_1.push_back(pixel2cam(keypoint_1[m.queryIdx].pt, K));//相机坐标1
pts_2.push_back(pixel2cam(keypoint_2[m.trainIdx].pt, K));//相机坐标2
}
Mat pts_4d;
// pts_4d表示世界坐标（齐次的）
cv::triangulatePoints(T1, T2, pts_1, pts_2, pts_4d);
// 转换成非齐次坐标
for (int i = 0; i < pts_4d.cols; i++) {
Mat x = pts_4d.col(i);
x /= x.at<float>(3, 0); // 归一化
Point3d p(
x.at<float>(0, 0),
x.at<float>(1, 0),
x.at<float>(2, 0)
);
points.push_back(p);
}
}
Point2f pixel2cam(const Point2d &p, const Mat &K) {
return Point2f
(
(p.x - K.at<double>(0, 2)) / K.at<double>(0, 0),
(p.y - K.at<double>(1, 2)) / K.at<double>(1, 1)
);
}
5.4 结果
六、PnP法（3D->2D）：已知特征点世界坐标及其投影位置估计相机位姿
6.1 理论部分
- PnP：已知特征点世界坐标及其投影位置估计相机位姿
- 使用场景：
1）在单目视觉里程计中，先用对极几何初始化，转化成3D->2D，再使用PnP；
2）在双目/RGB-D中，直接使用PnP估计相机运动 - PnP问题求解方法？
1）线性变换：P3P(3对点)，直接线性变换DLT(6对点)，EPnP，UPnP等
2）非线性变换：构建最小二乘法并迭代求解，即光束法平差(Bundle Adjustment,BA)法
注意：线性变换，是先求解相机位姿，再估计空间点位置；而非线性优化，将它们都看作优化变量，放在一起优化。 - 直接线性变换DLT：已知3D空间点及其投影位置估计相机位姿。
- P3P：已知3对3D空间点及其投影位置，可以求出3个空间点在相机坐标下的3D坐标。从而将3D->2D问题转化成3D->3D问题，利用ICP求解相机位姿。
- BA：这里利用第六讲非线性优化的知识，采用G-N、L-M方法求解相机位姿。针对不同问题，雅可比和海塞阵的方式不同，本节实验针对的是的是课本185-187页。
6.2 CMakeLists.txt
参见2.1.1
6.3 代码展示
按照课本给出的问题（其中曲线函数和雅可比矩阵参见P187），这里使用3种方法求解PnP：
1）OpenCV自带的PnP求解法
2）自己写的高斯牛顿的PnP
3）利用g2o实现BA的PnP
注意：代码一定要结合第六讲理解
#include <iostream>
#include <opencv2/core/core.hpp>
#include <opencv2/features2d/features2d.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/calib3d/calib3d.hpp>
#include <Eigen/Core>
#include <g2o/core/base_vertex.h>
#include <g2o/core/base_unary_edge.h>
#include <g2o/core/sparse_optimizer.h>
#include <g2o/core/block_solver.h>
#include <g2o/core/solver.h>
#include <g2o/core/optimization_algorithm_gauss_newton.h>
#include <g2o/solvers/dense/linear_solver_dense.h>
#include <sophus/se3.hpp>
#include <chrono>
#include "opencv2/imgcodecs/legacy/constants_c.h"
using namespace std;
using namespace cv;
void find_feature_matches(
const Mat &img_1, const Mat &img_2,
std::vector<KeyPoint> &keypoints_1,
std::vector<KeyPoint> &keypoints_2,
std::vector<DMatch> &matches);
// 像素坐标转相机归一化坐标
Point2d pixel2cam(const Point2d &p, const Mat &K);
// 通过g2o实现BA
typedef vector<Eigen::Vector2d, Eigen::aligned_allocator<Eigen::Vector2d>> VecVector2d;
typedef vector<Eigen::Vector3d, Eigen::aligned_allocator<Eigen::Vector3d>> VecVector3d;
void bundleAdjustmentG2O(
const VecVector3d &points_3d,
const VecVector2d &points_2d,
const Mat &K,
Sophus::SE3d &pose
);
// 通过高斯牛顿实现BA
void bundleAdjustmentGaussNewton(
const VecVector3d &points_3d,
const VecVector2d &points_2d,
const Mat &K,
Sophus::SE3d &pose
);
int main(int argc, char **argv) {
// 1. 读取图像
Mat img_1 = imread("/home/robot/桌面/slambook2-master/ch7/1.png", cv::IMREAD_COLOR);
Mat img_2 = imread("/home/robot/桌面/slambook2-master/ch7/2.png", cv::IMREAD_COLOR);
assert(img_1.data && img_2.data && "Can not load images!");
// 2. 特征点提取与匹配
vector<KeyPoint> keypoints_1, keypoints_2;
vector<DMatch> matches;
find_feature_matches(img_1, img_2, keypoints_1, keypoints_2, matches);
cout << "一共找到了" << matches.size() << "组匹配点" << endl;
// 3. 建立3D和2D点
Mat d1 = imread("/home/robot/桌面/slambook2-master/ch7/1_depth.png", CV_LOAD_IMAGE_UNCHANGED); // 深度图为16位无符号数，单通道图像
Mat K = (Mat_<double>(3, 3) << 520.9, 0, 325.1, 0, 521.0, 249.7, 0, 0, 1);
vector<Point3f> pts_3d;
vector<Point2f> pts_2d;
for (DMatch m:matches) {
ushort d = d1.ptr<unsigned short>(int(keypoints_1[m.queryIdx].pt.y))[int(keypoints_1[m.queryIdx].pt.x)];
if (d == 0) // bad depth
continue;
float dd = d / 5000.0;
Point2d p1 = pixel2cam(keypoints_1[m.queryIdx].pt, K);//归一化坐标
pts_3d.push_back(Point3f(p1.x * dd, p1.y * dd, dd));//世界坐标
pts_2d.push_back(keypoints_2[m.trainIdx].pt);//像素坐标
}
cout << "3d-2d pairs: " << pts_3d.size() << endl;
// 4.1 使用OpenCV的PnP求解
chrono::steady_clock::time_point t1 = chrono::steady_clock::now();
Mat r, t;
// solvePnP(pts_3d, pts_2d, K, Mat(), r, t, false, SOLVEPNP_EPNP); //使用 EPnP
solvePnP(pts_3d, pts_2d, K, Mat(), r, t, false); // 调用OpenCV 的 PnP 求解，可选择EPNP，DLS等方法
Mat R;
cv::Rodrigues(r, R); // r为旋转向量形式，用Rodrigues公式转换为矩阵
chrono::steady_clock::time_point t2 = chrono::steady_clock::now();
chrono::duration<double> time_used = chrono::duration_cast<chrono::duration<double>>(t2 - t1);
cout << "solve pnp in opencv cost time: " << time_used.count() << " seconds." << endl;
cout << "R=" << endl << R << endl;
cout << "t=" << endl << t << endl;
VecVector3d pts_3d_eigen;
VecVector2d pts_2d_eigen;
for (size_t i = 0; i < pts_3d.size(); ++i) {
pts_3d_eigen.push_back(Eigen::Vector3d(pts_3d[i].x, pts_3d[i].y, pts_3d[i].z));
pts_2d_eigen.push_back(Eigen::Vector2d(pts_2d[i].x, pts_2d[i].y));
}
// 4.2 使用高斯牛顿的PnP求解
cout << "calling bundle adjustment by gauss newton" << endl;
Sophus::SE3d pose_gn;// 位姿
t1 = chrono::steady_clock::now();
bundleAdjustmentGaussNewton(pts_3d_eigen, pts_2d_eigen, K, pose_gn);// 使用高斯牛顿的BA
t2 = chrono::steady_clock::now();
time_used = chrono::duration_cast<chrono::duration<double>>(t2 - t1);
cout << "solve pnp by gauss newton cost time: " << time_used.count() << " seconds." << endl;
// 4.3 使用g2o的PnP求解
cout << "calling bundle adjustment by g2o" << endl;
Sophus::SE3d pose_g2o;
t1 = chrono::steady_clock::now();
bundleAdjustmentG2O(pts_3d_eigen, pts_2d_eigen, K, pose_g2o); //使用g2o的BA
t2 = chrono::steady_clock::now();
time_used = chrono::duration_cast<chrono::duration<double>>(t2 - t1);
cout << "solve pnp by g2o cost time: " << time_used.count() << " seconds." << endl;
return 0;
}
void find_feature_matches(const Mat &img_1, const Mat &img_2,
std::vector<KeyPoint> &keypoints_1,
std::vector<KeyPoint> &keypoints_2,
std::vector<DMatch> &matches) {
//-- 初始化
Mat descriptors_1, descriptors_2;
// used in OpenCV3
Ptr<FeatureDetector> detector = ORB::create();
Ptr<DescriptorExtractor> descriptor = ORB::create();
// use this if you are in OpenCV2
// Ptr<FeatureDetector> detector = FeatureDetector::create ( "ORB" );
// Ptr<DescriptorExtractor> descriptor = DescriptorExtractor::create ( "ORB" );
Ptr<DescriptorMatcher> matcher = DescriptorMatcher::create("BruteForce-Hamming");
//-- 第一步:检测 Oriented FAST 角点位置
detector->detect(img_1, keypoints_1);
detector->detect(img_2, keypoints_2);
//-- 第二步:根据角点位置计算 BRIEF 描述子
descriptor->compute(img_1, keypoints_1, descriptors_1);
descriptor->compute(img_2, keypoints_2, descriptors_2);
//-- 第三步:对两幅图像中的BRIEF描述子进行匹配，使用 Hamming 距离
vector<DMatch> match;
// BFMatcher matcher ( NORM_HAMMING );
matcher->match(descriptors_1, descriptors_2, match);
//-- 第四步:匹配点对筛选
double min_dist = 10000, max_dist = 0;
//找出所有匹配之间的最小距离和最大距离, 即是最相似的和最不相似的两组点之间的距离
for (int i = 0; i < descriptors_1.rows; i++) {
double dist = match[i].distance;
if (dist < min_dist) min_dist = dist;
if (dist > max_dist) max_dist = dist;
}
printf("-- Max dist : %f \n", max_dist);
printf("-- Min dist : %f \n", min_dist);
//当描述子之间的距离大于两倍的最小距离时,即认为匹配有误.但有时候最小距离会非常小,设置一个经验值30作为下限.
for (int i = 0; i < descriptors_1.rows; i++) {
if (match[i].distance <= max(2 * min_dist, 30.0)) {
matches.push_back(match[i]);
}
}
}
Point2d pixel2cam(const Point2d &p, const Mat &K) {
return Point2d
(
(p.x - K.at<double>(0, 2)) / K.at<double>(0, 0),
(p.y - K.at<double>(1, 2)) / K.at<double>(1, 1)
);
}
void bundleAdjustmentGaussNewton(
const VecVector3d &points_3d,
const VecVector2d &points_2d,
const Mat &K,
Sophus::SE3d &pose) {
typedef Eigen::Matrix<double, 6, 1> Vector6d;
const int iterations = 10;
double cost = 0, lastCost = 0;
// 内参
double fx = K.at<double>(0, 0);
double fy = K.at<double>(1, 1);
double cx = K.at<double>(0, 2);
double cy = K.at<double>(1, 2);
for (int iter = 0; iter < iterations; iter++) {// 高斯牛顿迭代次数
// 1）构建Hx=b
Eigen::Matrix<double, 6, 6> H = Eigen::Matrix<double, 6, 6>::Zero();
Vector6d b = Vector6d::Zero();
cost = 0;
// 2）计算成本compute cost
for (int i = 0; i < points_3d.size(); i++) {
Eigen::Vector3d pc = pose * points_3d[i];
double inv_z = 1.0 / pc[2];
double inv_z2 = inv_z * inv_z;
Eigen::Vector2d proj(fx * pc[0] / pc[2] + cx, fy * pc[1] / pc[2] + cy);
Eigen::Vector2d e = points_2d[i] - proj;//误差
cost += e.squaredNorm();
// 计算2✖6的雅可比，公式P187
Eigen::Matrix<double, 2, 6> J;
J << -fx * inv_z,
0,
fx * pc[0] * inv_z2,
fx * pc[0] * pc[1] * inv_z2,
-fx - fx * pc[0] * pc[0] * inv_z2,
fx * pc[1] * inv_z,
0,
-fy * inv_z,
fy * pc[1] * inv_z2,
fy + fy * pc[1] * pc[1] * inv_z2,
-fy * pc[0] * pc[1] * inv_z2,
-fy * pc[0] * inv_z;
// 利用高斯牛顿公式 P129
H += J.transpose() * J;
b += -J.transpose() * e;
}
// 3）Hx=b求出dx
Vector6d dx;
dx = H.ldlt().solve(b);
if (isnan(dx[0])) {
cout << "result is nan!" << endl;
break;
}
if (iter > 0 && cost >= lastCost) {
// cost increase, update is not good
cout << "cost: " << cost << ", last cost: " << lastCost << endl;
break;
}
// 4）利用误差估计更新位姿
pose = Sophus::SE3d::exp(dx) * pose;
lastCost = cost;
cout << "iteration " << iter << " cost=" << std::setprecision(12) << cost << endl;
if (dx.norm() < 1e-6) {
// converge
break;
}
}
cout << "pose by g-n: \n" << pose.matrix() << endl;
}
// 这里实现曲线模型的顶点（位姿）的更新
class VertexPose : public g2o::BaseVertex<6, Sophus::SE3d> {
public:
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
virtual void setToOriginImpl() override {
_estimate = Sophus::SE3d();
}
// 使用SE3的左乘法更新
virtual void oplusImpl(const double *update) override {
Eigen::Matrix<double, 6, 1> update_eigen;
update_eigen << update[0], update[1], update[2], update[3], update[4], update[5];
_estimate = Sophus::SE3d::exp(update_eigen) * _estimate;
}
virtual bool read(istream &in) override {}
virtual bool write(ostream &out) const override {}
};
// 这里实现曲线边的误差计算
class EdgeProjection : public g2o::BaseUnaryEdge<2, Eigen::Vector2d, VertexPose> {
public:
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
EdgeProjection(const Eigen::Vector3d &pos, const Eigen::Matrix3d &K) : _pos3d(pos), _K(K) {}
// 计算曲线模型误差
virtual void computeError() override {
const VertexPose *v = static_cast<VertexPose *> (_vertices[0]);
Sophus::SE3d T = v->estimate();
Eigen::Vector3d pos_pixel = _K * (T * _pos3d);
pos_pixel /= pos_pixel[2];
// w = y - exp(K*T*P/S) P185
_error = _measurement - pos_pixel.head<2>();
}
// 计算雅可比矩阵P187
virtual void linearizeOplus() override {
const VertexPose *v = static_cast<VertexPose *> (_vertices[0]);
Sophus::SE3d T = v->estimate();
Eigen::Vector3d pos_cam = T * _pos3d;
double fx = _K(0, 0);
double fy = _K(1, 1);
double cx = _K(0, 2);
double cy = _K(1, 2);
double X = pos_cam[0];
double Y = pos_cam[1];
double Z = pos_cam[2];
double Z2 = Z * Z;
_jacobianOplusXi
<< -fx / Z, 0, fx * X / Z2, fx * X * Y / Z2, -fx - fx * X * X / Z2, fx * Y / Z,
0, -fy / Z, fy * Y / (Z * Z), fy + fy * Y * Y / Z2, -fy * X * Y / Z2, -fy * X / Z;
}
virtual bool read(istream &in) override {}
virtual bool write(ostream &out) const override {}
private:
Eigen::Vector3d _pos3d;
Eigen::Matrix3d _K;
};
void bundleAdjustmentG2O(
const VecVector3d &points_3d,
const VecVector2d &points_2d,
const Mat &K,
Sophus::SE3d &pose) {
// 构建图优化，先设定g2o
typedef g2o::BlockSolver<g2o::BlockSolverTraits<6, 3>> BlockSolverType; // pose is 6, landmark is 3
typedef g2o::LinearSolverDense<BlockSolverType::PoseMatrixType> LinearSolverType; // 线性求解器类型
// 梯度下降方法，可以从GN, LM, DogLeg 中选
auto solver = new g2o::OptimizationAlgorithmGaussNewton(
g2o::make_unique<BlockSolverType>(g2o::make_unique<LinearSolverType>()));
g2o::SparseOptimizer optimizer; // 图模型
optimizer.setAlgorithm(solver); // 设置求解器
optimizer.setVerbose(true); // 打开调试输出
// 往图中增加顶点
VertexPose *vertex_pose = new VertexPose(); // camera vertex_pose
vertex_pose->setId(0);
vertex_pose->setEstimate(Sophus::SE3d());
optimizer.addVertex(vertex_pose);
// 内参K
Eigen::Matrix3d K_eigen;
K_eigen <<
K.at<double>(0, 0), K.at<double>(0, 1), K.at<double>(0, 2),
K.at<double>(1, 0), K.at<double>(1, 1), K.at<double>(1, 2),
K.at<double>(2, 0), K.at<double>(2, 1), K.at<double>(2, 2);
// 往图中增加边
int index = 1;
for (size_t i = 0; i < points_2d.size(); ++i) {
auto p2d = points_2d[i];
auto p3d = points_3d[i];
EdgeProjection *edge = new EdgeProjection(p3d, K_eigen);
edge->setId(index);
edge->setVertex(0, vertex_pose);
edge->setMeasurement(p2d);
edge->setInformation(Eigen::Matrix2d::Identity());
optimizer.addEdge(edge);
index++;
}
// 执行优化
chrono::steady_clock::time_point t1 = chrono::steady_clock::now();
optimizer.setVerbose(true);
optimizer.initializeOptimization();
optimizer.optimize(10);
chrono::steady_clock::time_point t2 = chrono::steady_clock::now();
chrono::duration<double> time_used = chrono::duration_cast<chrono::duration<double>>(t2 - t1);
cout << "optimization costs time: " << time_used.count() << " seconds." << endl;
// 输出优化值
cout << "pose estimated by g2o =\n" << vertex_pose->estimate().matrix() << endl;
pose = vertex_pose->estimate();
}
6.4 结果
-- Max dist : 94.000000
-- Min dist : 4.000000
一共找到了79组匹配点
3d-2d pairs: 75
solve pnp in opencv cost time: 0.000473166 seconds.
R=
[0.9979059095501289, -0.05091940089111062, 0.03988747043647115;
0.04981866254254162, 0.9983623157438141, 0.02812094175381178;
-0.04125404886071617, -0.02607491352889358, 0.9988083912027663]
t=
[-0.1267821389556796;
-0.008439496817594587;
0.06034935748886031]
calling bundle adjustment by gauss newton
iteration 0 cost=40517.7576706
iteration 1 cost=410.547029116
iteration 2 cost=299.76468142
iteration 3 cost=299.763574327
pose by g-n:
0.997905909549 -0.0509194008562 0.0398874705187 -0.126782139096
0.049818662505 0.998362315745 0.0281209417649 -0.00843949683874
-0.0412540489424 -0.0260749135374 0.998808391199 0.0603493575229
0 0 0 1
solve pnp by gauss newton cost time: 6.6108e-05 seconds.
calling bundle adjustment by g2o
iteration= 0 chi2= 410.547029 time= 1.6601e-05 cumTime= 1.6601e-05 edges= 75 schur= 0
iteration= 1 chi2= 299.764681 time= 8.169e-06 cumTime= 2.477e-05 edges= 75 schur= 0
iteration= 2 chi2= 299.763574 time= 7.401e-06 cumTime= 3.2171e-05 edges= 75 schur= 0
iteration= 3 chi2= 299.763574 time= 7.201e-06 cumTime= 3.9372e-05 edges= 75 schur= 0
iteration= 4 chi2= 299.763574 time= 7.187e-06 cumTime= 4.6559e-05 edges= 75 schur= 0
iteration= 5 chi2= 299.763574 time= 7.188e-06 cumTime= 5.3747e-05 edges= 75 schur= 0
iteration= 6 chi2= 299.763574 time= 7.181e-06 cumTime= 6.0928e-05 edges= 75 schur= 0
iteration= 7 chi2= 299.763574 time= 7.151e-06 cumTime= 6.8079e-05 edges= 75 schur= 0
iteration= 8 chi2= 299.763574 time= 7.162e-06 cumTime= 7.5241e-05 edges= 75 schur= 0
iteration= 9 chi2= 299.763574 time= 7.129e-06 cumTime= 8.237e-05 edges= 75 schur= 0
optimization costs time: 0.000256374 seconds.
pose estimated by g2o =
0.99790590955 -0.0509194008911 0.0398874704367 -0.126782138956
0.0498186625425 0.998362315744 0.0281209417542 -0.00843949681823
-0.0412540488609 -0.0260749135293 0.998808391203 0.0603493574888
0 0 0 1
solve pnp by g2o cost time: 0.000329383 seconds.
七、ICP法（3D->3D）：根据一组匹配好了的3D点（世界坐标）估计其位姿
7.1 理论部分
- 迭代最近点法ICP：已经有一组匹配好了的3D点估计其位姿。迭代最近点指，认为距离最近的两个点为同一个。
- ICP的求解方式？
1）线性变换：SVD法 P197
2）非线性优化：BA P198
7.2 CMakeLists.txt
参见2.1.1
7.3 代码展示
按照课本给出的问题P196-199，这里使用2种方法求解ICP：
1）SVD法
2）利用g2o的BA法
注意：代码一定要结合第六讲理解
#include <iostream>
#include <opencv2/core/core.hpp>
#include <opencv2/features2d/features2d.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/calib3d/calib3d.hpp>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <Eigen/Geometry>
#include <Eigen/SVD>
#include <g2o/core/base_vertex.h>
#include <g2o/core/base_unary_edge.h>
#include <g2o/core/block_solver.h>
#include <g2o/core/optimization_algorithm_gauss_newton.h>
#include <g2o/core/optimization_algorithm_levenberg.h>
#include <g2o/solvers/dense/linear_solver_dense.h>
#include <chrono>
#include <sophus/se3.hpp>
#include "opencv2/imgcodecs/legacy/constants_c.h"
using namespace std;
using namespace cv;
void find_feature_matches(
const Mat &img_1, const Mat &img_2,
std::vector<KeyPoint> &keypoints_1,
std::vector<KeyPoint> &keypoints_2,
std::vector<DMatch> &matches);
// 像素坐标转相机归一化坐标
Point2d pixel2cam(const Point2d &p, const Mat &K);
void pose_estimation_3d3d(
const vector<Point3f> &pts1,
const vector<Point3f> &pts2,
Mat &R, Mat &t
);
void bundleAdjustment(
const vector<Point3f> &points_3d,
const vector<Point3f> &points_2d,
Mat &R, Mat &t
);
/// 这里实现曲线模型的顶点（位姿）的更新
class VertexPose : public g2o::BaseVertex<6, Sophus::SE3d> {
public:
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
virtual void setToOriginImpl() override {
_estimate = Sophus::SE3d();
}
/// 使用SE3的左乘法更新
virtual void oplusImpl(const double *update) override {
Eigen::Matrix<double, 6, 1> update_eigen;
update_eigen << update[0], update[1], update[2], update[3], update[4], update[5];
_estimate = Sophus::SE3d::exp(update_eigen) * _estimate;
}
virtual bool read(istream &in) override {}
virtual bool write(ostream &out) const override {}
};
/// 这里实现曲线边的误差计算
class EdgeProjectXYZRGBDPoseOnly : public g2o::BaseUnaryEdge<3, Eigen::Vector3d, VertexPose> {
public:
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
EdgeProjectXYZRGBDPoseOnly(const Eigen::Vector3d &point) : _point(point) {}
// 计算曲线模型误差
virtual void computeError() override {
const VertexPose *pose = static_cast<const VertexPose *> ( _vertices[0] );
// w = y - (R * p + t) P196 这里pose->estimate将R和t合并了
_error = _measurement - pose->estimate() * _point;
}
// 计算雅可比矩阵P199
virtual void linearizeOplus() override {
VertexPose *pose = static_cast<VertexPose *>(_vertices[0]);
Sophus::SE3d T = pose->estimate();
Eigen::Vector3d xyz_trans = T * _point;
_jacobianOplusXi.block<3, 3>(0, 0) = -Eigen::Matrix3d::Identity();
_jacobianOplusXi.block<3, 3>(0, 3) = Sophus::SO3d::hat(xyz_trans);
}
bool read(istream &in) {}
bool write(ostream &out) const {}
protected:
Eigen::Vector3d _point;
};
int main(int argc, char **argv) {
// 1. 读取图像
Mat img_1 = imread("/home/robot/桌面/slambook2-master/ch7/1.png", cv::IMREAD_COLOR);
Mat img_2 = imread("/home/robot/桌面/slambook2-master/ch7/2.png", cv::IMREAD_COLOR);
// 2. 特征点的提取与匹配
vector<KeyPoint> keypoints_1, keypoints_2;
vector<DMatch> matches;
find_feature_matches(img_1, img_2, keypoints_1, keypoints_2, matches);
cout << "一共找到了" << matches.size() << "组匹配点" << endl;
// 3. 建立两个3D点
Mat depth1 = imread("/home/robot/桌面/slambook2-master/ch7/1_depth.png", CV_LOAD_IMAGE_UNCHANGED); // 深度图为16位无符号数，单通道图像
Mat depth2 = imread("/home/robot/桌面/slambook2-master/ch7/2_depth.png", CV_LOAD_IMAGE_UNCHANGED); // 深度图为16位无符号数，单通道图像
Mat K = (Mat_<double>(3, 3) << 520.9, 0, 325.1, 0, 521.0, 249.7, 0, 0, 1);
vector<Point3f> pts1, pts2;
for (DMatch m:matches) {
ushort d1 = depth1.ptr<unsigned short>(int(keypoints_1[m.queryIdx].pt.y))[int(keypoints_1[m.queryIdx].pt.x)];
ushort d2 = depth2.ptr<unsigned short>(int(keypoints_2[m.trainIdx].pt.y))[int(keypoints_2[m.trainIdx].pt.x)];
if (d1 == 0 || d2 == 0) // bad depth
continue;
Point2d p1 = pixel2cam(keypoints_1[m.queryIdx].pt, K);// 像素坐标->归一化坐标
Point2d p2 = pixel2cam(keypoints_2[m.trainIdx].pt, K);// 像素坐标->归一化坐标
float dd1 = float(d1) / 5000.0;
float dd2 = float(d2) / 5000.0;
pts1.push_back(Point3f(p1.x * dd1, p1.y * dd1, dd1));// 归一化坐标 * 深度-> 世界坐标
pts2.push_back(Point3f(p2.x * dd2, p2.y * dd2, dd2));// 归一化坐标 * 深度-> 世界坐标
}
// 4. SVD方法：算法参见课本P197
cout << "3d-3d pairs: " << pts1.size() << endl;
Mat R, t;
pose_estimation_3d3d(pts1, pts2, R, t);
cout << "ICP via SVD results: " << endl;
cout << "R = " << R << endl;
cout << "t = " << t << endl;
cout << "R_inv = " << R.t() << endl;
cout << "t_inv = " << -R.t() * t << endl;
// 5. BA方法：：算法参见课本P198
cout << "calling bundle adjustment" << endl;
bundleAdjustment(pts1, pts2, R, t);
// 6. 验证 p1 = R * p2 + t
for (int i = 0; i < 5; i++) {
cout << "p1 = " << pts1[i] << endl;
cout << "p2 = " << pts2[i] << endl;
cout << "P1 = (R*p2+t) = " <<
R * (Mat_<double>(3, 1) << pts2[i].x, pts2[i].y, pts2[i].z) + t
<< endl;
cout << endl;
}
return 0;
}
void find_feature_matches(const Mat &img_1, const Mat &img_2,
std::vector<KeyPoint> &keypoints_1,
std::vector<KeyPoint> &keypoints_2,
std::vector<DMatch> &matches) {
//-- 初始化
Mat descriptors_1, descriptors_2;
// used in OpenCV3
Ptr<FeatureDetector> detector = ORB::create();
Ptr<DescriptorExtractor> descriptor = ORB::create();
// use this if you are in OpenCV2
// Ptr<FeatureDetector> detector = FeatureDetector::create ( "ORB" );
// Ptr<DescriptorExtractor> descriptor = DescriptorExtractor::create ( "ORB" );
Ptr<DescriptorMatcher> matcher = DescriptorMatcher::create("BruteForce-Hamming");
//-- 第一步:检测 Oriented FAST 角点位置
detector->detect(img_1, keypoints_1);
detector->detect(img_2, keypoints_2);
//-- 第二步:根据角点位置计算 BRIEF 描述子
descriptor->compute(img_1, keypoints_1, descriptors_1);
descriptor->compute(img_2, keypoints_2, descriptors_2);
//-- 第三步:对两幅图像中的BRIEF描述子进行匹配，使用 Hamming 距离
vector<DMatch> match;
// BFMatcher matcher ( NORM_HAMMING );
matcher->match(descriptors_1, descriptors_2, match);
//-- 第四步:匹配点对筛选
double min_dist = 10000, max_dist = 0;
//找出所有匹配之间的最小距离和最大距离, 即是最相似的和最不相似的两组点之间的距离
for (int i = 0; i < descriptors_1.rows; i++) {
double dist = match[i].distance;
if (dist < min_dist) min_dist = dist;
if (dist > max_dist) max_dist = dist;
}
printf("-- Max dist : %f \n", max_dist);
printf("-- Min dist : %f \n", min_dist);
//当描述子之间的距离大于两倍的最小距离时,即认为匹配有误.但有时候最小距离会非常小,设置一个经验值30作为下限.
for (int i = 0; i < descriptors_1.rows; i++) {
if (match[i].distance <= max(2 * min_dist, 30.0)) {
matches.push_back(match[i]);
}
}
}
Point2d pixel2cam(const Point2d &p, const Mat &K) {
return Point2d(
(p.x - K.at<double>(0, 2)) / K.at<double>(0, 0),
(p.y - K.at<double>(1, 2)) / K.at<double>(1, 1)
);
}
void pose_estimation_3d3d(const vector<Point3f> &pts1,
const vector<Point3f> &pts2,
Mat &R, Mat &t) {
// 1. 定义两组质心
Point3f p1, p2;
int N = pts1.size();
for (int i = 0; i < N; i++) {
p1 += pts1[i];
p2 += pts2[i];
}
p1 = Point3f(Vec3f(p1) / N);
p2 = Point3f(Vec3f(p2) / N);
// 2. 定义两组去质心坐标
vector<Point3f> q1(N), q2(N);
for (int i = 0; i < N; i++) {
q1[i] = pts1[i] - p1;
q2[i] = pts2[i] - p2;
}
// 3. 计算W=q1*q2^T
Eigen::Matrix3d W = Eigen::Matrix3d::Zero();
for (int i = 0; i < N; i++) {
W += Eigen::Vector3d(q1[i].x, q1[i].y, q1[i].z) * Eigen::Vector3d(q2[i].x, q2[i].y, q2[i].z).transpose();
}
cout << "W=" << W << endl;
// 4. 使用SVD分解W=U*对角阵*V^T
Eigen::JacobiSVD<Eigen::Matrix3d> svd(W, Eigen::ComputeFullU | Eigen::ComputeFullV);
Eigen::Matrix3d U = svd.matrixU();
Eigen::Matrix3d V = svd.matrixV();
cout << "U=" << U << endl;
cout << "V=" << V << endl;
// 5. 求出R=U*V^T
Eigen::Matrix3d R_ = U * (V.transpose());
if (R_.determinant() < 0) {
R_ = -R_;
}
// 6. 求出t
Eigen::Vector3d t_ = Eigen::Vector3d(p1.x, p1.y, p1.z) - R_ * Eigen::Vector3d(p2.x, p2.y, p2.z);
// 7. 将R和t合并
R = (Mat_<double>(3, 3) <<
R_(0, 0), R_(0, 1), R_(0, 2),
R_(1, 0), R_(1, 1), R_(1, 2),
R_(2, 0), R_(2, 1), R_(2, 2)
);
t = (Mat_<double>(3, 1) << t_(0, 0), t_(1, 0), t_(2, 0));
}
void bundleAdjustment(const vector<Point3f> &pts1, const vector<Point3f> &pts2, Mat &R, Mat &t) {
// 构建图优化，先设定g2o
typedef g2o::BlockSolverX BlockSolverType;
typedef g2o::LinearSolverDense<BlockSolverType::PoseMatrixType> LinearSolverType; // 线性求解器类型
// 梯度下降方法，可以从GN, LM, DogLeg 中选
auto solver = new g2o::OptimizationAlgorithmLevenberg(
g2o::make_unique<BlockSolverType>(g2o::make_unique<LinearSolverType>()));
g2o::SparseOptimizer optimizer; // 图模型
optimizer.setAlgorithm(solver); // 设置求解器
optimizer.setVerbose(true); // 打开调试输出
// vertex
VertexPose *pose = new VertexPose(); // camera pose
pose->setId(0);
pose->setEstimate(Sophus::SE3d());
optimizer.addVertex(pose);
// edges
for (size_t i = 0; i < pts1.size(); i++) {
EdgeProjectXYZRGBDPoseOnly *edge = new EdgeProjectXYZRGBDPoseOnly(
Eigen::Vector3d(pts2[i].x, pts2[i].y, pts2[i].z));
edge->setVertex(0, pose);
edge->setMeasurement(Eigen::Vector3d(
pts1[i].x, pts1[i].y, pts1[i].z));
edge->setInformation(Eigen::Matrix3d::Identity());
optimizer.addEdge(edge);
}
chrono::steady_clock::time_point t1 = chrono::steady_clock::now();
optimizer.initializeOptimization();
optimizer.optimize(10);
chrono::steady_clock::time_point t2 = chrono::steady_clock::now();
chrono::duration<double> time_used = chrono::duration_cast<chrono::duration<double>>(t2 - t1);
cout << "optimization costs time: " << time_used.count() << " seconds." << endl;
cout << endl << "after optimization:" << endl;
cout << "T=\n" << pose->estimate().matrix() << endl;
// convert to cv::Mat
Eigen::Matrix3d R_ = pose->estimate().rotationMatrix();
Eigen::Vector3d t_ = pose->estimate().translation();
R = (Mat_<double>(3, 3) <<
R_(0, 0), R_(0, 1), R_(0, 2),
R_(1, 0), R_(1, 1), R_(1, 2),
R_(2, 0), R_(2, 1), R_(2, 2)
);
t = (Mat_<double>(3, 1) << t_(0, 0), t_(1, 0), t_(2, 0));
}
7.4 结果
-- Max dist : 94.000000
-- Min dist : 4.000000
一共找到了79组匹配点
3d-3d pairs: 72
W= 10.871 -1.01948 2.54771
-2.16033 3.85307 -5.77742
3.94738 -5.79979 9.62203
U= 0.558087 -0.829399 -0.0252034
-0.428009 -0.313755 0.847565
0.710878 0.462228 0.530093
V= 0.617887 -0.784771 -0.0484806
-0.399894 -0.366747 0.839989
0.676979 0.499631 0.540434
ICP via SVD results:
R = [0.9969452351705235, 0.0598334759429696, -0.05020112774999549;
-0.05932607556034211, 0.9981719680327525, 0.01153858709846634;
0.05079975225724825, -0.008525103530306, 0.9986724727258676]
t = [0.1441598281917405;
-0.06667849447794799;
-0.03009747343724256]
R_inv = [0.9969452351705235, -0.05932607556034211, 0.05079975225724825;
0.0598334759429696, 0.9981719680327525, -0.008525103530306;
-0.05020112774999549, 0.01153858709846634, 0.9986724727258676]
t_inv = [-0.1461462830262246;
0.0576744363694081;
0.03806387978797152]
calling bundle adjustment
iteration= 0 chi2= 1.816112 time= 2.4661e-05 cumTime= 2.4661e-05 edges= 72 schur= 0 lambda= 0.000758 levenbergIter= 1
iteration= 1 chi2= 1.815514 time= 1.8223e-05 cumTime= 4.2884e-05 edges= 72 schur= 0 lambda= 0.000505 levenbergIter= 1
iteration= 2 chi2= 1.815514 time= 1.7451e-05 cumTime= 6.0335e-05 edges= 72 schur= 0 lambda= 0.000337 levenbergIter= 1
iteration= 3 chi2= 1.815514 time= 1.7568e-05 cumTime= 7.7903e-05 edges= 72 schur= 0 lambda= 0.000225 levenbergIter= 1
iteration= 4 chi2= 1.815514 time= 1.7873e-05 cumTime= 9.5776e-05 edges= 72 schur= 0 lambda= 0.000150 levenbergIter= 1
optimization costs time: 0.00031751 seconds.
after optimization:
iteration= 5 chi2= 1.815514 time= 1.7842e-05 cumTime= 0.000113618 edges= 72 schur= 0 lambda= 0.000299 levenbergIter= 1
T=
0.996945 0.0598335 -0.0502011 0.14416
-0.0593261 0.998172 0.0115386 -0.0666785
0.0507998 -0.0085251 0.998672 -0.0300979
0 0 0 1
p1 = [-0.243698, -0.117719, 1.5848]
p2 = [-0.297211, -0.0956614, 1.6558]
P1 = (R*p2+t) = [-0.2409901495364604;
-0.1254270500587826;
1.609221205029395]
p1 = [0.402045, -0.341821, 2.2068]
p2 = [0.378811, -0.262859, 2.2196]
P1 = (R*p2+t) = [0.3946591022539743;
-0.3259188829495218;
2.20803983035825]
p1 = [-0.522843, -0.214436, 1.4956]
p2 = [-0.58581, -0.208584, 1.6052]
P1 = (R*p2+t) = [-0.532923946912698;
-0.2216052393093164;
1.54499035805527]
p1 = [-0.627753, 0.160186, 1.3396]
p2 = [-0.709645, 0.159033, 1.4212]
P1 = (R*p2+t) = [-0.6251478068660965;
0.1505624195985039;
1.351809862638435]
p1 = [0.594266, -0.0256024, 1.5332]
p2 = [0.514795, 0.0391393, 1.5332]
P1 = (R*p2+t) = [0.5827556962439571;
-0.04046060384335358;
1.526884519595548]
952
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
