如何使用Ceres Solver库拟合二阶多元多项式（3D曲面）

问题说明 #

你需要拟合的多项式形式为：

z = a + b*x + c*y + d*x*y + e*x^2*y + f*y^2*x + g*x^2*y^2

共7个待求系数a~g，以下是不需要复杂模板编程的实现步骤，直接套用即可。 #

实现步骤 #

1. 定义残差计算结构体 #

用Ceres提供的自动微分能力，不需要手动推导雅可比，仅需要定义残差计算逻辑：

#include <ceres/ceres.h>
#include <iostream>

// 残差计算 functor
struct ResidualFunctor {
    ResidualFunctor(double x, double y, double z) : x_(x), y_(y), z_(z) {}

    // 模板化的()运算符，供自动微分调用，不需要自行修改
    template <typename T>
    bool operator()(const T* const a, const T* const b, const T* const c,
                    const T* const d, const T* const e, const T* const f,
                    const T* const g, T* residual) const {
        // 计算残差：观测值z减去多项式预测值
        residual[0] = T(z_) - (a[0] + b[0]*T(x_) + c[0]*T(y_) + d[0]*T(x_)*T(y_)
                            + e[0]*T(x_)*T(x_)*T(y_) + f[0]*T(y_)*T(y_)*T(x_)
                            + g[0]*T(x_)*T(x_)*T(y_)*T(y_));
        return true;
    }

private:
    const double x_; // 观测的x坐标（pupil-glints纵向坐标）
    const double y_; // 观测的y坐标（pupil-glints横向坐标）
    const double z_; // 观测的z坐标（给定点纵向坐标）
};

2. 构建优化问题并求解 #

你只需要把你实验采集到的(x,y,z)观测数据替换到下面的示例中即可：

int main(int argc, char** argv) {
    // 待求解参数，顺序对应a,b,c,d,e,f,g，初始值可以设为0或者你有先验的初始估计
    double params[7] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};

    ceres::Problem problem;
    // ---------------------- 替换为你的观测数据 start ----------------------
    // 示例：模拟5组观测数据，你需要换成自己实验采集的所有数据
    double observations[5][3] = {
        {1.0, 2.0, 5.1},
        {2.0, 1.0, 4.8},
        {3.0, 2.0, 12.3},
        {1.0, 3.0, 8.9},
        {2.0, 2.0, 7.2}
    };
    int obs_num = sizeof(observations)/sizeof(observations[0]);
    for (int i = 0; i < obs_num; ++i) {
        double x = observations[i][0];
        double y = observations[i][1];
        double z = observations[i][2];
        // 给每一组观测添加残差块
        problem.AddResidualBlock(
            // 自动微分模板参数：<残差类型, 残差维度, 7个参数各自的维度>，不需要修改
            new ceres::AutoDiffCostFunction<ResidualFunctor, 1, 1,1,1,1,1,1,1>(
                new ResidualFunctor(x, y, z)
            ),
            nullptr, // 如果有外点可以换成鲁棒核函数：new ceres::HuberLoss(1.0)
            &params[0], &params[1], &params[2],
            &params[3], &params[4], &params[5], &params[6]
        );
    }
    // ---------------------- 替换为你的观测数据 end ----------------------

    // 配置求解器
    ceres::Solver::Options options;
    options.max_num_iterations = 100; // 最大迭代次数，该场景下足够使用
    options.linear_solver_type = ceres::DENSE_QR; // 小维度问题用DENSE_QR性能最优
    options.minimizer_progress_to_stdout = true; // 不需要打印迭代过程可以设为false

    ceres::Solver::Summary summary;
    ceres::Solve(options, &problem, &summary);

    // 输出结果
    std::cout << "迭代结果摘要：" << summary.BriefReport() << std::endl;
    std::cout << "拟合得到的系数：" << std::endl;
    std::cout << "a = " << params[0] << std::endl;
    std::cout << "b = " << params[1] << std::endl;
    std::cout << "c = " << params[2] << std::endl;
    std::cout << "d = " << params[3] << std::endl;
    std::cout << "e = " << params[4] << std::endl;
    std::cout << "f = " << params[5] << std::endl;
    std::cout << "g = " << params[6] << std::endl;

    return 0;
}

3. 编译配置 #

你的CMakeLists.txt需要添加Ceres和Eigen的依赖，参考配置：

cmake_minimum_required(VERSION 3.10)
project(poly_fit)

find_package(Ceres REQUIRED)
include_directories(${CERES_INCLUDE_DIRS})

add_executable(poly_fit main.cpp)
target_link_libraries(poly_fit ${CERES_LIBRARIES})

注意事项 #

• 如果你采集的观测数据存在明显异常点，可以把AddResidualBlock的第二个参数从nullptr换成new ceres::HuberLoss(1.0)，使用鲁棒核函数降低外点的干扰，阈值可以根据你的z值范围调整。
• 你当前的多项式本质是线性最小二乘问题，也可以直接用Eigen的矩阵分解求解，但上面的Ceres实现结构通用性更强，后续如果需要修改为更复杂的非线性模型不需要大改代码。
• 整个实现不需要你手动编写模板逻辑，直接套用上面的结构替换观测数据即可运行。

内容的提问来源于stack exchange，提问作者PHỦ Đỗ Tường