为了防止忘记，特转载至此。本文方法来源是《最小二乘法直线拟合: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; }
下面是运行结果：
从上图可以看出来，红色直线受噪点影响有所偏离，而绿色直线较好地贴合大部分的样本。因此给噪点赋予小的权重可以优化结果。实际应用中，可能采用的策略是第一次拟合时所有点的权重相同；第二次拟合时根据每个点到直线的距离赋予不同的权重来优化结果。