直线一般式拟合直线

为了防止忘记，特转载至此。本文方法来源是《最小二乘法直线拟合:Ax+By+C=0 - 会飞的大象会飞的大象 (whudj.cn)》。用一次函数${ y=kx+b }$形式拟合直线非常简单，直接带入最小二乘法公式就行了。而用直线一般式${ ax+by+c=0 }$拟合由于不是线性方程组则需要一些求解技巧。这里不再重复原文内容，而是在原文基础上更进一步。考虑实际情况，在机器视觉应用中从图像中提取的点往往包含一定的噪声，所以通常采取多次加权拟合以消除噪点对结果的影响。下面将给出加权拟合的公式。拟合直线的加权的误差函数如下：

$${ \begin{equation} e=\sum_{i=1}^{N}w_{i}(ax_{i}+by_{i}+c)^{2},满足a^{2}+b^{2}=1 \end{equation} }$$

对c求偏导数并令其为0：

$${ \frac{\partial e}{\partial c}=2\sum_{i=1}^{N}w_{i}(ax_{i}+by_{i}+c)=0 }$$

解得：

$${ \begin{equation} c=-\left( a\frac{ \sum_{i=1}^{N}w_{i}x_{i} }{ \sum_{i=1}^{N}w_{i} }+b\frac{ \sum_{i=1}^{N}w_{i}y_{i} }{ \sum_{i=1}^{N}w_{i} } \right) \end{equation} }$$

把上式代入(1)式，有：

$${ \begin{align*} e &= \sum_{j=1}^{N} w_{j}\left( ax_{j}+by_{j} - a\frac{ \sum_{i=1}^{N}w_{i}x_{i} }{ \sum_{i=1}^{N}w_{i} }-b\frac{ \sum_{i=1}^{N}w_{i}y_{i} }{ \sum_{i=1}^{N}w_{i} } \right)^{2} \\  &= \sum_{j=1}^{N} \left[ a \sqrt{w_{j}} \left( x_{j}-\frac{ \sum_{i=1}^{N}w_{i}x_{i} }{ \sum_{i=1}^{N}w_{i} } \right) + b \sqrt{w_{j}} \left( y_{j}-\frac{ \sum_{i=1}^{N}w_{i}y_{i} }{ \sum_{i=1}^{N}w_{i} } \right) \right]^{2} \end{align*} }$$

现在令${ p_{j}=\sqrt{w_{j}} \left( x_{j}-\frac{ \sum_{i=1}^{N}w_{i}x_{i} }{ \sum_{i=1}^{N}w_{i} } \right),q_{j}=\sqrt{w_{j}} \left( y_{j}-\frac{ \sum_{i=1}^{N}w_{i}y_{i} }{ \sum_{i=1}^{N}w_{i} } \right) }$则上式可简写成：

$${ \begin{align*} e &=\sum_{j=1}^{N} \left( a p_{j} + b q_{j} \right)^{2} \\  &=\begin{pmatrix} a & b \end{pmatrix}\underbrace{\begin{pmatrix} \sum_{j=1}^{N} p_{j}^{2} & \sum_{j=1}^{N} p_{j}q_{j} \\ \sum_{j=1}^{N} p_{j}q_{j} & \sum_{j=1}^{N} q_{j}^{2} \end{pmatrix}}_{S}\begin{pmatrix} a \\ b \end{pmatrix} \end{align*} }$$

接下来的内容就跟不加权的方法相同了。矩阵S的最小特征值对应的特征向量里的两个数就是系数a,b的值。至于为什么可以搜索“二次型求最小值”相关内容。调用OpenCV的cv::eigen(...)函数可得到矩阵的特征值与特征向量。c的值用式(2)求。对应的代码如下。程序运行环境是VS2017和OpenCV430：

Point3f fitLine(const vector<Point2f>& points, const vector<float>& weights)
{
    float swx = 0;
    float swy = 0;
    float sw = 0;
    int count = (int)points.size();
    for (int i = 0; i < count; i++)
    {
        swx += weights[i] * points[i].x;
        swy += weights[i] * points[i].y;
        sw += weights[i];
    }
    float p = 0;
    float q = 0;
    Matx22f pq = Matx22f::zeros();
    for (int i = 0; i < count; i++)
    {
        float pi = sqrtf(weights[i]) * (points[i].x - (swx / sw));
        float qi = sqrtf(weights[i]) * (points[i].y - (swy / sw));
        pq(0, 0) += pi * pi;
        pq(0, 1) += pi * qi;
        pq(1, 0) += pi * qi;
        pq(1, 1) += qi * qi;
    }
    Mat eval, evec;
    eigen(pq, eval, evec);
    float a, b, c;
    a = evec.at<float>(1, 0);
    b = evec.at<float>(1, 1);
    c = -(a * swx + b * swy) / sw;
    return Point3f(a, b, c);
}

int main()
{
    vector<Point2f> points;
    points.push_back(Point2f(10, -0.3f));
    points.push_back(Point2f(20.3f, 10));
    points.push_back(Point2f(30, 18.0f));
    points.push_back(Point2f(80, 71.0f));
    points.push_back(Point2f(110, 100));
    points.push_back(Point2f(215.0f, 200));
    points.push_back(Point2f(315.0f, 299));
    points.push_back(Point2f(425.0f, 412));
    points.push_back(Point2f(565.4f, 560));
    points.push_back(Point2f(410, 440));
    points.push_back(Point2f(59, 20));

    vector<float> weights = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 };

    Point3f coes = fitLine(points, weights);
    Mat image(600, 600, CV_8UC3, Scalar(0, 0, 0));
    for (auto& item : points)
    {
        image.at<Vec3b>(item) = Vec3b::all(255);
    }
    Point2f p1(0,0), p2(600,0);
    p1.y = -(coes.x * p1.x + coes.z) / coes.y;
    p2.y = -(coes.x * p2.x + coes.z) / coes.y;
    line(image, p1, p2, Scalar(127, 127, 255)); /* 红色 */

    weights[9] = 0.2f;
    weights[10] = 0.2f;
    coes = fitLine(points, weights);
    p1.y = -(coes.x * p1.x + coes.z) / coes.y;
    p2.y = -(coes.x * p2.x + coes.z) / coes.y;
    line(image, p1, p2, Scalar(127, 255, 127)); /* 绿色 */
    imshow("拟合直线示意图", image);
    waitKey(-1);
    return 0;
}

下面是运行结果：

从上图可以看出来，红色直线受噪点影响有所偏离，而绿色直线较好地贴合大部分的样本。因此给噪点赋予小的权重可以优化结果。实际应用中，可能采用的策略是第一次拟合时所有点的权重相同；第二次拟合时根据每个点到直线的距离赋予不同的权重来优化结果。