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

来源：互联网发布：sigmaplot linux 编辑：程序博客网时间：2024/06/06 03:12
参考了文章:
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-3）式可推知，的取值范围为[-d/2,d/2]，为方便计算三线性插值下标从0开始，公式(3) 加了d/2，使(x”,y")的范围为[0,d]，而d=4. 因此在水平X轴和垂直Y轴的单位为1的间隔下，共有16个种子点区域（见左上图的绿色格子），每个区域的中心点分别为，bin(0,0),b(0,1),...,b(3,3)（见左上图的蓝色格子在每个绿色格子的中心交叉点为bin的坐标）。根据(x”,y")的值确定其最邻近的4个中心点。
   
为什么直方图长度histlen = (d+2)*(d+2)*(n+2)?
因为计算的是每个Pixel对每一个蓝色端点(包含绿色格子的中心点)的影响，所以每一行一共6个端点，每一个列也有6个端点。也就是d+2个点。

如左上图，若(x”,y")距离X轴和Y轴的距离分别为（dr，dc），则另外的几个点分别为 （dr，1-dc）,（1-dr，dc）,（1-dr，1-dc）；

相邻4个中心点的距离：
b1中心点的距离为：(dr，dc)
b2中心点的距离为：(1-dr，dc)
b3中心点的距离为：(dr，1-dc)
b4中心点的距离为：(1-dr，1-dc)

根据权重与距离成反比的原则计算邻近的四个种子点分配的到权重：
即：用该点对角线方向的矩形面积为该点的权重值
则分配到
               b1的值为： W * (1-dr)*(1-dc)
b2的值为： W * dr*(1-dc)
               b3的值为： W * (1-dr)*dc
b4的值为： W * dr*dc  
考虑到方向，见右上图。
角度转化为比例分量，为do，那么距离其邻近方向的距离为(1-do)。
则点(x”,y") 的W分配的权重如下：

方向2与方向1相邻，比如方向1是[0° 45°) 方向2就是[45° 90°)

分配给b1在方向1上的权重为：  W * (1-dr)*(1-dc)*(1-do) ,  方向2上的权重为：W * (1-dr)*(1-dc) *do
分配给b2在方向1上的权重为：  W * dr*(1-dc) *(1-do) ,       方向2上的权重为：W * dr*(1-dc) *do
分配给b3在方向1上的权重为：  W * (1-dr)*dc *(1-do) ,       方向2上的权重为：W * (1-dr)*dc *do
分配给b4在方向1上的权重为：  W * dr*dc*(1-do) , 方向2上的权重为：W*dr*dc*do

以上为对三线性插值的分步解释。也可以使用公式


进行直接计算。k,m,n 为0或1.

2. 为什么进行三线性插值？
答：为了减少一个梯度幅值从一个格子漂移(shift)到另一个格子引起的描述子突变，需要对梯度值做三线性插值。也就是根据三维坐标计算距离周围格子的距离，按距离的倒数计算权重，将梯度幅值按权重分配到临近的格子里。

3. 附zddhub的一幅图
 
       * 注意： 在(3)式进行计算(x'',y'')时，会减去0.5, 主要是为了平衡之前计算radius时多加的0.5.


    //计算梯度直方图  
    for( k = 0; k < len; k++ )  
    {  
        float rbin = RBin[k];
float cbin = CBin[k];  
//(Ori[k] - ori)重新调整角度,这个时候算是真正的翻转
        float obin = (Ori[k] - ori)*bins_per_rad;  
        float mag = Mag[k]*W[k];  
  
//因为rbin > -1, 所以r0 = cvFloor(rbin)有可能是-1
        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;  
  
//根据以上公式可得:
// v_rco000 =(1-rbin)*(1-cbin)*(1-obin)
// v_rco001 =(1-rbin)*(1-cbin)*obin

//v_rco010 =(1-rbin)*cbin*(1-obin)
//v_rco011 =(1-rbin)*cbin*obin

//v_rco100 =rbin*(1-cbin)*(1-obin)
//v_rco101 =rbin*(1-cbin)*obin

// v_rco110 =rbin*cbin*(1-obin)
// v_rco111 =rbin*cbin*obin
// 因为r0,c0的范围是[-1, d)，所以要先+1
        int idx = ((r0+1)*(d+2) + c0+1)*(n+2) + o0;  
        hist[idx] += v_rco000;   ==> b1关键点
        hist[idx+1] += v_rco001;  
        hist[idx+(n+2)] += v_rco010;   ==> b3关键点
        hist[idx+(n+3)] += v_rco011;  
        hist[idx+(d+2)*(n+2)] += v_rco100; ==> b2关键点  
        hist[idx+(d+2)*(n+2)+1] += v_rco101;  
        hist[idx+(d+3)*(n+2)] += v_rco110;  ==> b4关键点 
        hist[idx+(d+3)*(n+2)+1] += v_rco111;  
    }  

    for( i = 0; i < d; i++ )  
        for( j = 0; j < d; j++ )  
        {  
//i,j各+1,因为存储的边界点不需要统计，只需要统计绿色框中蓝色端点
            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];  
        } 
    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++ )  
        {  
            // Calculate sample's histogram array coords rotated relative to ori.  
            // Subtract 0.5 so samples that fall e.g. in the center of row 1 (i.e.  
            // r_rot = 1.5) have full weight placed in row 1 after interpolation.  
        //根据旋转后的坐标,得到旋转前的坐标
 float c_rot = j * cos_t - i * sin_t;  
            float r_rot = j * sin_t + i * cos_t;  
  
//由于cos_t 和 sin_t除过了hist_width，
//所以这时旋转后的坐标的 c_rot 和 r_rot 的单位直接就是种子点了。
//下面的 rbin 和 cbin 计算的就是种子点的编号。
//种子点的编号从(-1,-1)到(d,d)。减去0.5是为了方便差值时计算权重。  
  
//rbin = r_rot + 1.5
//cbin = c_rot + 1.5
//调整原点位置,如上图的点(0,0)。
            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;  
  
//确保旋转前的坐标在有效范围内,然后计算梯度幅值和角度。
//(-1,d)的范围,其实就是上图的蓝色框框里的所有pixel。
//旋转的时候,不会影响梯度幅值,所以可以先计算出来。
            if( rbin > -1 && rbin < d && cbin > -1 && cbin < d &&  
                r > 0 && r < rows - 1 && c > 0 && c < cols - 1 )  
//超出窗边界或图像边界的点不予考虑。
//注意这里 rbin 和 cbin 的下线都是 -1,也就是说它们可能为负  
            {  
                float dx = (float)(img.at<sift_wt>(r, c+1) - img.at<sift_wt>(r, c-1));  
                float dy = (float)(img.at<sift_wt>(r-1, c) - img.at<sift_wt>(r+1, 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;  
//计算角度
    fastAtan2(Y, X, Ori, len, true); 
//计算幅值
    magnitude(X, Y, Mag, len);  
    exp(W, W, len);  
0 0