直线一般式拟合直线
为了防止忘记,特转载至此。本文方法来源是《最小二乘法直线拟合: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{equation} \begin{split} 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}}_{ \mathbf{S} }\begin{pmatrix} a \\ b \end{pmatrix} \end{split} \end{equation} }$$
接下来的内容就跟不加权的方法相同了。矩阵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; }
下面是运行结果:
从上图可以看出来,红色直线受噪点影响有所偏离,而绿色直线较好地贴合大部分的样本。因此给噪点赋予小的权重可以优化结果。实际应用中,可能采用的策略是第一次拟合时所有点的权重相同;第二次拟合时根据每个点到直线的距离赋予不同的权重来优化结果。
附,关于(3)式中矩阵S的最小特征值对应的单位特征向量是${ a,b }$的解的说明
为了求e的最小值,我们对(3)式中的e构造拉格朗日乘子:
$${ e=\lambda \left( a^{2}+b^{2}-1 \right) + \sum_{j=1}^{N} \left( a p_{j} + b q_{j} \right)^{2} }$$
对各个变量求偏导数:
$${ \left\{ \begin{matrix} \frac{\partial e}{\partial a}=2 \lambda a + 2 \sum_{j=1}^{N}p_{j}\left( ap_{j}+bq_{j} \right)=0 \\ \frac{\partial e}{\partial b}=2 \lambda b + 2 \sum_{j=1}^{N}q_{j}\left( ap_{j}+bq_{j} \right)=0 \\ \frac{\partial e}{\partial \lambda}=a^{2}+b^{2}-1 = 0 \end{matrix} \right. }$$
整理得:
$${ \left\{ \begin{matrix} a\sum_{j=1}^{N}p_{j}^{2}+b\sum_{j=1}^{N}p_{j}q_{j} + \lambda a =0 \\ a\sum_{j=1}^{N}p_{j}q_{j}+b\sum_{j=1}^{N}q_{j}^{2} + \lambda b =0 \\ a^{2}+b^{2}-1 = 0 \end{matrix} \right. }$$
观察可知,上述方程组可写成如下矩阵形式。对线性代数还有印象的话,可以看出来下面的矩阵表达式其实就是矩阵特征值和特征向量的概念(教科书上的${ \mathbf{A} \mathbf{x}=\mu \mathbf{x} }$)。
$${ \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}}_{\mathbf{S}}\begin{pmatrix} a \\ b \end{pmatrix}=-\lambda \begin{pmatrix} a \\ b \end{pmatrix},需要\begin{pmatrix} a & b \end{pmatrix}\begin{pmatrix} a \\ b \end{pmatrix}=1 }$$
所以上式的解${ \begin{pmatrix} a \\ b \end{pmatrix} }$就是矩阵S的单位特征向量。我们知道,矩阵S一般有两个特征值和特征向量。那么哪个才是让e的值最小的特征向量呢?答案是特征值较小的那个特征向量。下面给出说明。由于恰巧上式中的矩阵S就是(3)式中的矩阵S,所以我们可以把(3)式抄下来并进行变换。下面用${ \lambda^{'} }$表示与${ \begin{pmatrix} a \\ b \end{pmatrix} }$对应的特征值:
$${ \begin{align*} e&=\begin{pmatrix} a & b \end{pmatrix} \mathbf{S}\begin{pmatrix} a \\ b \end{pmatrix} \\ &=\begin{pmatrix} a & b \end{pmatrix} \lambda^{'}\begin{pmatrix} a \\ b \end{pmatrix} \\ &=\lambda^{'} \begin{pmatrix} a & b \end{pmatrix}\begin{pmatrix} a \\ b \end{pmatrix} \\ &=\lambda^{'} \end{align*} }$$
很明显,最小的特征值${ \lambda^{'} }$使e取得最小值。