简介面向Matlab初、中级使用者以及需要做椭圆曲线拟合的科研、工程人员这份代码以最小二乘法为框架借助广义矩阵特征值分解求解椭圆参数可将离散采样点拟合为一般式椭圆方程并得到aX^2bXYcY^2dXeYf0的系数结果。压缩包十分轻量仅527B共包含2个文件一个可运行的m脚本和一个txt文本说明前者实现椭圆拟合算法后者介绍输入坐标与输出格式。实际使用时将采样点的x、y坐标替换进脚本运行后即可获得椭圆方程各项系数便于后续绘图与误差检验。目前已有641人学习下载可作为理解特征值分解在曲线拟合中应用的基础范例。通过该资源读者能同时掌握最小二乘拟合流程、广义特征值求解思路并可迁移到椭圆检测、图像边缘点拟合等常见任务。1. 椭圆拟合绕不开的特征值问题拿到一堆离散点要拟合椭圆时最常见的错误不是算法选错而是直接把六个系数做普通最小二乘结果得到一条双曲线或者长短轴互换、倾角偏了 90 度。原因在于椭圆约束本质上是一个带不等式的二次型问题而所有靠谱解法最后都会回到同一个动作对 2x2 或 6x6 矩阵做特征分解。命名里带 eigen 的椭圆拟合程序核心也就是这件事。这篇文章把曲线拟合特征抽取的完整链路拆开讲先建立椭圆参数与特征值、特征向量的对应关系再给出可在 MATLAB 里直接跑的拟合程序最后讨论由特征值排序引发的工程坑。适合正在做图像测量、相机标定、点云后处理或者只是想把椭圆拟合代码写得比现有工具更稳的 MATLAB 使用者。2. 二次型与椭圆参数特征分解决定倾角和半轴椭圆在平面上的通用表达式是二次曲线a x^2 b x y c y^2 d x e y f 0它不一定表示椭圆只有当判别式满足 b^2 - 4ac 0 时才是椭圆。把二次项部分写成矩阵形式A [ a b/2 b/2 c ]二次型 x^T A x 的特征值和特征向量直接告诉你椭圆主轴方向和半轴长度。这是整个拟合程序的数学地基。2.1 二次型的特征分解与几何参数对应关系对实对称矩阵 A 做特征分解 A V Λ V^T其中 V 的列是特征向量Λ 是对角矩阵。椭圆方程里A 的主特征向量就是长轴方向对应的特征值决定该方向的“弯曲程度”。更直观的理解是椭圆可以看成高斯分布的等密线协方差矩阵做特征分解得到的主成分方向和椭圆主轴方向一致这就是为什么 PCA 和椭圆拟合经常写在同一段代码里。设中心为 (cx, cy)把坐标平移到中心后二次式变为λ1 u^2 λ2 v^2 F0其中 u、v 是沿着特征向量方向的坐标F0 是平移后的常数项。由此得到几何参数计算方式说明长轴长度 asqrt(abs(F0 / λ_min))特征值绝对值小的方向是长轴短轴长度 bsqrt(abs(F0 / λ_max))特征值绝对值大的方向是短轴倾角 φatan2(V(2,i), V(1,i))i 是选定主轴的列索引中心 (cx, cy)-0.5 * A^(-1) * [d; e]由一次项系数决定表里的问号是需要警惕的地方如果直接用 λ 的绝对值会掩盖 A 中有一个负特征值的情况。椭圆约束下平移后的常数 F0 和两个特征值符号相反所以更稳妥的写法是a sqrt(abs(F0 / lambda(1)))别在意符号先把约束条件满足即可。2.2 用 MATLAB 生成一个已知椭圆验证特征分解要验证拟合程序第一步是生成已知参数的椭圆点。下面这段代码生成中心为 (2, 3)、长轴为 5、短轴为 3、倾角为 30 度的椭圆点集% 生成已知椭圆点集用于后续拟合验证 a_true 5; b_true 3; phi_true pi/6; cx_true 2; cy_true 3; n 300; t linspace(0, 2*pi, n); % 单位圆先拉伸成椭圆再旋转 P [a_true * cos(t), b_true * sin(t)]; R [cos(phi_true), -sin(phi_true); sin(phi_true), cos(phi_true)]; pts P * R [cx_true, cy_true]; plot(pts(:,1), pts(:,2), .); axis equal; grid on;这里用旋转矩阵 R 对单位圆周上的点做线性变换等价于协方差矩阵R * diag([a^2, b^2]) * R。特征值分解这个协方差矩阵得到的特征向量就是主轴的向量特征值开根号就是半轴长度。理解了这条线后面拟合出的参数你才能看明白。2.3 从二次型系数反解几何参数的子函数实际拟合得到的是六个系数不是直接给中心、半轴和角度。所以需要一个转换函数function [cx, cy, a, b, phi] quad_to_ellipse(v) % v [a; b; c; d; e; f]对应一般二次曲线方程 A [v(1), v(2)/2; v(2)/2, v(3)]; lin [v(4); v(5)]; fval v(6); % 中心 -0.5 * A^{-1} * lin cx -0.5 * (A \ lin); cx cx(1); cy cx(2); % 这只是示意实际要分别取两个分量实际代码里应当写成center -0.5 * (A \ lin); cx center(1); cy center(2); % 平移后的常数项 F0 fval - center * A * center; [V, D] eig(A); lambda diag(D); % 特征值按绝对值降序排列保证长轴短轴不颠倒 [~, idx] sort(abs(lambda), descend); a sqrt(abs(F0 / lambda(idx(1)))); b sqrt(abs(F0 / lambda(idx(2)))); phi atan2(V(2, idx(1)), V(1, idx(1))); end注意eig返回的特征向量每列已经归一化但特征向量的方向可能翻转 180 度。若 φ 出现在某个数据里是 30 度另一段数据是 210 度说明不是拟合错误而是特征向量取了反向。后续统一用mod(phi, pi)归一化即可。3. MATLAB 最小二乘椭圆拟合程序从 SVD 到约束特征解拟合程序的输入是一组坐标点 (x_i, y_i)目标是求六个系数使每个点代入方程后的残差最小。这个问题写成矩阵形式后解法的选择直接决定结果是不是椭圆。3.1 齐次最小二乘的解不能直接调 pinv对 n 个点构造设计矩阵D [x.^2, x.*y, y.^2, x, y, ones(n, 1)]; v [a; b; c; d; e; f];目标是让 D * v 接近零向量。注意右边是零不是某个观测向量所以不能写v D \ (-y)这种形式。标准做法是求 D 的最小奇异值对应的右奇异向量也就是约束 ||v|| 1 下的最小二乘解。用 MATLAB 实现非常短[~, ~, V] svd(D, 0); v V(:, end);这步的含义是六个系数被限制在单位球面上避免全零解。SVD 保证即便数据只覆盖一小段弧代进去的方差也是最小的只是结果不一定是椭圆。3.2 判别式校验和 Fitzgibbon 约束解自由 SVD 的解可能是双曲线或抛物线此时判别式 b^2 - 4ac 0。工程上碰到这种情况常见做法是切到约束最小二乘最经典的是 Fitzgibbon 提出的约束条件 4ac - b^2 1。约束矩阵 C 是 6x6 对称矩阵只在前三行三列有值C zeros(6); C(1,3) 2; C(3,1) 2; % 对应 4ac 项 C(2,2) -1; % 对应 -b^2 项然后求解广义特征值问题 S v λ C v其中 S D*D。完整函数可以这样写function [v, ok] fit_quadratic_constrained(x, y) x x(:); y y(:); n numel(x); D [x.^2, x.*y, y.^2, x, y, ones(n,1)]; % 第一步自由最小二乘 [~, ~, V] svd(D, 0); v V(:, end); % 判别式校验b^2 - 4ac 0 才可能是椭圆 if v(2)^2 - 4*v(1)*v(3) 0 ok true; return; end % 第二步Fitzgibbon 约束4ac - b^2 1 S D * D; C zeros(6); C(1,3) 2; C(3,1) 2; C(2,2) -1; [Vg, Lg] eig(S, C); evals real(diag(Lg)); % 取最小正特征值对应特征向量 pos find(evals eps isfinite(evals)); [~, j] min(evals(pos)); v Vg(:, pos(j)); v v / norm(v); % 最后还是校验一次 ok v(2)^2 - 4*v(1)*v(3) 0; end代码里eig(S, C)返回的广义特征值满足 S v C v λ对应广义瑞利商 v^T S v / v^T C v。因为约束是 v^T C v 1取最小正 λ 的向量就是使二次残差最小的椭圆解。选特征值的逻辑如果写成找最大负特征值多半是约束矩阵的符号取反了排查时先检查 C 矩阵里 C(2,2) 的符号。3.3 从六系数到椭圆参数的完整调用流程把前面的子函数串起来一个最小可用的椭圆拟合调用是x pts(:,1); y pts(:,2); [v, ok] fit_quadratic_constrained(x, y); if ~ok error(数据退化无法构成椭圆); end [cx, cy, a, b, phi] quad_to_ellipse(v); fprintf(中心: %.3f, %.3f\n, cx, cy); fprintf(半轴: %.3f, %.3f\n, a, b); fprintf(倾角: %.3f 度\n, rad2deg(mod(phi, pi)));mod(phi, pi)会把角度归一到 0 到 π 之间避免特征向量方向翻转造成角度跳变。若你发现长轴短轴互换多半是quad_to_ellipse里对特征值排序时用了sort的默认升序而特征值有负有正时绝对值排序更重要。前面给的代码已经用sort(abs(lambda), descend)锁定长轴索引。3.4 参数选择的三个工程细节首先是弧段覆盖不足。数据只覆盖椭圆三分之一周长时自由 SVD 的解经常直接退化约束解能强行给出一个椭圆但外推区域误差极大。这种情况不要盲目相信拟合结果至少要看中心是否落在数据点密集区域附近。其次是数值尺度。如果坐标值很大比如像素坐标在数千量级D 矩阵中 x^2 项和常数项差 10^6 倍SVD 或广义特征分解都可能遇到数值警告。常见做法是先做坐标归一化把点平移到以数据重心为原点再除以整体尺度最后把拟合出的中心换算回原坐标系。最后是自由度冗余。六个系数乘以任意非零常数表示同一条曲线所以不需要对系数做额外归一化但比较两组系数时不能直接做差要先归一化到相同尺度或者干脆比较几何参数。4. 几何距离迭代拟合与特征值排序的三个坑代数距离拟合只最小化二次曲线方程的值当点分布不均匀或者噪声与坐标值相关时拟合结果会偏向大尺度区域。更精细的做法是迭代加权最小二乘让每次迭代按点到椭圆的近似几何距离重新分配权重。4.1 近似几何距离的梯度公式点 (x,y) 到曲线 F(x,y) 0 的近似距离可以用一阶泰勒展开r ≈ |F(x,y)| / sqrt(Fx^2 Fy^2)其中 Fx、Fy 是偏导数。对二次曲线Fx 2ax by d Fy bx 2cy e这个公式比直接算点到椭圆的正交距离快一个数量级而且实现简单适合在拟合循环里反复调用。4.2 用 IRWLS 实现抗离群点拟合在上一章函数基础上加上权重矩阵的迭代就构成一个鲁棒拟合function [v, w] fit_ellipse_irls(x, y, maxiter) x x(:); y y(:); n numel(x); D [x.^2, x.*y, y.^2, x, y, ones(n,1)]; w ones(n,1); for iter 1:maxiter Dw D .* sqrt(w); [~, ~, V] svd(Dw, 0); v V(:, end); % 近似几何距离 Fx 2*v(1)*x v(2)*y v(4); Fy v(2)*x 2*v(3)*y v(5); denom sqrt(Fx.^2 Fy.^2) eps; r abs(D*v) ./ denom; % 权重用中位数绝对偏差做尺度估计 mad_s 1.4826 * median(abs(r - median(r))); w 1 ./ (r 0.3 * mad_s); if iter 1 norm(v - v_prev, inf) 1e-10 break; end v_prev v; end end权重公式里的 0.3 是收缩系数防止权值过大造成震荡。MAD中位数绝对偏差系数 1.4826 使它在高斯噪声下等价于标准差对离群点不敏感。迭代 10 轮通常足够收敛如果数据里离群点超过 30%建议先目视检查后再决定是否换数据源而不是无限制加大迭代次数。4.3 特征值排序、角度周期与半轴互换用eig提取椭圆参数时三个问题几乎每台机器上都会遇到。第一是特征值不是天然有序的。eig(A)返回的 D 矩阵对角元素可能按升序也可能按降序取决于底层 LAPACK 实现。必须显式排序并且按绝对值排否则 a、b 会互换。第二是特征向量方向翻转。同一个主轴特征向量 (cosφ, sinφ) 和 (-cosφ, -sinφ) 都合法导致角度相差 π。归一化角度前统一用mod(phi, pi)处理让输出落在 [0, π) 区间。第三是近圆退化。当长短轴之比接近 1 时角度对噪声极其敏感。半轴上 1% 的误差可能让倾角偏 20 度。这时应在输出里附带形状比 a/b由调用方决定是否该信任倾角。4.4 用数据协方差特征值快速验算拟合质量拟合完成后可以先不画图直接对原始点做一次 PCACov cov(x, y); [~, Latent] eig(Cov); [evals, ~] sort(diag(Latent), descend); shape_ratio sqrt(evals(2) / evals(1));这个 shape_ratio 是点云分布的长短轴之比虽然不能直接当椭圆半轴比但对判断退化很有效。如果数据点沿圆弧分布这个值会接近 1说明点云在方向上没有强烈偏好这时拟合出的椭圆角度基本没有意义如果拟合程序返回的 a/b 与此量级差异巨大通常是参数提取或约束分支哪里出了问题。5. 合成数据自检与批处理时保留的几个验证技巧椭圆拟合程序的验收不能只看一张图画得圆不圆。建议准备一个自检脚本先生成已知参数的椭圆再叠加噪声做污染最后比较恢复的参数与真实值的误差。5.1 闭合验证的脚本骨架% 生成已知椭圆并加噪声 t linspace(0, 2*pi, 200); P 5*[cos(t), 2*sin(t)] * R(pi/4) [1, 2]; P P 0.1 * randn(200, 2); % 高斯噪声 [v, ok] fit_quadratic_constrained(P(:,1), P(:,2)); [cx, cy, a, b, phi] quad_to_ellipse(v);自检时记录不同情况下的恢复误差建议至少覆盖四种数据形态见下表。数据形态验证重点推荐配置完整椭圆 360 度基础精度自由 SVD 即可iter 可省只覆盖 120 度弧段约束分支走 Fitzgibbon 分支检查中心漂移加入 20% 离群点鲁棒性IRWLS 权重迭代 10 轮长短轴比接近 1角度稳定性输出附上 a/b提示角度不可靠5.2 批处理时的参数一致性约定批量处理多张图或其他数据源时把参数写成一致的输出约定能省大量时间cx, cy, a, b, phi_deg, a_over_b, ok_flag其中 phi_deg 统一取 [0, 180) 范围a 恒大于或等于 ba_over_b 作为可信度参考。写入表格的时候a 由特征值绝对值排序决定角度由特征向量方向归一化决定两者绑定避免不同程序段之间出现半轴互换。5.3 最后一个小技巧用 Cholesky 分解替代部分特征运算如果数据点数量很大且只需要中心坐标不必每次都做完整特征分解。对散布矩阵 S D*D 做 Cholesky 分解并求解约束线性方程可以得到相同的中心但数值上更快更稳。完整的六参数拟合特征提取按 eigen 的方式做即可日常验证时把svd的结果直接画成长轴短轴覆盖在原点上肉眼确认角度和椭圆弧段是否贴合比任何参数指标都直接。本文还有配套的精品资源点击获取
