Ceres Solver从零开始手把手教学使用

2023-11-17

目录

一 、简介

二、安装

三、介绍

 四、Hello Word!

五、导数

       1 数值导数

        2解析求导

六、实践 

        Powell函数


一 、简介

        笔者已经半年没有更新新的内容了,最近学习视觉SLAM的过程中发现自己之前学习的库基础不够扎实。Ceres作为一个优化库在视觉SLAM中使用非常广泛,笔者会跟随着Ceres Solver的官方教程一直学习下来,并且会在专栏中持续更新,如需转载请声明出处

二、安装

        虽然说是从零开始教学Cere Solver的使用,但是基本linux和git指令还是需要一定基础。具体的安装步骤如下,需要不同的安装版本可以去GitHub上找到。

        1、依赖项

                Ceres Solver 2.2 需要完全符合 C++17 的编译器。

                Cmake需要3.10版本以上和Eigen3.3以上

                注意官方维护的ubuntu版本必须在18.04以上。老版本的ubuntu可能安装使用会出现不兼容。

# CMake
sudo apt-get install cmake
# google-glog + gflags
sudo apt-get install libgoogle-glog-dev libgflags-dev
# Use ATLAS for BLAS & LAPACK
sudo apt-get install libatlas-base-dev
# Eigen3
sudo apt-get install libeigen3-dev
# SuiteSparse (optional)
sudo apt-get install libsuitesparse-dev

        2、安装

                Ceres是使用Cmake进行安装的,具体的步骤如下:

tar zxf ceres-solver-2.1.0.tar.gz
mkdir ceres-bin
cd ceres-bin
cmake ../ceres-solver-2.1.0
make -j3 
make test
# Optionally install Ceres, it can also be exported using CMake which
# allows Ceres to be used without requiring installation, see the documentation
# for the EXPORT_BUILD_DIR option for more information.
make install

           安装完成之后可以测试自己是否安装成功,可以使用Ceres自带的华盛顿大学的 BAL 数据集进行测试运行看是否安装成功。

bin/simple_bundle_adjuster ../ceres-solver-2.1.0/data/problem-16-22106-pre.txt

     这将使用 DENSE_SCHUR 线性求解器运行 Ceres 最多 10 次迭代。Ceres还可以这Window、Mac、ios和Android安装。具体的安装请参看官方文档。本专栏只针对运行在linux上。

三、介绍

        Ceres可以解决以下边界约束鲁棒化的非线性最小二乘问题:

 这种形式的问题在工程领域运用的非常多比如拟合统计曲线到计算机视觉的图片构建等。

        接下来笔者将带大家一起学习Ceres Solver求解, ρi(‖fi(xi1,...,xik)‖2) 被称为ResidualBlock为残差块, fi(⋅) 为损失函数。 [xi1,...,xik]就是我们需要优化的变量,我们将它定位ParameterBlock。接下来为了理解,举一个例子。在视觉SLAM中,平移有3个分量,旋转使用四元数有4个分量,比如我们最小化重投影误差中需要优化的就是位置和位姿。lj和 uj是参数块的界限xj,当然,lj和 uj都可以取负无穷和正无穷,并且优化的参数也可以是一个,这样可以得到一个我们熟知的非线性最小化问题:

 四、Hello Word!

        学习的开始,我们将和编程一样从一个最简单的Hello Word开始!

        首先,我们先考虑以下函数的最小值:

         从式子可以看出最小值计算x=10的时候,用这样一个简单的例子来讲解Ceres。

        1、首先我们定义一个仿函数来评估函数f(x)=10−x。

                这里我先简单介绍一下什么是仿函数,这里需要一定的C++的基础。仿函数就是在类中对()运算符号进行重载。为什么要使用仿函数主要是为了方便维护类的成员函数和运行效率更高。笔者接触这个也很少,在ORB-SLAM中也同样遇到过仿函数,有空对仿函数专门进行更新。

struct CostFunctor {
   template <typename T>
   bool operator()(const T* const x, T* residual) const {
     residual[0] = 10.0 - x[0];
     return true;
   }
};

  这里解释一下输入的const T* const x, const T*x表示 *x指向的对象不可修改,T* const x表示指针不能修改,const T* const x就是两者的结合指向的对象和指针都不能进行修改。这里我们使用了模板函数,这也是C++的基本知识在这里就不详细介绍了。这里的x变量就是我们需要输入数据,x可以是n维向量。residual则是我们误差(或者残差)也就是损失函数。接下来看主函数:

int main(int, char**) {
    
    // 设置初始值
    double ininal_num=5.0;
    double x=ininal_num;

    // 构建问题Problem有两个重要的成员函数:
    // Problem::AddResidalBlock() and Problem::AddParameterBlock()
    ceres::Problem problem;

    //设置损失函数,这里使用自动求导计算雅克比矩阵以我们自己定义的CostFunctor为类型在进行初始化
    ceres::CostFunction* cost_functor=
    new ceres::AutoDiffCostFunction<CostFunctor,1,1>(new CostFunctor);
    // 添加残差
    problem.AddResidualBlock(cost_functor,nullptr,&x);

    // 配置求解器
    ceres::Solver::Options options;
    // options.max_num_iterations=20; //可以设置最大迭代次数
    options.linear_solver_type=ceres::DENSE_QR;
    options.minimizer_progress_to_stdout=true;
    ceres::Solver::Summary summary;
    //开始求解
    ceres::Solve(options,&problem,&summary);
    std::cout<<summary.BriefReport()<<"\n";
    std::cout<<"x:"<<ininal_num<<"->"<<x<<"\n";
    return 0;

}

五、导数

        和大多数的优化库一样,Ceres能够在任意参数值下估计目标函数中每个项的值和导数。上面我们就使用到了自动求导。现在我们讨论两种解析导数和数值导数。

       1 数值导数

                在某些情况下,是不能写出模板损失函数,比如损失函数使用了其他库中的函数,这样就不能使用自动求导AutoDiffCostFunction。在这种情况下我们可以使用数值导数。数值导数也和自动求导类似,首先定义一个仿函数,之后在传递给NumericDiffCostFunctor。例如:

struct NumericDiffCostFunctor {
  bool operator()(const double* const x, double* residual) const {
    residual[0] = 10.0 - x[0];
    return true;
  }
};
//添加到problem中

CostFunction* cost_function =
  new NumericDiffCostFunction<NumericDiffCostFunctor, ceres::CENTRAL, 1, 1>(
      new NumericDiffCostFunctor);
problem.AddResidualBlock(cost_function, nullptr, &x);

基本的构造和AutoDiffCostFunction相同,就是多了一个参数的输入。这个多出来参数用于计算数值导数的有限差分种类,这个之后讲解。在刚开始学习的过程官方也推荐我们最好使用自动求导

        2解析求导

                解析求导就是在以封闭式计算导数则不能使用自动求导,这是官方所说的有点难懂。简单来说就是你定义了一个抽象函数或者说这个函数中包含了更多函数时,Ceres就不能帮你自动求导了。这种情况就是自己提供残差和雅克比计算了。我们以10-x函数为例子(只是为了方便理解,这个函数可以使用自动求导解决):

class QuadraticCostFunction : public ceres::SizedCostFunction<1, 1> {
 public:
  virtual ~QuadraticCostFunction() {}
  virtual bool Evaluate(double const* const* parameters,
                        double* residuals,
                        double** jacobians) const {
    const double x = parameters[0][0];
    residuals[0] = 10 - x;

    //计算雅克比矩阵
    if (jacobians != nullptr && jacobians[0] != nullptr) {
      jacobians[0][0] = -1;
    }
    return true;
  }
};

我们重构虚函数Evaluate,第一个输入为参数数组,第二个为残差,第三个为雅克比矩阵。Evaluate检查雅克比矩阵是否为空,如果是则用残差函数的导数填充。本例子中因为是线性函数,所以是数值。

官方推荐我们如果我们有能力管理自己的雅克比计算,就使解析导数。否则还是优先自动求导和数值求导。

        Ceres提供了更多的导数计算,还有更加高级的计算DynamicAutoDiffCostFunction和CostFunctionToFunctor在之后我会进一步补充。

六、实践 

        可能你现在还一头雾水,不太理解上面代码的详细功能。接下来我将带大家完成几个例子加深对代码的理解。

        Powell函数

现在我们来考虑一个复一定的问题。powell函数x=[x1,x2,x3,x4]是F(x)中的四个参数有以下等式:

 可以看出我们有四个残差,我们希望1/2‖F(x)‖2是最小的。

首先我们从f4这个式子开始:

 首先还是一样需要定义目标函数中的项进行评估:

struct F1 {
  template <typename T>
  bool operator()(const T* const x1, const T* const x2, T* residual) const {
    // f1 = x1 + 10 * x2;
    residual[0] = x1[0] + 10.0 * x2[0];
    return true;
  }
};
struct F2 {
  template <typename T>
  bool operator()(const T* const x3, const T* const x4, T* residual) const {
    // f2 = sqrt(5) (x3 - x4)
    residual[0] = sqrt(5.0) * (x3[0] - x4[0]);
    return true;
  }
};
struct F3 {
  template <typename T>
  bool operator()(const T* const x2, const T* const x3, T* residual) const {
    // f3 = (x2 - 2 x3)^2
    residual[0] = (x2[0] - 2.0 * x3[0]) * (x2[0] - 2.0 * x3[0]);
    return true;
  }
};
struct F4 {
  template <typename T>
  bool operator()(const T* const x1, const T* const x4, T* residual) const {
    // f4 = sqrt(10) (x1 - x4)^2
    residual[0] = sqrt(10.0) * (x1[0] - x4[0]) * (x1[0] - x4[0]);
    return true;
  }
};

之后的我们开始构造问题:

 double x1 = 3.0;
  double x2 = -1.0;
  double x3 = 0.0;
  double x4 = 1.0;
  Problem problem;
     //添加残差并且使用自动求导
  problem.AddResidualBlock(
      new AutoDiffCostFunction<F1, 1, 1, 1>(new F1), nullptr, &x1, &x2);
  problem.AddResidualBlock(
      new AutoDiffCostFunction<F2, 1, 1, 1>(new F2), nullptr, &x3, &x4);
  problem.AddResidualBlock(
      new AutoDiffCostFunction<F3, 1, 1, 1>(new F3), nullptr, &x2, &x3);
  problem.AddResidualBlock(
      new AutoDiffCostFunction<F4, 1, 1, 1>(new F4), nullptr, &x1, &x4);

之后的优化和之前一样配置就可以了:

Solver::Options options;
options.max_num_iterations = 100;
options.linear_solver_type = ceres::DENSE_QR;
options.minimizer_progress_to_stdout = true;
Solver::Summary summary;
Solve(options, &problem, &summary);
std::cout << summary.FullReport() << "\n";

从我们自身进行观察我们就可以看出最优解是x1=x2=x3=x4=0 ,可以看出最终的优化结果也是是否接近。

 

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

Ceres Solver从零开始手把手教学使用 的相关文章

随机推荐

  • DDR2 DDR3的区别

    DDR2 DDR3的区别 功耗进一步减少 DDR2内存的默认电压为1 8V 而DDR3内存的默认电压只有1 5V 因此内存的功耗更小 发热量也相应地会减少 值得一提的是 DDR3内存还新增了温度监控 采用了ASR Automatic sel
  • Unity5.x运行场景直接卡死的问题

    今天遇到一个很奇葩的问题 就是电脑中Unity5 x版本的都不能运行场景了 包括新建空的工程也不行 重装unity软件重装VS环境也不行 神奇的是2017 2018 2019版本的都没有问题 并且界面也没有任何报错 爬各种论坛谷歌都没有找到
  • 山海演武传·黄道·第一卷 雏龙惊蛰 第二十二 ~ 二十四章 真龙之剑·星墟列将...

    我是第一次 请你 请你温柔一点 少女一边娇喘着 一边将稚嫩的红唇紧贴在男子耳边 樱桃小嘴盈溢出如兰香气 这样子 人家骑在上面 她紧紧地依偎在某个男子身上 窈窕的身躯与丰盈的酥胸 伴随着男子身体晃动而滑上滑下 起伏不定 啊 不要晃的那样厉害
  • ios appStore上架审核通过后,appStore搜索不到该应用

    问题描述 前两天上架一款ios App 周一到公司看审核已经通过了 去appStore上搜索一直搜索不到 ios appStore connect点击提示该商品在中国大陆没有上架 解决方法 通过app store connect 最下面的联
  • vue3.0 + element Plus实现页面中引入弹窗组件及defineExpose用法

    1 在需要弹窗显示的页面引入你所写的弹窗组件 index html
  • AI革命:AI+算力,霸主即将诞生!

    随着人工智能技术的飞速发展 AI 算力 的结合应用已成为科技行业的热点话题 甚至诞生出 AI 算力 最强龙头 的网络热门等式 该组合不仅可以提高计算效率 还可以为各行各业带来更强大的数据处理和分析能力 从而推动创新和增长 那么对于这个时下的
  • 二维多孔介质图像的粒度分布研究(Matlab代码实现)

    欢迎来到本博客 博主优势 博客内容尽量做到思维缜密 逻辑清晰 为了方便读者 座右铭 行百里者 半于九十 本文目录如下 目录 1 概述 2 运行结果 3 参考文献 4 Matlab代码实现 1 概述 使用流域分割算法对岩石二维二值图像进行粒度
  • Scala深入浅出——从Java到Scala

    本文适合有一定Java基础的 并想系统学习Scala的小伙伴借鉴学习 文章有大量实例 建议自己跑一遍 Scala深入浅出 从Java到Scala Scala 一 介绍 1 什么是Scala 2 特点 3 安装 二 Scala特点 三 sca
  • SecureCRT9.1高亮配色设置

    参考 http zh cjh com qita 1623 html https download csdn net download qq 45698138 88310255 spm 1001 2014 3001 5503 1 创建文件co
  • fork的例子

    以下是下列代码的头文件 forks c Examples of Unix process control include
  • Ruoyi-cloud集成Sa-Token SSO单点登录

    文章目录 服务端 客户端前端 客户端后端 https github com dromara Sa Token Sa Token SSO 模式三 修改本地hosts 127 0 0 1 sa sso server com 127 0 0 1
  • ionic3代码压缩和apk优化

    我们在做ionic打包的时候 通常执行这条命令 ionic cordova build android release prod 使用这个命令生成的apk是ionic项目导出的最优化的apk 但是如果还想继续压缩 那么还可以借助Androi
  • Unity 空气墙Shader

    废话不多说 先上效果图 具体代码如下 Shader Hidden AirWall Properties Color Color Color 1 1 1 1 颜色 Interval Interval float 10 间隔 SubShader
  • springmvc注解和参数传递

    一 SpringMVC注解入门 创建web项目 在springmvc的配置文件中指定注解驱动 配置扫描器 Xml代码 收藏代码
  • FFmpeg 实战指南

    文章目录 表达式 滤镜效果 zoompan 中心视距由远及近 中心视距由近及远 水平视距从左到右 水平视距从右到左 垂直视距从上到下 垂直视距从下到上 rotate 顺时针旋转 PI 6 弧度 逆时针旋转 PI 6 弧度 顺时针旋转 45
  • 【Flink】处理函数Process

    目录 处理函数 基本处理函数 ProcessFunction 处理函数的功能 ProcessFunction解析 处理函数的分类 按键分区处理函数 KeyedProcessFunction 定时器Timer 和定时服务 TimerServi
  • 几种css炫酷背景欣赏

    这里为大家带来几种表现惊人的css背景效果 纯css表现效果 有桌布效果 星空效果 心形效果 砖墙效果等 请欣赏 background radial gradient rgba 255 255 255 0 0 rgba 255 255 25
  • 2020-10-29 org.apache.commons.lang3.StringUtils

    public static void TestStr null 和 操作 判断是否Null 或者 System out println StringUtils isEmpty null System out println StringUt
  • 基于神经网络的模式识别

    一 项目设计的目的 通过构建BP网络和离散Hopfield网络模式识别实例 输出稳定结果 二 相关原理知识介绍 BP学习算法是通过反向学习过程使误差最小 其算法过程从输出节点开始 反向地向第一隐含层 即最接近输入层的隐含层 传播由总误差引起
  • Ceres Solver从零开始手把手教学使用

    目录 一 简介 二 安装 三 介绍 四 Hello Word 五 导数 1 数值导数 2解析求导 六 实践 Powell函数 一 简介 笔者已经半年没有更新新的内容了 最近学习视觉SLAM的过程中发现自己之前学习的库基础不够扎实 Ceres