SIFT算法详解与代码解析

Sift算法详解及代码解析

 

学了SIFT也有1个半月了，真的是坎坷不断，也因为我可能接触图像邻域时间不长，有很多相关知识要学习，直至今日，才把SIFT大致弄明白。但还有些细节值得去深究，我先把我个人对SIFT算法的理解分享给大家。如果有什么错误，欢迎大家指正。

要理解SIFT算法，首先要清楚你要干什么。SIFT的确可以做很多东西，比如说目标识别，图像区域匹配，又或者是三维视觉。但是对于每个不同的用途，需要SIFT算法里面的功能也是不一样的，也就是SIFT可以和别的算法组合使用，要看具体的用途。

SIFT算法从主体上可以分为两个部分：

1.检测特征点。

2.邻域匹配。

我就将从这两点开始往细讲。。。

首先，SIFT和另外的检测特征点方法不同，它是利用了和LOG很像的一种检测方法—DOG，也就是高斯差分金字塔。为的是让程序更加高效的执行，具体的我下面会讲到。而为什么检测特征点要用DOG（differfence of gaussian）呢，其实我们完全可以用HARRI或是FAST算法来代替，本身DOG是很慢的，如果不优化，几乎是不可能使程序达到实时运算的效果，但是DOG的魅力在哪里呢，DOG引入了一个尺度的概念，什么是尺度呢，就比如说：我要匹配图像，两个物体，一个很近，一个又很远，那两张图片就处在了不同的尺度。其实你可以试一试，拍出这样的照片，然后把远的物体放大，你会发现远的物体模糊了，没错这就是尺度，就是一种模糊，这种模糊是模拟人的眼睛看物体的远近，在实际匹配时，如果两个物体在多个尺度都能匹配上，就说明两个物体是同一个。小尺度来检测纹理，大尺度来检测形状，这就是SIFT算法为什么要选用DOG的原因，SIFT是一种物体匹配算法。如果你想做纹理匹配，其实也可以换用HARRIS角点检测，但是DOG检测出来的特征点比HARRIS稳定很多。这也就决定了到时候匹配特征点的数量问题了。

 

 

一、检测特征点

检测特征点我们可以分为三个步骤：

1.建立高斯金字塔

2.建立DOG差分金字塔

3.对特征点进行筛选与精确

 

什么是高斯金字塔？

实际上高斯金字塔就是在模拟人的眼球由远及近。高斯金字塔的顶端是最大的尺度，也就是最远的距离，而底部就是最小尺度，也就是最近距离。

 

怎么体现尺度这个概念？

已被证实，高斯模糊是尺度变换的为一核心，也就是高斯滤波器或者叫高斯模板卷积，我们对一张图像高斯模糊，实际上就是在变大这个图像是尺度。高斯金字塔就是一层一层不断地高斯模糊，不断地加大sigma值，sigma值越大，图像就越模糊。

 

金字塔图样：

 

由图可知，一个金字塔有很多层（oct），而每一层又有很多小层（inte），这些层是怎么来的呢？

首先，我们要对一张图像进行预处理，也就是去除噪声。

利用这张图片作为最底oct的第一个inte，往后每一个inte都是前一个inte进行高斯模糊得到的图像，为了尺度连续，我们让sigma=pow（k，inte），k=pow(2，inte/（G-3）)对每个oct第一张图进行模糊，这样就能保证尺度连续，也就是sigma是连续的。

有人问了为什么这么做，我后面会详细的讲（实际上在做的过程中，我在这个地方卡了很久，按照这个sigma方法求出来的特征点效果并不是很好，后来查看opencv源码发现源码中的sigma是通过前一张图像的sigma与前一张sigma*k的平方和来计算的，并且对前一张图像进行高斯模糊，我的代码示例会有注释 ）而第二层的第一个inte是由上一层的倒数第三张图片降采样得到的。

 

这里有个降采样的概念：降采样就是隔点采样，即对一个图像的奇数列或偶数列，奇数行或偶数行进行提取，得到一张有原图1/4大小的图片。

降采样图片实例：

 

 

降采样代码示例：

CvMat *GDGM(CvMat *bf)
{
	CvMat *af=cvCreateMat(bf->rows/2,bf->cols/2,CV_64FC1);
	int m,n;
	for (int i=0;i<bf->rows-1;i=i+2)
	{
		for (int j=0;j<bf->cols-1;j=j+2)
		{
			m=i/2;
			n=j/2;
			*(af->data.db+m*af->step/8+n)=*(bf->data.db+i*bf->step/8+j);
		}
	}
	return af;


}

高斯公式：

 

高斯模糊：

X-Xi，Y-Yj为周围像素点到中心点的距离。高斯模板为n*n，其中n为奇数，这样就能保证能找到中心点，这样就不用进行双线性插值，并且加权准确性得到保证。将周围每个像素点的n*n的邻域内像素点的灰度值进行高斯加权：将每个邻域像素灰度值乘上相应的高斯权值。边缘像素点的处理方式可以参考opencv的做法。

然后对高斯模板进行归一化，也就是让所有权值加起来为1，保证模糊后的亮度不变性；

 

Lowe建议选择5*5的高斯模板，因为在距离中心点3以外的权值过小，可以忽略，并且加一层带来的时间消耗过大，5*5时平均性能最好。

 计算高斯模板：

CvMat *GussianModel(int row,double stdless)
{
	double x,y;
	double gs;
	int arow=row/2;
	CvMat *gsm=cvCreateMat(row,row,CV_64FC1);
	for (int i=0;i<arow+1;i++)
	{
		for(int j=arow;j<row;j++)
		{
			/*x=j-2.0;
			y=2-i;*/
			gs=(1/(2*3.14*stdless*stdless))*(pow(2.72,-((j-arow)*(j-arow)+(i-arow)*(i-arow))/(2*stdless*stdless)));
			cout<<gs<<" "<<endl;
			*(gsm->data.db+i*gsm->step/8+j)=gs;

		}
	}
	for (int i=0;i<arow+1;i++)
	{
		for (int j=0;j<arow;j++)
		{
			*(gsm->data.db+i*gsm->step/8+j)=*(gsm->data.db+i*gsm->step/8+(4-j));
		}
	}
	for (int i=arow+1;i<row;i++)
	{
		for (int j=0;j<arow;j++)
		{
			*(gsm->data.db+i*gsm->step/8+j)=*(gsm->data.db+j*gsm->step/8+i);
		}
	}
	for (int i=arow;i<row;i++)
	{
		for (int j=arow;j<row;j++)
		{
			*(gsm->data.db+i*gsm->step/8+j)=*(gsm->data.db+(4-i)*gsm->step/8+j);
		}
	}
	cout<<"未归一化之前高斯模板----------------------------------------"<<endl;
	MatOUT(gsm,row,row,1);
	cout<<endl;
	return gsm;

}


CvMat *Gussian2one(CvMat *gsm,int row,int col)
{
	

	double ad=0;
	for (int i=0;i<row;i++)
	{
		for (int j=0;j<col;j++)
		{
			ad+=(*(gsm->data.db+i*gsm->step/8+j));
		}
	}
	for (int i=0;i<row;i++)
	{
		for (int j=0;j<col;j++)
		{
			double k=(*(gsm->data.db+i*gsm->step/8+j));
			(*(gsm->data.db+i*gsm->step/8+j))=k/ad;
		}
	}
	cout<<"归一化之后的高斯模板----------------------------------------"<<endl;
	MatOUT(gsm,row,row,1);
	cout<<endl;
	return gsm;


}

为了更加快速的计算高斯图像，可以将n*n的模板简化为两个n*1的模板，时间复杂度降低了n倍。因为在SIFT算法中，高斯模糊将会占用整个程序的大部分时间，当然也可以直接调用opencv的函数：cvSmooth，参数可以自行百度。

 

有了这些基础知识，那么高斯金字塔也就不难完成了。

 

我们首先要建立一个数据结构，class，struct都没问题。其中包括oct层，每个oct层里面又包括了inte层，每个inte层为一张图像，也就是CvMat *img；

在每个oct层中的inte图像都为第一张inte向后一次模糊产生的。那我们要几个oct层呢？

512*512的图像大概需要6~8层，这样能够保证在各个尺度识别的准确性（其实要是只是做双目视觉的匹配的话，其实1~2层就行，平行摄像头一般不存在两物体过大的尺度差距）。

每个oct层需要建立6个inte层，这个6就是上面所写的inte层间模糊公式的G参数；

为什么要G-3呢，因为差分时候会少首位两张图片，差分后的第一张图片不能计入DOG。

 数据结构示例：

class interval
{
public:
	CvMat *Gussianimg;
};
class octave
{
public:
	interval inter[10];
	
};

class delta
{
public:
	octave oct[20];
	int longoct;
	int longinte;
};

当我们建立好了每个oct的inte时，高斯金字塔也就完成了，恭喜，SIFT的第一步也就完成了，不要着急，DOG就在下面......

 

我们对高斯金字塔中的inte层进行差分，什么叫差分呢？

 

图像矩阵的差分：

两个图像矩阵对应的像素点相减，放进第三张图中，第三张图就是差分出的图像，但是要注意一点：如果减出来是负的，就要归一到0~255之间，比如减出来是-3，那放入第三张图像的灰度值为252。其实在灰度值放入时将值强制转化为unsigned char（UCHAR）就可以啦。

 

当差分完以后，DOG金字塔就形成了，就是下面这样。

 

DOG图像：

 

 

生成DOG：

delta DOG(delta Delta)
{
	int xrow,xcol;
	delta DogDelta;
	for (int i=0;i<Delta.longoct;i++)
	{
		for(int j=0;j<Delta.longinte-1;j++)
		{
			xrow=Delta.oct[i].inter[j].Gussianimg->rows;
			xcol=Delta.oct[i].inter[j].Gussianimg->cols;
			DogDelta.oct[i].inter[j].Gussianimg=cvCreateMat(xrow,xcol,CV_64FC1);
		}
	}

	for (int i=0;i<Delta.longoct;i++)
	{
		for(int j=0;j<Delta.longinte-1;j++)
		{
			CvMat *bf1=Delta.oct[i].inter[j+1].Gussianimg;
			CvMat *bf2=Delta.oct[i].inter[j].Gussianimg;
			CvMat *dst=DogDelta.oct[i].inter[j].Gussianimg;

			for(int i = 0; i<bf1->height;i++)
				for(int j=0;j<bf1->width;j++)
				{
					uchar temp = *(bf1->data.db+i*bf1->step/8+j);
					uchar temp_sec =*(bf2->data.db+i*bf2->step/8+j);
					*(dst->data.db+i*dst->step/8+j)=(uchar)(temp-temp_sec);
				}
			Mat bf=cvarrToMat(Delta.oct[i].inter[j].Gussianimg,true);
			Mat sub=af-bf;
			CvMat tempcopy=sub;
			cvCopy(&tempcopy,DogDelta.oct[i].inter[j].Gussianimg);*/



		}
	}
	DogDelta.longoct=Delta.longoct;
	DogDelta.longinte=Delta.longinte-1;
	return DogDelta;


}

为什么可以用差分来完成对特征点的检测呢？

在lowe论文中有提到

这个式子可以近似的看做Gaussian—Laplace算子。

 

从函数拟合来看 DOG的极值幅度虽然小于Laplace算子，但是这并不影响我们对极值特征点的检测，这些检测出的极值点是初步的，而且有非常大的边缘响应，边缘特征点其实是很不利于匹配的，因为弧度相似在两张图片中可能就会有几个像素的偏差，这样的偏差在做立体视觉会产生很大的误差。所以我们要尽量选择出稳定的特征点，比如角点，还有一些极值非常突出的点。

 

第三步：选取与精确关键点

首先我们要从DOG图像中选取我们的关键点，选取的原则是找到特征极值非常突出的点。怎么找呢？

 

从DOG金字塔的每个oct中的第三个inte层开始，将每个点和它所在inte周围8个邻域点以及上下两个inte层周围的各9个点，一共26个点比较，如果比他们都大或者小，将这个点选取为初步特征点。

 

进一步筛选关键点：将上一步的特征点进行曲率检测，如果曲率过小的点也就是边缘点，剔除。

我们用HESSIAN矩阵完成这个操作

 

Dxx Dxy Dyy分别为二阶求偏导

Dx Dy为一阶求偏导

图像矩阵中的具体操作为：

偏导其实就是微分，微分在大矩阵中体现的就是差分，那么Dx,Dy的操作就是

Dx=F（x+1，y）-F（x-1，y）

Dy=F（x，y+1）-F（x，y-1）

二阶则是

 

 

 

求出HESSIAN矩阵后，求出矩阵的迹和行列式

 

检测Tr^2与Det的比值，令alpha=r*beta

 

最后比较他们之间的值，lowe建议r=10，可以去除大部分曲率响应过小的点。

 

最后精确关键点的位置，这些关键点都是十分稳定的关键点，稳定的意思就是在光照变化及其他环境因素下依然可以检测到的点。拥有十分高的鲁棒性。

精确位置我们要用到泰勒公式去拟合。

我们都知道，极值点的求导等于零，在图像中表现为求x，y偏导都等于零。

 

 

其中，fn(x)表示f(x)的n阶导数，多项式称为函数f(x)在a处的泰勒展开式，剩余的Rn(x)是泰勒公式的余项，是(x-a)^n的高阶无穷小。

我们利用对最后的f（x）求导=0反求deltaX（x的偏移量）来实现对关键点的精确。

实际操作为

 

我们只取前三项就可以实现精确化。

令D（x）=0

 

 

 

 

 

这些公式并不够清晰明了说明怎么列式

列式应该如下：

 

 

 精确特征点：

void FixPoint(CvMat* spoint,int *num,delta dogdelta)
{
	int xoct=dogdelta.longoct;
	int xinte=dogdelta.longinte;
	int i,j;
	int m,n;
	int xrow;
	int xcol;
	int x,y;
	int h=1;
	double dxx,dxy,dyy,dx,dy;
	double extrmx,dex;
	double extrmy,dey;
	double sdey;
	double D_XY;
	double tdey;
	double tr,det;
	double inr=10;      //====曲率因子
	double fininr,fink;
	CvMat *del=cvCreateMat(xoct,10000,CV_64FC1);
	int dele[10];
	int numk=0;
	//=======计算一阶导与二阶导=======
	for (i=0;i<xoct;i++)
	{
		m=2*i;
		CvMat *fixing=cvCreateMat(dogdelta.oct[i].inter[1].Gussianimg->rows,dogdelta.oct[i].inter[1].Gussianimg->cols,CV_64FC1);
		fixing=cvCloneMat(dogdelta.oct[i].inter[1].Gussianimg);
		xrow=fixing->rows;
  		xcol=fixing->cols;
		dele[i]=0;
		
		for (n=0;n<*(num+i);n++)
		{
			x=*(spoint->data.db+m*spoint->step/8+n); 
			y=*(spoint->data.db+(m+1)*spoint->step/8+n);

			dx=(*(fixing->data.db+x*fixing->step/8+y+1)-*(fixing->data.db+x*fixing->step/8+y-1))/(double)(2.0*h);
			dy=(*(fixing->data.db+(x+1)*fixing->step/8+y)-*(fixing->data.db+(x-1)*fixing->step/8+y))/(double)(2.0*h);
			dxx=(*(fixing->data.db+x*fixing->step/8+y+1)+*(fixing->data.db+x*fixing->step/8+y-1)-*(fixing->data.db+x*fixing->step/8+y)*2)/(double)(pow(h,2.0));
			dyy=(*(fixing->data.db+(x+1)*fixing->step/8+y)+*(fixing->data.db+(x-1)*fixing->step/8+y)-*(fixing->data.db+x*fixing->step/8+y)*2)/(double)(pow(h,2.0));
  			dxy=((*(fixing->data.db+(x+1)*fixing->step/8+y+1)+*(fixing->data.db+(x-1)*fixing->step/8+y-1))-(*(fixing->data.db+(x-1)*fixing->step/8+y+1)+*(fixing->data.db+(x+1)*fixing->step/8+y-1)))/(double)(4.0*pow(h,2.0));

			//==========消除曲率响应=============
			tr=dxx+dyy;
			det=dxx*dyy-dxy*dxy;
			fininr=(inr+1)*(inr+1)/inr;
			fink=tr*tr/det;
			if (fink>=fininr)
			{
				*(del->data.db+i*del->step/8+dele[i])=n;
				dele[i]++;
			}
			//避免做检测时有特征点邻域溢出====
			int l=50;
			if (x>xrow-l*1.4)
			{
				*(del->data.db+i*del->step/8+dele[i])=n;
				dele[i]++;
			}

			if (x<l*1.4)
			{
				*(del->data.db+i*del->step/8+dele[i])=n;
				dele[i]++;
			}

			if (y>xcol-l*1.4)
			{
				*(del->data.db+i*del->step/8+dele[i])=n;
				dele[i]++;
			}

			if (y<l*1.4)
			{
				*(del->data.db+i*del->step/8+dele[i])=n;
				dele[i]++;
			}
			//===========================

			if (dxy!=0) //&& dx!=0 && dy!=0)
			{
				dex=(dy-dx*dyy/dxy)/(dxx*dyy/dxy-dxy);
				dey=(-dx-dex*dxx)/(dxy);
				sdey=(-dy-dex*dxy)/(dyy);
				tdey=(dy-dx*dxy/dxx)/(dxy*dxy/dxx-dyy);
				D_XY=dex*dx+dey*dy+0.5*dex*dex*dxx+0.5*dey*dey*dyy+dex*dey*dxy;


				if (abs(D_XY)<0.03)
				{
					*(del->data.db+i*del->step/8+dele[i])=n;
					dele[i]++;
				}
	
				*(spoint->data.db+m*spoint->step/8+n)=*(spoint->data.db+m*spoint->step/8+n)+dex;
				*(spoint->data.db+(m+1)*spoint->step/8+n)=*(spoint->data.db+(m+1)*spoint->step/8+n)+dey;
				



			}
			else
			{
	/*			*(del->data.db+i*del->step/8+dele[i])=n;
				dele[i]++;
*/
			}
		}
		

	}


	DelPoint(spoint,del,num,dele);//删除不要的特征点
	

	//for (int o=0;o<xoct;o++)
	//{
	//	MatOUT(spoint,2,*(num+o),1);
	//}




}

现在已经找到了足够多稳定的特征点了，这样SIFT的第一大步就完成了！！！！！！！

 

 

 

第二大步就是如何将两幅图的特征点匹配起来！

 

怎么去让两幅图中的特征点匹配起来呢，lowe想到了用邻域匹配的方法。

我们知道每个特征点周围的纹理都是不一样的，也就是特征点周围的每个点的梯度方向都是独特的。

梯度方向：灰度值下降最快的方向，一般会垂直于边缘。

那我们其实只需要将每个特征点周围邻域的剃度值描述出来就可以了。

描述时我们要借助梯度直方图：就是将梯度（角度值）作为横轴，而梯度模值作为纵轴。

 

Lowe用的是128维向量进行描述，在特征点周围选用长为16*16的正方形，当然也可以根据尺度的实际情况来选取长度，尺度越大，长度应该越长。

我们将这16*16的邻域分成4*4的区块，一个区块就称为一个种子点，每个种子点内又包含了4*4的像素点，我们对每个种子点都做一个8维直方图，也就是45°为一个直方格，所以最后得到的向量就是4*4*8=128维。

但是我们不能一开始就先描述，因为如果要检测的物体发生了旋转，那么他周围的纹理也会发生旋转，并且是同步的，所以我们要找到这个旋转的角度，让统计128维直方图时候也根据这个角度做旋转，这样不管怎么旋转，特征点周围的纹理描述出来都是一致的。这就是旋转不变性！

梯度与角度公式为：

 

这个角度其实就是特征点的梯度，但是为了剔除其他因素的影响并增强鲁棒性，我们选取特征点周围8*8的正方形邻域做统计，做出36维梯度直方图，也就是10°为一个直方格。

 

每个直方格的幅值就是邻域内所有在这个角度的特征点的梯度模值之和，距离特征点越远的邻域，所加到幅值的权值越小，可以利用高斯加权，但权值不能过大，大概在0.2~0.4即可。

最后统计出的直方图为了防止突变，lowe建议再对直方图高斯滤波。

最后，幅值最高的那个角度就是该特征点的主方向（MAIN DIRECTION）。

为了增强鲁棒性，我们还将最大幅值80%以上的角度作为辅助方向，即复制这个特征点，赋值上辅助方向。

这样，每个特征描述向量所需要旋转的角度就确定下来了。

生成MAIN_DIRECTION:

CvMat *hi(delta Delta,CvMat *spoint,int *num)
{
	//循环变量===
	int i,j,l,k,p;
	//========

	int xoct=Delta.longoct;//delta的层数

	CvMat *main_direction=cvCreateMat(xoct,100000,CV_64FC1);//=========================存放每个特征点的主方向 序号对应为spoint的序号


	//创建圆形邻域===
	int neibor[2][300];//======================储存每个邻域点的x y 调用时需要用中心点加上这些值
	int neibor_num=0;//======================储存邻域内一共有多少点
	for (i=0;i<17;i++)//======================8*2+1个像素点为长和宽的正方形里面寻找半径为8的圆
	{
		for (j=0;j<17;j++)
		{
			int nx=-8+i;
			int ny=-8+j;
			int dx=abs(nx);
			int dy=abs(ny);
			double modxy=sqrt(dx*dx+dy*dy);
			if (modxy<=8)
			{
				neibor[0][neibor_num]=nx;
				neibor[1][neibor_num]=ny;
				neibor_num++;
     		}

		}
	}
	//===========



	//对圆形邻域内的点进行分块 分成36块===
	int neiborbin[36][100];//======================储存每个扇形柱内点的序号，对应的是neibor[2][300]，即要寻找坐标，用序号去在neibor二维数组用对应的序号查找
	int neiborbin_num[36];//====================储存每个扇形柱内点的数量
	for (i=0;i<36;i++)
	{
		neiborbin_num[i]=0;
	}

	for (i=0;i<neibor_num;i++)
	{
		int nx=neibor[0][i];
		int ny=neibor[1][i];
		double k=(double)(ny)/(double)(nx);
		int sitabin;

		if ((ny==0)&&(nx==0))//====================中心点不计入
		{
			/*sitabin=0;
			neiborbin[sitabin][neiborbin_num[sitabin]]=i;*/
		}
		else if(ny==0)
		{
			if (nx>0)
			{
				sitabin=0;
				neiborbin[sitabin][neiborbin_num[sitabin]]=i;
				neiborbin_num[sitabin]++;
			} 
			else
			{
				sitabin=180/10;
				neiborbin[sitabin][neiborbin_num[sitabin]]=i;
				neiborbin_num[sitabin]++;
			}

		}
		else if(nx==0)
		{
			if (ny>0)
			{
				sitabin=90/10;
				neiborbin[sitabin][neiborbin_num[sitabin]]=i;
				neiborbin_num[sitabin]++;
			} 
			else
			{
				sitabin=270/10;
				neiborbin[sitabin][neiborbin_num[sitabin]]=i;
				neiborbin_num[sitabin]++;
			}

		}
		else
		{
			if (k>0)
			{
				if (nx>0)
				{
					double sita=atan(k)*180/3.14159;
					if (sita<0)
					{
						sita=sita+360;
					}
					sitabin=(int)(sita/10);
					neiborbin[sitabin][neiborbin_num[sitabin]]=i;
					neiborbin_num[sitabin]++;
				}
				else
				{
					double sita=atan(k)*180/3.14159+180;
					sitabin=(int)(sita/10);
					neiborbin[sitabin][neiborbin_num[sitabin]]=i;
					neiborbin_num[sitabin]++;
				}
			}

			if (k<0)
			{
				if (nx>0)
				{
					double sita=atan(k)*180/3.14159;
					if (sita<0)
					{
						sita=sita+360;
					}
					sitabin=(int)(sita/10);
					neiborbin[sitabin][neiborbin_num[sitabin]]=i;
					neiborbin_num[sitabin]++;
				} 
				else
				{
					double sita=atan(k)*180/3.14159+180;
					sitabin=(int)(sita/10);
					neiborbin[sitabin][neiborbin_num[sitabin]]=i;
					neiborbin_num[sitabin]++;
				}
			}

		}


	}

	for(i=0;i<36;i++)
	{
		int x=neiborbin_num[i];
		cout<<x<<" ";
	}

	//==================================================

	for (i=0;i<xoct;i++)
	{

		CvMat *NowMat=Delta.oct[i].inter[0].Gussianimg;//===========================需要处理的图像矩阵
		int point_num=*(num+i);
		for(j=0;j<point_num;j++)
		{
			double doux=*(spoint->data.db+(2*i)*spoint->step/8+j);
			double douy=*(spoint->data.db+(2*i+1)*spoint->step/8+j);
			int x=cvRound(doux);//======================================特征点的坐标x   整形坐标
			int y=cvRound(douy);//======================================特征点的坐标y   整形坐标


			//=======================================计算主方向与80%的辅助方向添加

			double direction[36];//==========================存放角度与模值
			for (l=0;l<36;l++)//============================初始化direction
			{
				direction[l]=0;
			}


			for (l=0;l<neibor_num;l++)
			{
				int dx=neibor[0][l];//==========================邻域中相对主坐标的相对坐标
				int dy=neibor[1][l];

				int nx=x+dx;//=============================围绕主坐标的邻域坐标
				int ny=y+dy;

				double mod;//=============================邻域坐标的模值
				double sita;//==============================邻域坐标的角度梯度

				mod=Mod(NowMat,nx,ny);//gussian_value(dx,dy,0.5);
				sita=Grad_sita(NowMat,nx,ny);

				int sitabin;//==============================存放每个角度在direction里面的位置
				sitabin=(int)sita/10;
				
				direction[sitabin]+=mod;
			}

			int topbin=0;//================================存放最大值的角度柱
			double topmod=direction[0];//=============================存放最大值的总模值
			int secbin=0;//================================存放辅助值的角度柱
			double secmod=direction[0];//=============================存放辅助值的总模值

			for (l=0;l<36;l++)//=============================找主方向
			{
				if (direction[l]>topmod)
				{
					topmod=direction[l];
					topbin=l;
				}

			}

			//for (l=0;l<36;l++)//==============================找辅助方向
			//{
			//	if (l==topbin)
			//	{
			//	}
			//	else if (direction[l]>secmod)
			//	{
			//		secmod=direction[l];
			//		secbin=l;
			//	}

			//}
		
			//if ((double)(secmod)>=(double)(0.8*topmod))
			//{
			//	*(spoint->data.db+(2*i)*spoint->step/8+*(num+i))=doux;
			//	*(spoint->data.db+(2*i+1)*spoint->step/8+*(num+i))=douy;
			//	*(main_direction->data.db+i*main_direction->step/8+*(num+i))=secbin*10;
			//	*(num+i)=*(num+i)+1;
			//	
			//}


			*(main_direction->data.db+i*main_direction->step/8+j)=topbin*10;


			//==============================================================主方向梯度计算完成

		}
	}







	return main_direction;
}
	

大概效果就是下面这样：

 

 

 

下一步我们就可以进行描述了：

先找出每个特征点的邻域坐标，邻域需要旋转，16*16的正方形旋转至45°时就会出现最大边界16*sqrt（2），再加上我们的XY轴的一个像素，最后我们需要确定的特征邻域为16*sqrt（2）+1。

然后我们要对这个大正方形的所有像素做旋转，求出旋转以后原像素点在新坐标轴下的坐标。旋转公式为：

 

 

 

得出的坐标就是原像素点在新坐标系下的坐标。

我们需要利用这些坐标上的像素灰度值来求出新坐标系中每个种子点的8维特征向量。

因为这些带有灰度信息的点并不是均匀分布在每个种子点里面，所以我们要对所有带有灰度信息的点进行三线性插值。

为什么说是三线性插值，因为每个种子点对x，y，sita都存在不均匀分布，所以我们首先要对这三个变量进行插值。

 

线性插值：

 

首先度Q12，Q22线性插值求出R2，就是按比例求出R2的灰度值，同理求出R1。这是单线性插值。利用R1，R2求出p这就是双线性插值，三线性插值也一样。

这里给出双线性插值的公式：

f(x,y)=f(0,0)(1-x)(1-y)+f(0,1)(1-x)y+f(1,1)xy+f(1,0)x(1-y)

最后在加到相应的直方块时，为了降低距离远的邻域点可能发生不可预测事件的几率，我们利用高斯权值来降低远邻域点的权重，lowe建议sigma选用0.5。

生成的向量如下：

 

描述特征点：

for (i=0;i<xoct;i++)
{

	CvMat *NowMat=Delta.oct[i].inter[0].Gussianimg;//===========================需要处理的图像矩阵
	IplImage *ax=cvCreateImage(cvGetSize(NowMat),IPL_DEPTH_8U,1);
	cvConvert(NowMat,ax);
	int point_num=*(num+i);
	for(j=0;j<point_num;j++)
	{

		/*Mat Matnow=cvarrToMat(ax,true);*/
		float ft[128];
		double x=*(spoint->data.db+2*i*spoint->step/8+j);
		double y=*(spoint->data.db+(2*i+1)*spoint->step/8+j);
		Point2f p(y,x);
		float direct=*(main_direction->data.db+i*main_direction->step/8+j);
		
		//cvNamedWindow("a");
		//cvShowImage("a",ax);
		//cvWaitKey(0);
		/*namedWindow("as");
		imshow("as",Matnow);
		waitKey(0);*/
		calcSIFTDescriptor(NowMat,p,direct,ft);
		for (k=0;k<128;k++)
		{
			CvMat *z=F->oct_imgFeature[i].p_Feature[j].Feature_Vector;
			*(z->data.db+k)=ft[k];
		}
		

	}
}

void calcSIFTDescriptor( CvMat* img, Point2f ptf, float ori, float* dst )
						   
{
	float SIFT_DESCR_MAG_THR=0.2;
	float SIFT_INT_DESCR_FCTR=512.0;
	int d=4;
	int n=8;
    Point pt(cvRound(ptf.x), cvRound(ptf.y));
	//计算余弦，正弦，CV_PI/180:将角度值转化为幅度值
    float cos_t = cosf(ori*(float)(CV_PI/180));
    float sin_t = sinf(ori*(float)(CV_PI/180));
    float bins_per_rad = n / 360.f;
    float exp_scale = -1.f/(d * d * 0.5f);
    float hist_width = 4.0;  
	                                      
	// 计算图像区域半径mσ(d+1)/2*sqrt(2)
	
    int radius = cvRound(hist_width * 1.4142135623730951f * (d + 1) * 0.5f);
    cos_t /= hist_width;
    sin_t /= hist_width;

    int i, j, k, len = (radius*2+1)*(radius*2+1), histlen = (d+2)*(d+2)*(n+2);
    int rows = img->rows, cols = img->cols;

    AutoBuffer<float> buf(len*6 + histlen);
    float *X = buf, *Y = X + len, *Mag = Y, *Ori = Mag + len, *W = Ori + len;
    float *RBin = W + len, *CBin = RBin + len, *hist = CBin + len;

	//初始化直方图
    for( i = 0; i < d+2; i++ )
    {
        for( j = 0; j < d+2; j++ )
            for( k = 0; k < n+2; k++ )
                hist[(i*(d+2) + j)*(n+2) + k] = 0.;
    }

	//计算采样区域点坐标旋转
    for( i = -radius, k = 0; i <= radius; i++ )
        for( j = -radius; j <= radius; j++ )
        {
           
            float c_rot = j * cos_t - i * sin_t;
            float r_rot = j * sin_t + i * cos_t;
            float rbin = r_rot + d/2 - 0.5f;
            float cbin = c_rot + d/2 - 0.5f;
            int r = pt.y + i, c = pt.x + j;

            if( rbin > -1 && rbin < d && cbin > -1 && cbin < d &&
               r > 0 && r < rows - 1 && c > 0 && c < cols - 1 )
            {
                
				float dx=(float)(*(img->data.db+r*img->step/8+c+1)-*(img->data.db+r*img->step/8+c-1));
				float dy=(float)(*(img->data.db+(r-1)*img->step/8+c)-*(img->data.db+(r+1)*img->step/8+c));
                X[k] = dx; Y[k] = dy; RBin[k] = rbin; CBin[k] = cbin;
                W[k] = (c_rot * c_rot + r_rot * r_rot)*exp_scale;
                k++;
            }
        }

	len = k;
	for (int i=0;i<len;i++)
	{
		
		Ori[i]=fastAtan2(Y[i], X[i]);
		Mag[i]=sqrt(X[i]*X[i]+Y[i]*Y[i]) ;
		W[i]=exp(W[i]);

	}


 
	
	//计算梯度直方图
    for( k = 0; k < len; k++ )
    {
        float rbin = RBin[k], cbin = CBin[k];
        float obin = (Ori[k] - ori)*bins_per_rad;
        float mag = Mag[k]*W[k];

        int r0 = cvFloor( rbin );
        int c0 = cvFloor( cbin );
        int o0 = cvFloor( obin );
        rbin -= r0;
        cbin -= c0;
        obin -= o0;

        //n为SIFT_DESCR_HIST_BINS：8，即将360°分为8个区间
		if( o0 < 0 )
            o0 += n;
        if( o0 >= n )
            o0 -= n;
		

        // histogram update using tri-linear interpolation
		// 双线性插值
        float v_r1 = mag*rbin, v_r0 = mag - v_r1;
        float v_rc11 = v_r1*cbin, v_rc10 = v_r1 - v_rc11;
        float v_rc01 = v_r0*cbin, v_rc00 = v_r0 - v_rc01;
        float v_rco111 = v_rc11*obin, v_rco110 = v_rc11 - v_rco111;
        float v_rco101 = v_rc10*obin, v_rco100 = v_rc10 - v_rco101;
        float v_rco011 = v_rc01*obin, v_rco010 = v_rc01 - v_rco011;
        float v_rco001 = v_rc00*obin, v_rco000 = v_rc00 - v_rco001;

        int idx = ((r0+1)*(d+2) + c0+1)*(n+2) + o0;
        hist[idx] += v_rco000;
        hist[idx+1] += v_rco001;
        hist[idx+(n+2)] += v_rco010;
        hist[idx+(n+3)] += v_rco011;
        hist[idx+(d+2)*(n+2)] += v_rco100;
        hist[idx+(d+2)*(n+2)+1] += v_rco101;
        hist[idx+(d+3)*(n+2)] += v_rco110;
        hist[idx+(d+3)*(n+2)+1] += v_rco111;
    }
	
    for( i = 0; i < d; i++ )
        for( j = 0; j < d; j++ )
        {
            int idx = ((i+1)*(d+2) + (j+1))*(n+2);
            hist[idx] += hist[idx+n];
            hist[idx+1] += hist[idx+n+1];
            for( k = 0; k < n; k++ )
                dst[(i*d + j)*n + k] = hist[idx+k];
        }

    float nrm2 = 0;
    len = d*d*n;
    for( k = 0; k < len; k++ )
        nrm2 += dst[k]*dst[k];
    float thr = sqrt(nrm2)*SIFT_DESCR_MAG_THR;
    for( i = 0, nrm2 = 0; i < k; i++ )
    {
        float val = min(dst[i], thr);
        dst[i] = val;
        nrm2 += val*val;
    }
    nrm2 = SIFT_INT_DESCR_FCTR/max(sqrt(nrm2), FLT_EPSILON);
    for( k = 0; k < len; k++ )
    {
        dst[k] = saturate_cast<uchar>(dst[k]*nrm2);
    }
}

到此，SIFT算法的生成就算完了，剩下的就是怎么把两个向量给匹配起来。

 

我们用到的概念是二范数，也就是欧几里得距离。

欧几里德距离：向量的自然长度（即该点到原点的距离）。

匹配其实就是找出欧式距离最近的，也就是匹配度最高的向量；

O_Distance=sqrt(pow(m11-m21,2)+pow(m12-m22,2)+pow(m13-m23,2)+..........+pow(m1128-m2128,2)));

 

计算欧式距离：

CvMat *O_dis_fix(Feature *F1,Feature *F2,int F1_oct,int F2_oct)
{
	int i,j,k,l;
	int F1_num=F1->oct_imgFeature[F1_oct].Feature_num;
	int F2_num=F2->oct_imgFeature[F2_oct].Feature_num;

	//O_dist存放计算好的欧式距离 行为F1的每个特征点 列为每个F1特征点对于与F2特征点的欧式距离
	CvMat *O_dist=cvCreateMat(F1_num,F2_num,CV_64FC1);
	double O=0;
	double temp=0;

	for (i=0;i<F1_num;i++)
	{
		for (j=0;j<F2_num;j++)
		{
			temp=0;
			CvMat *p1=F1->oct_imgFeature[F1_oct].p_Feature[i].Feature_Vector;
			CvMat *p2=F2->oct_imgFeature[F2_oct].p_Feature[j].Feature_Vector;

			for (k=0;k<128;k++)
			{
				double lp1=*(p1->data.db+k);
				double lp2=*(p2->data.db+k);
				temp+=pow(lp1-lp2,2.0);

			}

			O=sqrt(temp);
			*(O_dist->data.db+i*O_dist->step/8+j)=O;

		}
	}
	return O_dist;


}

但是因为每张图片的特殊性，我们不能将欧式距离的大小作为是否匹配成功的条件，Lowe提出了一种ratio比值限制法的解决方案。

我们把ratio比值定义为：次短欧式距离/最短欧式距离；

通常误匹配点的ratio会很大，因为误匹配点的最短距离一般会混淆于其他欧式距离，通俗的讲就是欧式距离短的不够明显。而正确的匹配点欧式距离会很突出，正好可以匹配上。这样我们就可以在不同的图像中匹配到最稳定的特征点了。

Lowe建议ratio取0.8，实际取值在0.6左右即可，如果是做三维视觉的匹配最好还是在0.4，这样能保证精确度，其次也可以加入一些限制条件：比如平行视觉就可以加入特征点平行的条件，或者是两幅图互匹配结果相同防止多个点对应于一个点的条件等。

至此，SIFT就已经全部完成。

欧式匹配法：

for (i=0;i<F1->oct_imgFeature[oct_F1].Feature_num;i++)
{
	comp1=0;
	comp2=0;
	for (j=0;j<F2->oct_imgFeature[oct_F2].Feature_num;j++)
	{
		odis1=*(dist->data.db+i*dist->step/8+j);
		if (j==0)
		{
			comp1=odis1;
			reco1=0;
		}
		else if(odis1<comp1)
		{
			comp1=odis1;
			reco1=j;

		}
	}
	for (j=0;j<F2->oct_imgFeature[oct_F2].Feature_num;j++)
	{
		odis2=*(dist->data.db+i*dist->step/8+j);
		if (j==0)
		{
			comp2=odis2;
			reco2=0;
		}
		else if(j==reco1)
		{
		}
		else if(odis2<comp2)
		{
			comp2=odis2;
			reco2=j;

		}
	}
		
	//if (comp<0.05)
	//{
	//	*(point_pos1->data.db+0*point_pos1->step/8+numx1)=i;
	//	*(point_pos1->data.db+1*point_pos1->step/8+numx1)=reco;
	//	numx1++;
	//}

	double ratio=comp1/comp2;
	if (ratio<_ratio)
	{
		/*if (comp1<)
		{
		}*/
		*(point_pos->data.db+0*point_pos->step/8+numx)=i;
		*(point_pos->data.db+1*point_pos->step/8+numx)=reco1;
		numx++;
	}
	//*(point_pos->data.db+0*point_pos->step/8+numx)=i;
	//*(point_pos->data.db+1*point_pos->step/8+numx)=reco1;
	//numx++;
		
		
}

自己做的SIFT程序匹配样张：

感谢David G.Lowe的SIFT算法论文

参考资料：

1.百度百科 http://baike.baidu.com/link?url=mveWSNZgmkbVESc2XOJviBHLapQ6wtE4SpOsvoz9bhwN3dA6mpQ0OPvkYHcdm0B7R4WsTmQ944J4-BiAQWiDLq

2.z大博文 http://blog.csdn.net/zddblog/article/details/7521424

3.ROB HESS维护SIFT库 http://blogs.oregonstate.edu/hess/code/sift/

4.R大博文 http://blog.csdn.net/abcjennifer/article/details/7639681

5.源码解析 http://blog.csdn.net/xiaowei_cqu/article/details/8113565