第8章 关于继黄黑色p站张宇考研数学之后，蓝色p站也可以学习了这码事儿——工程数学：最小二乘法在机器视觉中的应用（实习初步学习笔记）

　　……竟然凑够了三十人hhhh

　　 那大家就来看H文写手讲数学吧。

　　 哎呦沃日你大坝，word文档复制粘贴到p站上不会直接把图导进来，离谱……

　　 先修课程：高等数学（基础）、线性代数（重要！）、计算方法（重要！）、机器视觉与人工智能一类的视觉算法课（不是很必要）。

　　 这节课我们讲《最小二乘法在计算机视觉的一些应用》，全文分为三个部分，第一部分讲：什么是最小二乘法？第二第三部分讲最小二乘法的两个应用。

　　 最小二乘法：背景？性质？用途？

　　 很多时候，我们在工程中会遇到解多元一次方程组的问题，比如我们要解一个三元一次方程组，那很明显，我们需要三组互相线性不相关的方程，这样才能恰好解出方程的解X1、X2、X3。也就是说，当我们要解一个n元一次方程组的时候，一般来说是需要n组互相线性不相关的方程的。

　　 这个性质，我们用矩阵来表示：

　　 A(n*n) * X(n*1) = B(n*1)

　　 其中X就是我们要求的解向量。

　　 但是在工程中，我们遇到的问题往往是：我们所列出的方程，每一组都存在误差，比如模拟到数字的转化误差（比如电表量电压的时候，只能显示小数点后一位，但是实际上的电压可能小数点后面还有好几位）、数据类型转化产生的误差（如double到int，6.5->6，这0.5的误差就出来了）……等等等等。

　　 那这样一来，我拿着这些带着误差的值去求解，那当然是铁求错啊！于是乎，工程师们就想了个办法：既然反正错了，不如我多用几组数据去求解，这样去求一个方差最小的解出来，尽量满足计算要求。

　　 好家伙，这样一来，多组数据共同进入运算，互相去消除误差，听起来很美妙对吧？Wouldn’t that be lovely?

　　 我说婷婷。

　　 比如我们用十组（n=10）数据去求解一个六元一次方程吧。

　　 我们再把刚才的矩阵表达式带一遍，会是怎么样呢？

　　 A(10*10) * X(6*1) = B(10*1)

　　 发现问题了没有！矩阵维度根本不一致！这个矩阵乘法等式根本不成立！

　　 那有些同学就问了啊，A它不是一个10*6的矩阵吗？怎么到你这就成了10*10了呢？

　　 哎我说这讲课不是胡说，编课不是瞎编。我们再写一下矩阵表达式：

　　 A(10*6) * X(6*1) = B(10*1)

　　 好！矩阵等式至少在维度上成立了！很有精神！大功告成？

　　 我说婷婷。

　　 10组数据根本解不出来6个未知数，这是一个超定问题（已知的量的个数超过求解需要的量的个数）。

　　 那怎么办？我总不能一下子扔掉后四个已知量吧？且不说这前边的六位兄弟卖了后面的四个兄弟，他们的心里冷~冰~冰~，还有就是我们一开始的目的：多组数据共同进入运算，互相去消除误差没有达到，可以说没有充分地利用已知条件。

　　 那这个时候就要到最小二乘法大显神威了：我们在上一条等式的两边乘以A的转置：

　　 AT(6*6)*A(10*6) * X(6*1) = AT(6*6)*B(10*1)

　　 化简、合并：

　　 A\\u0027(6*6)*X(6*1)=B\\u0027(6*1)

　　 这样就化为了多元一次方程组的解的一般形式了。

　　 应用1：相机标定

　　 相机标定（Camera Calibration），是给定一个真实坐标系，然后将像素坐标转化到这个坐标系上的问题。我在实习中遇到的问题是：给一张含有DN码和棋盘格点的图片（一定是平拍的，不会倾斜，因此不涉及相机倾角的问题，因此可以不用张正友标定发那么复杂的方法），读DN码得到DN码的方向和DN码中心的真实坐标，并且依据方向建立基础坐标系、并利用多个棋盘的格点去确定更精确的图像坐标系到真实坐标的变换矩阵，最后再用这个变换矩阵去实现相机坐标系到真实坐标系的坐标变换。

　　 如图。这个大蓝圈儿的中心就是用真实坐标（13,21）做反变换得到的相机坐标。别的圈儿是我中间debug的时候画的。（二维码的中心坐标是（10,20））

　　 别的一些对应角点、角点排序等问题并不是我们今天的主题，我们不过多赘述，我们今天只讲opencv下最小二乘法的应用。

　　 我们三步走：

　　 第一步，根据实际模型填写矩阵A、B：

　　 第二步：A和B同时乘以A的转置，得A’、B’。

　　 第三步：用A’·X=B’解得X

　　 那么最小二乘法应用完了，最核心的问题就解决了，剩下就是将abcdef装回矩阵，然后用

　　 去做坐标变换罢了。

　　 应用2：小弧度圆的拟合

　　 给定一张图片，拍摄角度不超过10°的一个圆弧，对应的半径在像素上是从至少1000到约2600不等，圆心的位置完全不确定，只知道圆弧上的一组最广分布在±2°上的点和它们不超过±3°的对应旋转关系。

　　 由于第二个任务目前只是在让我做一个初步的算法，还没有具体到实际的可视化，因此不能像第上一个应用那样去给大家直接展示成果，很多时候只能讲理论。

　　 我们假设：原圆心为（x0，y0），一共有a个点，和它们对应的b个旋转关系。（这里的原圆心是（500，900），每个亚像素的噪声在x和y向的都是方差为0.5/2.58的AWGN（Addictive White Gaussian Noise），即以99%概率分布在±0.5，由于这个噪声对于小弧度的圆来说实在过于致命，后来我又用最小二乘法拟合了一条补正的二次曲线，这才有了后面比较精准的解算结果，不过后面的过程我们不多赘述，重点讲第一步这个“粗略”圆心是怎么获得的）。

　　 话不多说，三板斧开砍。

　　 第一步，根据实际模型填写矩阵A、B：

　　 这里的实际模型比第一例还要复杂很多，甚至矩阵维度都不是固定的——如果有a个点，A的维度将会是（a*b, a+2）。这里以a=3，b=3为例。

　　 注意这里建模的时候，建立在“三点隶属同一圆心”这一前提上是十分重要的。Halcon的拟合方法我也用了，它并没有这样的功能，只能用一个点去不停地旋。这种拟合出来的圆心是十分离谱的，完全不能用的。

　　 第二步：A和B同时乘以A的转置，得A’、B’。

　　 第三步：用A’·X=B’解得X

　　 然后，我们取X的前两项，乘以-0.5，就是圆心的坐标了。

　　 总结：

　　 最小二乘法三步走：填矩阵、两边乘转置、解得。说白了就是这么简单。

　　 绝了，经典不务正业，经典一周五千字给金主小刀刮痧但摸鱼却摸的很起劲23333