基于通用椭圆方程的最小二乘拟合：如何添加F系数？

1. 扩展M矩阵的结构 #

原来针对每个数据点(xi, yi)，M矩阵的行是[xi², xi*yi, yi², xi, yi]，现在要在每行末尾加一个1，对应通用方程里的F项，变成[xi², xi*yi, yi², xi, yi, 1]。

2. 保持目标向量不变 #

方程右侧还是1，所以目标向量依然是全1的列向量，每个数据点对应的等式为：Axi² + Bxi yi + Cyi² + Dxi + Eyi + F = 1。

3. 调整参数求解维度 #

原来求解的是5个参数[A,B,C,D,E]，现在扩展为6个参数[A,B,C,D,E,F]，通过最小二乘求解线性方程组M^T M θ = M^T b即可，其中θ就是这6个参数组成的向量。

4. 代码修改示例（Python） #

假设原代码是这样的：

import numpy as np

def fit_ellipse(x, y):
    M = np.column_stack([x**2, x*y, y**2, x, y])
    b = np.ones_like(x)
    params, _, _, _ = np.linalg.lstsq(M, b, rcond=None)
    A, B, C, D, E = params
    return A, B, C, D, E

修改后支持F系数的版本：

import numpy as np

def fit_ellipse_full(x, y):
    # 新增全1列对应F参数
    M = np.column_stack([x**2, x*y, y**2, x, y, np.ones_like(x)])
    b = np.ones_like(x)
    params, _, _, _ = np.linalg.lstsq(M, b, rcond=None)
    A, B, C, D, E, F = params
    return A, B, C, D, E, F

5. 注意椭圆有效性验证 #

拟合完成后要检查椭圆的必要条件：B² - 4AC < 0，如果不满足，说明拟合结果是双曲线或抛物线，此时可能需要加入约束条件的最小二乘方法来强制得到椭圆。

内容的提问来源于stack exchange，提问作者D.Gordon