SIFT 三线性差值原理与代码分析

参考了文章:

http://blog.csdn.net/fzthao/article/details/62424271

Jie Pro 在进行特征描述时，讲的非常详细。但未对三线性插值进行阐述。我也是花了好久的时间才慢慢搞懂。有错之处，请指出。

1. 首先需要的几个已知量包括：

     将关键点附近的邻域划分为d*d(Lowe建议d=4)个子区域，每个子区域做为一个种子点，每个种子点有8个方向。每个子区域的大小与关键点方向分配时相同，即每个区域有个子像素。

（1）

（2）

（3）

（4） 

其中（3）式在三线性插值作用巨大，待下面具体解释。 （4）式主要对所有圆形区域内的点进行高斯加权，m为梯度值.

1.     //计算梯度直方图  
2.     for( k = 0; k < len; k++ )  
3.     {  
4.         float rbin = RBin[k];
5. float cbin = CBin[k];  
6. //(Ori[k] - ori)重新调整角度,这个时候算是真正的翻转
7.         float obin = (Ori[k] - ori)*bins_per_rad;  
8.         float mag = Mag[k]*W[k];  
9.   
10. //因为rbin > -1, 所以r0 = cvFloor(rbin)有可能是-1
11.         int r0 = cvFloor( rbin );  
12.         int c0 = cvFloor( cbin );  
13.         int o0 = cvFloor( obin );  
14.         rbin -= r0;  
15.         cbin -= c0;  
16.         obin -= o0;  
17.   
18.         //n为SIFT_DESCR_HIST_BINS：8，即将360°分为8个区间  
19.         if( o0 < 0 )  
20.             o0 += n;  
21.         if( o0 >= n )  
22.             o0 -= n;  
23.           
24.   
25.         // histogram update using tri-linear interpolation  
26.         // 三线性插值  
27.         float v_r1 = mag*rbin, v_r0 = mag - v_r1;  
28.         float v_rc11 = v_r1*cbin, v_rc10 = v_r1 - v_rc11;  
29.         float v_rc01 = v_r0*cbin, v_rc00 = v_r0 - v_rc01;  
30.         float v_rco111 = v_rc11*obin, v_rco110 = v_rc11 - v_rco111;  
31.         float v_rco101 = v_rc10*obin, v_rco100 = v_rc10 - v_rco101;  
32.         float v_rco011 = v_rc01*obin, v_rco010 = v_rc01 - v_rco011;  
33.         float v_rco001 = v_rc00*obin, v_rco000 = v_rc00 - v_rco001;  
34.   
35. //根据以上公式可得:
36. // v_rco000 =(1-rbin)*(1-cbin)*(1-obin)
37. // v_rco001 =(1-rbin)*(1-cbin)*obin
38. 
    
39. //v_rco010 =(1-rbin)*cbin*(1-obin)
40. //v_rco011 =(1-rbin)*cbin*obin
41. 
    
42. //v_rco100 =rbin*(1-cbin)*(1-obin)
43. //v_rco101 =rbin*(1-cbin)*obin
44. 
    
45. // v_rco110 =rbin*cbin*(1-obin)
46. // v_rco111 =rbin*cbin*obin
48. // 因为r0,c0的范围是[-1, d)，所以要先+1
49.         int idx = ((r0+1)*(d+2) + c0+1)*(n+2) + o0;  
50.         hist[idx] += v_rco000;   ==> b1关键点
51.         hist[idx+1] += v_rco001;  
52.         hist[idx+(n+2)] += v_rco010;   ==> b3关键点
53.         hist[idx+(n+3)] += v_rco011;  
54.         hist[idx+(d+2)*(n+2)] += v_rco100; ==> b2关键点  
55.         hist[idx+(d+2)*(n+2)+1] += v_rco101;  
56.         hist[idx+(d+3)*(n+2)] += v_rco110;  ==> b4关键点 
57.         hist[idx+(d+3)*(n+2)+1] += v_rco111;  
58.     }  

1.     for( i = 0; i < d; i++ )  
2.         for( j = 0; j < d; j++ )  
3.         {  
4. //i,j各+1,因为存储的边界点不需要统计，只需要统计绿色框中蓝色端点
5.             int idx = ((i+1)*(d+2) + (j+1))*(n+2);  
6.             hist[idx] += hist[idx+n];  
7.             hist[idx+1] += hist[idx+n+1];  
8.             for( k = 0; k < n; k++ )  
9.                 dst[(i*d + j)*n + k] = hist[idx+k];  
10.         } 

1.     {  
2.         for( j = 0; j < d+2; j++ )  
3.             for( k = 0; k < n+2; k++ )  
4.                 hist[(i*(d+2) + j)*(n+2) + k] = 0.;  
5.     }  
6.   
7.     for( i = -radius, k = 0; i <= radius; i++ )  
8.         for( j = -radius; j <= radius; j++ )  
9.         {  
10.             // Calculate sample's histogram array coords rotated relative to ori.  
11.             // Subtract 0.5 so samples that fall e.g. in the center of row 1 (i.e.  
12.             // r_rot = 1.5) have full weight placed in row 1 after interpolation.  
13.         //根据旋转后的坐标,得到旋转前的坐标
14.  float c_rot = j * cos_t - i * sin_t;  
15.             float r_rot = j * sin_t + i * cos_t;  
16.   
17. //由于cos_t 和 sin_t除过了hist_width，
18. //所以这时旋转后的坐标的 c_rot 和 r_rot 的单位直接就是种子点了。
19. //下面的 rbin 和 cbin 计算的就是种子点的编号。
20. //种子点的编号从(-1,-1)到(d,d)。减去0.5是为了方便差值时计算权重。  
21.   
22. //rbin = r_rot + 1.5
23. //cbin = c_rot + 1.5
24. //调整原点位置,如上图的点(0,0)。
25.             float rbin = r_rot + d/2 - 0.5f;  
26.             float cbin = c_rot + d/2 - 0.5f;  
27.             int r = pt.y + i, c = pt.x + j;  
28.   
29. //确保旋转前的坐标在有效范围内,然后计算梯度幅值和角度。
30. //(-1,d)的范围,其实就是上图的蓝色框框里的所有pixel。
31. //旋转的时候,不会影响梯度幅值,所以可以先计算出来。
32.             if( rbin > -1 && rbin < d && cbin > -1 && cbin < d &&  
33.                 r > 0 && r < rows - 1 && c > 0 && c < cols - 1 )  
34. //超出窗边界或图像边界的点不予考虑。
35. //注意这里 rbin 和 cbin 的下线都是 -1,也就是说它们可能为负  
36.             {  
37.                 float dx = (float)(img.at<sift_wt>(r, c+1) - img.at<sift_wt>(r, c-1));  
38.                 float dy = (float)(img.at<sift_wt>(r-1, c) - img.at<sift_wt>(r+1, c));  
39.                 X[k] = dx; Y[k] = dy; RBin[k] = rbin; CBin[k] = cbin;  
40.                 W[k] = (c_rot * c_rot + r_rot * r_rot)*exp_scale;  
41.                 k++;  
42.             }  
43.         }  
44.   
45.     len = k;  
46. //计算角度
47.     fastAtan2(Y, X, Ori, len, true); 
48. //计算幅值
49.     magnitude(X, Y, Mag, len);  
50.     exp(W, W, len);