第一篇 特征提取以及openvslam中的相關實現詳解

尺度空間

參考尺度空間理論

金字塔

當用一個機器視覺系統分析未知場景時，計算機沒有辦法預先知道圖像中物體尺度，因此，我們需要同時考慮圖像在多尺度下的描述，獲知感興趣物體的最佳尺度。所以在很多時候，我們會在將圖像構建為一系列不同尺度的圖像集，在不同的尺度中去檢測我們感興趣的特征。比如：在Harr特征檢測人臉的時候，因為我們並不知道圖像中人臉的尺寸，所以需要生成一個不同大小的圖像組成的金字塔，掃描其中每一幅圖像來尋找可能的人臉。

圖像金字塔化的一般步驟：首先，圖像經過一個低通濾波器進行平滑（這個步驟會使圖像變模糊，好像模仿人的視覺中遠處的物體沒有近處的清晰的原理），然后，對這個平滑后的圖像進行抽樣（一般抽樣比例在水平和豎直方向上都為1/2），從而得到一系列的縮小的圖像。

數字圖像的濾波可以簡單的理解為：將圖像中各個像素點的值變更為該像素點及其相鄰區域像素點與一個濾波模板進行相乘后的值。圖像金字塔中通常我們采用高斯模板對圖像進行平滑。
高斯模板是如何生成的呢？

來看下opencv中的函數cv::getGaussianKernel()

\[G_i= \alpha *e^{-(i-( \texttt{ksize} -1)/2)^2/(2* \texttt{sigma}^2)},i=0,...,ksize-1 \]

\(\alpha\) 是尺度系數為了保證\(\sum_i G_i=1\), ksize就是模板的大小，模板通常是ksizeksize的矩形。使用這個函數我們得到的是ksize1的向量，由於gaussianBlur是一種可分離濾波器, 我們可以分別獲取x方向和y方向的高斯核模板，先對行濾波,再對列濾波的方式進行濾波。

cv::GaussianBlur(blurred_image, blurred_image, cv::Size(7, 7), 2, 2, cv::BORDER_REFLECT_101);

具體的計算分析可以看考手撕OpenCV源碼之GaussianBlur

為什么用高斯核

圖像的金字塔化能高效地（計算效率也較高）對圖像進行多尺度的表達，但它缺乏堅實的理論基礎，不能分析圖像中物體的各種尺度（雖然我們有小貓的金字塔圖像，我們還是不知道原圖像內小貓的大小）。
信號的尺度空間剛提出是就是通過一系列單參數、寬度遞增的高斯濾波器將原始信號濾波得到到組低頻信號。那么一個很明顯的疑問是：除了高斯濾波之外，其他帶有參數t的低通濾波器是否也可以用來生成一個尺度空間。

雖然很多研究者從可分性、旋轉不變性、因果性等特性推出高斯濾波器是建立線性尺度空間的最優濾波器。然后在數字圖像處理中，需要對核函數進行采樣，離散的高斯函數並不滿足連續高斯函數的的一些優良的性質。所以后來出現了一些非線性的濾波器組來建立尺度空間，如B樣條核函數。

使用高斯濾波器對圖像進行尺度空間金塔塔圖的構建，讓這個尺度空間具有下面的性質：
1）加權平均和有限孔徑效應
信號在尺度t上的表達可以看成是原信號在空間上的一系列加權平均，權重就是具有不同尺度參數的高斯核。信號在尺度t上的表達也對應於用一個無方向性的孔徑函數（特征長度為\(\sigma=\sqrt{t}\)）來觀測信號的結果。這時候信號中特征長度小於σ的精細結構會被抑制（理解為一維信號上小於σ的波動會被平滑掉）。

2）層疊平滑
也叫高斯核族的半群（Semi-Group）性質：兩個高斯核的卷積等同於另外一個不同核參數的高斯核卷積。

\[g(\mu,\sigma_1)\ast g(\mu,\sigma_2)=g(\mu,\sqrt{\sigma_1^2+\sigma_2^2}) \]

這個性質的意思就是說不同的高斯核對圖像的平滑是連續的。

3）局部極值遞性
這個特征可以從人眼的視覺原理去理解，人在看一件物體時，離得越遠，物體的細節看到的越少，細節特征是在減少的。高斯核對圖像進行濾波具有壓制局部細節的性質。

4）尺度伸縮不變性。
這里只是一個公式推導的問題，對原來的信號加一個變換函數，對變換后的信號再進行高斯核的尺度空間生成，新的信號的極值點等特征是不變的。

圖像的矩特征

參考圖像的矩特征

1. 矩的概念

圖像識別的一個核心問題是圖像的特征提取，簡單描述即為用一組簡單的數據（圖像描述量）來描述整個圖像，這組數據越簡單越有代表性越好。良好的特征不受光線、噪點、幾何形變的干擾。圖像識別發展幾十年，不斷有新的特征提出，而圖像不變矩就是其中一個。
矩是概率與統計中的一個概念，是隨機變量的一種數字特征。設X為隨機變量，c為常數，k為正整數。則量\(E[(x-c)^k]\)稱為X關於c點的k階矩。
比較重要的有兩種情況：

1. c=0。這時\(a_k=E(X^k)\)稱為X的k階原點矩

2. \(c=E(X)\)。這時\(\mu_k=E[(X-EX)^k]\)稱為X的k階中心矩。

一階原點矩就是期望。一階中心矩\(\mu_1=0\)，二階中心矩\(\mu_2\)就是X的方差\(Var(X)\)。在統計學上，高於4階的矩極少使用。\(\mu_3\)可以去衡量分布是否有偏。\(\mu_4\)可以去衡量分布（密度）在均值附近的陡峭程度如何。

針對於一幅圖像，我們把像素的坐標看成是一個二維隨機變量(X,Y)，那么一幅灰度圖像可以用二維灰度密度函數來表示，因此可以用矩來描述灰度圖像的特征。

不變矩(Invariant Moments)是一處高度濃縮的圖像特征，具有平移、灰度、尺度、旋轉不變性。M.K.Hu在1961年首先提出了不變矩的概念。1979年M.R.Teague根據正交多項式理論提出了Zernike矩。下面主要介紹這兩種矩特征的算法原理與實現。

2. Hu矩

一幅M×N的數字圖像\(f(i,j)\)，其p+q階幾何矩\(m_pq\)和中心矩\(\mu_{pq}\)為：

\[m_{pq}=\sum_{i=1}^M\sum_{j=1}^Ni^pj^qf(i,j) \\ \mu_{pq}=\sum_{i=1}^M\sum_{j=1}^N(i-\bar{i})^p(j-\bar{j})^qf(i,j) \]

其中\(f(i,j)\)為圖像在坐標點(i,j)處的灰度值。\(\bar{i}=\frac{m_{10}}{m_{00}},\bar{j}=\frac{m_{01}}{m_{00}}\)

若將\(m_{00}\)看作是圖像的灰度質量，則\((\bar{i},\bar{j})\)為圖像的質心坐標，那么中心矩\(\mu_{pq}\)反映的是圖像灰度相對於其灰度質心的分布情況。可以用幾何矩來表示中心矩，0~3階中心矩與幾何矩的關系如下：
下圖中\(\bar{x}=\bar{i},\bar{y}=\bar{j}\)

為了消除圖像比例變化帶來的影響，定義規格化中心矩如下：

\[\eta_{pq}=\frac{\mu_{pq}}{\mu_{00}^\gamma},(\gamma=1+\frac{p+q}{2},p+q \geqslant 2) \]

利用二階和三階規格中心矩可以導出下面7個不變矩組\((\Phi_1~\Phi_7)\)，它們在圖像平移、旋轉和比例變化時保持不變:

\[\begin{array}{l} \Phi_1= \eta _{20}+ \eta _{02} \\ \Phi_2=( \eta _{20}- \eta _{02})^{2}+4 \eta _{11}^{2} \\ \Phi_3=( \eta _{30}-3 \eta _{12})^{2}+ (3 \eta _{21}- \eta _{03})^{2} \\ \Phi_4=( \eta _{30}+ \eta _{12})^{2}+ ( \eta _{21}+ \eta _{03})^{2} \\ \Phi_5=( \eta _{30}-3 \eta _{12})( \eta _{30}+ \eta _{12})[( \eta _{30}+ \eta _{12})^{2}-3( \eta _{21}+ \eta _{03})^{2}]+(3 \eta _{21}- \eta _{03})( \eta _{21}+ \eta _{03})[3( \eta _{30}+ \eta _{12})^{2}-( \eta _{21}+ \eta _{03})^{2}] \\ \Phi_6 =( \eta _{20}- \eta _{02})[( \eta _{30}+ \eta _{12})^{2}- ( \eta _{21}+ \eta _{03})^{2}]+4 \eta _{11}( \eta _{30}+ \eta _{12})( \eta _{21}+ \eta _{03}) \\ \Phi_7=(3 \eta _{21}- \eta _{03})( \eta _{21}+ \eta _{03})[3( \eta _{30}+ \eta _{12})^{2}-( \eta _{21}+ \eta _{03})^{2}]-( \eta _{30}-3 \eta _{12})( \eta _{21}+ \eta _{03})[3( \eta _{30}+ \eta _{12})^{2}-( \eta _{21}+ \eta _{03})^{2}] \\ \end{array} \]

使用opencv可以方便的計算：

int main(int argc, char** argv) 
{ 
    Mat image = imread(argv[1]);  
    cvtColor(image, image, CV_BGR2GRAY); 
    Moments mts = moments(image); 
    double hu[7]; 
    HuMoments(mts, hu); 
    for (int i=0; i<7; i++) 
    { 
        cout << log(abs(hu[i])) <<endl; 
    } 
   return 0; 
}

可以看到圖像雖然發生了各種變化，但是矩的值基本保持不變。

Harris角點

參考Harris角點

基本原理

人眼對角點的識別通常是在一個局部的小區域或小窗口完成的。如果在各個方向上移動這個特征的小窗口，窗口內區域的灰度發生了較大的變化，那么就認為在窗口內遇到了角點。如果這個特定的窗口在圖像各個方向上移動時，窗口內圖像的灰度沒有發生變化，那么窗口內就不存在角點；如果窗口在某一個方向移動時，窗口內圖像的灰度發生了較大的變化，而在另一些方向上沒有發生變化，那么，窗口內的圖像可能就是一條直線的線段。

對於圖像\(I(x,y)\)，當在點\((x,y)\)處平移\((\Delta x,\Delta y)\)后的自相似性，可以通過自相關函數給出：

\[c(x,y;\Delta x,\Delta y) = \sum_{(u,v)\in W(x,y)}w(u,v)(I(u,v) – I(u+\Delta x,v+\Delta y))^2 \]

其中，\(W(x,y)\)是以點(x,y)為中心的窗口，\(w(u,v)\)為加權函數，它既可是常數，也可以是高斯加權函數。

根據泰勒展開，對圖像\(I(x,y)\)在平移\((\Delta x,\Delta y)\)后進行一階近似：

\[I(u+\Delta x,v+\Delta y) = I(u,v)+I_x(u,v)\Delta x+I_y(u,v)\Delta y+O(\Delta x^2,\Delta y^2)\approx I(u,v)+I_x(u,v)\Delta x+I_y(u,v)\Delta y \]

其中，\(I_x,I_y\)是圖像\(I(x,y)\)的偏導數，這樣的話，自相關函數則可以簡化為：

\[c(x,y;\Delta x,\Delta y)\approx \sum_w (I_x(u,v)\Delta x+I_y(u,v)\Delta y)^2=[\Delta x,\Delta y]M(x,y)\begin{bmatrix}\Delta x \\ \Delta y\end{bmatrix} \]

其中

\[M(x,y)=\sum_w \begin{bmatrix}I_x(x,y)^2&I_x(x,y)I_y(x,y) \\ I_x(x,y)I_y(x,y)&I_y(x,y)^2\end{bmatrix} = \begin{bmatrix}\sum_w I_x(x,y)^2&\sum_w I_x(x,y)I_y(x,y) \\\sum_w I_x(x,y)I_y(x,y)&\sum_w I_y(x,y)^2\end{bmatrix}=\begin{bmatrix}A&C\\C&B\end{bmatrix} \]

也就是說圖像\(I(x,y)\)在點(x,y)處平移\((\Delta x,\Delta y)\)后的自相關函數可以近似為二項函數：

\[c(x,y;\Delta x,\Delta y)\approx A\Delta x^2+2C\Delta x\Delta y+B\Delta y^2 \]

其中

\[A=\sum_w I_x^2, B=\sum_w I_y^2,C=\sum_w I_x I_y \]

二次項函數本質上就是一個橢圓函數。橢圓的扁率和尺寸是由\(M(x,y)\)的特征值\(\lambda_1、\lambda_2\)決定的，橢圓的方向是由\(M(x,y)\)的特征矢量決定的，如下圖所示，橢圓方程為：

\[[\Delta x,\Delta y]M(x,y)\begin{bmatrix}\Delta x \\ \Delta y\end{bmatrix} = 1 \]

橢圓函數特征值與圖像中的角點、直線（邊緣）和平面之間的關系如下圖所示。共可分為三種情況：

• 圖像中的直線。一個特征值大，另一個特征值小，λ1≫λ2或λ2≫λ1。自相關函數值在某一方向上大，在其他方向上小。
• 圖像中的平面。兩個特征值都小，且近似相等；自相關函數數值在各個方向上都小。
• 圖像中的角點。兩個特征值都大，且近似相等，自相關函數在所有方向都增大。
  

根據二次項函數特征值的計算公式，我們可以求M(x,y)矩陣的特征值。但是Harris給出的角點差別方法並不需要計算具體的特征值，而是計算一個角點響應值R來判斷角點。R的計算公式為：

\[R=det \boldsymbol{M} - \alpha(trace\boldsymbol{M})^2 \]

式中，\(det \boldsymbol{M}\)為矩陣\(\boldsymbol{M}=\begin{bmatrix}A&B\\B&C\end{bmatrix}\)的行列式；\(trace\boldsymbol{M}\)為矩陣M的直跡；\(\alpha\)為經常常數，取值范圍為0.04~0.06。事實上，特征是隱含在\(det \boldsymbol{M}\)和\(det \boldsymbol{M}\)中，因為，

\[det\boldsymbol{M} = \lambda_1\lambda_2=AC-B^2 \\ trace\boldsymbol{M}=\lambda_2+\lambda_2=A+C \]

特性

1. 增大α的值，將減小角點響應值R，降低角點檢測的靈性，減少被檢測角點的數量；減小α值，將增大角點響應值R，增加角點檢測的靈敏性，增加被檢測角點的數量。

2. Harris角點檢測算子對亮度和對比度的變化不敏感;
   這是因為在進行Harris角點檢測時，使用了微分算子對圖像進行微分運算，而微分運算對圖像密度的拉升或收縮和對亮度的抬高或下降不敏感。換言之，對亮度和對比度的仿射變換並不改變Harris響應的極值點出現的位置，但是，由於閾值的選擇，可能會影響角點檢測的數量。
   

3. Harris角點檢測算子具有旋轉不變性
   Harris角點檢測算子使用的是角點附近的區域灰度二階矩矩陣。而二階矩矩陣可以表示成一個橢圓，橢圓的長短軸正是二階矩矩陣特征值平方根的倒數。當特征橢圓轉動時，特征值並不發生變化，所以判斷角點響應值R也不發生變化，由此說明Harris角點檢測算子具有旋轉不變性。

4. Harris角點檢測算子不具有尺度不變性

FAST

參考FAST特征點檢測

Features from Accelerated Segment Test

1. 從圖片中選取一個像素P，下面我們將判斷它是否是一個特征點。我們首先把它的亮度值設為\(I_p\)。
2. 設定一個合適的閾值t。
3. 考慮以該像素點為中心的一個半徑等於3像素的離散化的Bresenham圓，這個圓的邊界上有16個像素。

4. 如果在這個大小為16個像素的圓上有n個連續的像素點，它們的像素值都比\(I_p+t\)大，或者都比\(I_p−t\)小，那么它就是一個角點。（如圖1中的白色虛線所示）。n的值可以設置為12或者9。

非極大值抑制

從鄰近的位置選取了多個特征點是另一個問題，我們可以使用Non-Maximal Suppression來解決。

1. 為每一個檢測到的特征點計算它的響應大小（score function）V。這里V定義為點p和它周圍16個像素點的絕對偏差的和。
2. 考慮兩個相鄰的特征點，並比較它們的V值。
3. V值較低的點將會被刪除。

BRIEF

Binary Robust Independent Elementary Features
參考BRIEF 特征描述子

BRIEF提供了一種計算二值串的捷徑，而並不需要去計算一個類似於SIFT的特征描述子。它需要先平滑圖像，然后在特征點周圍選擇一個Patch，在這個Patch內通過一種選定的方法來挑選出來\(n_d\)個點對。然后對於每一個點對\((p,q)\)，我們來比較這兩個點的亮度值，如果\(I(p)>I(q)\)則這個點對生成了二值串中一個的值為1，如果\(I(p)<I(q)\)，則對應在二值串中的值為-1，否則為0。所有\(n_d\)個點對，都進行比較之間，我們就生成了一個\(n_d\)長的二進制串。

設我們在特征點的鄰域塊大小為\(S*S\)內選擇\(n_d\)個點對\((p,q)\)，實驗中測試了5種采樣方法：

1）在圖像塊內平均采樣；

2）p和q都符合\((0,\frac{1}{25}S^2)\)的高斯分布；

3）p符合\((0,\frac{1}{25}S^2)\)的高斯分布，而q符合\((0,\frac{1}{100}S^2)\)的高斯分布；

4）在空間量化極坐標下的離散位置隨機采樣

5）把p固定為(0,0)，q在周圍平均采樣

第二種效果比較好;

ORB

Oriented FAST and Rotated BRIEF
參考ORB特征點檢測
ORB特征是將FAST特征點的檢測方法與BRIEF特征描述子結合起來，並在它們原來的基礎上做了改進與優化。

旋轉不變性

我們知道FAST特征點是沒有尺度不變性的，所以我們可以通過構建高斯金字塔，然后在每一層金字塔圖像上檢測角點，來實現尺度不變性。那么，對於局部不變性，我們還差一個問題沒有解決，就是FAST特征點不具有方向，ORB的論文中提出了一種利用灰度質心法來解決這個問題，灰度質心法假設角點的灰度與質心之間存在一個偏移，這個向量可以用於表示一個方向。對於任意一個特征點p來說，我們定義p的鄰域像素的矩為：

\[m_{pq} = \sum_{x,y}x^py^qI(x,y) \]

其中\(I(x,y)\)為點(x,y)處的灰度值。那么我們可以得到圖像的質心為：

\[C = \left(\frac{m_{10}}{m_{00}},\frac{m_{01}}{m_{00}}\right) \]

那么特征點與質心的夾角定義為FAST特征點的方向：

\[\theta = arctan(m_{01},m_{10}) \]

為了提高方法的旋轉不變性，需要確保x和y在半徑為r的圓形區域內，即\(x,y\in [-r,r]\)，r等於鄰域半徑。

描述子

ORB選擇了BRIEF作為特征描述方法，但是我們知道BRIEF是沒有旋轉不變性的，所以我們需要給BRIEF加上旋轉不變性，把這種方法稱為“Steer BREIF”。對於任何一個特征點來說，它的BRIEF描述子是一個長度為n的二值碼串，這個二值串是由特征點周圍n個點對（2n個點）生成的，現在我們將這2n個點\((x_i,y_i),i=1,2,⋯,2n\)組成一個矩陣S

\[S =\begin{pmatrix}x_1&x_2&\cdots&x_{2n}\\y_1&y_2&\cdots&y_{2n}\end{pmatrix} \]

Calonder建議為每個塊的旋轉和投影集合分別計算BRIEF描述子，但代價昂貴。ORB中采用了一個更有效的方法：使用鄰域方向θ和對應的旋轉矩陣\(R_θ\)，構建S的一個校正版本\(S_θ\)

\[S_{\theta} = R_{\theta}S \]

其中:

\[R_{\theta} = \begin{bmatrix}cos\theta & sin\theta \\ –sin\theta &cos\theta\end{bmatrix} \]

實際上，我們可以把角度離散化，即把360度分為12份，每一份是30度，然后我們對這個12個角度分別求得一個Sθ，這樣我們就創建了一個查找表，對於每一個θ，我們只需查表即可快速得到它的點對的集合\(S_θ\)。

解決描述子的區分性

BRIEF令人驚喜的特性之一是：對於n維的二值串的每個比特征位，所有特征點在該位上的值都滿足一個均值接近於0.5，而方差很大的高斯分布。方差越大，說明區分性越強，那么不同特征點的描述子就表現出來越大差異性，對匹配來說不容易誤配。但是當我們把BRIEF沿着特征點的方向調整為Steered BRIEF時，均值就漂移到一個更加分散式的模式。可以理解為有方向性的角點關鍵點對二值串則展現了一個更加均衡的表現。而且論文中提到經過PCA對各個特征向量進行分析，得知Steered BRIEF的方差很小，判別性小，各個成分之間相關性較大。
為了減少Steered BRIEF方差的虧損，並減少二進制碼串之間的相關性，ORB使用了一種學習的方法來選擇一個較小的點對集合。方法如下：
首先建立一個大約300k關鍵點的測試集，這些關鍵點來自於PASCAL2006集中的圖像。
對於這300k個關鍵點中的每一個特征點，考慮它的31×31的鄰域，我們將在這個鄰域內找一些點對。不同於BRIEF中要先對這個Patch內的點做平滑，再用以Patch中心為原點的高斯分布選擇點對的方法。ORB為了去除某些噪聲點的干擾，選擇了一個5×5大小的區域的平均灰度來代替原來一個單點的灰度，這里5×5區域內圖像平均灰度的計算可以用積分圖的方法。我們知道31×31的Patch里共有N=(31−5+1)×(31−5+1)個這種子窗口，那么我們要N個子窗口中選擇2個子窗口的話，共有\(C^2_N\)種方法。所以，對於300k中的每一個特征點，我們都可以從它的31×31大小的鄰域中提取出一個很長的二進制串，長度為\(M = C_N^2\)，表示為

\[binArray = [p_1,p_2,\cdots,p_M],p_i\in\{0,1\} \]

那么當300k個關鍵點全部進行上面的提取之后，我們就得到了一個300k×M的矩陣，矩陣中的每個元素值為0或1。
進行貪婪搜索：從T中把排在第一的那個列放到R中，T中就沒有這個點對了測試結果了。然后把T中的排下一個的列與R中的所有元素比較，計算它們的相關性，如果相關超過了某一事先設定好的閾值，就扔了它，否則就把它放到R里面。重復上面的步驟，只到R中有256個列向量為止。如果把T全找完也，也沒有找到256個，那么，可以把相關的閾值調高一些，再重試一遍。這樣，我們就得到了256個點對。上面這個過程我們稱它為rBRIEF。

openvslam特征提取

結合test/openvslam/feature/orb_extractor.cc代碼，逐步分析特征提取過程。openvslam提取的特征的一個顯著特征就是十分均勻，下文會詳細分析實現的方法。

初始化

首先定義orb參數

const auto params = feature::orb_params();

// -- src/openvslam/feature/orb_params.h

// 最多關鍵點的個數
unsigned int max_num_keypts_ = 2000;
// 圖像金字塔縮放比例
float scale_factor_ = 1.2;
// 圖像金字塔層數
unsigned int num_levels_ = 8;
// FAST腳點檢測的灰度差門限
unsigned int ini_fast_thr_ = 20;
// 如果上面的門限沒有檢測到角點，改用更小的門限值
unsigned int min_fast_thr = 7;

然后定義特征提取器

auto extractor = feature::orb_extractor(params);

此時做一些初始化

1. 計算圖像金塔每一層的縮放系數(以及倒數)、系數的平方（以及倒數）;
   

2. 根據設定的最大特征點數目計算圖像金字塔中每層對應的特征點數;

\[N*\frac{1-a}{1-a^n} + N*\frac{1-a}{1-a^n}*a + \cdots+ N*\frac{1-a}{1-a^n}*a^{n-1}\approx N \\ a=\frac{1}{scale}，N=最大特征點數，n=金字塔層數\]

這里層與層之間的特征點數滿足縮放因子，但是實際上應該滿足平方的關系，但是這樣的數列不好找。

3. 計算方向
   

提取特征點

std::vector<cv::KeyPoint> keypts;
cv::Mat desc;
extractor.extract(img, mask, keypts, desc);

1. 准備每層圖像金字塔的圖像;

//直接resize
cv::resize(image_pyramid_.at(level - 1), image_pyramid_.at(level), size, 0, 0, cv::INTER_LINEAR);

2. mask初始化;

用戶可以自定義任意個數的矩形mask區域，不進行特征點提取;也可以使用更細致的image mask;

//[x_min / cols, x_max / cols, y_min / rows, y_max / rows]
std::vector<std::vector<float>> mask_rects_;

void orb_extractor::create_rectangle_mask(const unsigned int cols, const unsigned int rows) {
    // rect_mask_和原圖有同樣的尺寸，初始填充數值為255
    if (rect_mask_.empty()) {
        rect_mask_ = cv::Mat(rows, cols, CV_8UC1, cv::Scalar(255));
    }
    // 根據mask_rects_得到mask區域
    for (const auto& mask_rect : orb_params_.mask_rects_) {
        // draw black rectangle
        const unsigned int x_min = std::round(cols * mask_rect.at(0));
        const unsigned int x_max = std::round(cols * mask_rect.at(1));
        const unsigned int y_min = std::round(rows * mask_rect.at(2));
        const unsigned int y_max = std::round(rows * mask_rect.at(3));
        // 將mask區域數值改為0
        cv::rectangle(rect_mask_, cv::Point2i(x_min, y_min), cv::Point2i(x_max, y_max), cv::Scalar(0), -1, cv::LINE_AA);
    }
}

3. 計算fast keypoints

//BRIEF orientation
static constexpr unsigned int fast_patch_size_ = 31;
//half size of FAST patch
static constexpr int fast_half_patch_size_ = fast_patch_size_ / 2;
//size of maximum ORB patch radius
static constexpr unsigned int orb_patch_radius_ = 19;

// 珊格化有效圖像區域
constexpr unsigned int cell_size = 64;
// 珊格化之后每個珊格與相鄰的珊格有一定的重疊區域
constexpr unsigned int overlap = 6;

程序流程：

遍歷圖像金字塔每一層
    根據orb_patch_radius_計算出該層的邊界，從而確定該層有效搜索區域的HxW
    按照cell_size對該有效區域珊格化，珊格中有任何在mask內，則不對該珊格進行特征點計算，剔除該珊格
    對珊格進行FAST檢測，首先使用較大的門限ini_fast_thr_，如果沒有檢測到FAST關鍵點，再使用較小的門限min_fast_thr進行檢測
    將檢測到的特征點確認不再mask中后添加到該層特征點集合中
    通過樹形結構篩選該層的特征點
    調整使用cv::FAST計算的特征點中的patch直徑，原始是7，設置為31（rBRIEF描述子的patch size）* scale_factor(金字塔層數的比例因子) 以及 _octave值
    計算FAST特征點的方向，和論文中一致，計算圖像的矩

通過樹形結構篩選特征點
    初始化根結點
    迭代展開其子節點
        一個父節點對應展開4個子葉節點 將一個矩形等分為4等份，每份為被一個子葉節點掌握，將父節點中的特征點分配代各自的子葉節點中
        展開的葉節點數 >= 該層的目標特征點數 或者 葉節點數不再發生變化則跳出迭代
    在將每個節點中響應值最大（特征最好）的特征點當作該節點上的特征點
    由於每個子葉節點均是由上層等分展開的，所得到的特征點在圖像平面上分布就十分均勻

• 關鍵點
  使用cv::FAST提取的特征點默認cv.FAST_FEATURE_DETECTOR_TYPE_9_16與論文中一致，即patch直徑為7個像素，沒有方向和金字塔層數，需要后面自己計算填寫。

/*
_pt: 關鍵點的像素坐標(x,y)
_size: 關鍵點patch直徑，像素點
_angle: 關鍵點的方向
_response： 響應值，越大越好
_octave： 關鍵點位於圖像金字塔的層數
_class_id： 關鍵點的類別ID
*/
KeyPoint(Point2f _pt, float _size, float _angle=-1, float _response=0, int _octave=0, int _class_id=-1);

• 初始化根結點
  根結點的數量與寬高比有關：
  一個根結點掌管寬高1：1的區域，通常假設寬度>高度

const unsigned int num_initial_nodes = std::round(static_cast<double>(max_x - min_x) / (max_y - min_y));

// 節點的結構
class orb_extractor_node {
public:
    orb_extractor_node() = default;

    // 一個父節點對應展開4個子葉節點 將一個矩形等分為4等份，每份為被一個子葉節點掌握
    std::array<orb_extractor_node, 4> divide_node();

    // 以該結點為父節點的特征點
    std::vector<cv::KeyPoint> keypts_;
    // 該節點掌控的區域，矩形的對角點
    cv::Point2i pt_begin_, pt_end_;

    std::list<orb_extractor_node>::iterator iter_;

    bool is_leaf_node_ = false;
};


// nodes是鏈表結構，為了方便給指定的node中添加特征點，使用vector保存每個node的指針
std::list<orb_extractor_node> nodes;
std::vector<orb_extractor_node*> initial_nodes;


// 用於判斷（x,y）是否為mash點
auto is_in_mask = [&mask](const unsigned int y, const unsigned int x, const float scale_factor) {
    return mask.at<unsigned char>(y * scale_factor, x * scale_factor) == 0;
};

• 使用openMP 加速計算:
  並行程序塊，同時只能有一個線程能訪問該並行程序塊

#pragma omp critical [(name)] //[]表示名字可選
{
    ...
}

4.計算描述子

和論文中一致，作者借鑒opencv中實現的方法。

其中點對兒的取法來自opencv

計算描述子
    圖像進行高斯平滑
    計算描述子以及方向，每個關鍵點生成一個32字節的描述子

5. 統一特征點位置

將所有圖像金字塔中的特征點位置統一到原始圖像中的相應的位置；

可優化內容

1. orb_extractor_node::divide_node() 中每個子節點預留的空間都是該父節點的特征點數，是不是太多了？

問題

1. 在提取FAST特征點的時候沒有對圖像使用高斯平滑，再作下采樣；計算描述子的時候采用了。
2. u_max_ 似乎沒有用；