SLAM练习题(十一)—— G2O实战

2023-11-19

SLAM 学习笔记

写在前面的话:

算是一点小小的感悟吧(* ̄︶ ̄)
估计位姿的方法有线性方法和非线性方法,线性方法就是特征点法中的2D-2D的对极约束,3D-2D的PnP问题。非线性方法有BA优化,它将位姿的估计问题转换成了一个误差关于优化量的最小二乘极值问题,它通过雅克比矩阵指导迭代方向,获得最优的位姿。需要注意的是g2o是一个图优化库,是一个非线性优化库,在使用线性方式求位姿的时候用不到这个库,一般是在使用线性方法(对极约束,PnP……)获得一个误差的估计值后,再用g2o构建图优化问题,设置顶点和边,对位姿进行优化。
但是如果不先求得一个位姿的估计值,直接使用g20构建图优化,将位姿的初始值设为0(或随便设一个初值),进行优化,这样也是可以估计出位姿的,但是迭代次数肯定比先得到一个粗略的位姿再优化要多的多,而且迭代次数与初始值的选取也有很大关系。

特征点法BA优化

分析

以下题目来自计算机视觉life从零开始一起学习SLAM系列
题目: 给定一组世界坐标系下的3D点(p3d.txt)以及它在相机中对应的坐标(p2d.txt),以及相机的内参矩阵。使用
bundle adjustment 方法(g2o库实现)来估计相机的位姿T。初始位姿T为单位矩阵。

(其实就是对3D-2D使用PnP方法得到位姿后,用g2o对该位姿进行优化。)
重投影误差:
在这里插入图片描述
图优化模型:
在这里插入图片描述

代码实现

#include <vector>
#include <fstream>
#include <iostream>
#include <opencv2/core/core.hpp>

#include <Eigen/Core>
#include <Eigen/Geometry>
#include <Eigen/Dense>
#include <g2o/core/base_vertex.h>
#include <g2o/core/base_unary_edge.h>
#include <g2o/core/block_solver.h>
#include <g2o/core/optimization_algorithm_levenberg.h>
#include <g2o/solvers/csparse/linear_solver_csparse.h>
#include <g2o/types/sba/types_six_dof_expmap.h>

using namespace Eigen;

using namespace cv;
using namespace std;


string p3d_file = "./p3d.txt";
string p2d_file = "./p2d.txt";

void bundleAdjustment (
        const vector<Point3f> points_3d,
        const vector<Point2f> points_2d,
        Mat& K );

int main(int argc, char **argv) {


    vector< Point3f > p3d;
    vector< Point2f > p2d;

    Mat K = ( Mat_<double> ( 3,3 ) << 520.9, 0, 325.1, 0, 521.0, 249.7, 0, 0, 1 );

    // 导入3D点和对应的2D点

    ifstream fp3d(p3d_file);
    if (!fp3d){
        cout<< "No p3d.text file" << endl;
        return -1;
    }
    else {
        while (!fp3d.eof()){
            double pt3[3] = {0};
            for (auto &p:pt3) {
                fp3d >> p;
            }
            p3d.push_back(Point3f(pt3[0],pt3[1],pt3[2]));
        }
    }
    ifstream fp2d(p2d_file);
    if (!fp2d){
        cout<< "No p2d.text file" << endl;
        return -1;
    }
    else {
        while (!fp2d.eof()){
            double pt2[2] = {0};
            for (auto &p:pt2) {
                fp2d >> p;
            }
            Point2f p2(pt2[0],pt2[1]);
            p2d.push_back(p2);
        }
    }

    assert(p3d.size() == p2d.size());

    int iterations = 100;
    double cost = 0, lastCost = 0;
    int nPoints = p3d.size();
    cout << "points: " << nPoints << endl;

    bundleAdjustment ( p3d, p2d, K );
    return 0;
}

void bundleAdjustment (
        const vector< Point3f > points_3d,
        const vector< Point2f > points_2d,
        Mat& K   )
{
    // creat g2o
    // new g2o version. Ref:https://www.cnblogs.com/xueyuanaichiyu/p/7921382.html

    // typedef g2o::BlockSolver< g2o::BlockSolverTraits<6,3> > Block;  // pose 维度为 6, landmark 维度为 3
    // 第1步:创建一个线性求解器LinearSolver
    // Block::LinearSolverType* linearSolver = new g2o::LinearSolverCSparse<Block::PoseMatrixType>();

    // 第2步:创建BlockSolver。并用上面定义的线性求解器初始化
    // Block* solver_ptr = new Block (  std::unique_ptr<Block::LinearSolverType>(linearSolver) );

    // 第3步:创建总求解器solver。并从GN, LM, DogLeg 中选一个,再用上述块求解器BlockSolver初始化
    // g2o::OptimizationAlgorithmLevenberg* solver = new g2o::OptimizationAlgorithmLevenberg ( std::unique_ptr<Block>(solver_ptr) );

    // 第4步:创建稀疏优化器
    // g2o::SparseOptimizer optimizer;
    // optimizer.setAlgorithm ( solver );

   // old g2o version
    typedef g2o::BlockSolver<g2o::BlockSolverTraits<6,3>> Block;    //pose维度为6 landmark 维度为3
    Block::LinearSolverType* linearSolver = new g2o::LinearSolverCSparse<Block::PoseMatrixType>();  // 线性求解方程器
    Block* solver_ptr = new Block(linearSolver);    // 矩阵求解器
    g2o::OptimizationAlgorithmLevenberg* solver = new g2o::OptimizationAlgorithmLevenberg(solver_ptr);
    g2o::SparseOptimizer optimizer;
    optimizer.setAlgorithm(solver);

    // 第5步:定义图的顶点和边。并添加到SparseOptimizer中
	
	// ----------------------开始你的代码:设置并添加顶点,初始位姿为单位矩阵
    g2o::VertexSE3Expmap* pose = new g2o::VertexSE3Expmap();    // camera pose
    pose->setId(0);
    pose->setEstimate(g2o::SE3Quat(Eigen::Matrix3d::Identity(), Eigen::Vector3d::Zero()));
    optimizer.addVertex(pose);

    int index =1;
    for (const Point3f p:points_3d)	// 设置空间点 顶点
    {
        g2o::VertexSBAPointXYZ* point = new g2o::VertexSBAPointXYZ();
        point->setId(index++);
        point->setEstimate(Eigen::Vector3d(p.x,p.y,p.z));
        point->setMarginalized(true);	// 空间点必须要边缘化,以进行加速
        optimizer.addVertex(point);
    }

	// ----------------------结束你的代码
	
    // 设置相机内参
    g2o::CameraParameters* camera = new g2o::CameraParameters (
            K.at<double> ( 0,0 ), Eigen::Vector2d ( K.at<double> ( 0,2 ), K.at<double> ( 1,2 ) ), 0);
    camera->setId ( 0 );	// 相机参数 的索引为0
    optimizer.addParameter ( camera );	

    // 设置边
    index = 1;	// 边的id是从1开始的 因为空间点顶点的索引是从1 开始的,方便下面的索引
    for ( const Point2f p:points_2d )
    {
        g2o::EdgeProjectXYZ2UV* edge = new g2o::EdgeProjectXYZ2UV();
        edge->setId ( index );
        edge->setVertex ( 0, dynamic_cast<g2o::VertexSBAPointXYZ*> ( optimizer.vertex ( index ) ) );	// 第一个顶点: 空间点
        edge->setVertex ( 1, pose );	// 第二个顶点:位姿
        edge->setMeasurement ( Eigen::Vector2d ( p.x, p.y ) );  //设置观测值
        edge->setParameterId ( 0,0 );	// 相机的参数
        edge->setInformation ( Eigen::Matrix2d::Identity() );	// 信息矩阵
        optimizer.addEdge ( edge );
        index++;
    }

    // 第6步:设置优化参数,开始执行优化
    optimizer.setVerbose ( true );	// 输出优化信息,包括迭代次数,cost time, 误差等等
    optimizer.initializeOptimization();
    optimizer.optimize ( 100 );

    // 输出优化结果
    cout<<endl<<"after optimization:"<<endl;
    cout<<"T="<<endl<<Eigen::Isometry3d ( pose->estimate() ).matrix() <<endl;
}

这里只考虑了一个位姿。

这是一个最小化重投影问题,使用了g2o 预定好的顶点和边,所以就不用自己再定义了

附CMakeLists.txt:

cmake_minimum_required( VERSION 2.8 )
project( BA )

# 添加c++ 11标准支持
set( CMAKE_CXX_FLAGS "-std=c++11" )

find_package( OpenCV REQUIRED )

# 添加g2o的依赖
# 因为g2o不是常用库,要添加它的findg2o.cmake文件
#set( G2O_ROOT /usr/local/include/g2o )
LIST( APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake_modules )
find_package( G2O REQUIRED )
find_package( CSparse REQUIRED )

# 添加头文件
include_directories( "/usr/include/eigen3"
        ${OpenCV_INCLUDE_DIRS}
        ${G2O_INCLUDE_DIRS}
        ${CSPARSE_INCLUDE_DIR}
        )

add_executable( BA BA-3Dto2D.cpp )

target_link_libraries( BA
        ${OpenCV_LIBS}
        ${CSPARSE_LIBRARY}
        g2o_core g2o_stuff g2o_types_sba g2o_csparse_extension
        )

雅克比矩阵

参考十四讲194页的推导过程:
在这里插入图片描述
误差相对于空间点P的导数:
在这里插入图片描述
稍有差别的是,g2o相机里用f统一描述fx,fy,并且李代数定义顺序不同(g2o是旋转在前,平移在后,书中是平移在前,旋转在后),所以矩阵前三列和后三列与书中的定义是相反的,此外都一致。

直接法BA优化

分析

题目: 用直接法Bundle Adjustment 估计相机位姿。
链接:https://pan.baidu.com/s/105a9WYr0pPOd691leyanjA 提取码:bjzi
知识点: g2o的边

直接法参考14讲上的介绍,主要是灰度不变原则,在这里我们使用一个4*4的窗口内像素的灰度值的误差构建优化目标函数。优化变量为位姿T,三维点P_i,误差为像素点的误差16×1;

在这里插入图片描述

代码实现

这次考虑了三个位姿,n个空间点。

g2o本身没有计算光度误差的边,我们需要自己定义。

/****************************
 * 题目:用直接法Bundle Adjustment 估计相机位姿,已经给定3张图,poses.txt,points.txt文件。
 *
 * 本程序学习目标:
 * 熟悉g2o库顶点和边的定义。
 * 
 * 本代码框架来自高翔,有改动
 * 公众号:计算机视觉life。发布于公众号旗下知识星球:从零开始学习SLAM
 * 时间:2019.03
****************************/

#include <iostream>

using namespace std;

#include <g2o/core/base_unary_edge.h>
#include <g2o/core/base_binary_edge.h>
#include <g2o/core/base_vertex.h>
#include <g2o/core/block_solver.h>
#include <g2o/core/optimization_algorithm_levenberg.h>
#include <g2o/solvers/dense/linear_solver_dense.h>
#include <g2o/core/robust_kernel.h>
#include <g2o/core/robust_kernel_impl.h>
#include <g2o/types/sba/types_six_dof_expmap.h>

#include <Eigen/Core>
#include <sophus/se3.h>
#include <opencv2/opencv.hpp>

#include <pangolin/pangolin.h>
#include <boost/format.hpp>
using namespace Sophus;
using namespace pangolin;
using namespace g2o;

#define DEBUG false

typedef vector<Sophus::SE3, Eigen::aligned_allocator<Sophus::SE3>> VecSE3;
typedef vector<Eigen::Vector3d, Eigen::aligned_allocator<Eigen::Vector3d>> VecVec3d;

// global variables
string pose_file = "./poses.txt";
string points_file = "./points.txt";

// intrinsics
float fx = 277.34;
float fy = 291.402;
float cx = 312.234;
float cy = 239.777;
// 判断当前像素点是否在当前帧内
inline bool bInImage(float u, float v, int w, int h)
{
    if(u>=0 && u<w && v>=0 && v<h)
        return true;
    else
        return false;
}
// bilinear interpolation 通过双线性差值获得浮点坐标对应插值后的像素值
inline float GetPixelValue(const cv::Mat &img, float x, float y) {
    uchar *data = &img.data[int(y) * img.step + int(x)];	// 取得该像素值的data位置 step为Mat中一行所占的字节数
    float xx = x - floor(x);
    float yy = y - floor(y);
    return float(
            (1 - xx) * (1 - yy) * data[0] +
            xx * (1 - yy) * data[1] +
            (1 - xx) * yy * data[img.step] +
            xx * yy * data[img.step + 1]
    );
}

// g2o vertex that use sophus::SE3 as pose
class VertexSophus : public g2o::BaseVertex<6, Sophus::SE3> {
public:
    EIGEN_MAKE_ALIGNED_OPERATOR_NEW

    VertexSophus() {}

    ~VertexSophus() {}

    bool read(std::istream &is) {}

    bool write(std::ostream &os) const {}

    virtual void setToOriginImpl() {
        _estimate = Sophus::SE3();
    }

    virtual void oplusImpl(const double *update_) {
        Eigen::Map<const Eigen::Matrix<double, 6, 1>> update(update_);
        setEstimate(Sophus::SE3::exp(update) * estimate());
    }
};

// TODO edge of projection error, implement it
// 16x1 error, which is the errors in patch
typedef Eigen::Matrix<double,16,1> Vector16d;
class EdgeDirectProjection : public g2o::BaseBinaryEdge<16, Vector16d, g2o::VertexSBAPointXYZ, VertexSophus> {
public:
    EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
	// *color是float类型的指针,指向16维的灰度值数组,也就是当前空间点(邻近窗口)的灰度信息,target是当前位姿对应的图片
    EdgeDirectProjection(float *color, cv::Mat &target) {
        this->origColor = color;
        this->targetImg = target;
        this->w = targetImg.cols;
        this->h = targetImg.rows;

    }

    ~EdgeDirectProjection() {}

    virtual void computeError() override {

        // ----------------  开始你的代码 ----------------------//
		// 参考十四讲中直接法BA部分
		const VertexSBAPointXYZ* vertexPw = static_cast<const VertexSBAPointXYZ*>(vertex(0));// 要注意索引值 只针对当前边的两个顶点
        const VertexSophus* vertexTcw = static_cast<const VertexSophus*>(vertex(1));
        Vector3d Pc = vertexTcw->estimate()*vertexPw->estimate();
        float u = Pc[0]/Pc[2]*fx + cx;
        float v = Pc[1]/Pc[2]*fy + cy;
        if(!bInImage(u-3, v-3, w, h)||!bInImage(u+2,v+2,w,h))
        {
            // 一般g2o默认只处理level为0的边,对于不再当前帧的路标,将其level设为1
            this->setLevel(1);
            for(int n =0;n<16;n++)
            {
                _error[n]=0;// 初始化误差为0
            }
        }else{
            for(int i =-2;i<2;i++){
                for(int j = -2;j<2;j++){
                    int num = 4*i + j + 10;// 0 1 2 3 4…… 15
                    // error = 真实观测值 - 计算观测值 ,u i 设置的妙啊
                    _error[num]=origColor[num]-GetPixelValue(targetImg, u+i, v+j);
                }
            }
        }
        // ----------------  结束你的代码 ----------------------//
    }
	
	// 下面函数不用自己写了,但是要能看懂
    virtual void linearizeOplus() override
    {
        if(level()==1)
        {
            _jacobianOplusXi = Matrix<double,16,3>::Zero();
            _jacobianOplusXj = Matrix<double,16,6>::Zero();
            return;
        }
        const VertexSBAPointXYZ* vertexPw = static_cast<const VertexSBAPointXYZ* >(vertex(0));
        const VertexSophus* vertexTcw = static_cast<const VertexSophus* >(vertex(1));
        Vector3d Pc = vertexTcw->estimate() * vertexPw->estimate();
        float x = Pc[0];
        float y = Pc[1];
        float z = Pc[2];
        float inv_z = 1.0/z;
        float inv_z2 = inv_z * inv_z;
        float u = x * inv_z * fx + cx;
        float v = y * inv_z * fy + cy;
		// 投影方程关于相机坐标系下三维点的导数
        Matrix<double,2,3> J_Puv_Pc;
        J_Puv_Pc(0,0) = fx * inv_z;
        J_Puv_Pc(0,1) = 0;
        J_Puv_Pc(0,2) = -fx * x * inv_z2;
        J_Puv_Pc(1,0) = 0;
        J_Puv_Pc(1,1) = fy * inv_z;
        J_Puv_Pc(1,2) = -fy * y * inv_z2;

		// 相机坐标系下三维点关于位姿扰动量的导数
        Matrix<double,3,6> J_Pc_kesi = Matrix<double,3,6>::Zero();
        J_Pc_kesi(0,0) = 1;
        J_Pc_kesi(0,4) = z;
        J_Pc_kesi(0,5) = -y;
        J_Pc_kesi(1,1) = 1;
        J_Pc_kesi(1,3) = -z;
        J_Pc_kesi(1,5) = x;
        J_Pc_kesi(2,2) = 1;
        J_Pc_kesi(2,3) = y;
        J_Pc_kesi(2,4) = -x;

        Matrix<double,1,2> J_I_Puv; // 像素梯度
        for(int i = -2; i<2; i++)
            for(int j = -2; j<2; j++) {
                int num = 4 * i + j + 10;
                J_I_Puv(0,0) = (GetPixelValue(targetImg,u+i+1,v+j) - GetPixelValue(targetImg,u+i-1,v+j))/2;
                J_I_Puv(0,1) = (GetPixelValue(targetImg,u+i,v+j+1) - GetPixelValue(targetImg,u+i,v+j-1))/2;
                // 光度误差对空间点的导数
                _jacobianOplusXi.block<1,3>(num,0) = -J_I_Puv * J_Puv_Pc * vertexTcw->estimate().rotation_matrix();
                // 光度误差对位姿的导数
                _jacobianOplusXj.block<1,6>(num,0) = -J_I_Puv * J_Puv_Pc * J_Pc_kesi;
            }
    }
    virtual bool read(istream &in) {}

    virtual bool write(ostream &out) const {}

private:
    cv::Mat targetImg;  // the target image
    float *origColor = nullptr;   // 16 floats, the color of this point
    int w;
    int h;
};

// plot the poses and points for you, need pangolin
void Draw(const VecSE3 &poses, const VecVec3d &points);

int main(int argc, char **argv) {

    // read poses and points
    VecSE3 poses;
    VecVec3d points;
    ifstream fin(pose_file);

    while (!fin.eof()) {
        double timestamp = 0;
        fin >> timestamp;
        if (timestamp == 0) break;
        double data[7];
        for (auto &d: data) fin >> d;
        poses.push_back(Sophus::SE3(
                Eigen::Quaterniond(data[6], data[3], data[4], data[5]),
                Eigen::Vector3d(data[0], data[1], data[2])
        ));
        if (!fin.good()) break;
    }
    fin.close();

    vector<float *> color;
    fin.open(points_file);
    while (!fin.eof()) {
        double xyz[3] = {0};
        for (int i = 0; i < 3; i++) fin >> xyz[i];
        if (xyz[0] == 0) break;
        points.push_back(Eigen::Vector3d(xyz[0], xyz[1], xyz[2]));
        float *c = new float[16];
        for (int i = 0; i < 16; i++) fin >> c[i];
        color.push_back(c);

        if (fin.good() == false) break;
    }
    fin.close();

    cout << "poses: " << poses.size() << ", points: " << points.size() << ", color: "<<color.size()<<endl;

    // read images
    vector<cv::Mat> images;
    boost::format fmt("../%d.png");
    for (int i = 0; i < 3; i++) {
        images.push_back(cv::imread((fmt % i).str(), 0));
    }
    cout<<"images: "<<images.size()<<endl;

    // typedef g2o::BlockSolver< g2o::BlockSolverTraits<6,3> > Block;  // pose 维度为 6, landmark 维度为 3
    // 第1步:创建一个线性求解器LinearSolver
    // Block::LinearSolverType* linearSolver = new g2o::LinearSolverDense<Block::PoseMatrixType>();

    // 第2步:创建BlockSolver。并用上面定义的线性求解器初始化
    // Block* solver_ptr = new Block (  std::unique_ptr<Block::LinearSolverType>(linearSolver) );

    // 第3步:创建总求解器solver。并从GN, LM, DogLeg 中选一个,再用上述块求解器BlockSolver初始化
    // g2o::OptimizationAlgorithmLevenberg* solver = new g2o::OptimizationAlgorithmLevenberg ( std::unique_ptr<Block>(solver_ptr) );

    // 第4步:创建稀疏优化器
    // g2o::SparseOptimizer optimizer;
    // optimizer.setAlgorithm ( solver );

   // old g2o version
   typedef g2o::BlockSolver< g2o::BlockSolverTraits<6,3> > Block;  // pose 维度为 6, landmark 维度为 3
   Block::LinearSolverType* linearSolver = new g2o::LinearSolverDense<Block::PoseMatrixType>(); // 线性方程求解器
   Block* solver_ptr = new Block ( linearSolver );     // 矩阵块求解器
   g2o::OptimizationAlgorithmLevenberg* solver = new g2o::OptimizationAlgorithmLevenberg ( solver_ptr );
   g2o::SparseOptimizer optimizer;
   optimizer.setAlgorithm ( solver );

    // 第5步:添加顶点和边
    // ----------------  开始你的代码 ----------------------//
	// 参考十四讲中直接法BA部分
	for (int i = 0; i < points.size(); i++)
    {
        VertexSBAPointXYZ* vertexPw = new VertexSBAPointXYZ();
        vertexPw->setId(i);
        vertexPw->setEstimate(points[i]);
        vertexPw->setMarginalized(true);
        optimizer.addVertex(vertexPw);
    }
    for (int i = 0; i < poses.size(); i++)
    {
        VertexSophus* vertexTcw  = new VertexSophus();
        vertexTcw->setId(i+points.size());
        vertexTcw->setEstimate(poses[i]);
        optimizer.addVertex(vertexTcw);
    }
    
    for (int i = 0; i < poses.size(); i++)
    {
        for (int j = 0; j < points.size(); j++)
        {	
            // color里面存的是第j个点周围4*4=16个像素点的灰度值的指针
            EdgeDirectProjection* edge = new EdgeDirectProjection(color[j], images[i]);
            edge->setVertex(0, dynamic_cast<VertexSBAPointXYZ*>(optimizer.vertex(j)));
            edge->setVertex(1, dynamic_cast<VertexSophus*>(optimizer.vertex(i+points.size())));
            edge->setInformation(Matrix<double,16,16>::Identity());
            RobustKernelHuber* rk = new RobustKernelHuber;	// 鲁棒核函数
            rk->setDelta(1.0);
            edge->setRobustKernel(rk);
            optimizer.addEdge(edge);
        }
        
    }
    
    
    // ----------------  结束你的代码 ----------------------//

    // 第6步:执行优化
    optimizer.initializeOptimization(0);
    optimizer.optimize(200);

    // 从optimizer中获取结果
    for(int c = 0; c < poses.size(); c++)
        for(int p = 0; p < points.size(); p++) {
            Vector3d Pw = dynamic_cast<VertexSBAPointXYZ*>(optimizer.vertex(p))->estimate();
            points[p] = Pw;
            SE3 Tcw = dynamic_cast<VertexSophus*>(optimizer.vertex(c + points.size()))->estimate();
            poses[c] = Tcw;
        }


    // plot the optimized points and poses
    Draw(poses, points);

    // delete color data
    for (auto &c: color) delete[] c;
    return 0;
}

// 画图函数,无需关注
void Draw(const VecSE3 &poses, const VecVec3d &points) {
    if (poses.empty() || points.empty()) {
        cerr << "parameter is empty!" << endl;
        return;
    }
    cout<<"Draw poses: "<<poses.size()<<", points: "<<points.size()<<endl;
    // create pangolin window and plot the trajectory
    pangolin::CreateWindowAndBind("Trajectory Viewer", 1024, 768);
    glEnable(GL_DEPTH_TEST);
    glEnable(GL_BLEND);
    glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);

    pangolin::OpenGlRenderState s_cam(
            pangolin::ProjectionMatrix(1024, 768, 500, 500, 512, 389, 0.1, 1000),
            pangolin::ModelViewLookAt(0, -0.1, -1.8, 0, 0, 0, 0.0, -1.0, 0.0)
    );

    pangolin::View &d_cam = pangolin::CreateDisplay()
            .SetBounds(0.0, 1.0, pangolin::Attach::Pix(175), 1.0, -1024.0f / 768.0f)
            .SetHandler(new pangolin::Handler3D(s_cam));


    while (pangolin::ShouldQuit() == false) {
        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

        d_cam.Activate(s_cam);
        glClearColor(0.0f, 0.0f, 0.0f, 0.0f);

        // draw poses
        float sz = 0.1;
        int width = 640, height = 480;
        for (auto &Tcw: poses) {
            glPushMatrix();
            Sophus::Matrix4f m = Tcw.inverse().matrix().cast<float>();
            glMultMatrixf((GLfloat *) m.data());
            glColor3f(1, 0, 0);
            glLineWidth(2);
            glBegin(GL_LINES);
            glVertex3f(0, 0, 0);
            glVertex3f(sz * (0 - cx) / fx, sz * (0 - cy) / fy, sz);
            glVertex3f(0, 0, 0);
            glVertex3f(sz * (0 - cx) / fx, sz * (height - 1 - cy) / fy, sz);
            glVertex3f(0, 0, 0);
            glVertex3f(sz * (width - 1 - cx) / fx, sz * (height - 1 - cy) / fy, sz);
            glVertex3f(0, 0, 0);
            glVertex3f(sz * (width - 1 - cx) / fx, sz * (0 - cy) / fy, sz);
            glVertex3f(sz * (width - 1 - cx) / fx, sz * (0 - cy) / fy, sz);
            glVertex3f(sz * (width - 1 - cx) / fx, sz * (height - 1 - cy) / fy, sz);
            glVertex3f(sz * (width - 1 - cx) / fx, sz * (height - 1 - cy) / fy, sz);
            glVertex3f(sz * (0 - cx) / fx, sz * (height - 1 - cy) / fy, sz);
            glVertex3f(sz * (0 - cx) / fx, sz * (height - 1 - cy) / fy, sz);
            glVertex3f(sz * (0 - cx) / fx, sz * (0 - cy) / fy, sz);
            glVertex3f(sz * (0 - cx) / fx, sz * (0 - cy) / fy, sz);
            glVertex3f(sz * (width - 1 - cx) / fx, sz * (0 - cy) / fy, sz);
            glEnd();
            glPopMatrix();
        }

        // points
        glPointSize(2);
        glBegin(GL_POINTS);
        for (size_t i = 0; i < points.size(); i++) {
            glColor3f(0.0, points[i][2]/4, 1.0-points[i][2]/4);
            glVertex3d(points[i][0], points[i][1], points[i][2]);
        }
        glEnd();

        pangolin::FinishFrame();
        usleep(5000);   // sleep 5 ms
    }
}


附CMakeLists.txt

cmake_minimum_required( VERSION 2.8 )
project( directBA )

# 添加c++ 11标准支持
set( CMAKE_CXX_FLAGS "-std=c++11 -O3" )
list( APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake_modules )
set(G2O_LIBS g2o_cli g2o_ext_freeglut_minimal g2o_simulator g2o_solver_slam2d_linear
        g2o_types_icp g2o_types_slam2d g2o_types_sba g2o_types_slam3d g2o_core g2o_interface
        g2o_solver_csparse g2o_solver_structure_only g2o_csparse_extension g2o_opengl_helper g2o_solver_dense
        g2o_stuff g2o_types_sclam2d g2o_parser g2o_solver_pcg g2o_types_data g2o_types_sim3 cxsparse )
# 添加cmake模块以使用g2o库
list( APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake_modules )
find_package( G2O REQUIRED)
find_package( Eigen3 REQUIRED)
find_package( OpenCV REQUIRED )beijign
find_package( Sophus REQUIRED )
find_package( Pangolin REQUIRED)

# 添加头文件
include_directories( "/usr/include/eigen3" )
include_directories( ${OpenCV_INCLUDE_DIRS} )
include_directories( ${Sophus_INCLUDE_DIRS} )
include_directories( ${Pangolin_INCLUDE_DIRS})
include_directories( ${G2O_INCLUDE_DIRS})

add_executable( directBA directBA.cpp )

target_link_libraries( directBA ${OpenCV_LIBS} ${G2O_LIBS})
target_link_libraries( directBA ${Sophus_LIBRARIES} )
target_link_libraries( directBA   /home/yan/3rdparty/Sophus/build/libSophus.so )
target_link_libraries( directBA ${Pangolin_LIBRARIES})

讨论

双线性插值

通过双线性差值获得浮点坐标对应插值后的像素值

背景:把空间点point(x,y,z)经过投影计算得到的图像坐标(u,v)(即下面函数的输入参数x,y)不一定是整数,我们需要对该浮点数进行插值处理。

// bilinear interpolation 通过双线性差值获得浮点坐标对应插值后的像素值
inline float GetPixelValue(const cv::Mat &img, float x, float y) {
    uchar *data = &img.data[int(y) * img.step + int(x)];	// 取得该像素值的data位置 step为Mat中一行所占的字节数
    float xx = x - floor(x);	// xx 为x的小数部分
    float yy = y - floor(y);	// yy为y的小数部分。
    return float(
            (1 - xx) * (1 - yy) * data[0] +
            xx * (1 - yy) * data[1] +
            (1 - xx) * yy * data[img.step] +
            xx * yy * data[img.step + 1]
    );
}

在图像的仿射变换(图像的放大缩小)中,很多地方需要用到插值运算,常见的插值运算包括最邻近插值,双线性插值,双三次插值,兰索思插值等方法,OpenCV提供了很多方法,其中,双线性插值由于折中的插值效果和运算速度,用的最多。

双线性插值算法:

对于一个目的像素,设置坐标通过反向变换得到的浮点坐标为(i+u,j+v) (其中i、j均为浮点坐标的整数部分u、v为浮点坐标的小数部分,是取值[0,1)区间的浮点数),则这个像素得值 f(i+u,j+v) 可由原图像中坐标为 (i,j)、(i+1,j)、(i,j+1)、(i+1,j+1)所对应的周围四个像素的值决定,即:f(i+u,j+v) = (1-u)(1-v)f(i,j) + (1-u)vf(i,j+1) + u(1-v)f(i+1,j) + uvf(i+1,j+1),其中f(i,j)表示源图像(i,j)处的的像素值,以此类推。

对应于程序里面:

data[0]为f(i,j),data[1]为f(i,j+1),data[img.step]为f(i+1,j), data[img.step+1]为f(i+1,j+1)

具体推导过程可参考:图像处理之双线性插值法

鲁棒核函数huber

参考维基百科:huber
在这里插入图片描述

当误差e大于某个阈值后,函数增长由二次形式变成了一次形式,相当于限制了梯度的最大值,同时Huber核函数又是光滑的,可以很方便的求导。参考Slam十四讲255页。

雅克比矩阵

推导过程参考14讲194页。

首先是光度误差关于位姿扰动量的导数,这里的q相当于特征点法中的P’,表示相机坐标系下的点。


J = − ∂ I 2 ∂ u ∂ u ∂ q ∂ q ∂ δ ξ = − ∂ I 2 ∂ u ∂ u ∂ δ ξ \large J=-\frac{\partial{I_2}}{\partial{u}}\frac{\partial{u}}{\partial{q}}\frac{\partial{q}}{\partial{\delta\xi}}=-\frac{\partial{I_2}}{\partial{u}}\frac{\partial{u}}{\partial{\delta\xi}} J=uI2quδξq=uI2δξu

  • ∂ I 2 ∂ u \LARGE \frac{\partial{I_2}}{\partial{u}} uI2是u处的像素梯度
    在上面 的代码中,使用了中值定理来计算梯度:
J_I_Puv(0,0) = (GetPixelValue(targetImg,u+i+1,v+j) - GetPixelValue(targetImg,u+i-1,v+j))/2;
                J_I_Puv(0,1) = (GetPixelValue(targetImg,u+i,v+j+1) - GetPixelValue(targetImg,u+i,v+j-1))/2;

也就是利用u处的前一个点和后一个点的像素均值来代替u处的像素梯度。

  • ∂ u ∂ q \LARGE \frac{\partial{u}}{\partial{q}} qu为投影方程关于相机坐标系下的三维点的导数

也就是像素坐标 [ u v ] \begin{bmatrix}u \\v \end{bmatrix} [uv]关于相机坐标系下三维点q [ X ’ Y ′ Z ′ ] \begin{bmatrix} X’&Y'&Z'\end{bmatrix} [XYZ]的导数:
∂ u ∂ q = [ f x Z ′ 0 − f x X ′ Z ′ 2 0 f y Z ′ − f y Y ′ Z ′ 2 ] \large \frac{\partial u}{\partial q}=\begin{bmatrix}\frac{f_x}{Z'}&0&-\frac{f_xX'}{Z'^2}\\ 0&\frac{f_y}{Z'}&-\frac{f_yY'}{Z'^2}\end{bmatrix} qu=Zfx00ZfyZ2fxXZ2fyY

  • ∂ q ∂ δ ξ \LARGE \frac{\partial q}{\partial \delta\xi} δξq为变换后的三维点(相机坐标系下的点)对位姿扰动量的导数。

在求左乘扰动的时候,已经证明:
∂ q ∂ δ ξ = [ I , − q ∧ ] \large \frac{\partial q}{\partial \delta \xi}=\begin{bmatrix}I,&-q^{\wedge}\end{bmatrix} δξq=[Iq]
其中
q ∧ = [ 0 − Z ′ Y ′ Z ′ 0 − X ′ − Y ′ X ′ 0 ] \large q^{\wedge}=\begin{bmatrix}0&-Z'&Y'\\ Z'&0&-X'\\ -Y'&X'&0 \end{bmatrix} q=0ZYZ0XYX0
后两项只与三维点P有关,而与图像无关,经常将这两项合并在一起:
∂ u ∂ δ ξ = [ f x Z ′ 0 − f x X ′ Z ′ 2 − f x X ′ Y ′ Z ′ 2 f x + f x X ′ 2 Z ′ 2 − f x Y ′ Z ′ 0 f y Z ′ − f y Y ′ Z ′ 2 − f y − f y Y ′ 2 Z ′ 2 f y X ′ Y ′ Z ′ 2 f y X ′ Z ′ ] \large \frac{\partial{u}}{\partial{\delta\xi}}=\begin{bmatrix}\frac{f_x}{Z'}&0&-\frac{f_xX'}{Z'^2} &-\frac{f_xX'Y'}{Z'^2}&f_x+\frac{f_xX'^2}{Z'^2}&-\frac{f_xY'}{Z'}\\ 0&\frac{f_y}{Z'}&-\frac{f_yY'}{Z'^2}&-f_y-\frac{f_yY'^2}{Z'^2}&\frac{f_yX'Y'}{Z'^2}&\frac{f_yX'}{Z'}\end{bmatrix} δξu=Zfx00ZfyZ2fxXZ2fyYZ2fxXYfyZ2fyY2fx+Z2fxX2Z2fyXYZfxYZfyX
所以光度误差相对于李代数的雅克比矩阵为:
J = − ∂ I 2 ∂ u ∂ u ∂ δ ξ \large J=-\frac{\partial{I_2}}{\partial{u}}\frac{\partial{u}}{\partial{\delta\xi}} J=uI2δξu
光度误差关于空间点的导数

和特征点法中的求法差不多:
∂ I ∂ P = ∂ I ∂ u ∂ u ∂ q ∂ q ∂ P \large \frac{\partial I}{\partial P}=\frac{\partial I}{\partial u}\frac{\partial u}{\partial q}\frac{\partial q}{\partial P} PI=uIquPq

前两项分别是像素梯度和投影坐标关于Pc的偏导数,第三项是Pc关于Pw的偏导数,又因为 q = e x p ( ξ ∧ ) P = R P + t q=exp(\xi^{\wedge})P=RP+t q=exp(ξ)P=RP+t,所以q对P的求导后只剩下R,于是:
∂ I ∂ P = ∂ I ∂ u ∂ u ∂ q R \large \frac{\partial I}{\partial P}=\frac{\partial I}{\partial u}\frac{\partial u}{\partial q}R PI=uIquR
在代码中为:

_jacobianOplusXi.block<1,3>(num,0) = -J_I_Puv * J_Puv_Pc * vertexTcw->estimate().rotation_matrix();
                _jacobianOplusXj.block<1,6>(num,0) = -J_I_Puv * J_Puv_Pc * J_Pc_kesi;

_jacobianOplusXi为光度误差对空间点的导数,里面存放16行1*3的矩阵。

_jacobianOplusXj光度误差对位姿(扰动量)的导数,里面存放16行1*6的矩阵。

参考:
视觉SLAM14讲
从零开始一起学习SLAM | 掌握g2o顶点编程套路

从零开始一起学习SLAM | 掌握g2o边的代码套路

本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)

SLAM练习题(十一)—— G2O实战 的相关文章

  • 给定焦点和相机主中心,如何创建透视投影矩阵

    我设法使用 OpenCV 获取相机的内在和外在参数 因此我有 fx fy cx 和 cy 我还有屏幕 图像的宽度和高度 但是如何根据这些参数创建 OpenGL 透视投影矩阵呢 glFrustum 展示了如何在给定 Z 近 Z 远以及图像宽度
  • 如何在 OpenCV 中绘制图像的 3D 直方图

    更新 我找到更多例子 我现在可以做到 我可以在 3d 中绘制多个直方图吗 https stackoverflow com questions 35210337 can i plot several histograms in 3d 我知道这
  • 如何忽略 OpenCV python 中的内部黑色轮廓?

    在上图中 我只需要以下轮廓 8 7 0 2 4 里面每隔一个轮廓都是黑色的空盒子 有没有办法使用 cv2 RETR TREE 自动提取此类轮廓 cv2 RETR EXTERNAL 将忽略我实际需要的 2 和 4 gt gt heirarch
  • 在 python 中使用 opencv 进行 Blob 过滤

    需要从图像中检测红色并根据屏幕尺寸获取坐标 使用掩模获取图像中具有红色的部分 将其转换为 BW 对其应用高斯滤波器 最终图像具有较小的主体 我需要将其删除并获取其余部分的坐标 我尝试了 SimpleBlobDetector 但没有帮助 这是
  • Android Studio、OpenCV、C++、jni、NDK - 无法配置?

    我正在尝试使用 jni 和 Android Studio NDK 在 Android Studio 中进行 OpenCV C 调用 这些是我一直在工作的参考资料 https github com opencv opencv tree mas
  • opencv中以下代码行的含义是什么?

    该代码行的含义是什么以及如何将此代码转换为 javacv gray Scalar all 255 这是与此代码行相关的完整代码 Mat src imread in jpg gray cvtColor src gray CV BGR2GRAY
  • 使用额外模块构建 opencv 错误 ocv_download

    我尝试使用 Extramodules 构建 openCV 如果我尝试使用 cmake 配置项目文件 Windows 10 上的 vc14 x32 则会出现错误 opencv contrib master modules xfeatures2
  • 如何在客户端-服务器模型(OpenCV C++)中通过套接字传输 cv::VideoCapture 帧?

    我通过套接字翻译视频 我在另一端看到视频流 但无法接收视频 我的视频文件是空的 我认为问题可能是接收端转换视频格式错误 我应该如何在服务器端对视频进行编码以及如何在客户端对视频进行解码 我使用 Linux Ubuntu Mate OpenC
  • 无法保存从网络摄像头捕获的图像(OpenCV 2.3 的 imwrite 编译错误)

    我正在使用 OpenCV 2 3 制作简单的网络摄像头程序 但遇到了编译错误 任何想法将不胜感激 编译后 我在 imwrite 处收到以下错误 在下面代码的 read 函数中 这个样本 https code ros org svn open
  • 如何在Python中删除图像的背景

    我有一个包含全角人类图像的数据集 我想删除这些图像中的所有背景 只留下全角人物 我的问题 有没有Python代码可以做到这一点 我是否需要每次都指定人员对象的坐标 这是使用 Python OpenCV 的一种方法 读取输入 转换为灰色 阈值
  • 使用 Python 和 OpenCV 改善图像歪斜校正

    我生成的用于检测和纠正偏差的代码给了我不一致的结果 我目前正在开展一个项目 该项目利用图像上的 OCR 文本提取 通过 Python 和 OpenCV 因此如果需要准确的结果 消除倾斜是关键 我的代码使用cv2 minAreaRect来检测
  • C++/OpenCV - 用于视频稳定的卡尔曼滤波器

    我尝试使用卡尔曼滤波器稳定视频以进行平滑 但我有一些问题 每次 我都有两帧 一帧是当前帧 另一帧是当前帧 这是我的工作流程 计算 goodFeaturesToTrack 使用 calcOpticalFlowPyrLK 计算光流 只保留优点
  • 使用python opencv从zip加载图像

    我能够成功地从 zip 加载图像 with zipfile ZipFile test zip r as zfile data zfile read test jpg how to open this using imread or imde
  • 在白色背景上将透明 PNG 保存为 JPEG

    假设我有一张 BGRA 图像numpy数组看起来非常像这样 233 228 230 128 233 228 230 128 233 228 230 0 164 160 159 65 199 197 196 65 255 255 254 12
  • 在 Mac OS X 上使用 pip 安装 pyopencv

    我正在尝试在 OS X Mountain Lion 中使用 pip 安装 pyopencv 但导入 setuptools 失败 以下是我的作品 setuptools中的 库 是什么 我以前没见过 我已经通过 homebrew 和其他东西安装
  • 将图像从 CV_64F 转换为 CV_8U

    我想转换类型的图像CV 64FC1 to CV 8UC1在 Python 中使用 OpenCV 在 C 中 使用convertTo函数 我们可以使用以下代码片段轻松转换图像类型 image convertTo image CV 8UC1 我
  • 使用opencv在灰度图像中绘制彩色矩形

    是否可以使用opencv在灰度图像中绘制彩色矩形 我尝试了几种方法 但整个图像要么变成灰度 要么变成 RGB 您不能获得混合的灰色和彩色图像 你可以看看有没有办法在同一图像opencv C 上同时拥有灰度和RGB像素 https stack
  • OpenCV Python - 如何实现RANSAC来检测直线?

    我正在尝试检测包含道路的图像上的线条 使用高斯平滑和 Canny 边缘检测 我在尝试实现 RANSAC 时遇到了困难 我基本上不知道如何去做 我可以获得一个粗略的想法或一个带有实现 RANSAC 的随机图像 foo png 的工作代码吗 N
  • 使用高斯混合模型进行皮肤检测

    我正在根据以下进行皮肤检测算法本文 http www cc gatech edu rehg Papers SkinDetect IJCV lowres pdf 第 21 页有两个模型 高斯皮肤混合模型和非皮肤颜色模型 第一个皮肤检测模型效果
  • 将 openCV 矩阵转换为向量

    看起来很容易 毕竟我们知道 std 或 openCV 向量可以轻松转换为矩阵 如下所示 vector

随机推荐

  • vue聊天页面在进入之后信息滑动到底部位置

    这是需要实现的目标 怎么做到进入到页面之后就滑动到底部 借助两个属性 scrollHeigh 该属性是指 元素中内容 的高度 图中的意思就是全部信息所占用的总高度 scrollTop 指的是 元素中的内容 超出 元素上边界 的那部分的高度
  • Python数据可视化详解(3/5)--------patch,饼图,柱状图和点图

    水平或垂直的单条柱状图 如图 上代码 import matplotlib pyplot as plt import numpy as np fig axes plt subplots 2 1 x 1 2 3 4 5 6 data 5 4 1
  • Redis高级客户端Lettuce详解

    前提 Lettuce是一个Redis的Java驱动包 初识她的时候是使用RedisTemplate的时候遇到点问题Debug到底层的一些源码 发现spring data redis的驱动包在某个版本之后替换为Lettuce Lettuce翻
  • 测试用例的基本要素 && properties配置文件 && 测试用例的基本要素 && SpringMVC背景知识 && 按照开发阶段划分测试类型

    第 1 题 简答题 题目名称 设计测试用例 登陆页面 题目内容 请你针对登陆场景来设计测试用例 记得用脑图哦 第 2 题 简答题 题目名称 测试用例设计 淘宝购物车 题目内容 设计淘宝购物车的测试用例 尽可能多的来设计测试用例 若想不出来了
  • DataGrip配置设定

    提示 文章写完后 目录可以自动生成 如何生成可参考右边的帮助文档 文章目录 前言 一 下载安装 二 使用步骤 1 配置成中文 2 导入旧的设定 3 常用操作 总结 前言 提示 这里可以添加本文要记录的大概内容 更换电脑 重装DataGrip
  • java通过word模板生成导出word

    java通过word模板生成导出word 生成word方法类 public class BokeWordUtils 根据模板生成新word文档 判断表格是需要替换还是需要插入 判断逻辑有 为替换 表格无 为插入 param inputUrl
  • MYSQL数据库转化成H2数据库

    mysql数据库数据导成CSV文件 H2 数据库导入CSV文件 代码如下 insert into tableName select from csvread c import data csv
  • DBCP 应用的总结(一)

    DBCP 应用的总结 1 Jar包下载 首先 下载必须的jar包 dbcp包 http jakarta apache org commons dbcp pool包 http jakarta apache org commons pool 如
  • POJ 2676 Sudoku 数独(dfs)

    POJ 2676 Sudoku 数独 dfs Sudoku is a very simple task A square table with 9 rows and 9 columns is divided to 9 smaller squ
  • Java jvm 内存溢出分析

    1 如何分析jvm内存溢出呢 我们经常用visualVm监控Jvm的内存 cpu 线程的使用情况 通常可以根据内存不断增长来判断内存是否存在不释放 但是我们不可能时时盯着去看 这里涉及jvm堆内存配置 堆内存参数配置和调优会在其他章节编写
  • 逆波兰表达式求值(C语言实现)

    实验项目 从文本文件输入任意一个语法正确的 中缀 表达式 显示并保存该表达式 利用栈结构 把上述 中缀 表达式转换成后缀表达式 并显示栈的状态变化过程和所得到的后缀表达式 利用栈结构 对上述后缀表达式进行求值 并显示栈的状态变化过程和最终结
  • 【mcuclub】数据存储AT24C02

    一 实物图 二 原理图 编号 名称 功能 1 A0 地址输入 A2 A1和A0是器件地址输入引脚 24C02 32 64使用A2 A1和AO输入引脚作为硬件地址 总线上可同时级联8个24C02 32 64器件 24C04使用A2和A1输入引
  • 身体语言密码

    本书来自www bookdown com cn免费txt小说下载站 更多更新免费电子书请关注www bookdown com cn 掀起职场白领流行新风尚 身体语言密码 全文 一个无心的眼神 一个不经意的微笑 一个细微的小动作 就可能决定了
  • CentOS 8部署Gitlab

    Gitlab简单部署 1 环境和准备工作 操作系统 CentOS6 或者7 8测试都是没有问题的 IP地址 192 168 106 129 配置 2核 最少4G 2 建立git用户 useradd git passwd git passwd
  • python摔倒检测,跌倒检测openpose站立行为检测

    python摔倒检测 跌倒检测openpose站立行为检测 import cv2 import numpy as np from torch import from numpy jit from modules keypoints impo
  • 利用fsck修复/dev/mapper/VolGroup-lv_root: unexpected inconsistency; run fsck manually

    fsck file system check 用来检查和维护不一致的文件系统 若系统掉电或磁盘发生问题 可利用fsck命令对文件系统进行检查 如果系统启动的时候提示 dev mapper VolGroup lv root unexpecte
  • 【已解决】mysql报错error: Found option without preceding group in config file: D:\mysql-5.7.25-win32\my.ini

    之前mysql数据库出现问题 因此卸载重安顺便写个教程 顺便给大家一个小建议 如果软件出现什么莫名其妙的错误 半个小时都没有解决 不如直接重装 因为你不知道哪里出了问题 重装比找出问题简单多了 mysql从官网下载安装包直接解压即可 注意最
  • MySQL - 新增数据时想要得到新增数据的自增ID,怎么办???

    我们在实际逻辑编写时 往往会有 想要使用我们刚刚插入数据的自增ID 来做另外一件事情 这时候我们是可以通过查询刚才插入数据的某个字段作为条件 来查出这条数据的ID 但是总感觉太麻烦效率不高 当然 我们通过SQL层面也是有类似的办法 例如 S
  • Spring有哪些注入方式呢?

    转自 Spring有哪些注入方式呢 下文笔者讲述Spring注入方式分享 如下所示 Spring中有以下注入方式 1 setter属性注入 2 构造方法注入 3 lombok注入 例 1 setter属性注入 Controller publ
  • SLAM练习题(十一)—— G2O实战

    SLAM 学习笔记 写在前面的话 算是一点小小的感悟吧 估计位姿的方法有线性方法和非线性方法 线性方法就是特征点法中的2D 2D的对极约束 3D 2D的PnP问题 非线性方法有BA优化 它将位姿的估计问题转换成了一个误差关于优化量的最小二乘