新手从零开始，相似图像匹配SIFT算法（三），完结版

时隔半个月，终于可以提笔写这篇从零开始学sift算法的博文了！

本文适合对理论有一定了解了的童鞋，帮你搞定java代码实现。如果不清楚理论，可以结合这篇博文来一起看 http://blog.csdn.net/zddblog/article/details/7521424

经过再三折腾，突然回头一看，发现SIFT并没有想象的那么难，也没有想象的那么强大（这里不指那些改进的sift）！我自己是完全用java语言写的，没有用opencv，或者metlab等工具，虽然过程比较纠结，但对sift的理解也算较为深刻的了！

其实sift的匹配效果受到初始sigma（也就是高斯模糊半径）以及后面的描述子的维度（可以自己增加或减少维度）影响较大，当然各种阈值也是影响匹配结果的！我在上一篇博文里得到的那么差的匹配结果就是因为描述子没有取好。

下面贴一下我的结果（没有剔除错误匹配的点，第一张仅仅是调了匹配的阈值），其实网上贴出来的图片大多都是给别人看的，其实也有蛮多错误匹配点的，只是很多都进行了错误点剔除），自认为也算不错了。在博文最后贴上了java的源代码：

前面的两篇的博文已写了一些参考理论博文或者文章，这里不多扯了！

         这里主要总结一下我从小白到写完sift的全过程需要特别注意的地方，也给自己写代码的童鞋一些参考：

 1、高斯金字塔：高斯金字塔的构建过程中，每一层的高斯模糊核其实都是一样的，如下：

double[] sig=new double[6];sig[0]=baseSigma;for(int i=1;i<gaussS;i++){double preSigma=baseSigma*Math.pow(2, (double)(i-1)/s);double nextSigma=preSigma*Math.pow(2, (double)1/s);sig[i]=Math.sqrt(nextSigma*nextSigma-preSigma*preSigma);}

因为每层的第一张图像是由上一层金字塔的第四张降采样而来的，所以用公式计算出的尺度刚好是上一层的第一张的2倍，即2*baseSigma。

2、计算获取特征点：获取dog图像的极值点并不难实现，但是对于新手来说，后面的去除低对比度和边缘响应点这一步，相信很多新手会焦头烂额，至少我曾经是。

/** * 过滤关键点，得到更加稳定的特征点——————————去除对比度低（方差）和边缘点（hessian矩阵去边缘点） * @param dogPyramid * @param keyPoints * @param r  在Lowe的文章中，取r＝10。图4.2右侧为消除边缘响应后的关键点分布图。 * @param contrast 对比度阈值 * @return */private static HashMap<Integer, List<MyPoint>>  filterPoints(HashMap<Integer,double[][]> dogPyramid,HashMap<Integer, List<MyPoint>> keyPoints,int r,double contrast ){HashMap<Integer, List<MyPoint>> resultMap=new HashMap<Integer, List<MyPoint>>();         //Set<Integer> gaussIndex=gaussPyramid.keySet();Set<Integer> pSet=keyPoints.keySet();//int length=pSet.size()-1;for(int i:pSet){List<MyPoint> points=keyPoints.get(i);List<MyPoint> resultPoints=new ArrayList<MyPoint>();///获取对应的dog图像double[][] gaussImage=dogPyramid.get(i);for(MyPoint mp:points){int x=mp.getX();int y=mp.getY();////对极值点进行精确定位/*mp=adjustLocation(dogPyramid, mp, 6, 6);if(null==mp){///如果精确定位返回null，则抛弃该点continue;}x=mp.getX();y=mp.getY();*/double xy00=gaussImage[y-1][x-1];double xy01=gaussImage[y-1][x];double xy02=gaussImage[y-1][x+1];double xy10=gaussImage[y][x-1];double xy11=gaussImage[y][x];double xy12=gaussImage[y][x+1];double xy20=gaussImage[y+1][x-1];double xy21=gaussImage[y+1][x];double xy22=gaussImage[y+1][x+1];double dxx=xy10+xy12-2*xy11;double dyy=xy01+xy21-2*xy11;double dxy=(xy22-xy20)-(xy02-xy00);///hessian矩阵的对角线值和行列式double trH=dxx+dyy;double detH=dxx*dyy;//-dxy*dxy;///邻域的均值double avg=(xy00+xy01+xy02+xy10+xy11+xy12+xy20+xy21+xy22)/9;///领域方差double DX=(xy00-avg)*(xy00-avg)+(xy01-avg)*(xy01-avg)+(xy02-avg)*(xy02-avg)+(xy10-avg)*(xy10-avg)+(xy11-avg)*(xy11-avg)+(xy12-avg)*(xy12-avg)+(xy20-avg)*(xy20-avg)+(xy21-avg)*(xy21-avg)+(xy22-avg)*(xy22-avg);DX=DX/9;double threshold=(double)(r+1)*(r+1)/r;if((detH>0&&(trH*trH)/detH<threshold)&&(DX>=contrast)){///主曲率小于阈值，则不是需要剔除的边缘响应点;方差大于0.03的为高对比度resultPoints.add(mp);}//System.out.println(DX);}///end of inner forresultMap.put(i, resultPoints);}return resultMap;}

在原文里有一段是精确特征点到亚像素，我在代码里写了，但是没有采用,因为过滤之后反而导致我的结果不好了（或许写的哪里有误，希望看了后交流）。

3、关键点描述：和lowe的原文计算方式一样，统计领域内的梯度方向，并让每个梯度模值高斯核的对应值，相当于加权平均获取模值；但是请注意，关键点的方向（比如0~10°内）并不是取0或者10°，而是计算落在该方向范围内的所有点的方向角度的加权平均数！这一点很重要！也是与原文不一样的。

double s_oct=baseSigma*Math.pow(2, (double)mp.getS()/nLayer);double sigma=1.5*s_oct;////高斯模糊核double[][] gtemplate=getGaussTemplate2D(sigma);////用于对模值进行高斯加权int radius=(int) (3*sigma);///领域采样半径int diameter=2*radius;int gtCenter=gtemplate.length/2;if(x>=diameter&&x<width-diameter&&y>=diameter&&y<height-diameter){///sigma=9for(int j=0;j<=2*radius;j++ ){for(int i=0;i<=2*radius;i++){/*if((j==radius)&&(i==radius)){continue;}*/double ty2=image[y-radius+j+1][x-radius+i];double ty1=image[y-radius+j-1][x-radius+i];double tx2=image[y-radius+j][x-radius+i+1];double tx1=image[y-radius+j][x-radius+i-1];//关梯度模值double tM=Math.sqrt((ty2-ty1)*(ty2-ty1)+(tx2-tx1)*(tx2-tx1));//梯度方向，转为成0~360°之间的角度值double tTheta=Math.atan2(ty2-ty1, tx2-tx1)*(180/Math.PI)+180;int section=(int) (tTheta/9);if(360-Math.abs(tTheta)<0.0001){///如果角度为360°，则和零一样算在第一个一区间内section=0;}keyTM[section]=keyTM[section]+tM*gtemplate[gtCenter-radius+j][gtCenter-radius+i];keyAngle[section]=keyAngle[section]+tTheta*gtemplate[gtCenter-radius+j][gtCenter-radius+i];////按比例对主方向产生角度贡献;angleRatio[section]+=gtemplate[gtCenter-radius+j][gtCenter-radius+i];}}for(int key=0;key<keyTM.length;key++){if(keyTM[key]>max){max=keyTM[key];index=key;}}theta=keyAngle[index]/angleRatio[index];}///对关键如果有多个辅方向，就复制成多个关键点for(int key=0;key<keyTM.length;key++){if(keyTM[key]>max*0.8){///大于最大值得80%都作为主方向之一，复制成一个关键点theta=keyAngle[key]/angleRatio[key];////获得辅方向//System.out.println("theta:"+theta);///计算每个代数圆内的梯度方向分布if(x>=largestDistance+1&&x<width-1-largestDistance&&y>=largestDistance+1&&y<height-1-largestDistance){int secNum=15;//每个代数圆内多少扇形int secAngle=360/secNum;///多大的角为一个扇形double[] grads=new double[secNum*(largestDistance/2)];///保存多维向量的数组double sum=0;for(int j=y-largestDistance;j<=y+largestDistance;j++){for(int i=x-largestDistance;i<=x+largestDistance;i++){if((j==y)&&(i==x)){continue;}d

4、特征sift描述子：原文是采用一个几何原，并计算统计16x16圆内的梯度，最后得到128维特征向量。本文采用了棋盘距离（如图）作为采样半径，半径大小依次为2、4、6、8；经过试验发现最大半径为8效果最好；在每一个小的采样领域内统计梯度分布，这里我选取了15个bin，每个bin  24° 。最后，每个点的方向都是直接减去特征点的方向的，相当于计算一个相对方向，从而得到旋转不变性：

5、特征点之间的特征向量比较，和lowe的理论一样，采用最近距离比上次进距离，如果比值小于一定的阈值则认为匹配上了。这里提一点，如果要在图片上画线，是要把特征点映射到原图上的，否则不用映射回去。

最后附上完整java代码，下载加入到eclipse中即可运行，当然要修改图片的路径！欢迎交流！（PS：代码里有部分数字图像工具类可能 没用到，懒得删除了，就直接打包了）。

 不知道怎么上传附件，直接去csdn下吧，不要积分：http://download.csdn.net/detail/abcd_d_/7107483