图像特征学习笔记（harris\sift\全景图像拼接）更新ing

最新推荐文章于 2023-01-16 18:54:36 发布

cc_sensii

最新推荐文章于 2023-01-16 18:54:36 发布

阅读量222

点赞数

本文链接：https://blog.csdn.net/cc_sensi/article/details/118102293

版权

harris

对于图像 $l\left ( x,y \right )$ ，当在点 $\left ( x,y \right )$ 处平移 $\left ( \Delta x,\Delta y \right )$ 后的自相似性：

$c\left ( x,y;\Delta x,\Delta y \right ) = \sum_{(u,v)\in W(x,y)}^{}w(u,v)(I(u,v)-I(u+\Delta x,v+\Delta y))^{2}$

W(x,y)是以点(x,y)为中心的窗口，即可是常数，也可以是高斯加权函数

基于泰勒展开，对图像 $l\left ( x,y \right )$ 在平移 $\left ( \Delta x,\Delta y \right )$ 后进行一阶相似

$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$

其中，Ix,Iy是I(x,y)的偏导数

近似可得：

$c(x,y;\Delta x,\Delta y) \approx \sum_{w}^{}(I_x(u,v)\Delta x + I_y(u,v)\Delta y)^2 = [\Delta x,\Delta y]M(x,y)\begin{bmatrix} \Delta x\\ \Delta y \end{bmatrix}$

其中M：

$M(x,y) = \sum_{w}^{} \begin{bmatrix} 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{bmatrix} = \begin{bmatrix} A &C \\ C & B \end{bmatrix}$

= $\begin{bmatrix} \lambda _1 & \\ & \lambda _2 \end{bmatrix}$

化简可得：

$c(x,y;\Delta x, \Delta y) \approx A\Delta x^2 + 2C\Delta x\Delta y + B\Delta y^2$

$c = \lambda _{1}\Delta x^2+ \lambda _{2}\Delta y^2$

sift

1 高斯模糊，产生一组图片，同一组图片，大小相等，模糊系数不同

2.降采样，参数为k，每行每列每隔k个点取一个点组，成一幅图像。论文中k=2。

for (o = 0; o < octvs; o++)//金字塔组数为octvs，
	for (i = 0; i < intvls + 3; i++)//每一组有intvls + 3 层，intvls一般为3
		{
		if (o == 0 && i == 0)//如果是第一组第1层
			gauss_pyr[o][i] = cvCloneImage(base);//base 为原始灰度图像经过升采样或降采样得到的图像
	/* base of new octvave is halved image from end of previous octave */
		else if (i == 0)//建立非第一组的第1层
			gauss_pyr[o][i] = downsample(gauss_pyr[o - 1][intvls]);//降采样图像
		/* blur the current octave's last image to create the next one */
		else//建立非第一组的非第1层
			{
			gauss_pyr[o][i] = cvCreateImage(cvGetSize(gauss_pyr[o][i - 1]),IPL_DEPTH_32F, 1);
			cvSmooth(gauss_pyr[o][i - 1], gauss_pyr[o][i],CV_GAUSSIAN, 0, 0, sig[i], sig[i]);// sig[i]为模糊系数
			}//cvSmooth 为平滑处理函数，也即模糊处理。CV_GAUSSIAN 为选用高斯函数对图像模糊
	return gauss_pyr;//返回建好的金字塔

3.高斯差分，生成高斯差分金字塔，1中每两个相邻图相减

for (o = 0; o < octvs; o++)//octvs为高斯金字塔组数
	for (i = 0; i < intvls + 2; i++)//因为相减，故高斯金字塔中每组有（intvls + 2）层图像
		{
		dog_pyr[o][i] = cvCreateImage(cvGetSize(gauss_pyr[o][i]),IPL_DEPTH_32F, 1);
		cvSub(gauss_pyr[o][i + 1], gauss_pyr[o][i], dog_pyr[o][i], NULL);//cvSub为opencv内置相减函数
		}
	return dog_pyr;//返回高斯差分金字塔

4.极值点检测，选取某点与同一张图周围八个点，和相邻图同位置十八个点相比较，若该点最大或最小，则选取该点为极值点

if (val > 0)//极大值检测{
		for (i = -1; i <= 1; i++)
		for (j = -1; j <= 1; j++)
		for (k = -1; k <= 1; k++)
		if (val < pixval32f(dog_pyr[octv][intvl + i], r + j, c + k))//pixval32f为提取图像像素位置上的灰度值
			return 0;}
		else /* check for minimum */
		{
		for (i = -1; i <= 1; i++)
		for (j = -1; j <= 1; j++)
		for (k = -1; k <= 1; k++)
		if (val > pixval32f(dog_pyr[octv][intvl + i], r + j, c + k))//r c为图像的行数和列数，dog_pyr为高斯差分图
			return 0;

5.关键点定位，4中方法检测到的极值点是离散空间的极值点，下面通过拟合三维二次函数来精确确定关键点的位置和尺度

在检测到的极值点 $X_{0}\left ( x_0,y_0,\sigma _0 \right )^T$ 进行三元二阶泰勒展开

对f(X)求导，并令导函数为零，得到X^，带回f(X)，得到极值f(X^)

其中, X^代表相对插值中心的偏移量，当它在任一维度上的偏移量大于0.5时（即x或y或 σ），意味着插值中心已经偏移到它的邻近点上，所以必须改变当前关键点的位置。同时在新的位置上反复插值直到收敛；也有可能超出所设定的迭代次数或者超出图像边界的范围，此时这样的点应该删除，在Lowe论文中进行了5次迭代。另外，过小的点易受噪声的干扰而变得不稳定，所以将小于某个经验值(Lowe论文中使用0.03)的极值点删除。同时，在此过程中获取特征点的精确位置(原位置加上拟合的偏移量)以及尺度(σ)。

6.消除边界响应，Hessian矩阵

一个定义不好的高斯差分算子的极值在横跨边缘的地方有较大的主曲率，而在垂直边缘的方向有较小的主曲率。DOG算子会产生较强的边缘响应，需要剔除不稳定的边缘响应点。

$H = \begin{bmatrix} D_{xx} &D_{xy} \\ D_{xy} & D_{yy} \end{bmatrix}$

Tr(H)为矩阵的迹，Det(H)为矩阵的行列式，令α=rβ则

公式(r+1)^2/r的值在两个特征值相等时最小，随着r的增大而增大。值越大，说明两个特征值的比值越大，即在某一个方向的梯度值越大，而在另一个方向的梯度值越小，而边缘恰恰就是这种情况。所以为了剔除边缘响应点，需要让该比值小于一定的阈值，因此，为了检测主曲率是否在某域值r下，只需检测:

论文中r=10

while (i < SIFT_MAX_INTERP_STEPS)//SIFT_MAX_INTERP_STEPS=5为最大迭代次数，避免长时迭代
		{
		interp_step(dog_pyr, octv, intvl, r, c, &xi, &xr, &xc);// 泰勒展开拟合，xi,xr,xc依次为x、y、σ方向偏移量, 
		if (ABS(xi) < 0.5  &&  ABS(xr) < 0.5  &&  ABS(xc) < 0.5)//如果当前偏移量绝对值中的每个值均小于0.5，退出迭代
			break;
		c += cvRound(xc);//计算行坐标，cvRound 为四舍五入。
		r += cvRound(xr);
		intvl += cvRound(xi);
		if (intvl < 1 ||//不在计算的图像层中
			intvl > intvls ||//高斯差分每组的层数为intvls 
			c < SIFT_IMG_BORDER ||//靠近图像边缘5个像素的区域不做检测，SIFT_IMG_BORDER=5,
			r < SIFT_IMG_BORDER ||
			c >= dog_pyr[octv][0]->width - SIFT_IMG_BORDER ||//靠近图像边缘5个像素的区域不做检测
			r >= dog_pyr[octv][0]->height - SIFT_IMG_BORDER)
			{
			return NULL;
			}
		i++;//迭代计数
		}
static void interp_step(IplImage*** dog_pyr, int octv, int intvl, int r, int c,double* xi, double* xr, double* xc)
	{
	CvMat* dD, *H, *H_inv, X;
	double x[3] = { 0 };
	dD = deriv_3D(dog_pyr, octv, intvl, r, c);//一阶偏导数
	H = hessian_3D(dog_pyr, octv, intvl, r, c);//Hessian 矩阵即二阶导数组成的矩阵
	H_inv = cvCreateMat(3, 3, CV_64FC1);
	cvInvert(H, H_inv, CV_SVD);//求Hessian矩阵的逆矩阵
	cvInitMatHeader(&X, 3, 1, CV_64FC1, x, CV_AUTOSTEP);
	cvGEMM(H_inv, dD, -1, NULL, 0, &X, 0); //cvGEMM为矩阵乘法，//第一个矩阵的系数；//H_inv、dD第一二个矩阵//-1矩阵前的常数//X结果矩阵
	cvReleaseMat(&dD); 
	cvReleaseMat(&H); 
	cvReleaseMat(&H_inv);
	*xi = x[2]; 
	*xr = x[1];
	*xc = x[0];
	}

7.关键点方向匹配

对于在 DOG 金字塔中检测出的关键点，采集其所在高斯金字塔图像3σ领域窗口内像素的梯度和方向分布特征。梯度的模值和方向如下：

保留峰值大于主方向峰值80％的方向作为该关键点的辅方向


#define interp_hist_peak( l, c, r ) ( 0.5 * ((l)-(r)) / ((l) - 2.0*(c) + (r)) )//插值计算式，l为左侧柱子值，r为左侧柱子值
static void add_good_ori_features(CvSeq* features, double* hist, int n,
								  double mag_thr, struct feature* feat)//精确主方向及辅方向
	{
	struct feature* new_feat;
	double bin, PI2 = CV_PI * 2.0;//CV_PI=pi
	int l, r, i;
	for (i = 0; i < n; i++)// 直方图有n=36个小柱子
		{
		l = (i == 0) ? n - 1 : i - 1;//把小柱子看成是循环的，角度的取值为0-360即一个圆周
		r = (i + 1) % n;
		//只对小柱子的值大于等于主峰80%且此小柱子比左右两边小柱子都高的柱子进行抛物线插值
		if (hist[i] > hist[l] && hist[i] > hist[r] && hist[i] >= mag_thr)// mag_thr为>=80%的最高峰值
			{
			bin = i + interp_hist_peak(hist[l], hist[i], hist[r]);//interp_hist_peak 插值函数
			bin = (bin < 0) ? n + bin : (bin >= n) ? bin - n : bin;//角度取值约束在0-360之间，且是连续循环的
			new_feat = clone_feature(feat);//幅值特征点
			new_feat->ori = ((PI2 * bin) / n) - CV_PI;//？
			cvSeqPush(features, new_feat);
			free(new_feat);
			}

8.特征点描述符

已经知道关键点的位置、尺度、方向，为每个关键点建立一个描述符，使其不随各种变化而改变，比如光照变化、视角变化等等。并且描述符应该有较高的独特性，以便于提高特征点正确匹配的概率。

SIFT描述子是关键点邻域高斯图像梯度统计结果的一种表示。通过对关键点周围图像区域分块，计算块内梯度直方图，生成具有独特性的向量，这个向量是该区域图像信息的一种抽象，具有唯一性。

将关键点附近的区域划分为d*d个子区域，每个子区域作为一个种子点，每个种子点有8个方向。论文建议描述子使用在关键点尺度空间内4*4的窗口中计算的8个方向的梯度信息，共4*4*8=128维向量表征。

1.确定描述子所需的图像区域

特征描述子与特征点所在的尺度有关，因此，对梯度的求取应在特征点对应的高斯图像上进行。考虑到实际计算时，需要采用三线性插值，以及旋转因素，实际计算所需的图像区域半径为：

2. 将坐标轴旋转为关键点的方向，以确保旋转不变性。

3. 将邻域内的采样点分配到对应的子区域内，将子区域内的梯度值分配到8个方向上，计算其权值。

旋转后的采样点坐标在半径为radius的圆内被分配到d x d的子区域，计算影响子区域的采样点的梯度和方向 (x' , y') 分配到8个方向上。

旋转后的采样点落在子区域的下标为

$\binom{x''}{y''}= \frac{1}{3\sigma -oct}\binom{x'}{y'}+\frac{d}{2}$

论文建议子区域的像素的梯度大小按的高斯加权计算，即

$w = m(a+x,b+y)*e^{-\frac{(x')^2+(y')^2}{2*(0.5d)^2}}$

其中a，b为关键点在高斯金字塔图像中的位置坐标。

4.插值计算每个种子点八个方向的梯度

将所得采样点在子区域中的下标(x'',y'')(图中蓝色窗口内红色点)线性插值，计算其对每个种子点的贡献。如图中的红色点，落在第0行和第1行之间，对这两行都有贡献。对第0行第3列种子点的贡献因子为dr，对第1行第3列的贡献因子为1-dr，同理，对邻近两列的贡献因子为dc和1-dc，对邻近两个方向的贡献因子为do和1-do。则最终累加在每个方向上的梯度大小为：

其中k，m，n为0或为1。

5.归一化处理

如上统计的4*4*8=128个梯度信息即为该关键点的特征向量。特征向量形成后，为了去除光照变化的影响，需要对它们进行归一化处理，对于图像灰度值整体漂移，图像各点的梯度是邻域像素相减得到，所以也能去除。得到的描述子向量为，归一化后的特征向量为则

6. 描述子向量门限。

非线性光照，相机饱和度变化对造成某些方向的梯度值过大，而对方向的影响微弱。因此设置门限值(向量归一化后，一般取0.2)截断较大的梯度值。然后，再进行一次归一化处理，提高特征的鉴别性。

7. 按特征点的尺度对特征描述向量进行排序。

简单地说，

其实是一个128维的向量，即128个数

区域，分成4×4个子区域

在每个子区域统计八个方向的梯度

统计每个子区域八个方向的梯度长度（通过高斯加权）

把128个数按顺序写出来，就是我们的描述符

描述子三线性插值源码分析

static void interp_hist_entry(double*** hist, double rbin, double cbin,
								  double obin, double mag, int d, int n)
		{
		double d_r, d_c, d_o, v_r, v_c, v_o;
		double** row, *h;
		int r0, c0, o0, rb, cb, ob, r, c, o;
		r0 = cvFloor(rbin);//向下取整
		c0 = cvFloor(cbin);
		o0 = cvFloor(obin);
		d_r = rbin - r0;//小数余项
		d_c = cbin - c0;
		d_o = obin - o0;
		for (r = 0; r <= 1; r++)//沿着行方向不超过1个单位长度
			{
			rb = r0 + r;
			if (rb >= 0 && rb < d)//如果此刻还在真正的描述子区间内
				{
				v_r = mag * ((r == 0) ? 1.0 - d_r : d_r);//d_r = rbin - r0;
				row = hist[rb];
				for (c = 0; c <= 1; c++)//沿着行方向不超过1个单位长度
					{
					cb = c0 + c;
					if (cb >= 0 && cb < d)
						{
						v_c = v_r * ((c == 0) ? 1.0 - d_c : d_c);
						h = row[cb];
						for (o = 0; o <= 1; o++)//沿着直方图方向不超过1个单位角度 
							{ 
							ob = (o0 + o) % n；//n=8，8个小柱子
							v_o = v_c * ((o == 0) ? 1.0 - d_o : d_o);
							h[ob] += v_o;
							}
						}
					}
				}