拟合圆的梯度下降法例子

最近研究了一下梯度下降法，所以写了个拟合圆的方法。拟合圆属于非线性拟合。网上的最小二乘法拟合圆公式并不是误差的平方，而是4次方（为了去掉公式里的开方）。一般可以先用网上的公式得到一个初始解，然后再用梯度下降法继续求精。梯度下降法的公式推导如下。定义误差函数：

$${e=\sum_{i}^{}\left [ \sqrt{ \left (x _{i}-x _{0} \right ) ^{2} + \left (y _{i}-y _{0} \right ) ^{2} } - r \right ] ^{2}}$$

分别对${x_{0},y_{0},r}$求偏导数可得：

$${ \begin{align} \frac{\partial e}{\partial x _{0}}&=\sum_{i}^{} 2f_{i} \left(x_{0},y_{0},r \right)\cdot \frac{1}{2}\cdot \frac{1}{g_{i} \left(x_{0},y_{0},r \right)}\cdot -2 \left(x_{i}-x_{0} \right) \\ \frac{\partial e}{\partial y _{0}}&=\sum_{i}^{} 2f_{i} \left(x_{0},y_{0},r \right)\cdot \frac{1}{2}\cdot \frac{1}{g_{i} \left(x_{0},y_{0},r \right)}\cdot -2 \left(y_{i}-y_{0} \right) \\ \frac{\partial e}{\partial r}&=\sum_{i}^{} 2f_{i} \left(x_{0},y_{0},r \right)\cdot -1 \end{align}}$$

其中：

$${ \begin{align*} g_{i} \left( x_{0},y_{0},r \right ) &= \sqrt{ \left (x _{i}-x _{0} \right ) ^{2} + \left (y _{i}-y _{0} \right ) ^{2} } \\ f_{i} \left(x_{0},y_{0},r \right ) &= g_{i} \left( x_{0},y_{0},r \right )-r \end{align*}}$$

最终迭代公式为：

$${ \left(x_{0},y_{0},r \right)_{n+1} = \left(x_{0},y_{0},r \right)_{n}-\left ( \frac{\partial e}{\partial x _{0}},\frac{\partial e}{\partial y _{0}},\frac{\partial e}{\partial r} \right ) \cdot C}$$

下述代码基于VS2017、Qt5.9和OpenCV430，通过了验证。代码中为了加速收敛限制了迭代的最小步长为0.005。如果步长小于它则会放大步长，因此该程序的求解精度约为5‰。代码如下：

void main()
{
    vector<Point2f> points = { { -1, -1 }, { 1, 1 }, { 1, -1 } };
    Matx13f resolve = Matx13f::zeros();
    for (int loop = 0; loop < 50; loop++)
    {
        float dex = 0;
        float dey = 0;
        float der = 0;
        for (auto pt : points)
        {
            float dx = pt.x - resolve(0);
            float dy = pt.y - resolve(1);
            float g = sqrtf(dx * dx + dy * dy);
            float f = g - resolve(2);
            dex += -2 * f / g * dx;
            dey += -2 * f / g * dy;
            der += -2 * f;
        }
        dex /= points.size();
        dey /= points.size();
        der /= points.size();
        dex *= 0.1f;
        dey *= 0.1f;
        der *= 0.1f;
        float maxv = std::max({ fabs(dex), fabs(dey), fabs(der) });
        if (maxv < 0.005f)
        {
            dex = dex / maxv * 0.005f;
            dey = dey / maxv * 0.005f;
            der = der / maxv * 0.005f;
        }
        resolve(0) -= dex;
        resolve(1) -= dey;
        resolve(2) -= der;
        qDebug() << dex << dey << der;
    }
    qDebug() << resolve(0) << resolve(1) << resolve(2) << "-<";
}

 附，在博文中插入公式的方法：

1. 在博客设置里添加支持数学公式（首页->管理->选项）；
2. 在Online LaTeX Equation Editor - create, integrate and download (codecogs.com)编辑数学公式，需要学习语法；
3. 把数学公式复制到博文中，并且用${<你的公式>}$或$${<你的公式>}$$包裹起来。注意公式内部不能换行；
4. 发表博文就能看到效果了。