Ceres::AutoDiffCostFunction拟合曲线

Ceres::AutoDiffCostFunction拟合曲线

基本思路：通过生成带有高斯噪声的数据，学习使用ceres拟合曲线。
结论：AutoDiffCostFunction只能拟合一些低阶线性和非线性的曲线。

参考：
http://www.ceres-solver.org/nnls_tutorial.html#curve-fitting
//https://stackoverflow.com/questions/32889309/adding-gaussian-noise

生成数据

在ground truth上加上噪声作为原始数据，曲线方程为$y=e^{3x}+2sin(x)$

#include<random>
\\...
    const double mean = 0.0;
    const double stddev = 0.5;
    std::default_random_engine generator;
    std::normal_distribution<double> dist(mean, stddev);
\\...
	double y_=exp(3.0*x_)+2.0*sin(x_)+dist(generator);

构建ceres求解器

构建以ceres::AutoDiffCostFunction为求解方式的ceres求解器主要包含两个步骤：构建ceres::Problem对象、构建目标成本函数ceres::CostFunction

1）构建ceres::CostFunction

第一步是定义结构体，在结构体中重载括号运算符“()”，在重载函数中定义残差计算方式，实际上就是变量之间应该是什么样的关系，比如说要拟合的曲线的参数是$y=e^{ax}+bsin(x)$中的$a$和$b$的值，那么结构体中的重载函数就直接写成$residual=y-[e^{ax}+bsin(x)]$：

struct curve_fitting
{
    curve_fitting(const double x,const double y):x_(x),y_(y){

    }
    template<typename T>
    bool operator()(const T* const a,const T* const b,T* residual)const{
        residual[0]=y_-exp(a[0]*x_)-b[0]*sin(x_);
        return true;
    }
private:
    const double x_;
    const double y_;
};

在这个地方我尝试了一些高阶的方程，发现都没有办法求解，可能是ceres::AutoDiffCostFunction只能求解一些低阶的方程。

2）构建ceres::Problem

有了函数对象，就可以结合函数对象和数据构建ceres的Problem类。

Problem problem;

    for(int i=0;i<x_num;i++){
        CostFunction * costf = new AutoDiffCostFunction<curve_fitting,1,1,1>(
                                    new curve_fitting(data_x.at(i),data_y.at(i)));
        problem.AddResidualBlock(costf,nullptr,&a,&b);
    }

这里就解释了为什么在 1） 里面写结构体的时候需要定义两个私有成员变量x_ 和 y_，因为一共需要往Problem里面放x_num（x_num是数据的个数）个残差（residual），每一组残差由一组数据得到，那么x_ 和 y_ 就分别存放某一个残差对应的数据。
上述代码可以理解为ceres::AutoDiffCostFunction是ceres::CostFunction的一个派生类

CostFunction * costf = new AutoDiffCostFunction<curve_fitting,1,1,1>(
                                    new curve_fitting(data_x.at(i),data_y.at(i)));

这部分代码可以理解为用 1) 中定义好的结构体作为模板，实例化一个AutoDiffCostFunction类，用AutoDiffCostFunction类实例化一个CostFunction类。

CostFunction * costf = new AutoDiffCostFunction<curve_fitting,1,1,1>(
									^__实例化一个AutoDiffCostFunction类，以curve_fitting为模板，第一个1代表残差项个数为1，第二个1代表待优化变量1的个数为1，第三个1代表待优化变量2的个数为1
                                    new curve_fitting(data_x.at(i),data_y.at(i)));
                                    ^__实例化一个curve_fitting类，构造函数列表参数为data_x.at(i)和data_y.at(i)

向ceres::Problem中放入构建好的CostFunction

problem.AddResidualBlock(costf,nullptr,&a,&b);

求解

求解过程主要的步骤是调用ceres::Solve方法，往其中传入求解对象（ceres::Problem）和求解过程中的一些配置，比如要用什么方法求解，求解过程中允许的最大迭代次数是多少，求解过程中是不是详细输出每一步的求解结果等等。

    Solver::Options option;
    option.minimizer_progress_to_stdout=true;
    option.max_num_iterations=100;
    option.linear_solver_type;ceres::DENSE_QR;
    
    Solver::Summary summary;

    Solve(option,&problem,&summary);

结果

代码

#include<ceres/ceres.h>
#include<iostream>
#include <random>
#include<fstream>
using ceres::AutoDiffCostFunction;
using ceres::CostFunction;
using ceres::Problem;
using ceres::Solve;
using ceres::Solver;
static int x_num=0;
std::string prefix_path{"/home/user/C++Project/ceres_use/data/"};
void generate_seq(std::vector<double>& x,std::vector<double>& y){
    const double mean = 0.0;
    const double stddev = 0.5;
    std::default_random_engine generator;
    std::normal_distribution<double> dist(mean, stddev);
    //https://stackoverflow.com/questions/32889309/adding-gaussian-noise
    std::ofstream file;
    file.open(prefix_path+"data.txt");
    if(file.is_open()){
        for(double x_=-1.;x_<1.;x_=x_+0.01){
            x_num++;
            //double y_=(exp(3.0*pow(x_,3))-1.)/(x_-sin(x_))+dist(generator);
            double y_=exp(3.0*x_)+2.0*sin(x_)+dist(generator);
            x.push_back(x_);
            y.push_back(y_);
            file<<std::to_string(x_);
            file<<" ";
            file<<std::to_string(y_);
            file<<"\n";
        }
    }
    file.close();
}
struct curve_fitting
{
    curve_fitting(const double x,const double y):x_(x),y_(y){

    }
    template<typename T>
    bool operator()(const T* const a,const T* const b,T* residual)const{
        //residual[0]=y_*(b[0]*x_-sin(x_))-(exp(a[0]*pow(x_,3))-1.);
        //residual[0]=y_-(exp(a[0]*pow(x_,3))-1.0)/(b[0]*x_-sin(x_));
        residual[0]=y_-exp(a[0]*x_)-b[0]*sin(x_);
        return true;
    }
private:
    const double x_;
    const double y_;
};

int main(int argc,char** argv){
    std::vector<double> data_x;
    std::vector<double> data_y;
    double a{0.0};
    double b{0.0};
    generate_seq(data_x,data_y);
    std::cout<<"x_num => "<<x_num<<std::endl;

    Problem problem;

    for(int i=0;i<x_num;i++){
        CostFunction * costf = new AutoDiffCostFunction<curve_fitting,1,1,1>(
                                    new curve_fitting(data_x.at(i),data_y.at(i)));
        problem.AddResidualBlock(costf,nullptr,&a,&b);
    }

    Solver::Options option;
    option.minimizer_progress_to_stdout=true;
    option.max_num_iterations=100;
    option.linear_solver_type;ceres::DENSE_QR;
    
    Solver::Summary summary;

    Solve(option,&problem,&summary);

    std::cout<<summary.BriefReport()<<std::endl;
    std::cout<<"final a => "<<a<<std::endl<<"final b => "<<b<<std::endl;
    return 0;
}

CMakeLists.txt

cmake_minimum_required(VERSION 3.5)

project(main)

find_package(Ceres REQUIRED)


add_executable(curfitting curve_fitting.cpp)
target_link_libraries(curfitting Ceres::ceres)