企业官网建设流程全解析

简介：面向Matlab初、中级使用者以及需要做椭圆曲线拟合的科研、工程人员，这份代码以最小二乘法为框架，借助广义矩阵特征值分解求解椭圆参数，可将离散采样点拟合为一般式椭圆方程，并得到aX^2+bXY+cY^2+dX+eY+f=0的系数结果。压缩包十分轻量，仅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 是平移后的常数项。由此得到：

表里的问号是需要警惕的地方：如果直接用 λ 的绝对值，会掩盖 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 用数据协方差特征值快速验算拟合质量

拟合完成后，可以先不画图，直接对原始点做一次 PCA：

Cov = 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);

自检时记录不同情况下的恢复误差，建议至少覆盖四种数据形态，见下表。

5.2 批处理时的参数一致性约定

批量处理多张图或其他数据源时，把参数写成一致的输出约定能省大量时间：

cx, cy, a, b, phi_deg, a_over_b, ok_flag

其中 phi_deg 统一取 [0, 180) 范围，a 恒大于或等于 b，a_over_b 作为可信度参考。写入表格的时候，a 由特征值绝对值排序决定，角度由特征向量方向归一化决定，两者绑定，避免不同程序段之间出现半轴互换。

5.3 最后一个小技巧：用 Cholesky 分解替代部分特征运算

如果数据点数量很大且只需要中心坐标，不必每次都做完整特征分解。对散布矩阵 S = D'*D 做 Cholesky 分解并求解约束线性方程可以得到相同的中心，但数值上更快更稳。完整的六参数拟合特征提取按 eigen 的方式做即可，日常验证时把svd的结果直接画成长轴短轴覆盖在原点上，肉眼确认角度和椭圆弧段是否贴合，比任何参数指标都直接。

本文还有配套的精品资源，点击获取