青铜时代 —— 图像处理

• 研究图像中的变量，(宽，高，通道，像素值)
• 图像处理就是如何怎样操作这些变量

1 坐标维变换

图像拼接

• 完整的拼接包括：特征点检测匹配和筛选，计算仿射变换和透射变换矩阵，进行拼接（计算出roi，用一个大的矩阵盛放）
• opencv中的代码：

Stitcher X = Stitcher::createDefault(true);
Stitcher::Status s = X.stitch(imags, resMat)
if (s != Stitcher::OK) // ERROR

图像变换

• 仿射变换/透射变换 参考
  
  
  
• LogPolar变换

图像采样/插值

• 最近邻插值，双线性插值（默认设置），双三次插值

最近邻插值

• 在四邻像素中，将距离待求像素最近的邻像素灰度赋给待求像素。

双线性插值算法

三次插值法

• 双立方插值法中，我们选取的是最近的16个像素点作为计算目标图像B(X,Y)处像素值的参数。如图所示：
  
• 如下图：当a取不同值时可以用来逼近不同的样条函数（常用值-0.5, -0.75）
  

图像金字塔

2 通道维变换

图像灰度化 & 颜色通道变换

• cvtColoer(A,B,COLOR_X,dstCn);
• COLOR_X: BGR / RGB / BGRA / RGBA / GRAY / BGR555 / BGR565 / XYZ / YCrCb / Lab / HSV / HLS / Luv

Lab

• 由一个亮度通道（channel）和两个颜色通道组成的。在Lab颜色空间中，每个颜色用L、a、b三个数字表示，各个分量的含义是这样的：
  • L*代表亮度
  • a*代表从绿色到红色的分量
  • b*代表从蓝色到黄色的分量
• Lab有个很好的特性——设备无关（device-independent）。也就是说，在给定了颜色空间白点（white point）（下图中表示了一种颜色空间的白点）之后，这个颜色空间就能明确地确定各个颜色是如何被创建和显示的，与使用的显示介质没有关系
  
• Lab色域很大，有一大部分已经超出了人类视觉范围，也就不能称之为“颜色”了
• Lab定义的是相对于白点的颜色，只有定义完白点是什么颜色（比如定义为CIE standard illuminant D50），我们才能知道其他的颜色

Luv

• L是亮度，u和v是色度坐标。对于一般的图像，u和v*的取值范围为-100到+100，亮度为0到100。
• 它的计算公式可以由CIEXYZ通过非线性计算得到。
• CIE xyY空间是由 XYZ值导出的空间， Y 是颜色的明度或亮度。x和y是CIE xy色度坐标，它们是所有三个三色刺激值 X、Y 和 Z 的函数所规范化的三个值中的两个：
• -XYZ (CIE 1931)
  
• Luv (CIE 1976)
  
  

HSV

• 色相、饱和度和亮度
• hue,saturation,value

• hue表示颜色的相位角，取值范围是0—360;S表示颜色的饱和度；S为一比例值，范围从0到1，它表示成所选颜色的纯度和该颜色最大的纯度之间的比率，通俗点讲,S表示的是某种颜色的“纯度”, S取值越大，表示色彩越纯，取值越小，表示色彩越灰。V表示色彩的明亮程度，范围从0到1。V等于0表示圆锥的底部定点，也就是黑色，V等于1表示圆锥的顶面，当V=1并且S=0时表示纯白色；

YCbCr

• Y就是所谓的流明（luminance），Y’表示光的浓度且为非线性，使用伽马修正（gamma correction）编码处理，而CB和CR则为蓝色和红色的浓度偏移量成分
• YCbCr不是一种绝对色彩空间，是YUV压缩和偏移的版本。YCbCr的Y与YUV中的Y含义一致，Cb和Cr与UV同样都指色彩，Cb指蓝色色度，Cr指红色色度，

3 像素值变换

3.1 单像素值的像素值变换

阈值化 > 二值图

• Otsu算法
  • 思路：最大类间方差
  • 列出灰度直方图 + 穷举搜索 + 类间方差最大
  • 记T为前景与背景的分割圖值, 前景点数占图像比例为$w_0$，平均灰度为 $u_0$ 背景点数占图像比例为 $w_1$ 平均灰度为$u_1$，图像的总平均灰度为$u$，前景和背景图象的方差则有：
    $$\begin{array}{c}u=w_0\times u_0+w_1\times u_1\\ g=w_0\times {\left(u_0-u\right)}^2+w_1\times {\left(u_1-u\right)}^2\\ g=w_0\times w_1\times {\left(u_0-u_1\right)}^2\\ g=\frac{w_0}{1-w_0}\times {\left(u_0-u\right)}^2\end{array}$$
• 自适应阈值
  • 阈值不是全局的，是在一个块窗口内取均值或加高斯窗的均值作为该点的阈值。利用了局部信息！

Gamma校正

线性/对数变换

灰度分层/拉伸

对比度拉伸

积分图

傅里叶变换

$$\begin{gathered} F(u,v)=\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}f(x,y)e^{-j2\pi \left(\frac{\mathrm{u}\mathrm{x}}{M}+\frac{vy}{N}\right)} \\ f(x,y)=\frac{1}{MN}\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}F(u,v)e^{j2\pi \left(\frac{\mathrm{u}\mathrm{x}}{\mathrm{M}}+\frac{vy}{N}\right)} \end{gathered}$$

3.2 全局统计信息的像素值变换

直方图

• 计算：cv::calcHist()

均衡化

• Histogram Equalizatio
• void cv::equalizeHist(InputArray src, OutputArray dst)
• 原始图像的灰度直方图从比较集中的某个灰度区间变成在全部灰度范围内的均匀分布
• 直方图均衡化就是对图像进行非线性拉伸，重新分配图像像素值，使一定灰度范围内的像素数量大致相同
• 变换函数要求是单值单调增加
  
  
  
  
  

规定化（匹配）

• Histogram Normalization/Matching
• 指对一副图像进行变换，使其直方图与另一幅图像的直方图或特定函数形式的直方图进行匹配。

流程

• 对图像A和B进行均衡化操作
• 得到A、B的直方图，并分别计算累积直方图
• 得到B到A的累计直方图的映射关系，如下图得到规定结果：
  

直方图距离

• CompareHist
  

反向投影

模板匹配

共生矩阵

• 取图像N×N中任意一点 （x，y）及偏离它的另一点 （x+a，y+b），设该点对的灰度值为 （g1，g2）。令点（x，y） 在整个画面上移动，则会得到各种 （g1，g2）值，设灰度值的级数为 k，则（g1，g2） 的组合共有 k 的平方种。对于整个图像画面，统计出每一种 （g1，g2）值出现的次数，然后排列成一个方阵，再用（g1，g2） 出现的总次数将它们归一化为出现的概率P（g1，g2） ，这样的方阵称为灰度共生矩阵。

• 共生矩阵常用的特征包含ASM能量、对比度、熵、自相关，假如我们另G为计算得出的灰度共生矩阵，那么这些特征可以作如下的表示。

ASM 能量

• angular second moment，即角二节矩，共生矩阵G每个元素的平方之和。，如果G的元素分布比较集中（即特定空间关系下像素点对的灰度值对出现概率大），那么特征ASM能量的值就比较大；如果G的元素分布比较均匀，那么得到的ASM能量的值就比较小。所以ASM反映的图像的灰度分布均匀程度以及纹理的粗细程度。
  $$ASM=\sum _{i=1}^k\sum _{j=1}^k(G(i,j){)}^2$$

对比度

• 对比度，直接反映了某个像素值及其领域像素值的亮度的对比情况。如果偏离对角线的元素有较大值，即图像亮度值变化很快，则CON会有较大取值，这也符合对比度的定义。其中 。反映了图像的清晰度和纹理沟纹深浅的程度。纹理沟纹越深，其对比度越大，视觉效果越清晰；反之，对比度小，则沟纹浅，效果模糊。灰度差即对比度大的象素对越多，这个值越大。灰度公生矩阵中远离对角线的元素值越大，CON越大
  $$CON=\sum _{n=0}^{k-1}{n}^2\left|\sum _{|i-j=n}G(i,j)\right|$$

逆差矩

• 逆差矩也称反差分矩阵，灰度共生矩阵对角元素有较大值，IDM就会取较大的值。因此连续灰度的图像会有较大IDM值。逆差矩反映图像纹理的同质性，度量图像纹理局部变化的多少。其值大则说明图像纹理的不同区域间缺少变化，局部非常均匀。
  $$IDM=\sum _{i=1}^k\sum _{j=1}^k\frac{G(i,j)}{1+(i-j{)}^2}$$

熵

• 熵，熵越大表示信息越混乱。纹理信息也是信息的一种，而且纹理信息是图像的一个随机性的度量。所以这里的熵反映了图像中的纹理复杂程度（不均匀），熵值越大纹理越复杂。
  $$ENT=-\sum _{i=1}^k\sum _{j=1}^{k}G(i,j)\log G(i,j)$$

自相关

• 自相关，自相关反应了图像纹理的一致性。如果图像中有水平方向纹理，则水平方向矩阵的COR大于其余矩阵的COR值。它度量空间灰度共生矩阵元素在行或列方向上的相似程度，因此，相关值大小反映了图像中局部灰度相关性。当矩阵元素值均匀相等时，相关值就大;相反，如果矩阵像元值相差很大则相关值小。
  $$COR=\sum _{i=1}^k\sum _{j=1}^k\frac{(ij)G(i,j)-u_i{u}_j}{s_i{s}_j}$$
  
• 当距离差分值（a,b）取四种值的时候，可以综合得到向量 : ASM1,CON1,IDM1,ENT1,COR1,…,ASM4,CON4,IDM4,ENT4,COR4

图像复原

3.3 局部邻域信息的像素值变换

形态学

图像平滑

图像去噪

图像稠密特征

HOG

• Gamma校正和灰度化
  
• 计算图像像素梯度图 (Sobel算子)
• 计算梯度直方图:联合一个8×8的小格子内部一些像素，计算其梯度幅度和梯度方向的统计直方图，这样一来就可以以这个梯度直方图来代替原本庞大的矩阵。每个像素有一个梯度幅度和梯度方向两个取值，那么一个小格子一共有8×8×2=128个取值；[0,180]分为9组
• Block归一化：HOG在选取8×8为一个单元格的基础之上，再以2×2个单元格为一组，称为block。作者提出要对block进行归一化，由于每个单元cell有9个向量，2×2个单元格则有36个向量，需要对这36个向量进行归一化
• HOG特征描述：每一个16×16大小的block将会得到36大小的vector，一个64×128大小的图像，按照上图的方式提取block，将会有7个水平位置和15个竖直位可以取得，所以一共有7×15=105。图片HOG特征的总维度36×105=3780！

fHOG

• 建立像素级特征映射：计算梯度幅值和方向！可以使用对方向敏感 B1（0-360度）也可以使用对方向不敏感 B2（0-180度）
  $$\begin{array}{l}{B}_1(x,y)=\mathrm{round}\left(\frac{p\theta (x,y)}{2\pi }\right)\mathrm{m}\mathrm{o}\mathrm{d}p\\ B_2(x,y)=\mathrm{round}\left(\frac{p\theta (x,y)}{\pi }\right)\mathrm{m}\mathrm{o}\mathrm{d}p\end{array}$$
• 空间聚合：增加36维！利用三线插值进行处理，这种处理也是合理的。能进一步较小混叠效应，通过插值可以使得每个像素对其周围的四个cell的特征向量产生贡献！
  
• 归一化和截断。四个归一化因子，考虑中间的cell5，绿色的四个方式！
  
  
• FHOG PCA降维：作者经过实验，发现：（PCA）由前11个主特征向量定义的线性子空间基本包含了hog特征的所有信息，并且用降维之后的特征在他们的任务（目标检测）中取得了和用36维特征一样的结果。然后作者采用了另一种方法：13维特征有一个比较简单的解释，9个方向特征以及反映cell周围区域能量的四个特征。取得一样的结果！
  
• 最终得到31维的特征向量：27个在不同归一化因子上的累加和（列和）以及4个在不同方向上的累加和（行和），27维中包含了27个bin通道，其中18个对方向敏感，9个对方向不敏感，4维分别捕获了当前cell周围4个cell组成的梯度能量。如上图所示！

LBP

4 点检测

4.1 角点检测

Moravec

• 以该点为中心的5x5窗口，计算其灰度总和与右移，下移，右下移和左下移的灰度总和的最小差值作为角点响应值CRF，之后接人工阈值判断
• 问题：对边缘响应很强 + 方向局限性(没有旋转不变性) + 对噪声敏感
• Opencv里面没有实现！

Harris

• Opencv里面的goodFeaturesToTrack ( 或 cornerHarris() ) 实现了Harris和Shi Tomasi角点检测算法
• 计算(x,y)平移后的自相关性：
  $$c(x,y;\Delta x,\Delta y)=\sum _{(u,v)\in W(x,y)}w(u,v)(I(u,v)-I(u+\Delta x,v+\Delta y){)}^2$$
• W是以(x,y)为中心的窗，w为加权函数，可以常数可以高斯加权函数
• $I$的一阶近似：
  $$\begin{gathered} I(u+\Delta x,v+\Delta y)=I(u,v)+I_x(u,v)\Delta x+I_y(u,v)\Delta y+O\left(\Delta x^2,\Delta y^2\right) \\ \approx I(u,v)+I_x(u,v)\Delta x+I_y(u,v)\Delta y \end{gathered}$$
• 即：
  $$\begin{gathered} c(x,y;\Delta x,\Delta y)\approx \sum _w{\left(I_x(u,v)\Delta x+I_y(u,v)\Delta y\right)}^2= \\ [\Delta x,\Delta y]M(x,y)\left[\begin{array}{c}\Delta x\\ \Delta y\end{array}\right] \\ M(x,y)=\sum _w\left[\begin{array}{cc}{I}_x(x,y{)}^2& I_x(x,y)I_y(x,y)\\ I_x(x,y)I_y(x,y)& I_y(x,y{)}^2\end{array}\right]= \\ \left[\begin{array}{cc}{\sum }_w{I}_x(x,y{)}^2& {\sum }_w{I}_x(x,y)I_y(x,y)\\ {\sum }_w{I}_x(x,y)I_y(x,y)& {\sum }_w{I}_y(x,y{)}^2\end{array}\right]=\left[\begin{array}{cc}A& C\\ C& B\end{array}\right] \end{gathered}$$
• 即二次项函数可以二次表示，就是一个椭圆函数！其长短轴由M的特征值决定！
• 分析特征值的大小：① 一个大一个小，就是直线 ② 都小，就是图像中的平面（变化小） ③ 都大，就是角点！
• 定义角点响应值：$R=\det M-\alpha (\mathrm{trace}M{)}^2$；$\alpha =[0.04,0.06]$(和检测出的点个数成反比！)
• 流程：设置窗口和加权窗 + 计算窗口内的梯度 + 计算窗口内各点的M矩阵 + 计算每个点的响应值进行阈值化 + 在3x3和5x5窗口进行最大值抑制！
  • Harris对亮度和对比度变化不太敏感
  • Harris具有旋转不变性
  • Harris不具有尺度不变性，cornerHairrs()
• 多尺度Harris角点检测：
  $$\begin{gathered} \mathbf{M}=\mu \left(x,{\sigma }_I,{\sigma }_D\right)= \\ {\sigma }_D^{2}g\left({\sigma }_I\right)\otimes \left[\begin{array}{cc}{L}_x^2\left(x,{\sigma }_D\right)& L_x{L}_y\left(x,{\sigma }_D\right)\\ L_x{L}_y\left(x,{\sigma }_D\right)& L_y^2\left(x,{\sigma }_D\right)\end{array}\right] \end{gathered}$$
  • 即先用高斯核g进行进行加窗，而${\sigma }_D$为微分尺度或局部尺度, 它是决定角点附近微分偵变化的变量。

Shi Tomasi

• Good Features to Track
• Harris进行改进：Shi Tomasi 发现，角点的稳定性其实和矩阵 M 的较小特征值有关，于是直接用较小的那个特征值作为分数。这样就不用调整k值了 $R=\min \left({\lambda }_1,{\lambda }_2\right)$，之后也是阈值化！

Fast

void cv::FAST	(	
InputArray 	image,
std::vector< KeyPoint > & 	keypoints,
int 	threshold,
bool 	nonmaxSuppression,
int 	type 
)	

• 定义一个阈值。计算p1、p9、p5、p13与中心p的像素差，若它们的绝对值有至少3个超过阈值，则当做候选角点，再进行下一步考察；否则，不可能是角点；
• 若p是候选点，则计算p1到p16这16个点与中心p的像素差，若它们有至少连续9个超过阈值，则是角点；否则，不可能是角点。
• 对图像进行非极大值抑制：计算特征点出的FAST得分值（即score值，也即s值），判断以特征点p为中心的一个邻域（如3x3或5x5）！
  $$Score=\max \{\begin{array}{l}\sum (\text{ pixel values }-p)\text{ if (value }-p)>t\\ \sum (p-\text{ pixel values })\text{ if }(p-\text{ value })>t\end{array}$$

4.2 关键点检测和描述

sift

static Ptr<SIFT> cv::SIFT::create	(	int 	nfeatures = 0,
int 	nOctaveLayers = 3,
double 	contrastThreshold = 0.04,
double 	edgeThreshold = 10,
double 	sigma = 1.6 
)	

• 流程：
  • 生成高斯金字塔，并由此生成差分高斯金字塔（DoG），该过程完成尺度空间的构建，以便后续的空间极值点检测。
  • 稳定关键点的精确定位、去除不稳定特征点
  • 稳定关键点信息分配
  • 关键点匹配

高斯图像金字塔

• 高斯窗口$(6\sigma +1{)}^2$
• Tony Lindeberg指出尺度规范化的LoG(Laplacion of Gaussian)算子具有真正的尺度不变性，Lowe使用高斯差分金字塔近似LoG算子，在尺度空间检测稳定的关键点
• 尺度空间：高斯模糊 + 降采样
• 金字塔每层多张图像(Octave)，金字塔每层只有一组图像，组数和金字塔层数相等，使用公式(3-3)计算，每组含有多张(也叫层Interval)图像

第0组(即第-1组)	0 1 2 3 4 5
第1组	6 7 8 9 10 11
第2组	第一张图根据第一组中索引为9的图片降采样

空间极值点检测

• 每组需S+3层图像
  • 由于要在相邻尺度进行比较
  • 构建尺度空间需确定的参数： —尺度空间坐标$\sigma$；O组(octave)数；S组内层数
    $$\sigma (o,s)={\sigma }_0{2}^{o+\frac{s}{S}}o\in [0,\dots ,O-1],s\in [0,\dots ,S+2]$$
• 关键点定位
  
  $$D(X)=D+\frac{\partial D^T}{\partial X}X+\frac{1}{2}{X}^T\frac{{\partial }^{2}D}{\partial X^2}X$$
  • 求导=0，得到极值点的偏移量：$\hat{X}=-\frac{{\partial }^2{D}^{-1}}{\partial X^2}\frac{\partial D}{\partial X}$
  • 其中，$\hat{X}=(x,y,\sigma {)}^T$是相对于插值中心的偏移量！当大于0.5则代表需要改变点位置，反复插值迭代到收敛！（如5次）
  • $|D(x)|$小的点容易受噪声影响，将小于0.03的点删掉！
• 消除边缘响应
  $$\begin{gathered} H=\left[\begin{array}{ll}{D}_{xx}& D_{xy}\\ D_{xy}& D_{xy}\end{array}\right] \\ \frac{\mathrm{Tr}(H{)}^2}{\mathrm{Det}(H)}=\frac{(\alpha +\beta {)}^2}{\alpha \beta }=\frac{(r\beta +\beta {)}^2}{r{\beta }^2}=\frac{(r+1{)}^2}{r} \end{gathered}$$
  • 取r=10，小于$11^2/10$的点认为是边缘点，删掉！

关键点方向分配

• 对于在DOG金字塔中检测出的关键点点，采集其所在高斯金字塔图像$3\sigma$邻域窗口内像素的梯度和方向分布特征。
• 梯度模进行高斯加权，$\sigma =1.5{\sigma }_{oct}$，按照$3\sigma$原则，窗口尺寸$4.5{\sigma }_{oct}$
• 梯度直方图将0~360度的方向范围分为36个柱(bins)，其中每柱10度
  
  • 以直方图中最大值作为该关键点的主方向。为了增强匹配的鲁棒性，只保留峰值大于主方向峰值80％的方向作为该关键点的辅方向

关键点特征描述

• 确定计算描述子所需的图像区域(4x4)
  $$radius=\frac{3{\sigma }_{-}oct\times \sqrt{2}\times (d+1)}{2}$$
  
• 将坐标轴旋转为关键点的方向，以确保旋转不变性
  $$\begin{gathered} \left(\begin{array}{l}{x}^{\prime }\\ y^{\prime }\end{array}\right)=\left(\begin{array}{cc}\cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{array}\right)\left(\begin{array}{l}x\\ y\end{array}\right)(x,y\in [-radius,radius]) \\ \left(\begin{array}{l}{x}^{\prime \prime }\\ y^{\prime \prime }\end{array}\right)=\frac{1}{3{\sigma }_{-}\mathrm{oct}}\left(\begin{array}{l}{x}^{\prime }\\ y^{\prime }\end{array}\right)+\frac{d}{2} \end{gathered}$$
  
• 将邻域内的采样点分配到对应的子区域内，将子区域内的梯度值分配到8个方向上，计算其权值：owe建议子区域的像素的梯度大小按$\sigma =0.5d$的高斯加权计算：$w=m(a+x,b+y)e^{-\frac{{\left(x^{\prime }\right)}^2+{\left(y^{\prime }\right)}^2}{2\times (0.5d{)}^2}}$
• 插值计算每个种子点八个方向的梯度：双线性插值
  weight $=w^{\ast }d{r}^k\ast (1-dr{)}^{1-k}\ast d{c}^{m\ast }(1-dc{)}^{1-m}\ast d{o}^n\ast (1-do{)}^{1-n}$ ；{0,1}
  
• 特征向量归一化处理：对4x4x8=128维向量进行归一化！
• 描述子向量门限：截断大于0.2的梯度值
• 按特征点的尺度对特征描述向量进行排序
• 缺点：实时性不高 + 特征点有时候小 + 边缘光滑无法准确提取特征点！

surf (Speeded Up Robust Features)

• cv::xfeatures2d::SURF

Hessian检测

$$\mathcal{H}(\mathbf{x},\sigma )=\left[\begin{array}{ll}{L}_{xx}(\mathbf{x},\sigma )& L_{xy}(\mathbf{x},\sigma )\\ L_{xy}(\mathbf{x},\sigma )& L_{yy}(\mathbf{x},\sigma )\end{array}\right]$$

• $L_{xx}$表示图像经过高斯二阶LoG卷积得到，离散化近似（盒函数）：
  
• 那么就可以采用积分图像加速Hessian矩阵的计算！得到近似值$D_{xx}$等，计算指标：$\det \left({\mathcal{H}}_{\text{approx }}\right)=D_{xx}{D}_{yy}-{\left(0.9D_{xy}\right)}^2$，0.9是为了补偿近似的误差！

SURF的尺度空间

• SURF采用的盒函数+积分图，不需要不断的高斯模糊核子采样，而是采用不同尺寸的核函数+采用 来计算图像金字塔！
• 99，1515，2121，2727 二倍增大
• 左图是sift算法，其是图像大小减少，而模板不变（这里只是指每组间，组内层之间还是要变的）。而SURF算法（右图）刚好相反，其是图像大小不变，而模板大小扩大
  

兴趣点主方向获得

• SURF方法则是通过计算其在x,y方向上的haar-wavelet响应，这是在兴趣点周围一个6s半径大小的圆形区域内。当然小波变换的大小也同尺度参数s有关，其步长为s，其大小为4s
  
• 60度

SURF描述子

• SURF也是通过建立兴趣点附近区域内的信息来作为描述子的
• SURF首先在兴趣点附近建立一个20s大小的方形区域，为了获得旋转不变性，同sift算法一样，我们需要将其先旋转到主方向，然后再将方形区域划分成16个（44）子域。对每个子域（其大小为5s5s）我们计算25（5*5）个空间归一化的采样点的Haar小波响应dx和dy
  
• 之后我们将每个子区域（共4*4）的dx,dy相加，因此每个区域都有一个描述子（如下式），为了增加鲁棒性，我们可以给描述子再添加高斯权重（尺度为3.3s，以兴趣点为中心）
• 所以最后在所有的16个子区域内的四位描述子结合，将得到该兴趣点的64位描述子
  
• 也建立128位的SURF描述子，其将原来小波的结果再细分，比如dx的和将根据dy的符号，分成了两类，所以此时每个子区域内都有8个分量，SURF-128有非常好效果

快速索引匹配

• 我们发现兴趣点其Laplacian（Hessian矩阵的迹）的符号（即正负）可以将区分相同对比形状的不同区域，如黑暗背景下的白斑和相反的情况。
• 只需要为其中一个建立描述子，而给另一个索引，而在匹配过程中，只要比较一个描述子，就能确定两个位置
  

brief

• cv::BriefDescriptorExtractor()
• BRIEF算法计算出来的是一个二进制串的特征描述符
• n取128、256或512，opencv默认为256

流程

• 特征点SxS的区域内选取点对的方法，和中心点比较灰度大小，进行01赋值！
  • 在图像块内平均采样；
  • p和q都符合$(0,S^2/25)$的高斯分布；论文认为这个好 √
  • p符合$(0,S^2/25)$的高斯分布，而q符合$(0,S^2/100)$的高斯分布；
  • 在空间量化极坐标下的离散位置随机采样；
  • 把p固定为$(0,0)$，q在周围平均采样
    
• 旋转大于30°后，BRIEF的匹配率快速降到0左右

orb

• Oriented FAST and Rotated BRIEF

oFast特征提取

• 粗提取：和原Fast特征提取一致！
• 机器学习的方法筛选最优特征点. 使用ID3决策树训练筛选最优的FAST特征点
• 非极大值抑制去除局部较密集特征点。响应计算，特征点P和其周围16个特征点偏差的绝对值和
• 特征点的尺度不变形。比例因子scaleFactor（opencv默认为1.2）和金字塔的层数nlevels（Opencv默认为8）
• 特征点的旋转不变性。ORB算法提出使用矩（moment）法来确定FAST特征点的方向。也就是说通过矩来计算特征点以r为半径范围内的质心，特征点坐标到质心形成一个向量作为该特征点的方向

rBRIEF特征描述

• steered BRIEF
  • 在BRIEF特征描述的基础上加入旋转因子改进的
  • 将brief采用的点对旋转至主方向（oFast提供）
• rBRIEF
  • sBRIFE降低了描述符的可区分性
    
• rBRIEF
  • 为了解决描述子的可区分性和相关性的问题，ORB论文中没有使用5种方法中的任意一种，而是使用统计学习的方法来重新选择点对集合。
  • 其31x31邻域
  • 这里在对图像进行高斯平滑之后，使用邻域中的某个点的5x5邻域灰度平均值来代替某个点对的值，进而比较点对的大小。这样特征值更加具备抗噪性。另外可以使用积分图像加快求取5x5邻域灰度平均值的速度。
  • 31x31的邻域内共有(31-5+1)x(31-5+1)=729个这样的子窗口
  • 取点对的方法共有M=265356种
  • 1）在300k特征点的每个31x31邻域内按M种方法取点对，比较点对大小，形成一个300kxM的二进制矩阵Q。矩阵的每一列代表300k个点按某种取法得到的二进制数。
  • 2）对Q矩阵的每一列求取平均值，按照平均值到0.5的距离大小重新对Q矩阵的列向量排序，形成矩阵T。
  • 3）将T的第一列向量放到R中。
  • 4）取T的下一列向量和R中的所有列向量计算相关性，如果相关系数小于设定的阈值，则将T中的该列向量移至R中。
  • 5）按照4）的方式不断进行操作，直到R中的向量数量为256。通过这种方法就选取了这256种取点对的方法。这就是rBRIEF算法。

brisk

• 《BRISK:Binary Robust Invariant Scalable Keypoints》
• 速度比较：SIFT<SURF<BRISK<FREAK<ORB

特征点检测

• BRISK算法主要利用FAST9-16进行特征点检测
• 尺度空间：n个octave层（用ci表示）和n个intra-octave层（用di表示），文章中n=4，i={0,1,…,n-1}
  • octave层的产生：c0层就是img原图像，c1层是c0层的2倍下采样，c2层是c1层的2倍下采样，以此类推
  • intra-octave层的产生：d0层是img的1.5倍下采样，d1层是d0层的2倍下采样（即img的2*1.5倍下采样），d2层是d1层的2倍下采样，以此类推
  • $t\left(c_i\right)=2^i,t\left(d_i\right)=2^i\cdot 1.5$
    
• 特征点检测：这8张图进行FAST9-16角点检测，得到具有角点信息的8张图，对原图像img进行一次FAST5-8角点检测（当做d(-1)层，虚拟层），总共会得到9幅有角点信息的图像
• 非极大值抑制：特征点在位置空间（8邻域点）和尺度空间（上下层2x9个点），共26个邻域点的FAST的得分值要最大，否则不能当做特征点；此时得到的极值点还比较粗糙，需要进一步精确定位
• 亚像素插值：类似SIFT，上下两层进行插值！

特征点描述

• 均匀采样模式：以特征点为中心，构建不同半径的同心圆，在每个圆上获取一定数目的等间隔采样点（所有采样点包括特征点，一共N=60），由于这种邻域采样模式会引起混叠效应，所以需要对同心圆上的采样点进行高斯滤波。 蓝圈是采样中心，红圈表示滤波半径！
  
• 局部梯度计算：两两一对，有$N(N-1)/2$组合方法
  $$\mathbf{g}\left({\mathbf{p}}_i,{\mathbf{p}}_j\right)=\left({\mathbf{p}}_j-{\mathbf{p}}_i\right)\cdot \frac{I\left({\mathbf{p}}_j,{\sigma }_j\right)-I\left({\mathbf{p}}_i,{\sigma }_i\right)}{{\Vert {\mathbf{p}}_j-{\mathbf{p}}_i\Vert }^2}$$
• 短距离点对子集、长距离点对子集（L个）
  S = { ( p i , p j ) ∈ A ∣ ∥ p j − p i ∥ < δ max ⁡ } ⊆ A L = { ( p i , p j ) ∈ A ∣ ∥ p j − p i ∥ > δ min ⁡ } ⊆ A δ max ⁡ = 9.75 t , δ min ⁡ = 13.67 t \begin{array}{l} \mathcal{S}=\left\{\left(\mathbf{p}_{i}, \mathbf{p}_{j}\right) \in \mathcal{A} |\left\|\mathbf{p}_{j}-\mathbf{p}_{i}\right\|<\delta_{\max }\right\} \subseteq \mathcal{A} \\ \mathcal{L}=\left\{\left(\mathbf{p}_{i}, \mathbf{p}_{j}\right) \in \mathcal{A} |\left\|\mathbf{p}_{j}-\mathbf{p}_{i}\right\|>\delta_{\min }\right\} \subseteq \mathcal{A} \end{array}\\ \delta_{\max }=9.75 t, \delta_{\min }=13.67 t S={(pi​,pj​)∈A∣∥pj​−pi​∥<δmax​}⊆AL={(pi​,pj​)∈A∣∥pj​−pi​∥>δmin​}⊆A​δmax​=9.75t,δmin​=13.67t
• 特征点的主方向(此处用了长距离子集)
  $$\begin{array}{c}\mathbf{g}=\left(\begin{array}{c}{g}_x\\ g_y\end{array}\right)=\frac{1}{L}\cdot {\sum }_{\left({\mathbf{p}}_i,{\mathbf{p}}_j\right)\in \mathcal{L}}\mathbf{g}\left({\mathbf{p}}_i,{\mathbf{p}}_j\right)\\ \alpha =\arctan 2\left(g_y,g_x\right)\end{array}$$
• BRISK描述子是二进制的特征，由采样点集合可得到N(N-1)/2对采样点对，就可以得到N(N-1)/2个距离的集合，考虑其中短距离子集中的512个短距离点对，进行二进制编码：
  $$\begin{array}{l}b=\{\begin{array}{ll}1,& I\left({\mathrm{p}}_j^{\alpha },{\sigma }_j\right)>I\left({\mathrm{p}}_i^{\alpha },{\sigma }_i\right)\\ 0,& \text{ otherwise }\end{array}\\ \forall \left({\mathrm{p}}_i^{\alpha },{\mathrm{p}}_j^{\alpha }\right)\in \mathcal{S}\end{array}$$
• 其中, $I\left({\mathbf{p}}_j^{\alpha },{\sigma }_j\right)$ 表示经过旋转$\alpha$角度后的新的采样点。由此可得到512Bit的二进制编码, 也就是64个字节(BRISK64)

freak

• 《FREAK: Fast Retina Keypoint》
• cv::xfeatures2d::FREAK
• 与BRISK算法非常相似，个人觉得就是在BRISK算法上的改进

采样模式

特征描述

• 二进制
  $$\begin{array}{c}F={\sum }_{0\le a<N}{2}^{a}T\left(P_a\right)\\ T\left(P_a\right)=\{\begin{array}{ll}1& \text{ if }(I\left(P_a^{r_1}\right)-I\left(P_a^{r_2}\right)>0\\ 0& \text{ otherwise }\end{array}\end{array}$$
• $I\left(P_a^{r_1}\right)$ 表示采样点对Pa中前一个采样点的像素值，$I\left(P_a^{r_2}\right)$ 表示后一个采样点的像素值
• 降维
• 1、建立矩阵D，D的每一行是一个FREAK二进制描述符，即每一行有N个元素；在上图的采样模式中公有43个采样点，可以产生N=43*(43-1)/2=903个采样点对，也就是说矩阵D有903行列（修改于：2015-11-15）；
• 2、对矩阵D的每一列计算其均值，由于D中元素都是0/1分布的，均值在0.5附近说明该列具有高的方差；
• 3、每一列都有一个均值，以均值最接近0.5的排在第一位，均值离0.5越远的排在越靠后，对列进行排序；
• 4、选取前512列作为最终的二进制描述符。(也可以是256、128、64、32等)

特征方向

• 有43个采样点，可产生43*(43-1)/2=903个采样点对，FREAK算法选取其中45个长的、对称的采样点对来提取特征点的方向
  
• O表示局部梯度信息，M表示采样点对个数，G表示采样点对集合，Po表示采样点对的位置
  $$O=\frac{1}{M}\sum _{P_o\in G}\left(I\left(P_o^{r_1}\right)-I\left(P_o^{r_2}\right)\right)\frac{P_o^{r_1}-P_o^{r_2}}{\Vert P_o^{r_1}-P_o^{r_2}}$$
• 和BRISK一样，计算出主方向！

特征匹配

• 在特征描述中，我们得到了512bit的二进制描述符，该描述符的列是高方差>低方差的排列，而高方差表征了模糊信息，低方差表征了细节信息，与人眼视网膜相似，人眼先处理的是模糊信息，再处理细节信息。
• 选取前128bit即16bytes进行匹配（异或），若两个待匹配的特征点前16bytes距离小于设定的阈值，则再用剩余的位信息进行匹配。这种方法可以剔除掉90%的不相关匹配点。
• 注意：这里的16bytes的选取是建立在并行处理技术（SIMD）上的，并行处理技术处理16bytes与处理1bytes的时间相同；也就是说，16bytes并不是固定的，如果你的并行处理技术能处理32bytes与处理1bytes的时间相同的话，那么你也可以选取前32bytes

特征点过滤

• cv::KeyPointsFilter
• 四个方法：
  • runByImageBorder - 删除图片之外的！
  • runByKeypointSize - 删除关键点尺寸小于minSz，大于maxSz的点
  • runByPixelsMask - 删除掩膜的0区域
  • removeDuplicated - 删除重复
  • retainBest - 选择top-k

特征匹配

• 暴力匹配 BFMatcher
• 快速近似最近邻FlannBasedMatcher

5 几何检测

5.1 基于二值图的几何检测

• Hough变换
• 线检测
• LSD快速直线检测
• 圆检测

5.2 基于非二值图的边缘检测

• 边缘产生的原因：深度不同；所属平面不同；材料导致表面反射光不同；外部场景光不同
• 一阶微分方法，二阶微分方法；混合方法

梯度

• 微分算子Sobel/laplace/Roberts/Prewitt
• Sobel【1 2 1】；Scharr【3 10 3】；laplace【-4】或【-8】
• Prewitt【1 1 1】；Roberts【0 1】或【1 0】
  1）Robert算子对陡峭的低噪声图像效果较好，尤其是边缘正负45度较多的图像，但定位准确率较差；
  2）Prewitt算子对灰度渐变的图像边缘提取效果较好，而没有考虑相邻点的距离远近对当前像素点的影响；
  3）Sobel算子考虑了综合因素，对噪声较多的图像处理效果更好。
  4）Laplacian算子对噪声比较敏感，由于其算法可能会出现双像素边界，常用来判断边缘像素位于图像的明区或暗区，很少用于边缘检测；

canny算子

• 使用高斯滤波器，以平滑图像，滤除噪声。
• 计算图像中每个像素点的梯度强度和方向。Sobel算子
• 应用非极大值（Non-Maximum Suppression）抑制，以消除边缘检测带来的杂散响应。
  • 将当前像素的梯度强度与沿正负梯度方向上的两个像素进行比较。
  • 如果当前像素的梯度强度与另外两个像素相比最大，则该像素点保留为边缘点，否则该像素点将被抑制
    
    
    
• 应用双阈值（Double-Threshold）检测来确定真实的和潜在的边缘。如果边缘像素的梯度值高于高阈值，则将其标记为强边缘像素；如果边缘像素的梯度值小于高阈值并且大于低阈值，则将其标记为弱边缘像素；如果边缘像素的梯度值小于低阈值，则会被抑制。
• 通过抑制孤立的弱边缘最终完成边缘检测。为了跟踪边缘连接，通过查看弱边缘像素及其8个邻域像素，只要其中一个为强边缘像素，则该弱边缘点就可以保留为真实的边缘。

6 形状检测

6.1 基于二值图的形状检测

• 轮廓检测
• 凸包检测
• 边界框
• 距离变换 + 骨架提取

6.2 基于非二值图的区域检测分割

• 最大极值稳定区域MSER提取
• Blob检测器
• 分水岭分割
• FloodFill算法
• MeanShift算法
• Crabcut / Graphcut

7 图像场景识别(分类)

• 词袋模型

8 图像目标检测

• 滑动窗 + 特征工程 + 分类器

9 图片序列(视频)

HDR

• 使用32位浮点数存储 RBG和表示光照亮度变化的A
• 把多图像和曝光值结合，能够全面表示自然场景的高动态范围
• HDR高动态范围的图像生成，需要结合图像的不同曝光率+相机响应曲线 + 压缩算法完成HDR图像生成
• 彩色RGB图像需要对三个通道分别恢复标定相机的响应曲线，计算量大，颜色失真
• Mitsunage等提出使用N次多项式拟合响应曲线，使用曝光度比的估计值来计算多项式系数，然后用多项式重新估计曝光度比，迭代至收敛得到响应函数和精确曝光度比！后序进一步改进！
• Opencv的代码：

void loadExposureSeq(String path, vector<Mat>& images, vector<float>& times)
{
	// 获取当前路径
	path = path + std::string("/");
	// 获取图像相机曝光信息
	ifstream list_file((path + "list.txt").c_str());
	string name;
	float val;
	//  遍历完成图像源信息输入
	while (list_file >> name >> val) {
		Mat img = imread(path + name);
		images.push_back(img);
		times.push_back(1 / val);
	}
	list_file.close();
}
vector<Mat> images;
vector<float> times;
// 获取图像序列
string strPath = "../images/HDR";
loadExposureSeq(strPath, images, times);
// 创建Calibrate对象
Mat response;
Ptr<CalibrateDebevec> calibrate = createCalibrateDebevec();
// 估计相机响应，校准算法来估计逆CRF所有256像素值
calibrate->process(images, response, times);
// 构建HDR图像
Mat hdr;
Ptr<MergeDebevec> merge_debevec = createMergeDebevec();
merge_debevec->process(images, hdr, times, response);
// 图像映射 利用双边滤波 Gamma变换
Mat ldr;
Ptr<TonemapDurand> tonemap = createTonemapDurand(2.2f);
tonemap->process(hdr, ldr);
// 执行曝光融合
Mat fusion;
Ptr<MergeMertens> merge_mertens = createMergeMertens();
merge_mertens->process(images, fusion);
cv::imshow("hdr", hdr);
// 结果输出
imwrite("fusion.png", fusion * 255);
imwrite("ldr.png", ldr * 255);
imwrite("hdr.hdr", hdr);

光流

• 计算机视觉–光流法(optical flow)简介

LK光流

背景提取

目标跟踪

10 立体视觉

相机标定

• 本质矩阵E（Essential Matrix）：反映【空间一点P的像点】在【不同视角摄像机】下【摄像机坐标系】中的表示之间的关系。
• 基础矩阵F（Fundamental Matrix）：反映【空间一点P的像素点】在【不同视角摄像机】下【图像坐标系】中的表示之间的关系
• 双目摄像头内参：焦距 + 扭曲值 + 光点中心

投影

三维视觉