新手從零開始，相似影象匹配SIFT演算法（三），完結版
時隔半個月，終於可以提筆寫這篇從零開始學sift演算法的博文了！
經過再三折騰，突然回頭一看，發現SIFT並沒有想象的那麼難，也沒有想象的那麼強大（這裡不指那些改進的sift）！我自己是完全用java語言寫的，沒有用opencv，或者metlab等工具，雖然過程比較糾結，但對sift的理解也算較為深刻的了！
其實sift的匹配效果受到初始sigma（也就是高斯模糊半徑）以及後面的描述子的維度（可以自己增加或減少維度）影響較大，當然各種閾值也是影響匹配結果的！我在上一篇博文裡得到的那麼差的匹配結果就是因為描述子沒有取好。
下面貼一下我的結果（沒有剔除錯誤匹配的點，第一張僅僅是調了匹配的閾值），其實網上貼出來的圖片大多都是給別人看的，其實也有蠻多錯誤匹配點的，只是很多都進行了錯誤點剔除），自認為也算不錯了。在博文最後貼上了java的原始碼：
前面的兩篇的博文已寫了一些參考理論博文或者文章，這裡不多扯了！
這裡主要總結一下我從小白到寫完sift的全過程需要特別注意的地方，也給自己寫程式碼的童鞋一些參考：
1、高斯金字塔：高斯金字塔的構建過程中，每一層的高斯模糊核其實都是一樣的，如下：
因為每層的第一張影象是由上一層金字塔的第四張降取樣而來的，所以用公式計算出的尺度剛好是上一層的第一張的2倍，即2*baseSigma。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、計算獲取特徵點：獲取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 for resultMap.put(i, resultPoints); } return resultMap; }
在原文裡有一段是精確特徵點到亞畫素，我在程式碼裡寫了，但是沒有采用,因為過濾之後反而導致我的結果不好了（或許寫的哪裡有誤，希望看了後交流）。
3、關鍵點描述：和lowe的原文計算方式一樣，統計
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=9
for(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;
}
d4、特徵sift描述子：原文是採用一個幾何原，並計算統計16x16圓內的梯度，最後得到128維特徵向量。本文采用了棋盤距離（如圖）作為取樣半徑，半徑大小依次為2、4、6、8；經過試驗發現最大半徑為8效果最好；在每一個小的取樣領域內統計梯度分佈，這裡我選取了15個bin，每個bin 24° 。最後，每個點的方向都是直接減去特徵點的方向的，相當於計算一個相對方向，從而得到旋轉不變性：
5、特徵點之間的特徵向量比較，和lowe的理論一樣，採用最近距離比上次進距離，如果比值小於一定的閾值則認為匹配上了。這裡提一點，如果要在圖片上畫線，是要把特徵點對映到原圖上的，否則不用映射回去。
最後附上完整java程式碼，下載加入到eclipse中即可執行，當然要修改圖片的路徑！歡迎交流！（PS：程式碼裡有部分數字影象工具類可能 沒用到，懶得刪除了，就直接打包了）。