1. 包含必要的头文件

#include <opencv2/core.hpp>
#include <opencv2/imgproc.hpp>
#include <opencv2/highgui.hpp>
#include <vector>
#include <cmath>

2. 定义生成椭圆点的函数

编写一个函数，接受椭圆的中心坐标、长轴半径、短轴半径、旋转角度和点的数量作为参数，生成椭圆上的点并存储到 std::vector<cv::Point2d> 中：

void generateEllipsePoints(double xc, double yc, double a, double b, double theta, int num_points, std::vector<cv::Point2d>& points) {points.clear(); // 清空之前的点double angle_step = 2 * M_PI / num_points; // 每个点的角度间隔for (int i = 0; i < num_points; ++i) {double angle = i * angle_step;// 计算椭圆上点的坐标double x = xc + a * cos(angle) * cos(theta) - b * sin(angle) * sin(theta);double y = yc + a * cos(angle) * sin(theta) + b * sin(angle) * cos(theta);points.emplace_back(x, y);}
}

3. 定义误差函数

编写一个误差函数，用于计算每个点到拟合椭圆的距离：

• 计算当前参数下，所有点到椭圆的误差。
• p 是当前优化的参数（5 个值）。
• x 是计算得到的误差（x_r^2 / a^2 + y_r^2 / b^2 - 1）。
• data 是传入的离散点数据。

void ellipse_residuals(double *p, double *x, int m, int n, void *data) {double xc = p[0];  // 椭圆中心 xdouble yc = p[1];  // 椭圆中心 ydouble a = p[2];   // 长轴半径double b = p[3];   // 短轴半径double theta = p[4];  // 旋转角度double cos_t = cos(theta);double sin_t = sin(theta);double *points = (double *)data;for (int i = 0; i < n; ++i) {double x_i = points[2 * i];     // x 坐标double y_i = points[2 * i + 1]; // y 坐标// 计算旋转后的坐标double x_r = (x_i - xc) * cos_t + (y_i - yc) * sin_t;double y_r = -(x_i - xc) * sin_t + (y_i - yc) * cos_t;// 计算残差x[i] = (x_r * x_r) / (a * a) + (y_r * y_r) / (b * b) - 1;}
}

4. 使用 levmar 进行椭圆拟合

• 定义 points[] 存储数据点，每两个值是一对 (x, y)。
• p[] 存储初始参数，包括椭圆中心、长短轴、旋转角度。
• opts[] 是 levmar 的优化控制参数。
• dlevmar_dif 进行优化，nullptr 代表不需要目标函数值。
• 代码最终输出 xc, yc, a, b, theta，即拟合出的椭圆参数。

int main() {double xc = 0.0; // 椭圆中心 x 坐标double yc = 0.0; // 椭圆中心 y 坐标double a = 5.0;   // 长轴半径double b = 3.0;   // 短轴半径double theta = M_PI / 4; // 旋转角度 45 度int num_points = 100;     // 生成 100 个点std::vector<cv::Point2d> ellipse_points;generateEllipsePoints(xc, yc, a, b, theta, num_points, ellipse_points);// 将 std::vector<cv::Point2d> 转换为 double 数组double* points = new double[2 * num_points];for (int i = 0; i < num_points; ++i) {points[2 * i] = ellipse_points[i].x;points[2 * i + 1] = ellipse_points[i].y;}// 初始猜测参数 [xc, yc, a, b, theta]double p[5] = {0.0, 0.0, 5.0, 3.0, M_PI / 4};// 设置优化参数double opts[LM_OPTS_SZ], info[LM_INFO_SZ];opts[0] = LM_INIT_MU;  // 初始阻尼因子opts[1] = 1e-10;       // 终止条件：误差变化opts[2] = 1e-10;       // 终止条件：参数变化opts[3] = 1e-10;       // 终止条件：梯度opts[4] = LM_DIFF_DELTA; // 数值微分步长// 运行优化int m = 5;  // 参数数目int n = num_points;  // 观测值数目（残差数）dlevmar_dif(ellipse_residuals, p, nullptr, m, n, 1000, opts, info, nullptr, nullptr, points);// 输出拟合结果std::cout << "拟合椭圆参数:" << std::endl;std::cout << "xc = " << p[0] << ", yc = " << p[1] << std::endl;std::cout << "a = " << p[2] << ", b = " << p[3] << std::endl;std::cout << "theta = " << p[4] * 180 / M_PI << "°" << std::endl;// 清理内存delete[] points;return 0;
}