使用Ceres进行曲线拟合

安装Ceres

github

clone 下来之后 用cmake 编译 安装

安装完成后

在/usr/local/include/ceres 下找到 Ceres 的头文件

在/usr/local/lib/下找到名为 libceres.a 的库文件。

在程序中使用Ceres CMakeLists.txt配置

注意，这个配置不好，会报错
直接给出 正确配置

#与Ceres 相关的

set(CMAKE_BUILD_TYPE "Release")
add_compile_options(-std=c++14)
set(CMAKE_CXX_FLAGS_RELEASE "-O3 -Wall -g")

# eigen
find_package(Eigen3 REQUIRED)
include_directories( "/usr/local/include/eigen3/" )

file(GLOB_RECURSE Opencv2.7_LIB "/usr/lib/python2.7/config-x86_64-linux-gnu/*.so")
file(GLOB_RECURSE Creres_LIB "/usr/local/lib/libceres.a")
file(GLOB_RECURSE cholmod_LIB "/usr/lib/x86_64-linux-gnu/libcholmod.so")
file(GLOB_RECURSE lapack_LIB "/usr/lib/lapack/liblapack.so")
file(GLOB_RECURSE f77blas_LIB "/usr/lib/libf77blas.so")
file(GLOB_RECURSE cxsparse_LIB "/usr/lib/x86_64-linux-gnu/libcxsparse.so")
file(GLOB_RECURSE glog_LIB "/usr/lib/x86_64-linux-gnu/libglog.so")
file(GLOB_RECURSE pthread_LIB "/usr/lib/x86_64-linux-gnu/libpthread.so")
file(GLOB_RECURSE Ceres_LIB  ${Creres_LIB} ${cholmod_LIB} ${lapack_LIB} ${f77blas_LIB} ${cxsparse_LIB} ${glog_LIB} ${pthread_LIB})

add_executable(Ceres_Curve_Fitting ceres_curve_fitting.cpp)
target_link_libraries(Ceres_Curve_Fitting ${Ceres_LIB})

进行曲线拟合

假设 曲线的 方程 为

其中 a, b, c 为曲线的参数,w 为高斯噪声

明显这是一个非线性模型

如果有N个 x，y的对应数据，根据该数据求 a、b、c对应的参数，就是曲线拟合的过程。

可以转化为下面的最小二乘问题

公式翻译过来，就是如果选取 a、b、c三个数，使得测量的y和模型计算的y的 差平方最小

在程序中

首先 就要根据 模型生成 x、y的真值。
然后添加 w的噪声 。
最后使用Ceres求解 从带噪声的 数据中 拟合出 三个参数

Code

#include <iostream>
#include <opencv2/core/core.hpp>
#include <ceres/ceres.h>
#include <chrono>


using namespace std;




//代价函数的计算模型
struct CURVE_FITTING_COST
{
    CURVE_FITTING_COST(double x,double y):_x(x),_y(y){}

    //残差的计算
    template <typename T>

    bool operator()(const T* const abc, T* residual )const
    {
        // y-exp(ax^2+bx+c)
        residual[0] = T ( _y ) - ceres::exp ( abc[0]*T ( _x ) *T ( _x ) + abc[1]*T ( _x ) + abc[2] );
        return true;
    }

    const double _x, _y; // x,y 数据
};



int main(int argc,char** argv)
{
    double a=1.0,b=2.0,c=1.0;//真是参数值

    int N=100;  // 数据个数

    double w_sigma =1.0;   //噪声的sigma值

    cv::RNG rng;//  OPENCV 随机数产生器

    double abc[3]={0,0,0};  //a、b、c参数的估计值

    vector<double>  x_data,y_data;//数据

    cout<<"generating data: "<<endl;

    for(int i=0;i<N;i++)//循环产生N个数据
    {
        double x = i/100.0;//x取值 0-1
        x_data.push_back(x);//记录x
        y_data.push_back(exp(a*x*x+b*x+c)+rng.gaussian(w_sigma));//记录添加了噪声产生的y

        cout<<x_data[i]<<" "<<y_data[i]<<endl;
    }

    //构建最小二乘问题
    ceres::Problem problem;
    for(int i=0;i<N;i++)
    {
        /*向问题中添加误差项*/
        problem.AddResidualBlock(
            /*使用自动求导，模板参数：误差类型，输出纬度，输入纬度，数值参照前面struct中写法*/
            new ceres::AutoDiffCostFunction<CURVE_FITTING_COST,1,3>(
                new CURVE_FITTING_COST(x_data[i],y_data[i])
            ),
            nullptr,//  核函数,这里不使用,为空
            abc//待估计参数
        );
    }

    //配置求解器
    ceres::Solver::Options options;//这里有很多配置项可以填
    options.linear_solver_type = ceres::DENSE_QR;//增量方程如何求解
    options.minimizer_progress_to_stdout = true;//输出到cout

    ceres::Solver::Summary summary;//优化信息
    chrono::steady_clock::time_point t1 = chrono::steady_clock::now();//记录时间
    
    /*开始优化*/
    ceres::Solve(options,&problem,&summary);

    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 time cost = "<<time_used.count()<<"seconds"<<endl;

    //输出结果
    cout<<summary.BriefReport()<<endl;
    cout<<"estimated a,b,c = ";
    for ( auto a:abc ) cout<<a<<" ";
    cout<<endl;


    return 0;
}