BRISK: Binary Robust Invariant Scalable Keypoints


注意:本文含有一些數學公式,如果chrome不能看見公式的話請用IE打開網站
1.特征點提取
 
特征點提取有以下幾個步驟:
a.尺度空間金字塔結構的構造
和SIFT類似,尺度空間金字塔是由不同的尺度構成,相互連續的兩個尺度之間由Octave構成. 我們令t表示尺度,它們之間的計算關系如下:
t({c_i}) = {2^i},t({d_i}) = {2^i} \cdot 1.5
圖像的大小為(width, height),舉個例子:
width,     height             scale1-octave1
(2/3)width, (2/3)height           scale1-octave2
(1/2)width, (1/2)height           scale2-octave1
(1/3)width, (1/3)height           scale2-octave2
(1/4)width, (1/4)height           scale3-octave1
(1/6)width, (1/6)height           scale3-octave2
 
其中 {d_0}{c_0}下采樣得到. 不同octave之間的采樣關系為2/3,不同尺度之間的采樣關系為1/2. 關於 構造的代碼:
void BriskScaleSpace::constructPyramid(const cv::Mat& image){

    // set correct size:
    pyramid_.clear();

    // fill the pyramid:
    pyramid_.push_back(BriskLayer(image.clone()));
    if(layers_>1){
        pyramid_.push_back(BriskLayer(pyramid_.back(),BriskLayer::CommonParams::TWOTHIRDSAMPLE));
    }
    const int octaves2=layers_;

    for(uint8_t i=2; i<octaves2; i+=2){
        pyramid_.push_back(BriskLayer(pyramid_[i-2],BriskLayer::CommonParams::HALFSAMPLE));
        pyramid_.push_back(BriskLayer(pyramid_[i-1],BriskLayer::CommonParams::HALFSAMPLE));
    }
}

 

b.通過閾值選取合適的關鍵點
關鍵點檢測是通過FAST的算法進行的. FAST和AGAST在特征點檢測中提供不同的模板。BRISK算法大部分使用FAST9-16提取特征點。首先,FAST9-16應用於每一個octave和每一層intra-octave,取相同的閾值T來分辨潛在的興趣區域。然后,對興趣區域中的點進行非極大值抑制:
1,問題點需要滿足最大值條件,就是在同層中八鄰域中的FAST score s最大。s定義為最大閾值是考慮到此點是圖像角點。
2,同層和上下層的scores 應該都比此點的score s 小。檢查等大小的方形patch內部:邊長選擇2像素。既然相鄰層是不同離散化的,就需要在patch邊緣進行插值。
3.由於確定一個特征點需要除本層以外的上下兩層,但c0層是最底層,故需虛擬有一個d-1層,但這層不使用FAST9-16,而使用FAST5-8
非極大值抑制的代碼:
  1 __inline__ bool BriskScaleSpace::isMax2D(const uint8_t layer,
  2         const int x_layer, const int y_layer){
  3     const cv::Mat& scores = pyramid_[layer].scores();
  4     const int scorescols = scores.cols;
  5     uchar* data=scores.data + y_layer*scorescols + x_layer;
  6     // decision tree:
  7     const uchar center = (*data);
  8     data--;
  9     const uchar s_10=*data;  //1
 10     if(center<s_10) return false;
 11     data+=2;
 12     const uchar s10=*data;    //2
 13     if(center<s10) return false;
 14     data-=(scorescols+1);
 15     const uchar s0_1=*data;    //3
 16     if(center<s0_1) return false;
 17     data+=2*scorescols;
 18     const uchar s01=*data;    //4
 19     if(center<s01) return false;
 20     data--;
 21     const uchar s_11=*data;    //5
 22     if(center<s_11) return false;
 23     data+=2;
 24     const uchar s11=*data;    //6
 25     if(center<s11) return false;
 26     data-=2*scorescols;
 27     const uchar s1_1=*data;    //7
 28     if(center<s1_1) return false;
 29     data-=2;
 30     const uchar s_1_1=*data;//8
 31     if(center<s_1_1) return false;
 32     
 33     /*8   3   7
 34       1   0   2
 35       5   4   6*/
 36 
 37     // reject neighbor maxima
 38     std::vector<int> delta;
 39     // put together a list of 2d-offsets to where the maximum is also reached
 40     if(center==s_1_1) { //8
 41         delta.push_back(-1);
 42         delta.push_back(-1);
 43     }
 44     if(center==s0_1) {    //3
 45         delta.push_back(0);
 46         delta.push_back(-1);
 47     }
 48     if(center==s1_1) {    //7
 49         delta.push_back(1);
 50         delta.push_back(-1);
 51     }
 52     if(center==s_10) {    //1
 53         delta.push_back(-1);
 54         delta.push_back(0);
 55     }
 56     if(center==s10) {    //2
 57         delta.push_back(1);
 58         delta.push_back(0);
 59     }
 60     if(center==s_11) {    //5
 61         delta.push_back(-1);
 62         delta.push_back(1);
 63     }
 64     if(center==s01) {    //4
 65         delta.push_back(0);
 66         delta.push_back(1);
 67     }
 68     if(center==s11) {    //6
 69         delta.push_back(1);
 70         delta.push_back(1);
 71     }
 72     const unsigned int deltasize=delta.size();
 73     if(deltasize!=0){
 74         // in this case, we have to analyze the situation more carefully:
 75         // the values are gaussian blurred and then we really decide
 76         data=scores.data + y_layer*scorescols + x_layer;
 77         int smoothedcenter=4*center+2*(s_10+s10+s0_1+s01)+s_1_1+s1_1+s_11+s11;
 78         for(unsigned int i=0; i<deltasize;i+=2){
 79             //這里把左上角作為中心點進行平滑不知道是何意?
 80             data=scores.data + (y_layer-1+delta[i+1])*scorescols + x_layer+delta[i]-1;
 81             int othercenter=*data;
 82             data++;
 83             othercenter+=2*(*data);
 84             data++;
 85             othercenter+=*data;
 86             data+=scorescols;
 87             othercenter+=2*(*data);
 88             data--;
 89             othercenter+=4*(*data);
 90             data--;
 91             othercenter+=2*(*data);
 92             data+=scorescols;
 93             othercenter+=*data;
 94             data++;
 95             othercenter+=2*(*data);
 96             data++;
 97             othercenter+=*data;
 98             if(othercenter>smoothedcenter) return false;
 99         }
100     }
101     return true;
102 }

 

 
c.去除不符合條件的關鍵點
 
 
2.特征點描述
和SIFT類似,尺度空間金字塔是由不同的尺度構成,相互連續的兩個尺度之間由Octave構成. 我們令t表示尺度,它們之間的計算關系如下:
t({c_i}) = {2^i},t({d_i}) = {2^i} \cdot 1.5
其中 {d_0}{c_0}下采樣得到. 不同octave之間的采樣關系為2/3,不同尺度之間的采樣關系為1/2.
 
 
b.通過閾值選取合適的關鍵點
關鍵點檢測是通過FAST的算法進行的. FAST和AGAST在特征點檢測中提供不同的模板。BRISK算法大部分使用FAST9-16提取特征點。首先,FAST9-16應用於每一個octave和每一層intra-octave,取相同的閾值T來分辨潛在的興趣區域。然后,對興趣區域中的點進行非極大值抑制:
1,問題點需要滿足最大值條件,就是在同層中八鄰域中的FAST score s最大。s定義為最大閾值是考慮到此點是圖像角點。
2,同層和上下層的scores 應該都比此點的score s 小。檢查等大小的方形patch內部:邊長選擇2像素。既然相鄰層是不同離散化的,就需要在patch邊緣進行插值。
3.由於確定一個特征點需要除本層以外的上下兩層,但c0層是最底層,故需虛擬有一個d-1層,但這層不使用FAST9-16,而使用FAST5-8
 
c.去除不符合條件的關鍵點
 
2.特征點描述
BRISK算法在每個模式設置了60個點。
小的藍色的圓表示在patch中的采樣位置;大的紅色虛線圓半徑為  \sigma 對應於用來平滑采樣點亮度的高斯核的標准差。上圖是尺度t=1的模式.
為了避免混疊效果,我們對在模式中的采樣點Pi應用了高斯平滑. 標准差 {\sigma _i}正比於每個采樣點對應於各自中心的距離.
然后對60個點兩兩選取組成點對。一共是(60-1)*60/2個點對。 計算局部梯度:
g({p_i},{p_j}) = ({p_j} - {p_i}) \cdot \frac{{I({p_j},{\sigma _j}) - ({p_i},{\sigma _i})}}{{{{\left\| {{p_j} - {p_i}} \right\|}^2}}}
所有點對的集合為:
A = \{ ({p_i},{p_j}) \in {\mathbb{R}^2} \times {\mathbb{R}^2}|i < N \wedge j < i \wedge i,j \in \mathbb{N}\}
 
S = \{ ({p_i},{p_j}) \in A|\left\| {{p_i} - {p_j}} \right\| < {\delta _{\max }}\} \subseteq A
 
L = \{ ({p_i},{p_j}) \in A|\left\| {{p_i} - {p_j}} \right\| > {\delta _{\min }}\} \subseteq A
 
閾值距離的選取: {\delta _{\max }} = 9.75t,{\delta _{\min }} = 13.67t
迭代整個L上的點對,這樣就可以估計出關鍵點k模式方向的整體特征:
g = \left( \begin{array}{l} {g_x}\ {g_y} \end{array} \right) = \frac{1}{L} \cdot \sum\limits_{({p_i},{p_j}) \in L} {g({p_i},{p_j})}
長距離的點對都參與了運算,基於本地梯度互相抵消的假說,所以全局梯度的計算是不必要的。這一點同時也被距離變量閾值 {\delta _{\min }}的實驗確認了
 
3.創建描述子
對於旋轉和尺度歸一化的描述子的建立,BRISK使用了關鍵點周圍的抽樣點旋轉 \alpha = \arctan 2({g_x},{g_y})角度作為模式。和BRIEF類似,BIRSK的描述子也是一個包含512個比特位的向量,每個描述子由短距離點對 (p_j^\alpha ,p_i^\alpha ) \in S兩兩進行比較產生的,上標alpha表示旋轉的模式。
 
b = \left\langle \begin{array}{l} 1,I(p_j^\alpha ,{\sigma _j}) > (p_i^\alpha ,{\sigma _i})\ 0,otherwise \end{array} \right.,\forall (p_j^\alpha ,p_i^\alpha ) \in S
 
 
與BRIEF不同的地方是,BRIEF只是進行亮度比較,除了預設尺度和預先對樣本模式的旋轉之外,BRISK和BRIEF有着根本的區別:
a.BRISK使用固定的樣本模式點,而且是以R為半徑圍繞關鍵點周圍的圓進行均勻取樣。因此特定的高斯核平滑不會突然地扭曲亮度內容的信息(模糊鄰近的兩個采樣點的亮度,從而保證亮度平滑過渡)
b.與兩兩組成的點對相比,BRISK顯著的減少了采樣點的數量(例如,單個的樣本點參與了更多的比較),限制了亮度查找表的復雜度
c.這里的比較是受到空間的限制的,所以亮度的改變僅僅只是需要局部一致性就可以了。
 
 
程序中用到的算法:
 
1.利用least square進行曲線擬合中的參數計算
二次曲線的標准公式如下:
y = a + bx + c{x^2}
給定數據:
({x_1},{y_1}),({x_2},{y_2}),({x_3},{y_3}),.......,({x_n},{y_n})
 
{y_\lambda }函數在 {x_1}的理論值,然后有:   e = {y_1} - {y_\lambda }
 
\begin{array}{l} {e_1} = {y_1} - (a + b{x_1} + cx_1^2)\ e_1^2 = {({y_1} - a - b{x_1} - cx_1^2)^2} \end{array}
 
根據最小二乘定理,當下列偏導數等於0時使得S最小.
\frac{{\partial S}}{{\partial u}} = 0,\frac{{\partial S}}{{\partial b}} = 0,and\frac{{\partial S}}{{\partial c}} = 0
 
用方程組的形式表示出來:
\sum y = na + b\sum x + c\sum {{x^2}}
 
\sum {xy} = a\sum x + b{\sum x ^2} + c\sum {{x^3}}
 
 
\sum {{x^2}y} = a{\sum x ^2} + b{\sum x ^3} + c\sum {{x^4}}
 
這里有個比較疑惑的地方, refine1D, refine1D_1, refine1D_2這三個函數的矩陣是選取什么樣的尺度初值算出來的?有知道朋友可以說說。
 
2.least square二次曲面擬合的參數計算
二次曲面的標准公式如下:
z = f(x,y) = {c_1}{x^2} + {c_2}xy + {c_3}{y^2} + {c_4}x + {c_5}y + {c_6} = C \cdot Q(x,y)
 
這里 C = ({c_1},{c_2},{c_3},{c_4},{c_5},{c_6}),Q(x,y) = ({x^2},xy,{y^2},x,y,1),選擇C使得平方誤差最小:
 
E(C) = \sum\limits_{i = 1}^m {{{(C \cdot {Q_i} - {z_i})}^2}}, 其中 {Q_i} = Q({x_i},{y_i})
 
當梯度E為0向量的時候上式取得最小值。
 
\nabla E = 2\sum\limits_{i = 1}^m {(C \cdot {Q_i} - {z_i}){Q_i} = 0}
 
也可以寫成含有6個未知變量的6個方程組的形式:
(\sum\limits_{i = 1}^m {{Q_i}{Q_i}^T} )C = \sum\limits_{i = 1}^m {{z_i}{Q_i}}
 
{Q_i}Q_i^T是6X1的矩陣 {Q_i}和1X6的矩陣 Q_i^T的乘積.
 
6x6的矩陣
A = (\sum\limits_{i = 1}^m {{Q_i}{Q_i}^T} )
6x1向量
B = \sum\limits_{i = 1}^m {{z_i}{Q_i}}
 
所以AC=B.
 
單個元素的例子:
s({x^3}y) = \sum\limits_{i = 1}^m {x_i^3{y_i}}
3.平滑函數
曲面擬合采用的是線性曲面擬合
注意,程序中的積分路線是45度角度路徑進行積分的。
4.曲面插值
smoothedIntensity,BriskLayer::value這兩個函數中使用了曲面插值算法。
 
5.重采樣
重采樣使用SSE指令對采樣進行加速。
 
 
 
 
reference:
1.<Curve Fitting and Solution of Equation>
2.<Least Squares Fitting of Data>      David Eberly
3.<BRISK: Binary Robust Invariant Scalable Keypoints> Stefan Leutenegger, Margarita Chli and Roland Y. Siegwart


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM